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