aboutsummaryrefslogtreecommitdiff
path: root/m
diff options
context:
space:
mode:
Diffstat (limited to 'm')
-rw-r--r--m/McLachlan_BSSN.m156
1 files changed, 120 insertions, 36 deletions
diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m
index 107ecd4..30cda03 100644
--- a/m/McLachlan_BSSN.m
+++ b/m/McLachlan_BSSN.m
@@ -110,7 +110,7 @@ Map [DefineTensor,
Rt, Rphi, gK,
T00, T0, T, rho, S,
x, y, z, r,
- epsdiss}];
+ epsdiss, Tet, Zet, dampk1k2}];
(* NOTE: It seems as if Lie[.,.] did not take these tensor weights
into account. Presumably, CD[.,.] and CDt[.,.] don't do this either. *)
@@ -191,7 +191,8 @@ evolvedGroups =
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" ],
+ SetGroupName [CreateGroupFromTensor [Tet ], prefix <> "Zenergy" ]};
evaluatedGroups =
{SetGroupName [CreateGroupFromTensor [H ], prefix <> "Ham"],
SetGroupName [CreateGroupFromTensor [M[la] ], prefix <> "mom"],
@@ -240,7 +241,8 @@ initialCalc =
alpha -> 1,
A -> 0,
beta[ua] -> 0,
- B[ua] -> 0
+ B[ua] -> 0,
+ Tet -> 0
}
};
@@ -367,7 +369,8 @@ initGammaCalc =
{
Xt[ua] -> 0,
A -> 0,
- B[ua] -> 0
+ B[ua] -> 0,
+ Tet -> 0
}
};
@@ -414,12 +417,13 @@ convertToADMBaseDtLapseShiftCalc =
(* see RHS *)
(*
admdtalpha -> - harmonicF alpha^harmonicN
- ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK)
+ ((1 - LapseAdvectionCoeff) A
+ + LapseAdvectionCoeff (trK - 2 Tet))
+ LapseAdvectionCoeff beta[ua] PDu[alpha,la],
*)
admdtalpha -> - harmonicF alpha^harmonicN
(+ LapseACoeff A
- + (1 - LapseACoeff) trK)
+ + (1 - LapseACoeff) (trK - 2 Tet))
+ LapseAdvectionCoeff Upwind[beta[ua], alpha, la],
admdtbeta[ua] -> IfThen[harmonicShift,
- 1/2 gtu[ua,uj] phi alpha
@@ -455,11 +459,12 @@ convertToADMBaseDtLapseShiftBoundaryCalc =
(* see RHS, but omit derivatives near the boundary *)
(*
admdtalpha -> - harmonicF alpha^harmonicN
- ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK),
+ ((1 - LapseAdvectionCoeff) A
+ + LapseAdvectionCoeff (trK - 2 Tet)),
*)
admdtalpha -> - harmonicF alpha^harmonicN
(+ LapseACoeff A
- + (1 - LapseACoeff) trK),
+ + (1 - LapseACoeff) (trK - 2 Tet)),
admdtbeta[ua] -> IfThen[harmonicShift,
0,
(* else *)
@@ -491,11 +496,12 @@ convertToADMBaseFakeDtLapseShiftCalc =
conditions) *)
(*
admdtalpha -> - harmonicF alpha^harmonicN
- ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK),
+ ((1 - LapseAdvectionCoeff) A
+ + LapseAdvectionCoeff (trK - 2 Tet)),
*)
admdtalpha -> - harmonicF alpha^harmonicN
(+ LapseACoeff A
- + (1 - LapseACoeff) trK),
+ + (1 - LapseACoeff) (trK - 2 Tet)),
admdtbeta[ua] -> IfThen[harmonicShift,
0,
(* else *)
@@ -528,7 +534,8 @@ evolCalc =
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]},
+ rho, S[la], trS, fac1, fac2, dottrK, dotXt[ua],
+ epsdiss[ua], dotTet, Zet[ua], dampk1k2},
Equations ->
{
dir[ua] -> Sign[beta[ua]],
@@ -546,7 +553,21 @@ 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]],
+ 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) *)
+ 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,
+
(* 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]
@@ -570,15 +591,13 @@ evolCalc =
Atm[ua,lb] -> gtu[ua,uc] At[lc,lb],
Atu[ua,ub] -> Atm[ua,lc] gtu[ub,uc],
- e4phi -> IfThen [conformalMethod, 1/phi^2, Exp[4 phi]],
- em4phi -> 1 / e4phi,
- g[la,lb] -> e4phi gt[la,lb],
- detg -> detgExpr,
- (* gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]], *)
- gu[ua,ub] -> em4phi gtu[ua,ub],
-
- R[la,lb] -> Rt[la,lb] + Rphi[la,lb],
-
+ (* 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],
+
(* Matter terms *)
(* rho = n^a n^b T_ab *)
@@ -599,11 +618,16 @@ evolCalc =
+ IfThen [conformalMethod, -1/3 phi, 1/6] PD[beta[ua],la],
(* PRD 62, 044034 (2000), eqn. (9) *)
- dot[gt[la,lb]] -> - 2 alpha At[la,lb]
+ (* 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])
+ 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 *)
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]
@@ -612,29 +636,52 @@ 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]
(* 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,
+
+ dot[Tet] -> dotTet,
+
(* 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 *)
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
(* 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) *)
+ (* Added Zet terms in the Ricci and Tet 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 At[la,lb] - 2 At[la,lc] Atm[uc,lb])
+ + alpha ((trK - 2 Tet) 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)),
@@ -643,18 +690,19 @@ evolCalc =
(* dot[alpha] -> - harmonicF alpha^harmonicN A + Lie[alpha, beta], *)
(*
dot[alpha] -> - harmonicF alpha^harmonicN (
- (1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK)
+ (1 - LapseAdvectionCoeff) A
+ + LapseAdvectionCoeff (trK - 2 Tet))
+ LapseAdvectionCoeff beta[ua] PDu[alpha,la],
- dot[A] -> (1 - LapseAdvectionCoeff) (dottrK - AlphaDriver A),
+ dot[A] -> (1 - LapseAdvectionCoeff) (dottrK - 2 dotTet - AlphaDriver A),
*)
dot[alpha] -> - harmonicF alpha^harmonicN
(+ LapseACoeff A
- + (1 - LapseACoeff) trK),
+ + (1 - LapseACoeff) (trK - 2 Tet)),
- dot[A] -> + LapseACoeff (dottrK - AlphaDriver A),
+ dot[A] -> + LapseACoeff (dottrK - 2 dotTet - AlphaDriver A),
eta -> etaExpr,
@@ -702,6 +750,8 @@ advectCalc =
dot[Xt[ui]] -> dot[Xt[ui]] + Upwind[beta[uj], Xt[ui], lj],
+ dot[Tet] -> dot[Tet] + Upwind[beta[ua], Tet, 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],
@@ -724,6 +774,7 @@ evolCalc1 = PartialCalculation[evolCalc, "1",
ConditionalOnKeyword -> {"RHS_calculation", "split"}
},
{
+ dot[Tet],
dot[trK],
dot[phi],
dot[gt[la,lb]],
@@ -755,7 +806,7 @@ dissCalc =
epsdiss[ua] -> EpsDiss,
Sequence@@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], trK, Tet, At[la,lb], alpha, A, beta[ua], B[ua]}}]
}
};
@@ -787,6 +838,7 @@ RHSStaticBoundaryCalc =
{
dot[phi] -> 0,
dot[gt[la,lb]] -> 0,
+ dot[Tet] -> 0,
dot[trK] -> 0,
dot[At[la,lb]] -> 0,
dot[Xt[ua]] -> 0,
@@ -809,6 +861,7 @@ initRHSCalc =
{
dot[phi] -> 0,
dot[gt[la,lb]] -> 0,
+ dot[Tet] -> 0,
dot[trK] -> 0,
dot[At[la,lb]] -> 0,
dot[Xt[ua]] -> 0,
@@ -848,6 +901,7 @@ 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],
@@ -899,6 +953,7 @@ boundaryCalc =
{
phi -> IfThen [conformalMethod, 1, 0],
gt[la,lb] -> KD[la,lb],
+ Tet -> 0,
trK -> 0,
At[la,lb] -> 0,
Xt[ua] -> 0,
@@ -925,7 +980,7 @@ constraintsCalc =
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},
+ rho, S[la], fac1, fac2, Zet[ua]},
Equations ->
{
detgt -> 1 (* detgtExpr *),
@@ -942,7 +997,17 @@ 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]],
+ 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]),
+
(* 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]
@@ -985,13 +1050,6 @@ constraintsCalc =
+ 4 cdphi[li] cdphi[lj]
- 4 gt[li,lj] gtu[ul,un] cdphi[ln] cdphi[ll],
- e4phi -> IfThen [conformalMethod, 1/phi^2, Exp[4 phi]],
- em4phi -> 1 / e4phi,
- g[la,lb] -> e4phi gt[la,lb],
- (* detg -> detgExpr, *)
- (* gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]], *)
- detg -> e4phi^3,
- gu[ua,ub] -> em4phi gtu[ua,ub],
(* 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 *)
@@ -999,7 +1057,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]),
- R[la,lb] -> Rt[la,lb] + Rphi[la,lb],
+ (* 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],
+
trR -> gu[ua,ub] R[la,lb],
(* K[la,lb] -> e4phi At[la,lb] + (1/3) g[la,lb] trK, *)
@@ -1252,6 +1316,26 @@ 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
}
};
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-29 09:52:24 +0000