Use correct one-sided stencils at boundaries.

Some reformatting, some commented out code removed.

1 files changed, 38 insertions, 50 deletions

diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m
index 4c3cfd0..75be42c 100644
--- a/m/McLachlan_BSSN.m
+++ b/m/McLachlan_BSSN.m
@@ -35,9 +35,10 @@ derivatives =
   PDstandardNth[i_,j_] -> StandardCenteredDifferenceOperator[1,derivOrder/2,i] *
                           StandardCenteredDifferenceOperator[1,derivOrder/2,j],
   
-(* PD: These comes from my mathematica notebook: Upwind-Kranc-Convert.nb that
-converts upwinding finite differencing operators generated by
-StandardUpwindDifferenceOperator into this form *)
+(* 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] - 
@@ -94,21 +95,15 @@ StandardUpwindDifferenceOperator into this form *)
                         140*shift[3]^(3*dir[3]) - 30*shift[3]^(4*dir[3]) +
                          3*shift[3]^(5*dir[3])))/(840*spacing[3])],
 
-(*  PDupwindpNth[i_]      -> StandardUpwindDifferenceOperator[1,derivOrder/2-1,derivOrder/2+1,i],
-  PDupwindmNth[i_]      -> StandardUpwindDifferenceOperator[1,derivOrder/2+1,derivOrder/2-1,i], *)
-
-(*
-  PDPlus    [i_] -> (+1) (-1 + shift[i]^(+1)) / spacing[i],
-  PDMinus   [i_] -> (-1) (-1 + shift[i]^(-1)) / spacing[i]
-*)
-  PDPlus    [i_] -> StandardUpwindDifferenceOperator[1,0,derivOrder/2,i],
-  PDMinus   [i_] -> StandardUpwindDifferenceOperator[1,derivOrder/2,0,i]
+  (* 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]
 };
 
-FD = PDstandardNth;
+FD  = PDstandardNth;
 FDu = PDupwindNth;
-(*FDp = PDupwindpNth;
-FDm = PDupwindmNth; *)
+FDo = PDonesided;
 
 ResetJacobians;
 If [useGlobalDerivs,
@@ -117,17 +112,9 @@ If [useGlobalDerivs,
 If [useGlobalDerivs,
     DefineJacobian[PDu, FDu, J, dJ],
     DefineJacobian[PDu, FDu, KD, Zero3]];
-(*If [useGlobalDerivs,
-    DefineJacobian[PDp, FDp, J, dJ],
-    DefineJacobian[PDp, FDp, KD, Zero3]];
 If [useGlobalDerivs,
-    DefineJacobian[PDm, FDm, J, dJ],
-    DefineJacobian[PDm, FDm, KD, Zero3]]; *)
-
-(*PD = FD;
-PDu = FDu;*)
-(*PDp = FDp;
-PDm = FDm; *)
+    DefineJacobian[PDo, FDo, J, dJ],
+    DefineJacobian[PDo, FDo, KD, Zero3]];
 
 
 
@@ -339,11 +326,12 @@ convertFromADMBaseGammaCalc =
   Where -> Interior,
   (* should not sync Gamma, since boundary conditions and
      synchronisation are applied later anyway *)
-  Shorthands -> {detgt, gtu[ua,ub], Gt[ua,lb,lc], dir[ua]},
+  Shorthands -> {dir[ua],
+                 detgt, gtu[ua,ub], Gt[ua,lb,lc]},
   Equations -> 
   {
-
     dir[ua] -> Sign[beta[ua]],
+    
     detgt        -> 1 (* detgtExpr *),
     gtu[ua,ub]   -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
     Gt[ua,lb,lc] -> 1/2 gtu[ua,ud]
@@ -370,11 +358,13 @@ convertToADMBaseCalc =
   Name -> BSSN <> "_convertToADMBase",
   Schedule -> {"IN " <> BSSN <> "_convertToADMBaseGroup"},
   Where -> Interior,
-  Shorthands -> {e4phi, g[la,lb], K[la,lb], dir[ua]},
+  Shorthands -> {dir[ua],
+                 e4phi, g[la,lb], K[la,lb]},
   Equations -> 
   {
-    dir[ua]  -> Sign[beta[ua]],
-    e4phi    -> xp[4 phi],
+    dir[ua] -> Sign[beta[ua]],
+    
+    e4phi    -> Exp[4 phi],
     g[la,lb] -> e4phi gt[la,lb],
     gxx      -> g11,
     gxy      -> g12,
@@ -487,16 +477,18 @@ evolCalc =
   Name -> BSSN <> "_RHS",
   Schedule -> {"IN " <> BSSN <> "_evolCalcGroup"},
   Where -> Interior,
-  Shorthands -> {detgt, ddetgt[la], gtu[ua,ub],
+  Shorthands -> {dir[ua],
+                 detgt, ddetgt[la], gtu[ua,ub],
                  dgtu[ua,ub,lc], ddgtu[ua,ub,lc,ld], Gt[ua,lb,lc],
                  Xtn[ua], Rt[la,lb], Rphi[la,lb], R[la,lb],
                  Atm[ua,lb], Atu[ua,ub],
                  e4phi, em4phi, g[la,lb], detg,
                  ddetg[la], gu[ua,ub], G[ua,lb,lc], Ats[la,lb], trAts,
-                 T00, T0[la], T[la,lb], rho, S[la], trS, dir[ua]},
+                 T00, T0[la], T[la,lb], rho, S[la], trS},
   Equations -> 
   {
     dir[ua] -> Sign[beta[ua]],
+    
     detgt        -> 1 (* detgtExpr *),
     ddetgt[la]   -> 0 (* ddetgtExpr[la] *),
     
@@ -638,11 +630,14 @@ RHSBoundaryCalc =
   Schedule -> {"IN MoL_RHSBoundaries"},
   ConditionalOnKeyword -> {"my_rhs_boundary_condition", "radiative"},
   Where -> Boundary,
-  Shorthands -> {detgt, gtu[ua,ub], em4phi, gu[ua,ub],
+  Shorthands -> {dir[ua],
+                 detgt, gtu[ua,ub], em4phi, gu[ua,ub],
                  nn[la], nu[ua], nlen, nlen2, su[ua],
                  vg},
   Equations -> 
   {
+    dir[ua] -> Sign[normal[ua]],
+    
     detgt      -> 1 (* detgtExpr *),
     gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
     em4phi     -> Exp[-4 phi],
@@ -656,18 +651,15 @@ RHSBoundaryCalc =
     
     vg -> Sqrt[harmonicF],
     
-    (* TODO: Use PDPlus at lower boundary *)
-    (* TODO: Use standard PD operators in directions which are not a
-       boundary *)
-    dot[phi]       -> - vg su[uc] PDMinus[phi      ,lc],
-    dot[gt[la,lb]] -> -    su[uc] PDMinus[gt[la,lb],lc],
-    dot[trK]       -> - vg su[uc] PDMinus[trK      ,lc],
-    dot[At[la,lb]] -> -    su[uc] PDMinus[At[la,lb],lc],
-    dot[Xt[ua]]    -> -    su[uc] PDMinus[Xt[ua]   ,lc],
-    dot[alpha]     -> - vg su[uc] PDMinus[alpha    ,lc],
-    dot[A]         -> - vg su[uc] PDMinus[A        ,lc],
-    dot[beta[ua]]  -> -    su[uc] PDMinus[beta[ua] ,lc],
-    dot[B[ua]]     -> -    su[uc] PDMinus[B[ua]    ,lc]
+    dot[phi]       -> - vg su[uc] PDo[phi      ,lc],
+    dot[gt[la,lb]] -> -    su[uc] PDo[gt[la,lb],lc],
+    dot[trK]       -> - vg su[uc] PDo[trK      ,lc],
+    dot[At[la,lb]] -> -    su[uc] PDo[At[la,lb],lc],
+    dot[Xt[ua]]    -> -    su[uc] PDo[Xt[ua]   ,lc],
+    dot[alpha]     -> - vg su[uc] PDo[alpha    ,lc],
+    dot[A]         -> - vg su[uc] PDo[A        ,lc],
+    dot[beta[ua]]  -> -    su[uc] PDo[beta[ua] ,lc],
+    dot[B[ua]]     -> -    su[uc] PDo[B[ua]    ,lc]
   }
 };
 
@@ -1012,9 +1004,5 @@ CreateKrancThornTT [groups, ".", BSSN,
    (matter seems cheap; it should be always enabled) *)
 
 createCode[4, False, 3, 0];
-(* createCode[4, False, 3, 0]; 
-createCode[4, True, 3, 0];
-createCode[6, False, 3, 0];
-createCode[8, False, 3, 0];
 createCode[4, False, 3, 1];
-createCode[4, True,  3, 0]; *)
+createCode[4, True,  3, 0];