Peano
Loading...
Searching...
No Matches
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
7template <class Shortcuts, int basisSize, int numberOfVariables, int numberOfParameters, typename T>
8void 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::computeParameters
void computeParameters(const double *Q, double &rho, double &cp, double cs[], int direction)
Definition riemannsolverAnisotropic.h:25
Numerics::getVelocities
void getVelocities(const T *Q, T &vx, T &vy, T &vz)
Definition riemannsolverRoutines.h:20
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::computeTractions
void computeTractions(const T *Q, const T *n, T &Tx, T &Ty, T &Tz)
Definition riemannsolverRoutines.h:6
Numerics::getNormals
void getNormals(const T *Q, int direction, T &norm, T *n)
Definition cartesianRoutines.h:5
euler.j
j
Definition euler.py:95