Peano 4
Loading...
Searching...
No Matches
LOH1OnGPU.cpp
Go to the documentation of this file.
1#include "LOH1OnGPU.h"
3
4
5tarch::logging::Log examples::exahype2::loh1::LOH1OnGPU::_log( "examples::exahype2::loh1::LOH1OnGPU" );
6
8 double * __restrict__ Q,
11 double t,
12 double dt
13)
14{
15 // initial conditions
16 if (tarch::la::equals(t,0.0) )
17 {
18 double center_curve[3] = {0.0,0.0,2.0};
19 bool upperLayer = volumeX(2) < 1.0;
20
21 double rho = upperLayer ? 4.0 : 6.0;
22 double cp = upperLayer ? 2.0 : 3.464;
23 double cs = upperLayer ? 2.6 : 2.7;
24
25 constexpr static int v_ = 0;
26 constexpr static int sigma_ = 3;
27 constexpr static int rho_ = 9;
28 constexpr static int cp_ = 10;
29 constexpr static int cs_ = 11;
30 constexpr static int alpha_ = 12;
31
32 Q[rho_ ] = rho;
33 Q[cp_ ] = cp;
34 Q[cs_ ] = cs;
35 Q[alpha_] = 1.0;
36
37 double radius = 0.1;
38
39 Q[ v_ + 0 ] = std::exp(-( (volumeX[0]-center_curve[0])*(volumeX[0]-center_curve[0])+
40 (volumeX[1]-center_curve[1])*(volumeX[1]-center_curve[1])+
41 (volumeX[2]-center_curve[2])*(volumeX[2]-center_curve[2]))/radius) * radius;
42 Q[ v_ + 1 ] = Q[ v_ + 0 ];
43 Q[ v_ + 2 ] = Q[ v_ + 1 ];
44
45 for(int i= 0 ; i < 6 ; i++) Q[ sigma_ + i ] = 0.0;
46 }
47 else
48 {
49 // other stuff
50 }
51 logTraceOut( "adjustSolution(...)" );
52}
53
54#if defined(GPUOffloadingOMP)
55#pragma omp declare target
56#endif
58 const double * __restrict__ Qinside,
59 double * __restrict__ Qoutside,
62 double t,
63 int normal
64)
65{
66 std::fill_n(Qoutside, 13, 0);
67
68 constexpr static int rho_ = 9;
69 constexpr static int cp_ = 10;
70 constexpr static int cs_ = 11;
71 constexpr static int alpha_ = 12;
72
73 Qoutside[ alpha_ ] = Qinside[ alpha_ ];
74 Qoutside[ rho_ ] = Qinside[ rho_ ];
75 Qoutside[ cp_ ] = Qinside[ cp_ ];
76 Qoutside[ cs_ ] = Qinside[ cs_ ];
77}
78#if defined(GPUOffloadingOMP)
79#pragma omp end declare target
80#endif
81
82
83#if defined(GPUOffloadingOMP)
84#pragma omp declare target
85#endif
87 const double * __restrict__ Q, // Q[9+4],
90 double t,
91 int normal
92)
93{
94 constexpr static int cp_ = 10;
95 constexpr static int cs_ = 11;
96 double cp = Q[cp_];
97 double cs = Q[cs_];
98
99/*
100 lambda[0] = -cp;
101 lambda[1] = cp;
102 lambda[2] = -cs;
103 lambda[3] = cp; I think this is a bug!
104 lambda[4] = 0.0;
105 lambda[5] = 0.0;
106 lambda[6] = 0.0;
107 lambda[7] = 0.0;
108 lambda[8] = 0.0;
109 lambda[9] = 0.0;
110*/
111
112 double result = 0.0;
113 result = std::max( result, -cp );
114 result = std::max( result, cp );
115 result = std::max( result, -cs );
116 result = std::max( result, cs );
117 return result;
118}
119#if defined(GPUOffloadingOMP)
120#pragma omp end declare target
121#endif
122
123
124
125#if defined(GPUOffloadingOMP)
126#pragma omp declare target
127#endif
129 const double * __restrict__ Q, // Q[9+4],
130 const double * __restrict__ deltaQ, // [9+4]
131 const tarch::la::Vector<Dimensions,double>& faceCentre,
133 double t,
134 int normal,
135 double * __restrict__ BgradQ // BgradQ[9]
136)
137{
138
139 constexpr static int v_ = 0;
140 constexpr static int sigma_ = 3;
141 constexpr static int rho_ = 9;
142 constexpr static int cp_ = 10;
143 constexpr static int cs_ = 11;
144 constexpr static int alpha_ = 12;
145
146 double rho = Q[rho_];
147 double cp = Q[cp_];
148 double cs = Q[cs_];
149 double alpha = Q[alpha_];
150
151 double mu = rho*cs*cs;
152 double lambda = rho*cp*cp-2*mu;
153 double rho_i = 1.0/rho;
154
155
156 double alpha_i = alpha/(alpha*alpha + 1.0e-13 * (1.0-alpha));
157
158 double u= Q[v_ + 0] * alpha_i;
159 double v= Q[v_ + 1] * alpha_i;
160 double w= Q[v_ + 2] * alpha_i;
161
162
163 switch ( normal ) {
164 //dx
165 case 0:
166 {
167 double alpha_x = deltaQ[alpha_];
168 BgradQ[v_ + 0] = -rho_i*(deltaQ[sigma_ + 0] + 2*Q[sigma_ + 0]*alpha_x);
169 BgradQ[v_ + 1] = -rho_i*(deltaQ[sigma_ + 3] + 2*Q[sigma_ + 3]*alpha_x);
170 BgradQ[v_ + 2] = -rho_i*(deltaQ[sigma_ + 4] + 2*Q[sigma_ + 4]*alpha_x);
171
172 BgradQ[sigma_ + 0] = -(lambda+2*mu) * ( deltaQ[v_ + 0] - u * alpha_x);
173 BgradQ[sigma_ + 1] = -lambda * ( deltaQ[v_ + 0] - u * alpha_x);
174 BgradQ[sigma_ + 2] = -lambda * ( deltaQ[v_ + 0] - u * alpha_x);
175 BgradQ[sigma_ + 3] = -mu * ( deltaQ[v_ + 1] - v * alpha_x);
176 BgradQ[sigma_ + 4] = -mu * ( deltaQ[v_ + 2] - w * alpha_x);
177 BgradQ[sigma_ + 5] = 0.0;
178 }
179 break;
180
182 case 1:
183 {
184 double alpha_y = deltaQ[alpha_];
185 BgradQ[v_ + 0] = -rho_i*(deltaQ[sigma_ + 3] + 2*Q[sigma_ + 3]*alpha_y);
186 BgradQ[v_ + 1] = -rho_i*(deltaQ[sigma_ + 1] + 2*Q[sigma_ + 1]*alpha_y);
187 BgradQ[v_ + 2] = -rho_i*(deltaQ[sigma_ + 5] + 2*Q[sigma_ + 5]*alpha_y);
188
189 BgradQ[sigma_ + 0] = -lambda * ( deltaQ[v_ + 1] - v * alpha_y);
190 BgradQ[sigma_ + 1] = -(lambda+2*mu) * ( deltaQ[v_ + 1] - v * alpha_y);
191 BgradQ[sigma_ + 2] = -lambda * ( deltaQ[v_ + 1] - v * alpha_y);
192 BgradQ[sigma_ + 3] = -mu * ( deltaQ[v_ + 0] - u * alpha_y);
193 BgradQ[sigma_ + 4] = 0.0;
194 BgradQ[sigma_ + 5] = -mu * ( deltaQ[v_ + 2] - w * alpha_y);
195 }
196 break;
199 case 2:
200 {
201 double alpha_z = deltaQ[alpha_];
202 BgradQ[v_ + 0] = -rho_i*(deltaQ[sigma_ + 4] + 2*Q[sigma_ + 4]*alpha_z);
203 BgradQ[v_ + 1] = -rho_i*(deltaQ[sigma_ + 5] + 2*Q[sigma_ + 5]*alpha_z);
204 BgradQ[v_ + 2] = -rho_i*(deltaQ[sigma_ + 2] + 2*Q[sigma_ + 2]*alpha_z);
205 BgradQ[sigma_ + 0] = -lambda * ( deltaQ[v_ + 2] - w * alpha_z);
206 BgradQ[sigma_ + 1] = -lambda * ( deltaQ[v_ + 2] - w * alpha_z);
207 BgradQ[sigma_ + 2] = -(lambda+2*mu) * ( deltaQ[v_ + 2] - w * alpha_z);
208 BgradQ[sigma_ + 3] =
209 BgradQ[sigma_ + 4] = -mu * ( deltaQ[v_ + 0] - u * alpha_z);
210 BgradQ[sigma_ + 5] = -mu * ( deltaQ[v_ + 1] - v * alpha_z);
211 }
212 break;
213
214 default:
215 break;
216 }
217}
218#if defined(GPUOffloadingOMP)
219#pragma omp end declare target
220#endif
#define logTraceOut(methodName)
Definition Log.h:379
virtual void boundaryConditions(const double *__restrict__ Qinside, double *__restrict__ Qoutside, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, int normal) override
Apply boundary conditions.
Definition LOH1OnGPU.cpp:57
static double maxEigenvalue(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, int normal)
Definition LOH1OnGPU.cpp:86
static tarch::logging::Log _log
Definition LOH1OnGPU.h:26
static void nonconservativeProduct(const double *__restrict__ Q, const double *__restrict__ deltaQ, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, int normal, double *__restrict__ BgradQ)
virtual void adjustSolution(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &volumeCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, double dt) override
Feel free to change the solution in the particular finite volume.
Definition LOH1OnGPU.cpp:7
Log Device.
Definition Log.h:516
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