Peano
Loading...
Searching...
No Matches
riemannsolver_isotropic.h
Go to the documentation of this file.
1// #include "../../../ExaHyPE/kernels/KernelUtils.h"
2// #include "exahype2/aderdg/kernels/Utils/KernelUtils.h"
3
4namespace Numerics {
5
6 // Getter routines
7 template <class Shortcuts>
8 inline void compute_parameters(const double* Q, double& rho, double& cp, double& cs, double& mu, double& lam) {
9 Shortcuts s;
10 rho = Q[s.rho];
11 cp = Q[s.cp];
12 cs = Q[s.cs];
13 mu = cs * cs * rho;
14 lam = rho * cp * cp - 2.0 * mu;
15 }
16
17 inline void riemannSolver_boundary(
18 int faceIndex,
19 double r,
20 double vn,
21 double vm,
22 double vl,
23 double Tn,
24 double Tm,
25 double Tl,
26 double zp,
27 double zs,
28 double& vn_hat,
29 double& vm_hat,
30 double& vl_hat,
31 double& Tn_hat,
32 double& Tm_hat,
33 double& Tl_hat
34 ) {
35 // if (faceIndex % 2 == 0) { //left face
36 if (faceIndex < Dimensions) { //left face
37 riemannSolver_BC0(vn, Tn, zp, r, vn_hat, Tn_hat);
38 riemannSolver_BC0(vm, Tm, zs, r, vm_hat, Tm_hat);
39 riemannSolver_BC0(vl, Tl, zs, r, vl_hat, Tl_hat);
40 }
41
42 // if (faceIndex % 2 == 1) { //right face
43 else { //right face
44 riemannSolver_BCn(vn, Tn, zp, r, vn_hat, Tn_hat);
45 riemannSolver_BCn(vm, Tm, zs, r, vm_hat, Tm_hat);
46 riemannSolver_BCn(vl, Tl, zs, r, vl_hat, Tl_hat);
47 }
48 }
49
50 template <class Shortcuts, int basisSize, int numberOfVariables, int numberOfData, int surface = 2>
51 void riemannSolver(
52 double* FL,
53 double* FR,
54 const double* const QL,
55 const double* const QR,
56 const double dt,
57 const int direction,
58 bool isBoundaryFace,
59 int faceIndex
60 ) {
61
62 kernels::idx3 idx_QLR(basisSize, basisSize, numberOfData);
63 kernels::idx3 idx_FLR(basisSize, basisSize, numberOfVariables);
64
65 double FLn, FLm, FLl, FRn, FRm, FRl;
66 double FLx, FLy, FLz, FRx, FRy, FRz;
67 double FL_n, FL_m, FL_l, FR_n, FR_m, FR_l;
68 double FL_x, FL_y, FL_z, FR_x, FR_y, FR_z;
69
70 Shortcuts s;
71 for (int i = 0; i < basisSize; i++) {
72#pragma simd
73 for (int j = 0; j < basisSize; j++) {
74
75 const double* Q_p = QR + idx_QLR(i, j, 0);
76 const double* Q_m = QL + idx_QLR(i, j, 0);
77
78 double* F_m = FL + idx_FLR(i, j, 0);
79 double* F_p = FR + idx_FLR(i, j, 0);
80 double rho_m, cp_m, cs_m, mu_m, lam_m;
81 double rho_p, cp_p, cs_p, mu_p, lam_p;
82
83 compute_parameters<Shortcuts>(Q_m, rho_m, cp_m, cs_m, mu_m, lam_m);
84 compute_parameters<Shortcuts>(Q_p, rho_p, cp_p, cs_p, mu_p, lam_p);
85
86 double n_m[3], m_m[3], l_m[3];
87 double n_p[3], m_p[3], l_p[3];
88 double norm_p, norm_m;
89
90 get_normals<Shortcuts>(Q_m, direction, norm_m, n_m);
91 get_normals<Shortcuts>(Q_p, direction, norm_p, n_p);
92
93 double Tx_m, Ty_m, Tz_m;
94 double Tx_p, Ty_p, Tz_p;
95 compute_tractions<Shortcuts>(Q_p, n_p, Tx_p, Ty_p, Tz_p);
96 compute_tractions<Shortcuts>(Q_m, n_m, Tx_m, Ty_m, Tz_m);
97
98 double vx_m, vy_m, vz_m;
99 double vx_p, vy_p, vz_p;
100 get_velocities<Shortcuts>(Q_p, vx_p, vy_p, vz_p);
101 get_velocities<Shortcuts>(Q_m, vx_m, vy_m, vz_m);
102
103 create_local_basis(n_p, m_p, l_p);
104 create_local_basis(n_m, m_m, l_m);
105
106 double Tn_m, Tm_m, Tl_m;
107 double Tn_p, Tm_p, Tl_p;
108
109 // rotate fields into l, m, n basis
110 rotate_into_orthogonal_basis(n_m, m_m, l_m, Tx_m, Ty_m, Tz_m, Tn_m, Tm_m, Tl_m);
111 rotate_into_orthogonal_basis(n_p, m_p, l_p, Tx_p, Ty_p, Tz_p, Tn_p, Tm_p, Tl_p);
112
113 double vn_m, vm_m, vl_m;
114 double vn_p, vm_p, vl_p;
115 rotate_into_orthogonal_basis(n_m, m_m, l_m, vx_m, vy_m, vz_m, vn_m, vm_m, vl_m);
116 rotate_into_orthogonal_basis(n_p, m_p, l_p, vx_p, vy_p, vz_p, vn_p, vm_p, vl_p);
117
118 // extract local s-wave and p-wave impedances
119 double zs_m = rho_m * cs_m;
120 double zs_p = rho_p * cs_p;
121
122 double zp_m = rho_m * cp_m;
123 double zp_p = rho_p * cp_p;
124
125 // impedance must be greater than zero !
126 assertion3(!(zp_p <= 0.0 || zp_m <= 0.0), "Impedance must be greater than zero !", zp_p, zp_m);
127
128 // generate interface data preserving the amplitude of the outgoing charactertritics
129 // and satisfying interface conditions exactly.
130 double vn_hat_p, vm_hat_p, vl_hat_p;
131 double Tn_hat_p, Tm_hat_p, Tl_hat_p;
132 double vn_hat_m, vm_hat_m, vl_hat_m;
133 double Tn_hat_m, Tm_hat_m, Tl_hat_m;
134
135 if (isBoundaryFace) {
136 // 0 absorbing 1 free surface
137 double r = faceIndex == surface ? 1 : 0;
138 // double r=1.;
139 riemannSolver_boundary(
140 faceIndex,
141 r,
142 vn_m,
143 vm_m,
144 vl_m,
145 Tn_m,
146 Tm_m,
147 Tl_m,
148 zp_m,
149 zs_m,
150 vn_hat_m,
151 vm_hat_m,
152 vl_hat_m,
153 Tn_hat_m,
154 Tm_hat_m,
155 Tl_hat_m
156 );
157 riemannSolver_boundary(
158 faceIndex,
159 r,
160 vn_p,
161 vm_p,
162 vl_p,
163 Tn_p,
164 Tm_p,
165 Tl_p,
166 zp_p,
167 zs_p,
168 vn_hat_p,
169 vm_hat_p,
170 vl_hat_p,
171 Tn_hat_p,
172 Tm_hat_p,
173 Tl_hat_p
174 );
175 } else {
176 riemannSolver_Nodal(vn_p, vn_m, Tn_p, Tn_m, zp_p, zp_m, vn_hat_p, vn_hat_m, Tn_hat_p, Tn_hat_m);
177 riemannSolver_Nodal(vm_p, vm_m, Tm_p, Tm_m, zs_p, zs_m, vm_hat_p, vm_hat_m, Tm_hat_p, Tm_hat_m);
178 riemannSolver_Nodal(vl_p, vl_m, Tl_p, Tl_m, zs_p, zs_m, vl_hat_p, vl_hat_m, Tl_hat_p, Tl_hat_m);
179 }
180
181 // generate fluctuations in the local basis coordinates: n, m, l
182 compute_fluctuations_left(zp_m, Tn_m, Tn_hat_m, vn_m, vn_hat_m, FLn);
183 compute_fluctuations_left(zs_m, Tm_m, Tm_hat_m, vm_m, vm_hat_m, FLm);
184 compute_fluctuations_left(zs_m, Tl_m, Tl_hat_m, vl_m, vl_hat_m, FLl);
185
186 compute_fluctuations_right(zp_p, Tn_p, Tn_hat_p, vn_p, vn_hat_p, FRn);
187 compute_fluctuations_right(zs_p, Tm_p, Tm_hat_p, vm_p, vm_hat_p, FRm);
188 compute_fluctuations_right(zs_p, Tl_p, Tl_hat_p, vl_p, vl_hat_p, FRl);
189
190 // Consider acoustic boundary
191 FL_n = FLn / zp_m;
192 if (zs_m > 0) {
193 FL_m = FLm / zs_m;
194 FL_l = FLl / zs_m;
195 } else {
196 FL_m = 0;
197 FL_l = 0;
198 }
199
200 FR_n = FRn / zp_p;
201 if (zs_p > 0) {
202 FR_m = FRm / zs_p;
203 FR_l = FRl / zs_p;
204 } else {
205 FR_m = 0;
206 FR_l = 0;
207 }
208
209 // rotate back to the physical coordinates x, y, z
210 rotate_into_physical_basis(n_m, m_m, l_m, FLn, FLm, FLl, FLx, FLy, FLz);
211 rotate_into_physical_basis(n_p, m_p, l_p, FRn, FRm, FRl, FRx, FRy, FRz);
212 rotate_into_physical_basis(n_m, m_m, l_m, FL_n, FL_m, FL_l, FL_x, FL_y, FL_z);
213 rotate_into_physical_basis(n_p, m_p, l_p, FR_n, FR_m, FR_l, FR_x, FR_y, FR_z);
214
215 // construct flux fluctuation vectors obeying the eigen structure of the PDE
216 // and choose physically motivated penalties such that we can prove
217 // numerical stability.
218 F_p[s.v + 0] = norm_p / rho_p * FRx;
219 F_m[s.v + 0] = norm_m / rho_m * FLx;
220
221 F_p[s.v + 1] = norm_p / rho_p * FRy;
222 F_m[s.v + 1] = norm_m / rho_m * FLy;
223
224 F_p[s.v + 2] = norm_p / rho_p * FRz;
225 F_m[s.v + 2] = norm_m / rho_m * FLz;
226
227 F_m[s.sigma + 0] = norm_m * ((2 * mu_m + lam_m) * n_m[0] * FL_x + lam_m * n_m[1] * FL_y + lam_m * n_m[2] * FL_z);
228 F_m[s.sigma + 1] = norm_m * ((2 * mu_m + lam_m) * n_m[1] * FL_y + lam_m * n_m[0] * FL_x + lam_m * n_m[2] * FL_z);
229 F_m[s.sigma + 2] = norm_m * ((2 * mu_m + lam_m) * n_m[2] * FL_z + lam_m * n_m[0] * FL_x + lam_m * n_m[1] * FL_y);
230
231 F_p[s.sigma + 0] = -norm_p
232 * ((2 * mu_p + lam_p) * n_p[0] * FR_x + lam_p * n_p[1] * FR_y + lam_p * n_p[2] * FR_z);
233 F_p[s.sigma + 1] = -norm_p
234 * ((2 * mu_p + lam_p) * n_p[1] * FR_y + lam_p * n_p[0] * FR_x + lam_p * n_p[2] * FR_z);
235 F_p[s.sigma + 2] = -norm_p
236 * ((2 * mu_p + lam_p) * n_p[2] * FR_z + lam_p * n_p[0] * FR_x + lam_p * n_p[1] * FR_y);
237
238 F_m[s.sigma + 3] = norm_m * mu_m * (n_m[1] * FL_x + n_m[0] * FL_y);
239 F_m[s.sigma + 4] = norm_m * mu_m * (n_m[2] * FL_x + n_m[0] * FL_z);
240 F_m[s.sigma + 5] = norm_m * mu_m * (n_m[2] * FL_y + n_m[1] * FL_z);
241
242 F_p[s.sigma + 3] = -norm_p * mu_p * (n_p[1] * FR_x + n_p[0] * FR_y);
243 F_p[s.sigma + 4] = -norm_p * mu_p * (n_p[2] * FR_x + n_p[0] * FR_z);
244 F_p[s.sigma + 5] = -norm_p * mu_p * (n_p[2] * FR_y + n_p[1] * FR_z);
245 }
246 }
247 }
248
249 /*
250 template<class Shortcuts,int basisSize, int numberOfVariables, int numberOfData, int surface = 2>
251 void riemannSolver_vect(double* FL,double* FR,const double* const QL,const double* const QR,const double dt,const
252 int direction, bool isBoundaryFace, int faceIndex){
253
254 Shortcuts s;
255 constexpr int num_nodes = basisSize*basisSize;
256
257 double j_m[num_nodes][3][3];
258 double j_p[num_nodes][3][3];
259
260 kernels::idx2 idx_normal(num_nodes,3);
261 double norm_p[num_nodes],norm_m[num_nodes];
262
263 for (int n = 0; n < num_nodes; n++) {
264 const double* Q_p = QR+idx_QLR(i,j,0);
265 const double* Q_m = QL+idx_QLR(i,j,0);
266
267 norm_p[n] = std::sqrt(Q_p[s.metric_derivative + 0 + direction * 3] * Q_p[s.metric_derivative + 0 + direction * 3]
268 + Q_p[s.metric_derivative + 1 + direction * 3] * Q_p[s.metric_derivative + 1 + direction * 3] +
269 Q_p[s.metric_derivative + 2 + direction * 3] * Q_p[s.metric_derivative + 2 + direction * 3]);
270 j_p[n][0][0] = Q_p[s.metric_derivative + 0 + direction * 3];
271 j_p[n][0][1] = Q_p[s.metric_derivative + 1 + direction * 3];
272 j_p[n][0][2] = Q_p[s.metric_derivative + 2 + direction * 3];
273
274 norm_m[n] = std::sqrt(Q_m[s.metric_derivative + 0 + direction * 3] * Q_m[s.metric_derivative + 0 + direction * 3]
275 + Q_m[s.metric_derivative + 1 + direction * 3] * Q_m[s.metric_derivative + 1 + direction * 3] +
276 Q_m[s.metric_derivative + 2 + direction * 3] * Q_m[s.metric_derivative + 2 + direction * 3]);
277
278 j_m[n][0][0] = Q_m[s.metric_derivative + 0 + direction * 3];
279 j_m[n][0][1] = Q_m[s.metric_derivative + 1 + direction * 3];
280 j_m[n][0][2] = Q_m[s.metric_derivative + 2 + direction * 3];
281 }
282
283
284 //compute tractions
285 double Tractions_arr[*3];
286 double* Tractions[3];
287 Tractions[0] = Tractions_arr;
288 Tractions[1] = Tractions_arr + basisSize * basisSize;
289 Tractions[2] = Tractions_arr + basisSize * basisSize * 2;
290
291 for (int i = 0; i < basisSize; i++) {
292 for (int j = 0; j < basisSize; j++) {
293 const double* Q_p = QR+idx_QLR(i,j,0);
294 const double* Q_m = QL+idx_QLR(i,j,0);
295
296
297
298
299 }
300 }
301 }
302 */
303} // namespace Numerics
assertion3
#define assertion3(expr, param0, param1, param2)
Definition Assertions.h:1020
Numerics::compute_fluctuations_left
void compute_fluctuations_left(double z, double T, double T_hat, double v, double v_hat, double &F)
Definition riemannsolver_routines.h:193
Numerics::get_stiffness_tensor
void get_stiffness_tensor(const double *Q, double &c11, double &c22, double &c33, double &c44, double &c55, double &c66, double &c12, double &c13, double &c23)
Definition anisotropic.h:4
Numerics::rotate_into_orthogonal_basis
void rotate_into_orthogonal_basis(double *n, double *m, double *l, double Tx, double Ty, double Tz, double &Tn, double &Tm, double &Tl)
Definition riemannsolver_routines.h:176
Numerics::riemannSolver_BC0
void riemannSolver_BC0(double v, double sigma, double z, double r, double &v_hat, double &sigma_hat)
Definition riemannsolver_routines.h:154
Numerics::riemannSolver_BCn
void riemannSolver_BCn(double v, double sigma, double z, double r, double &v_hat, double &sigma_hat)
Definition riemannsolver_routines.h:165
Numerics::riemannSolver_boundary
void riemannSolver_boundary(int faceIndex, double r, double vn, double vm, double vl, double Tn, double Tm, double Tl, double zp, double zs[2], double &vn_hat, double &vm_hat, double &vl_hat, double &Tn_hat, double &Tm_hat, double &Tl_hat)
Definition riemannsolver_anisotropic.h:32
Numerics::compute_parameters
void compute_parameters(const double *Q, double &rho, double &cp, double cs[], int direction)
Definition riemannsolver_anisotropic.h:7
Numerics::riemannSolver_Nodal
void riemannSolver_Nodal(double v_p, double v_m, double sigma_p, double sigma_m, double z_p, double z_m, double &v_hat_p, double &v_hat_m, double &sigma_hat_p, double &sigma_hat_m)
Definition riemannsolver_routines.h:97
Numerics::compute_fluctuations_right
void compute_fluctuations_right(double z, double T, double T_hat, double v, double v_hat, double &F)
Definition riemannsolver_routines.h:197
Numerics::riemannSolver
void riemannSolver(double *FL, double *FR, const double *const QL, const double *const QR, const double dt, const int direction, bool isBoundaryFace, int faceIndex)
Definition riemannsolver_anisotropic.h:66
Numerics::create_local_basis
void create_local_basis(double *n, double *m, double *l)
Definition riemannsolver_routines.h:61
Numerics::rotate_into_physical_basis
void rotate_into_physical_basis(double *n, double *m, double *l, double Fn, double Fm, double Fl, double &Fx, double &Fy, double &Fz)
Definition riemannsolver_routines.h:184