Correct schedule, in particular for checkpoint/recovery

Do not mark ADMBase variables for non-checkpointing if they have
multiple timelevels. (Variables with multiple timelevels must always
be checkpointed, because the past timelevels cannot be regenerated
after recovery.)

Finally remove all perl post-processing of the auto-generated code;
instead, use proper Kranc mechanisms.

Schedule the ADM constraints and ADM quantities after MoL_PostStep,
since this is where the ADMBase variables are set.

Schedule enforcing the BSSN constraints in the new schedule group
MoL_PostStepModify, since they should not be enforced after recovery.
(This would lead to inconsistencies at floating-point round-off
level.)

Regenerate all thorns.

1 files changed, 28 insertions, 4 deletions

diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m
index b05c8aa..927b894 100644
--- a/m/McLachlan_BSSN.m
+++ b/m/McLachlan_BSSN.m
@@ -808,7 +808,6 @@ evolCalc2 = PartialCalculation[evolCalc, "2",
     dot[At[la,lb]]
   }];
 
-
 dissCalc =
 {
   Name -> BSSN <> "_Dissipation",
@@ -826,6 +825,24 @@ dissCalc =
   }
 };
 
+dissCalcs =
+Table[
+{
+  Name -> BSSN <> "_Dissipation_" <> ToString[var /. {Tensor[n_,__] -> n}],
+  Schedule -> {"IN " <> BSSN <> "_evolCalcGroup " <>
+               "AFTER (" <> BSSN <> "_RHS1 " <> BSSN <> "_RHS2)"},
+  ConditionalOnKeyword -> {"apply_dissipation", "always"},
+  Where -> InteriorNoSync,
+  Shorthands -> {epsdiss[ua]},
+  Equations ->
+  {
+    epsdiss[ua] -> EpsDiss,
+    dot[var]    -> dot[var] + epsdiss[ux] PDdiss[var,lx]
+  }
+},
+  {var, {phi, gt[la,lb], Xt[ui], trK, At[la,lb], alpha, A, beta[ua], B[ua]}}
+];
+
 RHSStaticBoundaryCalc =
 {
   Name -> BSSN <> "_RHSStaticBoundary",
@@ -910,11 +927,14 @@ RHSRadiativeBoundaryCalc =
 enforceCalc =
 {
   Name -> BSSN <> "_enforce",
-  Schedule -> {"IN MoL_PostStep BEFORE " <> BSSN <> "_SelectBoundConds"},
+  Schedule -> {"IN MoL_PostStepModify"},
   Shorthands -> {detgt, gtu[ua,ub], trAt},
   Equations -> 
   {
-    (* Enforcing the constraints needs to be a projection, because it
+    (* The following comment is still interesting, but is not correct
+       any more since it is now scheduled in MoL_PostStepModify instead:
+
+       Enforcing the constraints needs to be a projection, because it
        is applied in MoL_PostStep and may thus be applied multiple
        times, not only during time evolution. Therefore detgt has to
        be calculated correctly, without assuming that det gt_ij = 1,
@@ -963,6 +983,7 @@ constraintsCalc =
 {
   Name -> BSSN <> "_constraints",
   Schedule -> Automatic,
+  After -> "MoL_PostStep",
   Where -> Interior,
   Shorthands -> {detgt, ddetgt[la], gtu[ua,ub],
                  Gt[ua,lb,lc], Gtl[la,lb,lc], Gtlu[la,lb,uc], Xtn[ua],
@@ -1293,6 +1314,7 @@ realParameters =
 (******************************************************************************)
 
 calculations =
+Join[
 {
   initialCalc,
   convertFromADMBaseCalc,
@@ -1314,7 +1336,9 @@ calculations =
   convertToADMBaseFakeDtLapseShiftCalc,
   (* constraintsCalc, *)
   constraintsCalc1, constraintsCalc2
-};
+},
+  {} (*dissCalcs*)
+];
 
 CreateKrancThornTT [groups, ".", BSSN,
   Calculations -> calculations,