Peano 4
Loading...
Searching...
No Matches
PatchwiseInsitu.cpph
Go to the documentation of this file.
1// This file is part of the ExaHyPE2 project. For conditions of distribution and
2// use, please see the copyright notice at www.peano-framework.org
3
4template <
5 class SolverType,
6 int NumberOfVolumesPerAxisInPatch,
7 int HaloSize,
8 int NumberOfUnknowns,
9 int NumberOfAuxiliaryVariables,
10 bool EvaluateFlux,
11 bool EvaluateNonconservativeProduct,
12 bool EvaluateSource,
13 bool EvaluateMaximumEigenvalueAfterTimeStep>
14void exahype2::fv::rusanov::timeStepWithRusanovPatchwiseInsituStateless(CellData& patchData, peano4::utils::LoopPlacement loopParallelism) {
15 static_assert(HaloSize == 1);
16
17 static tarch::logging::Log _log("exahype2::fv::rusanov");
18 logTraceIn("timeStepWithRusanovPatchwiseInsituStateless()");
19
20 const enumerator::AoSLexicographicEnumerator QInEnumerator(1, NumberOfVolumesPerAxisInPatch, HaloSize, NumberOfUnknowns, NumberOfAuxiliaryVariables);
21 const enumerator::AoSLexicographicEnumerator QOutEnumerator(1, NumberOfVolumesPerAxisInPatch, 0, NumberOfUnknowns, NumberOfAuxiliaryVariables);
22 const enumerator::SingleDoFEnumerator singleVolumeEnumerator(NumberOfUnknowns, 0);
23
24 for (int patchIndex = 0; patchIndex < patchData.numberOfCells; patchIndex++) {
25 // ====================================================
26 // Copy solution over and evaluate source (if required)
27 // ====================================================
28 if constexpr (EvaluateSource) {
29#if Dimensions == 2
30#pragma omp simd collapse(2)
31 for (int y = 0; y < NumberOfVolumesPerAxisInPatch; y++) {
32 for (int x = 0; x < NumberOfVolumesPerAxisInPatch; x++) {
33 loopbodies::copySolutionAndAddSourceTerm<SolverType>(
34 patchData.QIn[patchIndex],
35 QInEnumerator,
36 patchData.cellCentre[patchIndex],
37 patchData.cellSize[patchIndex],
38 patchIndex,
39 volumeIndex(x, y),
40 patchData.t[patchIndex],
41 patchData.dt[patchIndex],
42 patchData.QOut[patchIndex],
43 QOutEnumerator
44 );
45 }
46 }
47#else
48#pragma omp simd collapse(3)
49 for (int z = 0; z < NumberOfVolumesPerAxisInPatch; z++) {
50 for (int y = 0; y < NumberOfVolumesPerAxisInPatch; y++) {
51 for (int x = 0; x < NumberOfVolumesPerAxisInPatch; x++) {
52 loopbodies::copySolutionAndAddSourceTerm<SolverType>(
53 patchData.QIn[patchIndex],
54 QInEnumerator,
55 patchData.cellCentre[patchIndex],
56 patchData.cellSize[patchIndex],
57 patchIndex,
58 volumeIndex(x, y, z),
59 patchData.t[patchIndex],
60 patchData.dt[patchIndex],
61 patchData.QOut[patchIndex],
62 QOutEnumerator
63 );
64 }
65 }
66 }
67#endif
68 } else {
69#if Dimensions == 2
70#pragma omp simd collapse(3)
71 for (int y = 0; y < NumberOfVolumesPerAxisInPatch; y++) {
72 for (int x = 0; x < NumberOfVolumesPerAxisInPatch; x++) {
73 for (int unknown = 0; unknown < NumberOfUnknowns + NumberOfAuxiliaryVariables; unknown++) {
74 loopbodies::copySolution(patchData.QIn[patchIndex], QInEnumerator, patchIndex, volumeIndex(x, y), unknown, patchData.QOut[patchIndex], QOutEnumerator);
75 }
76 }
77 }
78#else
79#pragma omp simd collapse(4)
80 for (int z = 0; z < NumberOfVolumesPerAxisInPatch; z++) {
81 for (int y = 0; y < NumberOfVolumesPerAxisInPatch; y++) {
82 for (int x = 0; x < NumberOfVolumesPerAxisInPatch; x++) {
83 for (int unknown = 0; unknown < NumberOfUnknowns + NumberOfAuxiliaryVariables; unknown++) {
84 loopbodies::copySolution(patchData.QIn[patchIndex], QInEnumerator, patchIndex, volumeIndex(x, y, z), unknown, patchData.QOut[patchIndex], QOutEnumerator);
85 }
86 }
87 }
88 }
89#endif
90 }
91
92 // ====================================================
93 // Compute damping due to max eigenvalue
94 // ====================================================
95#if Dimensions == 2
96 for (int shift = 0; shift < 2; shift++) {
97 for (int x = -HaloSize + shift; x < NumberOfVolumesPerAxisInPatch + HaloSize - 1; x += 2) {
98#pragma omp simd
99 for (int y = 0; y < NumberOfVolumesPerAxisInPatch; y++) {
100 double maxEigenvalueL __attribute__((aligned(AlignmentOnHeap)));
101 double maxEigenvalueR __attribute__((aligned(AlignmentOnHeap)));
102 loopbodies::computeMaxEigenvalue<SolverType>(
103 patchData.QIn[patchIndex],
104 QInEnumerator,
105 patchData.cellCentre[patchIndex],
106 patchData.cellSize[patchIndex],
107 patchIndex,
108 volumeIndex(x, y),
109 patchData.t[patchIndex],
110 patchData.dt[patchIndex],
111 0, // normal
112 &maxEigenvalueL,
113 singleVolumeEnumerator
114 );
115 loopbodies::computeMaxEigenvalue<SolverType>(
116 patchData.QIn[patchIndex],
117 QInEnumerator,
118 patchData.cellCentre[patchIndex],
119 patchData.cellSize[patchIndex],
120 patchIndex,
121 volumeIndex(x + 1, y),
122 patchData.t[patchIndex],
123 patchData.dt[patchIndex],
124 0, // normal
125 &maxEigenvalueR,
126 singleVolumeEnumerator
127 );
128 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
129 loopbodies::updateSolutionWithEigenvalueDamping(
130 patchData.QIn[patchIndex],
131 QInEnumerator,
132 maxEigenvalueL,
133 maxEigenvalueR,
134 patchData.cellCentre[patchIndex],
135 patchData.cellSize[patchIndex],
136 patchIndex,
137 volumeIndex(x, y),
138 unknown,
139 patchData.dt[patchIndex],
140 0, // normal
141 patchData.QOut[patchIndex],
142 QOutEnumerator
143 );
144 }
145 }
146 }
147 }
148
149 for (int shift = 0; shift < 2; shift++) {
150 for (int y = -HaloSize + shift; y < NumberOfVolumesPerAxisInPatch + HaloSize - 1; y += 2) {
151#pragma omp simd
152 for (int x = 0; x < NumberOfVolumesPerAxisInPatch; x++) {
153 double maxEigenvalueL __attribute__((aligned(AlignmentOnHeap)));
154 double maxEigenvalueR __attribute__((aligned(AlignmentOnHeap)));
155 loopbodies::computeMaxEigenvalue<SolverType>(
156 patchData.QIn[patchIndex],
157 QInEnumerator,
158 patchData.cellCentre[patchIndex],
159 patchData.cellSize[patchIndex],
160 patchIndex,
161 volumeIndex(x, y),
162 patchData.t[patchIndex],
163 patchData.dt[patchIndex],
164 1, // normal
165 &maxEigenvalueL,
166 singleVolumeEnumerator
167 );
168 loopbodies::computeMaxEigenvalue<SolverType>(
169 patchData.QIn[patchIndex],
170 QInEnumerator,
171 patchData.cellCentre[patchIndex],
172 patchData.cellSize[patchIndex],
173 patchIndex,
174 volumeIndex(x, y + 1),
175 patchData.t[patchIndex],
176 patchData.dt[patchIndex],
177 1, // normal
178 &maxEigenvalueR,
179 singleVolumeEnumerator
180 );
181 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
182 loopbodies::updateSolutionWithEigenvalueDamping(
183 patchData.QIn[patchIndex],
184 QInEnumerator,
185 maxEigenvalueL,
186 maxEigenvalueR,
187 patchData.cellCentre[patchIndex],
188 patchData.cellSize[patchIndex],
189 patchIndex,
190 volumeIndex(x, y),
191 unknown,
192 patchData.dt[patchIndex],
193 1, // normal
194 patchData.QOut[patchIndex],
195 QOutEnumerator
196 );
197 }
198 }
199 }
200 }
201#else
202 for (int shift = 0; shift < 2; shift++) {
203 for (int x = -HaloSize + shift; x < NumberOfVolumesPerAxisInPatch + HaloSize - 1; x += 2) {
204#pragma omp simd collapse(2)
205 for (int y = 0; y < NumberOfVolumesPerAxisInPatch; y++) {
206 for (int z = 0; z < NumberOfVolumesPerAxisInPatch; z++) {
207 double maxEigenvalueL __attribute__((aligned(AlignmentOnHeap)));
208 double maxEigenvalueR __attribute__((aligned(AlignmentOnHeap)));
209 loopbodies::computeMaxEigenvalue<SolverType>(
210 patchData.QIn[patchIndex],
211 QInEnumerator,
212 patchData.cellCentre[patchIndex],
213 patchData.cellSize[patchIndex],
214 patchIndex,
215 volumeIndex(x, y, z),
216 patchData.t[patchIndex],
217 patchData.dt[patchIndex],
218 0, // normal
219 &maxEigenvalueL,
220 singleVolumeEnumerator
221 );
222 loopbodies::computeMaxEigenvalue<SolverType>(
223 patchData.QIn[patchIndex],
224 QInEnumerator,
225 patchData.cellCentre[patchIndex],
226 patchData.cellSize[patchIndex],
227 patchIndex,
228 volumeIndex(x + 1, y, z),
229 patchData.t[patchIndex],
230 patchData.dt[patchIndex],
231 0, // normal
232 &maxEigenvalueR,
233 singleVolumeEnumerator
234 );
235 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
236 loopbodies::updateSolutionWithEigenvalueDamping(
237 patchData.QIn[patchIndex],
238 QInEnumerator,
239 maxEigenvalueL,
240 maxEigenvalueR,
241 patchData.cellCentre[patchIndex],
242 patchData.cellSize[patchIndex],
243 patchIndex,
244 volumeIndex(x, y, z),
245 unknown,
246 patchData.dt[patchIndex],
247 0, // normal
248 patchData.QOut[patchIndex],
249 QOutEnumerator
250 );
251 }
252 }
253 }
254 }
255 }
256
257 for (int shift = 0; shift < 2; shift++) {
258 for (int y = -HaloSize + shift; y < NumberOfVolumesPerAxisInPatch + HaloSize - 1; y += 2) {
259#pragma omp simd collapse(2)
260 for (int x = 0; x < NumberOfVolumesPerAxisInPatch; x++) {
261 for (int z = 0; z < NumberOfVolumesPerAxisInPatch; z++) {
262 double maxEigenvalueL __attribute__((aligned(AlignmentOnHeap)));
263 double maxEigenvalueR __attribute__((aligned(AlignmentOnHeap)));
264 loopbodies::computeMaxEigenvalue<SolverType>(
265 patchData.QIn[patchIndex],
266 QInEnumerator,
267 patchData.cellCentre[patchIndex],
268 patchData.cellSize[patchIndex],
269 patchIndex,
270 volumeIndex(x, y, z),
271 patchData.t[patchIndex],
272 patchData.dt[patchIndex],
273 1, // normal
274 &maxEigenvalueL,
275 singleVolumeEnumerator
276 );
277 loopbodies::computeMaxEigenvalue<SolverType>(
278 patchData.QIn[patchIndex],
279 QInEnumerator,
280 patchData.cellCentre[patchIndex],
281 patchData.cellSize[patchIndex],
282 patchIndex,
283 volumeIndex(x, y + 1, z),
284 patchData.t[patchIndex],
285 patchData.dt[patchIndex],
286 1, // normal
287 &maxEigenvalueR,
288 singleVolumeEnumerator
289 );
290 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
291 loopbodies::updateSolutionWithEigenvalueDamping(
292 patchData.QIn[patchIndex],
293 QInEnumerator,
294 maxEigenvalueL,
295 maxEigenvalueR,
296 patchData.cellCentre[patchIndex],
297 patchData.cellSize[patchIndex],
298 patchIndex,
299 volumeIndex(x, y, z),
300 unknown,
301 patchData.dt[patchIndex],
302 1, // normal
303 patchData.QOut[patchIndex],
304 QOutEnumerator
305 );
306 }
307 }
308 }
309 }
310 }
311
312 for (int shift = 0; shift < 2; shift++) {
313 for (int z = -HaloSize + shift; z < NumberOfVolumesPerAxisInPatch + HaloSize - 1; z += 2) {
314#pragma omp simd collapse(2)
315 for (int x = 0; x < NumberOfVolumesPerAxisInPatch; x++) {
316 for (int y = 0; y < NumberOfVolumesPerAxisInPatch; y++) {
317 double maxEigenvalueL __attribute__((aligned(AlignmentOnHeap)));
318 double maxEigenvalueR __attribute__((aligned(AlignmentOnHeap)));
319 loopbodies::computeMaxEigenvalue<SolverType>(
320 patchData.QIn[patchIndex],
321 QInEnumerator,
322 patchData.cellCentre[patchIndex],
323 patchData.cellSize[patchIndex],
324 patchIndex,
325 volumeIndex(x, y, z),
326 patchData.t[patchIndex],
327 patchData.dt[patchIndex],
328 2, // normal
329 &maxEigenvalueL,
330 singleVolumeEnumerator
331 );
332 loopbodies::computeMaxEigenvalue<SolverType>(
333 patchData.QIn[patchIndex],
334 QInEnumerator,
335 patchData.cellCentre[patchIndex],
336 patchData.cellSize[patchIndex],
337 patchIndex,
338 volumeIndex(x, y, z + 1),
339 patchData.t[patchIndex],
340 patchData.dt[patchIndex],
341 2, // normal
342 &maxEigenvalueR,
343 singleVolumeEnumerator
344 );
345 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
346 loopbodies::updateSolutionWithEigenvalueDamping(
347 patchData.QIn[patchIndex],
348 QInEnumerator,
349 maxEigenvalueL,
350 maxEigenvalueR,
351 patchData.cellCentre[patchIndex],
352 patchData.cellSize[patchIndex],
353 patchIndex,
354 volumeIndex(x, y, z),
355 unknown,
356 patchData.dt[patchIndex],
357 2, // normal
358 patchData.QOut[patchIndex],
359 QOutEnumerator
360 );
361 }
362 }
363 }
364 }
365 }
366#endif
367
368 // ====================================================
369 // Normal (conservative) flux
370 // ====================================================
371 if constexpr (EvaluateFlux) {
372#if Dimensions == 2
373 for (int shift = 0; shift < 2; shift++) {
374 for (int x = -HaloSize + shift; x < NumberOfVolumesPerAxisInPatch + HaloSize - 1; x += 2) {
375#pragma omp simd
376 for (int y = 0; y < NumberOfVolumesPerAxisInPatch; y++) {
377 double fluxFL[NumberOfUnknowns] __attribute__((aligned(AlignmentOnHeap)));
378 double fluxFR[NumberOfUnknowns] __attribute__((aligned(AlignmentOnHeap)));
379 loopbodies::computeFlux<SolverType>(
380 patchData.QIn[patchIndex],
381 QInEnumerator,
382 patchData.cellCentre[patchIndex],
383 patchData.cellSize[patchIndex],
384 patchIndex,
385 volumeIndex(x, y),
386 patchData.t[patchIndex],
387 patchData.dt[patchIndex],
388 0, // normal
389 fluxFL,
390 singleVolumeEnumerator
391 );
392 loopbodies::computeFlux<SolverType>(
393 patchData.QIn[patchIndex],
394 QInEnumerator,
395 patchData.cellCentre[patchIndex],
396 patchData.cellSize[patchIndex],
397 patchIndex,
398 volumeIndex(x + 1, y),
399 patchData.t[patchIndex],
400 patchData.dt[patchIndex],
401 0, // normal
402 fluxFR,
403 singleVolumeEnumerator
404 );
405 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
406 loopbodies::updateSolutionWithFlux(
407 fluxFL,
408 fluxFR,
409 singleVolumeEnumerator,
410 patchData.cellCentre[patchIndex],
411 patchData.cellSize[patchIndex],
412 patchIndex,
413 volumeIndex(x, y),
414 unknown,
415 patchData.dt[patchIndex],
416 0, // normal
417 patchData.QOut[patchIndex],
418 QOutEnumerator
419 );
420 }
421 }
422 }
423 }
424
425 for (int shift = 0; shift < 2; shift++) {
426 for (int y = -HaloSize + shift; y < NumberOfVolumesPerAxisInPatch + HaloSize - 1; y += 2) {
427#pragma omp simd
428 for (int x = 0; x < NumberOfVolumesPerAxisInPatch; x++) {
429 double fluxFL[NumberOfUnknowns] __attribute__((aligned(AlignmentOnHeap)));
430 double fluxFR[NumberOfUnknowns] __attribute__((aligned(AlignmentOnHeap)));
431 loopbodies::computeFlux<SolverType>(
432 patchData.QIn[patchIndex],
433 QInEnumerator,
434 patchData.cellCentre[patchIndex],
435 patchData.cellSize[patchIndex],
436 patchIndex,
437 volumeIndex(x, y),
438 patchData.t[patchIndex],
439 patchData.dt[patchIndex],
440 1, // normal
441 fluxFL,
442 singleVolumeEnumerator
443 );
444 loopbodies::computeFlux<SolverType>(
445 patchData.QIn[patchIndex],
446 QInEnumerator,
447 patchData.cellCentre[patchIndex],
448 patchData.cellSize[patchIndex],
449 patchIndex,
450 volumeIndex(x, y + 1),
451 patchData.t[patchIndex],
452 patchData.dt[patchIndex],
453 1, // normal
454 fluxFR,
455 singleVolumeEnumerator
456 );
457 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
458 loopbodies::updateSolutionWithFlux(
459 fluxFL,
460 fluxFR,
461 singleVolumeEnumerator,
462 patchData.cellCentre[patchIndex],
463 patchData.cellSize[patchIndex],
464 patchIndex,
465 volumeIndex(x, y),
466 unknown,
467 patchData.dt[patchIndex],
468 1, // normal
469 patchData.QOut[patchIndex],
470 QOutEnumerator
471 );
472 }
473 }
474 }
475 }
476#else
477 for (int shift = 0; shift < 2; shift++) {
478 for (int x = HaloSize + shift; x < NumberOfVolumesPerAxisInPatch + HaloSize - 1; x += 2) {
479#pragma omp simd collapse(2)
480 for (int y = 0; y < NumberOfVolumesPerAxisInPatch; y++) {
481 for (int z = 0; z < NumberOfVolumesPerAxisInPatch; z++) {
482 double fluxFL[NumberOfUnknowns] __attribute__((aligned(AlignmentOnHeap)));
483 double fluxFR[NumberOfUnknowns] __attribute__((aligned(AlignmentOnHeap)));
484 loopbodies::computeFlux<SolverType>(
485 patchData.QIn[patchIndex],
486 QInEnumerator,
487 patchData.cellCentre[patchIndex],
488 patchData.cellSize[patchIndex],
489 patchIndex,
490 volumeIndex(x, y, z),
491 patchData.t[patchIndex],
492 patchData.dt[patchIndex],
493 0, // normal
494 fluxFL,
495 singleVolumeEnumerator
496 );
497 loopbodies::computeFlux<SolverType>(
498 patchData.QIn[patchIndex],
499 QInEnumerator,
500 patchData.cellCentre[patchIndex],
501 patchData.cellSize[patchIndex],
502 patchIndex,
503 volumeIndex(x + 1, y, z),
504 patchData.t[patchIndex],
505 patchData.dt[patchIndex],
506 0, // normal
507 fluxFR,
508 singleVolumeEnumerator
509 );
510 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
511 loopbodies::updateSolutionWithFlux(
512 fluxFL,
513 fluxFR,
514 singleVolumeEnumerator,
515 patchData.cellCentre[patchIndex],
516 patchData.cellSize[patchIndex],
517 patchIndex,
518 volumeIndex(x, y, z),
519 unknown,
520 patchData.dt[patchIndex],
521 0, // normal
522 patchData.QOut[patchIndex],
523 QOutEnumerator
524 );
525 }
526 }
527 }
528 }
529 }
530
531 for (int shift = 0; shift < 2; shift++) {
532 for (int y = -HaloSize + shift; y < NumberOfVolumesPerAxisInPatch + HaloSize - 1; y += 2) {
533#pragma omp simd collapse(2)
534 for (int x = 0; x < NumberOfVolumesPerAxisInPatch; x++) {
535 for (int z = 0; z < NumberOfVolumesPerAxisInPatch; z++) {
536 double fluxFL[NumberOfUnknowns] __attribute__((aligned(AlignmentOnHeap)));
537 double fluxFR[NumberOfUnknowns] __attribute__((aligned(AlignmentOnHeap)));
538 loopbodies::computeFlux<SolverType>(
539 patchData.QIn[patchIndex],
540 QInEnumerator,
541 patchData.cellCentre[patchIndex],
542 patchData.cellSize[patchIndex],
543 patchIndex,
544 volumeIndex(x, y, z),
545 patchData.t[patchIndex],
546 patchData.dt[patchIndex],
547 1, // normal
548 fluxFL,
549 singleVolumeEnumerator
550 );
551 loopbodies::computeFlux<SolverType>(
552 patchData.QIn[patchIndex],
553 QInEnumerator,
554 patchData.cellCentre[patchIndex],
555 patchData.cellSize[patchIndex],
556 patchIndex,
557 volumeIndex(x, y + 1, z),
558 patchData.t[patchIndex],
559 patchData.dt[patchIndex],
560 1, // normal
561 fluxFR,
562 singleVolumeEnumerator
563 );
564 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
565 loopbodies::updateSolutionWithFlux(
566 fluxFL,
567 fluxFR,
568 singleVolumeEnumerator,
569 patchData.cellCentre[patchIndex],
570 patchData.cellSize[patchIndex],
571 patchIndex,
572 volumeIndex(x, y, z),
573 unknown,
574 patchData.dt[patchIndex],
575 1, // normal
576 patchData.QOut[patchIndex],
577 QOutEnumerator
578 );
579 }
580 }
581 }
582 }
583 }
584
585 for (int shift = 0; shift < 2; shift++) {
586 for (int z = -HaloSize + shift; z < NumberOfVolumesPerAxisInPatch + HaloSize - 1; z += 2) {
587#pragma omp simd collapse(2)
588 for (int x = 0; x < NumberOfVolumesPerAxisInPatch; x++) {
589 for (int y = 0; y < NumberOfVolumesPerAxisInPatch; y++) {
590 double fluxFL[NumberOfUnknowns] __attribute__((aligned(AlignmentOnHeap)));
591 double fluxFR[NumberOfUnknowns] __attribute__((aligned(AlignmentOnHeap)));
592 loopbodies::computeFlux<SolverType>(
593 patchData.QIn[patchIndex],
594 QInEnumerator,
595 patchData.cellCentre[patchIndex],
596 patchData.cellSize[patchIndex],
597 patchIndex,
598 volumeIndex(x, y, z),
599 patchData.t[patchIndex],
600 patchData.dt[patchIndex],
601 2, // normal
602 fluxFL,
603 singleVolumeEnumerator
604 );
605 loopbodies::computeFlux<SolverType>(
606 patchData.QIn[patchIndex],
607 QInEnumerator,
608 patchData.cellCentre[patchIndex],
609 patchData.cellSize[patchIndex],
610 patchIndex,
611 volumeIndex(x, y, z + 1),
612 patchData.t[patchIndex],
613 patchData.dt[patchIndex],
614 2, // normal
615 fluxFR,
616 singleVolumeEnumerator
617 );
618 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
619 loopbodies::updateSolutionWithFlux(
620 fluxFL,
621 fluxFR,
622 singleVolumeEnumerator,
623 patchData.cellCentre[patchIndex],
624 patchData.cellSize[patchIndex],
625 patchIndex,
626 volumeIndex(x, y, z),
627 unknown,
628 patchData.dt[patchIndex],
629 2, // normal
630 patchData.QOut[patchIndex],
631 QOutEnumerator
632 );
633 }
634 }
635 }
636 }
637 }
638#endif
639 }
640
641 if constexpr (EvaluateNonconservativeProduct) {
642#if Dimensions == 2
643 for (int shift = 0; shift < 2; shift++) {
644 for (int x = -HaloSize + shift; x < NumberOfVolumesPerAxisInPatch + HaloSize - 1; x += 2) {
645#pragma omp simd
646 for (int y = 0; y < NumberOfVolumesPerAxisInPatch; y++) {
647 double fluxNCP[NumberOfUnknowns] __attribute__((aligned(AlignmentOnHeap)));
648 loopbodies::computeNonconservativeFlux<SolverType>(
649 patchData.QIn[patchIndex],
650 QInEnumerator,
651 patchData.cellCentre[patchIndex],
652 patchData.cellSize[patchIndex],
653 patchIndex,
654 volumeIndex(x, y),
655 patchData.t[patchIndex],
656 patchData.dt[patchIndex],
657 0, // normal
658 fluxNCP,
659 singleVolumeEnumerator
660 );
661 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
662 loopbodies::updateSolutionWithNonconservativeFlux(
663 fluxNCP,
664 singleVolumeEnumerator,
665 patchData.cellCentre[patchIndex],
666 patchData.cellSize[patchIndex],
667 patchIndex,
668 volumeIndex(x, y),
669 unknown,
670 patchData.dt[patchIndex],
671 0, // normal
672 patchData.QOut[patchIndex],
673 QOutEnumerator
674 );
675 }
676 }
677 }
678 }
679
680 for (int shift = 0; shift < 2; shift++) {
681 for (int y = -HaloSize + shift; y < NumberOfVolumesPerAxisInPatch + HaloSize - 1; y += 2) {
682#pragma omp simd
683 for (int x = 0; x < NumberOfVolumesPerAxisInPatch; x++) {
684 double fluxNCP[NumberOfUnknowns] __attribute__((aligned(AlignmentOnHeap)));
685 loopbodies::computeNonconservativeFlux<SolverType>(
686 patchData.QIn[patchIndex],
687 QInEnumerator,
688 patchData.cellCentre[patchIndex],
689 patchData.cellSize[patchIndex],
690 patchIndex,
691 volumeIndex(x, y),
692 patchData.t[patchIndex],
693 patchData.dt[patchIndex],
694 1, // normal
695 fluxNCP,
696 singleVolumeEnumerator
697 );
698 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
699 loopbodies::updateSolutionWithNonconservativeFlux(
700 fluxNCP,
701 singleVolumeEnumerator,
702 patchData.cellCentre[patchIndex],
703 patchData.cellSize[patchIndex],
704 patchIndex,
705 volumeIndex(x, y),
706 unknown,
707 patchData.dt[patchIndex],
708 1, // normal
709 patchData.QOut[patchIndex],
710 QOutEnumerator
711 );
712 }
713 }
714 }
715 }
716#else
717 for (int shift = 0; shift < 2; shift++) {
718 for (int x = -HaloSize + shift; x < NumberOfVolumesPerAxisInPatch + HaloSize - 1; x += 2) {
719#pragma omp simd collapse(2)
720 for (int y = 0; y < NumberOfVolumesPerAxisInPatch; y++) {
721 for (int z = 0; z < NumberOfVolumesPerAxisInPatch; z++) {
722 double fluxNCP[NumberOfUnknowns] __attribute__((aligned(AlignmentOnHeap)));
723 loopbodies::computeNonconservativeFlux<SolverType>(
724 patchData.QIn[patchIndex],
725 QInEnumerator,
726 patchData.cellCentre[patchIndex],
727 patchData.cellSize[patchIndex],
728 patchIndex,
729 volumeIndex(x, y, z),
730 patchData.t[patchIndex],
731 patchData.dt[patchIndex],
732 0, // normal
733 fluxNCP,
734 singleVolumeEnumerator
735 );
736 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
737 loopbodies::updateSolutionWithNonconservativeFlux(
738 fluxNCP,
739 singleVolumeEnumerator,
740 patchData.cellCentre[patchIndex],
741 patchData.cellSize[patchIndex],
742 patchIndex,
743 volumeIndex(x, y, z),
744 unknown,
745 patchData.dt[patchIndex],
746 0, // normal
747 patchData.QOut[patchIndex],
748 QOutEnumerator
749 );
750 }
751 }
752 }
753 }
754 }
755
756 for (int shift = 0; shift < 2; shift++) {
757 for (int y = -HaloSize + shift; y < NumberOfVolumesPerAxisInPatch + HaloSize - 1; y += 2) {
758#pragma omp simd collapse(2)
759 for (int x = 0; x < NumberOfVolumesPerAxisInPatch; x++) {
760 for (int z = 0; z < NumberOfVolumesPerAxisInPatch; z++) {
761 double fluxNCP[NumberOfUnknowns] __attribute__((aligned(AlignmentOnHeap)));
762 loopbodies::computeNonconservativeFlux<SolverType>(
763 patchData.QIn[patchIndex],
764 QInEnumerator,
765 patchData.cellCentre[patchIndex],
766 patchData.cellSize[patchIndex],
767 patchIndex,
768 volumeIndex(x, y, z),
769 patchData.t[patchIndex],
770 patchData.dt[patchIndex],
771 1, // normal
772 fluxNCP,
773 singleVolumeEnumerator
774 );
775 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
776 loopbodies::updateSolutionWithNonconservativeFlux(
777 fluxNCP,
778 singleVolumeEnumerator,
779 patchData.cellCentre[patchIndex],
780 patchData.cellSize[patchIndex],
781 patchIndex,
782 volumeIndex(x, y, z),
783 unknown,
784 patchData.dt[patchIndex],
785 1, // normal
786 patchData.QOut[patchIndex],
787 QOutEnumerator
788 );
789 }
790 }
791 }
792 }
793 }
794
795 for (int shift = 0; shift < 2; shift++) {
796 for (int z = -HaloSize + shift; z < NumberOfVolumesPerAxisInPatch + HaloSize - 1; z += 2) {
797#pragma omp simd collapse(2)
798 for (int x = 0; x < NumberOfVolumesPerAxisInPatch; x++) {
799 for (int y = 0; y < NumberOfVolumesPerAxisInPatch; y++) {
800 double fluxNCP[NumberOfUnknowns] __attribute__((aligned(AlignmentOnHeap)));
801 loopbodies::computeNonconservativeFlux<SolverType>(
802 patchData.QIn[patchIndex],
803 QInEnumerator,
804 patchData.cellCentre[patchIndex],
805 patchData.cellSize[patchIndex],
806 patchIndex,
807 volumeIndex(x, y, z),
808 patchData.t[patchIndex],
809 patchData.dt[patchIndex],
810 2, // normal
811 fluxNCP,
812 singleVolumeEnumerator
813 );
814 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
815 loopbodies::updateSolutionWithNonconservativeFlux(
816 fluxNCP,
817 singleVolumeEnumerator,
818 patchData.cellCentre[patchIndex],
819 patchData.cellSize[patchIndex],
820 patchIndex,
821 volumeIndex(x, y, z),
822 unknown,
823 patchData.dt[patchIndex],
824 2, // normal
825 patchData.QOut[patchIndex],
826 QOutEnumerator
827 );
828 }
829 }
830 }
831 }
832 }
833#endif
834 }
835
836 if constexpr (EvaluateMaximumEigenvalueAfterTimeStep) {
837 double newMaxEigenvalue = 0.0;
838#if Dimensions == 2
839#pragma omp simd collapse(2) reduction(max : newMaxEigenvalue)
840 for (int y = 0; y < NumberOfVolumesPerAxisInPatch; y++) {
841 for (int x = 0; x < NumberOfVolumesPerAxisInPatch; x++) {
842 newMaxEigenvalue = std::max(
843 newMaxEigenvalue,
844 loopbodies::reduceMaxEigenvalue<SolverType>(
845 patchData.QOut[patchIndex],
846 QOutEnumerator,
847 patchData.cellCentre[patchIndex],
848 patchData.cellSize[patchIndex],
849 patchIndex,
850 volumeIndex(x, y),
851 patchData.t[patchIndex],
852 patchData.dt[patchIndex]
853 )
854 );
855 }
856 }
857#else
858#pragma omp simd collapse(3) reduction(max : newMaxEigenvalue)
859 for (int z = 0; z < NumberOfVolumesPerAxisInPatch; z++) {
860 for (int y = 0; y < NumberOfVolumesPerAxisInPatch; y++) {
861 for (int x = 0; x < NumberOfVolumesPerAxisInPatch; x++) {
862 newMaxEigenvalue = std::max(
863 newMaxEigenvalue,
864 loopbodies::reduceMaxEigenvalue<SolverType>(
865 patchData.QOut[patchIndex],
866 QOutEnumerator,
867 patchData.cellCentre[patchIndex],
868 patchData.cellSize[patchIndex],
869 patchIndex,
870 volumeIndex(x, y, z),
871 patchData.t[patchIndex],
872 patchData.dt[patchIndex]
873 )
874 );
875 }
876 }
877 }
878#endif
879 patchData.maxEigenvalue[patchIndex] = newMaxEigenvalue;
880 }
881 }
882
883 logTraceOut("timeStepWithRusanovPatchwiseInsituStateless()");
884}
885
886
887template <
888 class SolverType,
889 int NumberOfVolumesPerAxisInPatch,
890 int HaloSize,
891 int NumberOfUnknowns,
892 int NumberOfAuxiliaryVariables,
893 bool EvaluateFlux,
894 bool EvaluateNonconservativeProduct,
895 bool EvaluateSource,
896 bool EvaluateMaximumEigenvalueAfterTimeStep>
897void exahype2::fv::rusanov::timeStepWithRusanovPatchwiseInsituStateless(
898 CellData& patchData, tarch::timing::Measurement& measurement, peano4::utils::LoopPlacement loopParallelism
899) {
900 tarch::timing::Watch watch("exahype2::fv::rusanov", "timeStepWithRusanovPatchwiseInsituStateless", false, true);
901 timeStepWithRusanovPatchwiseInsituStateless<
902 SolverType,
903 NumberOfVolumesPerAxisInPatch,
904 HaloSize,
905 NumberOfUnknowns,
906 NumberOfAuxiliaryVariables,
907 EvaluateFlux,
908 EvaluateNonconservativeProduct,
909 EvaluateSource,
910 EvaluateMaximumEigenvalueAfterTimeStep>(patchData, loopParallelism);
911 watch.stop();
912 measurement.setValue(watch.getCalendarTime());
913}
__attribute__
int __attribute__((optimize("O0"))) toolbox
Definition CompressedFloatingPointNumbers.cpp:30
HaloSize
static constexpr int HaloSize
Definition KernelBenchmarksFVRusanov-main.cpp:43
AlignmentOnHeap
#define AlignmentOnHeap
Definition LinuxAMD.h:23
logTraceOut
#define logTraceOut(methodName)
Definition Log.h:379
logTraceIn
#define logTraceIn(methodName)
Definition Log.h:369
tarch::logging::Log
Log Device.
Definition Log.h:516
tarch::timing::Measurement::setValue
void setValue(const double &value)
Set the value.
Definition Measurement.cpp:101
tarch::timing::Watch
A simple class that has to be included to measure the clock ticks required for an operation.
Definition Watch.h:45
exahype2::fv::rusanov::timeStepWithRusanovPatchwiseInsituStateless
KeywordToAvoidDuplicateSymbolsForInlinedFunctions void timeStepWithRusanovPatchwiseInsituStateless(CellData &patchData, tarch::timing::Measurement &measurement, peano4::utils::LoopPlacement loopParallelism=peano4::utils::LoopPlacement::Serial) InlineMethod
exahype2::volumeIndex
auto volumeIndex(Args... args)
Definition VolumeIndex.h:54
makeIC.z
tuple z
Definition makeIC.py:53
peano4::utils::LoopPlacement
LoopPlacement
Guide loop-level parallelism.
Definition Loop.h:60
_log
tarch::logging::Log _log("examples::unittests")