Regenerate code

2 files changed, 475 insertions, 0 deletions

diff --git a/m/PARSER/McLachlan_BSSN.kranc b/m/PARSER/McLachlan_BSSN.kranc
new file mode 100644
index 0000000..745bcc1
--- /dev/null
+++ b/m/PARSER/McLachlan_BSSN.kranc
@@ -0,0 +1,200 @@
+# This is an until-end-of-line comment
+
+THORN BSSN
+
+
+
+DERIVATIVE
+  PDstandardNth[i_]    -> StandardCenteredDifferenceOperator[1,derivOrder/2,i]
+
+DERIVATIVE
+  PDstandardNth[i_,i_] -> StandardCenteredDifferenceOperator[2,derivOrder/2,i]
+DERIVATIVE
+  PDstandardNth[i_,j_] -> StandardCenteredDifferenceOperator[1,derivOrder/2,i] *
+                          StandardCenteredDifferenceOperator[1,derivOrder/2,j]
+DERIVATIVE
+  PDdissipationNth[i_] ->
+    spacing[i]^(derivOrder+1) / 2^(derivOrder+2) *
+    StandardCenteredDifferenceOperator[derivOrder+2,derivOrder/2+1,i]
+
+JACOBIAN {PD, FD, J, dJ}
+
+
+
+TENSOR normal, tangentA, tangentB, dir
+TENSOR xx, rr, th, ph
+TENSOR J, dJ
+TENSOR admg, admK, admalpha, admdtalpha, admbeta, admdtbeta, H, M
+TENSOR g, detg, gu, G, R, trR, Km, trK, cdphi, cdphi2
+TENSOR phi, gt, At, Xt, Xtn, alpha, A, beta, B, Atm, Atu, trA, Ats, trAts
+TENSOR dottrK, dotXt
+TENSOR cXt, cS, cA
+TENSOR e4phi, em4phi, ddetg, detgt, gtu, ddetgt, dgtu, ddgtu, Gtl, Gtlu, Gt
+TENSOR Rt, Rphi, gK
+TENSOR T00, T0, T, rho, S
+TENSOR x, y, z, r
+TENSOR epsdiss
+
+
+
+SYMMETRIC admg[la,lb], admK[la,lb]
+SYMMETRIC g[la,lb], K[la,lb], R[la,lb], cdphi2[la,lb]
+SYMMETRIC gt[la,lb], At[la,lb], Ats[la,lb], Rt[la,lb], Rphi[la,lb], T[la,lb]
+
+SYMMETRIC {dJ[ua,lb,lc], lb, lc}
+SYMMETRIC {G[ua,lb,lc], lb, lc}
+SYMMETRIC {Gtl[la,lb,lc], lb, lc}
+SYMMETRIC {Gt[ua,lb,lc], lb, lc}
+SYMMETRIC {gK[la,lb,lc], la, lb}
+
+SYMMETRIC gu[ua,ub], gtu[ua,ub], Atu[ua,ub]
+SYMMETRIC {dgtu[ua,ub,lc], ua, ub}
+SYMMETRIC {ddgtu[ua,ub,lc,ld], ua, ub}
+SYMMETRIC {ddgtu[ua,ub,lc,ld], lc, ld}
+
+
+
+CONNECTION {CD, PD, G}
+CONNECTION {CDt, PD, Gt}
+
+
+
+GROUP {phi      , "log_confac"}
+GROUP {gt[la,lb], "metric"    }
+GROUP {Xt[ua   ], "Gamma"     }
+
+EXTRA_GROUP {"Grid::coordinates", {x, y, z, r}}
+EXTRA_GROUP {"ADMBase::metric",  {gxx, gxy, gxz, gyy, gyz, gzz}}
+EXTRA_GROUP {"ADMBase::curv",    {kxx, kxy, kxz, kyy, kyz, kzz}}
+
+
+
+DEFINE pi = N[Pi,40]
+
+DEFINE
+  detgExpr  = Det [MatrixOfComponents [g [la,lb]]]
+DEFINE
+  ddetgExpr[la_] =
+    Sum [D[Det[MatrixOfComponents[g[la, lb]]], X] PD[X, la],
+         {X, Union[Flatten[MatrixOfComponents[g[la, lb]]]]}]
+
+DEFINE
+  detgtExpr = Det [MatrixOfComponents [gt[la,lb]]]
+DEFINE
+  ddetgtExpr[la_] =
+    Sum [D[Det[MatrixOfComponents[gt[la, lb]]], X] PD[X, la],
+         {X, Union[Flatten[MatrixOfComponents[gt[la, lb]]]]}]
+
+
+
+CALCULATION "Minkowski"
+  Schedule: {"IN ADMBase_InitialData"}
+  ConditionalOnKeyword: {"my_initial_data", "Minkowski"}
+  BEGIN EQUATIONS
+    phi       -> IfThen[conformalMethod, 1, 0]
+    phi       -> conformalMethod ? 1 : 0
+    gt[la,lb] -> KD[la,lb]
+    trK       -> 0
+    At[la,lb] -> 0
+    Xt[ua]    -> 0
+    alpha     -> 1
+    A         -> 0
+    beta[ua]  -> 0
+    B[ua]     -> 0
+  END EQUATIONS
+END CALCULATION
+
+
+
+CALCULATION "convertFromADMBase"
+  Schedule: {"AT initial AFTER ADMBase_PostInitial"}
+  ConditionalOnKeyword: {"my_initial_data", "ADMBase"}
+  SHORTHAND g[la,lb], detg, gu[ua,ub], em4phi
+  BEGIN EQUATIONS
+    g[la,lb]  -> admg[la,lb]
+    detg      -> detgExpr
+    gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]]
+    
+    phi       -> IfThen [conformalMethod, detg^(-1/6), Log[detg]/12]
+    em4phi    -> IfThen [conformalMethod, phi^2, Exp[-4 phi]]
+    gt[la,lb] -> em4phi g[la,lb]
+    
+    trK       -> gu[ua,ub] admK[la,lb]
+    At[la,lb] -> em4phi (admK[la,lb] - (1/3) g[la,lb] trK)
+    
+    alpha     -> admalpha
+    
+    beta[ua]  -> admbeta[ua]
+  END EQUATIONS
+END CALCULATION
+
+
+
+CALCULATION convertFromADMBaseGammaCalc
+  Name: BSSN <> "_convertFromADMBaseGamma"
+  Schedule: {"AT initial AFTER " <> BSSN <> "_convertFromADMBase"}
+  ConditionalOnKeyword: {"my_initial_data", "ADMBase"}
+  # Do not synchronise right after this routine; instead, synchronise
+  # after extrapolating
+  Where: Interior
+  # Synchronise after this routine, so that the refinement boundaries
+  # are set correctly before extrapolating.  (We will need to
+  # synchronise again after extrapolating because extrapolation does
+  # not fill ghost zones, but this is irrelevant here.)
+  SHORTHAND dir[ua]
+  SHORTHAND detgt, gtu[ua,ub], Gt[ua,lb,lc], theta
+  BEGIN 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]
+                    (PD[gt[lb,ld],lc] + PD[gt[lc,ld],lb] - PD[gt[lb,lc],ld])
+    Xt[ua] -> gtu[ub,uc] Gt[ua,lb,lc]
+    
+    # If LapseACoeff=0, then A is not evolved, in the sense that it
+    # does not influence the time evolution of other variables.
+    A -> IfThen [LapseACoeff != 0
+                 1 / (- harmonicF alpha^harmonicN)
+                 (+ admdtalpha
+                  - LapseAdvectionCoeff beta[ua] PDua[alpha,la]
+                  - LapseAdvectionCoeff Abs[beta[ua]] PDus[alpha,la])
+                 0]
+    
+    theta -> thetaExpr
+    
+    # If ShiftBCoeff=0 or theta ShiftGammaCoeff=0, then B^i is not
+    # evolved, in the sense that it does not influence the time
+    # evolution of other variables.
+    B[ua] -> IfThen [ShiftGammaCoeff ShiftBCoeff != 0
+                     1 / (theta ShiftGammaCoeff)
+                     (+ admdtbeta[ua]
+                      - ShiftAdvectionCoeff beta[ub] PDua[beta[ua],lb]
+                      - ShiftAdvectionCoeff Abs[beta[ub]] PDus[beta[ua],lb])
+                     0]
+  END EQUATIONS
+END CALCULATION
+
+
+
+INHERITED_IMPLEMENTATION ADMBase, TmunuBase
+
+
+
+KEYWORD_PARAMETER "my_initial_data"
+  # Visibility: "restricted"
+  # Description: "ddd"
+  AllowedValues: {"ADMBase", "Minkowski"}
+  Default: "ADMBase"
+END KEYWORD_PARAMETER
+
+
+
+REAL_PARAMETER LapseACoeff
+  Description: "Whether to evolve A in time"
+  Default: 0
+END REAL_PARAMETER
+
+
+
+END THORN
diff --git a/m/PARSER/McLachlan_BSSN.m b/m/PARSER/McLachlan_BSSN.m
new file mode 100644
index 0000000..810edaa
--- /dev/null
+++ b/m/PARSER/McLachlan_BSSN.m
@@ -0,0 +1,275 @@
+$Path = Join[$Path, {"../../../repos/Kranc/Tools/CodeGen",
+                     "../../../repos/Kranc/Tools/MathematicaMisc"}];
+
+Get["KrancThorn`"];
+
+SetEnhancedTimes[False];
+SetSourceLanguage["C"];
+
+
+
+(******************************************************************************)
+(* Derivatives *)
+(******************************************************************************)
+
+KD = KroneckerDelta;
+
+partialDerivatives =
+{
+  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],
+  PDdissipationNth[i_] ->
+    spacing[i]^(derivOrder+1) / 2^(derivOrder+2) *
+    StandardCenteredDifferenceOperator[derivOrder+2,derivOrder/2+1,i],
+};
+
+ResetJacobians;
+
+DefineJacobian[PD, FD, J, dJ];
+
+
+
+(******************************************************************************)
+(* Tensors *)
+(******************************************************************************)
+
+DefineTensor[normal];
+DefineTensor[tangentA];
+DefineTensor[tangentB];
+DefineTensor[dir];
+
+AssertSymmetricIncreasing[admg[lalb]];
+AssertSymmetricIncreasing[admK[lalb]];
+AssertSymmetricIncreasing [dJ[ua,lb,lc], lb, lc];
+AssertSymmetricIncreasing [G[ua,lb,lc], lb, lc];
+AssertSymmetricIncreasing [Gtl[la,lb,lc], lb, lc];
+AssertSymmetricIncreasing [Gt[ua,lb,lc], lb, lc];
+AssertSymmetricIncreasing [gK[la,lb,lc], la, lb];
+
+DefineConnection [CD, PD, G];
+DefineConnection [CDt, PD, Gt];
+
+
+
+(******************************************************************************)
+(* Expressions *)
+(******************************************************************************)
+
+pi = N[Pi,40];
+
+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]]]]}];
+
+detgtExpr = Det [MatrixOfComponents [gt[la,lb]]];
+ddetgtExpr[la_] =
+  Sum [D[Det[MatrixOfComponents[gt[la, lb]]], X] PD[X, la],
+       {X, Union[Flatten[MatrixOfComponents[gt[la, lb]]]]}];
+
+
+
+(******************************************************************************)
+(* Groups *)
+(******************************************************************************)
+
+declaredGroups = {
+  SetGroupName [CreateGroupFromTensor [phi      ], "log_confac"],
+  SetGroupName [CreateGroupFromTensor [gt[la,lb]], "metric"    ],
+  SetGroupName [CreateGroupFromTensor [Xt[ua]   ], "Gamma"     ]
+};
+
+declaredGroupNames = Map[First, declaredGroups];
+
+extraGroups =
+  {{"Grid::coordinates", {x, y, z, r}},
+   {"ADMBase::metric",  {gxx, gxy, gxz, gyy, gyz, gzz}},
+   {"ADMBase::curv",    {kxx, kxy, kxz, kyy, kyz, kzz}}}
+};
+
+allGroups = Join[declaredGroups, extraGroups];
+
+
+
+(******************************************************************************)
+(* Initial data *)
+(******************************************************************************)
+
+initialCalc =
+{
+  Name -> BSSN <> "_Minkowski",
+  Schedule -> {"IN ADMBase_InitialData"},
+  ConditionalOnKeyword -> {"my_initial_data", "Minkowski"},
+  Equations -> 
+  {
+    phi       -> IfThen [conformalMethod, 1, 0],
+    gt[la,lb] -> KD[la,lb],
+    trK       -> 0,
+    At[la,lb] -> 0,
+    Xt[ua]    -> 0,
+    alpha     -> 1,
+    A         -> 0,
+    beta[ua]  -> 0,
+    B[ua]     -> 0
+  }
+};
+
+
+
+(******************************************************************************)
+(* Convert from ADMBase *)
+(******************************************************************************)
+
+convertFromADMBaseCalc =
+{
+  Name -> BSSN <> "_convertFromADMBase",
+  Schedule -> {"AT initial AFTER ADMBase_PostInitial"},
+  ConditionalOnKeyword -> {"my_initial_data", "ADMBase"},
+  Shorthands -> {g[la,lb], detg, gu[ua,ub], em4phi},
+  Equations -> 
+  {
+    g[la,lb]  -> admg[la,lb],
+    detg      -> detgExpr,
+    gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]],
+    
+    phi       -> IfThen [conformalMethod, detg^(-1/6), Log[detg]/12],
+    em4phi    -> IfThen [conformalMethod, phi^2, Exp[-4 phi]],
+    gt[la,lb] -> em4phi g[la,lb],
+    
+    trK       -> gu[ua,ub] admK[la,lb],
+    At[la,lb] -> em4phi (admK[la,lb] - (1/3) g[la,lb] trK),
+    
+    alpha     -> admalpha,
+    
+    beta[ua]  -> admbeta[ua]
+  }
+};
+
+convertFromADMBaseGammaCalc =
+{
+  Name -> BSSN <> "_convertFromADMBaseGamma",
+  Schedule -> {"AT initial AFTER " <> BSSN <> "_convertFromADMBase"},
+  ConditionalOnKeyword -> {"my_initial_data", "ADMBase"},
+  (*
+  Where -> InteriorNoSync,
+  *)
+  (* Do not synchronise right after this routine; instead, synchronise
+     after extrapolating *)
+  Where -> Interior,
+  (* Synchronise after this routine, so that the refinement boundaries
+     are set correctly before extrapolating.  (We will need to
+     synchronise again after extrapolating because extrapolation does
+     not fill ghost zones, but this is irrelevant here.)  *)
+  Shorthands -> {dir[ua],
+                 detgt, gtu[ua,ub], Gt[ua,lb,lc], theta},
+  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]
+                    (PD[gt[lb,ld],lc] + PD[gt[lc,ld],lb] - PD[gt[lb,lc],ld]),
+    Xt[ua] -> gtu[ub,uc] Gt[ua,lb,lc],
+    
+(*
+    A -> - admdtalpha / (harmonicF alpha^harmonicN) (LapseAdvectionCoeff - 1),
+*)
+    (* If LapseACoeff=0, then A is not evolved, in the sense that it
+       does not influence the time evolution of other variables.  *)
+    A -> IfThen [LapseACoeff != 0,
+                 1 / (- harmonicF alpha^harmonicN)
+                 (+ admdtalpha
+                  - LapseAdvectionCoeff beta[ua] PDua[alpha,la]
+                  - LapseAdvectionCoeff Abs[beta[ua]] PDus[alpha,la]),
+                 0],
+    
+    theta -> thetaExpr,
+    
+    (* If ShiftBCoeff=0 or theta ShiftGammaCoeff=0, then B^i is not
+       evolved, in the sense that it does not influence the time
+       evolution of other variables.  *)
+    B[ua] -> IfThen [ShiftGammaCoeff ShiftBCoeff != 0,
+                     1 / (theta ShiftGammaCoeff)
+                     (+ admdtbeta[ua]
+                      - ShiftAdvectionCoeff beta[ub] PDua[beta[ua],lb]
+                      - ShiftAdvectionCoeff Abs[beta[ub]] PDus[beta[ua],lb]),
+                     0]
+  }
+};
+
+
+
+(******************************************************************************)
+(* Implementations *)
+(******************************************************************************)
+
+inheritedImplementations = {"ADMBase", "TmunuBase"};
+
+
+
+(******************************************************************************)
+(* Parameters *)
+(******************************************************************************)
+
+keywordParameters =
+{
+  {
+    Name -> "my_initial_data",
+    (* Visibility -> "restricted", *)
+    (* Description -> "ddd", *)
+    AllowedValues -> {"ADMBase", "Minkowski"},
+    Default -> "ADMBase"
+  },
+  {
+    Name -> "my_initial_boundary_condition",
+    Visibility -> "restricted",
+    (* Description -> "ddd", *)
+    AllowedValues -> {"none"},
+    Default -> "none"
+  }
+};
+
+realParameters =
+{
+  {
+    Name -> LapseACoeff,
+    Description -> "Whether to evolve A in time",
+    Default -> 0
+  },
+  {
+    Name -> harmonicF,
+    Description -> "d/dt alpha = - f alpha^n K   (harmonic=1, 1+log=2)",
+    Default -> 1
+  }
+};
+
+
+
+(******************************************************************************)
+(* Construct the thorns *)
+(******************************************************************************)
+
+calculations =
+{
+  initialCalc,
+  convertFromADMBaseCalc,
+  convertFromADMBaseGammaCalc
+};
+
+CreateKrancThornTT [allGroups, ".", BSSN,
+  Calculations -> calculations,
+  DeclaredGroups -> declaredGroupNames,
+  PartialDerivatives -> derivatives,
+  EvolutionTimelevels -> evolutionTimelevels,
+  DefaultEvolutionTimelevels -> 3,
+  UseLoopControl -> True,
+  InheritedImplementations -> inheritedImplementations,
+  InheritedKeywordParameters -> inheritedKeywordParameters,
+  ExtendedKeywordParameters -> extendedKeywordParameters,
+  KeywordParameters -> keywordParameters,
+  IntParameters -> intParameters,
+  RealParameters -> realParameters
+];