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