From e574bf527f6e15bc0df239404df428caae93d797 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Fri, 6 Jan 2012 00:00:00 +0000 Subject: Some superficial changes: * Rename Tet -> Theta and Zet -> Z * Introduce an enum to select between CMBSSNphi, CMBSSNW and CMCCZ4. * Rename gamashift parameter to GammaShift --- m/McLachlan_BSSN.m | 281 +++++++++++++++++++++++++++++------------------------ 1 file changed, 155 insertions(+), 126 deletions(-) (limited to 'm') diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m index 6943a3e..e6dd50b 100644 --- a/m/McLachlan_BSSN.m +++ b/m/McLachlan_BSSN.m @@ -103,14 +103,15 @@ Map [DefineTensor, 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, alpha, A, beta, B, Atm, Atu, trA, Ats, trAts, + phi, gt, At, Xt, Xtn, Theta, Z, + alpha, A, beta, B, Atm, Atu, trA, Ats, trAts, dottrK, dotXt, cXt, cS, cA, e4phi, em4phi, ddetg, detgt, gtu, ddetgt, dgtu, ddgtu, Gtl, Gtlu, Gt, Rt, Rphi, gK, T00, T0, T, rho, S, x, y, z, r, - epsdiss, Tet, Zet, dampk1k2}]; + epsdiss}]; (* NOTE: It seems as if Lie[.,.] did not take these tensor weights into account. Presumably, CD[.,.] and CDt[.,.] don't do this either. *) @@ -159,6 +160,12 @@ 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 *) +CMBSSNphi = 0; +CMBSSNW = 1; +CMCCZ4 = 2; + detgExpr = Det [MatrixOfComponents [g [la,lb]]]; ddetgExpr[la_] = Sum [D[Det[MatrixOfComponents[g[la, lb]]], X] PD[X, la], @@ -182,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 [Tet ], prefix <> "Zenergy" ]}; + SetGroupName [CreateGroupFromTensor [B[ua] ], prefix <> "dtshift" ]}; evaluatedGroups = {SetGroupName [CreateGroupFromTensor [H ], prefix <> "Ham"], SetGroupName [CreateGroupFromTensor [M[la] ], prefix <> "mom"], @@ -229,16 +236,16 @@ initialCalc = ConditionalOnKeyword -> {"my_initial_data", "Minkowski"}, Equations -> { - phi -> IfThen [conformalMethod, 1, 0], + phi -> IfThen[conformalMethod==CMBSSNW, 1, 0], gt[la,lb] -> KD[la,lb], trK -> 0, At[la,lb] -> 0, Xt[ua] -> 0, + Theta -> 0, alpha -> 1, A -> 0, beta[ua] -> 0, - B[ua] -> 0, - Tet -> 0 + B[ua] -> 0 } }; @@ -285,13 +292,15 @@ convertFromADMBaseCalc = detg -> detgExpr, gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]], - phi -> IfThen [conformalMethod, detg^(-1/6), Log[detg]/12], - em4phi -> IfThen [conformalMethod, phi^2, Exp[-4 phi]], + phi -> IfThen[conformalMethod==CMBSSNW, detg^(-1/6), Log[detg]/12], + em4phi -> IfThen[conformalMethod==CMBSSNW, 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), + Theta -> 0, + alpha -> admalpha, beta[ua] -> admbeta[ua] @@ -365,8 +374,7 @@ initGammaCalc = { Xt[ua] -> 0, A -> 0, - B[ua] -> 0, - Tet -> 0 + B[ua] -> 0 } }; @@ -384,7 +392,7 @@ convertToADMBaseCalc = Shorthands -> {e4phi}, Equations -> { - e4phi -> IfThen [conformalMethod, 1/phi^2, Exp[4 phi]], + e4phi -> IfThen[conformalMethod==CMBSSNW, 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, @@ -413,13 +421,13 @@ convertToADMBaseDtLapseShiftCalc = (* see RHS *) (* admdtalpha -> - harmonicF alpha^harmonicN - ((1 - LapseAdvectionCoeff) A - + LapseAdvectionCoeff (trK - 2 Tet)) + ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK) + LapseAdvectionCoeff beta[ua] PDu[alpha,la], *) admdtalpha -> - harmonicF alpha^harmonicN (+ LapseACoeff A - + (1 - LapseACoeff) (trK - 2 Tet)) + + ((1 - LapseACoeff) + (trK - IfThen[conformalMethod==CMCCZ4, 2 Theta, 0]))) + LapseAdvectionCoeff Upwind[beta[ua], alpha, la], admdtbeta[ua] -> IfThen[harmonicShift, - 1/2 gtu[ua,uj] phi alpha @@ -455,12 +463,12 @@ convertToADMBaseDtLapseShiftBoundaryCalc = (* see RHS, but omit derivatives near the boundary *) (* admdtalpha -> - harmonicF alpha^harmonicN - ((1 - LapseAdvectionCoeff) A - + LapseAdvectionCoeff (trK - 2 Tet)), + ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK), *) admdtalpha -> - harmonicF alpha^harmonicN (+ LapseACoeff A - + (1 - LapseACoeff) (trK - 2 Tet)), + + ((1 - LapseACoeff) + (trK - IfThen[conformalMethod==CMCCZ4, 2 Theta, 0]))), admdtbeta[ua] -> IfThen[harmonicShift, 0, (* else *) @@ -492,12 +500,12 @@ convertToADMBaseFakeDtLapseShiftCalc = conditions) *) (* admdtalpha -> - harmonicF alpha^harmonicN - ((1 - LapseAdvectionCoeff) A - + LapseAdvectionCoeff (trK - 2 Tet)), + ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK), *) admdtalpha -> - harmonicF alpha^harmonicN (+ LapseACoeff A - + (1 - LapseACoeff) (trK - 2 Tet)), + + ((1 - LapseACoeff) + (trK - IfThen[conformalMethod==CMCCZ4, 2 Theta, 0]))), admdtbeta[ua] -> IfThen[harmonicShift, 0, (* else *) @@ -524,14 +532,14 @@ evolCalc = radiative boundary conditions. *) Where -> InteriorNoSync, Shorthands -> {dir[ua], - detgt, gtu[ua,ub], + detgt, gtu[ua,ub], Z[ua], 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], dotTet, Zet[ua], dampk1k2}, + rho, S[la], trS, fac1, fac2, dottrK, dotXt[ua], dotTheta, + epsdiss[ua]}, Equations -> { dir[ua] -> Sign[beta[ua]], @@ -549,7 +557,7 @@ evolCalc = used instead of Xt where no derivatives of Xt are taken *) Xtn[ui] -> gtu[uj,uk] Gt[ui,lj,lk], - e4phi -> IfThen [conformalMethod, 1/phi^2, Exp[4 phi]], + e4phi -> IfThen[conformalMethod==CMBSSNW, 1/phi^2, Exp[4 phi]], em4phi -> 1 / e4phi, g[la,lb] -> e4phi gt[la,lb], detg -> detgExpr, @@ -557,13 +565,13 @@ evolCalc = (* The Z quantities *) (* gr-qc:1106.2254 (2011), eqn. (23) *) - Zet[ud] -> (1/2) gu[ua,ud] (-PD[gt[la,lb],lc] gtu[ub,uc] + gt[la,lc] Xt[uc]), - - (* shorthand for the product of damping terms k1 and k2 in gr-qc:1106.2254 (2011) *) - dampk1k2 -> dampk1 dampk2, + Z[ud] -> IfThen[conformalMethod==CMCCZ4, + (1/2) gu[ua,ud] (- PD[gt[la,lb],lc] gtu[ub,uc] + + gt[la,lc] Xt[uc]), + 0], (* PRD 62, 044034 (2000), eqn. (18) *) - (* Added Zet term by changing Xtn to Xt *) + (* 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] @@ -573,9 +581,9 @@ evolCalc = + Gt[uk,lj,ll] Gtlu[li,lk,ul] + Gt[uk,li,ll] Gtlu[lk,lj,ul]), - fac1 -> IfThen [conformalMethod, -1/(2 phi), 1], + fac1 -> IfThen[conformalMethod==CMBSSNW, -1/(2 phi), 1], cdphi[la] -> fac1 CDt[phi,la], - fac2 -> IfThen [conformalMethod, 1/(2 phi^2), 0], + fac2 -> IfThen[conformalMethod==CMBSSNW, 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) *) @@ -587,13 +595,14 @@ evolCalc = Atm[ua,lb] -> gtu[ua,uc] At[lc,lb], Atu[ua,ub] -> Atm[ua,lc] gtu[ub,uc], - (* Added Zet terms *) - R[la,lb] -> Rt[la,lb] + Rphi[la,lb] + - (2/phi) (g[la,lc] Zet[uc] PD[phi,lb] + - g[lb,lc] Zet[uc] PD[phi,la] - - g[la,lb] Zet[uc] PD[phi,lc]) + - e4phi Zet[uc] PD[gt[la,lb],lc], - + R[la,lb] -> + Rt[la,lb] + Rphi[la,lb] + + IfThen[conformalMethod==CMCCZ4, + + (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], + (* Matter terms *) (* rho = n^a n^b T_ab *) @@ -610,20 +619,23 @@ evolCalc = (* PRD 62, 044034 (2000), eqn. (10) *) (* PRD 67 084023 (2003), eqn. (16) and (23) *) - dot[phi] -> IfThen [conformalMethod, 1/3 phi, -1/6] alpha trK - + IfThen [conformalMethod, -1/3 phi, 1/6] PD[beta[ua],la], + dot[phi] -> IfThen[conformalMethod==CMBSSNW, 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 insures that detg = 1 *) - dot[gt[la,lb]] -> - 2 alpha (At[la,lb] - - (1/3) At[lc,ld] gtu[uc,ud] gt[la,lb]) + (* removing trA from Aij ensures that detg = 1 *) + dot[gt[la,lb]] -> - 2 alpha (+ At[la,lb] + - IfThen[conformalMethod==CMCCZ4, + (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) *) - (* Added Zet terms by changing Xtn to Xt, and extra Zet and Tet terms *) + (* 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] @@ -632,34 +644,43 @@ 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] - + gamashift (2 e4phi) (- Zet[uj] PD[beta[ui],lj] - + (2/3) Zet[ui] PD[beta[uj],lj]) - - (4/3) alpha e4phi Zet[ui] trK - + 2 gtu[ui,uj] (alpha PD[Tet,lj] - Tet PD[alpha,lj]) - - 2 alpha e4phi dampk1 Zet[ui] + + IfThen[conformalMethod==CMCCZ4, + + IfThen[GammaShift, + 2 e4phi (- Z[uj] PD[beta[ui],lj] + + (2/3) Z[ui] PD[beta[uj],lj]), + 0] + - (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], - (* Equation for Theta *) (* gr-qc:1106.2254 (2011), eqn. (18) *) - dotTet -> - PD[alpha,la] Zet[ua] - + (1/2) alpha (gu[ua,ub] R[la,lb] - - Atm[ua,lb] Atm[ub,la] + (2/3) trK^2 - 2 trK Tet) - - dampk1 2 alpha Tet - dampk1k2 alpha Tet, + dotTheta -> IfThen[conformalMethod==CMCCZ4, + - 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], - dot[Tet] -> dotTet, + dot[Theta] -> dotTheta, (* PRD 62, 044034 (2000), eqn. (11) *) (* gr-qc:1106.2254 (2011), eqn. (17) *) - (* Added the rhs of Tet to K, because K_Z4 = K_BSSN + 2 Tet. *) - (* And the Z term, as it has to cancel with the one in Tet *) + (* 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) - + 2 dotTet + 2 PD[alpha,la] Zet[ua] - + dampk1 alpha Tet - dampk1k2 alpha Tet + + IfThen[conformalMethod==CMCCZ4, + + 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, @@ -667,17 +688,19 @@ evolCalc = (* PRD 62, 044034 (2000), eqn. (12) *) (* TODO: Should we use the Hamiltonian constraint to make Rij tracefree? *) (* gr-qc:1106.2254 (2011), eqn. (15) *) - (* Added Zet terms in the Ricci and Tet terms *) + (* 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 - 2 Tet) At[la,lb] - 2 At[la,lc] Atm[uc,lb]) + + alpha (+ ((trK - IfThen[conformalMethod==CMCCZ4, + 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] - (* damping term in trA, alternative to trA cleaning *) - (* - (dampA/3) alpha At[lc,ld] gtu[uc,ud] gt[la,lb] *) (* 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)), @@ -686,19 +709,22 @@ evolCalc = (* dot[alpha] -> - harmonicF alpha^harmonicN A + Lie[alpha, beta], *) (* dot[alpha] -> - harmonicF alpha^harmonicN ( - (1 - LapseAdvectionCoeff) A - + LapseAdvectionCoeff (trK - 2 Tet)) + (1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK) + LapseAdvectionCoeff beta[ua] PDu[alpha,la], - dot[A] -> (1 - LapseAdvectionCoeff) (dottrK - 2 dotTet - AlphaDriver A), + dot[A] -> (1 - LapseAdvectionCoeff) (dottrK - AlphaDriver A), *) dot[alpha] -> - harmonicF alpha^harmonicN (+ LapseACoeff A - + (1 - LapseACoeff) (trK + AlphaDriver (alpha - 1) - 2 Tet)), + + ((1 - LapseACoeff) + (+ trK - IfThen[conformalMethod==CMCCZ4, 2 Theta, 0] + + AlphaDriver (alpha - 1)))), - dot[A] -> + LapseACoeff (dottrK - 2 dotTet - AlphaDriver A), + dot[A] -> + (LapseACoeff + (+ dottrK - IfThen[conformalMethod==CMCCZ4, 2 dotTheta, 0] + - AlphaDriver A)), eta -> etaExpr, theta -> thetaExpr, @@ -745,7 +771,10 @@ advectCalc = dot[Xt[ui]] -> dot[Xt[ui]] + Upwind[beta[uj], Xt[ui], lj], - dot[Tet] -> dot[Tet] + Upwind[beta[ua], Tet, la], + dot[Theta] -> dot[Theta] + + IfThen[conformalMethod==CMCCZ4, + Upwind[beta[ua], Theta, la], + 0], dot[trK] -> dot[trK] + Upwind[beta[ua], trK, la], @@ -775,11 +804,11 @@ evolCalc1 = PartialCalculation[evolCalc, "1", ConditionalOnKeyword -> {"RHS_calculation", "split"} }, { - dot[Tet], - dot[trK], dot[phi], dot[gt[la,lb]], dot[Xt[ui]], + dot[Theta], + dot[trK], dot[alpha], dot[A], dot[beta[ua]], @@ -807,7 +836,8 @@ dissCalc = epsdiss[ua] -> EpsDiss, Sequence@@Table[ dot[var] -> dot[var] + epsdiss[ux] PDdiss[var,lx], - {var, {phi, gt[la,lb], Xt[ui], trK, Tet, At[la,lb], alpha, A, beta[ua], B[ua]}}] + {var, {phi, gt[la,lb], Xt[ui], Theta, trK, At[la,lb], + alpha, A, beta[ua], B[ua]}}] } }; @@ -826,7 +856,8 @@ Table[ dot[var] -> dot[var] + epsdiss[ux] PDdiss[var,lx] } }, - {var, {phi, gt[la,lb], Xt[ui], trK, At[la,lb], alpha, A, beta[ua], B[ua]}} + {var, {phi, gt[la,lb], Xt[ui], Theta, trK, At[la,lb], + alpha, A, beta[ua], B[ua]}} ]; RHSStaticBoundaryCalc = @@ -839,10 +870,10 @@ RHSStaticBoundaryCalc = { dot[phi] -> 0, dot[gt[la,lb]] -> 0, - dot[Tet] -> 0, 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, @@ -862,10 +893,10 @@ initRHSCalc = { dot[phi] -> 0, dot[gt[la,lb]] -> 0, - dot[Tet] -> 0, 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, @@ -889,7 +920,7 @@ RHSRadiativeBoundaryCalc = detgt -> 1 (* detgtExpr *), gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], - em4phi -> IfThen [conformalMethod, phi^2, Exp[-4 phi]], + em4phi -> IfThen[conformalMethod==CMBSSNW, phi^2, Exp[-4 phi]], gu[ua,ub] -> em4phi gtu[ua,ub], nn[la] -> Euc[la,lb] normal[ub], @@ -902,10 +933,10 @@ RHSRadiativeBoundaryCalc = dot[phi] -> - vg su[uc] PDo[phi ,lc], dot[gt[la,lb]] -> - su[uc] PDo[gt[la,lb],lc], - dot[Tet] -> - vg su[uc] PDo[Tet ,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[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], @@ -952,12 +983,12 @@ boundaryCalc = Where -> BoundaryWithGhosts, Equations -> { - phi -> IfThen [conformalMethod, 1, 0], + phi -> IfThen[conformalMethod==CMBSSNW, 1, 0], gt[la,lb] -> KD[la,lb], - Tet -> 0, trK -> 0, At[la,lb] -> 0, Xt[ua] -> 0, + Theta -> 0, alpha -> 1, A -> 0, beta[ua] -> 0, @@ -975,13 +1006,13 @@ constraintsCalc = Schedule -> Automatic, After -> "MoL_PostStep", Where -> Interior, - Shorthands -> {detgt, ddetgt[la], gtu[ua,ub], + 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, Zet[ua]}, + rho, S[la], fac1, fac2}, Equations -> { detgt -> 1 (* detgtExpr *), @@ -998,17 +1029,19 @@ constraintsCalc = used instead of Xt where no derivatives of Xt are taken *) Xtn[ui] -> gtu[uj,uk] Gt[ui,lj,lk], - e4phi -> IfThen [conformalMethod, 1/phi^2, Exp[4 phi]], + e4phi -> IfThen[conformalMethod==CMBSSNW, 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 *) - Zet[ud] -> (1/2) gu[ua,ud] (-PD[gt[la,lb],lc] gtu[ub,uc] + gt[la,lc] Xt[uc]), - + Z[ud] -> IfThen[conformalMethod==CMCCZ4, + (1/2) gu[ua,ud] (- PD[gt[la,lb],lc] gtu[ub,uc] + + gt[la,lc] Xt[uc]), + 0], + (* PRD 62, 044034 (2000), eqn. (18) *) - (* Added Zet 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] @@ -1040,9 +1073,9 @@ constraintsCalc = - PD[gt[la,lb],l1,l2] + PD[gt[l2,lb],l1,la]), *) - fac1 -> IfThen [conformalMethod, -1/(2 phi), 1], + fac1 -> IfThen[conformalMethod==CMBSSNW, -1/(2 phi), 1], cdphi[la] -> fac1 CDt[phi,la], - fac2 -> IfThen [conformalMethod, 1/(2 phi^2), 0], + fac2 -> IfThen[conformalMethod==CMBSSNW, 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) *) @@ -1058,13 +1091,13 @@ constraintsCalc = + 1/(2 detg) (+ KD[ua,lb] ddetg[lc] + KD[ua,lc] ddetg[lb] - (1/3) g[lb,lc] gu[ua,ud] ddetg[ld]), - (* Added Zet terms *) - R[la,lb] -> Rt[la,lb] + Rphi[la,lb] + - (2/phi) (g[la,lc] Zet[uc] PD[phi,lb] + - g[lb,lc] Zet[uc] PD[phi,la] - - g[la,lb] Zet[uc] PD[phi,lc]) + - e4phi Zet[uc] PD[gt[la,lb],lc], - + R[la,lb] -> + Rt[la,lb] + Rphi[la,lb] + + IfThen[conformalMethod==CMCCZ4, + + (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], trR -> gu[ua,ub] R[la,lb], (* K[la,lb] -> e4phi At[la,lb] + (1/3) g[la,lb] trK, *) @@ -1226,6 +1259,19 @@ keywordParameters = intParameters = { + { + Name -> conformalMethod, + Description -> "Treatment of conformal factor", + AllowedValues -> {{Value -> "0", Description -> "BSSN phi method"}, + {Value -> "1", Description -> "BSSN W method"}, + {Value -> "2", Description -> "CCZ4 method"}}, + Default -> 0 + }, + { + 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)", @@ -1236,28 +1282,31 @@ intParameters = Default -> 0 }, { - Name -> conformalMethod, - Description -> "Treatment of conformal factor", - AllowedValues -> {{Value -> "0", Description -> "phi method"}, - {Value -> "1", Description -> "W method"}}, + Name -> harmonicShift, + Description -> "Whether to use the harmonic shift", + AllowedValues -> {{Value -> "0", Description -> "Gamma driver shift"}, + {Value -> "1", Description -> "Harmonic shift"}}, 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 = { + { + Name -> dampk1, + Description -> "CCZ4 damping term 1 for Theta and Z", + Default -> 0 + }, + { + Name -> dampk2, + Description -> "CCZ4 damping term 2 for Theta and Z", + Default -> 0 + }, { Name -> LapseACoeff, Description -> "Whether to evolve A in time", @@ -1317,26 +1366,6 @@ realParameters = Description -> "Dissipation strength", AllowedValues -> {{Value -> "(0:*", Description -> "Positive"}}, Default -> 0 - }, - { - Name -> dampk1, - Description -> "Damping term 1 in Tet and Zet", - Default -> 0 - }, - { - Name -> dampk2, - Description -> "Damping term 2 in Tet and Zet", - Default -> 0 - }, - { - Name -> gamashift, - Description -> "Covariant shift term in Gamma", - Default -> 0 - }, - { - Name -> dampA, - Description -> "Damping term trA", - Default -> 0 } }; -- cgit v1.2.3