1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
|
SetEnhancedTimes[False];
SetSourceLanguage["C"];
prefix = "ML_";
thorn = prefix <> "Kretschmann";
Map [DefineTensor, {
g, gu, G, R, Km, trK, phi, gt, At, e4phi, Dphi, D2phi,
Gl, Riem, Riemu, K, Ku, CDK, CDKu, Kretsch, Dg, D2g
}
];
declaredGroups = {
SetGroupName[CreateGroupFromTensor[Kretsch], prefix <> "Kretschmann"]
};
declaredGroupNames = Map[First, declaredGroups];
extraGroups = {
{ "ML_BSSN::ML_metric", { gt11, gt12, gt13, gt22, gt23, gt33 } },
{ "ML_BSSN::ML_curv", { At11, At12, At13, At22, At23, At33 } },
{ "ML_BSSN::ML_trace_curv", { trK } },
{ "ML_BSSN::ML_log_confac", { phi } }
};
groups = Join[declaredGroups, extraGroups];
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]
};
PD = PDstandardNth;
Map[AssertSymmetricIncreasing,
{ g[la, lb], K[la, lb], R[la, lb],
gt[la, lb], At[la, lb], D2phi[la, lb]}];
Map [AssertSymmetricIncreasing,
{ gu[ua, ub], gtu[ua, ub], Ku[ua, ub]}];
AssertSymmetricIncreasing [G[ua,lb,lc], lb, lc];
AssertSymmetricIncreasing [Gl[la,lb,lc], lb, lc];
AssertSymmetricIncreasing [Dg[la, lb,lc], la, lb];
AssertSymmetricIncreasing [CDK[la, lb, lc], la, lb];
AssertSymmetricIncreasing [CDKu[ua, ub, uc], ua, ub];
AssertSymmetricIncreasing [D2g[la, lb, lc, ld], la, lb];
AssertSymmetricIncreasing [D2g[la, lb, lc, ld], lc, ld];
kretschmannCalc =
{
Name -> thorn <> "_kretschmann",
Schedule -> Automatic,
Where -> Interior,
After -> "MoL_PostStep",
Shorthands -> {term1, term2, term3,
e4phi, Dphi[la], D2phi[la, lb],
Dg[la, lb, lc], D2g[la, lb, lc, ld],
g[la, lb], gu[ua, ub],
Gl[la, lb, lc], G[ua, lb, lc],
Riem[la, lb, lc, ld], Riemu[ua, ub, uc, ud], R[la, lb],
K[la, lb], Km[ua, lb], Ku[ua, ub],
CDK[la, lb, lc], CDKu[ua, ub, uc]},
Equations -> {
e4phi -> IfThen[conformalMethod == 1, 1/phi^2, Exp[4 phi]],
Dphi[la] -> IfThen[conformalMethod == 1,
- PD[phi, la] / (2 phi ),
PD[phi, la]],
D2phi[la, lb] -> IfThen[conformalMethod == 1,
- 1/ (2 phi) (PD[phi, la, lb] -
PD[phi, la] PD[phi, lb] / phi),
PD[phi, la, lb]],
Dg[la, lb, lc] -> e4phi * (4 Dphi[lc] gt[la, lb] + PD[gt[la, lb], lc]),
D2g[la, lb, lc, ld] -> e4phi * (4 gt[la, lb] (4 Dphi[lc] Dphi[ld] +
D2phi[lc, ld]) +
4 Dphi[lc] PD[gt[la, lb], ld] +
4 Dphi[ld] PD[gt[la, lb], lc] +
PD[gt[la, lb], lc, ld]),
g[la, lb] -> e4phi gt[la, lb],
gu[ua, ub] -> MatrixInverse[g[ua, ub]],
Gl[la, lb, lc] -> 1/2 (Dg[la, lb, lc] + Dg[la, lc, lb] - Dg[lb, lc, la]),
G[ua, lb, lc] -> gu[ua, ud] Gl[ld, lb, lc],
Riem[la, lb, lc, ld] -> 1/2 (D2g[la, ld, lb, lc] + D2g[lb, lc, la, ld] -
D2g[lb, ld, la, lc] - D2g[la, lc, lb, ld]) +
Gl[le, la, ld] G[ue, lb, lc] -
Gl[le, la, lc] G[ue, lb, ld],
(*
* kranc properly zeroes the zero-from-antisymmetry terms from Riemann above,
* but fails to do so when raising indices;
* therefore indicate the (anti)symmetries explicitly.
* results in much faster code
*)
Riemu[ua, ub, uc, ud] -> gu[ua, ue] gu[ub, uf] gu[uc, ug] gu[ud, uh]
1/4 (Riem[le, lf, lg, lh] - Riem[lf, le, lg, lh] -
Riem[le, lf, lh, lg] + Riem[lf, le, lh, lg]),
R[la, lb] -> gu[uc, ud] Riem[lc, la, ld, lb],
K[la, lb] -> e4phi At[la, lb] + 1/3 trK g[la, lb],
Km[ua, lb] -> gu[ua, uc] K[lc, lb],
Ku[ua, ub] -> gu[ub, uc] Km[ua, lc],
CDK[la, lb, lc] -> e4phi * (4 Dphi[lc] At[la, lb] + PD[At[la, lb], lc]) +
1/3 (PD[trK, lc] g[la, lb] + trK Dg[la, lb, lc]) -
G[ud, la, lc] K[ld, lb] - G[ud, lb, lc] K[ld, la],
CDKu[ua, ub, uc] -> gu[ua, ud] gu[ub, ue] gu[uc, uf] CDK[ld, le, lf],
term1 -> (Riem[la, lb, lc, ld] + K[la, lc] K[lb, ld] - K[la, ld] K[lb, lc])
(Riemu[ua, ub, uc, ud] + Ku[ua, uc] Ku[ub, ud] - Ku[ua, ud] Ku[ub, uc]),
term2 -> 8 (CDK[la, lb, lc] CDKu[ua, uc, ub] - CDK[la, lb, lc] CDKu[ua, ub, uc]),
term3 -> 4 gu[ua, uc] gu[ub, ud]
(R[la, lb] + trK K[la, lb] - K[la, le] Km[ue, lb])
(R[lc, ld] + trK K[lc, ld] - K[lc, le] Km[ue, ld]),
Kretsch -> term1 + term2 + term3
}
};
intParameters = {
(* FIXME: pull this parameter value from ML_BSSN ? *)
{
Name -> conformalMethod,
Description -> "Treatment of conformal factor",
AllowedValues -> {{Value -> "0", Description -> "phi method"},
{Value -> "1", Description -> "W method"}},
Default -> 0
},
{
Name -> fdOrder,
Default -> 4,
AllowedValues -> {2, 4, 6, 8}
}
};
CreateKrancThornTT[ groups, ".", thorn,
Calculations -> { kretschmannCalc },
PartialDerivatives -> derivatives,
EvolutionTimelevels -> 3,
DeclaredGroups -> declaredGroupNames,
DefaultEvolutionTimelevels -> 3,
UseJacobian -> True,
UseLoopControl -> True,
UseVectors -> True,
IntParameters -> intParameters,
InheritedImplementations -> { "ML_BSSN" }
];
|
