Peano 4
Loading...
Searching...
No Matches
riemannSolverLOH.cpp
Go to the documentation of this file.
1
2#include "riemannSolverLOH.h"
3
4examples::exahype2::elastic::VariableShortcuts sh;
5
6template<typename T>
8 T* const FL_,
9 T* const FR_,
10 const T* const QL_,
11 const T* const QR_,
12 const double t,
13 const double dt,
15 const int normalNonZeroIndex,
16 bool isBoundaryFace,
17 int faceIndex){
18
19#ifdef OPT_KERNELS
20 T FL[converter::getFFaceGenArraySize()];
21 T FR[converter::getFFaceGenArraySize()];
22 T QL[converter::getQFaceGenArraySize()];
23 T QR[converter::getQFaceGenArraySize()];
24
25 converter::FFace_optimised2generic(FL_,FL);
26 converter::FFace_optimised2generic(FR_,FR);
27 converter::QFace_optimised2generic(QL_,QL);
28 converter::QFace_optimised2generic(QR_,QR);
29#else
30 T* FL=FL_;
31 T* FR=FR_;
32 const T* QL=QL_;
33 const T* QR=QR_;
34#endif
35
36 constexpr int numberOfVariables = examples::exahype2::elastic::Abstractloh::NumberOfVariables;
37 constexpr int numberOfVariables2 = numberOfVariables*numberOfVariables;
38 constexpr int numberOfParameters = examples::exahype2::elastic::Abstractloh::NumberOfParameters;
39 constexpr int numberOfData = numberOfVariables+numberOfParameters;
40 constexpr int basisSize = examples::exahype2::elastic::Abstractloh::Order+1;
41 constexpr int order = basisSize - 1;
42
43 kernels::idx3 idx_QLR(basisSize,basisSize,numberOfData);
44 kernels::idx3 idx_FLR(basisSize,basisSize,numberOfVariables);
45
46 T n[3]={0,0,0};
47 n[normalNonZeroIndex]=1;
48
49 T n_p[3]={0,0,0};
50 T n_m[3]={0,0,0};
51
52 T m_p[3]={0,0,0};
53 T m_m[3]={0,0,0};
54
55 T l_p[3]={0,0,0};
56 T l_m[3]={0,0,0};
57
58 T norm_p_qr = 1.0;
59 T norm_m_qr = 1.0;
60
61 T FLn, FLm, FLl, FRn,FRm,FRl;
62 T FL_n,FL_m,FL_l,FR_n,FR_m,FR_l;
63 T FLx,FLy,FLz,FRx,FRy,FRz;
64 T FL_x,FL_y,FL_z,FR_x,FR_y,FR_z;
65
66 for (int k = 0; k < 3; k++){
67
68 n_m[k] = n[k];
69 n_p[k] = n[k];
70 }
71
72 for (int i = 0; i < basisSize; i++) {
73 for (int j = 0; j < basisSize; j++) {
74 T rho_p=QR[idx_QLR(i,j,9)];
75 T c_p_p=QR[idx_QLR(i,j,10)];
76 T c_s_p=QR[idx_QLR(i,j,11)];
77 T mu_p=c_s_p*c_s_p*rho_p;
78 T lam_p = rho_p*c_p_p*c_p_p-2*mu_p;
79
80 T rho_m=QL[idx_QLR(i,j,9)];
81 T c_p_m=QL[idx_QLR(i,j,10)];
82 T c_s_m=QL[idx_QLR(i,j,11)];
83 T mu_m=c_s_m*c_s_m*rho_m;
84 T lam_m = rho_m*c_p_m*c_p_m-2*mu_m;
85
86
87 T Tx_m,Ty_m,Tz_m,Tx_p,Ty_p,Tz_p;
88 T vx_m,vy_m,vz_m,vx_p,vy_p,vz_p;
89
90 extract_tractions_and_particle_velocity(n_p,QR+idx_QLR(i,j,0),Tx_p,Ty_p,Tz_p,vx_p,vy_p,vz_p );
91 extract_tractions_and_particle_velocity(n_m,QL+idx_QLR(i,j,0),Tx_m,Ty_m,Tz_m,vx_m,vy_m,vz_m );
92
93 localBasis(n_p, m_p, l_p, 3);
94 localBasis(n_m, m_m, l_m, 3);
95
96 T Tn_m,Tm_m,Tl_m,vn_m,vm_m,vl_m;
97 T Tn_p,Tm_p,Tl_p,vn_p,vm_p,vl_p;
98
99 // rotate fields into l, m, n basis
100 rotate_into_orthogonal_basis(n_m,m_m,l_m,Tx_m,Ty_m,Tz_m,Tn_m,Tm_m,Tl_m);
101 rotate_into_orthogonal_basis(n_m,m_m,l_m,vx_m,vy_m,vz_m,vn_m,vm_m,vl_m);
102 rotate_into_orthogonal_basis(n_p,m_p,l_p,Tx_p,Ty_p,Tz_p,Tn_p,Tm_p,Tl_p);
103 rotate_into_orthogonal_basis(n_p,m_p,l_p,vx_p,vy_p,vz_p,vn_p,vm_p,vl_p);
104
105
106 // extract local s-wave and p-wave impedances
107 T zs_p=rho_p*c_s_p;
108 T zp_p=rho_p*c_p_p;
109 T zs_m=rho_m*c_s_m;
110 T zp_m=rho_m*c_p_m;
111
112 // impedance must be greater than zero !
113 if (zp_p <= 0.0 || zp_m <= 0.0){
114 std::cout << "At position: " << isBoundaryFace << "\n";
115 std::cout<<zs_p<<" "<<zp_p<<" "<<zs_m<<" "<<zp_m<<"\n";
116 std::cout<<" Impedance must be greater than zero ! "<< std::endl;
117 std::exit(-1);
118 }
119
120 // generate interface data preserving the amplitude of the outgoing charactertritics
121 // and satisfying interface conditions exactly.
122 T vn_hat_p,vm_hat_p,vl_hat_p,Tn_hat_p,Tm_hat_p,Tl_hat_p;
123 T vn_hat_m,vm_hat_m,vl_hat_m,Tn_hat_m,Tm_hat_m,Tl_hat_m;
124
125 if (isBoundaryFace) {
126 // 1 absorbing 0 free surface
127 T r= faceIndex==2 ? 1 : 0;
128 // T r=1.;
129 riemannSolver_boundary(faceIndex,r,vn_m,vm_m,vl_m,Tn_m,Tm_m,Tl_m,zp_m,zs_m,vn_hat_m,vm_hat_m,vl_hat_m,Tn_hat_m,Tm_hat_m,Tl_hat_m);
130 riemannSolver_boundary(faceIndex,r,vn_p,vm_p,vl_p,Tn_p,Tm_p,Tl_p,zp_p,zs_p,vn_hat_p,vm_hat_p,vl_hat_p,Tn_hat_p,Tm_hat_p,Tl_hat_p);
131 }else {
132 riemannSolver_Nodal(vn_p,vn_m, Tn_p, Tn_m, zp_p , zp_m, vn_hat_p , vn_hat_m, Tn_hat_p, Tn_hat_m);
133 riemannSolver_Nodal(vm_p,vm_m, Tm_p, Tm_m, zs_p , zs_m, vm_hat_p , vm_hat_m, Tm_hat_p, Tm_hat_m);
134 riemannSolver_Nodal(vl_p,vl_m, Tl_p, Tl_m, zs_p , zs_m, vl_hat_p , vl_hat_m, Tl_hat_p, Tl_hat_m);
135 }
136
137 //generate fluctuations in the local basis coordinates: n, m, l
138 generate_fluctuations_left(zp_m,Tn_m,Tn_hat_m,vn_m,vn_hat_m,FLn);
139 generate_fluctuations_left(zs_m,Tm_m,Tm_hat_m,vm_m,vm_hat_m,FLm);
140 generate_fluctuations_left(zs_m,Tl_m,Tl_hat_m,vl_m,vl_hat_m,FLl);
141
142 generate_fluctuations_right(zp_p,Tn_p,Tn_hat_p,vn_p,vn_hat_p,FRn);
143 generate_fluctuations_right(zs_p,Tm_p,Tm_hat_p,vm_p,vm_hat_p,FRm);
144 generate_fluctuations_right(zs_p,Tl_p,Tl_hat_p,vl_p,vl_hat_p,FRl);
145
146 FL_n = FLn/zp_m;
147 if(zs_m > 0){
148 FL_m = FLm/zs_m;
149 FL_l = FLl/zs_m;
150 }else{
151 FL_m=0;
152 FL_l=0;
153 }
154
155 FR_n = FRn/zp_p;
156 if(zs_p > 0){
157 FR_m = FRm/zs_p;
158 FR_l = FRl/zs_p;
159 }else{
160 FR_m=0;
161 FR_l=0;
162 }
163
164 // rotate back to the physical coordinates x, y, z
165 rotate_into_physical_basis(n_m,m_m,l_m,FLn,FLm,FLl,FLx,FLy,FLz);
166 rotate_into_physical_basis(n_p,m_p,l_p,FRn,FRm,FRl,FRx,FRy,FRz);
167 rotate_into_physical_basis(n_m,m_m,l_m,FL_n,FL_m,FL_l,FL_x,FL_y,FL_z);
168 rotate_into_physical_basis(n_p,m_p,l_p,FR_n,FR_m,FR_l,FR_x,FR_y,FR_z);
169
170 // construct flux fluctuation vectors obeying the eigen structure of the PDE
171 // and choose physically motivated penalties such that we can prove
172 // numerical stability.
173
174 FR[idx_FLR(i,j, 0)] = norm_p_qr/rho_p*FRx;
175 FL[idx_FLR(i,j, 0)] = norm_m_qr/rho_m*FLx;
176
177 FR[idx_FLR(i,j, 1)] = norm_p_qr/rho_p*FRy;
178 FL[idx_FLR(i,j, 1)] = norm_m_qr/rho_m*FLy;
179
180 FR[idx_FLR(i,j, 2)] = norm_p_qr/rho_p*FRz;
181 FL[idx_FLR(i,j, 2)] = norm_m_qr/rho_m*FLz;
182
183 FL[idx_FLR(i,j, 3)] = norm_m_qr*((2*mu_m+lam_m)*n_m[0]*FL_x+lam_m*n_m[1]*FL_y+lam_m*n_m[2]*FL_z);
184 FL[idx_FLR(i,j, 4)] = norm_m_qr*((2*mu_m+lam_m)*n_m[1]*FL_y+lam_m*n_m[0]*FL_x+lam_m*n_m[2]*FL_z);
185 FL[idx_FLR(i,j, 5)] = norm_m_qr*((2*mu_m+lam_m)*n_m[2]*FL_z+lam_m*n_m[0]*FL_x+lam_m*n_m[1]*FL_y);
186
187 FR[idx_FLR(i,j, 3)] = -norm_p_qr*((2*mu_p+lam_p)*n_p[0]*FR_x+lam_p*n_p[1]*FR_y+lam_p*n_p[2]*FR_z);
188 FR[idx_FLR(i,j, 4)] = -norm_p_qr*((2*mu_p+lam_p)*n_p[1]*FR_y+lam_p*n_p[0]*FR_x+lam_p*n_p[2]*FR_z);
189 FR[idx_FLR(i,j, 5)] = -norm_p_qr*((2*mu_p+lam_p)*n_p[2]*FR_z+lam_p*n_p[0]*FR_x+lam_p*n_p[1]*FR_y);
190
191 FL[idx_FLR(i,j, 6)] = norm_m_qr*mu_m*(n_m[1]*FL_x + n_m[0]*FL_y);
192 FL[idx_FLR(i,j, 7)] = norm_m_qr*mu_m*(n_m[2]*FL_x + n_m[0]*FL_z);
193 FL[idx_FLR(i,j, 8)] = norm_m_qr*mu_m*(n_m[2]*FL_y + n_m[1]*FL_z);
194
195 FR[idx_FLR(i,j, 6)] = -norm_p_qr*mu_p*(n_p[1]*FR_x + n_p[0]*FR_y);
196 FR[idx_FLR(i,j, 7)] = -norm_p_qr*mu_p*(n_p[2]*FR_x + n_p[0]*FR_z);
197 FR[idx_FLR(i,j, 8)] = -norm_p_qr*mu_p*(n_p[2]*FR_y + n_p[1]*FR_z);
198 }
199 }
200
201#ifdef OPT_KERNELS
202 converter::FFace_generic2optimised(FL,FL_);
203 converter::FFace_generic2optimised(FR,FR_);
204#endif
205
206}
207
208template void loh::riemannSolver::riemannSolver<float>(
209 float* const FL_,
210 float* const FR_,
211 const float* const QL_,
212 const float* const QR_,
213 const double t,
214 const double dt,
215 const tarch::la::Vector<Dimensions, double>& lengthScale,
216 const int normalNonZeroIndex,
217 bool isBoundaryFace,
218 int faceIndex);
219template void loh::riemannSolver::riemannSolver<double>(
220 double* const FL_,
221 double* const FR_,
222 const double* const QL_,
223 const double* const QR_,
224 const double t,
225 const double dt,
226 const tarch::la::Vector<Dimensions, double>& lengthScale,
227 const int normalNonZeroIndex,
228 bool isBoundaryFace,
229 int faceIndex);
230
231//Gram Schmidt orthonormalization
232template<typename T>
233void loh::riemannSolver::Gram_Schmidt(T* const y, T* const z){
234 T a_yz = y[0]*z[0] + y[1]*z[1] + y[2]*z[2];
235
236 for (int i = 0; i< 3; i++){
237 z[i] = z[i] - a_yz*y[i];
238 }
239
240 T norm_z = std::sqrt(z[0]*z[0] + z[1]*z[1] + z[2]*z[2]);
241
242 for (int i = 0; i< 3; i++){
243 z[i] = z[i]/norm_z;
244 }
245}
246
247template<typename T>
248void loh::riemannSolver::localBasis(T* const n, T* const m, T* const l, int d){
249 if (d == 2){
250 l[0] = 0.;
251 l[1] = 0.;
252 l[2] = 1.0;
253
254 m[0] = n[1]*l[2]-n[2]*l[1];
255 m[1] = -(n[0]*l[2]-n[2]*l[0]);
256 m[2] = n[0]*l[1]-n[1]*l[0];
257 }else if (d == 3){
258 T tol, diff_norm1, diff_norm2;
259 tol = 1e-12;
260 m[0] = 0.;
261 m[1] = 1.;
262 m[2] = 0.;
263
264 diff_norm1 = kernels::pow2(n[0]-m[0]) + kernels::pow2(n[1]-m[1]) + kernels::pow2(n[2]-m[2]);
265 diff_norm2 = kernels::pow2(n[0]+m[0]) + kernels::pow2(n[1]+m[1]) + kernels::pow2(n[2]+m[2]);
266
267 if (diff_norm1 >= kernels::pow2(tol) && diff_norm2 >= kernels::pow2(tol)){
268 Gram_Schmidt(n, m);
269 }else{
270 m[0] = 0.;
271 m[1] = 0.;
272 m[2] = 1.;
273 Gram_Schmidt(n, m);
274 }
275 l[0] = n[1]*m[2]-n[2]*m[1];
276 l[1] = -(n[0]*m[2]-n[2]*m[0]);
277 l[2] = n[0]*m[1]-n[1]*m[0];
278 }
279}
280
281
282
283template<typename T>
284void loh::riemannSolver::riemannSolver_Nodal(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){
285 T p=0;
286 T q=0;
287 T phi=0;
288 T v_hat=0;
289 T eta=0;
290
291 p=z_m*v_p + sigma_p;
292 q=z_p*v_m - sigma_m;
293
294 if(z_p > 0 && z_m > 0){
295 eta=(z_p*z_m)/(z_p+z_m);
296
297 phi= eta*(p/z_p - q/z_m);
298
299 sigma_hat_p=phi;
300 sigma_hat_m=phi;
301
302 v_hat_p=(q+phi)/z_m;
303 v_hat_m=(p-phi)/z_p;
304 }else if(z_p > 0){
305 sigma_hat_p=0;
306 sigma_hat_m=sigma_m;
307
308 v_hat_p=v_p;
309 v_hat_m=v_m;
310 }else if(z_m > 0){
311 sigma_hat_p=sigma_p;
312 sigma_hat_m=0;
313
314 v_hat_p=v_p;
315 v_hat_m=v_m;
316 }else{
317 sigma_hat_p=sigma_p;
318 sigma_hat_m=sigma_m;
319
320 v_hat_p=v_p;
321 v_hat_m=v_m;
322 }
323 }
324
325template<typename T>
326void loh::riemannSolver::riemannSolver_BC0(T v, T sigma, T z, T r, T& v_hat, T& sigma_hat){
327 T p = 0.5*(z*v + sigma);
328 if(z > 0){
329 v_hat = (1+r)/z*p;
330 sigma_hat = (1-r)*p;
331 }else{
332 v_hat = v;
333 sigma_hat = sigma;
334 }
335}
336
337template<typename T>
338void loh::riemannSolver::riemannSolver_BCn(T v,T sigma, T z, T r, T& v_hat, T& sigma_hat){
339 T q = 0.5*(z*v - sigma);
340 if(z > 0){
341 v_hat = (1+r)/z*q;
342 sigma_hat = -(1-r)*q;
343 }else{
344 v_hat = v;
345 sigma_hat = sigma;
346 }
347}
348
349
350
351template<typename T>
352void loh::riemannSolver::extract_tractions_and_particle_velocity(T* const n,const T* const Q, T& Tx,T& Ty,T& Tz,T& vx,T& vy,T& vz ){
353 T sigma_xx = Q[3];
354 T sigma_yy = Q[4];
355 T sigma_zz = Q[5];
356 T sigma_xy = Q[6];
357 T sigma_xz = Q[7];
358 T sigma_yz = Q[8];
359
360 Tx = n[0]*sigma_xx + n[1]*sigma_xy + n[2]*sigma_xz;
361 Ty = n[0]*sigma_xy + n[1]*sigma_yy + n[2]*sigma_yz;
362 Tz = n[0]*sigma_xz + n[1]*sigma_yz + n[2]*sigma_zz;
363
364 vx = Q[0];
365 vy = Q[1];
366 vz = Q[2];
367}
368
369template<typename T>
370void loh::riemannSolver::rotate_into_orthogonal_basis(T* const n,T* const m,T* const l, T Tx,T Ty,T Tz, T& Tn, T& Tm, T& Tl){
371 Tn= Tx*n[0] + Ty*n[1] + Tz*n[2];
372 Tm= Tx*m[0] + Ty*m[1] + Tz*m[2];
373 Tl= Tx*l[0] + Ty*l[1] + Tz*l[2];
374}
375
376template<typename T>
377void loh::riemannSolver::rotate_into_physical_basis(T* const n,T* const m,T* const l, T Fn,T Fm,T Fl, T& Fx, T& Fy, T& Fz){
378 Fx = n[0]*Fn + m[0]*Fm + l[0]*Fl;
379 Fy = n[1]*Fn + m[1]*Fm + l[1]*Fl;
380 Fz = n[2]*Fn + m[2]*Fm + l[2]*Fl;
381}
382
383template<typename T>
384void loh::riemannSolver::generate_fluctuations_left(T z, T myT,T T_hat,T v, T v_hat, T& F){
385 F = 0.5*(z*(v-v_hat) + (myT-T_hat));
386}
387
388template<typename T>
389void loh::riemannSolver::generate_fluctuations_right(T z, T myT,T T_hat,T v, T v_hat, T& F){
390 F = 0.5*(z*(v-v_hat) - (myT-T_hat));
391}
392
393template<typename T>
394void loh::riemannSolver::riemannSolver_boundary(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)
395{
396 if (faceIndex % 2 == 0) {
397 riemannSolver_BC0(vn, Tn, zp, r, vn_hat, Tn_hat);
398 riemannSolver_BC0(vm, Tm, zs, r, vm_hat, Tm_hat);
399 riemannSolver_BC0(vl, Tl, zs, r, vl_hat, Tl_hat);
400 }
401
402 if (faceIndex % 2 == 1) {
403 riemannSolver_BCn(vn, Tn, zp, r, vn_hat, Tn_hat);
404 riemannSolver_BCn(vm, Tm, zs, r, vm_hat, Tm_hat);
405 riemannSolver_BCn(vl, Tl, zs, r, vl_hat, Tl_hat);
406 }
407}
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ j
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ k
double pow2(const double &a)
Definition KernelUtils.h:30
void riemannSolver(T *const FL, T *const FR, const T *const QL, const T *const QR, const double t, const double dt, const tarch::la::Vector< Dimensions, double > &h, const int direction, bool isBoundaryFace, int faceIndex)
void rotate_into_orthogonal_basis(T *const n, T *const m, T *const l, T Tx, T Ty, T Tz, T &Tn, T &Tm, T &Tl)
void extract_tractions_and_particle_velocity(T *const n, const T *const Q, T &Tx, T &Ty, T &Tz, T &vx, T &vy, T &vz)
void riemannSolver_Nodal(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)
void generate_fluctuations_left(T z, T myT, T T_hat, T v, T v_hat, T &F)
void rotate_into_physical_basis(T *const n, T *const m, T *const l, T Fn, T Fm, T Fl, T &Fx, T &Fy, T &Fz)
void Gram_Schmidt(T *const y, T *const z)
void riemannSolver_boundary(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)
void riemannSolver_BC0(T v, T sigma, T z, T r, T &v_hat, T &sigma_hat)
void riemannSolver_BCn(T v, T sigma, T z, T r, T &v_hat, T &sigma_hat)
void generate_fluctuations_right(T z, T myT, T T_hat, T v, T v_hat, T &F)
void localBasis(T *const n, T *const m, T *const l, int d)
examples::exahype2::elastic::VariableShortcuts sh
Simple vector class.
Definition Vector.h:134