Peano
Loading...
Searching...
No Matches
GRMHD.py
Go to the documentation of this file.
1# This file is part of the ExaHyPE2 project. For conditions of distribution and
2# use, please see the copyright notice at www.peano-framework.org
3
4max_eigenvalue = r"""
5 // The eigenvalue represents the maximal velocity as a fraction of
6 // the speed of light, as such it can never be greater than 1, and
7 // returning 1 is safe as it an upper bound for the real eigenvalue
8 // return 1.0;
9
10 // PARAMETERS
11 constexpr double gamma = 4.0/3.0;
12
13 // Eigenvalues as given by Olindo
14 double V[NumberOfUnknowns];
15 //TODO: get whether this produced an error? make that an assertion?
16 GRMHD::PDE::PDECons2Prim(Q, V);
17
18 double rho = V[Shortcuts::rho];
19 double v_cov[3] = {
20 V[Shortcuts::vel + 0], V[Shortcuts::vel + 1], V[Shortcuts::vel + 2]
21 };
22 double p = V[Shortcuts::E];
23 double B_contr[3] = {
24 V[Shortcuts::B + 0], V[Shortcuts::B + 1], V[Shortcuts::B + 2]
25 };
26
27 double shift[3] = {
28 V[Shortcuts::shift + 0], V[Shortcuts::shift + 1], V[Shortcuts::shift + 2]
29 };
30 // double psi = V[Shortcuts::psi]; // unused variable
31 double lapse = V[Shortcuts::lapse];
32
33 double g_cov[3][3] = {
34 {V[Shortcuts::gij + 0], V[Shortcuts::gij + 1], V[Shortcuts::gij + 2]},
35 {V[Shortcuts::gij + 1], V[Shortcuts::gij + 3], V[Shortcuts::gij + 4]},
36 {V[Shortcuts::gij + 2], V[Shortcuts::gij + 4], V[Shortcuts::gij + 5]}
37 };
38
39 double g_contr[3][3];
40 double gp = GRMHD::PDE::matrixInverse3x3(g_cov, g_contr);
41 gp = std::sqrt(gp);
42 // double gm = 1./gp;
43
44 double v_contr[3];
45 GRMHD::PDE::matrixVectMult(g_contr, v_cov, v_contr);
46
47 double B_cov[3];
48 GRMHD::PDE::matrixVectMult(g_cov, B_contr, B_cov);
49
50 // evaluate useful quantities
51 double v2 = v_contr[0]*v_cov[0] + v_contr[1]*v_cov[1] + v_contr[2]*v_cov[2];
52 // double lf = 1.0/sqrt(1.0 - v2);
53 double lf2m = 1.0 - v2;
54 double b2 = B_contr[0]*B_cov[0] + B_contr[1]*B_cov[1] + B_contr[2]*B_cov[2];
55
56 double VdotB = v_contr[0]*B_cov[0] + v_contr[1]*B_cov[1] + v_contr[2]*B_cov[2];
57 double b2_4 = b2*lf2m + VdotB;
58 double gamma1 = gamma/(gamma-1.0);
59 double w = rho + gamma1*p + b2_4;
60 double cs2 = (gamma*p + b2_4)/w;
61
62 double u = v_contr[normal];
63 double sft = shift[normal];
64 double gg = g_contr[normal][normal];
65 double den = 1.0/(1.0 - v2*cs2);
66
67 double deltaU = sqrt( cs2*lf2m*( (1.0-v2*cs2)*gg - u*u*(1.0-cs2) ));
68
69 double maxEigenvalue = std::abs(lapse*u -sft);
70
71 maxEigenvalue = std::max( maxEigenvalue,
72 std::abs(
73 lapse*den*( u*(1.0-cs2) - deltaU ) - sft
74 ));
75
76 maxEigenvalue = std::max( maxEigenvalue,
77 std::abs(
78 lapse*den*( u*(1.0-cs2) + deltaU ) - sft
79 ));
80
81 maxEigenvalue = std::max( maxEigenvalue,
82 1.0
83 );
84
85 return maxEigenvalue;
86
87 //Safe mode: return also 'covariant' eigenvalues, should never be required
88 // double shift_cov[3];
89 // GRMHD::PDE::matrixVectMult(g_cov, shift, shift_cov);
90 // double vn2 = v_cov[normal];
91 // double sft2 = shift_cov[normal];
92 // double gg2 = g_cov[normal][normal];
93 // den = 1.0/(1.0 - v2*cs2);
94
95 // u = vn2;
96 // deltaU = sqrt( cs2*lf2m*( (1.0-v2*cs2)*gg2 - u*u*(1.0-cs2) ));
97
98 // maxEigenvalue = std::max( maxEigenvalue,
99 // std::abs(
100 // lapse*den*( u*(1.0-cs2) - deltaU ) - sft2
101 // ));
102
103 // maxEigenvalue = std::max( maxEigenvalue,
104 // std::abs(
105 // lapse*den*( u*(1.0-cs2) - deltaU ) - sft2
106 // ));
107
108 // maxEigenvalue = std::max( maxEigenvalue,
109 // std::abs(
110 // lapse*den*( u*(1.0-cs2) - deltaU ) - sft2
111 // ));
112"""
113
114flux = r"""
115 std::fill_n(F, NumberOfUnknowns, 0.0);
116
117 constexpr double gamma = 4.0/3.0;
118 double gamma1 = gamma/(gamma-1.0);
119
120 double V[NumberOfUnknowns];
121 //TODO: get whether this produced an error? make that an assertion?
122 GRMHD::PDE::PDECons2Prim(Q, V);
123
124 // int iErr;
125 // pdecons2prim_(V,Q,&iErr);
126
127 double rho = V[Shortcuts::rho];
128 double v_cov[3] = {
129 V[Shortcuts::vel + 0], V[Shortcuts::vel + 1], V[Shortcuts::vel + 2]
130 };
131 double p = V[Shortcuts::E];
132 double B_contr[3] = {
133 V[Shortcuts::B + 0], V[Shortcuts::B + 1], V[Shortcuts::B + 2]
134 };
135 double QB_contr[3] = {
136 Q[Shortcuts::B + 0], Q[Shortcuts::B + 1], Q[Shortcuts::B + 2]
137 };
138 double shift[3] = {
139 V[Shortcuts::shift + 0], V[Shortcuts::shift + 1], V[Shortcuts::shift + 2]
140 };
141 // double psi = V[Shortcuts::psi]; // unused variable
142 double lapse = V[Shortcuts::lapse];
143
144 double g_cov[3][3] = {
145 {V[Shortcuts::gij + 0], V[Shortcuts::gij + 1], V[Shortcuts::gij + 2]},
146 {V[Shortcuts::gij + 1], V[Shortcuts::gij + 3], V[Shortcuts::gij + 4]},
147 {V[Shortcuts::gij + 2], V[Shortcuts::gij + 4], V[Shortcuts::gij + 5]}
148 };
149
150 double g_contr[3][3];
151 double gp = GRMHD::PDE::matrixInverse3x3(g_cov, g_contr);
152 gp = std::sqrt(gp);
153 double gm = 1./gp;
154
155 double v_contr[3];
156 GRMHD::PDE::matrixVectMult(g_contr, v_cov, v_contr);
157
158 double S_contr[3];
159 GRMHD::PDE::matrixVectMult(g_contr, &Q[1], S_contr);
160
161 double BQ[3];
162 GRMHD::PDE::matrixVectMult(g_cov, QB_contr, BQ);
163
164 double B_cov[3];
165 GRMHD::PDE::matrixVectMult(g_cov, B_contr, B_cov);
166
167 // COVARIANT =>
168 double vxB_cov[3] = {
169 gp*(v_contr[1]*B_contr[2] - v_contr[2]*B_contr[1]),
170 gp*(v_contr[2]*B_contr[0] - v_contr[0]*B_contr[2]),
171 gp*(v_contr[0]*B_contr[1] - v_contr[1]*B_contr[0])
172 };
173
174 // CONTRAVARIANT =>
175 double vxB_contr[3] = {
176 gm*(v_cov[1]*B_cov[2] - v_cov[2]*B_cov[1]),
177 gm*(v_cov[2]*B_cov[0] - v_cov[0]*B_cov[2]),
178 gm*(v_cov[0]*B_cov[1] - v_cov[1]*B_cov[0])
179 };
180
181 double v2 = v_contr[0]*v_cov[0] + v_contr[1]*v_cov[1] + v_contr[2]*v_cov[2];
182 double e2 = vxB_contr[0]*vxB_cov[0] + vxB_contr[1]*vxB_cov[1] + vxB_contr[2]*vxB_cov[2];
183 double b2 = B_contr[0]*B_cov[0] + B_contr[1]*B_cov[1] + B_contr[2]*B_cov[2];
184
185 double uem = 0.5*(b2 + e2);
186 double lf = 1.0/sqrt(1.0 - v2);
187
188 double w = rho + gamma1*p;
189 double ww = w*lf*lf;
190 double wwx = ww*v_contr[0];
191 double wwy = ww*v_contr[1];
192 double wwz = ww*v_contr[2];
193
194 double Vtr[3] = {
195 lapse*v_contr[0]-shift[0], lapse*v_contr[1]-shift[1], lapse*v_contr[2]-shift[2]
196 };
197
198 double Fij[3][3];
199 for(int i=0; i<3; i++){
200 for(int j=0; j<3; j++){
201 Fij[i][j] = Vtr[i]*QB_contr[j]-Vtr[j]*QB_contr[i];
202 }
203 }
204
205 switch (normal) {
206 case 0:
207 F[0] = v_contr[0]*Q[0];
208 F[1] = gp*(wwx*v_cov[0] - vxB_contr[0]*vxB_cov[0] - B_contr[0]*B_cov[0] + p + uem);
209 F[2] = gp*(wwx*v_cov[1] - vxB_contr[0]*vxB_cov[1] - B_contr[0]*B_cov[1]);
210 F[3] = gp*(wwx*v_cov[2] - vxB_contr[0]*vxB_cov[2] - B_contr[0]*B_cov[2]);
211 F[4] = S_contr[0]-F[0];
212 F[5] = Fij[0][0];
213 F[6] = Fij[0][1];
214 F[7] = Fij[0][2];
215 break;
216 case 1:
217 F[0] = v_contr[1]*Q[0];
218 F[1] = gp*(wwy*v_cov[0] - vxB_contr[1]*vxB_cov[0] - B_contr[1]*B_cov[0]);
219 F[2] = gp*(wwy*v_cov[1] - vxB_contr[1]*vxB_cov[1] - B_contr[1]*B_cov[1] + p + uem);
220 F[3] = gp*(wwy*v_cov[2] - vxB_contr[1]*vxB_cov[2] - B_contr[1]*B_cov[2]);
221 F[4] = S_contr[1]-F[0];
222 F[5] = Fij[1][0];
223 F[6] = Fij[1][1];
224 F[7] = Fij[1][2];
225 break;
226 case 2:
227 F[0] = v_contr[2]*Q[0];
228 F[1] = gp*(wwz*v_cov[0] - vxB_contr[2]*vxB_cov[0] - B_contr[2]*B_cov[0]);
229 F[2] = gp*(wwz*v_cov[1] - vxB_contr[2]*vxB_cov[1] - B_contr[2]*B_cov[1]);
230 F[3] = gp*(wwz*v_cov[2] - vxB_contr[2]*vxB_cov[2] - B_contr[2]*B_cov[2] + p + uem);
231 F[4] = S_contr[2]-F[0];
232 F[5] = Fij[2][0];
233 F[6] = Fij[2][1];
234 F[7] = Fij[2][2];
235 break;
236 }
237
238 // Remember that Q(:) below contains already the factor gp, which is ok!
239 for(int i=0; i<5; i++){
240 F[i] = lapse*F[i] - shift[normal]*Q[i];
241 }
242"""
243
244ncp = r"""
245 std::fill_n(BTimesDeltaQ, NumberOfUnknowns, 0.0);
246
247 // PARAMETERS
248 constexpr double gamma = 4.0/3.0;
249 constexpr double DivCleaning_a = 1.0;
250
251 double lapse = Q[9];
252 double shift[3] = {
253 Q[Shortcuts::shift+0], Q[Shortcuts::shift+1], Q[Shortcuts::shift+2]
254 };
255
256 double g_cov[3][3] = {
257 {Q[13], Q[14], Q[15]},
258 {Q[14], Q[16], Q[17]},
259 {Q[15], Q[17], Q[18]}
260 };
261
262 double g_contr[3][3];
263 double gp = GRMHD::PDE::matrixInverse3x3(g_cov, g_contr);
264
265 gp = std::sqrt(gp);
266 double gm = 1./gp;
267
268 double Vc[NumberOfUnknowns];
269 //TODO: get whether this produced an error? make that an assertion?
270 GRMHD::PDE::PDECons2Prim(Q, Vc);
271
272 double gamma1 = gamma/(gamma-1.);
273
274 double rho = Vc[0];
275 double v_cov[3] = {
276 Vc[1], Vc[2], Vc[3]
277 };
278 double p = Vc[4];
279
280 double B_contr[3] = {
281 Vc[Shortcuts::B + 0], Vc[Shortcuts::B + 1], Vc[Shortcuts::B + 2]
282 };
283 // double QB_contr[3] = {
284 // Q[Shortcuts::B + 0], Q[Shortcuts::B + 1], Q[Shortcuts::B + 2]
285 // }; //unused variable
286
287 double B_cov[3];
288 GRMHD::PDE::matrixVectMult(g_cov, B_contr, B_cov);
289
290 // double psi = Vc[Shortcuts::psi];
291
292 double v_contr[3];
293 GRMHD::PDE::matrixVectMult(g_contr, v_cov, v_contr);
294
295 double S_contr[3];
296 GRMHD::PDE::matrixVectMult(g_contr, &Q[Shortcuts::vel], S_contr);
297
298 // COVARIANT =>
299 double vxB_cov[3] = {
300 gp*(v_contr[1]*B_contr[2] - v_contr[2]*B_contr[1]),
301 gp*(v_contr[2]*B_contr[0] - v_contr[0]*B_contr[2]),
302 gp*(v_contr[0]*B_contr[1] - v_contr[1]*B_contr[0])
303 };
304
305 // CONTRAVARIANT =>
306 double vxB_contr[3] = {
307 gm*(v_cov[1]*B_cov[2] - v_cov[2]*B_cov[1]),
308 gm*(v_cov[2]*B_cov[0] - v_cov[0]*B_cov[2]),
309 gm*(v_cov[0]*B_cov[1] - v_cov[1]*B_cov[0])
310 };
311
312 double v2 = v_contr[0]*v_cov[0] + v_contr[1]*v_cov[1] + v_contr[2]*v_cov[2];
313 double e2 = vxB_contr[0]*vxB_cov[0] + vxB_contr[1]*vxB_cov[1] + vxB_contr[2]*vxB_cov[2];
314 double b2 = B_contr[0]*B_cov[0] + B_contr[1]*B_cov[1] + B_contr[2]*B_cov[2];
315
316 double uem = 0.5*(b2 + e2);
317 double lf = 1.0/sqrt(1.0 - v2);
318
319 double w = rho + gamma1*p;
320 double ww = w*lf*lf;
321
322 double delta[3][3] = {
323 {1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}
324 };
325
326 //Consider unrolling these loops
327 int ccount = 0;
328 for(int i=0; i<3; i++){
329 double W_ij = ww*v_cov[i]*v_contr[normal]-vxB_cov[i]*vxB_contr[normal]-B_cov[i]*B_contr[normal]+(p+uem)*delta[i][normal];
330 BTimesDeltaQ[Shortcuts::vel + normal] -= Q[Shortcuts::vel+i]*deltaQ[Shortcuts::shift+i];
331 BTimesDeltaQ[Shortcuts::E] -= gp*W_ij*deltaQ[Shortcuts::shift+i];
332
333 for(int j=0; j<3; j++){
334 if(j>=i){
335 double Wim = ww*v_contr[i]*v_contr[j]-vxB_contr[i]*vxB_contr[j]-B_contr[i]*B_contr[j]+(p+uem)*g_contr[i][j];
336 if(i!=j){
337 Wim = 2.0*Wim;
338 }
339 BTimesDeltaQ[Shortcuts::vel+normal] -= 0.5*gp*lapse*Wim*deltaQ[Shortcuts::gij+ccount];
340 BTimesDeltaQ[Shortcuts::E] -= 0.5*gp*Wim*shift[normal]*deltaQ[Shortcuts::gij+ccount];
341 ccount += 1;
342 }
343 }
344 }
345
346 BTimesDeltaQ[Shortcuts::vel+normal] += (Q[Shortcuts::E]+Q[Shortcuts::rho])*deltaQ[Shortcuts::lapse];
347 BTimesDeltaQ[Shortcuts::E] += S_contr[normal]*deltaQ[Shortcuts::lapse];
348 BTimesDeltaQ[Shortcuts::B+0] = lapse*gp*g_contr[0][normal]*deltaQ[Shortcuts::psi]
349 - Q[Shortcuts::shift+0]*deltaQ[Shortcuts::B+normal];
350 BTimesDeltaQ[Shortcuts::B+1] = lapse*gp*g_contr[1][normal]*deltaQ[Shortcuts::psi]
351 - Q[Shortcuts::shift+1]*deltaQ[Shortcuts::B+normal];
352 BTimesDeltaQ[Shortcuts::B+2] = lapse*gp*g_contr[2][normal]*deltaQ[Shortcuts::psi]
353 - Q[Shortcuts::shift+2]*deltaQ[Shortcuts::B+normal];
354 BTimesDeltaQ[Shortcuts::psi] = gm*lapse*DivCleaning_a*DivCleaning_a*deltaQ[Shortcuts::B+normal]
355 - Q[Shortcuts::shift+normal]*deltaQ[Shortcuts::psi];
356"""