Move RHS evaluations around to optimize for 32Kb instruction cache machines.

Move RHS evaluations for trK and the lapse to RHS1 to better even out the
instruction count between RHS1 and RHS2 leading to better performance on
32Kb instruction cache machines. Also replace covariant derivatives with
respect to the physical metric with covariant derivatives with respect to
the conformal metric. This avoids having to compute the physical Christoffel
symbols. The performance on 64Kb instruction cache machines is unchanged.
Also get rid of some unused shorthands and generate code for 3 timelevels.

Signed-off-by: Peter Diener <diener@cct.lsu.edu>

1 files changed, 53 insertions, 38 deletions

diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m
index 9647ccf..877a879 100644
--- a/m/McLachlan_BSSN.m
+++ b/m/McLachlan_BSSN.m
@@ -447,7 +447,7 @@ evolCalc =
                  Gt[ua,lb,lc], Xtn[ua], Rt[la,lb], Rphi[la,lb], R[la,lb],
                  Atm[ua,lb], Atu[ua,ub],
                  e4phi, em4phi, cdphi[la], cdphi2[la,lb], g[la,lb], detg,
-                 gu[ua,ub], G[ua,lb,lc], Ats[la,lb], trAts, eta,
+                 gu[ua,ub], Ats[la,lb], trAts, eta,
                  rho, S[la], trS, fac1, fac2},
   Equations -> 
   {
@@ -494,9 +494,6 @@ evolCalc =
     detg        -> detgExpr,
     (* gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]], *)
     gu[ua,ub]   -> em4phi gtu[ua,ub],
-    G[ua,lb,lc] -> Gt[ua,lb,lc]
-                   + 2 (KD[ua,lb] cdphi[lc] + KD[ua,lc] cdphi[lb] 
-                        - gtu[ua,ud] gt[lb,lc] cdphi[ld]),
     
     R[la,lb] -> Rt[la,lb] + Rphi[la,lb],
     
@@ -540,7 +537,9 @@ evolCalc =
                       + addMatter (- 16 pi alpha gtu[ui,uj] S[lj]),
     
     (* PRD 62, 044034 (2000), eqn. (11) *)
-    dot[trK]       -> - gu[ua,ub] CD[alpha,la,lb]
+    dot[trK]       -> - em4phi ( gtu[ua,ub] ( PD[alpha,la,lb]
+                                + 2 cdphi[la] PD[alpha,lb] )
+                                - Xtn[ua] PD[alpha,la] )
                       + alpha (Atm[ua,lb] Atm[ub,la] + (1/3) trK^2)
                       + beta[ua] PDu[trK,la]
     (* Equation (4.21) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *)
@@ -548,7 +547,9 @@ evolCalc =
 
     (* PRD 62, 044034 (2000), eqn. (12) *)
     (* TODO: Should we use the Hamiltonian constraint to make Rij tracefree? *)
-    Ats[la,lb]     -> - CD[alpha,la,lb] + alpha R[la,lb],
+    Ats[la,lb]     -> - CDt[alpha,la,lb] +
+                      + 2 (PD[alpha,la] cdphi[lb] + PD[alpha,lb] cdphi[la] )
+                      + alpha R[la,lb],
     trAts          -> gu[ua,ub] Ats[la,lb],
     dot[At[la,lb]] -> + em4phi (+ Ats[la,lb] - (1/3) g[la,lb] trAts )
                       + alpha (trK At[la,lb] - 2 At[la,lc] Atm[uc,lb])
@@ -594,10 +595,10 @@ evol1Calc =
   Where -> InteriorNoSync,
   Shorthands -> {dir[ua],
                  detgt, gtu[ua,ub],
-                 Gt[ua,lb,lc], Xtn[ua], Rt[la,lb], Rphi[la,lb], R[la,lb],
+                 Gt[ua,lb,lc], Xtn[ua],
                  Atm[ua,lb], Atu[ua,ub],
-                 e4phi, em4phi, cdphi[la], cdphi2[la,lb], g[la,lb], detg, eta,
-                 rho, S[la], trS, fac1, fac2},
+                 e4phi, em4phi, cdphi[la], eta,
+                 rho, S[la], trS, fac1},
   Equations -> 
   {
     dir[ua] -> Sign[beta[ua]],
@@ -619,13 +620,35 @@ evol1Calc =
     Atm[ua,lb] -> gtu[ua,uc] At[lc,lb],
     Atu[ua,ub] -> Atm[ua,lc] gtu[ub,uc],
     
+    e4phi       -> IfThen [conformalMethod, 1/phi^2, Exp[4 phi]],
+    em4phi      -> 1 / e4phi,
+
     (* Matter terms *)
     
+    (* rho = n^a n^b T_ab *)
+    rho -> addMatter
+           (1/alpha^2 (T00 - 2 beta[ui] T0[li] + beta[ui] beta[uj] T[li,lj])),
+
+    (* trS = gamma^ij T_ij  *)
+    trS -> addMatter (em4phi gtu[ui,uj] T[li,lj]),
+
     (* S_i = -p^a_i n^b T_ab, where p^a_i = delta^a_i + n^a n_i *)
     S[li] -> addMatter (-1/alpha (T0[li] - beta[uj] T[li,lj])),
     
     (* RHS terms *)
     
+    (* PRD 62, 044034 (2000), eqn. (11) *)
+    (* Note, Covariant derivatives with respect to the physical metric
+       has been replaced with derivatives with respect to the
+       conformal metric *)
+    dot[trK]       -> - em4phi ( gtu[ua,ub] ( PD[alpha,la,lb]
+                                + 2 cdphi[la] PD[alpha,lb] )
+                                - Xtn[ua] PD[alpha,la] )
+                      + alpha (Atm[ua,lb] Atm[ub,la] + (1/3) trK^2)
+                      + beta[ua] PDu[trK,la]
+    (* Equation (4.21) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *)
+                      + addMatter (4 pi alpha (rho + trS)),
+
     (* PRD 62, 044034 (2000), eqn. (10) *)
     (* PRD 67 084023 (2003), eqn. (16) and (23) *)
     dot[phi]       -> IfThen [conformalMethod, 1/3 phi, -1/6] alpha trK
@@ -637,6 +660,7 @@ evol1Calc =
                       + beta[uc] PDu[gt[la,lb],lc]
                       + gt[la,lc] PD[beta[uc],lb] + gt[lb,lc] PD[beta[uc],la]
                       - (2/3) gt[la,lb] PD[beta[uc],lc],
+
     (* PRD 62, 044034 (2000), eqn. (20) *)
     (* PRD 67 084023 (2003), eqn (26) *)
     dot[Xt[ui]]    -> - 2 Atu[ui,uj] PD[alpha,lj]
@@ -654,6 +678,14 @@ evol1Calc =
     eta -> BetaDriver
            IfThen [r > SpatialBetaDriverRadius, SpatialBetaDriverRadius / r, 1],
     
+    (* dot[alpha] -> - harmonicF alpha^harmonicN trK, *)
+    (* dot[alpha] -> - harmonicF alpha^harmonicN A + Lie[alpha, beta], *)
+    dot[alpha] -> - harmonicF alpha^harmonicN (
+                    (1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK)
+                  + LapseAdvectionCoeff beta[ua] PDu[alpha,la],
+
+    dot[A]     -> (1 - LapseAdvectionCoeff) (dot[trK] - AlphaDriver A),
+
     (* dot[beta[ua]] -> eta Xt[ua], *)
     (* dot[beta[ua]] -> ShiftGammaCoeff alpha^ShiftAlphaPower B[ua], *)
     dot[beta[ua]] -> + ShiftGammaCoeff B[ua]
@@ -682,7 +714,7 @@ evol2Calc =
                  Rt[la,lb], Rphi[la,lb], R[la,lb],
                  Atm[ua,lb],
                  e4phi, em4phi, cdphi[la], cdphi2[la,lb], g[la,lb], detg,
-                 gu[ua,ub], G[ua,lb,lc], Ats[la,lb], trAts,
+                 gu[ua,ub], Ats[la,lb], trAts,
                  rho, S[la], trS, fac1, fac2},
   Equations -> 
   {
@@ -729,33 +761,24 @@ evol2Calc =
     g[la,lb]    -> e4phi gt[la,lb],
     (* gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]], *)
     gu[ua,ub]   -> em4phi gtu[ua,ub],
-    G[ua,lb,lc] -> Gt[ua,lb,lc]
-                   + 2 (KD[ua,lb] cdphi[lc] + KD[ua,lc] cdphi[lb] 
-                        - gtu[ua,ud] gt[lb,lc] cdphi[ld]),
     
     R[la,lb] -> Rt[la,lb] + Rphi[la,lb],
     
     (* Matter terms *)
     
-    (* rho = n^a n^b T_ab *)
-    rho -> addMatter
-           (1/alpha^2 (T00 - 2 beta[ui] T0[li] + beta[ui] beta[uj] T[li,lj])),
-    
     (* trS = gamma^ij T_ij  *)
     trS -> addMatter (gu[ui,uj] T[li,lj]),
     
     (* RHS terms *)
     
-    (* PRD 62, 044034 (2000), eqn. (11) *)
-    dot[trK]       -> - gu[ua,ub] CD[alpha,la,lb]
-                      + alpha (Atm[ua,lb] Atm[ub,la] + (1/3) trK^2)
-                      + beta[ua] PDu[trK,la]
-    (* Equation (4.21) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *)
-                      + addMatter (4 pi alpha (rho + trS)),
-
     (* PRD 62, 044034 (2000), eqn. (12) *)
     (* TODO: Should we use the Hamiltonian constraint to make Rij tracefree? *)
-    Ats[la,lb]     -> - CD[alpha,la,lb] + alpha R[la,lb],
+    (* Note, Covariant derivatives with respect to the physical metric
+       has been replaced with derivatives with respect to the
+       conformal metric *)
+    Ats[la,lb]     -> - CDt[alpha,la,lb] +
+                      + 2 (PD[alpha,la] cdphi[lb] + PD[alpha,lb] cdphi[la] )
+                      + alpha R[la,lb],
     trAts          -> gu[ua,ub] Ats[la,lb],
     dot[At[la,lb]] -> + em4phi (+ Ats[la,lb] - (1/3) g[la,lb] trAts )
                       + alpha (trK At[la,lb] - 2 At[la,lc] Atm[uc,lb])
@@ -764,15 +787,7 @@ evol2Calc =
                       - (2/3) At[la,lb] PD[beta[uc],lc]
     (* Equation (4.23) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *)
                       + addMatter (- em4phi alpha 8 pi
-                                     (T[la,lb] - (1/3) g[la,lb] trS)),
-    
-    (* dot[alpha] -> - harmonicF alpha^harmonicN trK, *)
-    (* dot[alpha] -> - harmonicF alpha^harmonicN A + Lie[alpha, beta], *)
-    dot[alpha] -> - harmonicF alpha^harmonicN (
-                    (1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK)
-                  + LapseAdvectionCoeff beta[ua] PDu[alpha,la],
-
-    dot[A]     -> (1 - LapseAdvectionCoeff) (dot[trK] - AlphaDriver A)
+                                     (T[la,lb] - (1/3) g[la,lb] trS))
   }
 };
 
@@ -1226,7 +1241,7 @@ CreateKrancThornTT [groups, ".", BSSN,
 (* matter: 0 or 1
    (matter seems cheap; it should be always enabled) *)
 
-createCode[2, False, 4, 1];
-createCode[4, False, 4, 1];
-createCode[4, True,  4, 1];
-createCode[8, False, 4, 1];
+createCode[2, False, 3, 1]; 
+createCode[4, False, 3, 1];
+createCode[4, True,  3, 1];
+createCode[8, False, 3, 1];