diff options
| author | Erik Schnetter <schnetter@cct.lsu.edu> | 2010-01-09 12:25:37 -0600 |
|---|---|---|
| committer | Erik Schnetter <schnetter@cct.lsu.edu> | 2010-01-09 12:25:37 -0600 |
| commit | e58418f219164cd87eaf03076a2671f0a1ea2dd4 (patch) | |
| tree | ac9b65d1b95ebaded0bdba88810663812c1c6d88 /m/McLachlan_BSSN_Peter.m | |
| parent | 5efd84e6066707cffd0b277802213a19cb806ff4 (diff) | |
Regenerate code
Diffstat (limited to 'm/McLachlan_BSSN_Peter.m')
| -rw-r--r-- | m/McLachlan_BSSN_Peter.m | 960 |
1 files changed, 960 insertions, 0 deletions
diff --git a/m/McLachlan_BSSN_Peter.m b/m/McLachlan_BSSN_Peter.m new file mode 100644 index 0000000..d72f5f0 --- /dev/null +++ b/m/McLachlan_BSSN_Peter.m @@ -0,0 +1,960 @@ +$Path = Join[$Path, {"~/Calpha/kranc/Tools/CodeGen",
+ "~/Calpha/kranc/Tools/MathematicaMisc"}];
+
+Get["KrancThorn`"];
+
+SetEnhancedTimes[False];
+SetSourceLanguage["C"];
+
+(******************************************************************************)
+(* Options *)
+(******************************************************************************)
+
+createCode[derivOrder_, useGlobalDerivs_, evolutionTimelevels_, addMatter_] :=
+Module[{},
+
+prefix = "ML_";
+suffix =
+ If [useGlobalDerivs, "_MP", ""] <>
+ If [derivOrder!=4, "_O" <> ToString[derivOrder], ""] <>
+ If [evolutionTimelevels!=3, "_TL" <> ToString[evolutionTimelevels], ""] <>
+ If [addMatter==1, "_M", ""];
+
+BSSN = prefix <> "BSSN" <> suffix;
+
+(******************************************************************************)
+(* Derivatives *)
+(******************************************************************************)
+
+KD = KroneckerDelta;
+
+derivatives =
+{
+ PDstandardNth[i_] -> StandardCenteredDifferenceOperator[1,derivOrder/2,i],
+ PDstandardNth[i_, i_] -> StandardCenteredDifferenceOperator[2,derivOrder/2,i],
+ PDstandardNth[i_, j_] -> StandardCenteredDifferenceOperator[1,derivOrder/2,i]
+ StandardCenteredDifferenceOperator[1,derivOrder/2,j],
+
+(* PD: These comes from my mathematica notebook: Upwind-Kranc-Convert.nb that
+converts upwinding finite differencing operators generated by
+StandardUpwindDifferenceOperator into this form *)
+
+ Switch[derivOrder,2,
+ PDupwindNth[1] -> (dir[1]*(-3 + 4*shift[1]^dir[1] -
+ shift[1]^(2*dir[1])))/(2*spacing[1]),
+ 4,
+ PDupwindNth[1] -> (dir[1]*(-10 - 3/shift[1]^dir[1] + 18*shift[1]^dir[1] -
+ 6*shift[1]^(2*dir[1]) +
+ shift[1]^(3*dir[1])))/(12*spacing[1]),
+ 6,
+ PDupwindNth[1] -> (dir[1]*(-35 + 2/shift[1]^(2*dir[1]) -
+ 24/shift[1]^dir[1] + 80*shift[1]^dir[1] -
+ 30*shift[1]^(2*dir[1]) + 8*shift[1]^(3*dir[1]) -
+ shift[1]^(4*dir[1])))/(60*spacing[1]),
+ 8,
+ PDupwindNth[1] -> (dir[1]*(-378 - 5/shift[1]^(3*dir[1]) +
+ 60/shift[1]^(2*dir[1]) - 420/shift[1]^dir[1] +
+ 1050*shift[1]^dir[1] - 420*shift[1]^(2*dir[1]) +
+ 140*shift[1]^(3*dir[1]) - 30*shift[1]^(4*dir[1]) +
+ 3*shift[1]^(5*dir[1])))/(840*spacing[1])],
+ Switch[derivOrder,2,
+ PDupwindNth[2] -> (dir[2]*(-3 + 4*shift[2]^dir[2] -
+ shift[2]^(2*dir[2])))/(2*spacing[2]),
+ 4,
+ PDupwindNth[2] -> (dir[2]*(-10 - 3/shift[2]^dir[2] + 18*shift[2]^dir[2] -
+ 6*shift[2]^(2*dir[2]) +
+ shift[2]^(3*dir[2])))/(12*spacing[2]),
+ 6,
+ PDupwindNth[2] -> (dir[2]*(-35 + 2/shift[2]^(2*dir[2]) -
+ 24/shift[2]^dir[2] + 80*shift[2]^dir[2] -
+ 30*shift[2]^(2*dir[2]) + 8*shift[2]^(3*dir[2]) -
+ shift[2]^(4*dir[2])))/(60*spacing[2]),
+ 8,
+ PDupwindNth[2] -> (dir[2]*(-378 - 5/shift[2]^(3*dir[2]) +
+ 60/shift[2]^(2*dir[2]) - 420/shift[2]^dir[2] +
+ 1050*shift[2]^dir[2] - 420*shift[2]^(2*dir[2]) +
+ 140*shift[2]^(3*dir[2]) - 30*shift[2]^(4*dir[2]) +
+ 3*shift[2]^(5*dir[2])))/(840*spacing[2])],
+ Switch[derivOrder,2,
+ PDupwindNth[3] -> (dir[3]*(-3 + 4*shift[3]^dir[3] -
+ shift[3]^(2*dir[3])))/(2*spacing[3]),
+ 4,
+ PDupwindNth[3] -> (dir[3]*(-10 - 3/shift[3]^dir[3] + 18*shift[3]^dir[3] -
+ 6*shift[3]^(2*dir[3]) +
+ shift[3]^(3*dir[3])))/(12*spacing[3]),
+ 6,
+ PDupwindNth[3] -> (dir[3]*(-35 + 2/shift[3]^(2*dir[3]) -
+ 24/shift[3]^dir[3] + 80*shift[3]^dir[3] -
+ 30*shift[3]^(2*dir[3]) + 8*shift[3]^(3*dir[3]) -
+ shift[3]^(4*dir[3])))/(60*spacing[3]),
+ 8,
+ PDupwindNth[3] -> (dir[3]*(-378 - 5/shift[3]^(3*dir[3]) +
+ 60/shift[3]^(2*dir[3]) - 420/shift[3]^dir[3] +
+ 1050*shift[3]^dir[3] - 420*shift[3]^(2*dir[3]) +
+ 140*shift[3]^(3*dir[3]) - 30*shift[3]^(4*dir[3]) +
+ 3*shift[3]^(5*dir[3])))/(840*spacing[3])],
+
+(* PDupwindpNth[i_] -> StandardUpwindDifferenceOperator[1,derivOrder/2-1,derivOrder/2+1,i],
+ PDupwindmNth[i_] -> StandardUpwindDifferenceOperator[1,derivOrder/2+1,derivOrder/2-1,i], *)
+
+(*
+ PDPlus [i_] -> (+1) (-1 + shift[i]^(+1)) / spacing[i],
+ PDMinus [i_] -> (-1) (-1 + shift[i]^(-1)) / spacing[i]
+*)
+ PDPlus [i_] -> StandardUpwindDifferenceOperator[1,0,derivOrder/2,i],
+ PDMinus [i_] -> StandardUpwindDifferenceOperator[1,derivOrder/2,0,i]
+};
+
+FD = PDstandardNth;
+FDu = PDupwindNth;
+(*FDp = PDupwindpNth;
+FDm = PDupwindmNth; *)
+
+ResetJacobians;
+If [useGlobalDerivs,
+ DefineJacobian[PD, FD, J, dJ],
+ DefineJacobian[PD, FD, KD, Zero3]];
+If [useGlobalDerivs,
+ DefineJacobian[PDu, FDu, J, dJ],
+ DefineJacobian[PDu, FDu, KD, Zero3]];
+(*If [useGlobalDerivs,
+ DefineJacobian[PDp, FDp, J, dJ],
+ DefineJacobian[PDp, FDp, KD, Zero3]];
+If [useGlobalDerivs,
+ DefineJacobian[PDm, FDm, J, dJ],
+ DefineJacobian[PDm, FDm, KD, Zero3]]; *)
+
+(*PD = FD;
+PDu = FDu; *)
+(*PDp = FDp;
+PDm = FDm; *)
+
+
+
+(******************************************************************************)
+(* 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,
+ J, dJ,
+ g, K, alpha, beta, H, M, detg, gu, G, R, trR, Km, trK,
+ phi, gt, At, Xt, Xtn, A, B, Atm, Atu, trA, Ats, trAts, cXt, cS, cA,
+ e4phi, em4phi, ddetg, detgt, gtu, ddetgt, dgtu, ddgtu, Gt, Rt, Rphi, gK,
+ T00, T0, T, rho, S,
+ Psi0re, Psi0im, Psi1re, Psi1im, Psi2re, Psi2im, Psi3re, Psi3im,
+ Psi4re, Psi4im,
+ er, eth, eph, mm1A, mm1L, mm1, mm2A, mm2B, mm2L, mm2,
+ ssA, ssB, ssC, ssL, ss, ss0, tt, ss0, kk, nn, kk0, nn0, mmre, mmim,
+ EE, BB}];
+
+(* 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,
+ {g[la,lb], K[la,lb], R[la,lb],
+ gt[la,lb], At[la,lb], Ats[la,lb], Rt[la,lb], Rphi[la,lb], T[la,lb]}];
+AssertSymmetricIncreasing [dJ[ua,lb,lc], lb, lc];
+AssertSymmetricIncreasing [G[ua,lb,lc], lb, lc];
+AssertSymmetricIncreasing [Gt[ua,lb,lc], lb, lc];
+AssertSymmetricIncreasing [gK[la,lb,lc], la, lb];
+Map [AssertSymmetricDecreasing, {gu[ua,ub], gtu[ua,ub], Atu[ua,ub]}];
+AssertSymmetricDecreasing [dgtu[ua,ub,lc], ua, ub];
+AssertSymmetricDecreasing [ddgtu[ua,ub,lc,ld], ua, ub];
+AssertSymmetricIncreasing [ddgtu[ua,ub,lc,ld], lc, ld];
+
+DefineConnection [CD, PD, G];
+DefineConnection [CDt, PD, Gt];
+
+Map [DefineTensor,
+ {gxx, gxy, gxz, gyy, gyz, gzz,
+ kxx, kxy, kxz, kyy, kyz, kzz,
+ alp,
+ dtalp,
+ betax, betay, betaz,
+ dtbetax, dtbetay, dtbetaz,
+ eTtt,
+ eTtx, eTty, eTtz,
+ eTxx, eTxy, eTxz, eTyy, eTyz, eTzz}];
+
+(******************************************************************************)
+(* Expressions *)
+(******************************************************************************)
+
+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]]]]}];
+
+pi = N[Pi,40];
+
+(******************************************************************************)
+(* Groups *)
+(******************************************************************************)
+
+SetGroupTimelevels[g_,tl_] = Join[g, {Timelevels -> tl}];
+
+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" ]};
+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 =
+ {{"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}},
+ {"Coordinates::jacobian", {J11, J12, J13, J21, J22, J23, J31, J32, J33}},
+ {"Coordinates::jacobian2", {dJ111, dJ112, dJ113, dJ122, dJ123, dJ133,
+ dJ211, dJ212, dJ213, dJ222, dJ223, dJ233,
+ dJ311, dJ312, dJ313, dJ322, dJ323, dJ333}}
+};
+
+
+
+groups = Join [declaredGroups, extraGroups];
+
+(******************************************************************************)
+(* Initial data *)
+(******************************************************************************)
+
+initialCalc =
+{
+ Name -> BSSN <> "_Minkowski",
+ Schedule -> {"IN ADMBase_InitialData"},
+ ConditionalOnKeyword -> {"my_initial_data", "Minkowski"},
+ Equations ->
+ {
+ phi -> 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
+ }
+};
+
+(******************************************************************************)
+(* Convert from ADMBase *)
+(******************************************************************************)
+
+convertFromADMBaseCalc =
+{
+ Name -> BSSN <> "_convertFromADMBase",
+ Schedule -> {"AT initial AFTER ADMBase_PostInitial"},
+ ConditionalOnKeyword -> {"my_initial_data", "ADMBase"},
+ Shorthands -> {g[la,lb], detg, gu[ua,ub], em4phi, K[la,lb]},
+ Equations ->
+ {
+ g11 -> gxx,
+ g12 -> gxy,
+ g13 -> gxz,
+ g22 -> gyy,
+ g23 -> gyz,
+ g33 -> gzz,
+
+ detg -> detgExpr,
+ gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]],
+
+ phi -> Log [detg] / 12,
+ em4phi -> Exp [-4 phi],
+ gt[la,lb] -> em4phi g[la,lb],
+
+ K11 -> kxx,
+ K12 -> kxy,
+ K13 -> kxz,
+ K22 -> kyy,
+ K23 -> kyz,
+ K33 -> kzz,
+
+ trK -> gu[ua,ub] K[la,lb],
+ At[la,lb] -> em4phi (K[la,lb] - (1/3) g[la,lb] trK),
+
+ alpha -> alp,
+
+ beta1 -> betax,
+ beta2 -> betay,
+ beta3 -> betaz
+ }
+};
+
+convertFromADMBaseGammaCalc =
+{
+ Name -> BSSN <> "_convertFromADMBaseGamma",
+ Schedule -> {"AT initial AFTER " <> BSSN <> "_convertFromADMBase"},
+ ConditionalOnKeyword -> {"my_initial_data", "ADMBase"},
+ Where -> Interior,
+ Shorthands -> {detgt, gtu[ua,ub], Gt[ua,lb,lc], dir[ua]},
+ 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 -> - dtalp / (harmonicF alpha^harmonicN) (LapseAdvectionCoeff - 1),
+
+ B1 -> 1/ShiftGammaCoeff
+ (dtbetax - ShiftAdvectionCoeff beta[ua] PDu[beta1,la]),
+ B2 -> 1/ShiftGammaCoeff
+ (dtbetay - ShiftAdvectionCoeff beta[ua] PDu[beta2,la]),
+ B3 -> 1/ShiftGammaCoeff
+ (dtbetaz - ShiftAdvectionCoeff beta[ua] PDu[beta3,la])
+ }
+};
+
+(******************************************************************************)
+(* Convert to ADMBase *)
+(******************************************************************************)
+
+convertToADMBaseCalc =
+{
+ Name -> BSSN <> "_convertToADMBase",
+ Schedule -> {"IN ML_BSSN_convertToADMBaseGroup"},
+ Where -> Interior,
+ Shorthands -> {e4phi, g[la,lb], K[la,lb], dir[ua]},
+ Equations ->
+ {
+ dir[ua] -> Sign[beta[ua]],
+ e4phi -> Exp [4 phi],
+ g[la,lb] -> e4phi gt[la,lb],
+ gxx -> g11,
+ gxy -> g12,
+ gxz -> g13,
+ gyy -> g22,
+ gyz -> g23,
+ gzz -> g33,
+ K[la,lb] -> e4phi At[la,lb] + (1/3) g[la,lb] trK,
+ kxx -> K11,
+ kxy -> K12,
+ kxz -> K13,
+ kyy -> K22,
+ kyz -> K23,
+ kzz -> K33,
+ alp -> alpha,
+ betax -> beta1,
+ betay -> beta2,
+ betaz -> beta3,
+ (* see RHS *)
+ dtalp -> - harmonicF alpha^harmonicN
+ ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK)
+ + LapseAdvectionCoeff beta[ua] PDu[alpha,la],
+ dtbetax -> + ShiftGammaCoeff B1
+ + ShiftAdvectionCoeff beta[ub] PDu[beta1,lb],
+ dtbetay -> + ShiftGammaCoeff B2
+ + ShiftAdvectionCoeff beta[ub] PDu[beta2,lb],
+ dtbetaz -> + ShiftGammaCoeff B3
+ + ShiftAdvectionCoeff beta[ub] PDu[beta3,lb]
+
+ }
+};
+
+boundaryCalcADMBase =
+{
+ Name -> BSSN <> "_ADMBaseBoundary",
+ Schedule -> {"IN ML_BSSN_convertToADMBaseGroup AFTER " <> BSSN <> "_convertToADMBase"},
+ ConditionalOnKeyword -> {"my_boundary_condition", "Minkowski"},
+ Where -> BoundaryWithGhosts,
+ Equations ->
+ {
+ gxx -> 1,
+ gxy -> 0,
+ gxz -> 0,
+ gyy -> 1,
+ gyz -> 0,
+ gzz -> 1,
+ kxx -> 0,
+ kxy -> 0,
+ kxz -> 0,
+ kyy -> 0,
+ kyz -> 0,
+ kzz -> 0,
+ alp -> 1,
+ dtalp -> 0,
+ betax -> 0,
+ betay -> 0,
+ betaz -> 0,
+ dtbetax -> 0,
+ dtbetay -> 0,
+ dtbetaz -> 0
+ }
+};
+
+(******************************************************************************)
+(* Evolution equations *)
+(******************************************************************************)
+
+evolCalc =
+{
+ Name -> BSSN <> "_RHS",
+ Schedule -> {"IN " <> BSSN <> "_evolCalcGroup"},
+ Where -> Interior,
+ Shorthands -> {detgt, ddetgt[la], gtu[ua,ub],
+ dgtu[ua,ub,lc], ddgtu[ua,ub,lc,ld], Gt[ua,lb,lc],
+ Xtn[ua], Rt[la,lb], Rphi[la,lb], R[la,lb],
+ Atm[ua,lb], Atu[ua,ub],
+ e4phi, em4phi, g[la,lb], detg,
+ ddetg[la], gu[ua,ub], G[ua,lb,lc], Ats[la,lb], trAts,
+ T00, T0[la], T[la,lb], rho, S[la], trS, dir[ua]},
+ Equations ->
+ {
+ dir[ua] -> Sign[beta[ua]],
+ detgt -> 1 (* detgtExpr *),
+ ddetgt[la] -> 0 (* ddetgtExpr[la] *),
+
+ (* This leads to simpler code... *)
+ gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
+ dgtu[ua,ub,lc] -> - gtu[ua,ud] gtu[ub,ue] PD[gt[ld,le],lc],
+ ddgtu[ua,ub,lc,ld] -> - dgtu[ua,ue,ld] gtu[ub,uf] PD[gt[le,lf],lc]
+ - gtu[ua,ue] dgtu[ub,uf,ld] PD[gt[le,lf],lc]
+ - gtu[ua,ue] gtu[ub,uf] PD[gt[le,lf],lc,ld],
+ Gt[ua,lb,lc] -> 1/2 gtu[ua,ud]
+ (PD[gt[lb,ld],lc] + PD[gt[lc,ld],lb] - PD[gt[lb,lc],ld]),
+
+ (* 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],
+
+ (* 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] gt[li,ln] Gt[un,lj,lk]
+ + (1/2) Xtn[uk] gt[lj,ln] Gt[un,li,lk]
+ + gtu[ul,um] (+ Gt[uk,ll,li] gt[lj,ln] Gt[un,lk,lm]
+ + Gt[uk,ll,lj] gt[li,ln] Gt[un,lk,lm]
+ + Gt[uk,li,lm] gt[lk,ln] Gt[un,ll,lj]),
+ (* PRD 62, 044034 (2000), eqn. (15) *)
+ Rphi[li,lj] -> - 2 CDt[phi,lj,li]
+ - 2 gt[li,lj] gtu[ul,un] CDt[phi,ll,ln]
+ + 4 CDt[phi,li] CDt[phi,lj]
+ - 4 gt[li,lj] gtu[ul,un] CDt[phi,ln] CDt[phi,ll],
+
+ Atm[ua,lb] -> gtu[ua,uc] At[lc,lb],
+ Atu[ua,ub] -> Atm[ua,lc] gtu[ub,uc],
+
+ e4phi -> 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],
+ G[ua,lb,lc] -> Gt[ua,lb,lc]
+ + 2 (KD[ua,lb] PD[phi,lc] + KD[ua,lc] PD[phi,lb]
+ - gtu[ua,ud] gt[lb,lc] PD[phi,ld]),
+
+ R[la,lb] -> Rt[la,lb] + Rphi[la,lb],
+
+ (* Matter terms *)
+
+ T00 -> addMatter eTtt,
+ T01 -> addMatter eTtx,
+ T02 -> addMatter eTty,
+ T03 -> addMatter eTtz,
+ T11 -> addMatter eTxx,
+ T12 -> addMatter eTxy,
+ T13 -> addMatter eTxz,
+ T22 -> addMatter eTyy,
+ T23 -> addMatter eTyz,
+ T33 -> addMatter eTzz,
+
+ (* 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 (gu[ui,uj] T[li,lj]),
+
+ (* RHS terms *)
+
+ (* PRD 62, 044034 (2000), eqn. (10) *)
+ (* PRD 67 084023 (2003), eqn. (16) and (23) *)
+ dot[phi] -> - (1/6) alpha trK
+ + beta[ua] PDu[phi,la]
+ + (1/6) PD[beta[ua],la],
+
+ (* PRD 62, 044034 (2000), eqn. (9) *)
+ dot[gt[la,lb]] -> - 2 alpha At[la,lb]
+ + beta[uc] PDu[gt[la,lb],lc]
+ + 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) *)
+ dot[Xt[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] PD[phi,lj])
+ + gtu[uj,ul] PD[beta[ui],lj,ll]
+ + (1/3) gtu[ui,uj] PD[beta[ul],lj,ll]
+ + beta[uj] PDu[Xt[ui],lj]
+ - Xtn[uj] PD[beta[ui],lj]
+ + (2/3) Xtn[ui] PD[beta[uj],lj]
+ (* Equation (4.28) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *)
+ + addMatter (- 16 pi alpha gtu[ui,uj] S[lj]),
+
+ (* PRD 62, 044034 (2000), eqn. (11) *)
+ dot[trK] -> - gu[ua,ub] CD[alpha,la,lb]
+ + alpha (Atm[ua,lb] Atm[ub,la] + (1/3) trK^2)
+ + beta[ua] PDu[trK,la]
+ (* Equation (4.21) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *)
+ + addMatter (4 pi alpha (rho + trS)),
+
+ (* PRD 62, 044034 (2000), eqn. (12) *)
+ (* TODO: use Hamiltonian constraint to make tracefree *)
+ Ats[la,lb] -> - CD[alpha,la,lb] + 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])
+ + beta[uc] PDu[At[la,lb],lc]
+ + 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) (dot[trK] - AlphaDriver A),
+ (* dot[beta[ua]] -> eta Xt[ua], *)
+ (* dot[beta[ua]] -> ShiftGammaCoeff alpha^ShiftAlphaPower B[ua], *)
+
+ dot[beta[ua]] -> + ShiftGammaCoeff B[ua]
+ + ShiftAdvectionCoeff beta[ub] PDu[beta[ua],lb],
+
+ dot[B[ua]] -> + dot[Xt[ua]] - BetaDriver B[ua]
+ + ShiftAdvectionCoeff beta[ub] ( PDu[B[ua],lb]
+ - PDu[Xt[ua],lb] )
+ }
+};
+
+RHSBoundaryCalc =
+{
+ Name -> BSSN <> "_RHSBoundary",
+ Schedule -> {"IN MoL_RHSBoundaries"},
+ ConditionalOnKeyword -> {"my_rhs_boundary_condition", "radiative"},
+ Where -> Boundary,
+ Shorthands -> {detgt, gtu[ua,ub], em4phi, gu[ua,ub],
+ nn[la], nu[ua], nlen, nlen2, su[ua],
+ vg},
+ Equations ->
+ {
+ detgt -> 1 (* detgtExpr *),
+ gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
+ em4phi -> Exp [-4 phi],
+ gu[ua,ub] -> em4phi gtu[ua,ub],
+
+ nn[la] -> normal[la],
+ 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] PDMinus[phi ,lc],
+ dot[gt[la,lb]] -> - su[uc] PDMinus[gt[la,lb],lc],
+ dot[trK] -> - vg su[uc] PDMinus[trK ,lc],
+ dot[At[la,lb]] -> - su[uc] PDMinus[At[la,lb],lc],
+ dot[Xt[ua]] -> - su[uc] PDMinus[Xt[ua] ,lc],
+ dot[alpha] -> - vg su[uc] PDMinus[alpha ,lc],
+ dot[A] -> - vg su[uc] PDMinus[A ,lc],
+ dot[beta[ua]] -> - su[uc] PDMinus[beta[ua] ,lc],
+ dot[B[ua]] -> - su[uc] PDMinus[B[ua] ,lc]
+ }
+};
+
+enforceCalc =
+{
+ Name -> BSSN <> "_enforce",
+ Schedule -> {"IN MoL_PostStep BEFORE " <> BSSN <> "_BoundConds"},
+ Shorthands -> {detgt, gtu[ua,ub], trAt},
+ Equations ->
+ {
+ 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, 10^(-10)] *)
+ }
+};
+
+(******************************************************************************)
+(* Boundary conditions *)
+(******************************************************************************)
+
+boundaryCalc =
+{
+ Name -> BSSN <> "_boundary",
+ Schedule -> {"IN MoL_PostStep"},
+ ConditionalOnKeyword -> {"my_boundary_condition", "Minkowski"},
+ Where -> BoundaryWithGhosts,
+ Equations ->
+ {
+ phi -> 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
+ }
+};
+
+(******************************************************************************)
+(* Constraint equations *)
+(******************************************************************************)
+
+constraintsCalc =
+{
+ Name -> BSSN <> "_constraints",
+ Schedule -> {"IN " <> BSSN <> "_constraintsCalcGroup"},
+ Where -> Interior,
+ Shorthands -> {detgt, ddetgt[la], gtu[ua,ub], Gt[ua,lb,lc], 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[la,lb],
+ gK[la,lb,lc],
+ T00, T0[la], T[la,lb], rho, S[la]},
+ Equations ->
+ {
+ detgt -> 1 (* detgtExpr *),
+ ddetgt[la] -> 0 (* ddetgtExpr[la] *),
+ 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]),
+
+ (* PRD 62, 044034 (2000), eqn. (18) *)
+ (* Note: This differs from the Goddard formulation,
+ which is e.g. described in PRD 70 (2004) 124025, eqn. (6).
+ Goddard has a Gamma^k replaced by its definition in terms
+ of metric derivatives. *)
+ 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) Xt[uk] gt[li,ln] Gt[un,lj,lk]
+ + (1/2) Xt[uk] gt[lj,ln] Gt[un,li,lk]
+ + gtu[ul,um] (+ Gt[uk,ll,li] gt[lj,ln] Gt[un,lk,lm]
+ + Gt[uk,ll,lj] gt[li,ln] Gt[un,lk,lm]
+ + Gt[uk,li,lm] gt[lk,ln] Gt[un,ll,lj]),
+
+ (* 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]),
+*)
+ (* PRD 62, 044034 (2000), eqn. (15) *)
+ Rphi[li,lj] -> - 2 CDt[phi,lj,li]
+ - 2 gt[li,lj] gtu[ul,un] CDt[phi,ll,ln]
+ + 4 CDt[phi,li] CDt[phi,lj]
+ - 4 gt[li,lj] gtu[ul,un] CDt[phi,ln] CDt[phi,ll],
+
+ e4phi -> 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 *)
+ 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],
+ 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 *)
+
+ T00 -> eTtt,
+ T01 -> eTtx,
+ T02 -> eTty,
+ T03 -> eTtz,
+ T11 -> eTxx,
+ T12 -> eTxy,
+ T13 -> eTxz,
+ T22 -> eTyy,
+ T23 -> eTyz,
+ T33 -> eTzz,
+
+ (* 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] PD[phi,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]
+ }
+};
+
+constraintsBoundaryCalc =
+{
+ Name -> BSSN <> "_constraints_boundary",
+ Schedule -> {"IN " <> BSSN <> "_constraintsCalcGroup AFTER " <> BSSN <> "_constraints"},
+ Where -> BoundaryWithGhosts,
+ Equations ->
+ {
+ H -> 0,
+ M[la] -> 0
+ }
+};
+
+(******************************************************************************)
+(* Implementations *)
+(******************************************************************************)
+
+inheritedImplementations =
+ Join[{"ADMBase"},
+ If [addMatter!=0, {"TmunuBase"}, {}],
+ If [useGlobalDerivs, {"Coordinates"}, {}]];
+
+(******************************************************************************)
+(* Parameters *)
+(******************************************************************************)
+
+inheritedKeywordParameters = {};
+
+extendedKeywordParameters =
+{
+ {
+ Name -> "ADMBase::evolution_method",
+ AllowedValues -> {BSSN}
+ },
+ {
+ Name -> "ADMBase::lapse_evolution_method",
+ AllowedValues -> {BSSN}
+ },
+ {
+ Name -> "ADMBase::shift_evolution_method",
+ AllowedValues -> {BSSN}
+ }
+};
+
+keywordParameters =
+{
+ {
+ Name -> "my_initial_data",
+ (* Visibility -> "restricted", *)
+ (* Description -> "ddd", *)
+ AllowedValues -> {"ADMBase", "Minkowski"},
+ Default -> "ADMBase"
+ },
+ {
+ Name -> "my_rhs_boundary_condition",
+ (* Visibility -> "restricted", *)
+ (* Description -> "ddd", *)
+ AllowedValues -> {"none", "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"
+ }
+};
+
+intParameters =
+{
+ {
+ Name -> harmonicN,
+ Description -> "d/dt alpha = - f alpha^n K (harmonic=2, 1+log=1)",
+ Default -> 2
+ },
+ {
+ Name -> ShiftAlphaPower,
+ Default -> 0
+ }
+};
+
+realParameters =
+{
+ {
+ Name -> harmonicF,
+ Description -> "d/dt alpha = - f alpha^n K (harmonic=1, 1+log=2)",
+ Default -> 1
+ },
+ {
+ Name -> AlphaDriver,
+ Default -> 0
+ },
+ {
+ Name -> ShiftGammaCoeff,
+ Default -> 0
+ },
+ {
+ Name -> BetaDriver,
+ Default -> 0
+ },
+ {
+ Name -> LapseAdvectionCoeff,
+ Description -> "Factor in front of the shift advection terms in 1+log",
+ Default -> 1
+ },
+ {
+ Name -> ShiftAdvectionCoeff,
+ Description -> "Factor in front of the shift advection terms in gamma driver",
+ Default -> 1
+ }
+};
+
+
+(******************************************************************************)
+(* Construct the thorns *)
+(******************************************************************************)
+
+calculations =
+{
+ initialCalc,
+ convertFromADMBaseCalc,
+ convertFromADMBaseGammaCalc,
+ evolCalc,
+ RHSBoundaryCalc,
+ enforceCalc,
+ boundaryCalc,
+ convertToADMBaseCalc,
+ boundaryCalcADMBase,
+ constraintsCalc,
+ constraintsBoundaryCalc
+};
+
+CreateKrancThornTT [groups, ".", BSSN,
+ Calculations -> calculations,
+ DeclaredGroups -> declaredGroupNames,
+ PartialDerivatives -> derivatives,
+ EvolutionTimelevels -> evolutionTimelevels,
+ UseLoopControl -> True,
+ InheritedImplementations -> inheritedImplementations,
+ InheritedKeywordParameters -> inheritedKeywordParameters,
+ ExtendedKeywordParameters -> extendedKeywordParameters,
+ KeywordParameters -> keywordParameters,
+ IntParameters -> intParameters,
+ RealParameters -> realParameters
+];
+
+];
+
+
+
+(******************************************************************************)
+(* Options *)
+(******************************************************************************)
+
+(* derivative order: 2, 4, 6, 8, ... *)
+(* useGlobalDerivs: 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) *)
+
+createCode[4, False, 3, 0];
+(*createCode[6, False, 3, 0];
+createCode[8, False, 3, 0];
+createCode[4, False, 3, 1];
+createCode[4, True, 3, 0]; *)
|