From e636585d71de3af8c961319e93ceb3bd9e31582f Mon Sep 17 00:00:00 2001 From: Ian Hinder Date: Thu, 5 Sep 2013 19:29:23 +0200 Subject: Tests: Add McLachlan.mt This exercises a large amount of code, but takes about 2 minutes to run, even with NoSimplify -> True. --- Tests/McLachlan.mt | 1459 ++++++++++++++++++++++++++++++++++++++++++++++++++++ Tests/RunTests.m | 3 +- 2 files changed, 1461 insertions(+), 1 deletion(-) create mode 100644 Tests/McLachlan.mt (limited to 'Tests') diff --git a/Tests/McLachlan.mt b/Tests/McLachlan.mt new file mode 100644 index 0000000..5eeea8d --- /dev/null +++ b/Tests/McLachlan.mt @@ -0,0 +1,1459 @@ + +(* Mathematica Test File *) + +SetEnhancedTimes[False]; +SetSourceLanguage["C"]; + +(******************************************************************************) +(* Options *) +(******************************************************************************) + +createCode[derivOrder_, useJacobian_, splitUpwindDerivs_, evolutionTimelevels_, addMatter_, formulation_] := +Module[{prefix, suffix, thorn}, + +prefix = "ML_"; +suffix = + "" + <> If [useJacobian, "_MP", ""] + <> If [derivOrder!=4, "_O" <> ToString[derivOrder], ""] + <> If [splitUpwindDerivs, "", "_UPW"] + (* <> 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, Ddetgt, + 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, em4phi, Ddetgt[la]}, + Equations -> + { + dir[ua] -> Sign[beta[ua]], + + detgt -> 1 (* detgtExpr *), + (* This leads to simpler code... *) + gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], + em4phi -> IfThen[conformalMethod==CMW, phi^2, Exp[-4 phi]], + + eta -> etaExpr, + theta -> thetaExpr, + + (* Ddetgt should be zero analytically, but we're not assuming it here. Change commenting to assume it.*) + Ddetgt[la] -> gtu[uk,ul] PD[gt[lk,ll],la], + (*Ddetgt[la] -> 0,*) + + (* 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] em4phi alpha + (- 2 alpha IfThen[conformalMethod==CMW,1/phi,-2] PD[phi,lj] + + 2 PD[alpha,lj] + + alpha (Ddetgt[lj] - 2 gtu[uk,ul] 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, + 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], Ddetgt[la]}, + 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, + + (* Ddetgt should be zero analytically, but we're not assuming it here. Change commenting to assume it.*) + Ddetgt[la] -> gtu[uk,ul] PD[gt[lk,ll],la], + (*Ddetgt[la] -> 0,*) + + (* dot[beta[ua]] -> eta Xt[ua], *) + (* dot[beta[ua]] -> ShiftGammaCoeff alpha^ShiftAlphaPower B[ua], *) + dot[beta[ua]] -> IfThen[harmonicShift, + - 1/2 gtu[ua,uj] em4phi alpha + (- 2 alpha IfThen[conformalMethod==CMW,1/phi,-2] PD[phi,lj] + + 2 PD[alpha,lj] + + alpha (Ddetgt[lj] - 2 gtu[uk,ul] 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 <> "_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, + 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]]. *) + } +}; + +evolCalc1 = PartialCalculation[evolCalc, "1", + { + ConditionalOnKeyword -> {"RHS_calculation", "split"} + }, + { + 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_calculation", "split"} + }, + { + 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", + 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" + } + +}; + +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, + initRHSCalc, + (* evol1Calc, evol2Calc, *) + RHSStaticBoundaryCalc, + (* RHSRadiativeBoundaryCalc, *) + enforceCalc, + boundaryCalc, + convertToADMBaseCalc, + convertToADMBaseDtLapseShiftCalc, + convertToADMBaseDtLapseShiftBoundaryCalc, + convertToADMBaseFakeDtLapseShiftCalc, + (* constraintsCalc, *) + constraintsCalc1, constraintsCalc2 +}, + {} (*dissCalcs*) +]; + +calculations = Map[Append[#, NoSimplify -> True] &, calculations]; + +CreateKrancThornTT [groups, "TestThorns", thorn, + Calculations -> calculations, + DeclaredGroups -> declaredGroupNames, + PartialDerivatives -> derivatives, + EvolutionTimelevels -> evolutionTimelevels, + DefaultEvolutionTimelevels -> 3, + UseJacobian -> True, + UseLoopControl -> True, + UseVectors -> True, + InheritedImplementations -> inheritedImplementations, + InheritedKeywordParameters -> inheritedKeywordParameters, + ExtendedKeywordParameters -> extendedKeywordParameters, + KeywordParameters -> keywordParameters, + IntParameters -> intParameters, + RealParameters -> realParameters +]; + +]; + +(****************************************************************) +(* McLachlan *) +(****************************************************************) + +Test[ + +(******************************************************************************) +(* Options *) +(******************************************************************************) + +(* These are the arguments to createCode: + - derivative order: 2, 4, 6, 8, ... + - useJacobian: False or True + - split upwind derivatives: 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 , 3, 1, "BSSN"]; +(* createCode[4, False, True , 3, 1, "CCZ4"]; *) + , + Null + , + TestID->"McLachlan" +] diff --git a/Tests/RunTests.m b/Tests/RunTests.m index 46fee00..5d749ac 100755 --- a/Tests/RunTests.m +++ b/Tests/RunTests.m @@ -21,7 +21,8 @@ Needs["KrancThorn`"]; SetDebugLevel[DebugQuiet]; tests = { - "Kranc" + "Kranc", + "McLachlan" }; (Print["\n"]; TestRun[#<>".mt", Loggers -> {VerbosePrintLogger[]}, TestRunTitle -> #]) & /@ tests; -- cgit v1.2.3