Peano
Loading...
Searching...
No Matches
PDE.h
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#pragma once
4
5namespace GRMHD {
6 namespace PDE {
7 using Shortcuts = tests::exahype2::aderdg::VariableShortcutsADERDGSolver;
8
9 // Computes the vector resulting from the multiplication of a 3x3 matrix with a
10 // vector of size 3
11 inline void matrixVectMult(const double M[][3], const double VIn[3], double VOut[3]) {
12 VOut[0] = M[0][0] * VIn[0] + M[0][1] * VIn[1] + M[0][2] * VIn[2];
13 VOut[1] = M[1][0] * VIn[0] + M[1][1] * VIn[1] + M[1][2] * VIn[2];
14 VOut[2] = M[2][0] * VIn[0] + M[2][1] * VIn[1] + M[2][2] * VIn[2];
15 }
16
17 // Computes the inverse of a 3x3 matrix and returns its determinant
18 inline double matrixInverse3x3(double M[][3], double iM[][3]) {
19 double det = M[0][0] * M[1][1] * M[2][2] - M[0][0] * M[1][2] * M[2][1] - M[1][0] * M[0][1] * M[2][2]
20 + M[1][0] * M[0][2] * M[2][1] + M[2][0] * M[0][1] * M[1][2]
21 - M[2][0] * M[0][2] * M[1][1]; // ComputeDet(M,3)
22
23 iM[0][0] = (M[1][1] * M[2][2] - M[1][2] * M[2][1]) / det;
24 iM[0][1] = (M[0][2] * M[2][1] - M[0][1] * M[2][2]) / det;
25 iM[0][2] = (M[0][1] * M[1][2] - M[0][2] * M[1][1]) / det;
26 iM[1][0] = (M[1][2] * M[2][0] - M[1][0] * M[2][2]) / det;
27 iM[1][1] = (M[0][0] * M[2][2] - M[0][2] * M[2][0]) / det;
28 iM[1][2] = (M[0][2] * M[1][0] - M[0][0] * M[1][2]) / det;
29 iM[2][0] = (M[1][0] * M[2][1] - M[1][1] * M[2][0]) / det;
30 iM[2][1] = (M[0][1] * M[2][0] - M[0][0] * M[2][1]) / det;
31 iM[2][2] = (M[0][0] * M[1][1] - M[0][1] * M[1][0]) / det;
32
33 return det;
34 }
35
36 inline void FUNC_C2P_RMHD1(
37 double x,
38 double* f,
39 double* df,
40 double gam,
41 double d,
42 double e,
43 double s2,
44 double b2,
45 double sb2,
46 double* w
47 ) {
48 // INTENT(IN) :: x,gam,d,e,s2,b2,sb2
49 // INTENT(OUT):: f,df,w
50
51 // ! This is the CONS2PRIM strategy adopted by Del Zanna et al. (2007) A&A, 473, 11-30
52 // ! and it corresponds to their choice 3 in Section 3.2
53 constexpr double tolerance = 1e-10;
54 constexpr double third = 1. / 3.;
55 constexpr double i27 = 1. / 27.;
56
57 double v2 = x;
58 double rho = d * sqrt(1. - v2);
59 double c3 = 1. - gam * (1. - v2);
60 double c2 = gam * rho + .5 * b2 * (1 + v2) - e;
61 double c0 = -0.5 * sb2;
62
63 double c2ic3 = c2 / c3;
64 if (std::abs(c0) < 1.0e-20) {
65 w[0] = -c2ic3;
66 } else {
67 double c2ic33 = c2ic3 * c2ic3 * c2ic3;
68 double q = c0 / c3 + 2.0 * i27 * c2ic33;
69 double p = -third * c2ic3 * c2ic3;
70 double sqrtdet = std::sqrt(0.25 * q * q + i27 * p * p * p);
71 w[0] = -third * c2ic3 + std::pow(-0.5 * q + sqrtdet, third) + std::pow(-0.5 * q - sqrtdet, third);
72 // w[0] = std::max(-c2/c3, std::pow(-c0/c3, third));
73 // for(int i=0; i<100; i++){
74 // double dw = -((c3*w[0] + c2)*w[0]*w[0] + c0)/((3*c3*w[0] + 2*c2)*w[0]);
75 // if(std::abs(dw/w[0])<tolerance){
76 // break;
77 // }
78 // w[0] = w[0] + dw;
79 // }
80 }
81
82 double dc3 = gam;
83 double dc2 = 0.5 * (b2 - gam * rho / (1.0 - v2));
84
85 double dlogw = (dc3 * w[0] + dc2) / (3.0 * c3 * w[0] + 2.0 * c2);
86 double wb = w[0] + b2;
87 double vb2 = sb2 / (w[0] * w[0]);
88 f[0] = wb * wb * v2 - (2.0 * w[0] + b2) * vb2 - s2;
89 df[0] = wb * (wb + 2.0 * dlogw * (w[0] * v2 + vb2));
90 }
91
92 inline void RTSAFE_C2P_RMHD1(
93 double* v2,
94 double* X1,
95 double* X2,
96 double* XACC,
97 double* w,
98 double gam,
99 double d,
100 double e,
101 double s2,
102 double b2,
103 double sb2,
104 bool* failed
105 ) {
106 // REAL, INTENT(INOUT) :: v2,X1,X2,XACC,w
107 // REAL, INTENT(IN) :: gam,d,e,s2,b2,sb2
108 // LOGICAL, INTENT(INOUT) :: FAILED
109
110 double xacc2 = 1e-60;
111 failed[0] = false;
112
113 double FL, FH, DF;
114
115 FUNC_C2P_RMHD1(*X1, &FL, &DF, gam, d, e, s2, b2, sb2, w);
116
117 if (std::abs(FL) < xacc2) {
118 v2[0] = X1[0];
119 return;
120 }
121
122 FUNC_C2P_RMHD1(*X2, &FH, &DF, gam, d, e, s2, b2, sb2, w);
123
124 if (std::abs(FH) < xacc2) {
125 v2[0] = X2[0];
126 return;
127 }
128
129 if (FL * FH > 0.) {
130 failed[0] = true;
131 v2[0] = 0.; // assign value even if <it fails
132 return;
133 }
134
135 double XL, XH, SWAP;
136
137 if (FL < 0.) {
138 XL = X1[0];
139 XH = X2[0];
140 } else {
141 XH = X1[0];
142 XL = X2[0];
143 SWAP = FL;
144 FL = FH;
145 FH = SWAP;
146 }
147
148 v2[0] = .5 * (X1[0] + X2[0]);
149 double dxOld = std::abs(X2[0] - X1[0]);
150 double dx = dxOld;
151
152 double F;
153 FUNC_C2P_RMHD1(v2[0], &F, &DF, gam, d, e, s2, b2, sb2, w);
154
155 constexpr int MaxIt = 400;
156 for (int j = 0; j < MaxIt; j++) {
157 if ((((v2[0] - XH) * DF - F) * ((v2[0] - XL) * DF - F) > 0.) or (std::abs(2 * F) > std::abs(dxOld * DF))) {
158 dxOld = dx;
159 dx = 0.5 * (XH - XL);
160 v2[0] = XL + dx;
161 if (XL == v2[0]) {
162 return;
163 }
164 } else {
165 dxOld = dx;
166 dx = F / DF;
167 double TEMP = v2[0];
168 v2[0] = v2[0] - dx;
169 if (TEMP == v2[0]) {
170 return;
171 }
172 }
173
174 if (std::abs(dx) < XACC[0]) {
175 return;
176 }
177
178 FUNC_C2P_RMHD1(v2[0], &F, &DF, gam, d, e, s2, b2, sb2, w);
179
180 if (F < 0) {
181 XL = v2[0];
182 FL = F;
183 } else {
184 XH = v2[0];
185 FH = F;
186 }
187
188 if (std::abs(F) < xacc2) {
189 return;
190 }
191 } // for j
192
193 failed[0] = true;
194 v2[0] = 0.;
195 return;
196 }
197
198 // transforms the PDE from its formulation in primitive forms to its formulation
199 // in conservative form
200 inline void PDEPrim2Cons(double* Q, const double* V) {
201 constexpr double gamma = 4.0 / 3.0;
202
203 double rho = V[0];
204 double v_cov[3] = {V[1], V[2], V[3]};
205 double p = V[4];
206 double B_contr[3] = {V[5], V[6], V[7]};
207 // double psi = V[8];
208 // double lapse = V[9];
209 // double shift[3] = {V[10], V[11], V[12]};
210
211 // NOTE: not in order as there are six values being assigned to a matrix with
212 // nine elements
213 double g_cov[3][3] = {{V[13], V[14], V[15]}, {V[14], V[16], V[17]}, {V[15], V[17], V[18]}};
214
215 double g_contr[3][3];
216 double gp = matrixInverse3x3(g_cov, g_contr);
217 gp = std::sqrt(gp);
218
219 double gm = 1. / gp;
220
221 double v_contr[3];
222 matrixVectMult(g_cov, v_cov, v_contr);
223
224 double B_cov[3];
225 matrixVectMult(g_cov, B_contr, B_cov);
226
227 // COVARIANT ==>
228 double vxB_cov[3] = {
229 gp * (v_contr[1] * B_contr[2] - v_contr[2] * B_contr[1]),
230 gp * (v_contr[2] * B_contr[0] - v_contr[0] * B_contr[2]),
231 gp * (v_contr[0] * B_contr[1] - v_contr[1] * B_contr[0])
232 };
233
234 // CONTRAVARIANT ==>
235 double vxB_contr[3] = {
236 gm * (v_cov[1] * B_cov[2] - v_cov[2] * B_cov[1]),
237 gm * (v_cov[2] * B_cov[0] - v_cov[0] * B_cov[2]),
238 gm * (v_cov[0] * B_cov[1] - v_cov[1] * B_cov[0])
239 };
240
241 double vb = v_cov[0] * B_contr[0] + v_cov[1] * B_contr[1] + v_cov[2] * B_contr[2];
242 double v2 = v_contr[0] * v_cov[0] + v_contr[1] * v_cov[1] + v_contr[2] * v_cov[2];
243 double e2 = vxB_contr[0] * vxB_cov[0] + vxB_contr[1] * vxB_cov[1] + vxB_contr[2] * vxB_cov[2];
244 double b2 = B_contr[0] * B_cov[0] + B_contr[1] * B_cov[1] + B_contr[2] * B_cov[2];
245
246 double uem = 0.5 * (b2 + e2);
247
248 double lf = 1.0 / sqrt(1.0 - v2);
249 double gamma1 = gamma / (gamma - 1.0);
250 double w = rho + gamma1 * p;
251 double ww = w * lf * lf;
252
253 // using Shortcuts = tests::exahype2::aderdg::VariableShortcutsADERDGSolver;
254 Q[Shortcuts::rho] = gp * rho * lf;
255 Q[Shortcuts::vel + 0] = gp * ww * v_cov[0] + b2 * v_cov[0] - vb * B_cov[0];
256 Q[Shortcuts::vel + 1] = gp * ww * v_cov[1] + b2 * v_cov[1] - vb * B_cov[1];
257 Q[Shortcuts::vel + 2] = gp * ww * v_cov[2] + b2 * v_cov[2] - vb * B_cov[2];
258 Q[Shortcuts::E] = gp * ww - p + uem - Q[Shortcuts::rho]; // we subtract PDE(Q[rho]))
259 Q[Shortcuts::B + 0] = gp * V[Shortcuts::B + 0];
260 Q[Shortcuts::B + 1] = gp * V[Shortcuts::B + 1];
261 Q[Shortcuts::B + 2] = gp * V[Shortcuts::B + 2];
262
263 // simple copy
264 Q[Shortcuts::psi] = V[Shortcuts::psi];
265 Q[Shortcuts::lapse] = V[Shortcuts::lapse];
266 Q[Shortcuts::shift + 0] = V[Shortcuts::shift + 0];
267 Q[Shortcuts::shift + 1] = V[Shortcuts::shift + 1];
268 Q[Shortcuts::shift + 2] = V[Shortcuts::shift + 2];
269 Q[Shortcuts::gij + 0] = V[Shortcuts::gij + 0];
270 Q[Shortcuts::gij + 1] = V[Shortcuts::gij + 1];
271 Q[Shortcuts::gij + 2] = V[Shortcuts::gij + 2];
272 Q[Shortcuts::gij + 3] = V[Shortcuts::gij + 3];
273 Q[Shortcuts::gij + 4] = V[Shortcuts::gij + 4];
274 Q[Shortcuts::gij + 5] = V[Shortcuts::gij + 5];
275 }
276
277 // transforms the PDE from its formulation in conservative form to its
278 // formulation in primitive form
279 inline void PDECons2Prim(const double* Q, double* V) {
280 // PARAMETERS
281 constexpr double gamma = 4.0 / 3.0;
282 constexpr double p_floor = 1.0e-16;
283 constexpr double rho_floor = 1.0e-10;
284
285 constexpr double igamma = 1.0 / gamma;
286 constexpr int NSTOV_kappa = 100;
287 // constexpr double NSTOV_p_atmo = 1e-15;
288 // constexpr double NSTOV_rho_atmo = std::pow(NSTOV_p_atmo/NSTOV_kappa, igamma);
289 constexpr double NSTOV_p_atmo = 1e-12;
290 constexpr double NSTOV_rho_atmo = 1e-11;
291
292 // constexpr double tol = 1.0e-12;
293
294 // using Shortcuts = tests::exahype2::aderdg::VariableShortcutsADERDGSolver;
295 double psi = Q[Shortcuts::psi];
296 double lapse = Q[Shortcuts::lapse];
297 double shift[3] = {Q[Shortcuts::shift + 0], Q[Shortcuts::shift + 1], Q[Shortcuts::shift + 2]};
298 double gammaij[6] = {
299 Q[Shortcuts::gij + 0],
300 Q[Shortcuts::gij + 1],
301 Q[Shortcuts::gij + 2],
302 Q[Shortcuts::gij + 3],
303 Q[Shortcuts::gij + 4],
304 Q[Shortcuts::gij + 5]
305 };
306 double g_cov[3][3] = {
307 // TODO: switch these to gij notation
308 {Q[13], Q[14], Q[15]},
309 {Q[14], Q[16], Q[17]},
310 {Q[15], Q[17], Q[18]}
311 };
312
313 // TODO: turn this into a proper assertion
314 // IF(SUM(g_cov**2).LT.1e-9) THEN
315 // print *, 'FATAL ERROR: g_cov = 0'
316
317 double g_contr[3][3];
318 double gp = matrixInverse3x3(g_cov, g_contr);
319
320 gp = std::sqrt(gp);
321 double gm = 1. / gp;
322
323 double Qloc[19];
324 for (int i = 0; i < 8; i++) {
325 Qloc[i] = gm * Q[i];
326 }
327
328 double dd = Qloc[0];
329
330 double B_contr[3] = {
331 Qloc[5 + 0],
332 Qloc[5 + 1],
333 Qloc[5 + 2] // wrong: magnetic field is contravariant
334 };
335 double sm_cov[3] = {Qloc[1], Qloc[2], Qloc[3]};
336
337 double gamma1 = gamma / (gamma - 1.0);
338 double gam = 1.0 / gamma1;
339
340 // Solve for p
341 bool failed = false;
342
343 double sm[3];
344 matrixVectMult(g_contr, sm_cov, sm);
345 double B_cov[3];
346 matrixVectMult(g_cov, B_contr, B_cov);
347
348 double s2 = sm_cov[0] * sm[0] + sm_cov[1] * sm[1] + sm_cov[2] * sm[2];
349 double b2 = B_contr[0] * B_cov[0] + B_contr[1] * B_cov[1] + B_contr[2] * B_cov[2];
350 double sb = sm_cov[0] * B_contr[0] + sm_cov[1] * B_contr[1] + sm_cov[2] * B_contr[2];
351 double sb2 = sb * sb;
352
353 // First option [Del Zanna et al. (2007) A&A, 473, 11-30 (method 3)]
354
355 double e = Qloc[4] + dd; // Q(5) = gamma^1/2 ( U - dd )
356 // double x1 = 0.;
357 double eps = 1.0e-10;
358 // double x2 = 1.0-eps; // !tol !
359 // double w=0;
360
361 // double x1, x2, tol, v2;
362 // CALL RTSAFE_C2P_RMHD1(v2,x1,x2,tol,gam,dd,e,s2,b2,sb2,w,FAILED)
363 // double v2 = RTSAFE_C2P_RMHD1(x1, x2, tol, gam, dd, e, s2, b2, sb2, &w, &failed);
364
365 double x1 = 0.;
366 double x2 = 1.0 - 1e-10;
367 double tol = 1.0e-18;
368 double w = 0;
369 double v2;
370 RTSAFE_C2P_RMHD1(&v2, &x1, &x2, &tol, &w, gam, dd, e, s2, b2, sb2, &failed);
371
372 double p, rho, den, vb;
373 double v_cov[3];
374 if (failed) {
375 // iErr = -1
376 p = NSTOV_p_atmo;
377 rho = NSTOV_rho_atmo;
378
379 v_cov[0] = 0.0;
380 v_cov[1] = 0.0;
381 v_cov[2] = 0.0;
382 } else {
383 den = 1.0 / (w + b2);
384 vb = sb / w;
385 rho = dd * sqrt(1. - v2);
386 v_cov[0] = (sm_cov[0] + vb * B_cov[0]) * den;
387 v_cov[1] = (sm_cov[1] + vb * B_cov[1]) * den;
388 v_cov[2] = (sm_cov[2] + vb * B_cov[2]) * den;
389 p = gam * (w * (1. - v2) - rho);
390 }
391
392 if (rho < NSTOV_rho_atmo * 1.01) {
393 v_cov[0] = 0.0;
394 v_cov[1] = 0.0;
395 v_cov[2] = 0.0;
396 v2 = 0.;
397 rho = NSTOV_rho_atmo;
398 p = gam * (w * (1. - v2) - rho);
399 }
400 if (p < NSTOV_p_atmo * 1.01) {
401 v_cov[0] = 0.0;
402 v_cov[1] = 0.0;
403 v_cov[2] = 0.0;
404 v2 = 0.;
405 p = gam * (w * (1. - v2) - rho);
406 p = NSTOV_p_atmo;
407 }
408 if (p < NSTOV_p_atmo) {
409 p = NSTOV_p_atmo;
410 }
411
412 V[0] = rho;
413 V[1] = v_cov[0];
414 V[2] = v_cov[1];
415 V[3] = v_cov[2];
416 V[4] = p;
417 V[5] = B_contr[0];
418 V[6] = B_contr[1];
419 V[7] = B_contr[2];
420 V[8] = psi;
421 V[9] = lapse;
422 V[10] = shift[0];
423 V[11] = shift[1];
424 V[12] = shift[2];
425
426 for (int i = 0; i < 6; i++) {
427 V[13 + i] = gammaij[i];
428 }
429 } // PDECons2Prim
430 } // namespace PDE
431} // namespace GRMHD
GRMHD::PDE::FUNC_C2P_RMHD1
void FUNC_C2P_RMHD1(double x, double *f, double *df, double gam, double d, double e, double s2, double b2, double sb2, double *w)
Definition PDE.h:36
GRMHD::PDE::PDECons2Prim
void PDECons2Prim(const double *Q, double *V)
Definition PDE.h:279
GRMHD::PDE::matrixVectMult
void matrixVectMult(const double M[][3], const double VIn[3], double VOut[3])
Definition PDE.h:11
GRMHD::PDE::matrixInverse3x3
double matrixInverse3x3(double M[][3], double iM[][3])
Definition PDE.h:18
GRMHD::PDE::Shortcuts
tests::exahype2::aderdg::VariableShortcutsADERDGSolver Shortcuts
Definition PDE.h:7
GRMHD::PDE::PDEPrim2Cons
void PDEPrim2Cons(double *Q, const double *V)
Definition PDE.h:200
GRMHD::PDE::RTSAFE_C2P_RMHD1
void RTSAFE_C2P_RMHD1(double *v2, double *X1, double *X2, double *XACC, double *w, double gam, double d, double e, double s2, double b2, double sb2, bool *failed)
Definition PDE.h:92
GRMHD
Definition GRMHD.py:1
PDE
Definition PDE.py:1