(* Mathematica Test File *) SetEnhancedTimes[False]; SetSourceLanguage["C"]; (******************************************************************************) (* Options *) (******************************************************************************) Options[createCode] = {"DGFE" -> False, "CaKernel" -> False}; createCode[derivOrder_, useJacobian_, splitUpwindDerivs_, evolutionTimelevels_, addMatter_, formulation_, vectorise_, opencl_, OptionsPattern[]] := Module[{prefix, suffix, thorn}, prefix = "ML_"; suffix = "" <> If [useJacobian, "_MP", ""] <> If [!vectorise, "_NoVec", ""] <> If [opencl, "_OpenCL", ""] <> If [OptionValue[DGFE], "_DGFE", ""] <> If [OptionValue[CaKernel], "_CaKernel", ""] <> 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])], 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])], 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])], (* 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} }, { 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 -> vectorise, UseOpenCL -> opencl, UseDGFE -> OptionValue[DGFE], UseCaKernel -> OptionValue[CaKernel], InheritedImplementations -> inheritedImplementations, InheritedKeywordParameters -> inheritedKeywordParameters, ExtendedKeywordParameters -> extendedKeywordParameters, KeywordParameters -> keywordParameters, IntParameters -> intParameters, RealParameters -> realParameters ]; ]; (****************************************************************) (* McLachlan *) (****************************************************************) (******************************************************************************) (* 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 - vectorise - opencl *) Test[ ClearAllTensors[]; createCode[4, False, True , 3, 1, "BSSN", True, False]; , Null , TestID->"McLachlanVec" ] Test[ ClearAllTensors[]; createCode[4, False, True , 3, 1, "BSSN", False, False]; , Null , TestID->"McLachlanNoVec" ] Test[ ClearAllTensors[]; createCode[4, False, True , 3, 1, "BSSN", True, True]; , Null , TestID->"McLachlanOpenCL" ] Test[ ClearAllTensors[]; createCode[4, False, True , 3, 1, "BSSN", True, False, DGFE -> True]; , Null , TestID->"McLachlanDGFE" ] Test[ ClearAllTensors[]; createCode[4, False, True , 3, 1, "BSSN", True, False, CaKernel -> True]; , Null , TestID->"McLachlanCaKernel" ]