Cadabra hangs on long expression

Question

Hi Troops, Here is a code that computes the Ricci tensor explicitly in terms of the metric gab and its inverse. That works fine. The second part of the code evaluates the Rab as a Sympy string (which I'll later use to generate C-code, that bit is not is this code).

The problem is that when Cadabra hits the last line it hangs. I let it run for over 7 hours with out a result (the memory usage also crept up slowly form 41Mbytes to about 115Mbytes).

I admit that I'm asking a lot here -- the expression for Rab has many thousands of terms. Am I asking too much from Cadabra? Or is there a better way to do this job?

Cheers, Leo

{a,b,c,d,e,f,i,j,k,l,m,n,o,p,q,r,s,t,u#}::Indices(position=independent,values={x,y,z}).
{x,y,z}::Coordinate.

\partial{#}::PartialDerivative;

g_{a b}::Symmetric;
g^{a b}::Symmetric;
g_{a}^{b}::KroneckerDelta.
g^{a}_{b}::KroneckerDelta.

g_{a b}::Depends(\partial{#});
g^{a b}::Depends(\partial{#});

KDelta := g_{a b} g^{b c} -> g_{a}^{c}.

Gamma := \Gamma^{a}_{b c} ->
         (1/2) g^{a e} (   \partial_{b}{g_{e c}}
                         + \partial_{c}{g_{b e}}
                         - \partial_{e}{g_{b c}}).

Riem := R^{a}_{b c d} ->
        \partial_{c}{\Gamma^{a}_{b d}} + \Gamma^{a}_{e c} \Gamma^{e}_{b d}
      - \partial_{d}{\Gamma^{a}_{b c}} - \Gamma^{a}_{e d} \Gamma^{e}_{b c}.

Ric := R_{a b} -> g^{c d} g_{a e} R^{e}_{c b d}.

Rab := R_{a b}.

substitute (Rab, Ric)
substitute (Rab, Riem)
substitute (Rab, Gamma)

distribute   (Rab)
product_rule (Rab)
distribute   (Rab)
substitute   (Rab, KDelta)
eliminate_kronecker (Rab)

# the Ricci tensor in terms of the metric
canonicalise (Rab);

# now start preparing to "evaluate" Rab

substitute (Rab, $\partial_{a b}{g_{c d}} -> dg_{c d a b}$)
substitute (Rab, $\partial_{a}{g_{b c}} -> dg_{b c a}$)
substitute (Rab, $\partial_{a}{g^{b c}} -> dg^{b c}_{a}$)

gup=[]
gdn=[]
Rdn=[]

for a in ["x","y","z"]:
   for b in ["x","y","z"]:
      gup.append("g^{"+a+" "+b+"} = gup"+a+b)
      gdn.append("g_{"+a+" "+b+"} = gdn"+a+b)
      Rdn.append("R_{"+a+" "+b+"} = R"+a+b)

gup_rl = Ex(','.join(gup))
gdn_rl = Ex(','.join(gdn))
Rdn_rl = Ex(','.join(Rdn))

gupd1=[]
gdnd1=[]

for a in ["x","y","z"]:
   for b in ["x","y","z"]:
      for c in ["x","y","z"]:
         gupd1.append("dg^{"+a+" "+b+"}_{"+c+"} = gup"+a+b+c)
         gdnd1.append("dg_{"+a+" "+b+" "+c+"} = gdn"+a+b+c)

gupd1_rl = Ex(','.join(gupd1))
gdnd1_rl = Ex(','.join(gdnd1))

gdnd2=[]

for a in ["x","y","z"]:
   for b in ["x","y","z"]:
      for c in ["x","y","z"]:
         for d in ["x","y","z"]:
            gdnd2.append("dg_{"+a+" "+b+" "+c+" "+d+"} = gdn"+a+b+c+d)

gdnd2_rl = Ex(','.join(gdnd2))

# the final expression for Rab but ...
#    Cadabra (or sympy) hangs (:

evaluate (Rab,gup_rl+gdn_rl+gupd1_rl+gdnd1_rl+gdnd2_rl,rhsonly=True);

Answer 1

With the current version on github, try

evaluate(Rab, simplify=False);

directly after the canonicalise(Rab) line. Runs in seconds (and it's not optimised yet in any form). Suggestions for further improvements are welcome.

Comments

Excellent news! Thanks so much :)

---

Hi Kasper,

I've tested the new code and it runs very quickly, great! When I coupled this to my C-code-export (from a previous post) my C-code contained lines like

x0 = UP(z);
x1 = g(x0, x0);

x8 = UP(x);
x9 = g(x8, x8);

I guess this is because the "evaluate (Rab,simplify=False)" command was placed before the rules to convert g_{a b} to gxx, gzy etc. So I removed the line "evaluate(Rab, simplify=False)" and added the ",simplify=False" to the evaluate commands after all substitutions had been made. This worked perfectly. I now have raw C-ode to evaluate the Ricci tensor in terms of the metric and its derivatives. Excellent. Thanks so much :).

I'm not sure at that this point what further improvements I would want, but I'm sure I'll think of something soon :)

Cheers, Leo

---

Very good. Can you perhaps say in a few lines what you want to do with this result? Some generic numerical evolution/evaluation of 3d spacetimes?

---

Yes, that's exactly what I'll be using the code for. In computational general relativity one of the key tasks is to compute estimates of the Riemann curvatures given the point values of the metric. This is not something anybody would reasonable code by hand (the number of terms in the full expansion is in the 10's of thousands). So being able to code this using Cadabra is a huge step. But also being able to start with an ode/pde and then use Cadabra/Python to automate the construction of the numerical code means I can easily explore alternative formulations of the field equations. This idea clearly applies to other ares of physics. Thanks again.

Answer 2

The problem here is that this makes zillions of round-trips through sympy with increasingly length expressions which cannot be simplified, but for which sympy still takes time to figure out that it can't do anything. The Cadabra side is fast.

What I should do is add a 'sympy=False' type of option to evaluate to avoid all these round trips. Will have a look at that, but it will take a few days as I'm very busy.

Comments

Thanks for that. I'm in no rush for a fix, plenty of things here to keep me busy. Cheers, Leo

---

In case someone missed it from the posts above: there is now a simplify=False option to prevent evaluate from going through sympy for every step of the computation. This speeds things up massively.