Peano 4
Loading...
Searching...
No Matches
NavierStokesSolver.cpp
Go to the documentation of this file.
2
4 "applications::exahype2:"
5 ":CompressibleNavierStok"
6 "es::NavierStokesSolver"
7);
8
22
23Scenario scenario = Scenario::CJDetonation;
24
26 double* __restrict__ Q,
27 const tarch::la::Vector<Dimensions, double>& volumeCentre,
29 bool gridIsConstructed
30) {
31 assertion(!std::isnan(Q[rho]));
32 assertion(!std::isnan(Q[u]));
33 assertion(!std::isnan(Q[v]));
34 assertion(!std::isnan(Q[e]));
35
36 switch (scenario) {
37
38 case Scenario::LidDrivenCavity: {
39 Q[rho] = 1.0;
40 Q[u] = 0.0;
41 Q[v] = 0.0;
42#if Dimensions == 3
43 Q[w] = 0.0;
44#endif
45 Q[e] = (100 / Gamma) / (Gamma - 1) + 0.5 * (Q[u] * Q[u] + Q[v] * Q[v]) / Q[rho];
46 Q[Z] = 0.0;
47 break;
48 }
49
50 case Scenario::PointExplosion: {
51#if Dimensions == 2
52 const tarch::la::Vector<2, double> circleCentre = {0.5, 0.5};
53 bool isInTheCentre = (tarch::la::norm2(volumeCentre - circleCentre) < 0.2);
54#endif
55#if Dimensions == 3
56 const tarch::la::Vector<3, double> circleCentre = {0.5, 0.5, 0.5};
57 bool isInTheCentre = (tarch::la::norm2(volumeCentre - circleCentre) < 0.2);
58#endif
59 Q[rho] = 1.0;
60 Q[u] = 0;
61 Q[v] = 0;
62#if Dimensions == 3
63 Q[w] = 0.0;
64#endif
65 Q[e] = isInTheCentre ? 1.00 : 1.01; // inner energy
66 break;
67 }
68
69 case Scenario::TaylorGreenVortex: {
70#if Dimensions == 2
71 const double pressure = 0.25 * (std::cos(4.0 * M_PI * volumeCentre[0]) + std::sin(4.0 * M_PI * volumeCentre[1]))
72 + 100 / Gamma;
73 Q[rho] = 1.0;
74 Q[u] = std::sin(2.0 * M_PI * volumeCentre[0]) * std::cos(2.0 * M_PI * volumeCentre[1]);
75 Q[v] = -std::cos(2.0 * M_PI * volumeCentre[0]) * std::sin(2.0 * M_PI * volumeCentre[1]);
76 Q[e] = pressure / (Gamma - 1) + 0.5 * (Q[u] * Q[u] + Q[v] * Q[v]) / Q[rho];
77#else
78 std::cout << "Taylor green vortex only works for 2D!\n";
79 exit(0);
80#endif
81 break;
82 }
83
84 case Scenario::ShockTube: {
85 double pressure;
86 if (volumeCentre[0] < 0.45) {
87 Q[rho] = 0.125;
88 pressure = 0.1;
89 } else {
90 Q[rho] = 1.0;
91 pressure = 1.0;
92 }
93 Q[u] = 0;
94 Q[v] = 0;
95#if Dimensions == 2
96 Q[e] = pressure / (Gamma - 1) + 0.5 * (Q[u] * Q[u] + Q[v] * Q[v]) / Q[rho];
97#else
98 Q[w] = 0;
99 Q[e] = pressure / (Gamma - 1) + 0.5 * (Q[u] * Q[u] + Q[v] * Q[v] + Q[w] * Q[w]) / Q[rho];
100#endif
101 break;
102 }
103
104 case Scenario::DoubleShockTube: {
105 double pressure;
106 if (volumeCentre[0] < 0.33 || volumeCentre[0] > 0.66) {
107 Q[rho] = 0.125;
108 pressure = 0.1;
109 } else {
110 Q[rho] = 1.0;
111 pressure = 1.0;
112 }
113 Q[u] = 0;
114 Q[v] = 0;
115#if Dimensions == 2
116 Q[e] = pressure / (Gamma - 1) + 0.5 * (Q[u] * Q[u] + Q[v] * Q[v]) / Q[rho];
117#else
118 Q[w] = 0;
119 Q[e] = pressure / (Gamma - 1) + 0.5 * (Q[u] * Q[u] + Q[v] * Q[v] + Q[w] * Q[w]) / Q[rho];
120#endif
121 break;
122 }
123
124 case Scenario::SmoothWave: {
125 const auto distX = volumeCentre[0] - 0.5;
126 const auto distY = volumeCentre[1] - 0.5;
127 const auto dist = distX * distX + distY * distY;
128 const double pressure = 1.0 - dist;
129
130 Q[rho] = 1.0 - dist;
131 Q[u] = 0.0;
132 Q[v] = 0.0;
133 Q[e] = pressure / (Gamma - 1) + 0.5 * (Q[u] * Q[u] + Q[v] * Q[v]) / Q[rho];
134#if Dimensions == 3
135 Q[w] = 0.0;
136 Q[e] += 0.5 * Q[w] * Q[w] / Q[rho];
137#endif
138 break;
139 }
140
141 case Scenario::CJDetonation: {
142
143#if Dimensions == 2
144 // B: Burned, U: Unburned
145 const auto rhoB = 1.4;
146 const auto rhoU = 0.887565;
147 const auto pB = 1.0;
148 const auto pU = 0.191709;
149 const auto ZB = 0.0;
150 const auto ZU = 1.0;
151
152 const auto alpha = std::atan2(volumeCentre[1], volumeCentre[0]); // Angle in polar coordinates
153 const auto uB = 0.0;
154 const auto uU = -0.577350 * std::cos(alpha);
155 const auto vB = 0.0;
156 const auto vU = -0.577350 * std::sin(alpha);
157
158 const auto distX = volumeCentre[0];
159 const auto distY = volumeCentre[1];
160 const auto distToCenter = std::sqrt(distX * distX + distY * distY);
161 const auto radiusBurned = 0.3;
162 const auto isBurned = distToCenter < radiusBurned;
163
164 if (isBurned) {
165 Q[rho] = rhoB;
166 Q[u] = uB * rhoB;
167 Q[v] = vB * rhoB;
168 Q[Z] = ZB * rhoB;
169 Q[e] = pB / (Gamma - 1) + 0.5 * (Q[u] * Q[u] + Q[v] * Q[v]) / Q[rho] + Q[Z];
170 } else {
171 Q[rho] = rhoU;
172 Q[u] = uU * rhoU;
173 Q[v] = vU * rhoU;
174 Q[Z] = ZU * rhoU;
175 Q[e] = pU / (Gamma - 1) + 0.5 * (Q[u] * Q[u] + Q[v] * Q[v]) / Q[rho] + Q[Z];
176 }
177 break;
178#else
179 std::cout << "CJ Detonation Wave Only works in 2D!\n";
180 exit(0);
181#endif
182 }
183 }
184}
185
187 const double* __restrict__ Qinside,
188 double* __restrict__ Qoutside,
191 double t,
192 int normal
193) {
194 assertion(!std::isnan(Qinside[rho]));
195 assertion(!std::isnan(Qinside[u]));
196 assertion(!std::isnan(Qinside[v]));
197 assertion(!std::isnan(Qinside[e]));
198
199 switch (scenario) {
200 case Scenario::LidDrivenCavity:
201 Qoutside[rho] = Qinside[rho];
202 Qoutside[u] = 0;
203 if (faceCentre[y] >= DomainSize[y] - std::numeric_limits<double>::epsilon()) { // TODO: DomainOffset missing here
204 Qoutside[u] = LidDrivenCavityVelocity;
205 }
206 Qoutside[v] = 0;
207#if Dimensions == 3
208 Qoutside[w] = 0;
209#endif
210 Qoutside[e] = Qinside[e];
211 break;
212
213 case Scenario::CJDetonation:
214 case Scenario::ShockTube:
215 case Scenario::DoubleShockTube:
216 case Scenario::PointExplosion:
217 for (int i = 0; i < NumberOfUnknowns + NumberOfAuxiliaryVariables; i++) {
218 Qoutside[i] = Qinside[i];
219 }
220 break;
221
222 case Scenario::SmoothWave:
223 for (int i = 0; i < NumberOfUnknowns + NumberOfAuxiliaryVariables; i++) {
224 Qoutside[i] = Qinside[i];
225 }
226 if (normal != 0) {
227 Qoutside[u + normal] = -Qinside[u + normal];
228 }
229 break;
230 case Scenario::TaylorGreenVortex:
231 const double pressure = 0.25 * std::exp(-4 * Viscosity * t)
232 * (std::cos(4.0 * M_PI * faceCentre[0]) + std::sin(4.0 * M_PI * faceCentre[1]))
233 + 100 / Gamma;
234
235 Qoutside[rho] = 1.0;
236 Qoutside[u] = std::exp(-2 * Viscosity * t) * std::sin(2.0 * M_PI * faceCentre[0])
237 * std::cos(2.0 * M_PI * faceCentre[1]);
238 Qoutside[v] = std::exp(-2 * Viscosity * t) * -std::cos(2.0 * M_PI * faceCentre[0])
239 * std::sin(2.0 * M_PI * faceCentre[1]);
240 Qoutside[e] = pressure / (Gamma - 1) + 0.5 * (Qoutside[u] * Qoutside[u] + Qoutside[v] * Qoutside[v]) / Qoutside[rho];
241 break;
242 }
243}
244
246 const double* __restrict__ Q,
249 double t,
250 double dt,
251 int normal
252) {
253 return maxEigenvalue(Q, faceCentre, volumeH, t, dt, normal, Offloadable::Yes);
254}
255
257 const double* __restrict__ Q,
260 double t,
261 double dt,
262 int normal,
263 double* __restrict__ F
264) {
265 flux(Q, faceCentre, volumeH, t, dt, normal, F, Offloadable::Yes);
266}
267
269 const double* __restrict__ Q,
270 const double* __restrict__ deltaQ,
273 double t,
274 double dt,
275 int normal,
276 double* __restrict__ BTimesDeltaQ
277) {
278 nonconservativeProduct(Q, deltaQ, faceCentre, volumeH, t, dt, normal, BTimesDeltaQ, Offloadable::Yes);
279}
280
282 const double* __restrict__ Q,
283 const tarch::la::Vector<Dimensions, double>& volumeCentre,
285 double t,
286 double dt,
287 double* __restrict__ S
288) {
289 if (scenario == Scenario::CJDetonation) {
290 const double irho = 1.0 / Q[rho];
291#if Dimensions == 2
292 const double pressure = (Gamma - 1) * (Q[e] - 0.5 * irho * (Q[u] * Q[u] + Q[v] * Q[v]) - Q[Z]);
293#else
294 const double pressure = (Gamma - 1) * (Q[e] - 0.5 * irho * (Q[u] * Q[u] + Q[v] * Q[v] + Q[w] * Q[w]) - Q[Z]);
295#endif
296 const double T = pressure * irho / R;
297
298 S[Z] = -Q[Z] * (T > IgnitionTemperature ? -1.0 / Timescale : 0.0);
299 }
300}
301
303 const double* __restrict__ Q,
304 const tarch::la::Vector<Dimensions, double>& volumeCentre,
306 double t
307) {
309
310 switch (scenario) {
311 case Scenario::LidDrivenCavity:
312 result = volumeCentre[y] > DomainSize[y] / 2 ? ::exahype2::RefinementCommand::Refine : result;
313 break;
314 }
315
316 return result;
317}
#define assertion(expr)
@ LidDrivenCavity
@ RayleighTaylorInstability
@ TaylorGreenVortex
@ PointExplosion
@ ChannelFlowWithObstacle
@ DoubleShockTube
virtual void sourceTerm(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &volumeCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, double dt, double *__restrict__ S) override
virtual void flux(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, double dt, int normal, double *__restrict__ F) override
virtual void initialCondition(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &volumeCentre, const tarch::la::Vector< Dimensions, double > &volumeH, bool gridIsConstructed) override
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
virtual::exahype2::RefinementCommand refinementCriterion(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &volumeCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t) override
virtual double maxEigenvalue(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, double dt, int normal) override
virtual 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, double dt, int normal, double *__restrict__ BTimesDeltaQ) override
Log Device.
Definition Log.h:516
Scalar norm2(const Vector< Size, Scalar > &vector)
Computes the 2-norm of the vector, i.e.
Simple vector class.
Definition Vector.h:134