McLachlan_BSSN.m: Make finite difference order a run-time parameter

1 files changed, 46 insertions, 142 deletions

diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m
index a999fb6..cc5ef93 100644
--- a/m/McLachlan_BSSN.m
+++ b/m/McLachlan_BSSN.m
@@ -33,156 +33,55 @@ KD = KroneckerDelta;
 
 derivatives =
 {
-  PDstandardNth[i_]    -> StandardCenteredDifferenceOperator[1,derivOrder/2,i],
-  PDstandardNth[i_,i_] -> StandardCenteredDifferenceOperator[2,derivOrder/2,i],
-  PDstandardNth[i_,j_] -> StandardCenteredDifferenceOperator[1,derivOrder/2,i] *
-                          StandardCenteredDifferenceOperator[1,derivOrder/2,j],
+  PDstandardNth[i_]    -> StandardCenteredDifferenceOperator[1,fdOrder/2,i],
+  PDstandardNth[i_,i_] -> StandardCenteredDifferenceOperator[2,fdOrder/2,i],
+  PDstandardNth[i_,j_] -> StandardCenteredDifferenceOperator[1,fdOrder/2,i] *
+                          StandardCenteredDifferenceOperator[1,fdOrder/2,j],
   PDdissipationNth[i_] ->
-    spacing[i]^(derivOrder+1) / 2^(derivOrder+2) *
-    StandardCenteredDifferenceOperator[derivOrder+2,derivOrder/2+1,i],
+    spacing[i]^(fdOrder+1) / 2^(fdOrder+2) *
+    StandardCenteredDifferenceOperator[fdOrder+2,fdOrder/2+1,i],
   
 (* PD: These come from my mathematica notebook
    "Upwind-Kranc-Convert.nb" that converts upwinding finite
    differencing operators generated by
    StandardUpwindDifferenceOperator into this form *)
 
-  Switch[derivOrder,
-    2,
-    PDupwindNth[1] -> (dir[1]*(-3 + 4*shift[1]^dir[1] - 
-                         shift[1]^(2*dir[1])))/(2*spacing[1]),
-    4,
-    PDupwindNth[1] -> (dir[1]*(-10 - 3/shift[1]^dir[1] + 18*shift[1]^dir[1] -
-                         6*shift[1]^(2*dir[1]) + 
-                         shift[1]^(3*dir[1])))/(12*spacing[1]),
-    6,
-    PDupwindNth[1] -> (dir[1]*(-35 + 2/shift[1]^(2*dir[1]) - 
-                         24/shift[1]^dir[1] + 80*shift[1]^dir[1] -
-                         30*shift[1]^(2*dir[1]) + 8*shift[1]^(3*dir[1]) -
-                         shift[1]^(4*dir[1])))/(60*spacing[1]),
-    8,
-    PDupwindNth[1] -> (dir[1]*(-378 - 5/shift[1]^(3*dir[1]) +
-                         60/shift[1]^(2*dir[1]) - 420/shift[1]^dir[1] +
-                         1050*shift[1]^dir[1] - 420*shift[1]^(2*dir[1]) + 
-                         140*shift[1]^(3*dir[1]) - 30*shift[1]^(4*dir[1]) +
-                         3*shift[1]^(5*dir[1])))/(840*spacing[1])],
-  Switch[derivOrder,
-    2,
-    PDupwindNth[2] -> (dir[2]*(-3 + 4*shift[2]^dir[2] - 
-                         shift[2]^(2*dir[2])))/(2*spacing[2]),
-    4,
-    PDupwindNth[2] -> (dir[2]*(-10 - 3/shift[2]^dir[2] + 18*shift[2]^dir[2] -
-                         6*shift[2]^(2*dir[2]) + 
-                         shift[2]^(3*dir[2])))/(12*spacing[2]),
-    6,
-    PDupwindNth[2] -> (dir[2]*(-35 + 2/shift[2]^(2*dir[2]) - 
-                         24/shift[2]^dir[2] + 80*shift[2]^dir[2] -
-                         30*shift[2]^(2*dir[2]) + 8*shift[2]^(3*dir[2]) -
-                         shift[2]^(4*dir[2])))/(60*spacing[2]),
-    8,
-    PDupwindNth[2] -> (dir[2]*(-378 - 5/shift[2]^(3*dir[2]) +
-                         60/shift[2]^(2*dir[2]) - 420/shift[2]^dir[2] +
-                         1050*shift[2]^dir[2] - 420*shift[2]^(2*dir[2]) + 
-                         140*shift[2]^(3*dir[2]) - 30*shift[2]^(4*dir[2]) +
-                         3*shift[2]^(5*dir[2])))/(840*spacing[2])],
-  Switch[derivOrder,
-    2,
-    PDupwindNth[3] -> (dir[3]*(-3 + 4*shift[3]^dir[3] - 
-                         shift[3]^(2*dir[3])))/(2*spacing[3]),
-    4,
-    PDupwindNth[3] -> (dir[3]*(-10 - 3/shift[3]^dir[3] + 18*shift[3]^dir[3] -
-                         6*shift[3]^(2*dir[3]) + 
-                         shift[3]^(3*dir[3])))/(12*spacing[3]),
-    6,
-    PDupwindNth[3] -> (dir[3]*(-35 + 2/shift[3]^(2*dir[3]) - 
-                         24/shift[3]^dir[3] + 80*shift[3]^dir[3] -
-                         30*shift[3]^(2*dir[3]) + 8*shift[3]^(3*dir[3]) -
-                         shift[3]^(4*dir[3])))/(60*spacing[3]),
-    8,
-    PDupwindNth[3] -> (dir[3]*(-378 - 5/shift[3]^(3*dir[3]) +
-                         60/shift[3]^(2*dir[3]) - 420/shift[3]^dir[3] +
-                         1050*shift[3]^dir[3] - 420*shift[3]^(2*dir[3]) + 
-                         140*shift[3]^(3*dir[3]) - 30*shift[3]^(4*dir[3]) +
-                         3*shift[3]^(5*dir[3])))/(840*spacing[3])],
-
-  (* TODO: make these higher order stencils *)
-  PDonesided[1] -> dir[1] (-1 + shift[1]^dir[1]) / spacing[1],
-  PDonesided[2] -> dir[2] (-1 + shift[2]^dir[2]) / spacing[2],
-  PDonesided[3] -> dir[3] (-1 + shift[3]^dir[3]) / spacing[3]
+  Sequence@@Flatten[Table[
+   {PDupwindNth[i] -> Switch[fdOrder,
+      2, (dir[i]*(-3 + 4*shift[i]^dir[i] - shift[i]^(2*dir[i])))/(2*spacing[i]),
+      4, (dir[i]*(-10 - 3/shift[i]^dir[i] + 18*shift[i]^dir[i] -
+          6*shift[i]^(2*dir[i]) + shift[i]^(3*dir[i])))/(12*spacing[i]),
+      6, (dir[i]*(-35 + 2/shift[i]^(2*dir[i]) - 24/shift[i]^dir[i] + 80*shift[i]^dir[i] -
+          30*shift[i]^(2*dir[i]) + 8*shift[i]^(3*dir[i]) - shift[i]^(4*dir[i])))/(60*spacing[i]),
+      8, (dir[i]*(-378 - 5/shift[i]^(3*dir[i]) + 60/shift[i]^(2*dir[i]) - 420/shift[i]^dir[i] +
+          1050*shift[i]^dir[i] - 420*shift[i]^(2*dir[i]) + 140*shift[i]^(3*dir[i]) - 30*shift[i]^(4*dir[i]) +
+          3*shift[i]^(5*dir[i])))/(840*spacing[i])],
+
+    PDupwindNthAnti[i] -> Switch[fdOrder,
+      2, (+1 shift[i]^(-2) -4 shift[i]^(-1) +0 shift[i]^( 0) +4 shift[i]^(+1) -1 shift[i]^(+2)) / (4 spacing[i]),
+      4, (-1 shift[i]^(-3) +6 shift[i]^(-2) -21 shift[i]^(-1 )+0 shift[i]^( 0) +21 shift[i]^(+1)
+          -6 shift[i]^(+2) +1 shift[i]^(+3)) / (24 spacing[i]),
+      6, (+1 shift[i]^(-4) -8 shift[i]^(-3) +32 shift[i]^(-2) -104 shift[i]^(-1) +0 shift[i]^( 0)
+          +104 shift[i]^(+1) -32 shift[i]^(+2) +8 shift[i]^(+3) -1 shift[i]^(+4)) / (120 spacing[i]),
+      8, (-3 shift[i]^(-5) +30 shift[i]^(-4) -145 shift[i]^(-3) +480 shift[i]^(-2) -1470 shift[i]^(-1)
+          +0 shift[i]^( 0) +1470 shift[i]^(+1) -480 shift[i]^(+2) +145 shift[i]^(+3) -30 shift[i]^(+4)
+          +3 shift[i]^(+5)) / (1680 spacing[i])],
+
+    PDupwindNthSymm[i] -> Switch[fdOrder,
+     2, (-1 shift[i]^(-2) +4 shift[i]^(-1) -6 shift[i]^( 0) +4 shift[i]^(+1) -1 shift[i]^(+2)) / (4 spacing[i]),
+     4, (+1 shift[i]^(-3) -6 shift[i]^(-2) +15 shift[i]^(-1) -20 shift[i]^( 0) +15 shift[i]^(+1)
+         -6 shift[i]^(+2) +1 shift[i]^(+3)) / (24 spacing[i]),
+     6, (-1 shift[i]^(-4) +8 shift[i]^(-3) - 28 shift[i]^(-2)+56 shift[i]^(-1)-70 shift[i]^( 0)
+         +56 shift[i]^(+1) -28 shift[i]^(+2) +8 shift[i]^(+3) -1 shift[i]^(+4)) / (120 spacing[i]),
+     8, (+3 shift[i]^(-5) -30 shift[i]^(-4) +135 shift[i]^(-3) -360 shift[i]^(-2) +630 shift[i]^(-1)
+         -756 shift[i]^( 0) +630 shift[i]^(+1) -360 shift[i]^(+2) +135 shift[i]^(+3) -30 shift[i]^(+4)
+         +3 shift[i]^(+5)) / (1680 spacing[i])],
+
+    (* TODO: make these higher order stencils *)
+    PDonesided[i] -> dir[i] (-1 + shift[i]^dir[i]) / spacing[i]} /. i->j, {j,1,3}],1]
 };
 
-derivatives = Join[derivatives,
-  Flatten[Table[
-  Switch[derivOrder,
-    2,
-    {PDupwindNthAnti[i] -> (  +1 shift[i]^(-2)
-                              -4 shift[i]^(-1)
-                              +0 shift[i]^( 0)
-                              +4 shift[i]^(+1)
-                              -1 shift[i]^(+2)) / (4 spacing[i]),
-     PDupwindNthSymm[i] -> (  -1 shift[i]^(-2)
-                              +4 shift[i]^(-1)
-                              -6 shift[i]^( 0)
-                              +4 shift[i]^(+1)
-                              -1 shift[i]^(+2)) / (4 spacing[i])},
-    4,
-    {PDupwindNthAnti[i] -> (  -1 shift[i]^(-3)
-                              +6 shift[i]^(-2)
-                             -21 shift[i]^(-1)
-                              +0 shift[i]^( 0)
-                             +21 shift[i]^(+1)
-                              -6 shift[i]^(+2)
-                              +1 shift[i]^(+3)) / (24 spacing[i]),
-     PDupwindNthSymm[i] -> (  +1 shift[i]^(-3)
-                              -6 shift[i]^(-2)
-                             +15 shift[i]^(-1)
-                             -20 shift[i]^( 0)
-                             +15 shift[i]^(+1)
-                              -6 shift[i]^(+2)
-                              +1 shift[i]^(+3)) / (24 spacing[i])},
-    6,
-    {PDupwindNthAnti[i] -> (  +1 shift[i]^(-4)
-                              -8 shift[i]^(-3)
-                             +32 shift[i]^(-2)
-                            -104 shift[i]^(-1)
-                              +0 shift[i]^( 0)
-                            +104 shift[i]^(+1)
-                             -32 shift[i]^(+2)
-                              +8 shift[i]^(+3)
-                              -1 shift[i]^(+4)) / (120 spacing[i]),
-     PDupwindNthSymm[i] -> (  -1 shift[i]^(-4)
-                              +8 shift[i]^(-3)
-                             -28 shift[i]^(-2)
-                             +56 shift[i]^(-1)
-                             -70 shift[i]^( 0)
-                             +56 shift[i]^(+1)
-                             -28 shift[i]^(+2)
-                              +8 shift[i]^(+3)
-                              -1 shift[i]^(+4)) / (120 spacing[i])},
-    8,
-    {PDupwindNthAnti[i] -> (   -3 shift[i]^(-5)
-                              +30 shift[i]^(-4)
-                             -145 shift[i]^(-3)
-                             +480 shift[i]^(-2)
-                            -1470 shift[i]^(-1)
-                               +0 shift[i]^( 0)
-                            +1470 shift[i]^(+1)
-                             -480 shift[i]^(+2)
-                             +145 shift[i]^(+3)
-                              -30 shift[i]^(+4)
-                               +3 shift[i]^(+5)) / (1680 spacing[i]),
-     PDupwindNthSymm[i] -> (   +3 shift[i]^(-5)
-                              -30 shift[i]^(-4)
-                             +135 shift[i]^(-3)
-                             -360 shift[i]^(-2)
-                             +630 shift[i]^(-1)
-                             -756 shift[i]^( 0)
-                             +630 shift[i]^(+1)
-                             -360 shift[i]^(+2)
-                             +135 shift[i]^(+3)
-                              -30 shift[i]^(+4)
-                               +3 shift[i]^(+5)) / (1680 spacing[i])}],
-  {i,1,3}]]
-];
+(* Put[derivatives, "derivatives.m"]; *)
 
 PD     = PDstandardNth;
 PDu    = PDupwindNth;
@@ -1242,6 +1141,11 @@ intParameters =
     AllowedValues -> {{Value -> "0", Description -> "phi method"},
                       {Value -> "1", Description -> "W method"}},
     Default -> 0
+  },
+  {
+    Name -> fdOrder,
+    Default -> derivOrder,
+    AllowedValues -> {2,4,6,8}
   }
 };
 
@@ -1346,7 +1250,7 @@ CreateKrancThornTT [groups, ".", BSSN,
   PartialDerivatives -> derivatives,
   EvolutionTimelevels -> evolutionTimelevels,
   DefaultEvolutionTimelevels -> 3,
-  UseJacobian -> useJacobian,
+  UseJacobian -> True,
   UseLoopControl -> True,
   UseVectors -> True,
   InheritedImplementations -> inheritedImplementations,