Peano 4
Loading...
Searching...
No Matches
ADERDGLOH1.cpp
Go to the documentation of this file.
1#include "ADERDGLOH1.h"
3
4
5tarch::logging::Log examples::exahype2::loh1::ADERDGLOH1::_log( "examples::exahype2::loh1::ADERDGLOH1" );
6
9 double Q[]
10) {
11 double center_curve[3];
12
13 center_curve[0] = 0.0;
14 center_curve[1] = 0.0;
15 center_curve[2] = 2.0;
16
17 // 1.0 is the parameter from the setup, but we need at least one voxel
18 bool layerWidth = std::max( getMinMeshSize() / (2*Order+1), 1.0 ); // TODO guess
19
20 bool upperLayer = x(2) <= layerWidth;
21
22 //double rho = upperLayer ? 2.67 : 2.6;
23 //double cp = upperLayer ? 6.0 : 4.0;
24 //double cs = upperLayer ? 3.343 : 2.0;
25 double rho = upperLayer ? 4.0 : 6.0;
26 double cp = upperLayer ? 2.0 : 3.464;
27 double cs = upperLayer ? 2.6 : 2.7;
28
29 Q[s.rho ] = rho;
30 Q[s.cp ] = cp;
31 Q[s.cs ] = cs;
32 Q[s.alpha] = 1.0;
33
34 double radius = 2.0*layerWidth;
35 //double height = 3.0;
36
37 Q[ s.v + 0 ] = std::exp(-((x[0]-center_curve[0])*(x[0]-center_curve[0])+
38 (x[1]-center_curve[1])*(x[1]-center_curve[1])+
39 (x[2]-center_curve[2])*(x[2]-center_curve[2]))/radius) * radius;
40 Q[ s.v + 1 ] = Q[ s.v + 0 ];
41 Q[ s.v + 2 ] = Q[ s.v + 1 ];
42
43
44 for(int i= 0 ; i < 6 ; i++){
45 Q[ s.sigma + i ] = 0.0;
46 }
47}
48
50 double * __restrict__ Q,
52 double t
53) {
54 logTraceInWith2Arguments( "adjustSolution(...)", x, t );
55 if (tarch::la::equals(t,0.0) ) {
56 prescribeGaussianWave(x,Q);
57 }
58 else {
59 // other stuff
60 }
61 logTraceOut( "adjustSolution(...)" );
62}
63
64
65
66
68 const double * __restrict__ Qinside, // Qinside[9+4]
69 double * __restrict__ Qoutside, // Qoutside[9+4]
71 double t,
72 int normal
73) {
74 logTraceInWith3Arguments( "boundaryConditions(...)", x, t, normal );
75 #ifdef Asserts
76 if(!std::isfinite(1.0/Qinside[s.rho])){
77 std::cout << "FV: rho_inv not finite in boundary values" << std::endl;
78 for(int i=0;i< Unknowns+NumberOfParameters; i++){
79 std::cout << Qinside[i] << std::endl;
80 }
81 }
82 #endif
83
84 assertion2(std::isfinite(1.0/Qinside[s.rho]),
85 "rho_i is not finite",Qinside[s.rho]);
86
87 std::fill_n(Qoutside, 13, 0);
88
89 Qoutside[ s.alpha ] = Qinside[ s.alpha ];
90 Qoutside[ s.rho ] = Qinside[ s.rho ];
91 Qoutside[ s.cp ] = Qinside[ s.cp ];
92 Qoutside[ s.cs ] = Qinside[ s.cs ];
93
94 assertion2(std::isfinite(1.0/Qoutside[s.rho]), "rho_i is not finite",Qoutside[s.rho]);
95
96 logTraceOut( "boundaryConditions(...)" );
97}
98
100 const double * __restrict__ Q, // Q[9+4],
102 double t,
103 int normal
104) {
105 logTraceInWith2Arguments( "maxEigenvalue(...)", x, t );
106 double cp = Q[s.cp];
107 double cs = Q[s.cs];
108
109 double result = 0.0;
110 result = std::max( result, std::abs(cp) );
111 result = std::max( result, std::abs(cs) );
112 logTraceOut( "maxEigenvalue(...)" );
113 return result;
114}
115
116
118 const double * __restrict__ Q, // [9+4],
119 const double * __restrict__ deltaQ, // [9+4],
121 double t,
122 int normal,
123 double * __restrict__ BgradQ // BgradQ[13]
124 ) {
125 logTraceIn( "nonconservativeProduct(...)" );
126
127 double rho = Q[s.rho];
128 double cp = Q[s.cp];
129 double cs = Q[s.cs];
130 double alpha = Q[s.alpha];
131
132 double mu = rho*cs*cs;
133 double lambda = rho*cp*cp-2*mu;
134 double rho_i = 1.0/rho;
135
136 assertion2(std::isfinite(rho_i),"BGradQ: rho inverse not finite",rho);
137
138 double alpha_i = alpha/(alpha*alpha + 1.0e-13 * (1.0-alpha));
139
140 double u= Q[s.v + 0] * alpha_i;
141 double v= Q[s.v + 1] * alpha_i;
142 double w= Q[s.v + 2] * alpha_i;
143
144 assertion(normal>=0);
145 assertion(normal<3);
146 switch ( normal ) {
147 //dx
148 case 0:
149 {
150 double alpha_x = deltaQ[s.alpha];
151 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)
152 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)
153 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)
154 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)
155 BgradQ[s.sigma + 1] = -lambda * ( deltaQ[s.v + 0] - u * alpha_x); // AQx(2) = - lam*Qx(7) + lam*u(1)*Qx(13)
156 BgradQ[s.sigma + 2] = -lambda * ( deltaQ[s.v + 0] - u * alpha_x); // AQx(3) = - lam*Qx(7) + lam*u(1)*Qx(13)
157 BgradQ[s.sigma + 3] = -mu * ( deltaQ[s.v + 1] - v * alpha_x); // AQx(4) = - mu *Qx(8) + mu *u(2)*Qx(13)
158 BgradQ[s.sigma + 4] = -mu * ( deltaQ[s.v + 2] - w * alpha_x); // AQx(6) = - mu *Qx(9) + mu *u(3)*Qx(13)
159 BgradQ[s.sigma + 5] = 0.0; // AQx(5) = 0.0
160 }
161 break;
162
163 //dy
164 case 1:
165 {
166 double alpha_y = deltaQ[s.alpha];
167 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)
168 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)
169 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)
170 BgradQ[s.sigma + 0] = -lambda * ( deltaQ[s.v + 1] - v * alpha_y); // BQy(1) = - lam*Qy(8) + lam*u(2)*Qy(13)
171 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)
172 BgradQ[s.sigma + 2] = -lambda * ( deltaQ[s.v + 1] - v * alpha_y); // BQy(3) = - lam*Qy(8) + lam*u(2)*Qy(13)
173 BgradQ[s.sigma + 3] = -mu * ( deltaQ[s.v + 0] - u * alpha_y); // BQy(4) = - mu *Qy(7) + mu *u(1)*Qy(13)
174 BgradQ[s.sigma + 4] = 0.0; // BQy(6) = 0.0
175 BgradQ[s.sigma + 5] = -mu * ( deltaQ[s.v + 2] - w * alpha_y); // BQy(5) = - mu *Qy(9) + mu *u(3)*Qy(13)
176 }
177 break;
178 // !
179 //dz
180 case 2:
181 {
182 double alpha_z = deltaQ[s.alpha];
183 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)
184 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)
185 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)
186 BgradQ[s.sigma + 0] = -lambda * ( deltaQ[s.v + 2] - w * alpha_z); // CQz(1) = - lam*Qz(9) + lam*u(3)*Qz(13)
187 BgradQ[s.sigma + 1] = -lambda * ( deltaQ[s.v + 2] - w * alpha_z); // CQz(2) = - lam*Qz(9) + lam*u(3)*Qz(13)
188 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)
189 BgradQ[s.sigma + 3] = 0.0; // CQz(4) = 0.0
190 BgradQ[s.sigma + 4] = -mu * ( deltaQ[s.v + 0] - u * alpha_z); // CQz(6) = - mu *Qz(7) + mu *u(1)*Qz(13)
191 BgradQ[s.sigma + 5] = -mu * ( deltaQ[s.v + 1] - v * alpha_z); // CQz(5) = - mu *Qz(8) + mu *u(2)*Qz(13)
192 }
193 break;
194
195 default:
196#ifdef Asserts
197 assertionMsg(false,"Normal direction/gradQ component must be in {0,1,2}.");
198#endif
199 break;
200 }
201
202#ifdef Asserts
203 assertion2(std::isfinite(BgradQ[s.v + 0]),"BGradQ: v not finite",BgradQ[s.v + 0]);
204 assertion2(std::isfinite(BgradQ[s.v + 1]),"BGradQ: u not finite",BgradQ[s.v + 1]);
205 assertion2(std::isfinite(BgradQ[s.v + 2]),"BGradQ: w not finite",BgradQ[s.v + 2]);
206 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);
207 assertion2(std::isfinite(BgradQ[s.sigma + 1]),"BGradQ: sigma_yy not finite",BgradQ[s.sigma + 1]);
208 assertion2(std::isfinite(BgradQ[s.sigma + 2]),"BGradQ: sigma_zz not finite",BgradQ[s.sigma + 2]);
209 assertion2(std::isfinite(BgradQ[s.sigma + 3]),"BGradQ: sigma_xy not finite",BgradQ[s.sigma + 3]);
210 assertion2(std::isfinite(BgradQ[s.sigma + 4]),"BGradQ: sigma_xz not finite",BgradQ[s.sigma + 4]);
211 assertion2(std::isfinite(BgradQ[s.sigma + 5]),"BGradQ: sigma_yz not finite",BgradQ[s.sigma + 5]);
212#endif
213 logTraceOut( "nonconservativeProduct(...)" );
214}
#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 logTraceInWith3Arguments(methodName, argument0, argument1, argument2)
Definition Log.h:372
#define logTraceIn(methodName)
Definition Log.h:369
#define logTraceInWith2Arguments(methodName, argument0, argument1)
Definition Log.h:371
double maxEigenvalue(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, int normal) override
static tarch::logging::Log _log
Definition ADERDGLOH1.h:26
void nonconservativeProduct(const double *__restrict__ Q, const double *__restrict__ deltaQ, const tarch::la::Vector< Dimensions, double > &x, double t, int normal, double *__restrict__ BgradQ) override
void prescribeGaussianWave(const tarch::la::Vector< Dimensions, double > &x, double Q[])
Definition ADERDGLOH1.cpp:7
void adjustSolution(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t) override
struct examples::exahype2::loh1::ADERDGLOH1::VariablesShortcuts s
virtual void boundaryConditions(const double *__restrict__ Qinside, double *__restrict__ Qoutside, const tarch::la::Vector< Dimensions, double > &x, double t, int normal) 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