Add thorn that evaluates the ADM constraints

1 files changed, 351 insertions, 0 deletions

diff --git a/m/McLachlan_ADMConstraints.m b/m/McLachlan_ADMConstraints.m
new file mode 100644
index 0000000..a3312a3
--- /dev/null
+++ b/m/McLachlan_ADMConstraints.m
@@ -0,0 +1,351 @@
+$Path = Join[$Path, {"../../../kranc/Tools/CodeGen",
+                     "../../../kranc/Tools/MathematicaMisc"}];
+
+Get["KrancThorn`"];
+
+SetEnhancedTimes[False];
+SetSourceLanguage["C"];
+
+(******************************************************************************)
+(* Options *)
+(******************************************************************************)
+
+createCode[derivOrder_, useGlobalDerivs_, evolutionTimelevels_, addMatter_] :=
+Module[{},
+
+prefix = "ML_";
+suffix =
+  ""
+  <> If [useGlobalDerivs, "_MP", ""]
+  <> If [derivOrder!=4, "_O" <> ToString[derivOrder], ""]
+  (* <> If [evolutionTimelevels!=3, "_TL" <> ToString[evolutionTimelevels], ""] *)
+  (* <> If [addMatter==1, "_M", ""] *)
+  ;
+
+ADMConstraints = prefix <> "ADMConstraints" <> suffix;
+
+(******************************************************************************)
+(* Derivatives *)
+(******************************************************************************)
+
+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],
+  
+(* 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]
+};
+
+FD  = PDstandardNth;
+FDu = PDupwindNth;
+FDo = PDonesided;
+
+ResetJacobians;
+If [useGlobalDerivs,
+    DefineJacobian[PD, FD, J, dJ],
+    DefineJacobian[PD, FD, KD, Zero3]];
+If [useGlobalDerivs,
+    DefineJacobian[PDu, FDu, J, dJ],
+    DefineJacobian[PDu, FDu, KD, Zero3]];
+If [useGlobalDerivs,
+    DefineJacobian[PDo, FDo, J, dJ],
+    DefineJacobian[PDo, FDo, KD, Zero3]];
+
+
+
+(******************************************************************************)
+(* Tensors *)
+(******************************************************************************)
+
+(* Register the tensor quantities with the TensorTools package *)
+Map [DefineTensor,
+     {normal, tangentA, tangentB, dir,
+      nn, nu, nlen, nlen2, su, vg,
+      J, dJ,
+      g, K, alpha, beta, detg, gu, G, R, trR, Km, trK,
+      H, M,
+      T00, T0, T, rho, S,
+      x, y, z, r}];
+
+Map [AssertSymmetricIncreasing,
+     {g[la,lb], K[la,lb], R[la,lb], T[la,lb]}];
+AssertSymmetricIncreasing [dJ[ua,lb,lc], lb, lc];
+AssertSymmetricIncreasing [G[ua,lb,lc], lb, lc];
+Map [AssertSymmetricDecreasing,
+     {gu[ua,ub]}];
+
+DefineConnection [CD, PD, G];
+
+(* Use the CartGrid3D variable names *)
+x1=x; x2=y; x3=z;
+
+(* Use the ADMBase variable names *)
+g11=gxx; g12=gxy; g22=gyy; g13=gxz; g23=gyz; g33=gzz;
+K11=kxx; K12=kxy; K22=kyy; K13=kxz; K23=kyz; K33=kzz;
+alpha=alp;
+beta1=betax; beta2=betay; beta3=betaz;
+
+Map [DefineTensor,
+     {eTtt,
+      eTtx, eTty, eTtz,
+      eTxx, eTxy, eTxz, eTyy, eTyz, eTzz}];
+
+(******************************************************************************)
+(* Expressions *)
+(******************************************************************************)
+
+detgExpr = Det [MatrixOfComponents [g [la,lb]]];
+ddetgExpr[la_] =
+  Sum [D[Det[MatrixOfComponents[g[la, lb]]], X] PD[X, la],
+       {X, Union[Flatten[MatrixOfComponents[g[la, lb]]]]}];
+
+pi = N[Pi,40]; 
+
+(******************************************************************************)
+(* Groups *)
+(******************************************************************************)
+
+SetGroupTimelevels[g_,tl_] = Join[g, {Timelevels -> tl}];
+
+evolvedGroups = {};
+evaluatedGroups =
+  {SetGroupName [CreateGroupFromTensor [H    ], prefix <> "Ham"],
+   SetGroupName [CreateGroupFromTensor [M[la]], prefix <> "mom"]};
+(*
+evolvedGroups =
+  {SetGroupName [CreateGroupFromTensor [H    ], prefix <> "Ham"],
+   SetGroupName [CreateGroupFromTensor [M[la]], prefix <> "mom"]};
+evaluatedGroups = {};
+*)
+
+declaredGroups = Join [evolvedGroups, evaluatedGroups];
+declaredGroups = Map [AddGroupExtra [#, Timelevels -> evolutionTimelevels] &, declaredGroups];
+
+declaredGroupNames = Map [First, declaredGroups];
+
+
+
+extraGroups =
+  {{"ADMBase::metric",  {gxx, gxy, gxz, gyy, gyz, gzz}},
+   {"ADMBase::curv",    {kxx, kxy, kxz, kyy, kyz, kzz}},
+   {"ADMBase::lapse",   {alp}},
+   {"ADMBase::dtlapse", {dtalp}},
+   {"ADMBase::shift",   {betax, betay, betaz}},
+   {"ADMBase::dtshift", {dtbetax, dtbetay, dtbetaz}},
+   {"Grid::coordinates", {x, y, z, r}},
+   {"TmunuBase::stress_energy_scalar", {eTtt}},
+   {"TmunuBase::stress_energy_vector", {eTtx, eTty, eTtz}},
+   {"TmunuBase::stress_energy_tensor", {eTxx, eTxy, eTxz, eTyy, eTyz, eTzz}},
+   {"Coordinates::jacobian",  {J11, J12, J13, J21, J22, J23, J31, J32, J33}},
+   {"Coordinates::jacobian2", {dJ111, dJ112, dJ113, dJ122, dJ123, dJ133,
+                               dJ211, dJ212, dJ213, dJ222, dJ223, dJ233,
+                               dJ311, dJ312, dJ313, dJ322, dJ323, dJ333}}
+};
+
+
+
+groups = Join [declaredGroups, extraGroups];
+
+(******************************************************************************)
+(* Constraint equations *)
+(******************************************************************************)
+
+ADMConstraintsCalc =
+{
+  Name -> ADMConstraints,
+  (*
+  Schedule -> {"AT analysis"},
+  *)
+  Schedule -> {"AT evol AFTER MoL_Evolution"},
+  Where -> Interior,
+  Shorthands -> {detg, gu[ua,ub], G[ua,lb,lc],
+                 R[la,lb], trR, Km[la,lb], trK,
+                 T00, T0[la], T[la,lb], rho, S[la]},
+  Equations -> 
+  {
+    detg -> detgExpr,
+    gu[ua,ub] -> 1/detg detgExpr MatrixInverse[g[ua,ub]],
+    G[ua,lb,lc] -> 1/2 gu[ua,ud]
+                   (PD[g[lb,ld],lc] + PD[g[lc,ld],lb] - PD[g[lb,lc],ld]),
+    R[la,lb] -> G[u1,l2,la] G[l1,lb,u2] - G[u1,la,lb] G[l1,l2,u2]
+                + 1/2 gu[u1,u2] (- PD[g[l1,l2],la,lb] + PD[g[l1,la],l2,lb]
+                                 - PD[g[la,lb],l1,l2] + PD[g[l2,lb],l1,la]),
+    trR -> R[la,lb] gu[ua,ub],
+
+    Km[ua,lb] -> gu[ua,uc] K[lc,lb],
+    trK -> Km[ua,la],
+
+    (* Matter terms *)
+
+    T00 -> eTtt,
+    T01 -> eTtx,
+    T02 -> eTty,
+    T03 -> eTtz,
+    T11 -> eTxx,
+    T12 -> eTxy,
+    T13 -> eTxz,
+    T22 -> eTyy,
+    T23 -> eTyz,
+    T33 -> eTzz,
+
+    (* rho = n^a n^b T_ab *)
+    rho -> 1/alpha^2 (T00 - 2 beta[ui] T0[li] + beta[ui] beta[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] -> -1/alpha (T0[li] - beta[uj] T[li,lj]),
+
+    (* ADM constraints *)
+
+    H     -> + trR - Km[ua,lb] Km[ub,la] + trK^2
+             - addMatter 16 pi rho,
+    M[la] -> + gu[ub,uc] (CD[K[lc,la], lb] - CD[K[lc,lb], la])
+             - addMatter 8 pi S[la]
+  }
+};
+
+ADMConstraintsBoundaryCalc =
+{
+  Name -> ADMConstraints <> "_boundary",
+  (*
+  Schedule -> {"AT analysis AFTER " <> ADMConstraints},
+  *)
+  Schedule -> {"AT evol AFTER MoL_Evolution AFTER " <> ADMConstraints},
+  Where -> BoundaryWithGhosts,
+  Equations -> 
+  {
+    H     -> 0,
+    M[la] -> 0
+  }
+};
+
+(******************************************************************************)
+(* Implementations *)
+(******************************************************************************)
+
+inheritedImplementations =
+  Join[{"ADMBase"},
+       If [addMatter!=0, {"TmunuBase"}, {}],
+       If [useGlobalDerivs, {"Coordinates"}, {}]];
+
+(******************************************************************************)
+(* Parameters *)
+(******************************************************************************)
+
+intParameters =
+{
+  {
+    Name -> useMatter,
+    Description -> "Add matter terms",
+    AllowedValues -> {{Value -> "0", Description -> "no matter"},
+                      {Value -> "1", Description -> "matter"}},
+    Default -> addMatter
+  }
+};
+
+(******************************************************************************)
+(* Construct the thorns *)
+(******************************************************************************)
+
+calculations =
+{
+  ADMConstraintsCalc,
+  ADMConstraintsBoundaryCalc
+};
+
+CreateKrancThornTT [groups, ".", ADMConstraints,
+  Calculations -> calculations,
+  DeclaredGroups -> declaredGroupNames,
+  PartialDerivatives -> derivatives,
+  EvolutionTimelevels -> evolutionTimelevels,
+  UseLoopControl -> True,
+  InheritedImplementations -> inheritedImplementations,
+  IntParameters -> intParameters
+];
+
+];
+
+
+
+(******************************************************************************)
+(* Options *)
+(******************************************************************************)
+
+(* derivative order: 2, 4, 6, 8, ... *)
+(* useGlobalDerivs: False or True *)
+(* timelevels: 2 or 3
+   (keep this at 3; this is better chosen with a run-time parameter) *)
+(* matter: 0 or 1
+   (matter seems cheap; it should be always enabled) *)
+
+createCode[4, False, 3, 1];
+createCode[4, True,  3, 1];