Peano
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 
7 namespace 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 }
ModeCalc.zs
float zs
Definition: ModeCalc.py:147
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::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: riemannsolverAnisotropic.h:83
Numerics::riemannSolverBC0
void riemannSolverBC0(T v, T sigma, T z, T r, T &v_hat, T &sigma_hat)
Definition: riemannsolverRoutines.h:129
euler.j
j
Definition: euler.py:95