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