Peano
Loading...
Searching...
No Matches
InitialData.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
5#include "PDE.h"
6
7namespace GRMHD {
8 namespace InitialData {
9 // Currently this is only a helper function for the following initial conditions
10 inline void metric(
11 const tarch::la::Vector<DIMENSIONS, double>& x,
12 double* lapse,
13 double* gp,
14 double* gm, // size 1
15 double* shift, // size 3
16 double** Kex,
17 double g_cov[][3],
18 double g_contr[][3], // all size 3x3
19 double* phi
20 ) {
21 constexpr double aom = 0.0;
22 constexpr double Mbh = 1.0;
23
24 phi[0] = 1.0;
25
26 // Rotating black hole in Kerr-Schild Cartesian coordinates. See De Felice & Clarke Sect. 11.4
27 double z2 = x[2] * x[2];
28 double aom2 = aom * aom;
29
30 double d_aom2 = 0.5 * (x[0] * x[0] + x[1] * x[1] + x[2] * x[2] - aom2);
31 // r = SQRT( (x**2 + y**2 + z**2 - aom2)/2.0 + SQRT(((x**2 + y**2 + z**2 - aom2)/2.0)**2 + z2*aom2))
32 double r = std::sqrt(d_aom2 + std::sqrt(d_aom2 * d_aom2 + z2 * aom2));
33 double r2 = r * r;
34
35 double HH = Mbh * r2 * r / (r2 * r2 + aom2 * z2);
36 double SS = 1.0 + 2.0 * HH;
37
38 double lx = (r * x[0] + aom * x[1]) / (r2 + aom2);
39 double ly = (r * x[1] - aom * x[0]) / (r2 + aom2);
40 double lz = x[2] / r;
41
42 lapse[0] = 1.0 / std::sqrt(SS);
43 shift[0] = 2.0 * HH / SS * lx;
44 shift[1] = 2.0 * HH / SS * ly;
45 shift[2] = 2.0 * HH / SS * lz;
46
47 g_cov[0][0] = 1.0 + 2.0 * HH * lx * lx;
48 g_cov[0][1] = 2.0 * HH * lx * ly;
49 g_cov[0][2] = 2.0 * HH * lx * lz;
50
51 g_cov[1][0] = 2.0 * HH * lx * ly;
52 g_cov[1][1] = 1.0 + 2.0 * HH * ly * ly;
53 g_cov[1][2] = 2.0 * HH * ly * lz;
54
55 g_cov[2][0] = 2.0 * HH * lx * lz;
56 g_cov[2][1] = 2.0 * HH * ly * lz;
57 g_cov[2][2] = 1.0 + 2.0 * HH * lz * lz;
58
59 g_contr[0][0] = (1.0 + 2.0 * HH * ly * ly + 2.0 * HH * lz * lz) / SS;
60 g_contr[0][1] = (-2.0 * HH * lx * ly) / SS;
61 g_contr[0][2] = (-2.0 * HH * lx * lz) / SS;
62
63 g_contr[1][0] = (-2.0 * HH * lx * ly) / SS;
64 g_contr[1][1] = (1.0 + 2.0 * HH * lx * lx + 2.0 * HH * lz * lz) / SS;
65 g_contr[1][2] = (-2.0 * HH * ly * lz) / SS;
66
67 g_contr[2][0] = (-2.0 * HH * lx * lz) / SS;
68 g_contr[2][1] = (-2.0 * HH * ly * lz) / SS;
69 g_contr[2][2] = (1.0 + 2.0 * HH * lx * lx + 2.0 * HH * ly * ly) / SS;
70
71 gp[0] = std::sqrt(SS);
72 gm[0] = 1.0 / gp[0];
73 }
74
75 // Initial conditions and analytical solution for the "Michel accretion test" described e.g.
76 // in https://doi.org/10.1016/j.cpc.2020.107251 or http://dx.doi.org/10.1007/BF00649949
77 // this test represents the spherical accretion onto a stationary black hole.
78 inline void initialAccretionDisc3D(
79 const tarch::la::Vector<DIMENSIONS, double>& x,
80 const double t,
81 double* const Q
82 ) {
83 // PARAMETERS
84 constexpr double gamma = 4.0 / 3.0;
85
86 double rhoc = 0.0625; // = 1/16, ICSphericalAccretionrhoc
87 double rc = 8.0; // ICSphericalAccretionrc
88
89 int MaxNewton = 50;
90
91 double ng = 1.0 / (gamma - 1.0);
92
93 double lapse, gp, gm, phi;
94 double shift[3];
95 double g_cov[3][3], g_contr[3][3];
96
97 metric(x, &lapse, &gp, &gm, shift, nullptr, g_cov, g_contr, &phi);
98
99 double r = sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
100
101 double V[19];
102 if (r < 0.8) {
103 // To avoid division by zero, never used for evolution or BC
104 using Shortcuts = tests::exahype2::aderdg::VariableShortcutsADERDGSolver;
105 V[Shortcuts::rho] = 1;
106 V[Shortcuts::vel + 0] = 0.;
107 V[Shortcuts::vel + 1] = 0.;
108 V[Shortcuts::vel + 2] = 0.;
109 V[Shortcuts::E] = 1.; // should be p?, likely has to do with primitives and conservative forms
110 V[Shortcuts::B + 0] = 0.;
111 V[Shortcuts::B + 1] = 0.;
112 V[Shortcuts::B + 2] = 0.;
113 V[Shortcuts::psi] = 0.;
114 V[Shortcuts::lapse] = 1.;
115 V[Shortcuts::shift + 0] = 0.;
116 V[Shortcuts::shift + 1] = 0.;
117 V[Shortcuts::shift + 2] = 0.;
118 V[Shortcuts::gij + 0] = 1.;
119 V[Shortcuts::gij + 1] = 0.;
120 V[Shortcuts::gij + 2] = 0.;
121 V[Shortcuts::gij + 3] = 1.;
122 V[Shortcuts::gij + 4] = 0.;
123 V[Shortcuts::gij + 5] = 1.;
124
125 PDE::PDEPrim2Cons(Q, V);
126
127 return;
128 }
129
130 double theta = acos(x[2] / r);
131 phi = atan2(x[1], x[0]);
132 double zz = 2.0 / r; // we are computing the solution at theta=pi/2
133 double betaru = zz / (1.0 + zz);
134 // double g_tt = zz - 1.0; // unused
135
136 double urc = sqrt(1.0 / (2.0 * rc));
137 double vc2 = urc * urc / (1.0 - 3.0 * urc * urc);
138 double tc = ng * vc2 / ((1.0 + ng) * (1.0 - ng * vc2));
139 // double pc = rhoc*tc; // unused
140
141 // double c1 = urc*tc**ng*rc**2;
142 // double c2 = (1.0 + ( 1.0 + ng)*tc)**2*(1.0 - 2.0/rc+urc*urc);
143 double c1 = urc * std::pow(tc, ng) * (rc * rc);
144 double c2 = (1.0 + (1.0 + ng) * tc) * (1.0 + (1.0 + ng) * tc) * (1.0 - 2.0 / rc + urc * urc);
145
146 double tt = tc;
147 double urr;
148 for (int i = 0; i < MaxNewton; i++) {
149 // urr = c1 / (r**2*tt**ng)
150 urr = c1 / (r * r * std::pow(tt, ng));
151 // f = (1.0 + (1.0 + ng)*tt)**2*(1.0 - 2.0/r + urr**2) - c2
152 double f = (1.0 + (1.0 + ng) * tt) * (1.0 + (1.0 + ng) * tt) * (1.0 - 2.0 / r + urr * urr) - c2;
153 // df = 2.0 * (1.0 + ng)*(1.0 + (1.0 + ng)*tt)*(1.0 - 2.0/r + urr**2) - 2.0*ng*urr**2/tt*(1.0 + (1.0 +
154 // ng)*tt)**2
155 double df = 2.0 * (1.0 + ng) * (1.0 + (1.0 + ng) * tt) * (1.0 - 2.0 / r + urr * urr)
156 - 2.0 * ng * urr * urr / tt * (1.0 + (1.0 + ng) * tt) * (1.0 + (1.0 + ng) * tt);
157 double dt = -f / df;
158 if (std::abs(dt) < 1.e-10) {
159 break;
160 }
161 tt = tt + dt;
162 }
163
164 double ut = (-zz * urr + sqrt(urr * urr - zz + 1.0)) / (zz - 1.0);
165 double LF = lapse * ut;
166 double vr = (urr / LF + betaru / lapse);
167
168 double vx = sin(theta) * cos(phi) * vr;
169 double vy = sin(theta) * sin(phi) * vr;
170 double vz = cos(theta) * vr;
171
172 double VV[3] = {vx, vy, vz};
173
174 // Convert to covariant velocities
175 double VV_cov[3];
176 PDE::matrixVectMult(g_cov, VV, VV_cov);
177
178 double rho = rhoc * std::pow(tt / tc, ng); // rho = rhoc*(tt/tc)**ng
179 double p = std::pow(rho, tt);
180
181 double gammaij[6] = {g_cov[0][0], g_cov[0][1], g_cov[0][2], g_cov[1][1], g_cov[1][2], g_cov[2][2]};
182
183 // !If hydro (Clasical) Michel accretion
184 // !BV(1:3) = 0.
185 // ! MHD mICHEl accretion
186 // ! detgamma= gammaij(1)*( gammaij(4)*gammaij(6)-gammaij(5)*gammaij(5)) &
187 // ! - gammaij(2)*( gammaij(2)*gammaij(6)-gammaij(5)*gammaij(3)) &
188 // ! + gammaij(3)*( gammaij(2)*gammaij(5)-gammaij(3)*gammaij(4))
189 // ! detgamma = sqrt(detgamma)
190 // ! B0 = (2.2688/M)*(sqrt(b^2/rho0)), here sqrt(b^2/rho0)=4.0
191
192 double B0 = 0;
193 double BV_contr[3] = {
194 // BV_contr = B0*x(1:3)/(sqrt(gp) * r*r*r)
195 B0 * x[0] / (sqrt(gp) * r * r * r),
196 B0 * x[1] / (sqrt(gp) * r * r * r),
197 B0 * x[2] / (sqrt(gp) * r * r * r)
198 };
199
200 double BV[3];
201 PDE::matrixVectMult(g_cov, BV_contr, BV);
202
203 using Shortcuts = tests::exahype2::aderdg::VariableShortcutsADERDGSolver;
204 V[Shortcuts::rho] = rho;
205 V[Shortcuts::vel + 0] = VV_cov[0];
206 V[Shortcuts::vel + 1] = VV_cov[1];
207 V[Shortcuts::vel + 2] = VV_cov[2];
208 V[Shortcuts::E] = p; // should be p?, likely has to do with primitives and conservative forms
209 V[Shortcuts::B + 0] = BV[0];
210 V[Shortcuts::B + 1] = BV[1];
211 V[Shortcuts::B + 2] = BV[2];
212 V[Shortcuts::psi] = 0.;
213 V[Shortcuts::lapse] = lapse;
214 V[Shortcuts::shift + 0] = shift[0];
215 V[Shortcuts::shift + 1] = shift[1];
216 V[Shortcuts::shift + 2] = shift[2];
217 V[Shortcuts::gij + 0] = gammaij[0];
218 V[Shortcuts::gij + 1] = gammaij[1];
219 V[Shortcuts::gij + 2] = gammaij[2];
220 V[Shortcuts::gij + 3] = gammaij[3];
221 V[Shortcuts::gij + 4] = gammaij[4];
222 V[Shortcuts::gij + 5] = gammaij[5];
223
224 PDE::PDEPrim2Cons(Q, V);
225 }
226
227 inline void AlfvenWaveTest(const tarch::la::Vector<DIMENSIONS, double>& x, const double t, double* const Q) {
228 // PARAMETERS
229 constexpr double gamma = 4.0 / 3.0;
230
231 constexpr double rho0 = 1.;
232 constexpr double p0 = 1.;
233 constexpr double eta = 1.;
234 constexpr double B0 = 1.;
235
236 double hh = 1.0 + gamma / (gamma - 1.) * p0 / rho0;
237 double tempaa = rho0 * hh + B0 * B0 * (1.0 + eta * eta);
238 double tempab = 2.0 * eta * B0 * B0 / tempaa;
239 double tempac = 0.5 * (1.0 + sqrt(1.0 - tempab * tempab));
240
241 double va2 = B0 * B0 / (tempaa * tempac);
242 double vax = sqrt(va2);
243
244 using Shortcuts = tests::exahype2::aderdg::VariableShortcutsADERDGSolver;
245 double V[19];
246
247 V[Shortcuts::rho] = rho0;
248
249 V[Shortcuts::B + 0] = B0;
250 V[Shortcuts::B + 1] = eta * B0 * cos(x[0] - vax * t);
251 V[Shortcuts::B + 2] = eta * B0 * sin(x[0] - vax * t);
252
253 V[Shortcuts::vel + 0] = 0.0;
254 V[Shortcuts::vel + 1] = -vax * V[Shortcuts::B + 1] / B0;
255 V[Shortcuts::vel + 2] = -vax * V[Shortcuts::B + 2] / B0;
256 V[Shortcuts::E] = p0;
257
258 V[Shortcuts::psi] = 0.;
259 V[Shortcuts::lapse] = 1.0;
260 V[Shortcuts::shift + 0] = 0.;
261 V[Shortcuts::shift + 1] = 0.;
262 V[Shortcuts::shift + 2] = 0.;
263
264 V[Shortcuts::gij + 0] = 1.;
265 V[Shortcuts::gij + 1] = 0.;
266 V[Shortcuts::gij + 2] = 0.;
267 V[Shortcuts::gij + 3] = 1.;
268 V[Shortcuts::gij + 4] = 0.;
269 V[Shortcuts::gij + 5] = 1.;
270
271 PDE::PDEPrim2Cons(Q, V);
272 }
273 } // namespace InitialData
274} // namespace GRMHD
GRMHD::InitialData::initialAccretionDisc3D
void initialAccretionDisc3D(const tarch::la::Vector< DIMENSIONS, double > &x, const double t, double *const Q)
Definition InitialData.h:78
GRMHD::InitialData::metric
void metric(const tarch::la::Vector< DIMENSIONS, double > &x, double *lapse, double *gp, double *gm, double *shift, double **Kex, double g_cov[][3], double g_contr[][3], double *phi)
Definition InitialData.h:10
GRMHD::InitialData::AlfvenWaveTest
void AlfvenWaveTest(const tarch::la::Vector< DIMENSIONS, double > &x, const double t, double *const Q)
Definition InitialData.h:227
GRMHD::PDE::matrixVectMult
void matrixVectMult(const double M[][3], const double VIn[3], double VOut[3])
Definition PDE.h:11
GRMHD::PDE::PDEPrim2Cons
void PDEPrim2Cons(double *Q, const double *V)
Definition PDE.h:200
GRMHD
Definition GRMHD.py:1