Peano 4
Loading...
Searching...
No Matches
NavierStokesSolver.h
Go to the documentation of this file.
1// This file is part of the Peano project. For conditions of distribution and
2// use, please see the copyright notice at www.peano-framework.org
3#pragma once
4
5#include "AbstractNavierStokesSolver.h"
6#include "peano4/utils/Loop.h"
8#include "tarch/logging/Log.h"
9
12 x = 0,
14#if Dimensions == 3
15 z
16#endif
17 };
18
19 enum Variables {
20 rho = 0,
23#if Dimensions == 3
24 w,
25#endif
27 Z
28 };
29
31 dudx = 0,
33#if Dimensions == 3
34 dudz,
35#endif
38#if Dimensions == 3
39 dvdz,
40 dwdx,
41 dwdy,
42 dwdz,
43#endif
45 dTdy
46#if Dimensions == 3
47 ,
48 dTdz
49#endif
50 };
51 class NavierStokesSolver;
52} // namespace applications::exahype2::CompressibleNavierStokes
53
55private:
57
58public:
59 virtual void initialCondition(
60 double* __restrict__ Q,
61 const tarch::la::Vector<Dimensions, double>& volumeCentre,
63 bool gridIsConstructed
64 ) override;
65
66 virtual void boundaryConditions(
67 const double* __restrict__ Qinside,
68 double* __restrict__ Qoutside,
71 double t,
72 int normal
73 ) override;
74
75 virtual ::exahype2::RefinementCommand refinementCriterion(
76 const double* __restrict__ Q,
77 const tarch::la::Vector<Dimensions, double>& volumeCentre,
79 double t
80 ) override;
81
82 virtual double maxEigenvalue(
83 const double* __restrict__ Q,
86 double t,
87 double dt,
88 int normal
89 ) override;
90
91 static inline double maxEigenvalue(
92 const double* __restrict__ Q,
95 double t,
96 double dt,
97 int normal,
98 Offloadable
99 ) {
100
101#ifdef GPUOffloadingOff
102 assertion(normal >= 0);
103 assertion(normal < Dimensions);
104 assertion(Q[rho] == Q[0]);
105 assertion(Q[rho] > 0);
106 assertion(Q[e] > 0);
107 assertion((v + normal) < (NumberOfUnknowns + NumberOfAuxiliaryVariables));
108 assertion(!std::isnan(Q[rho]));
109 assertion(!std::isnan(Q[u]));
110 assertion(!std::isnan(Q[v]));
111#if Dimensions == 3
112 assertion(!std::isnan(Q[w]));
113#endif
114 assertion(!std::isnan(Q[e]));
115#endif
116
117 const double irho = 1.0 / Q[rho];
118#if Dimensions == 2
119 const double pressure = (Gamma - 1) * (Q[e] - 0.5 * irho * (Q[u] * Q[u] + Q[v] * Q[v]) - Q[Z]);
120#else
121 const double pressure = (Gamma - 1) * (Q[e] - 0.5 * irho * (Q[u] * Q[u] + Q[v] * Q[v] + Q[w] * Q[w]) - Q[Z]);
122#endif
123 assertion(pressure > 0);
124
125 const double c_speedOfSound = std::sqrt(Gamma * pressure / Q[rho]);
126
127 const double max_lambda_convective = std::abs(Q[u + normal]) / Q[rho] + c_speedOfSound;
128 const double max_lambda_viscous = std::max(
129 4.0 / 3.0 * Viscosity / Q[rho], Gamma * Viscosity / (PrandtlNumber * Q[rho])
130 );
131 return std::max(max_lambda_convective, max_lambda_viscous);
132 }
133
134 virtual void flux(
135 const double* __restrict__ Q,
138 double t,
139 double dt,
140 int normal,
141 double* __restrict__ F
142 ) override;
143
144 static inline void flux(
145 const double* __restrict__ Q,
148 double t,
149 double dt,
150 int normal,
151 double* __restrict__ F,
152 Offloadable
153 ) {
154
155#ifdef GPUOffloadingOff
156 assertion3(normal >= 0, faceCentre, t, normal);
157 assertion3(normal < Dimensions, faceCentre, t, normal);
158 assertion(Q[rho] > 0);
159 assertion(Q[e] > 0);
160 assertion((v + normal) < (NumberOfUnknowns + NumberOfAuxiliaryVariables));
161 assertion(!std::isnan(Q[rho]));
162 assertion(!std::isnan(Q[u]));
163 assertion(!std::isnan(Q[v]));
164#if Dimensions == 3
165 assertion(!std::isnan(Q[w]));
166#endif
167 assertion(!std::isnan(Q[e]));
168 assertion(!std::isnan(Q[Z]));
169#endif
170
171 const double irho = 1.0 / Q[rho];
172#if Dimensions == 2
173 const double pressure = (Gamma - 1) * (Q[e] - 0.5 * irho * (Q[u] * Q[u] + Q[v] * Q[v]) - Q[Z]);
174#else
175 const double pressure = (Gamma - 1) * (Q[e] - 0.5 * irho * (Q[u] * Q[u] + Q[v] * Q[v] + Q[w] * Q[w]) - Q[Z]);
176#endif
177
178#ifdef GPUOffloadingOff
179 assertion(pressure > 0);
180#endif
181
182 F[rho] = Q[u + normal];
183 F[u] = irho * Q[u] * Q[u + normal];
184 F[v] = irho * Q[v] * Q[u + normal];
185#if Dimensions == 3
186 F[w] = irho * Q[w] * Q[u] + normal;
187#endif
188 F[u + normal] += pressure;
189 F[e] = irho * Q[u + normal] * (Q[e] + pressure);
190 F[Z] = irho * Q[u + normal] * Q[Z];
191 }
192
193 virtual void nonconservativeProduct(
194 const double* __restrict__ Q,
195 const double* __restrict__ deltaQ,
198 double t,
199 double dt,
200 int normal,
201 double* __restrict__ BTimesDeltaQ
202 ) override;
203
204 static inline void nonconservativeProduct(
205 const double* __restrict__ Q,
206 const double* __restrict__ deltaQ,
209 double t,
210 double dt,
211 int normal,
212 double* __restrict__ BTimesDeltaQ,
213 Offloadable
214 ) {
215
216#ifdef GPUOffloadingOff
217 assertion(normal >= 0);
218 assertion(normal < Dimensions);
219 assertion(Q[rho] > 0);
220 assertion(!std::isnan(Q[rho]));
221 assertion(!std::isnan(Q[u]));
222 assertion(!std::isnan(Q[v]));
223 assertion(!std::isnan(Q[e]));
224 assertion(!std::isnan(Q[Z]));
225 assertion(!std::isnan(Q[NumberOfUnknowns + dudx]));
226 assertion(!std::isnan(Q[NumberOfUnknowns + dudy]));
227 assertion(!std::isnan(Q[NumberOfUnknowns + dvdx]));
228 assertion(!std::isnan(Q[NumberOfUnknowns + dvdy]));
229 assertion(!std::isnan(Q[NumberOfUnknowns + dTdx]));
230 assertion(!std::isnan(Q[NumberOfUnknowns + dTdy]));
231#if Dimensions == 3
232 assertion(!std::isnan(Q[w]));
233 assertion(!std::isnan(Q[NumberOfUnknowns + dudz]));
234 assertion(!std::isnan(Q[NumberOfUnknowns + dvdz]));
235 assertion(!std::isnan(Q[NumberOfUnknowns + dwdx]));
236 assertion(!std::isnan(Q[NumberOfUnknowns + dwdy]));
237 assertion(!std::isnan(Q[NumberOfUnknowns + dwdz]));
238 assertion(!std::isnan(Q[NumberOfUnknowns + dTdz]));
239#endif
240#endif
241
242 double kappa = Viscosity * c_p / PrandtlNumber;
243 double c_v = c_p / Gamma;
244 double derivative[NumberOfAuxiliaryVariables];
245
246 auto irho2 = 1.0 / (Q[rho] * Q[rho]);
247 auto irho3 = irho2 / Q[rho];
248 auto a = -Q[e] * irho3 + (Q[u] * Q[u] + Q[v] * Q[v]) * irho3;
249 auto b = -Q[u] * irho2;
250 auto c = -Q[v] * irho2;
251
252 if (normal == x) {
253 derivative[dudx] = deltaQ[u];
254 derivative[dvdx] = deltaQ[v];
255 derivative[dudy] = Q[NumberOfUnknowns + dudy] - 0.5 * deltaQ[NumberOfUnknowns + dudy];
256 derivative[dvdy] = Q[NumberOfUnknowns + dvdy] - 0.5 * deltaQ[NumberOfUnknowns + dvdy];
257 derivative[dTdx]
258 = (a * deltaQ[rho] + b * derivative[dudx] + c * derivative[dvdx] - irho2 * (deltaQ[Z] - 2 * Q[Z] * deltaQ[rho]));
259 // derivative[dTdx] = Q[NumberOfUnknowns + dTdx] - 0.5 * deltaQ[NumberOfUnknowns + dTdx];
260
261#if Dimension == 3
262 derivative[dwdx] = deltaQ[w];
263 derivative[dudz] = Q[NumberOfUnknowns + dudz] - 0.5 * deltaQ[NumberOfUnknowns + dudz];
264 derivative[dwdz] = Q[NumberOfUnknowns + dwdz] - 0.5 * deltaQ[NumberOfUnknowns + dwdz];
265#endif
266
267 // -mu{2(du/dx) - 2/3 (du/dx + dv/dy)}
268 BTimesDeltaQ[u] = -Viscosity * (2.0 * derivative[dudx] - 2.0 / 3.0 * (derivative[dudx] + derivative[dvdy]
269#if Dimension == 3
270 +derivative[dwdz]
271#endif
272 )) / Q[rho];
273
274 // -mu{dv/dx + du/dy}
275 BTimesDeltaQ[v] = -Viscosity * (derivative[dvdx] + derivative[dudy]) / Q[rho];
276
277#if Dimension == 3
278 BTimesDeltaQ[w] = -Viscosity * (derivative[dudz] + derivative[dwdx]) / Q[rho];
279#endif
280 } else {
281#if Dimension == 3
282 if (normal == y) {
283 derivative[dwdy] = deltaQ[w];
284 derivative[dvdz] = Q[NumberOfUnknowns + dvdz] - 0.5 * deltaQ[NumberOfUnknowns + dvdz];
285 derivative[dwdz] = Q[NumberOfUnknowns + dwdz] - 0.5 * deltaQ[NumberOfUnknowns + dwdz];
286#endif
287
288 derivative[dudy] = deltaQ[u];
289 derivative[dvdy] = deltaQ[v];
290 derivative[dudx] = Q[NumberOfUnknowns + dudx] - 0.5 * deltaQ[NumberOfUnknowns + dudx];
291 derivative[dvdx] = Q[NumberOfUnknowns + dvdx] - 0.5 * deltaQ[NumberOfUnknowns + dvdx];
292
293 derivative[dTdy]
294 = (a * deltaQ[rho] + b * derivative[dudy] + c * derivative[dvdy] - irho2 * (deltaQ[Z] - 2 * Q[Z] * deltaQ[rho]));
295
296 // derivative[dTdy] = Q[NumberOfUnknowns + dTdy] - 0.5 * deltaQ[NumberOfUnknowns + dTdy];
297
298 // -mu{dv/dx + du/dy}
299 BTimesDeltaQ[u] = -Viscosity * (derivative[dvdx] + derivative[dudy]) / Q[rho];
300
301 // -mu{2(dv/dy) - 2/3 (du/dx + dv/dy)}
302 BTimesDeltaQ[v] = -Viscosity * (2.0 * derivative[dvdy] - 2.0 / 3.0 * (derivative[dudx] + derivative[dvdy]
303#if Dimension == 3
304 +derivative[dwdz]
305#endif
306 )) / Q[rho];
307#if Dimension == 3
308 BTimesDeltaQ[w] = -Viscosity * (derivative[dvdz] + derivative[dwdy]) / Q[rho];
309 } else {
310
311 derivative[dudz] = deltaQ[u];
312 derivative[dvdz] = deltaQ[v];
313 derivative[dwdz] = deltaQ[w];
314 derivative[dudx] = Q[NumberOfUnknowns + dudx] - 0.5 * deltaQ[NumberOfUnknowns + dudx];
315 derivative[dvdy] = Q[NumberOfUnknowns + dvdy] - 0.5 * deltaQ[NumberOfUnknowns + dvdy];
316 derivative[dudz] = Q[NumberOfUnknowns + dudz] - 0.5 * deltaQ[NumberOfUnknowns + dudz];
317 derivative[dwdy] = Q[NumberOfUnknowns + dwdy] - 0.5 * deltaQ[NumberOfUnknowns + dwdy];
318 derivative[dwdx] = Q[NumberOfUnknowns + dwdx] - 0.5 * deltaQ[NumberOfUnknowns + dwdx];
319 derivative[dTdz] = Q[NumberOfUnknowns + dTdz] - 0.5 * deltaQ[NumberOfUnknowns + dTdz];
320
321 BTimesDeltaQ[u] = -Viscosity * (derivative[dudz] + derivative[dwdx]) / Q[rho];
322 BTimesDeltaQ[v] = -Viscosity * (derivative[dvdz] + derivative[dwdy]) / Q[rho];
323 BTimesDeltaQ[w] = -Viscosity
324 * (2.0 * derivative[dwdz] - 2.0 / 3.0 * (derivative[dudx] + derivative[dvdy] + derivative[dwdz]))
325 / Q[rho];
326 }
327#endif
328 }
329 BTimesDeltaQ[rho] = 0;
330 BTimesDeltaQ[Z] = 0;
331#if Dimensions == 2
332 BTimesDeltaQ[e] = (Q[u] * BTimesDeltaQ[u] + Q[v] * BTimesDeltaQ[v]) / Q[rho] - kappa * derivative[dTdx + normal];
333#else
334 BTimesDeltaQ[e] = (Q[u] * BTimesDeltaQ[u] + Q[v] * BTimesDeltaQ[v] + Q[w] * BTimesDeltaQ[w]) / Q[rho]
335 - kappa * derivative[dTdx + normal];
336#endif
337 }
338
339
340 virtual void sourceTerm(
341 const double* __restrict__ Q,
342 const tarch::la::Vector<Dimensions, double>& volumeCentre,
344 double t,
345 double dt,
346 double* __restrict__ S
347 ) override;
348
349 static inline void extrapolateHalo(double* __restrict__ Q) {
350#if Dimensions == 2
351 tarch::la::Vector<Dimensions, int> topLeftCell = {0, NumberOfFiniteVolumesPerAxisPerPatch + 1};
353 NumberOfFiniteVolumesPerAxisPerPatch + 1, NumberOfFiniteVolumesPerAxisPerPatch + 1
354 };
355 tarch::la::Vector<Dimensions, int> bottomLeftCell = {0, 0};
356 tarch::la::Vector<Dimensions, int> bottomRightCell = {NumberOfFiniteVolumesPerAxisPerPatch + 1, 0};
357
358 tarch::la::Vector<Dimensions, int> topLeftCellN1 = topLeftCell;
359 topLeftCellN1(0)++;
360 tarch::la::Vector<Dimensions, int> topLeftCellN2 = topLeftCell;
361 topLeftCellN2(1)--;
362
363 tarch::la::Vector<Dimensions, int> topRightCellN1 = topRightCell;
364 topRightCellN1(0)--;
365 tarch::la::Vector<Dimensions, int> topRightCellN2 = topRightCell;
366 topRightCellN2(1)--;
367
368 tarch::la::Vector<Dimensions, int> bottomLeftCellN1 = bottomLeftCell;
369 bottomLeftCellN1(0)++;
370 tarch::la::Vector<Dimensions, int> bottomLeftCellN2 = bottomLeftCell;
371 bottomLeftCellN2(1)++;
372
373 tarch::la::Vector<Dimensions, int> bottomRightCellN1 = bottomRightCell;
374 bottomRightCellN1(0)--;
375 tarch::la::Vector<Dimensions, int> bottomRightCellN2 = bottomRightCell;
376 bottomRightCellN2(1)++;
377
378 const int topLeftCellSerialised = peano4::utils::dLinearised(topLeftCell, NumberOfFiniteVolumesPerAxisPerPatch + 2);
379 const int topRightCellSerialised = peano4::utils::dLinearised(
380 topRightCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
381 );
382
383 const int bottomLeftCellSerialised = peano4::utils::dLinearised(
384 bottomLeftCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
385 );
386 const int bottomRightCellSerialised = peano4::utils::dLinearised(
387 bottomRightCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
388 );
389
390 const int topLeftCellN1Serialised = peano4::utils::dLinearised(
391 topLeftCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
392 );
393 const int topLeftCellN2Serialised = peano4::utils::dLinearised(
394 topLeftCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
395 );
396
397 const int topRightCellN1Serialised = peano4::utils::dLinearised(
398 topRightCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
399 );
400 const int topRightCellN2Serialised = peano4::utils::dLinearised(
401 topRightCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
402 );
403
404 const int bottomLeftCellN1Serialised = peano4::utils::dLinearised(
405 bottomLeftCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
406 );
407 const int bottomLeftCellN2Serialised = peano4::utils::dLinearised(
408 bottomLeftCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
409 );
410
411 const int bottomRightCellN1Serialised = peano4::utils::dLinearised(
412 bottomRightCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
413 );
414 const int bottomRightCellN2Serialised = peano4::utils::dLinearised(
415 bottomRightCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
416 );
417
418 for (int i = 0; i < NumberOfUnknowns + NumberOfAuxiliaryVariables; i++) {
419
420 Q[(topLeftCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]
421 = 0.5
422 * (Q[(topLeftCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(topLeftCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]);
423 Q[(topRightCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]
424 = 0.5
425 * (Q[(topRightCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(topRightCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]);
426
427 Q[(bottomLeftCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]
428 = 0.5
429 * (Q[(bottomLeftCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(bottomLeftCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]);
430
431 Q[(bottomRightCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]
432 = 0.5
433 * (Q[(bottomRightCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(bottomRightCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]);
434 }
435#else
436 // Edges
437 for (int i = 1; i < NumberOfFiniteVolumesPerAxisPerPatch + 1; i++) {
438 tarch::la::Vector<Dimensions, int> backLeftCell = {0, i, 0};
439 tarch::la::Vector<Dimensions, int> backLeftCellN1 = {1, i, 0};
440 tarch::la::Vector<Dimensions, int> backLeftCellN2 = {0, i, 1};
441
442 const int backLeftCellSerialised = peano4::utils::dLinearised(
443 backLeftCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
444 );
445 const int backLeftCellN1Serialised = peano4::utils::dLinearised(
446 backLeftCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
447 );
448 const int backLeftCellN2Serialised = peano4::utils::dLinearised(
449 backLeftCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
450 );
451
452 tarch::la::Vector<Dimensions, int> backRightCell = {NumberOfFiniteVolumesPerAxisPerPatch + 1, i, 0};
453 tarch::la::Vector<Dimensions, int> backRightCellN1 = {NumberOfFiniteVolumesPerAxisPerPatch, i, 0};
454 tarch::la::Vector<Dimensions, int> backRightCellN2 = {NumberOfFiniteVolumesPerAxisPerPatch + 1, i, 1};
455
456 const int backRightCellSerialised = peano4::utils::dLinearised(
457 backRightCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
458 );
459 const int backRightCellN1Serialised = peano4::utils::dLinearised(
460 backRightCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
461 );
462 const int backRightCellN2Serialised = peano4::utils::dLinearised(
463 backRightCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
464 );
465
466 tarch::la::Vector<Dimensions, int> backTopCell = {i, NumberOfFiniteVolumesPerAxisPerPatch + 1, 0};
467 tarch::la::Vector<Dimensions, int> backTopCellN1 = {i, NumberOfFiniteVolumesPerAxisPerPatch + 1, 1};
468 tarch::la::Vector<Dimensions, int> backTopCellN2 = {i, NumberOfFiniteVolumesPerAxisPerPatch, 0};
469
470 const int backTopCellSerialised = peano4::utils::dLinearised(
471 backTopCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
472 );
473 const int backTopCellN1Serialised = peano4::utils::dLinearised(
474 backTopCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
475 );
476 const int backTopCellN2Serialised = peano4::utils::dLinearised(
477 backTopCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
478 );
479
480 tarch::la::Vector<Dimensions, int> backBottomCell = {i, 0, 0};
481 tarch::la::Vector<Dimensions, int> backBottomCellN1 = {i, 0, 1};
482 tarch::la::Vector<Dimensions, int> backBottomCellN2 = {i, 1, 0};
483
484 const int backBottomCellSerialised = peano4::utils::dLinearised(
485 backBottomCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
486 );
487 const int backBottomCellN1Serialised = peano4::utils::dLinearised(
488 backBottomCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
489 );
490 const int backBottomCellN2Serialised = peano4::utils::dLinearised(
491 backBottomCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
492 );
493
494 tarch::la::Vector<Dimensions, int> frontLeftCell = {0, i, NumberOfFiniteVolumesPerAxisPerPatch + 1};
495 tarch::la::Vector<Dimensions, int> frontLeftCellN1 = {1, i, NumberOfFiniteVolumesPerAxisPerPatch + 1};
496 tarch::la::Vector<Dimensions, int> frontLeftCellN2 = {0, i, NumberOfFiniteVolumesPerAxisPerPatch};
497
498 const int frontLeftCellSerialised = peano4::utils::dLinearised(
499 frontLeftCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
500 );
501 const int frontLeftCellN1Serialised = peano4::utils::dLinearised(
502 frontLeftCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
503 );
504 const int frontLeftCellN2Serialised = peano4::utils::dLinearised(
505 frontLeftCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
506 );
507
508 tarch::la::Vector<Dimensions, int> frontRightCell = {
509 NumberOfFiniteVolumesPerAxisPerPatch + 1, i, NumberOfFiniteVolumesPerAxisPerPatch + 1
510 };
511 tarch::la::Vector<Dimensions, int> frontRightCellN1 = {
512 NumberOfFiniteVolumesPerAxisPerPatch, i, NumberOfFiniteVolumesPerAxisPerPatch + 1
513 };
514 tarch::la::Vector<Dimensions, int> frontRightCellN2 = {
515 NumberOfFiniteVolumesPerAxisPerPatch + 1, i, NumberOfFiniteVolumesPerAxisPerPatch
516 };
517
518 const int frontRightCellSerialised = peano4::utils::dLinearised(
519 frontRightCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
520 );
521 const int frontRightCellN1Serialised = peano4::utils::dLinearised(
522 frontRightCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
523 );
524 const int frontRightCellN2Serialised = peano4::utils::dLinearised(
525 frontRightCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
526 );
527
529 i, NumberOfFiniteVolumesPerAxisPerPatch + 1, NumberOfFiniteVolumesPerAxisPerPatch + 1
530 };
531 tarch::la::Vector<Dimensions, int> frontTopCellN1 = {
532 i, NumberOfFiniteVolumesPerAxisPerPatch + 1, NumberOfFiniteVolumesPerAxisPerPatch
533 };
534 tarch::la::Vector<Dimensions, int> frontTopCellN2 = {
535 i, NumberOfFiniteVolumesPerAxisPerPatch, NumberOfFiniteVolumesPerAxisPerPatch + 1
536 };
537
538 const int frontTopCellSerialised = peano4::utils::dLinearised(
539 frontTopCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
540 );
541 const int frontTopCellN1Serialised = peano4::utils::dLinearised(
542 frontTopCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
543 );
544 const int frontTopCellN2Serialised = peano4::utils::dLinearised(
545 frontTopCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
546 );
547
548 tarch::la::Vector<Dimensions, int> frontBottomCell = {i, 0, NumberOfFiniteVolumesPerAxisPerPatch + 1};
549 tarch::la::Vector<Dimensions, int> frontBottomCellN1 = {i, 0, NumberOfFiniteVolumesPerAxisPerPatch};
550 tarch::la::Vector<Dimensions, int> frontBottomCellN2 = {i, 1, NumberOfFiniteVolumesPerAxisPerPatch + 1};
551
552 const int frontBottomCellSerialised = peano4::utils::dLinearised(
553 frontBottomCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
554 );
555 const int frontBottomCellN1Serialised = peano4::utils::dLinearised(
556 frontBottomCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
557 );
558 const int frontBottomCellN2Serialised = peano4::utils::dLinearised(
559 frontBottomCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
560 );
561
562 tarch::la::Vector<Dimensions, int> bottomLeftCell = {0, 0, i};
563 tarch::la::Vector<Dimensions, int> bottomLeftCellN1 = {1, 0, i};
564 tarch::la::Vector<Dimensions, int> bottomLeftCellN2 = {0, 1, i};
565
566 const int bottomLeftCellSerialised = peano4::utils::dLinearised(
567 bottomLeftCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
568 );
569 const int bottomLeftCellN1Serialised = peano4::utils::dLinearised(
570 bottomLeftCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
571 );
572 const int bottomLeftCellN2Serialised = peano4::utils::dLinearised(
573 bottomLeftCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
574 );
575
576 tarch::la::Vector<Dimensions, int> bottomRightCell = {NumberOfFiniteVolumesPerAxisPerPatch + 1, 0, i};
577 tarch::la::Vector<Dimensions, int> bottomRightCellN1 = {NumberOfFiniteVolumesPerAxisPerPatch, 0, i};
578 tarch::la::Vector<Dimensions, int> bottomRightCellN2 = {NumberOfFiniteVolumesPerAxisPerPatch + 1, 1, i};
579
580 const int bottomRightCellSerialised = peano4::utils::dLinearised(
581 bottomRightCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
582 );
583 const int bottomRightCellN1Serialised = peano4::utils::dLinearised(
584 bottomRightCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
585 );
586 const int bottomRightCellN2Serialised = peano4::utils::dLinearised(
587 bottomRightCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
588 );
589
590 tarch::la::Vector<Dimensions, int> topLeftCell = {0, NumberOfFiniteVolumesPerAxisPerPatch + 1, i};
591 tarch::la::Vector<Dimensions, int> topLeftCellN1 = {1, NumberOfFiniteVolumesPerAxisPerPatch + 1, i};
592 tarch::la::Vector<Dimensions, int> topLeftCellN2 = {0, NumberOfFiniteVolumesPerAxisPerPatch, i};
593
594 const int topLeftCellSerialised = peano4::utils::dLinearised(
595 topLeftCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
596 );
597 const int topLeftCellN1Serialised = peano4::utils::dLinearised(
598 topLeftCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
599 );
600 const int topLeftCellN2Serialised = peano4::utils::dLinearised(
601 topLeftCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
602 );
603
605 NumberOfFiniteVolumesPerAxisPerPatch + 1, NumberOfFiniteVolumesPerAxisPerPatch + 1, i
606 };
607 tarch::la::Vector<Dimensions, int> topRightCellN1 = {
608 NumberOfFiniteVolumesPerAxisPerPatch, NumberOfFiniteVolumesPerAxisPerPatch + 1, i
609 };
610 tarch::la::Vector<Dimensions, int> topRightCellN2 = {
611 NumberOfFiniteVolumesPerAxisPerPatch + 1, NumberOfFiniteVolumesPerAxisPerPatch, i
612 };
613
614 const int topRightCellSerialised = peano4::utils::dLinearised(
615 topRightCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
616 );
617 const int topRightCellN1Serialised = peano4::utils::dLinearised(
618 topRightCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
619 );
620 const int topRightCellN2Serialised = peano4::utils::dLinearised(
621 topRightCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
622 );
623
624 for (int j = 0; j < NumberOfAuxiliaryVariables + NumberOfUnknowns; j++) {
625 Q[(backLeftCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]
626 = 0.5
627 * (Q[(backLeftCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j] + Q[(backLeftCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]);
628
629 Q[(backRightCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]
630 = 0.5
631 * (Q[(backRightCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j] + Q[(backRightCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]);
632
633 Q[(backTopCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]
634 = 0.5
635 * (Q[(backTopCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j] + Q[(backTopCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]);
636
637 Q[(backBottomCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]
638 = 0.5
639 * (Q[(backBottomCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j] + Q[(backBottomCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]);
640
641 Q[(frontLeftCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]
642 = 0.5
643 * (Q[(frontLeftCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j] + Q[(frontLeftCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]);
644
645 Q[(frontRightCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]
646 = 0.5
647 * (Q[(frontRightCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j] + Q[(frontRightCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]);
648
649 Q[(frontTopCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]
650 = 0.5
651 * (Q[(frontTopCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j] + Q[(frontTopCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]);
652
653 Q[(frontBottomCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]
654 = 0.5
655 * (Q[(frontBottomCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j] + Q[(frontBottomCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]);
656
657 Q[(bottomLeftCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]
658 = 0.5
659 * (Q[(bottomLeftCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j] + Q[(bottomLeftCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]);
660
661 Q[(bottomRightCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]
662 = 0.5
663 * (Q[(bottomRightCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j] + Q[(bottomRightCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]);
664
665 Q[(topLeftCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]
666 = 0.5
667 * (Q[(topLeftCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j] + Q[(topLeftCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]);
668
669 Q[(topRightCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]
670 = 0.5
671 * (Q[(topRightCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j] + Q[(topRightCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + j]);
672 }
673 }
674
675 // Corners
676 tarch::la::Vector<Dimensions, int> topLeftBackCell = {0, NumberOfFiniteVolumesPerAxisPerPatch + 1, 0};
677 tarch::la::Vector<Dimensions, int> topRightBackCell = {
678 NumberOfFiniteVolumesPerAxisPerPatch + 1, NumberOfFiniteVolumesPerAxisPerPatch + 1, 0
679 };
680 tarch::la::Vector<Dimensions, int> bottomLeftBackCell = {0, 0, 0};
681 tarch::la::Vector<Dimensions, int> bottomRightBackCell = {NumberOfFiniteVolumesPerAxisPerPatch + 1, 0, 0};
682
683 tarch::la::Vector<Dimensions, int> topLeftFrontCell = {
684 0, NumberOfFiniteVolumesPerAxisPerPatch + 1, NumberOfFiniteVolumesPerAxisPerPatch + 1
685 };
686 tarch::la::Vector<Dimensions, int> topRightFrontCell = {
687 NumberOfFiniteVolumesPerAxisPerPatch + 1,
688 NumberOfFiniteVolumesPerAxisPerPatch + 1,
689 NumberOfFiniteVolumesPerAxisPerPatch + 1
690 };
691 tarch::la::Vector<Dimensions, int> bottomLeftFrontCell = {0, 0, NumberOfFiniteVolumesPerAxisPerPatch + 1};
692 tarch::la::Vector<Dimensions, int> bottomRightFrontCell = {
693 NumberOfFiniteVolumesPerAxisPerPatch + 1, 0, NumberOfFiniteVolumesPerAxisPerPatch + 1
694 };
695
696 tarch::la::Vector<Dimensions, int> topLeftBackCellN1 = topLeftBackCell;
697 topLeftBackCellN1(0)++;
698 tarch::la::Vector<Dimensions, int> topLeftBackCellN2 = topLeftBackCell;
699 topLeftBackCellN2(1)--;
700 tarch::la::Vector<Dimensions, int> topLeftBackCellN3 = topLeftBackCell;
701 topLeftBackCellN3(2)++;
702
703 tarch::la::Vector<Dimensions, int> topRightBackCellN1 = topRightBackCell;
704 topRightBackCellN1(0)--;
705 tarch::la::Vector<Dimensions, int> topRightBackCellN2 = topRightBackCell;
706 topRightBackCellN2(1)--;
707 tarch::la::Vector<Dimensions, int> topRightBackCellN3 = topRightBackCell;
708 topRightBackCellN3(2)++;
709
710 tarch::la::Vector<Dimensions, int> bottomLeftBackCellN1 = bottomLeftBackCell;
711 bottomLeftBackCellN1(0)++;
712 tarch::la::Vector<Dimensions, int> bottomLeftBackCellN2 = bottomLeftBackCell;
713 bottomLeftBackCellN2(1)++;
714 tarch::la::Vector<Dimensions, int> bottomLeftBackCellN3 = bottomLeftBackCell;
715 bottomLeftBackCellN3(2)++;
716
717 tarch::la::Vector<Dimensions, int> bottomRightBackCellN1 = bottomRightBackCell;
718 bottomRightBackCellN1(0)--;
719 tarch::la::Vector<Dimensions, int> bottomRightBackCellN2 = bottomRightBackCell;
720 bottomRightBackCellN2(1)++;
721 tarch::la::Vector<Dimensions, int> bottomRightBackCellN3 = bottomRightBackCell;
722 bottomRightBackCellN3(2)++;
723
724 tarch::la::Vector<Dimensions, int> topLeftFrontCellN1 = topLeftFrontCell;
725 topLeftFrontCellN1(0)++;
726 tarch::la::Vector<Dimensions, int> topLeftFrontCellN2 = topLeftFrontCell;
727 topLeftFrontCellN2(1)--;
728 tarch::la::Vector<Dimensions, int> topLeftFrontCellN3 = topLeftFrontCell;
729 topLeftFrontCellN3(2)--;
730
731 tarch::la::Vector<Dimensions, int> topRightFrontCellN1 = topRightFrontCell;
732 topRightFrontCellN1(0)--;
733 tarch::la::Vector<Dimensions, int> topRightFrontCellN2 = topRightFrontCell;
734 topRightFrontCellN2(1)--;
735 tarch::la::Vector<Dimensions, int> topRightFrontCellN3 = topRightFrontCell;
736 topRightFrontCellN3(2)--;
737
738 tarch::la::Vector<Dimensions, int> bottomLeftFrontCellN1 = {1, 0, NumberOfFiniteVolumesPerAxisPerPatch + 1};
739 tarch::la::Vector<Dimensions, int> bottomLeftFrontCellN2 = {0, 1, NumberOfFiniteVolumesPerAxisPerPatch + 1};
740 tarch::la::Vector<Dimensions, int> bottomLeftFrontCellN3 = {0, 0, NumberOfFiniteVolumesPerAxisPerPatch};
741
742 tarch::la::Vector<Dimensions, int> bottomRightFrontCellN1 = bottomRightFrontCell;
743 bottomRightFrontCellN1(0)--;
744 tarch::la::Vector<Dimensions, int> bottomRightFrontCellN2 = bottomRightFrontCell;
745 bottomRightFrontCellN2(1)++;
746 tarch::la::Vector<Dimensions, int> bottomRightFrontCellN3 = bottomRightFrontCell;
747 bottomRightFrontCellN3(2)--;
748
749 const int topLeftBackCellSerialised = peano4::utils::dLinearised(
750 topLeftBackCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
751 );
752 const int topRightBackCellSerialised = peano4::utils::dLinearised(
753 topRightBackCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
754 );
755 const int bottomLeftBackCellSerialised = peano4::utils::dLinearised(
756 bottomLeftBackCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
757 );
758 const int bottomRightBackCellSerialised = peano4::utils::dLinearised(
759 bottomRightBackCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
760 );
761
762 const int topLeftFrontCellSerialised = peano4::utils::dLinearised(
763 topLeftFrontCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
764 );
765 const int topRightFrontCellSerialised = peano4::utils::dLinearised(
766 topRightFrontCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
767 );
768 const int bottomLeftFrontCellSerialised = peano4::utils::dLinearised(
769 bottomLeftFrontCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
770 );
771 const int bottomRightFrontCellSerialised = peano4::utils::dLinearised(
772 bottomRightFrontCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
773 );
774
775 const int topLeftBackCellN1Serialised = peano4::utils::dLinearised(
776 topLeftBackCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
777 );
778 const int topLeftBackCellN2Serialised = peano4::utils::dLinearised(
779 topLeftBackCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
780 );
781 const int topLeftBackCellN3Serialised = peano4::utils::dLinearised(
782 topLeftBackCellN3, NumberOfFiniteVolumesPerAxisPerPatch + 2
783 );
784
785 const int topRightBackCellN1Serialised = peano4::utils::dLinearised(
786 topRightBackCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
787 );
788 const int topRightBackCellN2Serialised = peano4::utils::dLinearised(
789 topRightBackCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
790 );
791 const int topRightBackCellN3Serialised = peano4::utils::dLinearised(
792 topRightBackCellN3, NumberOfFiniteVolumesPerAxisPerPatch + 2
793 );
794
795 const int bottomLeftBackCellN1Serialised = peano4::utils::dLinearised(
796 bottomLeftBackCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
797 );
798 const int bottomLeftBackCellN2Serialised = peano4::utils::dLinearised(
799 bottomLeftBackCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
800 );
801 const int bottomLeftBackCellN3Serialised = peano4::utils::dLinearised(
802 bottomLeftBackCellN3, NumberOfFiniteVolumesPerAxisPerPatch + 2
803 );
804
805 const int bottomRightBackCellN1Serialised = peano4::utils::dLinearised(
806 bottomRightBackCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
807 );
808 const int bottomRightBackCellN2Serialised = peano4::utils::dLinearised(
809 bottomRightBackCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
810 );
811 const int bottomRightBackCellN3Serialised = peano4::utils::dLinearised(
812 bottomRightBackCellN3, NumberOfFiniteVolumesPerAxisPerPatch + 2
813 );
814
815 const int topLeftFrontCellN1Serialised = peano4::utils::dLinearised(
816 topLeftFrontCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
817 );
818 const int topLeftFrontCellN2Serialised = peano4::utils::dLinearised(
819 topLeftFrontCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
820 );
821 const int topLeftFrontCellN3Serialised = peano4::utils::dLinearised(
822 topLeftFrontCellN3, NumberOfFiniteVolumesPerAxisPerPatch + 2
823 );
824
825 const int topRightFrontCellN1Serialised = peano4::utils::dLinearised(
826 topRightFrontCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
827 );
828 const int topRightFrontCellN2Serialised = peano4::utils::dLinearised(
829 topRightFrontCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
830 );
831 const int topRightFrontCellN3Serialised = peano4::utils::dLinearised(
832 topRightFrontCellN3, NumberOfFiniteVolumesPerAxisPerPatch + 2
833 );
834
835 const int bottomLeftFrontCellN1Serialised = peano4::utils::dLinearised(
836 bottomLeftFrontCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
837 );
838 const int bottomLeftFrontCellN2Serialised = peano4::utils::dLinearised(
839 bottomLeftFrontCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
840 );
841 const int bottomLeftFrontCellN3Serialised = peano4::utils::dLinearised(
842 bottomLeftFrontCellN3, NumberOfFiniteVolumesPerAxisPerPatch + 2
843 );
844
845 const int bottomRightFrontCellN1Serialised = peano4::utils::dLinearised(
846 bottomRightFrontCellN1, NumberOfFiniteVolumesPerAxisPerPatch + 2
847 );
848 const int bottomRightFrontCellN2Serialised = peano4::utils::dLinearised(
849 bottomRightFrontCellN2, NumberOfFiniteVolumesPerAxisPerPatch + 2
850 );
851 const int bottomRightFrontCellN3Serialised = peano4::utils::dLinearised(
852 bottomRightFrontCellN3, NumberOfFiniteVolumesPerAxisPerPatch + 2
853 );
854
855 for (int i = 0; i < NumberOfUnknowns + NumberOfAuxiliaryVariables; i++) {
856
857 Q[(topLeftBackCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]
858 = 1.0 / 3.0
859 * (Q[(topLeftBackCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(topLeftBackCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(topLeftBackCellN3Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]);
860
861 Q[(topRightBackCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]
862 = 1.0 / 3.0
863 * (Q[(topRightBackCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(topRightBackCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(topRightBackCellN3Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]);
864
865 Q[(bottomLeftBackCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]
866 = 1.0 / 3.0
867 * (Q[(bottomLeftBackCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(bottomLeftBackCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(bottomLeftBackCellN3Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]);
868
869 Q[(bottomRightBackCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]
870 = 1.0 / 3.0
871 * (Q[(bottomRightBackCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(bottomRightBackCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(bottomRightBackCellN3Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]);
872
873 Q[(topLeftFrontCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]
874 = 1.0 / 3.0
875 * (Q[(topLeftFrontCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(topLeftFrontCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(topLeftFrontCellN3Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]);
876
877 Q[(topRightFrontCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]
878 = 1.0 / 3.0
879 * (Q[(topRightFrontCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(topRightFrontCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(topRightFrontCellN3Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]);
880
881 Q[(bottomLeftFrontCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]
882 = 1.0 / 3.0
883 * (Q[(bottomLeftFrontCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(bottomLeftFrontCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(bottomLeftFrontCellN3Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]);
884
885 Q[(bottomRightFrontCellSerialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]
886 = 0.5
887 * (Q[(bottomRightFrontCellN1Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(bottomRightFrontCellN2Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i] + Q[(bottomRightFrontCellN3Serialised) * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + i]);
888 }
889#endif
890 }
891
892 static inline void calculateDerivatives(double* __restrict__ Q) {
893
894#if Dimensions == 2
895
896 auto pressure = [](double* __restrict__ Q) {
897 return (Gamma - 1) * (Q[e] - 0.5 * (Q[u] * Q[u] + Q[v] * Q[v]) / Q[rho] - Q[Z]);
898 };
899
900 for (int x = 0; x < NumberOfFiniteVolumesPerAxisPerPatch + 1; x++)
901 for (int y = 0; y < NumberOfFiniteVolumesPerAxisPerPatch + 1; y++) {
902 const tarch::la::Vector<Dimensions, int> currentCell = {x, y};
903 const tarch::la::Vector<Dimensions, int> xNeighbourCell = {x + 1, y};
904 const tarch::la::Vector<Dimensions, int> yNeighbourCell = {x, y + 1};
905
906 const int currentCellSerialized = peano4::utils::dLinearised(
907 currentCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
908 );
909 const int xNeighbourCellSerialized = peano4::utils::dLinearised(
910 xNeighbourCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
911 );
912 const int yNeighbourCellSerialized = peano4::utils::dLinearised(
913 yNeighbourCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
914 );
915
916 // u_x
917 Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + NumberOfUnknowns]
918 = -Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + u]
919 + Q[xNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + u];
920 // u_y
921 Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + NumberOfUnknowns + 1]
922 = -Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 1]
923 + Q[yNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 1];
924 // v_x
925 Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + NumberOfUnknowns + 2]
926 = -Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 2]
927 + Q[xNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 2];
928 // v_y
929 Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + NumberOfUnknowns + 3]
930 = -Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 2]
931 + Q[yNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 2];
932 // T_x
933 Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + NumberOfUnknowns + 4]
934 = -pressure(&Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)]
935 ) / (R * Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)])
936 + pressure(&Q[xNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)]
937 ) / (R * Q[xNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)]);
938 // T_y
939 Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + NumberOfUnknowns + 4]
940 = -pressure(&Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)]
941 ) / (R * Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)])
942 + pressure(&Q[yNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)]
943 ) / (R * Q[yNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)]);
944 }
945#else
946
947 double pressure = {[](double* __restrict__ Q) {
948 return (Gamma - 1) * (Q[e] - 0.5 * (Q[u] * Q[u] + Q[v] * Q[v] + Q[w] * Q[w]) / Q[rho] - Q[Z]);
949 }};
950
951 for (int x = 0; x < NumberOfFiniteVolumesPerAxisPerPatch + 1; x++)
952 for (int y = 0; y < NumberOfFiniteVolumesPerAxisPerPatch + 1; y++)
953 for (int z = 0; z < NumberOfFiniteVolumesPerAxisPerPatch + 1; z++) {
954 const tarch::la::Vector<Dimensions, int> currentCell = {x, y, z};
955 const tarch::la::Vector<Dimensions, int> xNeighbourCell = {x + 1, y, z};
956 const tarch::la::Vector<Dimensions, int> yNeighbourCell = {x, y + 1, z};
957 const tarch::la::Vector<Dimensions, int> zNeighbourCell = {x, y, z + 1};
958
959 const int currentCellSerialized = peano4::utils::dLinearised(
960 currentCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
961 );
962 const int xNeighbourCellSerialized = peano4::utils::dLinearised(
963 xNeighbourCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
964 );
965 const int yNeighbourCellSerialized = peano4::utils::dLinearised(
966 yNeighbourCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
967 );
968 const int zNeighbourCellSerialized = peano4::utils::dLinearised(
969 zNeighbourCell, NumberOfFiniteVolumesPerAxisPerPatch + 2
970 );
971
972 // u_x
973 Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + NumberOfUnknowns]
974 = -Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 1]
975 + Q[xNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 1];
976 // u_y
977 Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + NumberOfUnknowns + 1]
978 = -Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 1]
979 + Q[yNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 1];
980 // u_z
981 Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + NumberOfUnknowns + 2]
982 = -Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 1]
983 + Q[zNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 1];
984 // v_x
985 Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + NumberOfUnknowns + 3]
986 = -Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 2]
987 + Q[xNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 2];
988 // v_y
989 Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + NumberOfUnknowns + 4]
990 = -Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 2]
991 + Q[yNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 2];
992 // v_z
993 Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + NumberOfUnknowns + 5]
994 = -Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 2]
995 + Q[zNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 2];
996 // w_x
997 Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + NumberOfUnknowns + 6]
998 = -Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 3]
999 + Q[xNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 3];
1000 // w_y
1001 Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + NumberOfUnknowns + 7]
1002 = -Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 3]
1003 + Q[yNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 3];
1004 // w_x
1005 Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + NumberOfUnknowns + 8]
1006 = -Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 3]
1007 + Q[zNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + 3];
1008 // T_x
1009 Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + NumberOfUnknowns + 9]
1010 = -pressure(&Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)]
1011 ) / (R * Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)])
1012 + pressure(&Q[xNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)]
1013 ) / (R * Q[xNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)]);
1014 // T_y
1015 Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + NumberOfUnknowns + 10]
1016 = -pressure(&Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)]
1017 ) / (R * Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)])
1018 + pressure(&Q[yNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)]
1019 ) / (R * Q[yNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)]);
1020 // T_z
1021 Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables) + NumberOfUnknowns + 11]
1022 = -pressure(&Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)]
1023 ) / (R * Q[currentCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)])
1024 + pressure(&Q[zNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)]
1025 ) / (R * Q[zNeighbourCellSerialized * (NumberOfUnknowns + NumberOfAuxiliaryVariables)]);
1026 }
1027#endif
1028 }
1029};
#define assertion3(expr, param0, param1, param2)
#define assertion(expr)
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
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
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, double dt, int normal, double *__restrict__ BTimesDeltaQ, Offloadable)
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
static 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, Offloadable)
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
static 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, Offloadable)
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
CPUGPUMethod int dLinearised(const tarch::la::Vector< Dimensions, int > &counter, int max)
Map d-dimensional vector onto integer index.
Definition Loop.cpp:106
Simple vector class.
Definition Vector.h:134