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