Peano
Loading...
Searching...
No Matches
riemannsolverPML.h
Go to the documentation of this file.
1#pragma once
2
3#include "./idx.h"
4#include "curvilinearRoutines.h"
5#include "riemannsolverRoutines.h"
6
7namespace Numerics{
8
9 //Getter routines
10 template <class Shortcuts, typename T>
11 inline void computeParameters(const T*Q, T& rho, T& cp, T& cs, T& mu, T& lam){
12 rho = Q[Shortcuts::rho];
13 cp = Q[Shortcuts::cp ];
14 cs = Q[Shortcuts::cs ];
15 mu = cs * cs * rho;
16 lam = rho* cp * cp - 2.0 * mu;
17 }
18
19 template<typename T>
20 inline void riemannSolverBoundary(int faceIndex,T r, T vn , T vm , T vl, T Tn , T Tm ,T Tl , T zp, T zs , T& vn_hat , T& vm_hat ,T& vl_hat , T& Tn_hat , T& Tm_hat ,T& Tl_hat)
21 {
22 // if (faceIndex % 2 == 0) { //left face
23 if (faceIndex < DIMENSIONS) { //left face
24 riemannSolverBC0(vn, Tn, zp, r, vn_hat, Tn_hat);
25 riemannSolverBC0(vm, Tm, zs, r, vm_hat, Tm_hat);
26 riemannSolverBC0(vl, Tl, zs, r, vl_hat, Tl_hat);
27 }
28
29 // if (faceIndex % 2 == 1) { //right face
30 else { //right face
31 riemannSolverBCn(vn, Tn, zp, r, vn_hat, Tn_hat);
32 riemannSolverBCn(vm, Tm, zs, r, vm_hat, Tm_hat);
33 riemannSolverBCn(vl, Tl, zs, r, vl_hat, Tl_hat);
34 }
35 }
36
37
38 template<class Shortcuts, typename T, int order, int numberOfVariables, int numberOfParameters , int surface = 2>
39 void riemannSolver(T* FL,T* FR,const T* const QL,const T* const QR,const double dt,const int direction, bool isBoundaryFace, int faceIndex){
40
41 constexpr int numberOfData = numberOfVariables+numberOfParameters;
42 constexpr int basisSize = order+1;
43
44 kernels::idx3 idx_QLR(basisSize,basisSize,numberOfData);
45 kernels::idx3 idx_FLR(basisSize,basisSize,numberOfVariables);
46
47 T FLn ,FLm ,FLl ,FRn ,FRm ,FRl;
48 T FLx ,FLy ,FLz ,FRx ,FRy ,FRz;
49 T FL_n,FL_m,FL_l,FR_n,FR_m,FR_l;
50 T FL_x,FL_y,FL_z,FR_x,FR_y,FR_z;
51
52 for (int i = 0; i < basisSize; i++) {
53 for (int j = 0; j < basisSize; j++) {
54 const T* Q_m = QL+idx_QLR(i,j,0);
55 const T* Q_p = QR+idx_QLR(i,j,0);
56
57 T dmp_pml_x_m = Q_m[Shortcuts::dmp_pml + 0];
58 T dmp_pml_y_m = Q_m[Shortcuts::dmp_pml + 1];
59 T dmp_pml_z_m = Q_m[Shortcuts::dmp_pml + 2];
60
61 T dmp_pml_x_p = Q_p[Shortcuts::dmp_pml + 0];
62 T dmp_pml_y_p = Q_p[Shortcuts::dmp_pml + 1];
63 T dmp_pml_z_p = Q_p[Shortcuts::dmp_pml + 2];
64
65 T* F_m = FL+idx_FLR(i,j,0);
66 T* F_p = FR+idx_FLR(i,j,0);
67 T rho_m,cp_m,cs_m,mu_m,lam_m;
68 T rho_p,cp_p,cs_p,mu_p,lam_p;
69
70 computeParameters<Shortcuts>(Q_m,rho_m,cp_m,cs_m,mu_m,lam_m);
71 computeParameters<Shortcuts>(Q_p,rho_p,cp_p,cs_p,mu_p,lam_p);
72
73 T n_m[3],m_m[3],l_m[3];
74 T n_p[3],m_p[3],l_p[3];
75 T norm_p,norm_m;
76
77 getNormals<Shortcuts>(Q_m,direction,norm_m,n_m);
78 getNormals<Shortcuts>(Q_p,direction,norm_p,n_p);
79
80
81 T Tx_m,Ty_m,Tz_m;
82 T Tx_p,Ty_p,Tz_p;
83 computeTractions<Shortcuts>(Q_p,n_p,Tx_p,Ty_p,Tz_p);
84 computeTractions<Shortcuts>(Q_m,n_m,Tx_m,Ty_m,Tz_m);
85
86 T vx_m,vy_m,vz_m;
87 T vx_p,vy_p,vz_p;
88 getVelocities<Shortcuts>(Q_p,vx_p,vy_p,vz_p);
89 getVelocities<Shortcuts>(Q_m,vx_m,vy_m,vz_m);
90
91 createLocalBasis(n_p, m_p, l_p);
92 createLocalBasis(n_m, m_m, l_m);
93
94 T Tn_m,Tm_m,Tl_m;
95 T Tn_p,Tm_p,Tl_p;
96
97 // rotate fields into l, m, n basis
98 rotateIntoOrthogonalBasis(n_m,m_m,l_m,Tx_m,Ty_m,Tz_m,Tn_m,Tm_m,Tl_m);
99 rotateIntoOrthogonalBasis(n_p,m_p,l_p,Tx_p,Ty_p,Tz_p,Tn_p,Tm_p,Tl_p);
100
101 T vn_m,vm_m,vl_m;
102 T vn_p,vm_p,vl_p;
103 rotateIntoOrthogonalBasis(n_m,m_m,l_m,vx_m,vy_m,vz_m,vn_m,vm_m,vl_m);
104 rotateIntoOrthogonalBasis(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 T zs_m=rho_m*cs_m;
109 T zs_p=rho_p*cs_p;
110
111 T zp_m=rho_m*cp_m;
112 T 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 T vn_hat_p,vm_hat_p,vl_hat_p;
120 T Tn_hat_p,Tm_hat_p,Tl_hat_p;
121 T vn_hat_m,vm_hat_m,vl_hat_m;
122 T Tn_hat_m,Tm_hat_m,Tl_hat_m;
123
124 if (isBoundaryFace) {
125 // 0 absorbing 1 free surface
126 T r= faceIndex == surface ? 1 : 0;
127 // double r=1.;
128 riemannSolverBoundary(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 riemannSolverBoundary(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 riemannSolverNodal(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 riemannSolverNodal(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 riemannSolverNodal(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 computeFluctuationsLeft(zp_m,
160 Tn_m,Tn_hat_m,
161 vn_m,vn_hat_m,
162 FLn);
163 computeFluctuationsLeft(zs_m,
164 Tm_m,Tm_hat_m,
165 vm_m,vm_hat_m,
166 FLm);
167 computeFluctuationsLeft(zs_m,
168 Tl_m,Tl_hat_m,
169 vl_m,vl_hat_m,
170 FLl);
171
172 computeFluctuationsRight(zp_p,
173 Tn_p,Tn_hat_p,
174 vn_p,vn_hat_p,
175 FRn);
176 computeFluctuationsRight(zs_p,
177 Tm_p,Tm_hat_p,
178 vm_p,vm_hat_p,
179 FRm);
180 computeFluctuationsRight(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 rotateIntoPhysicalBasis(n_m,m_m,l_m,
208 FLn,FLm,FLl,
209 FLx,FLy,FLz);
210 rotateIntoPhysicalBasis(n_p,m_p,l_p,
211 FRn,FRm,FRl,
212 FRx,FRy,FRz);
213 rotateIntoPhysicalBasis(n_m,m_m,l_m,
214 FL_n,FL_m,FL_l,
215 FL_x,FL_y,FL_z);
216 rotateIntoPhysicalBasis(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[Shortcuts::v + 0] = norm_p/rho_p*FRx;
224 F_m[Shortcuts::v + 0] = norm_m/rho_m*FLx;
225
226 F_p[Shortcuts::v + 1] = norm_p/rho_p*FRy;
227 F_m[Shortcuts::v + 1] = norm_m/rho_m*FLy;
228
229 F_p[Shortcuts::v + 2] = norm_p/rho_p*FRz;
230 F_m[Shortcuts::v + 2] = norm_m/rho_m*FLz;
231
232 F_m[Shortcuts::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[Shortcuts::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[Shortcuts::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[Shortcuts::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[Shortcuts::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[Shortcuts::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[Shortcuts::sigma + 3] = norm_m*mu_m*(n_m[1]*FL_x + n_m[0]*FL_y);
241 F_m[Shortcuts::sigma + 4] = norm_m*mu_m*(n_m[2]*FL_x + n_m[0]*FL_z);
242 F_m[Shortcuts::sigma + 5] = norm_m*mu_m*(n_m[2]*FL_y + n_m[1]*FL_z);
243
244 F_p[Shortcuts::sigma + 3] = -norm_p*mu_p*(n_p[1]*FR_x + n_p[0]*FR_y);
245 F_p[Shortcuts::sigma + 4] = -norm_p*mu_p*(n_p[2]*FR_x + n_p[0]*FR_z);
246 F_p[Shortcuts::sigma + 5] = -norm_p*mu_p*(n_p[2]*FR_y + n_p[1]*FR_z);
247
248 T norm_p_qr=norm_p;
249 T norm_m_qr=norm_m;
250 kernels::idx2 idx_pml(DIMENSIONS,9);
251
252 F_p[Shortcuts::pml + idx_pml(0,0)] = n_p[0]*dmp_pml_x_p*norm_p*FRx;
253 F_m[Shortcuts::pml + idx_pml(0,0)] = n_m[0]*dmp_pml_x_m*norm_m*FLx;
254
255 F_p[Shortcuts::pml + idx_pml(0,1)] = n_p[0]*dmp_pml_x_p*norm_p*FRy;
256 F_m[Shortcuts::pml + idx_pml(0,1)] = n_m[0]*dmp_pml_x_m*norm_m*FLy;
257
258 F_p[Shortcuts::pml + idx_pml(0,2)] = n_p[0]*dmp_pml_x_p*norm_p*FRz;
259 F_m[Shortcuts::pml + idx_pml(0,2)] = n_m[0]*dmp_pml_x_m*norm_m*FLz;
260
261 F_m[Shortcuts::pml + idx_pml(0,3)] = n_m[0]*dmp_pml_x_m*norm_m*n_m[0]*FL_x;
262 F_m[Shortcuts::pml + idx_pml(0,4)] = n_m[0]*dmp_pml_x_m*norm_m*n_m[1]*FL_y;
263 F_m[Shortcuts::pml + idx_pml(0,5)] = n_m[0]*dmp_pml_x_m*norm_m*n_m[2]*FL_z;
264
265 F_p[Shortcuts::pml + idx_pml(0,3)] = -n_p[0]*dmp_pml_x_p*norm_p*n_p[0]*FR_x;
266 F_p[Shortcuts::pml + idx_pml(0,4)] = -n_p[0]*dmp_pml_x_p*norm_p*n_p[1]*FR_y;
267 F_p[Shortcuts::pml + idx_pml(0,5)] = -n_p[0]*dmp_pml_x_p*norm_p*n_p[2]*FR_z;
268
269 F_m[Shortcuts::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[Shortcuts::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[Shortcuts::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[Shortcuts::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[Shortcuts::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[Shortcuts::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[Shortcuts::pml + idx_pml(1,0)] = n_p[1]*dmp_pml_y_p*norm_p*FRx;
279 F_m[Shortcuts::pml + idx_pml(1,0)] = n_m[1]*dmp_pml_y_m*norm_m*FLx;
280
281 F_p[Shortcuts::pml + idx_pml(1,1)] = n_p[1]*dmp_pml_y_p*norm_p*FRy;
282 F_m[Shortcuts::pml + idx_pml(1,1)] = n_m[1]*dmp_pml_y_m*norm_m*FLy;
283
284 F_p[Shortcuts::pml + idx_pml(1,2)] = n_p[1]*dmp_pml_y_p*norm_p*FRz;
285 F_m[Shortcuts::pml + idx_pml(1,2)] = n_m[1]*dmp_pml_y_m*norm_m*FLz;
286
287 F_m[Shortcuts::pml + idx_pml(1,3)] = n_m[1]*dmp_pml_y_m*norm_m*n_m[0]*FL_x;
288 F_m[Shortcuts::pml + idx_pml(1,4)] = n_m[1]*dmp_pml_y_m*norm_m*n_m[1]*FL_y;
289 F_m[Shortcuts::pml + idx_pml(1,5)] = n_m[1]*dmp_pml_y_m*norm_m*n_m[2]*FL_z;
290
291 F_p[Shortcuts::pml + idx_pml(1,3)] = -n_p[1]*dmp_pml_y_p*norm_p*n_p[0]*FR_x;
292 F_p[Shortcuts::pml + idx_pml(1,4)] = -n_p[1]*dmp_pml_y_p*norm_p*n_p[1]*FR_y;
293 F_p[Shortcuts::pml + idx_pml(1,5)] = -n_p[1]*dmp_pml_y_p*norm_p*n_p[2]*FR_z;
294
295 F_m[Shortcuts::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[Shortcuts::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[Shortcuts::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[Shortcuts::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[Shortcuts::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[Shortcuts::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[Shortcuts::pml + idx_pml(2,0)] = n_p[2]*dmp_pml_z_p*norm_p*FRx;
305 F_m[Shortcuts::pml + idx_pml(2,0)] = n_m[2]*dmp_pml_z_m*norm_m*FLx;
306
307 F_p[Shortcuts::pml + idx_pml(2,1)] = n_p[2]*dmp_pml_z_p*norm_p*FRy;
308 F_m[Shortcuts::pml + idx_pml(2,1)] = n_m[2]*dmp_pml_z_m*norm_m*FLy;
309
310 F_p[Shortcuts::pml + idx_pml(2,2)] = n_p[2]*dmp_pml_z_p*norm_p*FRz;
311 F_m[Shortcuts::pml + idx_pml(2,2)] = n_m[2]*dmp_pml_z_m*norm_m*FLz;
312
313 F_m[Shortcuts::pml + idx_pml(2,3)] = n_m[2]*dmp_pml_z_m*norm_m*n_m[0]*FL_x;
314 F_m[Shortcuts::pml + idx_pml(2,4)] = n_m[2]*dmp_pml_z_m*norm_m*n_m[1]*FL_y;
315 F_m[Shortcuts::pml + idx_pml(2,5)] = n_m[2]*dmp_pml_z_m*norm_m*n_m[2]*FL_z;
316
317 F_p[Shortcuts::pml + idx_pml(2,3)] = -n_p[2]*dmp_pml_z_p*norm_p*n_p[0]*FR_x;
318 F_p[Shortcuts::pml + idx_pml(2,4)] = -n_p[2]*dmp_pml_z_p*norm_p*n_p[1]*FR_y;
319 F_p[Shortcuts::pml + idx_pml(2,5)] = -n_p[2]*dmp_pml_z_p*norm_p*n_p[2]*FR_z;
320
321 F_m[Shortcuts::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[Shortcuts::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[Shortcuts::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[Shortcuts::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[Shortcuts::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[Shortcuts::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}
Numerics::computeFluctuationsLeft
void computeFluctuationsLeft(T z, T myT, T T_hat, T v, T v_hat, T &F)
Definition riemannsolverRoutines.h:172
Numerics::riemannSolverBoundary
void riemannSolverBoundary(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 riemannsolverAnisotropic.h:49
Numerics::rotateIntoOrthogonalBasis
void rotateIntoOrthogonalBasis(T *n, T *m, T *l, T Tx, T Ty, T Tz, T &Tn, T &Tm, T &Tl)
Definition riemannsolverRoutines.h:153
Numerics::riemannSolverNodal
void riemannSolverNodal(T v_p, T v_m, T sigma_p, T sigma_m, T z_p, T z_m, T &v_hat_p, T &v_hat_m, T &sigma_hat_p, T &sigma_hat_m)
Definition riemannsolverRoutines.h:81
Numerics::riemannSolverBCn
void riemannSolverBCn(T v, T sigma, T z, T r, T &v_hat, T &sigma_hat)
Definition riemannsolverRoutines.h:141
Numerics::createLocalBasis
void createLocalBasis(T *n, T *m, T *l)
Definition riemannsolverRoutines.h:44
Numerics::computeParameters
void computeParameters(const double *Q, double &rho, double &cp, double cs[], int direction)
Definition riemannsolverAnisotropic.h:25
Numerics::computeFluctuationsRight
void computeFluctuationsRight(T z, T myT, T T_hat, T v, T v_hat, T &F)
Definition riemannsolverRoutines.h:177
Numerics::rotateIntoPhysicalBasis
void rotateIntoPhysicalBasis(T *n, T *m, T *l, T Fn, T Fm, T Fl, T &Fx, T &Fy, T &Fz)
Definition riemannsolverRoutines.h:162
Numerics::riemannSolverBC0
void riemannSolverBC0(T v, T sigma, T z, T r, T &v_hat, T &sigma_hat)
Definition riemannsolverRoutines.h:129