mathematica blank
mathematica code SetEnhancedTimes[False];
mathematica code SetSourceLanguage["C"];
mathematica blank
mathematica comment (******************************************************************************)
mathematica comment (* Options *)
mathematica comment (******************************************************************************)
mathematica blank
mathematica code createCode[derivOrder_, useJacobian_, splitUpwindDerivs_, useVectors_, useOpenCL_, evolutionTimelevels_, addMatter_, formulation_] :=
mathematica code Module[{prefix, suffix, thorn},
mathematica blank
mathematica code prefix = "ML_";
mathematica code suffix =
mathematica code ""
mathematica code <> If [useJacobian, "_MP", ""]
mathematica code <> If [derivOrder!=4, "_O" <> ToString[derivOrder], ""]
mathematica code <> If [splitUpwindDerivs, "", "_UPW"]
mathematica code <> If [useVectors, "", "_NV"]
mathematica code <> If [useOpenCL, "_CL", ""]
mathematica comment (* <> If [evolutionTimelevels!=3, "_TL" <> ToString[evolutionTimelevels], ""] *)
mathematica comment (* <> If [addMatter==1, "_M", ""] *)
mathematica code ;
mathematica blank
mathematica code thorn = prefix <> formulation <> suffix;
mathematica blank
mathematica code SetAttributes[IfCCZ4, HoldAll];
mathematica code IfCCZ4[expr_, else_:Sequence[]] := If[formulation === "CCZ4", expr, Unevaluated[else]];
mathematica blank
mathematica comment (******************************************************************************)
mathematica comment (* Derivatives *)
mathematica comment (******************************************************************************)
mathematica blank
mathematica code KD = KroneckerDelta;
mathematica blank
mathematica code derivatives =
mathematica code {
mathematica code PDstandardNth[i_] -> StandardCenteredDifferenceOperator[1,fdOrder/2,i],
mathematica code PDstandardNth[i_,i_] -> StandardCenteredDifferenceOperator[2,fdOrder/2,i],
mathematica code PDstandardNth[i_,j_] -> StandardCenteredDifferenceOperator[1,fdOrder/2,i] *
mathematica code StandardCenteredDifferenceOperator[1,fdOrder/2,j],
mathematica code PDdissipationNth[i_] ->
mathematica code (-1)^(fdOrder/2) *
mathematica code spacing[i]^(fdOrder+1) / 2^(fdOrder+2) *
mathematica code StandardCenteredDifferenceOperator[fdOrder+2,fdOrder/2+1,i],
mathematica blank
mathematica comment (* PD: These come from my mathematica notebook
mathematica comment "Upwind-Kranc-Convert.nb" that converts upwinding finite
mathematica comment differencing operators generated by
mathematica comment StandardUpwindDifferenceOperator into this form *)
mathematica blank
mathematica code Sequence@@Flatten[Table[
mathematica code {PDupwindNth[i] -> Switch[fdOrder,
mathematica code 2, (dir[i]*(-3 + 4*shift[i]^dir[i] - shift[i]^(2*dir[i])))/(2*spacing[i]),
mathematica code 4, (dir[i]*(-10 - 3/shift[i]^dir[i] + 18*shift[i]^dir[i] -
mathematica code 6*shift[i]^(2*dir[i]) + shift[i]^(3*dir[i])))/(12*spacing[i]),
mathematica code 6, (dir[i]*(-35 + 2/shift[i]^(2*dir[i]) - 24/shift[i]^dir[i] + 80*shift[i]^dir[i] -
mathematica code 30*shift[i]^(2*dir[i]) + 8*shift[i]^(3*dir[i]) - shift[i]^(4*dir[i])))/(60*spacing[i]),
mathematica code 8, (dir[i]*(-378 - 5/shift[i]^(3*dir[i]) + 60/shift[i]^(2*dir[i]) - 420/shift[i]^dir[i] +
mathematica code 1050*shift[i]^dir[i] - 420*shift[i]^(2*dir[i]) + 140*shift[i]^(3*dir[i]) - 30*shift[i]^(4*dir[i]) +
mathematica code 3*shift[i]^(5*dir[i])))/(840*spacing[i])],
mathematica blank
mathematica code PDupwindNthAnti[i] -> Switch[fdOrder,
mathematica code 2, (+1 shift[i]^(-2) -4 shift[i]^(-1) +0 shift[i]^( 0) +4 shift[i]^(+1) -1 shift[i]^(+2)) / (4 spacing[i]),
mathematica code 4, (-1 shift[i]^(-3) +6 shift[i]^(-2) -21 shift[i]^(-1 )+0 shift[i]^( 0) +21 shift[i]^(+1)
mathematica code -6 shift[i]^(+2) +1 shift[i]^(+3)) / (24 spacing[i]),
mathematica code 6, (+1 shift[i]^(-4) -8 shift[i]^(-3) +32 shift[i]^(-2) -104 shift[i]^(-1) +0 shift[i]^( 0)
mathematica code +104 shift[i]^(+1) -32 shift[i]^(+2) +8 shift[i]^(+3) -1 shift[i]^(+4)) / (120 spacing[i]),
mathematica code 8, (-3 shift[i]^(-5) +30 shift[i]^(-4) -145 shift[i]^(-3) +480 shift[i]^(-2) -1470 shift[i]^(-1)
mathematica code +0 shift[i]^( 0) +1470 shift[i]^(+1) -480 shift[i]^(+2) +145 shift[i]^(+3) -30 shift[i]^(+4)
mathematica code +3 shift[i]^(+5)) / (1680 spacing[i])],
mathematica blank
mathematica code PDupwindNthSymm[i] -> Switch[fdOrder,
mathematica code 2, (-1 shift[i]^(-2) +4 shift[i]^(-1) -6 shift[i]^( 0) +4 shift[i]^(+1) -1 shift[i]^(+2)) / (4 spacing[i]),
mathematica code 4, (+1 shift[i]^(-3) -6 shift[i]^(-2) +15 shift[i]^(-1) -20 shift[i]^( 0) +15 shift[i]^(+1)
mathematica code -6 shift[i]^(+2) +1 shift[i]^(+3)) / (24 spacing[i]),
mathematica code 6, (-1 shift[i]^(-4) +8 shift[i]^(-3) - 28 shift[i]^(-2)+56 shift[i]^(-1)-70 shift[i]^( 0)
mathematica code +56 shift[i]^(+1) -28 shift[i]^(+2) +8 shift[i]^(+3) -1 shift[i]^(+4)) / (120 spacing[i]),
mathematica code 8, (+3 shift[i]^(-5) -30 shift[i]^(-4) +135 shift[i]^(-3) -360 shift[i]^(-2) +630 shift[i]^(-1)
mathematica code -756 shift[i]^( 0) +630 shift[i]^(+1) -360 shift[i]^(+2) +135 shift[i]^(+3) -30 shift[i]^(+4)
mathematica code +3 shift[i]^(+5)) / (1680 spacing[i])],
mathematica blank
mathematica comment (* TODO: make these higher order stencils *)
mathematica code PDonesided[i] -> dir[i] (-1 + shift[i]^dir[i]) / spacing[i]} /. i->j, {j,1,3}],1]
mathematica code };
mathematica blank
mathematica code PD = PDstandardNth;
mathematica code PDu = PDupwindNth;
mathematica code PDua = PDupwindNthAnti;
mathematica code PDus = PDupwindNthSymm;
mathematica comment (* PDo = PDonesided; *)
mathematica code PDdiss = PDdissipationNth;
mathematica blank
mathematica code If [splitUpwindDerivs,
mathematica code Upwind[dir_, var_, idx_] := dir PDua[var,idx] + Abs[dir] PDus[var,idx],
mathematica code Upwind[dir_, var_, idx_] := dir PDu[var,idx]];
mathematica blank
mathematica blank
mathematica blank
mathematica comment (******************************************************************************)
mathematica comment (* Tensors *)
mathematica comment (******************************************************************************)
mathematica blank
mathematica comment (* Register the tensor quantities with the TensorTools package *)
mathematica code Map [DefineTensor,
mathematica code {normal, tangentA, tangentB, dir,
mathematica code nn, nu, nlen, nlen2, su, vg,
mathematica code xx, rr, th, ph,
mathematica code admg, admK, admalpha, admdtalpha, admbeta, admdtbeta, H, M,
mathematica code g, detg, gu, G, R, trR, Km, trK, cdphi, cdphi2,
mathematica code phi, gt, At, Xt, Xtn, Theta, Z,
mathematica code alpha, A, beta, B, Atm, Atu, trA, Ats, trAts,
mathematica code dottrK, dotXt,
mathematica code cXt, cS, cA,
mathematica code e4phi, em4phi, ddetg, detgt, gtu, ddetgt, dgtu, ddgtu, Gtl, Gtlu, Gt,
mathematica code Rt, Rphi, gK,
mathematica code T00, T0, T, rho, S,
mathematica code x, y, z, r,
mathematica code epsdiss}];
mathematica blank
mathematica comment (* NOTE: It seems as if Lie[.,.] did not take these tensor weights
mathematica comment into account. Presumably, CD[.,.] and CDt[.,.] don't do this either. *)
mathematica code SetTensorAttribute[phi, TensorWeight, +1/6];
mathematica code SetTensorAttribute[gt, TensorWeight, -2/3];
mathematica code SetTensorAttribute[Xt, TensorWeight, +2/3];
mathematica code SetTensorAttribute[At, TensorWeight, -2/3];
mathematica code SetTensorAttribute[cXt, TensorWeight, +2/3];
mathematica code SetTensorAttribute[cS, TensorWeight, +2 ];
mathematica blank
mathematica code Map [AssertSymmetricIncreasing,
mathematica code {admg[la,lb], admK[la,lb], g[la,lb], K[la,lb], R[la,lb], cdphi2[la,lb],
mathematica code gt[la,lb], At[la,lb], Ats[la,lb], Rt[la,lb], Rphi[la,lb], T[la,lb]}];
mathematica code AssertSymmetricIncreasing [G[ua,lb,lc], lb, lc];
mathematica code AssertSymmetricIncreasing [Gtl[la,lb,lc], lb, lc];
mathematica code AssertSymmetricIncreasing [Gt[ua,lb,lc], lb, lc];
mathematica code AssertSymmetricIncreasing [gK[la,lb,lc], la, lb];
mathematica code Map [AssertSymmetricIncreasing,
mathematica code {gu[ua,ub], gtu[ua,ub], Atu[ua,ub]}];
mathematica code AssertSymmetricIncreasing [dgtu[ua,ub,lc], ua, ub];
mathematica code AssertSymmetricIncreasing [ddgtu[ua,ub,lc,ld], ua, ub];
mathematica code AssertSymmetricIncreasing [ddgtu[ua,ub,lc,ld], lc, ld];
mathematica blank
mathematica code DefineConnection [CD, PD, G];
mathematica code DefineConnection [CDt, PD, Gt];
mathematica blank
mathematica comment (* Use the CartGrid3D variable names *)
mathematica code x1=x; x2=y; x3=z;
mathematica blank
mathematica comment (* Use the ADMBase variable names *)
mathematica code admg11=gxx; admg12=gxy; admg22=gyy; admg13=gxz; admg23=gyz; admg33=gzz;
mathematica code admK11=kxx; admK12=kxy; admK22=kyy; admK13=kxz; admK23=kyz; admK33=kzz;
mathematica code admalpha=alp;
mathematica code admdtalpha=dtalp;
mathematica code admbeta1=betax; admbeta2=betay; admbeta3=betaz;
mathematica code admdtbeta1=dtbetax; admdtbeta2=dtbetay; admdtbeta3=dtbetaz;
mathematica blank
mathematica comment (* Use the TmunuBase variable names *)
mathematica code T00=eTtt;
mathematica code T01=eTtx; T02=eTty; T03=eTtz;
mathematica code T11=eTxx; T12=eTxy; T22=eTyy; T13=eTxz; T23=eTyz; T33=eTzz;
mathematica blank
mathematica blank
mathematica blank
mathematica comment (******************************************************************************)
mathematica comment (* Expressions *)
mathematica comment (******************************************************************************)
mathematica blank
mathematica comment (* enum constants for conformalMethod; these must be consistent
mathematica comment with the definition of the Cactus parameter conformalMethod *)
mathematica code CMphi = 0;
mathematica code CMW = 1;
mathematica blank
mathematica code detgExpr = Det [MatrixOfComponents [g [la,lb]]];
mathematica code ddetgExpr[la_] =
mathematica code Sum [D[Det[MatrixOfComponents[g[la, lb]]], X] PD[X, la],
mathematica code {X, Union[Flatten[MatrixOfComponents[g[la, lb]]]]}];
mathematica blank
mathematica code detgtExpr = Det [MatrixOfComponents [gt[la,lb]]];
mathematica code ddetgtExpr[la_] =
mathematica code Sum [D[Det[MatrixOfComponents[gt[la, lb]]], X] PD[X, la],
mathematica code {X, Union[Flatten[MatrixOfComponents[gt[la, lb]]]]}];
mathematica blank
mathematica code etaExpr = SpatialBetaDriverRadius / Max [r, SpatialBetaDriverRadius];
mathematica code thetaExpr = Min [Exp [1 - r / SpatialShiftGammaCoeffRadius], 1];
mathematica blank
mathematica blank
mathematica blank
mathematica comment (******************************************************************************)
mathematica comment (* Groups *)
mathematica comment (******************************************************************************)
mathematica blank
mathematica code evolvedGroups =
mathematica code {SetGroupName [CreateGroupFromTensor [phi ], prefix <> "log_confac"],
mathematica code SetGroupName [CreateGroupFromTensor [gt[la,lb]], prefix <> "metric" ],
mathematica code SetGroupName [CreateGroupFromTensor [Xt[ua] ], prefix <> "Gamma" ],
mathematica code SetGroupName [CreateGroupFromTensor [trK ], prefix <> "trace_curv"],
mathematica code SetGroupName [CreateGroupFromTensor [At[la,lb]], prefix <> "curv" ],
mathematica code SetGroupName [CreateGroupFromTensor [alpha ], prefix <> "lapse" ],
mathematica code SetGroupName [CreateGroupFromTensor [A ], prefix <> "dtlapse" ],
mathematica code SetGroupName [CreateGroupFromTensor [beta[ua] ], prefix <> "shift" ],
mathematica code SetGroupName [CreateGroupFromTensor [B[ua] ], prefix <> "dtshift" ],
mathematica code IfCCZ4[SetGroupName[CreateGroupFromTensor[Theta], prefix <> "Theta"]]};
mathematica code evaluatedGroups =
mathematica code {SetGroupName [CreateGroupFromTensor [H ], prefix <> "Ham"],
mathematica code SetGroupName [CreateGroupFromTensor [M[la] ], prefix <> "mom"],
mathematica code SetGroupName [CreateGroupFromTensor [cS ], prefix <> "cons_detg"],
mathematica code SetGroupName [CreateGroupFromTensor [cXt[ua]], prefix <> "cons_Gamma"],
mathematica code SetGroupName [CreateGroupFromTensor [cA ], prefix <> "cons_traceA"]};
mathematica blank
mathematica code declaredGroups = Join [evolvedGroups, evaluatedGroups];
mathematica code declaredGroupNames = Map [First, declaredGroups];
mathematica blank
mathematica blank
mathematica blank
mathematica code extraGroups =
mathematica code {{"grid::coordinates", {x, y, z, r}},
mathematica code {"ADMBase::metric", {gxx, gxy, gxz, gyy, gyz, gzz}},
mathematica code {"ADMBase::curv", {kxx, kxy, kxz, kyy, kyz, kzz}},
mathematica code {"ADMBase::lapse", {alp}},
mathematica code {"ADMBase::dtlapse", {dtalp}},
mathematica code {"ADMBase::shift", {betax, betay, betaz}},
mathematica code {"ADMBase::dtshift", {dtbetax, dtbetay, dtbetaz}},
mathematica code {"TmunuBase::stress_energy_scalar", {eTtt}},
mathematica code {"TmunuBase::stress_energy_vector", {eTtx, eTty, eTtz}},
mathematica code {"TmunuBase::stress_energy_tensor", {eTxx, eTxy, eTxz, eTyy, eTyz, eTzz}}
mathematica code };
mathematica blank
mathematica code groups = Join [declaredGroups, extraGroups];
mathematica blank
mathematica blank
mathematica blank
mathematica comment (******************************************************************************)
mathematica comment (* Initial data *)
mathematica comment (******************************************************************************)
mathematica blank
mathematica code initialCalc =
mathematica code {
mathematica code Name -> thorn <> "_Minkowski",
mathematica code Schedule -> {"IN ADMBase_InitialData"},
mathematica code ConditionalOnKeyword -> {"my_initial_data", "Minkowski"},
mathematica code Equations ->
mathematica code {
mathematica code phi -> IfThen[conformalMethod==CMW, 1, 0],
mathematica code gt[la,lb] -> KD[la,lb],
mathematica code trK -> 0,
mathematica code At[la,lb] -> 0,
mathematica code Xt[ua] -> 0,
mathematica code alpha -> 1,
mathematica code A -> 0,
mathematica code beta[ua] -> 0,
mathematica code B[ua] -> 0,
mathematica code IfCCZ4[Theta -> 0]
mathematica code }
mathematica code };
mathematica blank
mathematica blank
mathematica blank
mathematica comment (******************************************************************************)
mathematica comment (* Split a calculation *)
mathematica comment (******************************************************************************)
mathematica blank
mathematica code PartialCalculation[calc_, suffix_, updates_, evolVars_] :=
mathematica code Module[
mathematica code {name, calc1, replaces, calc2, vars, patterns, eqs, calc3},
mathematica comment (* Add suffix to name *)
mathematica code name = lookup[calc, Name] <> suffix;
mathematica code calc1 = mapReplace[calc, Name, name];
mathematica comment (* Replace some entries in the calculation *)
mathematica comment (* replaces = Map[Function[rule, mapReplace[#, rule[[1]], rule[[2]]]&], updates]; *)
mathematica code replaces = updates //. (lhs_ -> rhs_) -> (mapReplace[#, lhs, rhs]&);
mathematica code calc2 = Apply[Composition, replaces][calc1];
mathematica comment (* Remove unnecessary equations *)
mathematica code vars = Join[evolVars, lookup[calc2, Shorthands]];
mathematica code patterns = Replace[vars, { Tensor[n_,__] -> Tensor[n,__] ,
mathematica code dot[Tensor[n_,__]] -> dot[Tensor[n,__]]}, 1];
mathematica code eqs = FilterRules[lookup[calc, Equations], patterns];
mathematica code calc3 = mapReplace[calc2, Equations, eqs];
mathematica code calc3
mathematica code ];
mathematica blank
mathematica blank
mathematica blank
mathematica comment (******************************************************************************)
mathematica comment (* Convert from ADMBase *)
mathematica comment (******************************************************************************)
mathematica blank
mathematica code convertFromADMBaseCalc =
mathematica code {
mathematica code Name -> thorn <> "_convertFromADMBase",
mathematica code Schedule -> {"AT initial AFTER ADMBase_PostInitial"},
mathematica code ConditionalOnKeyword -> {"my_initial_data", "ADMBase"},
mathematica code Shorthands -> {g[la,lb], detg, gu[ua,ub], em4phi},
mathematica code Equations ->
mathematica code {
mathematica code g[la,lb] -> admg[la,lb],
mathematica code detg -> detgExpr,
mathematica code gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]],
mathematica blank
mathematica code phi -> IfThen[conformalMethod==CMW, detg^(-1/6), Log[detg]/12],
mathematica code em4phi -> IfThen[conformalMethod==CMW, phi^2, Exp[-4 phi]],
mathematica code gt[la,lb] -> em4phi g[la,lb],
mathematica blank
mathematica code trK -> gu[ua,ub] admK[la,lb],
mathematica code At[la,lb] -> em4phi (admK[la,lb] - (1/3) g[la,lb] trK),
mathematica blank
mathematica code alpha -> admalpha,
mathematica blank
mathematica code beta[ua] -> admbeta[ua],
mathematica blank
mathematica code IfCCZ4[Theta -> 0]
mathematica code }
mathematica code };
mathematica blank
mathematica code convertFromADMBaseGammaCalc =
mathematica code {
mathematica code Name -> thorn <> "_convertFromADMBaseGamma",
mathematica code Schedule -> {"AT initial AFTER " <> thorn <> "_convertFromADMBase"},
mathematica code ConditionalOnKeyword -> {"my_initial_data", "ADMBase"},
mathematica comment (*
mathematica comment Where -> InteriorNoSync,
mathematica comment *)
mathematica comment (* Do not synchronise right after this routine; instead, synchronise
mathematica comment after extrapolating *)
mathematica code Where -> Interior,
mathematica comment (* Synchronise after this routine, so that the refinement boundaries
mathematica comment are set correctly before extrapolating. (We will need to
mathematica comment synchronise again after extrapolating because extrapolation does
mathematica comment not fill ghost zones, but this is irrelevant here.) *)
mathematica code Shorthands -> {dir[ua],
mathematica code detgt, gtu[ua,ub], Gt[ua,lb,lc], theta},
mathematica code Equations ->
mathematica code {
mathematica code dir[ua] -> Sign[beta[ua]],
mathematica blank
mathematica code detgt -> 1 (* detgtExpr *),
mathematica code gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
mathematica code Gt[ua,lb,lc] -> 1/2 gtu[ua,ud]
mathematica code (PD[gt[lb,ld],lc] + PD[gt[lc,ld],lb] - PD[gt[lb,lc],ld]),
mathematica code Xt[ua] -> gtu[ub,uc] Gt[ua,lb,lc],
mathematica blank
mathematica comment (*
mathematica comment A -> - admdtalpha / (harmonicF alpha^harmonicN) (LapseAdvectionCoeff - 1),
mathematica comment *)
mathematica comment (* If LapseACoeff=0, then A is not evolved, in the sense that it
mathematica comment does not influence the time evolution of other variables. *)
mathematica code A -> IfThen [LapseACoeff != 0,
mathematica code 1 / (- harmonicF alpha^harmonicN)
mathematica code (+ admdtalpha
mathematica code - LapseAdvectionCoeff Upwind[beta[ua], alpha, la]),
mathematica code 0],
mathematica blank
mathematica code theta -> thetaExpr,
mathematica blank
mathematica comment (* If ShiftBCoeff=0 or theta ShiftGammaCoeff=0, then B^i is not
mathematica comment evolved, in the sense that it does not influence the time
mathematica comment evolution of other variables. *)
mathematica code B[ua] -> IfThen [ShiftGammaCoeff ShiftBCoeff != 0,
mathematica code 1 / (theta ShiftGammaCoeff)
mathematica code (+ admdtbeta[ua]
mathematica code - ShiftAdvectionCoeff Upwind[beta[ub], beta[ua], lb]),
mathematica code 0]
mathematica code }
mathematica code };
mathematica blank
mathematica comment (* Initialise the Gamma variables to 0. This is necessary with
mathematica comment multipatch because convertFromADMBaseGamma does not perform the
mathematica comment conversion in the boundary points, and the order in which symmetry
mathematica comment (interpatch) and outer boundary conditions is applied means that
mathematica comment points which are both interpatch and symmetry points are never
mathematica comment initialised. *)
mathematica code initGammaCalc =
mathematica code {
mathematica code Name -> thorn <> "_InitGamma",
mathematica code Schedule -> {"AT initial BEFORE " <> thorn <> "_convertFromADMBaseGamma"},
mathematica code ConditionalOnKeyword -> {"my_initial_data", "ADMBase"},
mathematica code Where -> Everywhere,
mathematica code Equations ->
mathematica code {
mathematica code Xt[ua] -> 0,
mathematica code A -> 0,
mathematica code B[ua] -> 0
mathematica code }
mathematica code };
mathematica blank
mathematica blank
mathematica blank
mathematica comment (******************************************************************************)
mathematica comment (* Convert to ADMBase *)
mathematica comment (******************************************************************************)
mathematica blank
mathematica code convertToADMBaseCalc =
mathematica code {
mathematica code Name -> thorn <> "_convertToADMBase",
mathematica code Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"},
mathematica code Where -> Everywhere,
mathematica code Shorthands -> {e4phi},
mathematica code Equations ->
mathematica code {
mathematica code e4phi -> IfThen[conformalMethod==CMW, 1/phi^2, Exp[4 phi]],
mathematica code admg[la,lb] -> e4phi gt[la,lb],
mathematica code admK[la,lb] -> e4phi At[la,lb] + (1/3) admg[la,lb] trK,
mathematica code admalpha -> alpha,
mathematica code admbeta[ua] -> beta[ua]
mathematica code }
mathematica code };
mathematica blank
mathematica code convertToADMBaseDtLapseShiftCalc =
mathematica code {
mathematica code Name -> thorn <> "_convertToADMBaseDtLapseShift",
mathematica code Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"},
mathematica code ConditionalOnKeyword -> {"dt_lapse_shift_method", "correct"},
mathematica code Where -> Interior,
mathematica code Shorthands -> {dir[ua], detgt, gtu[ua,ub], eta, theta},
mathematica code Equations ->
mathematica code {
mathematica code dir[ua] -> Sign[beta[ua]],
mathematica blank
mathematica code detgt -> 1 (* detgtExpr *),
mathematica comment (* This leads to simpler code... *)
mathematica code gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
mathematica blank
mathematica code eta -> etaExpr,
mathematica code theta -> thetaExpr,
mathematica blank
mathematica comment (* see RHS *)
mathematica comment (*
mathematica comment admdtalpha -> - harmonicF alpha^harmonicN
mathematica comment ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK)
mathematica comment + LapseAdvectionCoeff beta[ua] PDu[alpha,la],
mathematica comment *)
mathematica code admdtalpha -> - harmonicF alpha^harmonicN
mathematica code (+ LapseACoeff A
mathematica code + ((1 - LapseACoeff)
mathematica code (trK - IfCCZ4[2 Theta, 0])))
mathematica code + LapseAdvectionCoeff Upwind[beta[ua], alpha, la],
mathematica code admdtbeta[ua] -> IfThen[harmonicShift,
mathematica code - 1/2 gtu[ua,uj] phi alpha
mathematica code (- 2 alpha PD[phi,lj]
mathematica code + 2 phi PD[alpha,lj]
mathematica code + gtu[uk,ul] phi alpha
mathematica code (PD[gt[lk,ll],lj] - 2 PD[gt[lj,lk],ll])),
mathematica comment (* else *)
mathematica code + theta ShiftGammaCoeff
mathematica code (+ ShiftBCoeff B[ua]
mathematica code + (1 - ShiftBCoeff)
mathematica code (Xt[ua] - eta BetaDriver beta[ua]))]
mathematica code + ShiftAdvectionCoeff Upwind[beta[ub], beta[ua], lb]
mathematica code }
mathematica code };
mathematica blank
mathematica code convertToADMBaseDtLapseShiftBoundaryCalc =
mathematica code {
mathematica code Name -> thorn <> "_convertToADMBaseDtLapseShiftBoundary",
mathematica code Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"},
mathematica code ConditionalOnKeyword -> {"dt_lapse_shift_method", "correct"},
mathematica code Where -> BoundaryWithGhosts,
mathematica code Shorthands -> {detgt, gtu[ua,ub], eta, theta},
mathematica code Equations ->
mathematica code {
mathematica code detgt -> 1 (* detgtExpr *),
mathematica comment (* This leads to simpler code... *)
mathematica code gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
mathematica blank
mathematica code eta -> etaExpr,
mathematica code theta -> thetaExpr,
mathematica blank
mathematica comment (* see RHS, but omit derivatives near the boundary *)
mathematica comment (*
mathematica comment admdtalpha -> - harmonicF alpha^harmonicN
mathematica comment ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK),
mathematica comment *)
mathematica code admdtalpha -> - harmonicF alpha^harmonicN
mathematica code (+ LapseACoeff A
mathematica code + ((1 - LapseACoeff)
mathematica code (trK - IfCCZ4[2 Theta, 0]))),
mathematica code admdtbeta[ua] -> IfThen[harmonicShift,
mathematica code 0,
mathematica comment (* else *)
mathematica code + theta ShiftGammaCoeff
mathematica code (+ ShiftBCoeff B[ua]
mathematica code + (1 - ShiftBCoeff)
mathematica code (Xt[ua] - eta BetaDriver beta[ua]))]
mathematica code }
mathematica code };
mathematica blank
mathematica code convertToADMBaseFakeDtLapseShiftCalc =
mathematica code {
mathematica code Name -> thorn <> "_convertToADMBaseFakeDtLapseShift",
mathematica code Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"},
mathematica code ConditionalOnKeyword -> {"dt_lapse_shift_method", "noLapseShiftAdvection"},
mathematica code Where -> Everywhere,
mathematica code Shorthands -> {detgt, gtu[ua,ub], eta, theta},
mathematica code Equations ->
mathematica code {
mathematica code detgt -> 1 (* detgtExpr *),
mathematica comment (* This leads to simpler code... *)
mathematica code gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
mathematica blank
mathematica code eta -> etaExpr,
mathematica code theta -> thetaExpr,
mathematica blank
mathematica comment (* see RHS, but omit derivatives everywhere (which is wrong, but
mathematica comment faster, since it does not require synchronisation or boundary
mathematica comment conditions) *)
mathematica comment (*
mathematica comment admdtalpha -> - harmonicF alpha^harmonicN
mathematica comment ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK),
mathematica comment *)
mathematica code admdtalpha -> - harmonicF alpha^harmonicN
mathematica code (+ LapseACoeff A
mathematica code + ((1 - LapseACoeff)
mathematica code (trK - IfCCZ4[2 Theta, 0]))),
mathematica code admdtbeta[ua] -> IfThen[harmonicShift,
mathematica code 0,
mathematica comment (* else *)
mathematica code + theta ShiftGammaCoeff
mathematica code (+ ShiftBCoeff B[ua]
mathematica code + (1 - ShiftBCoeff)
mathematica code (Xt[ua] - eta BetaDriver beta[ua]))]
mathematica code }
mathematica code };
mathematica blank
mathematica comment (******************************************************************************)
mathematica comment (* Evolution equations *)
mathematica comment (******************************************************************************)
mathematica blank
mathematica code evolCalc =
mathematica code {
mathematica code Name -> thorn <> "_RHS",
mathematica code Schedule -> {"IN " <> thorn <> "_evolCalcGroup"},
mathematica comment (*
mathematica comment Where -> Interior,
mathematica comment *)
mathematica comment (* Synchronise the RHS grid functions after this routine, so that
mathematica comment the refinement boundaries are set correctly before applying the
mathematica comment radiative boundary conditions. *)
mathematica code Where -> InteriorNoSync,
mathematica code ConditionalOnKeyword -> {"RHS_split", "combined"},
mathematica code Shorthands -> {dir[ua],
mathematica code detgt, gtu[ua,ub],
mathematica code Gt[ua,lb,lc], Gtl[la,lb,lc], Gtlu[la,lb,uc], Xtn[ua],
mathematica code Rt[la,lb], Rphi[la,lb], R[la,lb],
mathematica code Atm[ua,lb], Atu[ua,ub],
mathematica code e4phi, em4phi, cdphi[la], cdphi2[la,lb], g[la,lb], detg,
mathematica code gu[ua,ub], Ats[la,lb], trAts, eta, theta,
mathematica code rho, S[la], trS, fac1, fac2, dottrK, dotXt[ua],
mathematica code epsdiss[ua], IfCCZ4[Z[ua]], IfCCZ4[dotTheta]},
mathematica code Equations ->
mathematica code {
mathematica code dir[ua] -> Sign[beta[ua]],
mathematica blank
mathematica code detgt -> 1 (* detgtExpr *),
mathematica blank
mathematica comment (* This leads to simpler code... *)
mathematica code gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
mathematica code Gtl[la,lb,lc] -> 1/2
mathematica code (PD[gt[lb,la],lc] + PD[gt[lc,la],lb] - PD[gt[lb,lc],la]),
mathematica code Gtlu[la,lb,uc] -> gtu[uc,ud] Gtl[la,lb,ld],
mathematica code Gt[ua,lb,lc] -> gtu[ua,ud] Gtl[ld,lb,lc],
mathematica blank
mathematica comment (* The conformal connection functions calculated from the conformal metric,
mathematica comment used instead of Xt where no derivatives of Xt are taken *)
mathematica code Xtn[ui] -> gtu[uj,uk] Gt[ui,lj,lk],
mathematica blank
mathematica code e4phi -> IfThen[conformalMethod==CMW, 1/phi^2, Exp[4 phi]],
mathematica code em4phi -> 1 / e4phi,
mathematica code g[la,lb] -> e4phi gt[la,lb],
mathematica code detg -> detgExpr,
mathematica code gu[ua,ub] -> em4phi gtu[ua,ub],
mathematica blank
mathematica comment (* The Z quantities *)
mathematica comment (* gr-qc:1106.2254 (2011), eqn. (23) *)
mathematica code IfCCZ4[
mathematica code Z[ud] -> (1/2) gu[ua,ud] (- PD[gt[la,lb],lc] gtu[ub,uc] + gt[la,lc] Xt[uc])
mathematica code ],
mathematica blank
mathematica comment (* PRD 62, 044034 (2000), eqn. (18) *)
mathematica comment (* Adding Z term by changing Xtn to Xt *)
mathematica code Rt[li,lj] -> - (1/2) gtu[ul,um] PD[gt[li,lj],ll,lm]
mathematica code + (1/2) gt[lk,li] PD[Xt[uk],lj]
mathematica code + (1/2) gt[lk,lj] PD[Xt[uk],li]
mathematica code + (1/2) Xtn[uk] Gtl[li,lj,lk]
mathematica code + (1/2) Xtn[uk] Gtl[lj,li,lk]
mathematica code + (+ Gt[uk,li,ll] Gtlu[lj,lk,ul]
mathematica code + Gt[uk,lj,ll] Gtlu[li,lk,ul]
mathematica code + Gt[uk,li,ll] Gtlu[lk,lj,ul]),
mathematica blank
mathematica code fac1 -> IfThen[conformalMethod==CMW, -1/(2 phi), 1],
mathematica code cdphi[la] -> fac1 CDt[phi,la],
mathematica code fac2 -> IfThen[conformalMethod==CMW, 1/(2 phi^2), 0],
mathematica code cdphi2[la,lb] -> fac1 CDt[phi,la,lb] + fac2 CDt[phi,la] CDt[phi,lb],
mathematica blank
mathematica comment (* PRD 62, 044034 (2000), eqn. (15) *)
mathematica code Rphi[li,lj] -> - 2 cdphi2[lj,li]
mathematica code - 2 gt[li,lj] gtu[ul,un] cdphi2[ll,ln]
mathematica code + 4 cdphi[li] cdphi[lj]
mathematica code - 4 gt[li,lj] gtu[ul,un] cdphi[ln] cdphi[ll],
mathematica blank
mathematica code Atm[ua,lb] -> gtu[ua,uc] At[lc,lb],
mathematica code Atu[ua,ub] -> Atm[ua,lc] gtu[ub,uc],
mathematica blank
mathematica code R[la,lb] -> Rt[la,lb] + Rphi[la,lb],
mathematica code IfCCZ4[
mathematica code R[la,lb] -> R[la,lb] + (2/phi) (+ g[la,lc] Z[uc] PD[phi,lb]
mathematica code + g[lb,lc] Z[uc] PD[phi,la] - g[la,lb] Z[uc] PD[phi,lc])
mathematica code + e4phi Z[uc] PD[gt[la,lb],lc]
mathematica code ],
mathematica blank
mathematica comment (* Matter terms *)
mathematica blank
mathematica comment (* rho = n^a n^b T_ab *)
mathematica code rho -> addMatter
mathematica code (1/alpha^2 (T00 - 2 beta[ui] T0[li] + beta[ui] beta[uj] T[li,lj])),
mathematica blank
mathematica comment (* S_i = -p^a_i n^b T_ab, where p^a_i = delta^a_i + n^a n_i *)
mathematica code S[li] -> addMatter (-1/alpha (T0[li] - beta[uj] T[li,lj])),
mathematica blank
mathematica comment (* trS = gamma^ij T_ij *)
mathematica code trS -> addMatter (em4phi gtu[ui,uj] T[li,lj]),
mathematica blank
mathematica comment (* RHS terms *)
mathematica blank
mathematica comment (* PRD 62, 044034 (2000), eqn. (10) *)
mathematica comment (* PRD 67 084023 (2003), eqn. (16) and (23) *)
mathematica code dot[phi] -> IfThen[conformalMethod==CMW, 1/3 phi, -1/6]
mathematica code (alpha trK - PD[beta[ua],la]),
mathematica blank
mathematica comment (* PRD 62, 044034 (2000), eqn. (9) *)
mathematica comment (* gr-qc:1106.2254 (2011), eqn. (14) *)
mathematica comment (* removing trA from Aij ensures that detg = 1 *)
mathematica code dot[gt[la,lb]] -> - 2 alpha (At[la,lb] - IfCCZ4[(1/3) At[lc,ld] gtu[uc,ud] gt[la,lb], 0])
mathematica code + gt[la,lc] PD[beta[uc],lb] + gt[lb,lc] PD[beta[uc],la]
mathematica code - (2/3) gt[la,lb] PD[beta[uc],lc],
mathematica comment (* PRD 62, 044034 (2000), eqn. (20) *)
mathematica comment (* PRD 67 084023 (2003), eqn (26) *)
mathematica comment (* gr-qc:1106.2254 (2011), eqn. (19) *)
mathematica comment (* Adding Z terms by changing Xtn to Xt,
mathematica comment also adding extra Z and Theta terms *)
mathematica code dotXt[ui] -> - 2 Atu[ui,uj] PD[alpha,lj]
mathematica code + 2 alpha (+ Gt[ui,lj,lk] Atu[uk,uj]
mathematica code - (2/3) gtu[ui,uj] PD[trK,lj]
mathematica code + 6 Atu[ui,uj] cdphi[lj])
mathematica code + gtu[uj,ul] PD[beta[ui],lj,ll]
mathematica code + (1/3) gtu[ui,uj] PD[beta[ul],lj,ll]
mathematica code - Xtn[uj] PD[beta[ui],lj]
mathematica code + (2/3) Xtn[ui] PD[beta[uj],lj]
mathematica code + IfCCZ4[
mathematica code + GammaShift 2 e4phi (- Z[uj] PD[beta[ui],lj]
mathematica code + (2/3) Z[ui] PD[beta[uj],lj])
mathematica code - (4/3) alpha e4phi Z[ui] trK
mathematica code + 2 gtu[ui,uj] (+ alpha PD[Theta,lj]
mathematica code - Theta PD[alpha,lj])
mathematica code - 2 alpha e4phi dampk1 Z[ui],
mathematica code 0]
mathematica comment (* Equation (4.28) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *)
mathematica code + addMatter (- 16 Pi alpha gtu[ui,uj] S[lj]),
mathematica code dot[Xt[ui]] -> dotXt[ui],
mathematica blank
mathematica comment (* gr-qc:1106.2254 (2011), eqn. (18) *)
mathematica code IfCCZ4[
mathematica code dotTheta ->
mathematica code - PD[alpha,la] Z[ua] - dampk1 (2 + dampk2) alpha Theta
mathematica code + (1/2) alpha (gu[ua,ub] R[la,lb] - Atm[ua,lb] Atm[ub,la] + (2/3) trK^2 - 2 trK Theta)
mathematica code + addMatter (- 8 Pi alpha rho)
mathematica code ],
mathematica blank
mathematica code IfCCZ4[
mathematica code dot[Theta] -> dotTheta
mathematica code ],
mathematica blank
mathematica comment (* PRD 62, 044034 (2000), eqn. (11) *)
mathematica comment (* gr-qc:1106.2254 (2011), eqn. (17) *)
mathematica comment (* Adding the RHS of Theta to K, because K_Z4 = K_BSSN + 2 Theta *)
mathematica comment (* Also adding the Z term, as it has to cancel with the one in Theta *)
mathematica code dottrK -> - em4phi ( gtu[ua,ub] ( PD[alpha,la,lb]
mathematica code + 2 cdphi[la] PD[alpha,lb] )
mathematica code - Xtn[ua] PD[alpha,la] )
mathematica code + alpha (Atm[ua,lb] Atm[ub,la] + (1/3) trK^2)
mathematica code + IfCCZ4[
mathematica code + 2 dotTheta + 2 PD[alpha,la] Z[ua]
mathematica code + dampk1 (1 - dampk2) alpha Theta,
mathematica code 0]
mathematica comment (* Equation (4.21) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *)
mathematica code + addMatter (4 Pi alpha (rho + trS)),
mathematica code dot[trK] -> dottrK,
mathematica blank
mathematica comment (* PRD 62, 044034 (2000), eqn. (12) *)
mathematica comment (* TODO: Should we use the Hamiltonian constraint to make Rij tracefree? *)
mathematica comment (* gr-qc:1106.2254 (2011), eqn. (15) *)
mathematica comment (* Adding Z terms in the Ricci and Theta terms *)
mathematica code Ats[la,lb] -> - CDt[alpha,la,lb] +
mathematica code + 2 (PD[alpha,la] cdphi[lb] + PD[alpha,lb] cdphi[la] )
mathematica code + alpha R[la,lb],
mathematica code trAts -> gu[ua,ub] Ats[la,lb],
mathematica code dot[At[la,lb]] -> + em4phi (+ Ats[la,lb] - (1/3) g[la,lb] trAts )
mathematica code + alpha (+ ((trK - IfCCZ4[2 Theta, 0])
mathematica code At[la,lb])
mathematica code - 2 At[la,lc] Atm[uc,lb])
mathematica code + At[la,lc] PD[beta[uc],lb] + At[lb,lc] PD[beta[uc],la]
mathematica code - (2/3) At[la,lb] PD[beta[uc],lc]
mathematica comment (* Equation (4.23) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *)
mathematica code + addMatter (- em4phi alpha 8 Pi
mathematica code (T[la,lb] - (1/3) g[la,lb] trS)),
mathematica blank
mathematica comment (* dot[alpha] -> - harmonicF alpha^harmonicN trK, *)
mathematica comment (* dot[alpha] -> - harmonicF alpha^harmonicN A + Lie[alpha, beta], *)
mathematica comment (*
mathematica comment dot[alpha] -> - harmonicF alpha^harmonicN (
mathematica comment (1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK)
mathematica comment + LapseAdvectionCoeff beta[ua] PDu[alpha,la],
mathematica comment dot[A] -> (1 - LapseAdvectionCoeff) (dottrK - AlphaDriver A),
mathematica comment *)
mathematica code dot[alpha] -> - harmonicF alpha^harmonicN
mathematica code (+ LapseACoeff A
mathematica code + ((1 - LapseACoeff)
mathematica code (+ trK - IfCCZ4[2 Theta, 0]
mathematica code + AlphaDriver (alpha - 1)))),
mathematica blank
mathematica code dot[A] -> + (LapseACoeff
mathematica code (+ dottrK - IfCCZ4[2 dotTheta, 0]
mathematica code - AlphaDriver A)),
mathematica blank
mathematica code eta -> etaExpr,
mathematica code theta -> thetaExpr,
mathematica blank
mathematica comment (* dot[beta[ua]] -> eta Xt[ua], *)
mathematica comment (* dot[beta[ua]] -> ShiftGammaCoeff alpha^ShiftAlphaPower B[ua], *)
mathematica code dot[beta[ua]] -> IfThen[harmonicShift,
mathematica code - 1/2 gtu[ua,uj] phi alpha
mathematica code (- 2 alpha PD[phi,lj]
mathematica code + 2 phi PD[alpha,lj]
mathematica code + gtu[uk,ul] phi alpha
mathematica code (PD[gt[lk,ll],lj] - 2 PD[gt[lj,lk],ll])),
mathematica comment (* else *)
mathematica code + theta ShiftGammaCoeff
mathematica code (+ ShiftBCoeff B[ua]
mathematica code + (1 - ShiftBCoeff)
mathematica code (Xt[ua] - eta BetaDriver beta[ua]))],
mathematica blank
mathematica code dot[B[ua]] -> + ShiftBCoeff (dotXt[ua] - eta BetaDriver B[ua])
mathematica comment (* Note that this dotXt[ua] is not yet \partial_t \tilde \Gamma^i, because the
mathematica comment advection term has not yet been added. It is actually
mathematica comment \partial_t \tilde \Gamma^i - \beta^j \partial_j \tilde \Gamma^i *)
mathematica code }
mathematica code };
mathematica blank
mathematica code advectCalc =
mathematica code {
mathematica code Name -> thorn <> "_Advect",
mathematica code Schedule -> {"IN " <> thorn <> "_evolCalcGroup " <>
mathematica code "AFTER (" <> thorn <> "_RHS " <> thorn <> "_RHS1 " <> thorn <> "_RHS2)"},
mathematica comment (*
mathematica comment Where -> Interior,
mathematica comment *)
mathematica comment (* Synchronise the RHS grid functions after this routine, so that
mathematica comment the refinement boundaries are set correctly before applying the
mathematica comment radiative boundary conditions. *)
mathematica code Where -> InteriorNoSync,
mathematica code ConditionalOnKeyword -> {"advection_split", "combined"},
mathematica code Shorthands -> {dir[ua]},
mathematica code Equations ->
mathematica code {
mathematica code dir[ua] -> Sign[beta[ua]],
mathematica blank
mathematica code dot[phi] -> dot[phi] + Upwind[beta[ua], phi, la],
mathematica blank
mathematica code dot[gt[la,lb]] -> dot[gt[la,lb]] + Upwind[beta[uc], gt[la,lb], lc],
mathematica blank
mathematica code dot[Xt[ui]] -> dot[Xt[ui]] + Upwind[beta[uj], Xt[ui], lj],
mathematica blank
mathematica code IfCCZ4[
mathematica code dot[Theta] -> dot[Theta] + Upwind[beta[ua], Theta, la]
mathematica code ],
mathematica blank
mathematica code dot[trK] -> dot[trK] + Upwind[beta[ua], trK, la],
mathematica blank
mathematica code dot[At[la,lb]] -> dot[At[la,lb]] + Upwind[beta[uc], At[la,lb], lc],
mathematica blank
mathematica code dot[alpha] -> dot[alpha]
mathematica code + LapseAdvectionCoeff Upwind[beta[ua], alpha, la],
mathematica blank
mathematica code dot[A] -> dot[A]
mathematica code + LapseACoeff (
mathematica code + LapseAdvectionCoeff Upwind[beta[ua], A, la]
mathematica code + (1 - LapseAdvectionCoeff) Upwind[beta[ua], trK, la]),
mathematica blank
mathematica code dot[beta[ua]] -> dot[beta[ua]]
mathematica code + ShiftAdvectionCoeff Upwind[beta[ub], beta[ua], lb],
mathematica blank
mathematica code dot[B[ua]] -> dot[B[ua]]
mathematica code + ShiftBCoeff (
mathematica code + ShiftAdvectionCoeff Upwind[beta[ub], B[ua], lb]
mathematica code + ((1 - ShiftAdvectionCoeff)
mathematica code Upwind[beta[ub], Xt[ua], lb]))
mathematica comment (* Note that the advection term \beta^j \partial_j \tilde \Gamma^i is not
mathematica comment subtracted here when ShiftAdvectionCoefficient == 1 because it was
mathematica comment implicitly subtracted before (see comment in previous calculation of
mathematica comment dot[B[ua]]. *)
mathematica code }
mathematica code };
mathematica blank
mathematica code varsNames = {
mathematica code {"phi", dot[phi]},
mathematica code {"gt", dot[gt[la,lb]]},
mathematica code {"Xt", dot[Xt[ui]]},
mathematica code {"trK", dot[trK]},
mathematica code {"At", dot[At[la,lb]]},
mathematica code {"alpha", dot[alpha]},
mathematica code {"A", dot[A]},
mathematica code {"beta", dot[beta[ua]]},
mathematica code {"B", dot[B[ua]]},
mathematica code IfCCZ4[{"Theta", dot[Theta]}]
mathematica code };
mathematica blank
mathematica code advectCalcs = Map[
mathematica code PartialCalculation[advectCalc, "_"<>ToString[First[#]],
mathematica code {ConditionalOnKeyword -> {"advection_split",
mathematica code "per variable"}},
mathematica code {Last[#]}]&,
mathematica code varsNames];
mathematica blank
mathematica code evolCalc1 = PartialCalculation[evolCalc, "1",
mathematica code {
mathematica code ConditionalOnKeyword -> {"RHS_split", "split At"}
mathematica code },
mathematica code {
mathematica code dot[phi],
mathematica code dot[gt[la,lb]],
mathematica code dot[Xt[ui]],
mathematica code dot[trK],
mathematica code dot[alpha],
mathematica code dot[A],
mathematica code dot[beta[ua]],
mathematica code dot[B[ua]],
mathematica code IfCCZ4[dot[Theta]]
mathematica code }];
mathematica blank
mathematica code evolCalc2 = PartialCalculation[evolCalc, "2",
mathematica code {
mathematica code ConditionalOnKeyword -> {"RHS_split", "split At"}
mathematica code },
mathematica code {
mathematica code dot[At[la,lb]]
mathematica code }];
mathematica blank
mathematica code dissCalc =
mathematica code {
mathematica code Name -> thorn <> "_Dissipation",
mathematica code Schedule -> {"IN " <> thorn <> "_evolCalcGroup " <>
mathematica code "AFTER (" <> thorn <> "_RHS1 " <> thorn <> "_RHS2)"},
mathematica code ConditionalOnKeyword -> {"apply_dissipation", "always"},
mathematica code Where -> InteriorNoSync,
mathematica code Shorthands -> {epsdiss[ua]},
mathematica code Equations ->
mathematica code {
mathematica code epsdiss[ua] -> EpsDiss,
mathematica code Sequence@@Table[
mathematica code dot[var] -> dot[var] + epsdiss[ux] PDdiss[var,lx],
mathematica code {var, {phi, gt[la,lb], Xt[ui], IfCCZ4[Theta], trK, At[la,lb],
mathematica code alpha, A, beta[ua], B[ua]}}]
mathematica code }
mathematica code };
mathematica blank
mathematica code dissCalcs =
mathematica code Table[
mathematica code {
mathematica code Name -> thorn <> "_Dissipation_" <> ToString[var /. {Tensor[n_,__] -> n}],
mathematica code Schedule -> {"IN " <> thorn <> "_evolCalcGroup " <>
mathematica code "AFTER (" <> thorn <> "_RHS1 " <> thorn <> "_RHS2)"},
mathematica code ConditionalOnKeyword -> {"apply_dissipation", "always"},
mathematica code Where -> InteriorNoSync,
mathematica code Shorthands -> {epsdiss[ua]},
mathematica code Equations ->
mathematica code {
mathematica code epsdiss[ua] -> EpsDiss,
mathematica code dot[var] -> dot[var] + epsdiss[ux] PDdiss[var,lx]
mathematica code }
mathematica code },
mathematica code {var, {phi, gt[la,lb], Xt[ui], IfCCZ4[Theta], trK, At[la,lb],
mathematica code alpha, A, beta[ua], B[ua]}}
mathematica code ];
mathematica blank
mathematica code RHSStaticBoundaryCalc =
mathematica code {
mathematica code Name -> thorn <> "_RHSStaticBoundary",
mathematica code Schedule -> {"IN MoL_CalcRHS"},
mathematica code ConditionalOnKeyword -> {"my_rhs_boundary_condition", "static"},
mathematica code Where -> Boundary,
mathematica code Equations ->
mathematica code {
mathematica code dot[phi] -> 0,
mathematica code dot[gt[la,lb]] -> 0,
mathematica code dot[trK] -> 0,
mathematica code dot[At[la,lb]] -> 0,
mathematica code dot[Xt[ua]] -> 0,
mathematica code dot[alpha] -> 0,
mathematica code dot[A] -> 0,
mathematica code dot[beta[ua]] -> 0,
mathematica code dot[B[ua]] -> 0,
mathematica code IfCCZ4[dot[Theta] -> 0]
mathematica code }
mathematica code };
mathematica blank
mathematica comment (* Initialise the RHS variables in analysis in case they are going to
mathematica comment be output - the noninterior points cannot be filled, so we define
mathematica comment them to be zero *)
mathematica code initRHSCalc =
mathematica code {
mathematica code Name -> thorn <> "_InitRHS",
mathematica code Schedule -> {"AT analysis BEFORE " <> thorn <> "_evolCalcGroup"},
mathematica code Where -> Everywhere,
mathematica code Equations ->
mathematica code {
mathematica code dot[phi] -> 0,
mathematica code dot[gt[la,lb]] -> 0,
mathematica code dot[trK] -> 0,
mathematica code dot[At[la,lb]] -> 0,
mathematica code dot[Xt[ua]] -> 0,
mathematica code dot[alpha] -> 0,
mathematica code dot[A] -> 0,
mathematica code dot[beta[ua]] -> 0,
mathematica code dot[B[ua]] -> 0,
mathematica code IfCCZ4[dot[Theta] -> 0]
mathematica code }
mathematica code };
mathematica blank
mathematica code RHSRadiativeBoundaryCalc =
mathematica code {
mathematica code Name -> thorn <> "_RHSRadiativeBoundary",
mathematica code Schedule -> {"IN MoL_CalcRHS"},
mathematica code ConditionalOnKeyword -> {"my_rhs_boundary_condition", "radiative"},
mathematica code Where -> Boundary,
mathematica code Shorthands -> {dir[ua],
mathematica code detgt, gtu[ua,ub], em4phi, gu[ua,ub],
mathematica code nn[la], nu[ua], nlen, nlen2, su[ua],
mathematica code vg},
mathematica code Equations ->
mathematica code {
mathematica code dir[ua] -> Sign[normal[ua]],
mathematica blank
mathematica code detgt -> 1 (* detgtExpr *),
mathematica code gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
mathematica code em4phi -> IfThen[conformalMethod==CMW, phi^2, Exp[-4 phi]],
mathematica code gu[ua,ub] -> em4phi gtu[ua,ub],
mathematica blank
mathematica code nn[la] -> Euc[la,lb] normal[ub],
mathematica code nu[ua] -> gu[ua,ub] nn[lb],
mathematica code nlen2 -> nu[ua] nn[la],
mathematica code nlen -> Sqrt[nlen2],
mathematica code su[ua] -> nu[ua] / nlen,
mathematica blank
mathematica code vg -> Sqrt[harmonicF],
mathematica blank
mathematica code dot[phi] -> - vg su[uc] PDo[phi ,lc],
mathematica code dot[gt[la,lb]] -> - su[uc] PDo[gt[la,lb],lc],
mathematica code dot[trK] -> - vg su[uc] PDo[trK ,lc],
mathematica code dot[At[la,lb]] -> - su[uc] PDo[At[la,lb],lc],
mathematica code dot[Xt[ua]] -> - su[uc] PDo[Xt[ua] ,lc],
mathematica code dot[alpha] -> - vg su[uc] PDo[alpha ,lc],
mathematica code dot[A] -> - vg su[uc] PDo[A ,lc],
mathematica code dot[beta[ua]] -> - su[uc] PDo[beta[ua] ,lc],
mathematica code dot[B[ua]] -> - su[uc] PDo[B[ua] ,lc],
mathematica code IfCCZ4[
mathematica code dot[Theta] -> - vg su[uc] PDo[Theta ,lc]
mathematica code ]
mathematica code }
mathematica code };
mathematica blank
mathematica code enforceCalc =
mathematica code {
mathematica code Name -> thorn <> "_enforce",
mathematica code Schedule -> {"IN MoL_PostStepModify"},
mathematica code Shorthands -> {detgt, gtu[ua,ub], trAt},
mathematica code Equations ->
mathematica code {
mathematica comment (* The following comment is still interesting, but is not correct
mathematica comment any more since it is now scheduled in MoL_PostStepModify instead:
mathematica blank
mathematica comment Enforcing the constraints needs to be a projection, because it
mathematica comment is applied in MoL_PostStep and may thus be applied multiple
mathematica comment times, not only during time evolution. Therefore detgt has to
mathematica comment be calculated correctly, without assuming that det gt_ij = 1,
mathematica comment which is not always the case (since we don't enforce it). On
mathematica comment the other hand, this may not be so important... *)
mathematica code detgt -> 1 (* detgtExpr *),
mathematica code gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
mathematica blank
mathematica code trAt -> gtu[ua,ub] At[la,lb],
mathematica blank
mathematica code At[la,lb] -> At[la,lb] - (1/3) gt[la,lb] trAt,
mathematica blank
mathematica code alpha -> Max[alpha, MinimumLapse]
mathematica code }
mathematica code };
mathematica blank
mathematica comment (******************************************************************************)
mathematica comment (* Boundary conditions *)
mathematica comment (******************************************************************************)
mathematica blank
mathematica code boundaryCalc =
mathematica code {
mathematica code Name -> thorn <> "_boundary",
mathematica code Schedule -> {"IN MoL_PostStep"},
mathematica code ConditionalOnKeyword -> {"my_boundary_condition", "Minkowski"},
mathematica code Where -> BoundaryWithGhosts,
mathematica code Equations ->
mathematica code {
mathematica code phi -> IfThen[conformalMethod==CMW, 1, 0],
mathematica code gt[la,lb] -> KD[la,lb],
mathematica code trK -> 0,
mathematica code At[la,lb] -> 0,
mathematica code Xt[ua] -> 0,
mathematica code alpha -> 1,
mathematica code A -> 0,
mathematica code beta[ua] -> 0,
mathematica code B[ua] -> 0,
mathematica code IfCCZ4[Theta -> 0]
mathematica code }
mathematica code };
mathematica blank
mathematica comment (******************************************************************************)
mathematica comment (* Constraint equations *)
mathematica comment (******************************************************************************)
mathematica blank
mathematica code constraintsCalc =
mathematica code {
mathematica code Name -> thorn <> "_constraints",
mathematica code Schedule -> Automatic,
mathematica code After -> "MoL_PostStep",
mathematica code Where -> Interior,
mathematica code Shorthands -> {detgt, ddetgt[la], gtu[ua,ub], Z[ua],
mathematica code Gt[ua,lb,lc], Gtl[la,lb,lc], Gtlu[la,lb,uc], Xtn[ua],
mathematica code e4phi, em4phi,
mathematica code g[la,lb], detg, gu[ua,ub], ddetg[la], G[ua,lb,lc],
mathematica code Rt[la,lb], Rphi[la,lb], R[la,lb], trR, Atm[ua,lb],
mathematica code gK[la,lb,lc], cdphi[la], cdphi2[la,lb],
mathematica code rho, S[la], fac1, fac2},
mathematica code Equations ->
mathematica code {
mathematica code detgt -> 1 (* detgtExpr *),
mathematica code ddetgt[la] -> 0 (* ddetgtExpr[la] *),
mathematica blank
mathematica comment (* This leads to simpler code... *)
mathematica code gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
mathematica code Gtl[la,lb,lc] -> 1/2
mathematica code (PD[gt[lb,la],lc] + PD[gt[lc,la],lb] - PD[gt[lb,lc],la]),
mathematica code Gtlu[la,lb,uc] -> gtu[uc,ud] Gtl[la,lb,ld],
mathematica code Gt[ua,lb,lc] -> gtu[ua,ud] Gtl[ld,lb,lc],
mathematica blank
mathematica comment (* The conformal connection functions calculated from the conformal metric,
mathematica comment used instead of Xt where no derivatives of Xt are taken *)
mathematica code Xtn[ui] -> gtu[uj,uk] Gt[ui,lj,lk],
mathematica blank
mathematica code e4phi -> IfThen[conformalMethod==CMW, 1/phi^2, Exp[4 phi]],
mathematica code em4phi -> 1 / e4phi,
mathematica code g[la,lb] -> e4phi gt[la,lb],
mathematica code detg -> e4phi^3,
mathematica code gu[ua,ub] -> em4phi gtu[ua,ub],
mathematica blank
mathematica comment (* The Z quantities *)
mathematica code IfCCZ4[
mathematica code Z[ud] -> (1/2) gu[ua,ud] (- PD[gt[la,lb],lc] gtu[ub,uc] + gt[la,lc] Xt[uc])
mathematica code ],
mathematica blank
mathematica comment (* PRD 62, 044034 (2000), eqn. (18) *)
mathematica code Rt[li,lj] -> - (1/2) gtu[ul,um] PD[gt[li,lj],ll,lm]
mathematica code + (1/2) gt[lk,li] PD[Xt[uk],lj]
mathematica code + (1/2) gt[lk,lj] PD[Xt[uk],li]
mathematica code + (1/2) Xtn[uk] Gtl[li,lj,lk]
mathematica code + (1/2) Xtn[uk] Gtl[lj,li,lk]
mathematica code + (+ Gt[uk,li,ll] Gtlu[lj,lk,ul]
mathematica code + Gt[uk,lj,ll] Gtlu[li,lk,ul]
mathematica code + Gt[uk,li,ll] Gtlu[lk,lj,ul]),
mathematica blank
mathematica comment (* From the long turducken paper.
mathematica comment This expression seems to give the same result as the one from 044034. *)
mathematica comment (* TODO: symmetrise correctly: (ij) = (1/2) [i+j] *)
mathematica comment (*
mathematica comment Rt[li,lj] -> - (1/2) gtu[uk,ul] PD[gt[li,lj],lk,ll]
mathematica comment + gt[lk,li] PD[Xt[uk],lj] + gt[lk,lj] PD[Xt[uk],li]
mathematica comment + gt[li,ln] Gt[un,lj,lk] gtu[um,ua] gtu[uk,ub] PD[gt[la,lb],lm]
mathematica comment + gt[lj,ln] Gt[un,li,lk] gtu[um,ua] gtu[uk,ub] PD[gt[la,lb],lm]
mathematica comment + gtu[ul,us] (+ 2 Gt[uk,ll,li] gt[lj,ln] Gt[un,lk,ls]
mathematica comment + 2 Gt[uk,ll,lj] gt[li,ln] Gt[un,lk,ls]
mathematica comment + Gt[uk,li,ls] gt[lk,ln] Gt[un,ll,lj]),
mathematica comment *)
mathematica blank
mathematica comment (* Below would be a straightforward calculation,
mathematica comment without taking any Gamma^i into account.
mathematica comment This expression gives a different answer! *)
mathematica comment (*
mathematica comment Rt[la,lb] -> + Gt[u1,l2,la] Gt[l1,lb,u2] - Gt[u1,la,lb] Gt[l1,l2,u2]
mathematica comment + 1/2 gtu[u1,u2] (- PD[gt[l1,l2],la,lb] + PD[gt[l1,la],l2,lb]
mathematica comment - PD[gt[la,lb],l1,l2] + PD[gt[l2,lb],l1,la]),
mathematica comment *)
mathematica blank
mathematica code fac1 -> IfThen[conformalMethod==CMW, -1/(2 phi), 1],
mathematica code cdphi[la] -> fac1 CDt[phi,la],
mathematica code fac2 -> IfThen[conformalMethod==CMW, 1/(2 phi^2), 0],
mathematica code cdphi2[la,lb] -> fac1 CDt[phi,la,lb] + fac2 CDt[phi,la] CDt[phi,lb],
mathematica blank
mathematica comment (* PRD 62, 044034 (2000), eqn. (15) *)
mathematica code Rphi[li,lj] -> - 2 cdphi2[lj,li]
mathematica code - 2 gt[li,lj] gtu[ul,un] cdphi2[ll,ln]
mathematica code + 4 cdphi[li] cdphi[lj]
mathematica code - 4 gt[li,lj] gtu[ul,un] cdphi[ln] cdphi[ll],
mathematica blank
mathematica comment (* ddetg[la] -> PD[e4phi detg,la], *)
mathematica code ddetg[la] -> e4phi ddetgt[la] + 4 detgt e4phi PD[phi,la],
mathematica comment (* TODO: check this equation, maybe simplify it by omitting ddetg *)
mathematica code G[ua,lb,lc] -> Gt[ua,lb,lc]
mathematica code + 1/(2 detg) (+ KD[ua,lb] ddetg[lc] + KD[ua,lc] ddetg[lb]
mathematica code - (1/3) g[lb,lc] gu[ua,ud] ddetg[ld]),
mathematica blank
mathematica code R[la,lb] -> + Rt[la,lb] + Rphi[la,lb],
mathematica blank
mathematica code IfCCZ4[
mathematica code R[la,lb] -> R[la, lb] + (2/phi) (+ g[la,lc] Z[uc] PD[phi,lb]
mathematica code + g[lb,lc] Z[uc] PD[phi,la] - g[la,lb] Z[uc] PD[phi,lc])
mathematica code + e4phi Z[uc] PD[gt[la,lb],lc]
mathematica code ],
mathematica blank
mathematica code trR -> gu[ua,ub] R[la,lb],
mathematica blank
mathematica comment (* K[la,lb] -> e4phi At[la,lb] + (1/3) g[la,lb] trK, *)
mathematica comment (* Km[ua,lb] -> gu[ua,uc] K[lc,lb], *)
mathematica code Atm[ua,lb] -> gtu[ua,uc] At[lc,lb],
mathematica blank
mathematica comment (* Matter terms *)
mathematica blank
mathematica comment (* rho = n^a n^b T_ab *)
mathematica code rho -> 1/alpha^2 (T00 - 2 beta[ui] T0[li] + beta[ui] beta[uj] T[li,lj]),
mathematica blank
mathematica comment (* S_i = -p^a_i n^b T_ab, where p^a_i = delta^a_i + n^a n_i *)
mathematica code S[li] -> -1/alpha (T0[li] - beta[uj] T[li,lj]),
mathematica blank
mathematica comment (* Constraints *)
mathematica blank
mathematica comment (* H -> trR - Km[ua,lb] Km[ub,la] + trK^2, *)
mathematica comment (* PRD 67, 084023 (2003), eqn. (19) *)
mathematica code H -> trR - Atm[ua,lb] Atm[ub,la] + (2/3) trK^2 - addMatter 16 Pi rho,
mathematica blank
mathematica comment (* gK[la,lb,lc] -> CD[K[la,lb],lc], *)
mathematica comment (* gK[la,lb,lc] -> + 4 e4phi PD[phi,lc] At[la,lb] + e4phi CD[At[la,lb],lc]
mathematica comment + (1/3) g[la,lb] PD[trK,lc],
mathematica comment M[la] -> gu[ub,uc] (gK[lc,la,lb] - gK[lc,lb,la]), *)
mathematica blank
mathematica code M[li] -> + gtu[uj,uk] (CDt[At[li,lj],lk] + 6 At[li,lj] cdphi[lk])
mathematica code - (2/3) PD[trK,li]
mathematica code - addMatter 8 Pi S[li],
mathematica comment (* TODO: use PRD 67, 084023 (2003), eqn. (20) *)
mathematica blank
mathematica comment (* det gamma-tilde *)
mathematica code cS -> Log[detgt],
mathematica blank
mathematica comment (* Gamma constraint *)
mathematica code cXt[ua] -> gtu[ub,uc] Gt[ua,lb,lc] - Xt[ua],
mathematica blank
mathematica comment (* trace A-tilde *)
mathematica code cA -> gtu[ua,ub] At[la,lb]
mathematica code }
mathematica code };
mathematica blank
mathematica code constraintsCalc1 = PartialCalculation[constraintsCalc, "1",
mathematica code {},
mathematica code {
mathematica code H
mathematica code }];
mathematica blank
mathematica code constraintsCalc2 = PartialCalculation[constraintsCalc, "2",
mathematica code {},
mathematica code {
mathematica code M[li],
mathematica code cS,
mathematica code cXt[ua],
mathematica code cA
mathematica code }];
mathematica blank
mathematica comment (******************************************************************************)
mathematica comment (* Implementations *)
mathematica comment (******************************************************************************)
mathematica blank
mathematica code inheritedImplementations =
mathematica code Join[{"ADMBase"},
mathematica code If [addMatter!=0, {"TmunuBase"}, {}]];
mathematica blank
mathematica comment (******************************************************************************)
mathematica comment (* Parameters *)
mathematica comment (******************************************************************************)
mathematica blank
mathematica code inheritedKeywordParameters = {};
mathematica blank
mathematica code extendedKeywordParameters =
mathematica code {
mathematica code {
mathematica code Name -> "ADMBase::evolution_method",
mathematica code AllowedValues -> {thorn}
mathematica code },
mathematica code {
mathematica code Name -> "ADMBase::lapse_evolution_method",
mathematica code AllowedValues -> {thorn}
mathematica code },
mathematica code {
mathematica code Name -> "ADMBase::shift_evolution_method",
mathematica code AllowedValues -> {thorn}
mathematica code },
mathematica code {
mathematica code Name -> "ADMBase::dtlapse_evolution_method",
mathematica code AllowedValues -> {thorn}
mathematica code },
mathematica code {
mathematica code Name -> "ADMBase::dtshift_evolution_method",
mathematica code AllowedValues -> {thorn}
mathematica code }
mathematica code };
mathematica blank
mathematica code keywordParameters =
mathematica code {
mathematica code {
mathematica code Name -> "my_initial_data",
mathematica comment (* Visibility -> "restricted", *)
mathematica comment (* Description -> "ddd", *)
mathematica code AllowedValues -> {"ADMBase", "Minkowski"},
mathematica code Default -> "ADMBase"
mathematica code },
mathematica code {
mathematica code Name -> "my_initial_boundary_condition",
mathematica code Visibility -> "restricted",
mathematica comment (* Description -> "ddd", *)
mathematica code AllowedValues -> {"none"},
mathematica code Default -> "none"
mathematica code },
mathematica code {
mathematica code Name -> "my_rhs_boundary_condition",
mathematica code Visibility -> "restricted",
mathematica comment (* Description -> "ddd", *)
mathematica code AllowedValues -> {"none", "static", "radiative"},
mathematica code Default -> "none"
mathematica code },
mathematica code {
mathematica code Name -> "my_boundary_condition",
mathematica comment (* Visibility -> "restricted", *)
mathematica comment (* Description -> "ddd", *)
mathematica code AllowedValues -> {"none", "Minkowski"},
mathematica code Default -> "none"
mathematica code },
mathematica code {
mathematica code Name -> "calculate_ADMBase_variables_at",
mathematica code Visibility -> "restricted",
mathematica comment (* Description -> "ddd", *)
mathematica code AllowedValues -> {"MoL_PostStep", "CCTK_EVOL", "CCTK_ANALYSIS"},
mathematica code Default -> "MoL_PostStep"
mathematica code },
mathematica code {
mathematica code Name -> "UseSpatialBetaDriver_UNUSED",
mathematica code Visibility -> "restricted",
mathematica comment (* Description -> "ddd", *)
mathematica code AllowedValues -> {"no", "yes"},
mathematica code Default -> "no"
mathematica code },
mathematica code {
mathematica code Name -> "dt_lapse_shift_method",
mathematica code Description -> "Treatment of ADMBase dtlapse and dtshift",
mathematica code AllowedValues -> {"correct",
mathematica code "noLapseShiftAdvection" (* omit lapse and shift advection terms (faster) *)
mathematica code },
mathematica code Default -> "correct"
mathematica code },
mathematica code {
mathematica code Name -> "apply_dissipation",
mathematica code Description -> "Whether to apply dissipation to the RHSs",
mathematica code AllowedValues -> {"always",
mathematica code "never" (* yes and no keyword values confuse Cactus, and Kranc
mathematica comment doesn't support boolean parameters *)
mathematica code },
mathematica code Default -> "never"
mathematica code },
mathematica code {
mathematica code Name -> "RHS_split",
mathematica code Description -> "How to split RHS calculation",
mathematica code AllowedValues -> {"combined",
mathematica code "split At"},
mathematica code Default -> "split At"
mathematica code },
mathematica code {
mathematica code Name -> "advection_split",
mathematica code Description -> "How to split advection calculation",
mathematica code AllowedValues -> {"combined",
mathematica code "per variable"},
mathematica code Default -> "combined"
mathematica code }
mathematica code };
mathematica blank
mathematica code intParameters =
mathematica code {
mathematica code {
mathematica code Name -> harmonicN,
mathematica code Description -> "d/dt alpha = - f alpha^n K (harmonic=2, 1+log=1)",
mathematica code Default -> 2
mathematica code },
mathematica code {
mathematica code Name -> ShiftAlphaPower,
mathematica code Default -> 0
mathematica code },
mathematica code {
mathematica code Name -> conformalMethod,
mathematica code Description -> "Treatment of conformal factor",
mathematica code AllowedValues -> {{Value -> "0", Description -> "phi method"},
mathematica code {Value -> "1", Description -> "W method"}},
mathematica code Default -> 0
mathematica code },
mathematica code {
mathematica code Name -> fdOrder,
mathematica code Default -> derivOrder,
mathematica code AllowedValues -> {2,4,6,8}
mathematica code },
mathematica code {
mathematica code Name -> harmonicShift,
mathematica code Description -> "Whether to use the harmonic shift",
mathematica code AllowedValues -> {{Value -> "0", Description -> "Gamma driver shift"},
mathematica code {Value -> "1", Description -> "Harmonic shift"}},
mathematica code Default -> 0
mathematica code }
mathematica code };
mathematica blank
mathematica code realParameters =
mathematica code {
mathematica code IfCCZ4[{
mathematica code Name -> GammaShift,
mathematica code Description -> "Covariant shift term in Gamma",
mathematica code Default -> 0.5
mathematica code }],
mathematica code IfCCZ4[{
mathematica code Name -> dampk1,
mathematica code Description -> "CCZ4 damping term 1 for Theta and Z",
mathematica code Default -> 0
mathematica code }],
mathematica code IfCCZ4[{
mathematica code Name -> dampk2,
mathematica code Description -> "CCZ4 damping term 2 for Theta and Z",
mathematica code Default -> 0
mathematica code }],
mathematica code {
mathematica code Name -> LapseACoeff,
mathematica code Description -> "Whether to evolve A in time",
mathematica code Default -> 0
mathematica code },
mathematica code {
mathematica code Name -> harmonicF,
mathematica code Description -> "d/dt alpha = - f alpha^n K (harmonic=1, 1+log=2)",
mathematica code Default -> 1
mathematica code },
mathematica code {
mathematica code Name -> AlphaDriver,
mathematica code Default -> 0
mathematica code },
mathematica code {
mathematica code Name -> ShiftBCoeff,
mathematica code Description -> "Whether to evolve B^i in time",
mathematica code Default -> 1
mathematica code },
mathematica code {
mathematica code Name -> ShiftGammaCoeff,
mathematica code Default -> 0
mathematica code },
mathematica code {
mathematica code Name -> BetaDriver,
mathematica code Default -> 0
mathematica code },
mathematica code {
mathematica code Name -> LapseAdvectionCoeff,
mathematica code Description -> "Factor in front of the lapse advection terms in 1+log",
mathematica code Default -> 1
mathematica code },
mathematica code {
mathematica code Name -> ShiftAdvectionCoeff,
mathematica code Description -> "Factor in front of the shift advection terms in gamma driver",
mathematica code Default -> 1
mathematica code },
mathematica code {
mathematica code Name -> MinimumLapse,
mathematica code Description -> "Minimum value of the lapse function",
mathematica code Default -> -1
mathematica code },
mathematica code {
mathematica code Name -> SpatialBetaDriverRadius,
mathematica code Description -> "Radius at which the BetaDriver starts to be reduced",
mathematica code AllowedValues -> {{Value -> "(0:*", Description -> "Positive"}},
mathematica code Default -> 10^12
mathematica code },
mathematica code {
mathematica code Name -> SpatialShiftGammaCoeffRadius,
mathematica code Description -> "Radius at which the ShiftGammaCoefficient starts to be reduced",
mathematica code AllowedValues -> {{Value -> "(0:*", Description -> "Positive"}},
mathematica code Default -> 10^12
mathematica code },
mathematica code {
mathematica code Name -> EpsDiss,
mathematica code Description -> "Dissipation strength",
mathematica code AllowedValues -> {{Value -> "(0:*", Description -> "Positive"}},
mathematica code Default -> 0
mathematica code }
mathematica code };
mathematica blank
mathematica comment (******************************************************************************)
mathematica comment (* Construct the thorns *)
mathematica comment (******************************************************************************)
mathematica blank
mathematica code calculations =
mathematica code Join[
mathematica code {
mathematica code initialCalc,
mathematica code convertFromADMBaseCalc,
mathematica code initGammaCalc,
mathematica code convertFromADMBaseGammaCalc,
mathematica code evolCalc,
mathematica code evolCalc1, evolCalc2,
mathematica code dissCalc,
mathematica code advectCalc,
mathematica comment (*advectCalcs,*)
mathematica code initRHSCalc,
mathematica code RHSStaticBoundaryCalc,
mathematica comment (* RHSRadiativeBoundaryCalc, *)
mathematica code enforceCalc,
mathematica code boundaryCalc,
mathematica code convertToADMBaseCalc,
mathematica code convertToADMBaseDtLapseShiftCalc,
mathematica code convertToADMBaseDtLapseShiftBoundaryCalc,
mathematica code convertToADMBaseFakeDtLapseShiftCalc,
mathematica comment (* constraintsCalc, *)
mathematica code constraintsCalc1, constraintsCalc2
mathematica code },
mathematica code advectCalcs
mathematica comment (*dissCalcs*)
mathematica code ];
mathematica blank
mathematica code CreateKrancThornTT [groups, ".", thorn,
mathematica code Calculations -> calculations,
mathematica code DeclaredGroups -> declaredGroupNames,
mathematica code PartialDerivatives -> derivatives,
mathematica code EvolutionTimelevels -> evolutionTimelevels,
mathematica code DefaultEvolutionTimelevels -> 3,
mathematica code UseJacobian -> True,
mathematica code UseLoopControl -> True,
mathematica code UseVectors -> useVectors,
mathematica code UseOpenCL -> useOpenCL,
mathematica code InheritedImplementations -> inheritedImplementations,
mathematica code InheritedKeywordParameters -> inheritedKeywordParameters,
mathematica code ExtendedKeywordParameters -> extendedKeywordParameters,
mathematica code KeywordParameters -> keywordParameters,
mathematica code IntParameters -> intParameters,
mathematica code RealParameters -> realParameters
mathematica code ];
mathematica blank
mathematica code ];
mathematica blank
mathematica blank
mathematica blank
mathematica comment (******************************************************************************)
mathematica comment (* Options *)
mathematica comment (******************************************************************************)
mathematica blank
mathematica comment (* These are the arguments to createComment:
mathematica comment - derivative order: 2, 4, 6, 8, ...
mathematica comment - useJacobian: False or True
mathematica comment - split upwind derivatives: False or True
mathematica comment - use vectorisation: False or True
mathematica comment - use OpenCL: False or True
mathematica comment - timelevels: 2 or 3
mathematica comment ## (keep this at 3; this is better chosen with a run-time parameter)
mathematica comment - matter: 0 or 1
mathematica comment ## (matter seems cheap; it should be always enabled)
mathematica comment - thorn base name
mathematica comment *)
mathematica blank
mathematica code createCode[4, False, True, True , False, 3, 1, "BSSN"];
mathematica code createCode[4, False, True, False, False, 3, 1, "BSSN"];
mathematica code createCode[4, False, True, True , True , 3, 1, "BSSN"];
mathematica blank
mathematica code createCode[4, False, True, True , False, 3, 1, "CCZ4"];
mathematica blank