Peano 4
Loading...
Searching...
No Matches
LOH1.cpp
Go to the documentation of this file.
1#include "LOH1.h"
2
3#include <algorithm>
4#include <iomanip>
5
6tarch::logging::Log examples::exahype2::loh1::LOH1::_log( "examples::exahype2::loh1::LOH1" );
7
8
9// TODO: No point sources supported
10//void examples::exahype2::loh1::LOH1::prescribeLOH1InitialData(
11// const tarch::la::Vector<Dimensions,double>& x,
12// double Q[]
13//) {
14// // initial conditions
15// constexpr idx_rho = 9;
16// if( x[1] < 1.0) {
17// Q[s.rho+0] = 2.6;
18// Q[s.cp] = 4.0;
19// Q[s.cs] = 2.0;
20// }else{
21// Q[s.rho+0] = 2.7;
22// Q[s.cp] = 6.0;
23// Q[s.cs] = 3.464;
24// }
25//}
26
29 double Q[]
30) {
31 double center_curve[3];
32
33 center_curve[0] = 0.0;
34 center_curve[1] = 0.0;
35 center_curve[2] = 2.0;
36
37 // 1.0 is the parameter from the setup, but we need at least one voxel
38 bool layerWidth = std::max( getMinMeshSize() / _NumberOfFiniteVolumesPerAxisPerPatch, 1.0 );
39
40 bool upperLayer = x(2) <= layerWidth;
41
42 //double rho = upperLayer ? 2.67 : 2.6;
43 //double cp = upperLayer ? 6.0 : 4.0;
44 //double cs = upperLayer ? 3.343 : 2.0;
45 double rho = upperLayer ? 4.0 : 6.0;
46 double cp = upperLayer ? 2.0 : 3.464;
47 double cs = upperLayer ? 2.6 : 2.7;
48
49 Q[s.rho ] = rho;
50 Q[s.cp ] = cp;
51 Q[s.cs ] = cs;
52 Q[s.alpha] = 1.0;
53
54 double radius = 2.0*layerWidth;
55 //double height = 3.0;
56
57 Q[ s.v + 0 ] = std::exp(-((x[0]-center_curve[0])*(x[0]-center_curve[0])+
58 (x[1]-center_curve[1])*(x[1]-center_curve[1])+
59 (x[2]-center_curve[2])*(x[2]-center_curve[2]))/radius) * radius;
60 Q[ s.v + 1 ] = Q[ s.v + 0 ];
61 Q[ s.v + 2 ] = Q[ s.v + 1 ];
62
63
64 for(int i= 0 ; i < 6 ; i++){
65 Q[ s.sigma + i ] = 0.0;
66 }
67}
68
69
70
72 double * __restrict__ Q, // [13]
75 double t,
76 double dt
77) {
78 logTraceInWith3Arguments( "adjustSolution(...)", x, h, t );
79
80 if (tarch::la::equals(t,0.0) ) {
81 prescribeGaussianWave(x,Q);
82 }
83
84 logTraceOut( "adjustSolution(...)" );
85}
86
87
89 const double* __restrict__ Q, // [13]
92 double t,
93 int normal
94) {
95 double cp = Q[s.cp];
96 double cs = Q[s.cs];
97
98/*
99 lambda[0] = -cp;
100 lambda[1] = cp;
101 lambda[2] = -cs;
102 lambda[3] = cp; I think this is a bug!
103 lambda[4] = 0.0;
104 lambda[5] = 0.0;
105 lambda[6] = 0.0;
106 lambda[7] = 0.0;
107 lambda[8] = 0.0;
108 lambda[9] = 0.0;
109*/
110
111 double result = 0.0;
112 result = std::max( result, -cp );
113 result = std::max( result, cp );
114 result = std::max( result, -cs );
115 result = std::max( result, cs );
116 return result;
117}
118
119
121 const double * __restrict__ Qinside,
122 double * __restrict__ Qoutside,
123 const tarch::la::Vector<Dimensions,double>& faceCentre,
125 double t,
126 int normal
127) {
128 logTraceInWith4Arguments( "boundaryConditions(...)", faceCentre, volumeH, t, normal );
129 std::fill_n(Qoutside, 13, 0);
130
131 #ifdef Asserts
132 if(!std::isfinite(1.0/Qinside[s.rho])){
133 std::cout << "FV: rho_inv not finite in boundary values" << std::endl;
134 for(int i=0;i< Unknowns+NumberOfParameters; i++){
135 std::cout << Qinside[i] << std::endl;
136 }
137 }
138 #endif
139
140 assertion2(std::isfinite(1.0/Qinside[s.rho]),
141 "rho_i is not finite",Qinside[s.rho]);
142
143 Qoutside[ s.alpha ] = Qinside[ s.alpha ];
144 Qoutside[ s.rho ] = Qinside[ s.rho ];
145 Qoutside[ s.cp ] = Qinside[ s.cp ];
146 Qoutside[ s.cs ] = Qinside[ s.cs ];
147
148 assertion2(std::isfinite(1.0/Qoutside[s.rho]), "rho_i is not finite",Qoutside[s.rho]);
149
150 logTraceOut( "boundaryConditions(...)" );
151}
152
153
155 const double * __restrict__ Q, // [9+4],
156 const double * __restrict__ deltaQ, // [9+4],
157 const tarch::la::Vector<Dimensions,double>& faceCentre,
159 double t,
160 int normal,
161 double * __restrict__ BgradQ // BgradQ[13] -> noe, nur 9
162) {
163 logTraceIn( "nonconservativeProduct(...)" );
164
165 double rho = Q[s.rho];
166 double cp = Q[s.cp];
167 double cs = Q[s.cs];
168 double alpha = Q[s.alpha];
169
170 double mu = rho*cs*cs;
171 double lambda = rho*cp*cp-2*mu;
172 double rho_i = 1.0/rho;
173
174 assertion2(std::isfinite(rho_i),"BGradQ: rho inverse not finite",rho);
175
176 double alpha_i = alpha/(alpha*alpha + 1.0e-13 * (1.0-alpha));
177
178 double u= Q[s.v + 0] * alpha_i;
179 double v= Q[s.v + 1] * alpha_i;
180 double w= Q[s.v + 2] * alpha_i;
181
182 assertion(normal>=0);
183 assertion(normal<3);
184 switch ( normal ) {
185 //dx
186 case 0:
187 {
188 double alpha_x = deltaQ[s.alpha];
189 BgradQ[s.v + 0] = -rho_i*(deltaQ[s.sigma + 0] + 2*Q[s.sigma + 0]*alpha_x); // AQx(7) = - irho * Qx(1) - 2*Q(1)*irho*Qx(13)
190 BgradQ[s.v + 1] = -rho_i*(deltaQ[s.sigma + 3] + 2*Q[s.sigma + 3]*alpha_x); // AQx(8) = - irho * Qx(4) - 2*Q(4)*irho*Qx(13)
191 BgradQ[s.v + 2] = -rho_i*(deltaQ[s.sigma + 4] + 2*Q[s.sigma + 4]*alpha_x); // AQx(9) = - irho * Qx(6) - 2*Q(6)*irho*Qx(13)
192 BgradQ[s.sigma + 0] =-(lambda+2*mu) * ( deltaQ[s.v + 0] - u * alpha_x); // AQx(1) = - (lam+2*mu)*Qx(7) + (lam+2*mu)*u(1)*Qx(13)
193 BgradQ[s.sigma + 1] = -lambda * ( deltaQ[s.v + 0] - u * alpha_x); // AQx(2) = - lam*Qx(7) + lam*u(1)*Qx(13)
194 BgradQ[s.sigma + 2] = -lambda * ( deltaQ[s.v + 0] - u * alpha_x); // AQx(3) = - lam*Qx(7) + lam*u(1)*Qx(13)
195 BgradQ[s.sigma + 3] = -mu * ( deltaQ[s.v + 1] - v * alpha_x); // AQx(4) = - mu *Qx(8) + mu *u(2)*Qx(13)
196 BgradQ[s.sigma + 4] = -mu * ( deltaQ[s.v + 2] - w * alpha_x); // AQx(6) = - mu *Qx(9) + mu *u(3)*Qx(13)
197 BgradQ[s.sigma + 5] = 0.0; // AQx(5) = 0.0
198 }
199 break;
200
201 //dy
202 case 1:
203 {
204 double alpha_y = deltaQ[s.alpha];
205 BgradQ[s.v + 0] = -rho_i*(deltaQ[s.sigma + 3] + 2*Q[s.sigma + 3]*alpha_y); // BQy(7) = - irho * Qy(4) - 2*Q(4)*irho*Qy(13)
206 BgradQ[s.v + 1] = -rho_i*(deltaQ[s.sigma + 1] + 2*Q[s.sigma + 1]*alpha_y); // BQy(8) = - irho * Qy(2) - 2*Q(2)*irho*Qy(13)
207 BgradQ[s.v + 2] = -rho_i*(deltaQ[s.sigma + 5] + 2*Q[s.sigma + 5]*alpha_y); // BQy(9) = - irho * Qy(5) - 2*Q(5)*irho*Qy(13)
208 BgradQ[s.sigma + 0] = -lambda * ( deltaQ[s.v + 1] - v * alpha_y); // BQy(1) = - lam*Qy(8) + lam*u(2)*Qy(13)
209 BgradQ[s.sigma + 1] = -(lambda+2*mu) * ( deltaQ[s.v + 1] - v * alpha_y); // BQy(2) = - (lam+2*mu)*Qy(8) + (lam+2*mu)*u(2)*Qy(13)
210 BgradQ[s.sigma + 2] = -lambda * ( deltaQ[s.v + 1] - v * alpha_y); // BQy(3) = - lam*Qy(8) + lam*u(2)*Qy(13)
211 BgradQ[s.sigma + 3] = -mu * ( deltaQ[s.v + 0] - u * alpha_y); // BQy(4) = - mu *Qy(7) + mu *u(1)*Qy(13)
212 BgradQ[s.sigma + 4] = 0.0; // BQy(6) = 0.0
213 BgradQ[s.sigma + 5] = -mu * ( deltaQ[s.v + 2] - w * alpha_y); // BQy(5) = - mu *Qy(9) + mu *u(3)*Qy(13)
214 }
215 break;
216 // !
217 //dz
218 case 2:
219 {
220 double alpha_z = deltaQ[s.alpha];
221 BgradQ[s.v + 0] = -rho_i*(deltaQ[s.sigma + 4] + 2*Q[s.sigma + 4]*alpha_z); // CQz(7) = - irho * Qz(6) - 2*Q(6)*irho*Qz(13)
222 BgradQ[s.v + 1] = -rho_i*(deltaQ[s.sigma + 5] + 2*Q[s.sigma + 5]*alpha_z); // CQz(8) = - irho * Qz(5) - 2*Q(5)*irho*Qz(13)
223 BgradQ[s.v + 2] = -rho_i*(deltaQ[s.sigma + 2] + 2*Q[s.sigma + 2]*alpha_z); // CQz(9) = - irho * Qz(3) - 2*Q(3)*irho*Qz(13)
224 BgradQ[s.sigma + 0] = -lambda * ( deltaQ[s.v + 2] - w * alpha_z); // CQz(1) = - lam*Qz(9) + lam*u(3)*Qz(13)
225 BgradQ[s.sigma + 1] = -lambda * ( deltaQ[s.v + 2] - w * alpha_z); // CQz(2) = - lam*Qz(9) + lam*u(3)*Qz(13)
226 BgradQ[s.sigma + 2] = -(lambda+2*mu) * ( deltaQ[s.v + 2] - w * alpha_z); // CQz(3) = - (lam+2*mu)*Qz(9) + (lam+2*mu)*u(3)*Qz(13)
227 BgradQ[s.sigma + 3] = 0.0; // CQz(4) = 0.0
228 BgradQ[s.sigma + 4] = -mu * ( deltaQ[s.v + 0] - u * alpha_z); // CQz(6) = - mu *Qz(7) + mu *u(1)*Qz(13)
229 BgradQ[s.sigma + 5] = -mu * ( deltaQ[s.v + 1] - v * alpha_z); // CQz(5) = - mu *Qz(8) + mu *u(2)*Qz(13)
230 }
231 break;
232
233 default:
234#ifdef Asserts
235 assertionMsg(false,"Normal direction/gradQ component must be in {0,1,2}.");
236#endif
237 break;
238 }
239
240#ifdef Asserts
241 assertion2(std::isfinite(BgradQ[s.v + 0]),"BGradQ: v not finite",BgradQ[s.v + 0]);
242 assertion2(std::isfinite(BgradQ[s.v + 1]),"BGradQ: u not finite",BgradQ[s.v + 1]);
243 assertion2(std::isfinite(BgradQ[s.v + 2]),"BGradQ: w not finite",BgradQ[s.v + 2]);
244 assertion10(std::isfinite(BgradQ[s.sigma + 0]),"BGradQ: sigma_xx not finite",BgradQ[s.sigma + 0],lambda,alpha_i,gradQ[idx(1,s.v + 1)],gradQ[idx(2,s.v + 2)],v,w,alpha_y,alpha_z);
245 assertion2(std::isfinite(BgradQ[s.sigma + 1]),"BGradQ: sigma_yy not finite",BgradQ[s.sigma + 1]);
246 assertion2(std::isfinite(BgradQ[s.sigma + 2]),"BGradQ: sigma_zz not finite",BgradQ[s.sigma + 2]);
247 assertion2(std::isfinite(BgradQ[s.sigma + 3]),"BGradQ: sigma_xy not finite",BgradQ[s.sigma + 3]);
248 assertion2(std::isfinite(BgradQ[s.sigma + 4]),"BGradQ: sigma_xz not finite",BgradQ[s.sigma + 4]);
249 assertion2(std::isfinite(BgradQ[s.sigma + 5]),"BGradQ: sigma_yz not finite",BgradQ[s.sigma + 5]);
250#endif
251 logTraceOut( "nonconservativeProduct(...)" );
252}
#define assertion10(expr, param0, param1, param2, param3, param4, param5, param6, param7, param8, param9)
#define assertion2(expr, param0, param1)
#define assertionMsg(expr, message)
#define assertion(expr)
#define logTraceOut(methodName)
Definition Log.h:379
#define logTraceInWith4Arguments(methodName, argument0, argument1, argument2, argument3)
Definition Log.h:373
#define logTraceInWith3Arguments(methodName, argument0, argument1, argument2)
Definition Log.h:372
#define logTraceIn(methodName)
Definition Log.h:369
static tarch::logging::Log _log
Definition LOH1.h:26
virtual void boundaryConditions(double Qinside[9+4], double Qoutside[9+4], const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, int normal) override
Apply boundary conditions.
virtual double maxEigenvalue(double Q[9+4], const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, int normal) override
Determine max eigenvalue over Jacobian in a given point with solution values (states) Q.
struct examples::exahype2::loh1::LOH1::VariablesShortcuts s
void adjustSolution(double *__restrict__ Q, double Q[9+4], const tarch::la::Vector< Dimensions, double > &volumeCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t) override
Feel free to change the solution in the particular finite volume.
void prescribeGaussianWave(const tarch::la::Vector< Dimensions, double > &x, double Q[])
Definition LOH1.cpp:27
void nonconservativeProduct(double Q[9+4], const double gradQ[9+4][Dimensions], const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, int normal, double BgradQ[9]) override
Log Device.
Definition Log.h:516
examples::exahype2::elastic::VariableShortcuts s
Definition loh.cpp:10
bool equals(const Matrix< Rows, Cols, Scalar > &lhs, const Matrix< Rows, Cols, Scalar > &rhs, const Scalar &tolerance=NUMERICAL_ZERO_DIFFERENCE)
Compares to matrices on equality by means of a numerical accuracy.
Simple vector class.
Definition Vector.h:134