-
Notifications
You must be signed in to change notification settings - Fork 38
/
ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear.py
266 lines (238 loc) · 14.2 KB
/
ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear.py
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
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
# This module performs the conversion between ADM
# spacetime variables in Spherical or Cartesian coordinates
# given as *numerical* expressions (i.e., ADM quantities
# are only given to floating-point precision; e.g., in the
# case of an initial data solver), to rescaled BSSN-in-curvilinear
# coordinate quantities, as defined in BSSN_RHSs.py
# Author: Zachariah B. Etienne
# zachetie **at** gmail **dot* com
# Step P1: Initialize core Python/NRPy+ modules
from outputC import lhrh,outCfunction,outputC # NRPy+: Core C code output module
import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends
import finite_difference as fin # NRPy+: Finite difference C code generation module
import grid as gri # NRPy+: Functions having to do with numerical grids
import indexedexp as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm # NRPy+: Reference metric support
import BSSN.BSSN_quantities as Bq # NRPy+: Computes useful BSSN quantities; e.g., gammabarUU & GammabarUDD needed below
import os, sys # Standard Python modules for multiplatform OS-level functions
def Convert_Spherical_or_Cartesian_ADM_to_BSSN_curvilinear(CoordType_in, ADM_input_function_name,
Ccodesdir = "BSSN", pointer_to_ID_inputs=False,loopopts=",oldloops"):
# The ADM & BSSN formalisms only work in 3D; they are 3+1 decompositions of Einstein's equations.
# To implement axisymmetry or spherical symmetry, simply set all spatial derivatives in
# the relevant angular directions to zero; DO NOT SET DIM TO ANYTHING BUT 3.
# Step 0: Set spatial dimension (must be 3 for BSSN)
DIM = 3
# Step 1: All ADM initial data quantities are now functions of xx0,xx1,xx2, but
# they are still in the Spherical or Cartesian basis. We can now directly apply
# Jacobian transformations to get them in the correct xx0,xx1,xx2 basis:
# All input quantities are in terms of r,th,ph or x,y,z. We want them in terms
# of xx0,xx1,xx2, so here we call sympify_integers__replace_rthph() to replace
# r,th,ph or x,y,z, respectively, with the appropriate functions of xx0,xx1,xx2
# as defined for this particular reference metric in reference_metric.py's
# xxSph[] or xx_to_Cart[], respectively:
# Define the input variables:
gammaSphorCartDD = ixp.declarerank2("gammaSphorCartDD", "sym01")
KSphorCartDD = ixp.declarerank2("KSphorCartDD", "sym01")
alphaSphorCart = sp.symbols("alphaSphorCart")
betaSphorCartU = ixp.declarerank1("betaSphorCartU")
BSphorCartU = ixp.declarerank1("BSphorCartU")
# Make sure that rfm.reference_metric() has been called.
# We'll need the variables it defines throughout this module.
if rfm.have_already_called_reference_metric_function == False:
print("Error. Called Convert_Spherical_ADM_to_BSSN_curvilinear() without")
print(" first setting up reference metric, by calling rfm.reference_metric().")
sys.exit(1)
r_th_ph_or_Cart_xyz_oID_xx = []
if CoordType_in == "Spherical":
r_th_ph_or_Cart_xyz_oID_xx = rfm.xxSph
elif CoordType_in == "Cartesian":
r_th_ph_or_Cart_xyz_oID_xx = rfm.xx_to_Cart
else:
print("Error: Can only convert ADM Cartesian or Spherical initial data to BSSN Curvilinear coords.")
sys.exit(1)
# Step 2: All ADM initial data quantities are now functions of xx0,xx1,xx2, but
# they are still in the Spherical or Cartesian basis. We can now directly apply
# Jacobian transformations to get them in the correct xx0,xx1,xx2 basis:
# alpha is a scalar, so no Jacobian transformation is necessary.
alpha = alphaSphorCart
Jac_dUSphorCart_dDrfmUD = ixp.zerorank2()
for i in range(DIM):
for j in range(DIM):
Jac_dUSphorCart_dDrfmUD[i][j] = sp.diff(r_th_ph_or_Cart_xyz_oID_xx[i], rfm.xx[j])
Jac_dUrfm_dDSphorCartUD, dummyDET = ixp.generic_matrix_inverter3x3(Jac_dUSphorCart_dDrfmUD)
betaU = ixp.zerorank1()
BU = ixp.zerorank1()
gammaDD = ixp.zerorank2()
KDD = ixp.zerorank2()
for i in range(DIM):
for j in range(DIM):
betaU[i] += Jac_dUrfm_dDSphorCartUD[i][j] * betaSphorCartU[j]
BU[i] += Jac_dUrfm_dDSphorCartUD[i][j] * BSphorCartU[j]
for k in range(DIM):
for l in range(DIM):
gammaDD[i][j] += Jac_dUSphorCart_dDrfmUD[k][i] * Jac_dUSphorCart_dDrfmUD[l][j] * \
gammaSphorCartDD[k][l]
KDD[i][j] += Jac_dUSphorCart_dDrfmUD[k][i] * Jac_dUSphorCart_dDrfmUD[l][j] * KSphorCartDD[k][l]
# Step 3: All ADM quantities were input into this function in the Spherical or Cartesian
# basis, as functions of r,th,ph or x,y,z, respectively. In Steps 1 and 2 above,
# we converted them to the xx0,xx1,xx2 basis, and as functions of xx0,xx1,xx2.
# Here we convert ADM quantities in the "rfm" basis to their BSSN Curvilinear
# counterparts, for all BSSN quantities *except* lambda^i:
import BSSN.BSSN_in_terms_of_ADM as BitoA
BitoA.gammabarDD_hDD(gammaDD)
BitoA.trK_AbarDD_aDD(gammaDD, KDD)
BitoA.cf_from_gammaDD(gammaDD)
BitoA.betU_vetU(betaU, BU)
hDD = BitoA.hDD
trK = BitoA.trK
aDD = BitoA.aDD
cf = BitoA.cf
vetU = BitoA.vetU
betU = BitoA.betU
# Step 4: Compute $\bar{\Lambda}^i$ (Eqs. 4 and 5 of
# [Baumgarte *et al.*](https://arxiv.org/pdf/1211.6632.pdf)),
# from finite-difference derivatives of rescaled metric
# quantities $h_{ij}$:
# \bar{\Lambda}^i = \bar{\gamma}^{jk}\left(\bar{\Gamma}^i_{jk} - \hat{\Gamma}^i_{jk}\right).
# The reference_metric.py module provides us with analytic expressions for
# $\hat{\Gamma}^i_{jk}$, so here we need only compute
# finite-difference expressions for $\bar{\Gamma}^i_{jk}$, based on
# the values for $h_{ij}$ provided in the initial data. Once
# $\bar{\Lambda}^i$ has been computed, we apply the usual rescaling
# procedure:
# \lambda^i = \bar{\Lambda}^i/\text{ReU[i]},
# and then output the result to a C file using the NRPy+
# finite-difference C output routine.
# We will need all BSSN gridfunctions to be defined, as well as
# expressions for gammabarDD_dD in terms of exact derivatives of
# the rescaling matrix and finite-difference derivatives of
# hDD's. This functionality is provided by BSSN.BSSN_unrescaled_and_barred_vars,
# which we call here to overwrite above definitions of gammabarDD,gammabarUU, etc.
Bq.gammabar__inverse_and_derivs() # Provides gammabarUU and GammabarUDD
gammabarUU = Bq.gammabarUU
GammabarUDD = Bq.GammabarUDD
# Next evaluate \bar{\Lambda}^i, based on GammabarUDD above and GammahatUDD
# (from the reference metric):
LambdabarU = ixp.zerorank1()
for i in range(DIM):
for j in range(DIM):
for k in range(DIM):
LambdabarU[i] += gammabarUU[j][k] * (GammabarUDD[i][j][k] - rfm.GammahatUDD[i][j][k])
# Finally apply rescaling:
# lambda^i = Lambdabar^i/\text{ReU[i]}
lambdaU = ixp.zerorank1()
for i in range(DIM):
lambdaU[i] = LambdabarU[i] / rfm.ReU[i]
if ADM_input_function_name == "DoNotOutputADMInputFunction":
return hDD,aDD,trK,vetU,betU,alpha,cf,lambdaU
# Step 5.A: Output files containing finite-differenced lambdas.
outCparams = "preindent=1,outCfileaccess=a,outCverbose=False,includebraces=False"
lambdaU_expressions = [lhrh(lhs=gri.gfaccess("in_gfs", "lambdaU0"), rhs=lambdaU[0]),
lhrh(lhs=gri.gfaccess("in_gfs", "lambdaU1"), rhs=lambdaU[1]),
lhrh(lhs=gri.gfaccess("in_gfs", "lambdaU2"), rhs=lambdaU[2])]
desc = "Output lambdaU[i] for BSSN, built using finite-difference derivatives."
name = "ID_BSSN_lambdas"
params = "const paramstruct *restrict params,REAL *restrict xx[3],REAL *restrict in_gfs"
preloop = ""
enableCparameters=True
if "oldloops" in loopopts:
params = "const int Nxx[3],const int Nxx_plus_2NGHOSTS[3],REAL *xx[3],const REAL dxx[3],REAL *in_gfs"
enableCparameters=False
preloop = """
const REAL invdx0 = 1.0/dxx[0];
const REAL invdx1 = 1.0/dxx[1];
const REAL invdx2 = 1.0/dxx[2];
"""
outCfunction(
outfile=os.path.join(Ccodesdir, name + ".h"), desc=desc, name=name, params=params,
preloop=preloop,
body=fin.FD_outputC("returnstring", lambdaU_expressions, outCparams),
loopopts="InteriorPoints,Read_xxs"+loopopts, enableCparameters=enableCparameters)
# Step 5: Output all ADM-to-BSSN expressions to a C function. This function
# must first call the ID_ADM_SphorCart() defined above. Using these
# Spherical or Cartesian data, it sets up all quantities needed for
# BSSNCurvilinear initial data, *except* $\lambda^i$, which must be
# computed from numerical data using finite-difference derivatives.
ID_inputs_param = "ID_inputs other_inputs,"
if pointer_to_ID_inputs == True:
ID_inputs_param = "ID_inputs *other_inputs,"
desc = "Write BSSN variables in terms of ADM variables at a given point xx0,xx1,xx2"
name = "ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs"
enableCparameters=True
params = "const paramstruct *restrict params, "
if "oldloops" in loopopts:
enableCparameters=False
params = ""
params += "const REAL xx0xx1xx2[3]," + ID_inputs_param + """
REAL *hDD00,REAL *hDD01,REAL *hDD02,REAL *hDD11,REAL *hDD12,REAL *hDD22,
REAL *aDD00,REAL *aDD01,REAL *aDD02,REAL *aDD11,REAL *aDD12,REAL *aDD22,
REAL *trK,
REAL *vetU0,REAL *vetU1,REAL *vetU2,
REAL *betU0,REAL *betU1,REAL *betU2,
REAL *alpha, REAL *cf"""
outCparams = "preindent=1,outCverbose=False,includebraces=False"
outCfunction(
outfile=os.path.join(Ccodesdir, name + ".h"), desc=desc, name=name, params=params,
body="""
REAL gammaSphorCartDD00,gammaSphorCartDD01,gammaSphorCartDD02,
gammaSphorCartDD11,gammaSphorCartDD12,gammaSphorCartDD22;
REAL KSphorCartDD00,KSphorCartDD01,KSphorCartDD02,
KSphorCartDD11,KSphorCartDD12,KSphorCartDD22;
REAL alphaSphorCart,betaSphorCartU0,betaSphorCartU1,betaSphorCartU2;
REAL BSphorCartU0,BSphorCartU1,BSphorCartU2;
const REAL xx0 = xx0xx1xx2[0];
const REAL xx1 = xx0xx1xx2[1];
const REAL xx2 = xx0xx1xx2[2];
REAL xyz_or_rthph[3];\n""" +
outputC(r_th_ph_or_Cart_xyz_oID_xx[0:3], ["xyz_or_rthph[0]", "xyz_or_rthph[1]", "xyz_or_rthph[2]"],
"returnstring",
outCparams + ",CSE_enable=False") + " " + ADM_input_function_name + """(xyz_or_rthph, other_inputs,
&gammaSphorCartDD00,&gammaSphorCartDD01,&gammaSphorCartDD02,
&gammaSphorCartDD11,&gammaSphorCartDD12,&gammaSphorCartDD22,
&KSphorCartDD00,&KSphorCartDD01,&KSphorCartDD02,
&KSphorCartDD11,&KSphorCartDD12,&KSphorCartDD22,
&alphaSphorCart,&betaSphorCartU0,&betaSphorCartU1,&betaSphorCartU2,
&BSphorCartU0,&BSphorCartU1,&BSphorCartU2);
// Next compute all rescaled BSSN curvilinear quantities:\n""" +
outputC([hDD[0][0], hDD[0][1], hDD[0][2], hDD[1][1], hDD[1][2], hDD[2][2],
aDD[0][0], aDD[0][1], aDD[0][2], aDD[1][1], aDD[1][2], aDD[2][2],
trK, vetU[0], vetU[1], vetU[2], betU[0], betU[1], betU[2],
alpha, cf],
["*hDD00", "*hDD01", "*hDD02", "*hDD11", "*hDD12", "*hDD22",
"*aDD00", "*aDD01", "*aDD02", "*aDD11", "*aDD12", "*aDD22",
"*trK", "*vetU0", "*vetU1", "*vetU2", "*betU0", "*betU1", "*betU2",
"*alpha", "*cf"], "returnstring", params=outCparams),
enableCparameters=enableCparameters)
# Step 5.a: Output the driver function for the above
# function ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs()
# Next write the driver function for ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs():
desc = """Driver function for ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs(),
which writes BSSN variables in terms of ADM variables at a given point xx0,xx1,xx2"""
name = "ID_BSSN__ALL_BUT_LAMBDAs"
params = "const paramstruct *restrict params,REAL *restrict xx[3]," + ID_inputs_param + "REAL *in_gfs"
enableCparameters = True
funccallparams = "params, "
idx3replace = "IDX3S"
idx4ptreplace = "IDX4ptS"
if "oldloops" in loopopts:
params = "const int Nxx_plus_2NGHOSTS[3],REAL *xx[3]," + ID_inputs_param + "REAL *in_gfs"
enableCparameters = False
funccallparams = ""
idx3replace = "IDX3"
idx4ptreplace = "IDX4pt"
outCfunction(
outfile=os.path.join(Ccodesdir, name + ".h"), desc=desc, name=name, params=params,
body="""
const int idx = IDX3(i0,i1,i2);
const REAL xx0xx1xx2[3] = {xx0,xx1,xx2};
ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs(""".replace("IDX3",idx3replace)+funccallparams+"""xx0xx1xx2,other_inputs,
&in_gfs[IDX4pt(HDD00GF,idx)],&in_gfs[IDX4pt(HDD01GF,idx)],&in_gfs[IDX4pt(HDD02GF,idx)],
&in_gfs[IDX4pt(HDD11GF,idx)],&in_gfs[IDX4pt(HDD12GF,idx)],&in_gfs[IDX4pt(HDD22GF,idx)],
&in_gfs[IDX4pt(ADD00GF,idx)],&in_gfs[IDX4pt(ADD01GF,idx)],&in_gfs[IDX4pt(ADD02GF,idx)],
&in_gfs[IDX4pt(ADD11GF,idx)],&in_gfs[IDX4pt(ADD12GF,idx)],&in_gfs[IDX4pt(ADD22GF,idx)],
&in_gfs[IDX4pt(TRKGF,idx)],
&in_gfs[IDX4pt(VETU0GF,idx)],&in_gfs[IDX4pt(VETU1GF,idx)],&in_gfs[IDX4pt(VETU2GF,idx)],
&in_gfs[IDX4pt(BETU0GF,idx)],&in_gfs[IDX4pt(BETU1GF,idx)],&in_gfs[IDX4pt(BETU2GF,idx)],
&in_gfs[IDX4pt(ALPHAGF,idx)],&in_gfs[IDX4pt(CFGF,idx)]);
""".replace("IDX4pt",idx4ptreplace),
loopopts="AllPoints,Read_xxs"+loopopts, enableCparameters=enableCparameters)