Peano
riemannsolverPMLDynamicRupture.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include "../Curvilinear/ContextDynamicRupture.h"
4 #include "riemannsolverPML.h"
5 #include "dynamicRupture.h"
6 
7 template <class Shortcuts, int basisSize, int numberOfVariables, int numberOfParameters, typename T>
8 void ContextDynamicRupture<Shortcuts, basisSize, numberOfVariables, numberOfParameters, T>::riemannSolver(
9  T* FL, T* FR,
10  const T* const QL, const T* const QR,
11  const double t, const double dt,
12  const tarch::la::Vector<DIMENSIONS, double>& cellSize,
13  const int direction,
14  bool isBoundaryFace,
15  int faceIndex,
16  int surface
17 ) {
18  constexpr int numberOfData = numberOfVariables+numberOfParameters;
19  // constexpr int basisSize = order+1;
20 
21  ::kernels::idx3 idx_QLR(basisSize,basisSize,numberOfData);
22  ::kernels::idx3 idx_FLR(basisSize,basisSize,numberOfVariables);
23 
24  //Checking whether the face is on a fault
25  int level = std::round(log(domainSize[0]/cellSize[0])/log(3.)) + 1;
26 
27  int elt_z = int(std::round( (QL[idx_QLR(0,0,Shortcuts::curve_grid+2)] -
28  this->domainOffset[2])/ this->max_dx)) * basisSize;
29  int elt_y = int(std::round((QL[idx_QLR(0,0,Shortcuts::curve_grid+1)] -
30  this->domainOffset[1])/ this->max_dx)) * basisSize;
31 
32  bool is_fault = (direction == 0);
33 
34  toolbox::curvi::Root* root = this->interface->getRoot();
35  toolbox::curvi::InnerNode* fault_node = static_cast<toolbox::curvi::InnerNode*>(root->getChild());
36 
37  toolbox::curvi::Coordinate fault_coords[2];
38  fault_node->getCoordinates(fault_coords);
39  toolbox::curvi::Coordinate fault_normal = fault_node->getFaceNormal();
40  T position = fault_node->getPosition();
41 
42  for (int i = 0; i < basisSize; i++) {
43  for (int j = 0; j < basisSize; j++) {
44  T eta = QL[idx_QLR(i,j,Shortcuts::curve_grid + (2-fault_normal))];
45  T xi = QL[idx_QLR(i,j,Shortcuts::curve_grid + (2-fault_coords[0]))];
46  T mu = QL[idx_QLR(i,j,Shortcuts::curve_grid + (2-fault_coords[1]))];
47 
48  T per_position = position + fault_node->evalPerturbation(xi,mu);
49 
50  is_fault = is_fault && (std::abs(eta - per_position) < cellSize[2-fault_normal] * 0.5);
51  }
52  }
53 
54  T FLn ,FLm ,FLl ,FRn ,FRm ,FRl;
55  T FLx ,FLy ,FLz ,FRx ,FRy ,FRz;
56  T FL_n,FL_m,FL_l,FR_n,FR_m,FR_l;
57  T FL_x,FL_y,FL_z,FR_x,FR_y,FR_z;
58 
59  for (int i = 0; i < basisSize; i++) {
60  for (int j = 0; j < basisSize; j++) {
61 
62  const T* Q_m = QL+idx_QLR(i,j,0);
63  const T* Q_p = QR+idx_QLR(i,j,0);
64 
65  T dmp_pml_x_m = Q_m[Shortcuts::dmp_pml + 0];
66  T dmp_pml_y_m = Q_m[Shortcuts::dmp_pml + 1];
67  T dmp_pml_z_m = Q_m[Shortcuts::dmp_pml + 2];
68 
69  T dmp_pml_x_p = Q_p[Shortcuts::dmp_pml + 0];
70  T dmp_pml_y_p = Q_p[Shortcuts::dmp_pml + 1];
71  T dmp_pml_z_p = Q_p[Shortcuts::dmp_pml + 2];
72 
73  T* F_m = FL + idx_FLR(i,j,0);
74  T* F_p = FR + idx_FLR(i,j,0);
75  T rho_m,cp_m,cs_m,mu_m,lam_m;
76  T rho_p,cp_p,cs_p,mu_p,lam_p;
77 
78  ::Numerics::computeParameters<Shortcuts>(Q_m,rho_m,cp_m,cs_m,mu_m,lam_m);
79  ::Numerics::computeParameters<Shortcuts>(Q_p,rho_p,cp_p,cs_p,mu_p,lam_p);
80 
81  T n_m[3],m_m[3],l_m[3];
82  T n_p[3],m_p[3],l_p[3];
83  T norm_p,norm_m;
84 
85  ::Numerics::getNormals<Shortcuts>(Q_m,direction,norm_m,n_m);
86  ::Numerics::getNormals<Shortcuts>(Q_p,direction,norm_p,n_p);
87 
88  T Tx_m,Ty_m,Tz_m;
89  T Tx_p,Ty_p,Tz_p;
90  ::Numerics::computeTractions<Shortcuts>(Q_p,n_p,Tx_p,Ty_p,Tz_p);
91  ::Numerics::computeTractions<Shortcuts>(Q_m,n_m,Tx_m,Ty_m,Tz_m);
92 
93  T vx_m,vy_m,vz_m;
94  T vx_p,vy_p,vz_p;
95  ::Numerics::getVelocities<Shortcuts>(Q_p,vx_p,vy_p,vz_p);
96  ::Numerics::getVelocities<Shortcuts>(Q_m,vx_m,vy_m,vz_m);
97 
98  ::Numerics::createLocalBasis(n_p, m_p, l_p);
99  ::Numerics::createLocalBasis(n_m, m_m, l_m);
100 
101  T Tn_m,Tm_m,Tl_m;
102  T Tn_p,Tm_p,Tl_p;
103 
104  // rotate fields into l, m, n basis
105  ::Numerics::rotateIntoOrthogonalBasis(n_m,m_m,l_m,Tx_m,Ty_m,Tz_m,Tn_m,Tm_m,Tl_m);
106  ::Numerics::rotateIntoOrthogonalBasis(n_p,m_p,l_p,Tx_p,Ty_p,Tz_p,Tn_p,Tm_p,Tl_p);
107 
108  T vn_m,vm_m,vl_m;
109  T vn_p,vm_p,vl_p;
110  ::Numerics::rotateIntoOrthogonalBasis(n_m,m_m,l_m,vx_m,vy_m,vz_m,vn_m,vm_m,vl_m);
111  ::Numerics::rotateIntoOrthogonalBasis(n_p,m_p,l_p,vx_p,vy_p,vz_p,vn_p,vm_p,vl_p);
112 
113 
114  // extract local s-wave and p-wave impedances
115  T zs_m=rho_m*cs_m;
116  T zs_p=rho_p*cs_p;
117 
118  T zp_m=rho_m*cp_m;
119  T zp_p=rho_p*cp_p;
120 
121  // impedance must be greater than zero !
122  assertion3(!(zp_p <= 0.0 || zp_m <= 0.0),"Impedance must be greater than zero !",zp_p,zs_p);
123 
124  // generate interface data preserving the amplitude of the outgoing charactertritics
125  // and satisfying interface conditions exactly.
126  T vn_hat_p,vm_hat_p,vl_hat_p;
127  T Tn_hat_p,Tm_hat_p,Tl_hat_p;
128  T vn_hat_m,vm_hat_m,vl_hat_m;
129  T Tn_hat_m,Tm_hat_m,Tl_hat_m;
130 
131  if(is_fault){
132 
133  T Sn_m,Sm_m,Sl_m,Sn_p,Sm_p,Sl_p;
134  T Sx_m,Sy_m,Sz_m,Sx_p,Sy_p,Sz_p;
135 
136  double x[3] = {QR[idx_QLR(i,j,Shortcuts::curve_grid + 0)],
137  QR[idx_QLR(i,j,Shortcuts::curve_grid + 1)],
138  QR[idx_QLR(i,j,Shortcuts::curve_grid + 2)]};
139 
140  Sx_p = QR[idx_QLR(i,j,Shortcuts::u + 0)];
141  Sy_p = QR[idx_QLR(i,j,Shortcuts::u + 1)];
142  Sz_p = QR[idx_QLR(i,j,Shortcuts::u + 2)];
143 
144  Sx_m = QL[idx_QLR(i,j,Shortcuts::u + 0)];
145  Sy_m = QL[idx_QLR(i,j,Shortcuts::u + 1)];
146  Sz_m = QL[idx_QLR(i,j,Shortcuts::u + 2)];
147 
148  // tarch::la::Vector<3,double> coords;
149  double coords[3] = {
150  QL[idx_QLR(i,j,Shortcuts::curve_grid + 0 )],
151  QL[idx_QLR(i,j,Shortcuts::curve_grid + 1 )],
152  QL[idx_QLR(i,j,Shortcuts::curve_grid + 2 )]
153  };
154 
155  ::Numerics::rotateIntoOrthogonalBasis(n_m, m_m, l_m, Sx_m, Sy_m, Sz_m, Sn_m, Sm_m, Sl_m);
156  ::Numerics::rotateIntoOrthogonalBasis(n_p, m_p, l_p, Sx_p, Sy_p, Sz_p, Sn_p, Sm_p, Sl_p);
157 
158  T S = std::sqrt((Sl_p- Sl_m)*(Sl_p- Sl_m)+(Sm_p- Sm_m)*(Sm_p- Sm_m));
159 
160  slipWeakeningFriction(
161  vn_p,vn_m, Tn_p,Tn_m, zp_p , zp_m, vn_hat_p , vn_hat_m, Tn_hat_p,Tn_hat_m, vm_p,vm_m,
162  Tm_p,Tm_m, zs_p,zs_m, vm_hat_p, vm_hat_m, Tm_hat_p,Tm_hat_m, vl_p,vl_m,Tl_p,Tl_m, zs_p,
163  zs_m, vl_hat_p , vl_hat_m, Tl_hat_p,Tl_hat_m, l_p, m_p, n_p,
164  // coords.data(),
165  coords,
166  S, t
167  );
168 
169  }
170  else if (isBoundaryFace) {
171  // 0 absorbing 1 free surface
172  T r= faceIndex == surface ? 1 : 0;
173 
174  ::Numerics::riemannSolverBoundary(faceIndex,r,
175  vn_m,vm_m,vl_m,
176  Tn_m,Tm_m,Tl_m,
177  zp_m,zs_m,
178  vn_hat_m,vm_hat_m,vl_hat_m,
179  Tn_hat_m,Tm_hat_m,Tl_hat_m);
180  ::Numerics::riemannSolverBoundary(faceIndex,r,
181  vn_p,vm_p,vl_p,
182  Tn_p,Tm_p,Tl_p,
183  zp_p,zs_p,
184  vn_hat_p,vm_hat_p,vl_hat_p,
185  Tn_hat_p,Tm_hat_p,Tl_hat_p);
186  }
187  else {
188  ::Numerics::riemannSolverNodal(vn_p, vn_m,
189  Tn_p, Tn_m,
190  zp_p , zp_m,
191  vn_hat_p , vn_hat_m,
192  Tn_hat_p, Tn_hat_m);
193  ::Numerics::riemannSolverNodal(vm_p, vm_m,
194  Tm_p, Tm_m,
195  zs_p , zs_m,
196  vm_hat_p , vm_hat_m,
197  Tm_hat_p, Tm_hat_m);
198  ::Numerics::riemannSolverNodal(vl_p, vl_m,
199  Tl_p, Tl_m,
200  zs_p , zs_m,
201  vl_hat_p , vl_hat_m,
202  Tl_hat_p, Tl_hat_m);
203  }
204 
205  //generate fluctuations in the local basis coordinates: n, m, l
206  ::Numerics::computeFluctuationsLeft(zp_m,
207  Tn_m,Tn_hat_m,
208  vn_m,vn_hat_m,
209  FLn);
210  ::Numerics::computeFluctuationsLeft(zs_m,
211  Tm_m,Tm_hat_m,
212  vm_m,vm_hat_m,
213  FLm);
214  ::Numerics::computeFluctuationsLeft(zs_m,
215  Tl_m,Tl_hat_m,
216  vl_m,vl_hat_m,
217  FLl);
218 
219  ::Numerics::computeFluctuationsRight(zp_p,
220  Tn_p,Tn_hat_p,
221  vn_p,vn_hat_p,
222  FRn);
223  ::Numerics::computeFluctuationsRight(zs_p,
224  Tm_p,Tm_hat_p,
225  vm_p,vm_hat_p,
226  FRm);
227  ::Numerics::computeFluctuationsRight(zs_p,
228  Tl_p,Tl_hat_p,
229  vl_p,vl_hat_p,
230  FRl);
231 
232  //Consider acoustic boundary
233  FL_n = FLn/zp_m;
234  if(zs_m > 0){
235  FL_m = FLm/zs_m;
236  FL_l = FLl/zs_m;
237  }else{
238  FL_m=0;
239  FL_l=0;
240  }
241 
242  FR_n = FRn/zp_p;
243  if(zs_p > 0){
244  FR_m = FRm/zs_p;
245  FR_l = FRl/zs_p;
246  }else{
247  FR_m=0;
248  FR_l=0;
249  }
250 
251  // rotate back to the physical coordinates x, y, z
252  ::Numerics::rotateIntoPhysicalBasis(n_m,m_m,l_m,
253  FLn,FLm,FLl,
254  FLx,FLy,FLz);
255  ::Numerics::rotateIntoPhysicalBasis(n_p,m_p,l_p,
256  FRn,FRm,FRl,
257  FRx,FRy,FRz);
258  ::Numerics::rotateIntoPhysicalBasis(n_m,m_m,l_m,
259  FL_n,FL_m,FL_l,
260  FL_x,FL_y,FL_z);
261  ::Numerics::rotateIntoPhysicalBasis(n_p,m_p,l_p,
262  FR_n,FR_m,FR_l,
263  FR_x,FR_y,FR_z);
264 
265  // construct flux fluctuation vectors obeying the eigen structure of the PDE
266  // and choose physically motivated penalties such that we can prove
267  // numerical stability.
268 
269  F_p[Shortcuts::v + 0] = norm_p/rho_p*FRx;
270  F_m[Shortcuts::v + 0] = norm_m/rho_m*FLx;
271 
272  F_p[Shortcuts::v + 1] = norm_p/rho_p*FRy;
273  F_m[Shortcuts::v + 1] = norm_m/rho_m*FLy;
274 
275  F_p[Shortcuts::v + 2] = norm_p/rho_p*FRz;
276  F_m[Shortcuts::v + 2] = norm_m/rho_m*FLz;
277 
278  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);
279  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);
280  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);
281 
282  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);
283  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);
284  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);
285 
286  F_m[Shortcuts::sigma + 3] = norm_m*mu_m*(n_m[1]*FL_x + n_m[0]*FL_y);
287  F_m[Shortcuts::sigma + 4] = norm_m*mu_m*(n_m[2]*FL_x + n_m[0]*FL_z);
288  F_m[Shortcuts::sigma + 5] = norm_m*mu_m*(n_m[2]*FL_y + n_m[1]*FL_z);
289 
290  F_p[Shortcuts::sigma + 3] = -norm_p*mu_p*(n_p[1]*FR_x + n_p[0]*FR_y);
291  F_p[Shortcuts::sigma + 4] = -norm_p*mu_p*(n_p[2]*FR_x + n_p[0]*FR_z);
292  F_p[Shortcuts::sigma + 5] = -norm_p*mu_p*(n_p[2]*FR_y + n_p[1]*FR_z);
293 
294  F_m[Shortcuts::u + 0] = 0;
295  F_m[Shortcuts::u + 1] = 0;
296  F_m[Shortcuts::u + 2] = 0;
297 
298  F_p[Shortcuts::u + 0] = 0;
299  F_p[Shortcuts::u + 1] = 0;
300  F_p[Shortcuts::u + 2] = 0;
301 
302  T norm_p_qr=norm_p;
303  T norm_m_qr=norm_m;
304  ::kernels::idx2 idx_pml(DIMENSIONS,9);
305 
306  F_p[Shortcuts::pml + idx_pml(0,0)] = n_p[0]*dmp_pml_x_p*norm_p*FRx;
307  F_m[Shortcuts::pml + idx_pml(0,0)] = n_m[0]*dmp_pml_x_m*norm_m*FLx;
308 
309  F_p[Shortcuts::pml + idx_pml(0,1)] = n_p[0]*dmp_pml_x_p*norm_p*FRy;
310  F_m[Shortcuts::pml + idx_pml(0,1)] = n_m[0]*dmp_pml_x_m*norm_m*FLy;
311 
312  F_p[Shortcuts::pml + idx_pml(0,2)] = n_p[0]*dmp_pml_x_p*norm_p*FRz;
313  F_m[Shortcuts::pml + idx_pml(0,2)] = n_m[0]*dmp_pml_x_m*norm_m*FLz;
314 
315  F_m[Shortcuts::pml + idx_pml(0,3)] = n_m[0]*dmp_pml_x_m*norm_m*n_m[0]*FL_x;
316  F_m[Shortcuts::pml + idx_pml(0,4)] = n_m[0]*dmp_pml_x_m*norm_m*n_m[1]*FL_y;
317  F_m[Shortcuts::pml + idx_pml(0,5)] = n_m[0]*dmp_pml_x_m*norm_m*n_m[2]*FL_z;
318 
319  F_p[Shortcuts::pml + idx_pml(0,3)] = -n_p[0]*dmp_pml_x_p*norm_p*n_p[0]*FR_x;
320  F_p[Shortcuts::pml + idx_pml(0,4)] = -n_p[0]*dmp_pml_x_p*norm_p*n_p[1]*FR_y;
321  F_p[Shortcuts::pml + idx_pml(0,5)] = -n_p[0]*dmp_pml_x_p*norm_p*n_p[2]*FR_z;
322 
323  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);
324  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);
325  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);
326 
327  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);
328  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);
329  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);
330 
331  // fs.pml + idx_pml(0,0)or y-direction auxiliary variables
332  F_p[Shortcuts::pml + idx_pml(1,0)] = n_p[1]*dmp_pml_y_p*norm_p*FRx;
333  F_m[Shortcuts::pml + idx_pml(1,0)] = n_m[1]*dmp_pml_y_m*norm_m*FLx;
334 
335  F_p[Shortcuts::pml + idx_pml(1,1)] = n_p[1]*dmp_pml_y_p*norm_p*FRy;
336  F_m[Shortcuts::pml + idx_pml(1,1)] = n_m[1]*dmp_pml_y_m*norm_m*FLy;
337 
338  F_p[Shortcuts::pml + idx_pml(1,2)] = n_p[1]*dmp_pml_y_p*norm_p*FRz;
339  F_m[Shortcuts::pml + idx_pml(1,2)] = n_m[1]*dmp_pml_y_m*norm_m*FLz;
340 
341  F_m[Shortcuts::pml + idx_pml(1,3)] = n_m[1]*dmp_pml_y_m*norm_m*n_m[0]*FL_x;
342  F_m[Shortcuts::pml + idx_pml(1,4)] = n_m[1]*dmp_pml_y_m*norm_m*n_m[1]*FL_y;
343  F_m[Shortcuts::pml + idx_pml(1,5)] = n_m[1]*dmp_pml_y_m*norm_m*n_m[2]*FL_z;
344 
345  F_p[Shortcuts::pml + idx_pml(1,3)] = -n_p[1]*dmp_pml_y_p*norm_p*n_p[0]*FR_x;
346  F_p[Shortcuts::pml + idx_pml(1,4)] = -n_p[1]*dmp_pml_y_p*norm_p*n_p[1]*FR_y;
347  F_p[Shortcuts::pml + idx_pml(1,5)] = -n_p[1]*dmp_pml_y_p*norm_p*n_p[2]*FR_z;
348 
349  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);
350  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);
351  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);
352 
353  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);
354  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);
355  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);
356 
357  // fs.pml + idx_pml(0,0)or z-direction auxiliary variables
358  F_p[Shortcuts::pml + idx_pml(2,0)] = n_p[2]*dmp_pml_z_p*norm_p*FRx;
359  F_m[Shortcuts::pml + idx_pml(2,0)] = n_m[2]*dmp_pml_z_m*norm_m*FLx;
360 
361  F_p[Shortcuts::pml + idx_pml(2,1)] = n_p[2]*dmp_pml_z_p*norm_p*FRy;
362  F_m[Shortcuts::pml + idx_pml(2,1)] = n_m[2]*dmp_pml_z_m*norm_m*FLy;
363 
364  F_p[Shortcuts::pml + idx_pml(2,2)] = n_p[2]*dmp_pml_z_p*norm_p*FRz;
365  F_m[Shortcuts::pml + idx_pml(2,2)] = n_m[2]*dmp_pml_z_m*norm_m*FLz;
366 
367  F_m[Shortcuts::pml + idx_pml(2,3)] = n_m[2]*dmp_pml_z_m*norm_m*n_m[0]*FL_x;
368  F_m[Shortcuts::pml + idx_pml(2,4)] = n_m[2]*dmp_pml_z_m*norm_m*n_m[1]*FL_y;
369  F_m[Shortcuts::pml + idx_pml(2,5)] = n_m[2]*dmp_pml_z_m*norm_m*n_m[2]*FL_z;
370 
371  F_p[Shortcuts::pml + idx_pml(2,3)] = -n_p[2]*dmp_pml_z_p*norm_p*n_p[0]*FR_x;
372  F_p[Shortcuts::pml + idx_pml(2,4)] = -n_p[2]*dmp_pml_z_p*norm_p*n_p[1]*FR_y;
373  F_p[Shortcuts::pml + idx_pml(2,5)] = -n_p[2]*dmp_pml_z_p*norm_p*n_p[2]*FR_z;
374 
375  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);
376  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);
377  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);
378 
379  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);
380  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);
381  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);
382  }
383  }
384 }
cellSize
const tarch::la::Vector< DIMENSIONS, double > cellSize
Definition: KernelBenchmarks-main.cpp:51
ContextDynamicRupture::riemannSolver
void riemannSolver(T *FL, T *FR, const T *const QL, const T *const QR, const double t, const double dt, const tarch::la::Vector< DIMENSIONS, double > &cellSize, const int direction, bool isBoundaryFace, int faceIndex, int surface=2)
Definition: riemannsolverDynamicRupture.h:8
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::createLocalBasis
void createLocalBasis(T *n, T *m, T *l)
Definition: riemannsolverRoutines.h:44
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
euler.u
u
Definition: euler.py:113
euler.j
j
Definition: euler.py:95