-
Notifications
You must be signed in to change notification settings - Fork 38
/
GiRaFFE_NRPy_C2P_P2C.py
133 lines (109 loc) · 6.4 KB
/
GiRaFFE_NRPy_C2P_P2C.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
from outputC import nrpyAbs # NRPy+: Core C code output module
import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends
import NRPy_param_funcs as par # NRPy+: Parameter interface
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 GRHD.equations as GRHD # NRPy+: Generate general relativistic hydrodynamics equations
import GRFFE.equations as GRFFE # NRPy+: Generate general relativisitic force-free electrodynamics equations
thismodule = __name__
# There are several C parameters that we will need in this module:
M_PI = par.Cparameters("#define",thismodule,["M_PI"], "")
GAMMA_SPEED_LIMIT = par.Cparameters("REAL",thismodule,"GAMMA_SPEED_LIMIT",10.0) # Default value based on
# IllinoisGRMHD.
# GiRaFFE default = 2000.0
# There are three fixes in this module; we might not always want to do all (or even any!) of them.
# So, we initialize some NRPy+ parameters to control this behavior
par.initialize_param(par.glb_param(type="bool", module=thismodule, parname="enforce_orthogonality_StildeD_BtildeU", defaultval=True))
par.initialize_param(par.glb_param(type="bool", module=thismodule, parname="enforce_speed_limit_StildeD", defaultval=True))
par.initialize_param(par.glb_param(type="bool", module=thismodule, parname="enforce_current_sheet_prescription", defaultval=True))
def GiRaFFE_NRPy_C2P(StildeD,BU,gammaDD,betaU,alpha):
GRHD.compute_sqrtgammaDET(gammaDD)
gammaUU,unusedgammadet = ixp.symm_matrix_inverter3x3(gammaDD)
BtildeU = ixp.zerorank1()
for i in range(3):
# \tilde{B}^i = B^i \sqrt{\gamma}
BtildeU[i] = GRHD.sqrtgammaDET*BU[i]
BtildeD = ixp.zerorank1()
for i in range(3):
for j in range(3):
BtildeD[j] += gammaDD[i][j]*BtildeU[i]
Btilde2 = sp.sympify(0)
for i in range(3):
Btilde2 += BtildeU[i]*BtildeD[i]
global outStildeD
outStildeD = StildeD
# Then, enforce the orthogonality:
if par.parval_from_str("enforce_orthogonality_StildeD_BtildeU"):
StimesB = sp.sympify(0)
for i in range(3):
StimesB += StildeD[i]*BtildeU[i]
for i in range(3):
# {\tilde S}_i = {\tilde S}_i - ({\tilde S}_j {\tilde B}^j) {\tilde B}_i/{\tilde B}^2
outStildeD[i] -= StimesB*BtildeD[i]/Btilde2
# Calculate \tilde{S}^2:
Stilde2 = sp.sympify(0)
for i in range(3):
for j in range(3):
Stilde2 += gammaUU[i][j]*outStildeD[i]*outStildeD[j]
# First we need to compute the factor f:
# f = \sqrt{(1-\Gamma_{\max}^{-2}){\tilde B}^4/(16 \pi^2 \gamma {\tilde S}^2)}
speed_limit_factor = sp.sqrt((sp.sympify(1)-GAMMA_SPEED_LIMIT**(-2.0))*Btilde2*Btilde2*sp.Rational(1,16)/\
(M_PI*M_PI*GRHD.sqrtgammaDET*GRHD.sqrtgammaDET*Stilde2))
import Min_Max_and_Piecewise_Expressions as noif
# Calculate B^2
B2 = sp.sympify(0)
for i in range(3):
for j in range(3):
B2 += gammaDD[i][j]*BU[i]*BU[j]
# Enforce the speed limit on StildeD:
if par.parval_from_str("enforce_speed_limit_StildeD"):
for i in range(3):
outStildeD[i] *= noif.min_noif(sp.sympify(1),speed_limit_factor)
global ValenciavU
ValenciavU = ixp.zerorank1()
# Recompute 3-velocity:
for i in range(3):
for j in range(3):
# \bar{v}^i = 4 \pi \gamma^{ij} {\tilde S}_j / (\sqrt{\gamma} B^2)
ValenciavU[i] += sp.sympify(4)*M_PI*gammaUU[i][j]*outStildeD[j]/(GRHD.sqrtgammaDET*B2)
# This number determines how far away (in grid points) we will apply the fix.
grid_points_from_z_plane = par.Cparameters("REAL",thismodule,"grid_points_from_z_plane",4.0)
if par.parval_from_str("enforce_current_sheet_prescription"):
# Calculate the drift velocity
driftvU = ixp.zerorank1()
for i in range(3):
driftvU[i] = alpha*ValenciavU[i] - betaU[i]
# The direct approach, used by the original GiRaFFE:
# v^z = -(\gamma_{xz} v^x + \gamma_{yz} v^y) / \gamma_{zz}
newdriftvU2 = -(gammaDD[0][2]*driftvU[0] + gammaDD[1][2]*driftvU[1])/gammaDD[2][2]
# Now that we have the z component, it's time to substitute its Valencia form in.
# Remember, we only do this if abs(z) < (k+0.01)*dz. Note that we add 0.01; this helps
# avoid floating point errors and division by zero. This is the same as abs(z) - (k+0.01)*dz<0
coord = nrpyAbs(rfm.xx[2])
bound =(grid_points_from_z_plane+sp.Rational(1,100))*gri.dxx[2]
ValenciavU[2] = noif.coord_leq_bound(coord,bound)*(newdriftvU2+betaU[2])/alpha \
+ noif.coord_greater_bound(coord,bound)*ValenciavU[2]
def GiRaFFE_NRPy_P2C(gammaDD,betaU,alpha, ValenciavU,BU, sqrt4pi):
# After recalculating the 3-velocity, we need to update the poynting flux:
# We'll reset the Valencia velocity, since this will be part of a second call to outCfunction.
# First compute stress-energy tensor T4UU and T4UD:
GRHD.compute_sqrtgammaDET(gammaDD)
# GRHD.u4U_in_terms_of_ValenciavU__rescale_ValenciavU_by_applying_speed_limit(alpha, betaU, gammaDD, ValenciavU)
R = sp.sympify(0)
for i in range(3):
for j in range(3):
R += gammaDD[i][j] * ValenciavU[i] * ValenciavU[j]
u4U_ito_ValenciavU = ixp.zerorank1(DIM=4)
u4U_ito_ValenciavU[0] = 1 / (alpha * sp.sqrt(1 - R))
# u^i = u^0 ( alpha v^i_{(n)} - beta^i ), where v^i_{(n)} is the Valencia 3-velocity
for i in range(3):
u4U_ito_ValenciavU[i + 1] = u4U_ito_ValenciavU[0] * (alpha * ValenciavU[i] - betaU[i])
GRFFE.compute_smallb4U_with_driftvU_for_FFE(gammaDD, betaU, alpha, u4U_ito_ValenciavU, BU, sqrt4pi)
GRFFE.compute_smallbsquared(gammaDD, betaU, alpha, GRFFE.smallb4_with_driftv_for_FFE_U)
GRFFE.compute_TEM4UU(gammaDD, betaU, alpha, GRFFE.smallb4_with_driftv_for_FFE_U, GRFFE.smallbsquared, u4U_ito_ValenciavU)
GRFFE.compute_TEM4UD(gammaDD, betaU, alpha, GRFFE.TEM4UU)
# Compute conservative variables in terms of primitive variables
GRHD.compute_S_tildeD(alpha, GRHD.sqrtgammaDET, GRFFE.TEM4UD)
global StildeD
StildeD = GRHD.S_tildeD