From a8794cd2227c87c7457ba51d040e863be2f981e7 Mon Sep 17 00:00:00 2001 From: Barry Wardell Date: Tue, 1 May 2012 17:09:05 +0100 Subject: Generate a separate thorn for CCZ4. This still uses the same Kranc script as the standard BSSN code, but having a separate thorn means that the BSSN is not negatively affected performance-wise. --- m/Makefile | 3 +- m/McLachlan_BSSN.m | 283 +++++++++++++++++++++++++---------------------------- 2 files changed, 136 insertions(+), 150 deletions(-) diff --git a/m/Makefile b/m/Makefile index 3af082a..cf8260b 100644 --- a/m/Makefile +++ b/m/Makefile @@ -13,13 +13,14 @@ McLachlan_ADM.out: McLachlan_ADM.m done McLachlan_BSSN.out: McLachlan_BSSN.m - rm -rf ML_BSSN* + rm -rf ML_BSSN* ML_CCZ4* ./runmath.sh $^ for thorn in ML_BSSN*; do \ ./copy-if-changed.sh $$thorn ../$$thorn && \ ./create-helper-thorn.sh $$thorn && \ ./copy-if-changed.sh $${thorn}_Helper ../$${thorn}_Helper; \ done + ./copy-if-changed.sh ML_CCZ4 ../ML_CCZ4 McLachlan_ADMConstraints.out: McLachlan_ADMConstraints.m rm -rf ML_ADMConstraints* diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m index 16d076a..0ce97b7 100644 --- a/m/McLachlan_BSSN.m +++ b/m/McLachlan_BSSN.m @@ -6,8 +6,8 @@ SetSourceLanguage["C"]; (* Options *) (******************************************************************************) -createCode[derivOrder_, useJacobian_, splitUpwindDerivs_, evolutionTimelevels_, addMatter_] := -Module[{}, +createCode[derivOrder_, useJacobian_, splitUpwindDerivs_, evolutionTimelevels_, addMatter_, formulation_] := +Module[{prefix, suffix, thorn}, prefix = "ML_"; suffix = @@ -19,9 +19,10 @@ suffix = (* <> If [addMatter==1, "_M", ""] *) ; -BSSN = prefix <> "BSSN" <> suffix; - +thorn = prefix <> formulation <> suffix; +SetAttributes[IfCCZ4, HoldAll]; +IfCCZ4[expr_, else_:Sequence[]] := If[formulation === "CCZ4", expr, Unevaluated[else]]; (******************************************************************************) (* Derivatives *) @@ -165,11 +166,6 @@ T11=eTxx; T12=eTxy; T22=eTyy; T13=eTxz; T23=eTyz; T33=eTzz; CMBSSNphi = 0; CMBSSNW = 1; -(* enum constants for formulation; these must be consistent - with the definition of the Cactus parameter formulation *) -FBSSN = 0; -FCCZ4 = 1; - detgExpr = Det [MatrixOfComponents [g [la,lb]]]; ddetgExpr[la_] = Sum [D[Det[MatrixOfComponents[g[la, lb]]], X] PD[X, la], @@ -193,13 +189,13 @@ evolvedGroups = {SetGroupName [CreateGroupFromTensor [phi ], prefix <> "log_confac"], SetGroupName [CreateGroupFromTensor [gt[la,lb]], prefix <> "metric" ], SetGroupName [CreateGroupFromTensor [Xt[ua] ], prefix <> "Gamma" ], - SetGroupName [CreateGroupFromTensor [Theta ], prefix <> "Theta" ], 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" ]}; + SetGroupName [CreateGroupFromTensor [B[ua] ], prefix <> "dtshift" ], + IfCCZ4[SetGroupName[CreateGroupFromTensor[Theta], prefix <> "Theta"]]}; evaluatedGroups = {SetGroupName [CreateGroupFromTensor [H ], prefix <> "Ham"], SetGroupName [CreateGroupFromTensor [M[la] ], prefix <> "mom"], @@ -235,7 +231,7 @@ groups = Join [declaredGroups, extraGroups]; initialCalc = { - Name -> BSSN <> "_Minkowski", + Name -> thorn <> "_Minkowski", Schedule -> {"IN ADMBase_InitialData"}, ConditionalOnKeyword -> {"my_initial_data", "Minkowski"}, Equations -> @@ -245,11 +241,11 @@ initialCalc = trK -> 0, At[la,lb] -> 0, Xt[ua] -> 0, - Theta -> 0, alpha -> 1, A -> 0, beta[ua] -> 0, - B[ua] -> 0 + B[ua] -> 0, + IfCCZ4[Theta -> 0] } }; @@ -286,7 +282,7 @@ Module[ convertFromADMBaseCalc = { - Name -> BSSN <> "_convertFromADMBase", + Name -> thorn <> "_convertFromADMBase", Schedule -> {"AT initial AFTER ADMBase_PostInitial"}, ConditionalOnKeyword -> {"my_initial_data", "ADMBase"}, Shorthands -> {g[la,lb], detg, gu[ua,ub], em4phi}, @@ -303,18 +299,18 @@ convertFromADMBaseCalc = trK -> gu[ua,ub] admK[la,lb], At[la,lb] -> em4phi (admK[la,lb] - (1/3) g[la,lb] trK), - Theta -> 0, - alpha -> admalpha, - beta[ua] -> admbeta[ua] + beta[ua] -> admbeta[ua], + + IfCCZ4[Theta -> 0] } }; convertFromADMBaseGammaCalc = { - Name -> BSSN <> "_convertFromADMBaseGamma", - Schedule -> {"AT initial AFTER " <> BSSN <> "_convertFromADMBase"}, + Name -> thorn <> "_convertFromADMBaseGamma", + Schedule -> {"AT initial AFTER " <> thorn <> "_convertFromADMBase"}, ConditionalOnKeyword -> {"my_initial_data", "ADMBase"}, (* Where -> InteriorNoSync, @@ -370,8 +366,8 @@ convertFromADMBaseGammaCalc = initialised. *) initGammaCalc = { - Name -> BSSN <> "_InitGamma", - Schedule -> {"AT initial BEFORE " <> BSSN <> "_convertFromADMBaseGamma"}, + Name -> thorn <> "_InitGamma", + Schedule -> {"AT initial BEFORE " <> thorn <> "_convertFromADMBaseGamma"}, ConditionalOnKeyword -> {"my_initial_data", "ADMBase"}, Where -> Everywhere, Equations -> @@ -390,13 +386,13 @@ initGammaCalc = convertToADMBaseCalc = { - Name -> BSSN <> "_convertToADMBase", - Schedule -> {"IN " <> BSSN <> "_convertToADMBaseGroup"}, + Name -> thorn <> "_convertToADMBase", + Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"}, Where -> Everywhere, Shorthands -> {e4phi}, Equations -> { - e4phi -> IfThen[conformalMethod==CMBSSNW, 1/phi^2, Exp[4 phi]], + e4phi -> IfThen[conformalMethod, 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, @@ -406,8 +402,8 @@ convertToADMBaseCalc = convertToADMBaseDtLapseShiftCalc = { - Name -> BSSN <> "_convertToADMBaseDtLapseShift", - Schedule -> {"IN " <> BSSN <> "_convertToADMBaseGroup"}, + Name -> thorn <> "_convertToADMBaseDtLapseShift", + Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"}, ConditionalOnKeyword -> {"dt_lapse_shift_method", "correct"}, Where -> Interior, Shorthands -> {dir[ua], detgt, gtu[ua,ub], eta, theta}, @@ -431,7 +427,7 @@ convertToADMBaseDtLapseShiftCalc = admdtalpha -> - harmonicF alpha^harmonicN (+ LapseACoeff A + ((1 - LapseACoeff) - (trK - IfThen[formulation==FCCZ4, 2 Theta, 0]))) + (trK - IfCCZ4[2 Theta, 0]))) + LapseAdvectionCoeff Upwind[beta[ua], alpha, la], admdtbeta[ua] -> IfThen[harmonicShift, - 1/2 gtu[ua,uj] phi alpha @@ -450,8 +446,8 @@ convertToADMBaseDtLapseShiftCalc = convertToADMBaseDtLapseShiftBoundaryCalc = { - Name -> BSSN <> "_convertToADMBaseDtLapseShiftBoundary", - Schedule -> {"IN " <> BSSN <> "_convertToADMBaseGroup"}, + Name -> thorn <> "_convertToADMBaseDtLapseShiftBoundary", + Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"}, ConditionalOnKeyword -> {"dt_lapse_shift_method", "correct"}, Where -> BoundaryWithGhosts, Shorthands -> {detgt, gtu[ua,ub], eta, theta}, @@ -472,7 +468,7 @@ convertToADMBaseDtLapseShiftBoundaryCalc = admdtalpha -> - harmonicF alpha^harmonicN (+ LapseACoeff A + ((1 - LapseACoeff) - (trK - IfThen[formulation==FCCZ4, 2 Theta, 0]))), + (trK - IfCCZ4[2 Theta, 0]))), admdtbeta[ua] -> IfThen[harmonicShift, 0, (* else *) @@ -485,8 +481,8 @@ convertToADMBaseDtLapseShiftBoundaryCalc = convertToADMBaseFakeDtLapseShiftCalc = { - Name -> BSSN <> "_convertToADMBaseFakeDtLapseShift", - Schedule -> {"IN " <> BSSN <> "_convertToADMBaseGroup"}, + Name -> thorn <> "_convertToADMBaseFakeDtLapseShift", + Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"}, ConditionalOnKeyword -> {"dt_lapse_shift_method", "noLapseShiftAdvection"}, Where -> Everywhere, Shorthands -> {detgt, gtu[ua,ub], eta, theta}, @@ -509,7 +505,7 @@ convertToADMBaseFakeDtLapseShiftCalc = admdtalpha -> - harmonicF alpha^harmonicN (+ LapseACoeff A + ((1 - LapseACoeff) - (trK - IfThen[formulation==FCCZ4, 2 Theta, 0]))), + (trK - IfCCZ4[2 Theta, 0]))), admdtbeta[ua] -> IfThen[harmonicShift, 0, (* else *) @@ -526,8 +522,8 @@ convertToADMBaseFakeDtLapseShiftCalc = evolCalc = { - Name -> BSSN <> "_RHS", - Schedule -> {"IN " <> BSSN <> "_evolCalcGroup"}, + Name -> thorn <> "_RHS", + Schedule -> {"IN " <> thorn <> "_evolCalcGroup"}, (* Where -> Interior, *) @@ -536,14 +532,14 @@ evolCalc = radiative boundary conditions. *) Where -> InteriorNoSync, Shorthands -> {dir[ua], - detgt, gtu[ua,ub], Z[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], dotTheta, - epsdiss[ua]}, + rho, S[la], trS, fac1, fac2, dottrK, dotXt[ua], + epsdiss[ua], IfCCZ4[Z[ua]], IfCCZ4[dotTheta]}, Equations -> { dir[ua] -> Sign[beta[ua]], @@ -569,10 +565,9 @@ evolCalc = (* The Z quantities *) (* gr-qc:1106.2254 (2011), eqn. (23) *) - Z[ud] -> IfThen[formulation==FCCZ4, - (1/2) gu[ua,ud] (- PD[gt[la,lb],lc] gtu[ub,uc] - + gt[la,lc] Xt[uc]), - 0], + 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 *) @@ -599,13 +594,12 @@ evolCalc = 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] - + IfThen[formulation==FCCZ4, - + (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], - 0], + 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 *) @@ -629,10 +623,7 @@ evolCalc = (* 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] - - IfThen[formulation==FCCZ4, - (1/3) At[lc,ld] gtu[uc,ud] gt[la,lb], - 0]) + 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) *) @@ -648,7 +639,7 @@ evolCalc = + (1/3) gtu[ui,uj] PD[beta[ul],lj,ll] - Xtn[uj] PD[beta[ui],lj] + (2/3) Xtn[ui] PD[beta[uj],lj] - + IfThen[formulation==FCCZ4, + + IfCCZ4[ + IfThen[GammaShift, 2 e4phi (- Z[uj] PD[beta[ui],lj] + (2/3) Z[ui] PD[beta[uj],lj]), @@ -657,21 +648,21 @@ evolCalc = + 2 gtu[ui,uj] (+ alpha PD[Theta,lj] - Theta PD[alpha,lj]) - 2 alpha e4phi dampk1 Z[ui], - 0] + 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) *) - dotTheta -> IfThen[formulation==FCCZ4, - - PD[alpha,la] Z[ua] - + (1/2) alpha (gu[ua,ub] R[la,lb] - - Atm[ua,lb] Atm[ub,la] - + (2/3) trK^2 - 2 trK Theta) - - dampk1 (2 + dampk2) alpha Theta, - 0], + 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) + ], - dot[Theta] -> dotTheta, + IfCCZ4[ + dot[Theta] -> dotTheta + ], (* PRD 62, 044034 (2000), eqn. (11) *) (* gr-qc:1106.2254 (2011), eqn. (17) *) @@ -681,7 +672,7 @@ evolCalc = + 2 cdphi[la] PD[alpha,lb] ) - Xtn[ua] PD[alpha,la] ) + alpha (Atm[ua,lb] Atm[ub,la] + (1/3) trK^2) - + IfThen[formulation==FCCZ4, + + IfCCZ4[ + 2 dotTheta + 2 PD[alpha,la] Z[ua] + dampk1 (1 - dampk2) alpha Theta, 0] @@ -698,9 +689,7 @@ evolCalc = + 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 - IfThen[formulation==FCCZ4, - 2 Theta, - 0]) + + 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] @@ -723,11 +712,11 @@ evolCalc = dot[alpha] -> - harmonicF alpha^harmonicN (+ LapseACoeff A + ((1 - LapseACoeff) - (+ trK - IfThen[formulation==FCCZ4, 2 Theta, 0] + (+ trK - IfCCZ4[2 Theta, 0] + AlphaDriver (alpha - 1)))), dot[A] -> + (LapseACoeff - (+ dottrK - IfThen[formulation==FCCZ4, 2 dotTheta, 0] + (+ dottrK - IfCCZ4[2 dotTheta, 0] - AlphaDriver A)), eta -> etaExpr, @@ -754,9 +743,9 @@ evolCalc = advectCalc = { - Name -> BSSN <> "_Advect", - Schedule -> {"IN " <> BSSN <> "_evolCalcGroup " <> - "AFTER (" <> BSSN <> "_RHS1 " <> BSSN <> "_RHS2)"}, + Name -> thorn <> "_Advect", + Schedule -> {"IN " <> thorn <> "_evolCalcGroup " <> + "AFTER (" <> thorn <> "_RHS1 " <> thorn <> "_RHS2)"}, (* Where -> Interior, *) @@ -775,10 +764,9 @@ advectCalc = dot[Xt[ui]] -> dot[Xt[ui]] + Upwind[beta[uj], Xt[ui], lj], - dot[Theta] -> dot[Theta] - + IfThen[formulation==FCCZ4, - Upwind[beta[ua], Theta, la], - 0], + IfCCZ4[ + dot[Theta] -> dot[Theta] + Upwind[beta[ua], Theta, la] + ], dot[trK] -> dot[trK] + Upwind[beta[ua], trK, la], @@ -811,12 +799,12 @@ evolCalc1 = PartialCalculation[evolCalc, "1", dot[phi], dot[gt[la,lb]], dot[Xt[ui]], - dot[Theta], dot[trK], dot[alpha], dot[A], dot[beta[ua]], - dot[B[ua]] + dot[B[ua]], + IfCCZ4[dot[Theta]] }]; evolCalc2 = PartialCalculation[evolCalc, "2", @@ -829,9 +817,9 @@ evolCalc2 = PartialCalculation[evolCalc, "2", dissCalc = { - Name -> BSSN <> "_Dissipation", - Schedule -> {"IN " <> BSSN <> "_evolCalcGroup " <> - "AFTER (" <> BSSN <> "_RHS1 " <> BSSN <> "_RHS2)"}, + Name -> thorn <> "_Dissipation", + Schedule -> {"IN " <> thorn <> "_evolCalcGroup " <> + "AFTER (" <> thorn <> "_RHS1 " <> thorn <> "_RHS2)"}, ConditionalOnKeyword -> {"apply_dissipation", "always"}, Where -> InteriorNoSync, Shorthands -> {epsdiss[ua]}, @@ -840,7 +828,7 @@ dissCalc = epsdiss[ua] -> EpsDiss, Sequence@@Table[ dot[var] -> dot[var] + epsdiss[ux] PDdiss[var,lx], - {var, {phi, gt[la,lb], Xt[ui], Theta, trK, At[la,lb], + {var, {phi, gt[la,lb], Xt[ui], IfCCZ4[Theta], trK, At[la,lb], alpha, A, beta[ua], B[ua]}}] } }; @@ -848,9 +836,9 @@ dissCalc = dissCalcs = Table[ { - Name -> BSSN <> "_Dissipation_" <> ToString[var /. {Tensor[n_,__] -> n}], - Schedule -> {"IN " <> BSSN <> "_evolCalcGroup " <> - "AFTER (" <> BSSN <> "_RHS1 " <> BSSN <> "_RHS2)"}, + 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]}, @@ -860,13 +848,13 @@ Table[ dot[var] -> dot[var] + epsdiss[ux] PDdiss[var,lx] } }, - {var, {phi, gt[la,lb], Xt[ui], Theta, trK, At[la,lb], + {var, {phi, gt[la,lb], Xt[ui], IfCCZ4[Theta], trK, At[la,lb], alpha, A, beta[ua], B[ua]}} ]; RHSStaticBoundaryCalc = { - Name -> BSSN <> "_RHSStaticBoundary", + Name -> thorn <> "_RHSStaticBoundary", Schedule -> {"IN MoL_CalcRHS"}, ConditionalOnKeyword -> {"my_rhs_boundary_condition", "static"}, Where -> Boundary, @@ -877,11 +865,11 @@ RHSStaticBoundaryCalc = dot[trK] -> 0, dot[At[la,lb]] -> 0, dot[Xt[ua]] -> 0, - dot[Theta] -> 0, dot[alpha] -> 0, dot[A] -> 0, dot[beta[ua]] -> 0, - dot[B[ua]] -> 0 + dot[B[ua]] -> 0, + IfCCZ4[dot[Theta] -> 0] } }; @@ -890,8 +878,8 @@ RHSStaticBoundaryCalc = them to be zero *) initRHSCalc = { - Name -> BSSN <> "_InitRHS", - Schedule -> {"AT analysis BEFORE " <> BSSN <> "_evolCalcGroup"}, + Name -> thorn <> "_InitRHS", + Schedule -> {"AT analysis BEFORE " <> thorn <> "_evolCalcGroup"}, Where -> Everywhere, Equations -> { @@ -900,17 +888,17 @@ initRHSCalc = dot[trK] -> 0, dot[At[la,lb]] -> 0, dot[Xt[ua]] -> 0, - dot[Theta] -> 0, dot[alpha] -> 0, dot[A] -> 0, dot[beta[ua]] -> 0, - dot[B[ua]] -> 0 + dot[B[ua]] -> 0, + IfCCZ4[dot[Theta] -> 0] } }; RHSRadiativeBoundaryCalc = { - Name -> BSSN <> "_RHSRadiativeBoundary", + Name -> thorn <> "_RHSRadiativeBoundary", Schedule -> {"IN MoL_CalcRHS"}, ConditionalOnKeyword -> {"my_rhs_boundary_condition", "radiative"}, Where -> Boundary, @@ -940,17 +928,19 @@ RHSRadiativeBoundaryCalc = 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[Theta] -> - vg su[uc] PDo[Theta ,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] + dot[B[ua]] -> - su[uc] PDo[B[ua] ,lc], + IfCCZ4[ + dot[Theta] -> - vg su[uc] PDo[Theta ,lc] + ] } }; enforceCalc = { - Name -> BSSN <> "_enforce", + Name -> thorn <> "_enforce", Schedule -> {"IN MoL_PostStepModify"}, Shorthands -> {detgt, gtu[ua,ub], trAt}, Equations -> @@ -981,7 +971,7 @@ enforceCalc = boundaryCalc = { - Name -> BSSN <> "_boundary", + Name -> thorn <> "_boundary", Schedule -> {"IN MoL_PostStep"}, ConditionalOnKeyword -> {"my_boundary_condition", "Minkowski"}, Where -> BoundaryWithGhosts, @@ -992,11 +982,11 @@ boundaryCalc = trK -> 0, At[la,lb] -> 0, Xt[ua] -> 0, - Theta -> 0, alpha -> 1, A -> 0, beta[ua] -> 0, - B[ua] -> 0 + B[ua] -> 0, + IfCCZ4[Theta -> 0] } }; @@ -1006,7 +996,7 @@ boundaryCalc = constraintsCalc = { - Name -> BSSN <> "_constraints", + Name -> thorn <> "_constraints", Schedule -> Automatic, After -> "MoL_PostStep", Where -> Interior, @@ -1040,10 +1030,9 @@ constraintsCalc = gu[ua,ub] -> em4phi gtu[ua,ub], (* The Z quantities *) - Z[ud] -> IfThen[formulation==FCCZ4, - (1/2) gu[ua,ud] (- PD[gt[la,lb],lc] gtu[ub,uc] - + gt[la,lc] Xt[uc]), - 0], + 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] @@ -1095,13 +1084,14 @@ constraintsCalc = + 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] - + IfThen[formulation==FCCZ4, - + (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], - 0], + R[la,lb] -> + Rt[la,lb] + Rphi[la,lb], + + IfCCZ4[ + 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, *) @@ -1177,23 +1167,23 @@ extendedKeywordParameters = { { Name -> "ADMBase::evolution_method", - AllowedValues -> {BSSN} + AllowedValues -> {thorn} }, { Name -> "ADMBase::lapse_evolution_method", - AllowedValues -> {BSSN} + AllowedValues -> {thorn} }, { Name -> "ADMBase::shift_evolution_method", - AllowedValues -> {BSSN} + AllowedValues -> {thorn} }, { Name -> "ADMBase::dtlapse_evolution_method", - AllowedValues -> {BSSN} + AllowedValues -> {thorn} }, { Name -> "ADMBase::dtshift_evolution_method", - AllowedValues -> {BSSN} + AllowedValues -> {thorn} } }; @@ -1263,25 +1253,11 @@ keywordParameters = intParameters = { - { - Name -> conformalMethod, - Description -> "Treatment of conformal factor", - AllowedValues -> {{Value -> "0", Description -> "phi method"}, - {Value -> "1", Description -> "W method"}}, - Default -> 0 - }, - { - Name -> formulation, - Description -> "Formulation of the equations", - AllowedValues -> {{Value -> "0", Description -> "BSSN"}, - {Value -> "1", Description -> "CCZ4"}}, - Default -> 0 - }, - { + IfCCZ4[{ Name -> GammaShift, Description -> "Covariant shift term in Gamma", Default -> 0 - }, + }], { Name -> harmonicN, Description -> "d/dt alpha = - f alpha^n K (harmonic=2, 1+log=1)", @@ -1292,31 +1268,38 @@ intParameters = Default -> 0 }, { - Name -> harmonicShift, - Description -> "Whether to use the harmonic shift", - AllowedValues -> {{Value -> "0", Description -> "Gamma driver shift"}, - {Value -> "1", Description -> "Harmonic shift"}}, + 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 -> 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", @@ -1410,7 +1393,7 @@ Join[ {} (*dissCalcs*) ]; -CreateKrancThornTT [groups, ".", BSSN, +CreateKrancThornTT [groups, ".", thorn, Calculations -> calculations, DeclaredGroups -> declaredGroupNames, PartialDerivatives -> derivatives, @@ -1443,9 +1426,11 @@ CreateKrancThornTT [groups, ".", BSSN, (* matter: 0 or 1 (matter seems cheap; it should be always enabled) *) -createCode[2, False, True , 3, 1]; -createCode[4, False, True , 3, 1]; -createCode[4, False, False, 3, 1]; -createCode[4, True , True , 3, 1]; -createCode[8, False, True , 3, 1]; -createCode[8, True , True , 3, 1]; +createCode[2, False, True , 3, 1, "BSSN"]; +createCode[4, False, True , 3, 1, "BSSN"]; +createCode[4, False, False, 3, 1, "BSSN"]; +createCode[4, True , True , 3, 1, "BSSN"]; +createCode[8, False, True , 3, 1, "BSSN"]; +createCode[8, True , True , 3, 1, "BSSN"]; + +createCode[4, False, True , 3, 1, "CCZ4"]; -- cgit v1.2.3