ianhinder / Kranc

A Mathematica package for generating code for solving time dependent partial differential equations
http://kranccode.org
GNU General Public License v2.0
29 stars 10 forks source link

Support more symmetries #81

Open eschnett opened 12 years ago

eschnett commented 12 years ago

To handle the Riemann tensor, at least support antisymmetric tensors and tensors with pairwise symmetric indices, as in

R_abcd = - R_bacd R_abcd = R_cdab

This should speed up WeylScal4 quite a bit, including re-generating and building this thorn.

barrywardell commented 12 years ago

I manually implemented the symmetries of Riemann in one of my own Kranc scripts using the following code:

R[i_, j_, k_, l_] /; i > j := -R[j, i, k, l];
R[i_, j_, k_, l_] /; i == j := 0;
R[i_, j_, k_, l_] /; k > l := -R[i, j, l, k];
R[i_, j_, k_, l_] /; k == l := 0;
R[i_, j_, k_, l_] /; i > k := R[k, l, i, j];
R[i_, j_, k_, l_] /; i == k && j > l := R[k, l, i, j];

This reduces the 4-Riemann tensor to 21 components. The elimination of the last one requires the Bianchi identity which is probably a lot more complicated to implement in a general way.

Another option would be to make more use of xTensor, which would support the symmetries natively.

eschnett commented 12 years ago

How can I use this in WeylScal4? Where do these replacements go?

ianhinder commented 12 years ago

I think these could go in WeylScal4.m after the existing symmetry call (AssertSymmetricDecreasing[gamma[ua,lb,lc], lb, lc]). Kranc expands tensor expressions R[la,lb,lc,ld] first as R[1,2,3,4] etc, then converts to symbols R1234. Once in the second form, if Mathematica has been given replacement rules such as those that Barry suggested, it will automatically use them. In fact, this is (almost) how symmetries are implemented internally. This will lead to the required simplification of the right hand sides of the expressions, but will leave duplicated entries for the R[la,lb,lc,ld] -> ... lines. I think that Kranc will eliminate these duplicates (within MakeExplicit of TensorTools) but I'm not 100% sure of this.

A more interesting question is whether or not these symmetries are appropriate, and for which tensors. In WeylScal4, we have R[la,lb,lc,ld] which is the 3-Riemann. This should definitely have all the symmetries applied. We also have R4p[li,lj,lk,ll] which is the 4-Riemann projected into the slice on all its indices. I think this should have the same symmetries. We also have Rojo[lj,ll] which is the 4-Riemann projected in the normal direction on two of its indices, and into the slice on the remaining two. We should work out which symmetries this resulting contraction has, and explicitly specify them in a manner similar to what Barry proposed.

We have test cases for WeylScal4, so it should be possible to make this change and verify that the results are correct.

eschnett commented 12 years ago

No, Kranc does in general not delete duplicate lines:

x -> x+1,
x -> x+1
ianhinder commented 12 years ago

It's not the duplicate lines in a calculation which would be deleted, it's the duplicate elements created in MakeExplicit in TensorTools. See line 477 on TensorTools.m (https://github.com/ianhinder/Kranc/blob/master/Tools/CodeGen/TensorTools.m):

MakeExplicit[l:List[Rule[_, _] ..]] :=
  Flatten[Map[removeDuplicatesFromMap, Map[MakeExplicit, l]],1];

From looking at the code, this should remove the rules in expanded lists of tensor assignments which have the same left hand side. In fact, this surely must happen, because it works for the symmetries supported by Kranc, and there is no other mechanism for removing these duplicates.

eschnett commented 12 years ago

It turns out that the method suggested above does not work, as this applies the symmetries also to the shorthands, leading to zeros and minus signs in the shorthands.

I believe this will require a Kranc modification instead.

ianhinder commented 12 years ago

This is probably because we are now dealing with antisymmetric tensors for the first time, so it's not just that duplicates need to be removed, but also negative components and zeroes. One solution would be to modify RemoveDuplicates in TensorTools.m to remove anything that was 0 or matches Times[-1, _Symbol]. At that point, the function should be renamed as RemoveDuplicatesFromSymmetries or something, and we should check that it isn't used anywhere where this behaviour would not work.

I also notice that there is a RemoveDuplicates defined in Differencing.m, and I'm not sure why these don't conflict!

Barry: was this not a problem for you when you did this?

barrywardell commented 12 years ago

In my particular case, I didn't have time/was too lazy to do things properly by modifying Kranc. Instead I avoided using the tensor support in Kranc altogether and did something like the following (I do not advocate actually doing this, but maybe it will help inform the best way to do things properly in Kranc)

(* Create a list of independent components of the 4-Riemann tensor *)
riemannComps = 
  Reap[Do[If[i!=k||l>=j, Sow[{i,j,k,l}]], {i,4}, {j,i+1,4}, {k,i,4}, {l,k+1,4}]][[2, 1]];

(* Start counting indices from 0, not 1 *)
riemannComps     -= 1;

(* Define certain components in terms of others using symmetries of the Riemann tensor *)
R[i_, j_, k_, l_] /; i > j := -R[j, i, k, l];
R[i_, j_, k_, l_] /; i == j := 0;
R[i_, j_, k_, l_] /; k > l := -R[i, j, l, k];
R[i_, j_, k_, l_] /; k == l := 0;
R[i_, j_, k_, l_] /; i > k := R[k, l, i, j];
R[i_, j_, k_, l_] /; i == k && j > l := R[k, l, i, j];

R[i_, j_, k_, l_] := Symbol["R"<>ToString[i]<>ToString[j]<>ToString[k]<>ToString[l]];

(* Create a group for all of the components *)
componentList[t_, l_List] := Map[Apply[t, #]&, l];

tensorGroups = 
  {...,
   {"Riemann", Flatten@componentList[R, riemannComps]},
   ...
  };

declaredGroupNames = Map[First, tensorGroups];

(* Calculate the Riemann tensor *)
calc = {
  ...,
  Equations -> Flatten@
    {
      ...,
      (* Riemann *)
      Map[With[{i=#[[1]], j=#[[2]], k=#[[3]], l=#[[4]]},
        R[i,j,k,l] -> 
          - 1/2 ddg[j,l,i,k] + 1/2 ddg[i,l,j,k]
          + 1/2 ddg[j,k,i,l] - 1/2 ddg[i,k,j,l]
          + Sum[Gd[j,k,m] G[i,l,m] - Gd[j,l,m] G[i,k,m], {m,0,3}]] &,
        riemannComps]
    }
  };
ianhinder commented 12 years ago

Erik, I can't get the effect that you saw. For some reason, my rules for R don't seem to have any effect. Can you send me a patch to WeylScal4.m to show me exactly how you did it?

ianhinder commented 12 years ago

You have to put the rules on Tensor[R, i, j, k, l] instead of R[i,j,k,l] because the latter is converted into the former. I have implemented this in WeylScal4 now (https://trac.einsteintoolkit.org/ticket/1036) and added code to Kranc to handle the zeroes and negatives from the antisymmetric tensor components. It would probably be good for Kranc to provide a friendly interface to implement symmetries like this, so I will leave this ticket open.