Peano
Loading...
Searching...
No Matches
LoopBodies.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 int NumberOfUnknowns,
6 int NumberOfAuxiliaryVariables,
7 class QInEnumeratorType,
8 class QOutEnumeratorType>
10 const double* __restrict__ QIn,
11 const QInEnumeratorType& QInEnumerator,
12 const SourceFunctor& sourceFunctor,
13 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
14 const ::tarch::la::Vector<Dimensions, double>& patchSize,
15 int patchIndex,
16 const ::tarch::la::Vector<Dimensions, int>& volumeIndex,
17 double t,
18 double dt,
19 double* __restrict__ QOut,
20 const QOutEnumeratorType& QOutEnumerator
21) {
22 double QInGathered[NumberOfUnknowns + NumberOfAuxiliaryVariables]{0.0};
23 double sourceGathered[NumberOfUnknowns]{0.0};
24
25 // gather
26#if defined(SharedOMP)
27#pragma omp simd
28#endif
29 for (int unknown = 0; unknown < NumberOfUnknowns + NumberOfAuxiliaryVariables;
30 unknown++) {
31 QInGathered[unknown] = QIn[QInEnumerator(patchIndex, volumeIndex, unknown)];
32 }
33
34 sourceFunctor(
35 QInGathered,
36 getVolumeCentre(
37 patchCentre,
38 patchSize,
39 QInEnumerator._numberOfDoFsPerAxisInCell,
40 volumeIndex
41 ),
42 getVolumeSize(patchSize, QInEnumerator._numberOfDoFsPerAxisInCell),
43 t,
44 dt,
45 sourceGathered
46 );
47
48 // scatter
49#if defined(SharedOMP)
50#pragma omp simd
51#endif
52 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
53 QOut[QOutEnumerator(
54 patchIndex,
55 volumeIndex,
56 unknown
57 )] = QIn[QInEnumerator(patchIndex, volumeIndex, unknown)]
58 + dt * sourceGathered[unknown];
60 QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)]
61 == QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)],
62 dt,
63 unknown,
64 volumeIndex,
65 QOutEnumerator(patchIndex, volumeIndex, unknown)
66 );
67 }
68#if defined(SharedOMP)
69#pragma omp simd
70#endif
71 for (int unknown = NumberOfUnknowns;
72 unknown < NumberOfUnknowns + NumberOfAuxiliaryVariables;
73 unknown++) {
74 QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)] = QIn
75 [QInEnumerator(patchIndex, volumeIndex, unknown)];
77 QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)]
78 == QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)],
79 dt,
80 unknown,
81 volumeIndex,
82 QOutEnumerator(patchIndex, volumeIndex, unknown)
83 );
84 }
85}
86
87
88template <int NumberOfUnknowns, int NumberOfAuxiliaryVariables>
90 const double* __restrict__ QIn,
91 const enumerator::AoSLexicographicEnumerator& QInEnumerator,
92 const SourceFunctor& sourceFunctor,
93 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
94 const ::tarch::la::Vector<Dimensions, double>& patchSize,
95 int patchIndex,
96 const ::tarch::la::Vector<Dimensions, int>& volumeIndex,
97 double t,
98 double dt,
99 double* __restrict__ QOut,
100 const enumerator::AoSLexicographicEnumerator& QOutEnumerator
101) {
102 sourceFunctor(
103 QIn + QInEnumerator(patchIndex, volumeIndex, 0),
104 getVolumeCentre(
105 patchCentre,
106 patchSize,
107 QInEnumerator._numberOfDoFsPerAxisInCell,
108 volumeIndex
109 ),
110 getVolumeSize(patchSize, QInEnumerator._numberOfDoFsPerAxisInCell),
111 t,
112 dt,
113 QOut + QOutEnumerator(patchIndex, volumeIndex, 0)
114 );
115
116#if defined(SharedOMP)
117#pragma omp simd
118#endif
119 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
120 QOut[QOutEnumerator(
121 patchIndex,
122 volumeIndex,
123 unknown
124 )] = dt * QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)]
125 + QIn[QInEnumerator(patchIndex, volumeIndex, unknown)];
127 QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)]
128 == QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)],
129 dt,
130 unknown,
131 volumeIndex,
132 QOutEnumerator(patchIndex, volumeIndex, unknown)
133 );
134 }
135#if defined(SharedOMP)
136#pragma omp simd
137#endif
138 for (int unknown = NumberOfUnknowns;
139 unknown < NumberOfUnknowns + NumberOfAuxiliaryVariables;
140 unknown++) {
141 QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)] = QIn
142 [QInEnumerator(patchIndex, volumeIndex, unknown)];
144 QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)]
145 == QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)],
146 dt,
147 unknown,
148 volumeIndex,
149 QOutEnumerator(patchIndex, volumeIndex, unknown)
150 );
151 }
152}
153
154
155template <class SolverType, class QInEnumeratorType, class QOutEnumeratorType>
157 const double* __restrict__ QIn,
158 const QInEnumeratorType& QInEnumerator,
159 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
160 const ::tarch::la::Vector<Dimensions, double>& patchSize,
161 int patchIndex,
162 const ::tarch::la::Vector<Dimensions, int>& volumeIndex,
163 double t,
164 double dt,
165 double* __restrict__ QOut,
166 const QOutEnumeratorType& QOutEnumerator
167) {
168 double QInGathered
169 [SolverType::NumberOfUnknowns + SolverType::NumberOfAuxiliaryVariables]{
170 0.0};
171 double sourceGathered[SolverType::NumberOfUnknowns]{0.0};
172
173 // gather
174#if defined(SharedOMP)
175#pragma omp simd
176#endif
177 for (int unknown = 0;
178 unknown
179 < SolverType::NumberOfUnknowns + SolverType::NumberOfAuxiliaryVariables;
180 unknown++) {
181 QInGathered[unknown] = QIn[QInEnumerator(patchIndex, volumeIndex, unknown)];
182 }
183
184 SolverType::sourceTerm(
185 QInGathered,
186 getVolumeCentre(
187 patchCentre,
188 patchSize,
189 QInEnumerator._numberOfDoFsPerAxisInCell,
190 volumeIndex
191 ),
192 getVolumeSize(patchSize, QInEnumerator._numberOfDoFsPerAxisInCell),
193 t,
194 dt,
195 sourceGathered,
196 Solver::Offloadable::Yes
197 );
198
199 // scatter
200#if defined(SharedOMP)
201#pragma omp simd
202#endif
203 for (int unknown = 0; unknown < SolverType::NumberOfUnknowns; unknown++) {
204 QOut[QOutEnumerator(
205 patchIndex,
206 volumeIndex,
207 unknown
208 )] = QIn[QInEnumerator(patchIndex, volumeIndex, unknown)]
209 + dt * sourceGathered[unknown];
210#if defined(GPUOffloadingOff)
212 QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)]
213 == QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)],
214 dt,
215 unknown,
216 volumeIndex,
217 QOutEnumerator(patchIndex, volumeIndex, unknown)
218 );
219#endif
220 }
221
222#if defined(SharedOMP)
223#pragma omp simd
224#endif
225 for (int unknown = SolverType::NumberOfUnknowns;
226 unknown
227 < SolverType::NumberOfUnknowns + SolverType::NumberOfAuxiliaryVariables;
228 unknown++) {
229 QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)] = QIn
230 [QInEnumerator(patchIndex, volumeIndex, unknown)];
231#if defined(GPUOffloadingOff)
233 QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)]
234 == QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)],
235 dt,
236 unknown,
237 volumeIndex,
238 QOutEnumerator(patchIndex, volumeIndex, unknown)
239 );
240#endif
241 }
242}
243
244
245template <class SolverType>
247 const double* __restrict__ QIn,
248 const enumerator::AoSLexicographicEnumerator& QInEnumerator,
249 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
250 const ::tarch::la::Vector<Dimensions, double>& patchSize,
251 int patchIndex,
252 const ::tarch::la::Vector<Dimensions, int>& volumeIndex,
253 double t,
254 double dt,
255 double* __restrict__ QOut,
256 const enumerator::AoSLexicographicEnumerator& QOutEnumerator
257) {
258 SolverType::sourceTerm(
259 QIn + QInEnumerator(patchIndex, volumeIndex, 0),
260 getVolumeCentre(
261 patchCentre,
262 patchSize,
263 QInEnumerator._numberOfDoFsPerAxisInCell,
264 volumeIndex
265 ),
266 getVolumeSize(patchSize, QInEnumerator._numberOfDoFsPerAxisInCell),
267 t,
268 dt,
269 QOut + QOutEnumerator(patchIndex, volumeIndex, 0),
270 Solver::Offloadable::Yes
271 );
272
273#if defined(SharedOMP)
274#pragma omp simd
275#endif
276 for (int unknown = 0; unknown < SolverType::NumberOfUnknowns; unknown++) {
277 QOut[QOutEnumerator(
278 patchIndex,
279 volumeIndex,
280 unknown
281 )] = dt * QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)]
282 + QIn[QInEnumerator(patchIndex, volumeIndex, unknown)];
283#if defined(GPUOffloadingOff)
285 QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)]
286 == QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)],
287 dt,
288 unknown,
289 volumeIndex,
290 QOutEnumerator(patchIndex, volumeIndex, unknown)
291 );
292#endif
293 }
294#if defined(SharedOMP)
295#pragma omp simd
296#endif
297 for (int unknown = SolverType::NumberOfUnknowns;
298 unknown
299 < SolverType::NumberOfUnknowns + SolverType::NumberOfAuxiliaryVariables;
300 unknown++) {
301 QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)] = QIn
302 [QInEnumerator(patchIndex, volumeIndex, unknown)];
303#if defined(GPUOffloadingOff)
305 QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)]
306 == QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)],
307 dt,
308 unknown,
309 volumeIndex,
310 QOutEnumerator(patchIndex, volumeIndex, unknown)
311 );
312#endif
313 }
314}
315
316
317template <
318 int NumberOfUnknowns,
319 int NumberOfAuxiliaryVariables,
320 class QInEnumeratorType,
321 class EigenvaluesEnumeratorType>
323 const double* __restrict__ QIn,
324 const QInEnumeratorType& QInEnumerator,
325 const EigenvaluesFunctor& eigenvaluesFunctor,
326 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
327 const ::tarch::la::Vector<Dimensions, double>& patchSize,
328 int patchIndex,
329 const ::tarch::la::Vector<Dimensions, int>& volumeIndex,
330 double t,
331 double dt,
332 int normal,
333 double* __restrict__ eigenvalues,
334 const EigenvaluesEnumeratorType& eigenvaluesEnumerator
335) {
336 double QInGathered[NumberOfUnknowns + NumberOfAuxiliaryVariables]{0.0};
337
338 // gather
339#if defined(SharedOMP)
340#pragma omp simd
341#endif
342 for (int unknown = 0; unknown < NumberOfUnknowns + NumberOfAuxiliaryVariables;
343 unknown++) {
344 QInGathered[unknown] = QIn[QInEnumerator(patchIndex, volumeIndex, unknown)];
345 }
346
347 eigenvaluesFunctor(
348 QInGathered,
349 getVolumeCentre(
350 patchCentre,
351 patchSize,
352 QInEnumerator._numberOfDoFsPerAxisInCell,
353 volumeIndex
354 ),
355 getVolumeSize(patchSize, QInEnumerator._numberOfDoFsPerAxisInCell),
356 t,
357 dt,
358 normal,
359 eigenvalues + eigenvaluesEnumerator(patchIndex, volumeIndex, 0)
360 );
361}
362
363
364template <int NumberOfUnknowns, int NumberOfAuxiliaryVariables>
366 const double* __restrict__ QIn,
367 const enumerator::AoSLexicographicEnumerator& QInEnumerator,
368 const EigenvaluesFunctor& eigenvaluesFunctor,
369 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
370 const ::tarch::la::Vector<Dimensions, double>& patchSize,
371 int patchIndex,
372 const ::tarch::la::Vector<Dimensions, int>& volumeIndex,
373 double t,
374 double dt,
375 int normal,
376 double* __restrict__ eigenvalues,
377 const enumerator::AoSLexicographicEnumerator& eigenvaluesEnumerator
378) {
379 eigenvaluesFunctor(
380 QIn + QInEnumerator(patchIndex, volumeIndex, 0),
381 getVolumeCentre(
382 patchCentre,
383 patchSize,
384 QInEnumerator._numberOfDoFsPerAxisInCell,
385 volumeIndex
386 ),
387 getVolumeSize(patchSize, QInEnumerator._numberOfDoFsPerAxisInCell),
388 t,
389 dt,
390 normal,
391 eigenvalues + eigenvaluesEnumerator(patchIndex, volumeIndex, 0)
392 );
393}
394
395
396template <
397 class SolverType,
398 class QInEnumeratorType,
399 class EigenvaluesEnumeratorType>
401 const double* __restrict__ QIn,
402 const QInEnumeratorType& QInEnumerator,
403 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
404 const ::tarch::la::Vector<Dimensions, double>& patchSize,
405 int patchIndex,
406 const ::tarch::la::Vector<Dimensions, int>& volumeIndex,
407 double t,
408 double dt,
409 int normal,
410 double* __restrict__ eigenvalues,
411 const EigenvaluesEnumeratorType& eigenvaluesEnumerator
412) {
413 double QInGathered
414 [SolverType::NumberOfUnknowns + SolverType::NumberOfAuxiliaryVariables]{
415 0.0};
416
417 // gather
418#if defined(SharedOMP)
419#pragma omp simd
420#endif
421 for (int unknown = 0;
422 unknown
423 < SolverType::NumberOfUnknowns + SolverType::NumberOfAuxiliaryVariables;
424 unknown++) {
425 QInGathered[unknown] = QIn[QInEnumerator(patchIndex, volumeIndex, unknown)];
426 }
427
428 SolverType::eigenvalues(
429 QInGathered,
430 getVolumeCentre(
431 patchCentre,
432 patchSize,
433 QInEnumerator._numberOfDoFsPerAxisInCell,
434 volumeIndex
435 ),
436 getVolumeSize(patchSize, QInEnumerator._numberOfDoFsPerAxisInCell),
437 t,
438 dt,
439 normal,
440 eigenvalues + eigenvaluesEnumerator(patchIndex, volumeIndex, 0),
441 Solver::Offloadable::Yes
442 );
443}
444
445
446template <class SolverType>
448 const double* __restrict__ QIn,
449 const enumerator::AoSLexicographicEnumerator& QInEnumerator,
450 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
451 const ::tarch::la::Vector<Dimensions, double>& patchSize,
452 int patchIndex,
453 const ::tarch::la::Vector<Dimensions, int>& volumeIndex,
454 double t,
455 double dt,
456 int normal,
457 double* __restrict__ eigenvalues,
458 const enumerator::AoSLexicographicEnumerator& eigenvaluesEnumerator
459) {
460 SolverType::eigenvalues(
461 QIn + QInEnumerator(patchIndex, volumeIndex, 0),
462 getVolumeCentre(
463 patchCentre,
464 patchSize,
465 QInEnumerator._numberOfDoFsPerAxisInCell,
466 volumeIndex
467 ),
468 getVolumeSize(patchSize, QInEnumerator._numberOfDoFsPerAxisInCell),
469 t,
470 dt,
471 normal,
472 eigenvalues + eigenvaluesEnumerator(patchIndex, volumeIndex, 0),
473 Solver::Offloadable::Yes
474 );
475}
476
477
478template <
479 int NumberOfUnknowns,
480 int NumberOfAuxiliaryVariables,
481 class QInEnumeratorType,
482 class FluxEnumeratorType>
484 const double* __restrict__ QIn,
485 const QInEnumeratorType& QInEnumerator,
486 const FluxFunctor& fluxFunctor,
487 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
488 const ::tarch::la::Vector<Dimensions, double>& patchSize,
489 int patchIndex,
490 const ::tarch::la::Vector<Dimensions, int>& volumeIndex,
491 double t,
492 double dt,
493 int normal,
494 double* __restrict__ flux,
495 const FluxEnumeratorType& fluxEnumerator
496) {
497 double QInGathered[NumberOfUnknowns + NumberOfAuxiliaryVariables]{0.0};
498 double fluxGathered[NumberOfUnknowns]{0.0};
499
500 // gather
501#if defined(SharedOMP)
502#pragma omp simd
503#endif
504 for (int unknown = 0; unknown < NumberOfUnknowns + NumberOfAuxiliaryVariables;
505 unknown++) {
506 QInGathered[unknown] = QIn[QInEnumerator(patchIndex, volumeIndex, unknown)];
507 }
508
509 fluxFunctor(
510 QInGathered,
511 getVolumeCentre(
512 patchCentre,
513 patchSize,
514 QInEnumerator._numberOfDoFsPerAxisInCell,
515 volumeIndex
516 ),
517 getVolumeSize(patchSize, QInEnumerator._numberOfDoFsPerAxisInCell),
518 t,
519 dt,
520 normal,
521 fluxGathered
522 );
523
524 // scatter
525#if defined(SharedOMP)
526#pragma omp simd
527#endif
528 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
529 flux[fluxEnumerator(patchIndex, volumeIndex, unknown)] = fluxGathered
530 [unknown];
531 }
532}
533
534
535template <int NumberOfUnknowns, int NumberOfAuxiliaryVariables>
537 const double* __restrict__ QIn,
538 const enumerator::AoSLexicographicEnumerator& QInEnumerator,
539 const FluxFunctor& fluxFunctor,
540 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
541 const ::tarch::la::Vector<Dimensions, double>& patchSize,
542 int patchIndex,
543 const ::tarch::la::Vector<Dimensions, int>& volumeIndex,
544 double t,
545 double dt,
546 int normal,
547 double* __restrict__ flux,
548 const enumerator::AoSLexicographicEnumerator& fluxEnumerator
549) {
550 fluxFunctor(
551 QIn + QInEnumerator(patchIndex, volumeIndex, 0),
552 getVolumeCentre(
553 patchCentre,
554 patchSize,
555 QInEnumerator._numberOfDoFsPerAxisInCell,
556 volumeIndex
557 ),
558 getVolumeSize(patchSize, QInEnumerator._numberOfDoFsPerAxisInCell),
559 t,
560 dt,
561 normal,
562 flux + fluxEnumerator(patchIndex, volumeIndex, 0)
563 );
564}
565
566
567template <class SolverType, class QInEnumeratorType, class FluxEnumeratorType>
569 const double* __restrict__ QIn,
570 const QInEnumeratorType& QInEnumerator,
571 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
572 const ::tarch::la::Vector<Dimensions, double>& patchSize,
573 int patchIndex,
574 const ::tarch::la::Vector<Dimensions, int>& volumeIndex,
575 double t,
576 double dt,
577 int normal,
578 double* __restrict__ flux,
579 const FluxEnumeratorType& fluxEnumerator
580) {
581 double QInGathered
582 [SolverType::NumberOfUnknowns + SolverType::NumberOfAuxiliaryVariables]{
583 0.0};
584 double fluxGathered[SolverType::NumberOfUnknowns]{0.0};
585
586 // gather
587#if defined(SharedOMP)
588#pragma omp simd
589#endif
590 for (int unknown = 0;
591 unknown
592 < SolverType::NumberOfUnknowns + SolverType::NumberOfAuxiliaryVariables;
593 unknown++) {
594 QInGathered[unknown] = QIn[QInEnumerator(patchIndex, volumeIndex, unknown)];
595 }
596
597 SolverType::flux(
598 QInGathered,
599 getVolumeCentre(
600 patchCentre,
601 patchSize,
602 QInEnumerator._numberOfDoFsPerAxisInCell,
603 volumeIndex
604 ),
605 getVolumeSize(patchSize, QInEnumerator._numberOfDoFsPerAxisInCell),
606 t,
607 dt,
608 normal,
609 fluxGathered,
610 Solver::Offloadable::Yes
611 );
612
613 // scatter
614#if defined(SharedOMP)
615#pragma omp simd
616#endif
617 for (int unknown = 0; unknown < SolverType::NumberOfUnknowns; unknown++) {
618 flux[fluxEnumerator(patchIndex, volumeIndex, unknown)] = fluxGathered
619 [unknown];
620 }
621}
622
623
624template <class SolverType, class FluxEnumeratorType>
626 const double* __restrict__ QIn,
627 const enumerator::AoSLexicographicEnumerator& QInEnumerator,
628 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
629 const ::tarch::la::Vector<Dimensions, double>& patchSize,
630 int patchIndex,
631 const ::tarch::la::Vector<Dimensions, int>& volumeIndex,
632 double t,
633 double dt,
634 int normal,
635 double* __restrict__ flux,
636 const FluxEnumeratorType& fluxEnumerator
637) {
638 double fluxGathered[SolverType::NumberOfUnknowns]{0.0};
639
640 SolverType::flux(
641 QIn + QInEnumerator(patchIndex, volumeIndex, 0),
642 getVolumeCentre(
643 patchCentre,
644 patchSize,
645 QInEnumerator._numberOfDoFsPerAxisInCell,
646 volumeIndex
647 ),
648 getVolumeSize(patchSize, QInEnumerator._numberOfDoFsPerAxisInCell),
649 t,
650 dt,
651 normal,
652 fluxGathered,
653 Solver::Offloadable::Yes
654 );
655
656 // scatter
657#if defined(SharedOMP)
658#pragma omp simd
659#endif
660 for (int unknown = 0; unknown < SolverType::NumberOfUnknowns; unknown++) {
661 flux[fluxEnumerator(patchIndex, volumeIndex, unknown)] = fluxGathered
662 [unknown];
663 }
664}
665
666
667template <class SolverType>
669 const double* __restrict__ QIn,
670 const enumerator::AoSLexicographicEnumerator& QInEnumerator,
671 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
672 const ::tarch::la::Vector<Dimensions, double>& patchSize,
673 int patchIndex,
674 const ::tarch::la::Vector<Dimensions, int>& volumeIndex,
675 double t,
676 double dt,
677 int normal,
678 double* __restrict__ flux,
679 const enumerator::AoSLexicographicEnumerator& fluxEnumerator
680) {
681 SolverType::flux(
682 QIn + QInEnumerator(patchIndex, volumeIndex, 0),
683 getVolumeCentre(
684 patchCentre,
685 patchSize,
686 QInEnumerator._numberOfDoFsPerAxisInCell,
687 volumeIndex
688 ),
689 getVolumeSize(patchSize, QInEnumerator._numberOfDoFsPerAxisInCell),
690 t,
691 dt,
692 normal,
693 flux + fluxEnumerator(patchIndex, volumeIndex, 0),
694 Solver::Offloadable::Yes
695 );
696}
697
698
699template <
700 int NumberOfUnknowns,
701 int NumberOfAuxiliaryVariables,
702 class QInEnumeratorType,
703 class NCPFaceEnumeratorType>
705 const double* __restrict__ QIn,
706 const QInEnumeratorType& QInEnumerator,
707 const NonconservativeProductFunctor& ncpFunctor,
708 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
709 const ::tarch::la::Vector<Dimensions, double>& patchSize,
710 int patchIndex,
711 const ::tarch::la::Vector<Dimensions, int>& volumeIndex,
712 double t,
713 double dt,
714 int normal,
715 double* __restrict__ ncp,
716 const NCPFaceEnumeratorType& ncpEnumerator
717) {
718 double QAverage[NumberOfUnknowns + NumberOfAuxiliaryVariables]{0.0};
719 double DeltaQ[NumberOfUnknowns + NumberOfAuxiliaryVariables]{0.0};
720 double ncpGathered[NumberOfUnknowns]{0.0};
721
722 const ::tarch::la::Vector<Dimensions, int> leftAdjacentVolume = volumeIndex;
724 rightAdjacentVolume(normal)++;
725
726 // gather
727#if defined(SharedOMP)
728#pragma omp simd
729#endif
730 for (int unknown = 0; unknown < NumberOfUnknowns + NumberOfAuxiliaryVariables;
731 unknown++) {
732 QAverage[unknown]
733 = 0.5 * QIn[QInEnumerator(patchIndex, leftAdjacentVolume, unknown)]
734 + 0.5 * QIn[QInEnumerator(patchIndex, rightAdjacentVolume, unknown)];
735 DeltaQ[unknown] = QIn[QInEnumerator(patchIndex, rightAdjacentVolume, unknown)]
736 - QIn
737 [QInEnumerator(patchIndex, leftAdjacentVolume, unknown)];
738 }
739
740 ncpFunctor(
741 QAverage,
742 DeltaQ,
743 getVolumeCentre(
744 patchCentre,
745 patchSize,
746 QInEnumerator._numberOfDoFsPerAxisInCell,
747 volumeIndex
748 ),
749 getVolumeSize(patchSize, QInEnumerator._numberOfDoFsPerAxisInCell),
750 t,
751 dt,
752 normal,
753 ncpGathered
754 );
755
756 // scatter
757#if defined(SharedOMP)
758#pragma omp simd
759#endif
760 for (int unknown = 0; unknown < NumberOfUnknowns; unknown++) {
761 ncp[ncpEnumerator(patchIndex, volumeIndex, unknown)] = ncpGathered[unknown];
762 }
763}
764
765
766template <class SolverType, class QInEnumeratorType, class NCPFaceEnumeratorType>
768 const double* __restrict__ QIn,
769 const QInEnumeratorType& QInEnumerator,
770 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
771 const ::tarch::la::Vector<Dimensions, double>& patchSize,
772 int patchIndex,
773 const ::tarch::la::Vector<Dimensions, int>& volumeIndex,
774 double t,
775 double dt,
776 int normal,
777 double* __restrict__ ncp,
778 const NCPFaceEnumeratorType& ncpEnumerator
779) {
780 double QAverage
781 [SolverType::NumberOfUnknowns + SolverType::NumberOfAuxiliaryVariables]{
782 0.0};
783 double DeltaQ
784 [SolverType::NumberOfUnknowns + SolverType::NumberOfAuxiliaryVariables]{
785 0.0};
786 double ncpGathered[SolverType::NumberOfUnknowns]{0.0};
787
788 const ::tarch::la::Vector<Dimensions, int> leftAdjacentVolume = volumeIndex;
790 rightAdjacentVolume(normal)++;
791
792 // gather
793#if defined(SharedOMP)
794#pragma omp simd
795#endif
796 for (int unknown = 0;
797 unknown
798 < SolverType::NumberOfUnknowns + SolverType::NumberOfAuxiliaryVariables;
799 unknown++) {
800 QAverage[unknown]
801 = 0.5 * QIn[QInEnumerator(patchIndex, leftAdjacentVolume, unknown)]
802 + 0.5 * QIn[QInEnumerator(patchIndex, rightAdjacentVolume, unknown)];
803 DeltaQ[unknown] = QIn[QInEnumerator(patchIndex, rightAdjacentVolume, unknown)]
804 - QIn
805 [QInEnumerator(patchIndex, leftAdjacentVolume, unknown)];
806 }
807
808 SolverType::nonconservativeProduct(
809 QAverage,
810 DeltaQ,
811 getVolumeCentre(
812 patchCentre,
813 patchSize,
814 QInEnumerator._numberOfDoFsPerAxisInCell,
815 volumeIndex
816 ),
817 getVolumeSize(patchSize, QInEnumerator._numberOfDoFsPerAxisInCell),
818 t,
819 dt,
820 normal,
821 ncpGathered,
822 Solver::Offloadable::Yes
823 );
824
825 // scatter
826#if defined(SharedOMP)
827#pragma omp simd
828#endif
829 for (int unknown = 0; unknown < SolverType::NumberOfUnknowns; unknown++) {
830 ncp[ncpEnumerator(patchIndex, volumeIndex, unknown)] = ncpGathered[unknown];
831 }
832}
833
834
835template <class NCPFaceEnumeratorType, class QOutEnumeratorType>
837 const double* __restrict__ ncpX,
838 const double* __restrict__ ncpY,
839 const double* __restrict__ ncpZ,
840 const NCPFaceEnumeratorType& ncpEnumerator,
841 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
842 const ::tarch::la::Vector<Dimensions, double>& patchSize,
843 int patchIndex,
844 const ::tarch::la::Vector<Dimensions, int>& volumeIndex,
845 int unknown,
846 double dt,
847 double* __restrict__ QOut,
848 const QOutEnumeratorType& QOutEnumerator
849) {
850 auto updateAlongOneCoordinateDirection =
851 [=](const double* __restrict__ ncp, int normal) {
853 const ::tarch::la::Vector<Dimensions, int> rightFace = volumeIndex;
854 leftFace(normal)--;
855
856 // The ncp is calculating (Q^+ - Q^-)B_i for the considered face, and B
857 // here is defined on the left of the equation. For ncpLeft, it is the
858 // ncp^+ of the left face, so it have a plus sign for ncpRight, it is the
859 // ncp^- of the right face, so it have a minus sign.
860 const double ncpLeft = +ncp[ncpEnumerator(patchIndex, leftFace, unknown)];
861 const double ncpRight = -ncp[ncpEnumerator(
862 patchIndex,
863 rightFace,
864 unknown
865 )];
866
867 // Left ncp have another minus sign as the normal vector direction for the
868 // left face goes against the flux.
869 QOut[QOutEnumerator(
870 patchIndex,
871 volumeIndex,
872 unknown
873 )] += dt / patchSize(normal) * QOutEnumerator._numberOfDoFsPerAxisInCell
874 * (-ncpLeft + ncpRight) * (0.5);
875
876#if defined(GPUOffloadingOff)
878 QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)]
879 == QOut[QOutEnumerator(patchIndex, volumeIndex, unknown)],
880 ncpLeft,
881 ncpRight,
882 dt,
883 normal,
884 unknown,
885 patchSize,
886 volumeIndex,
887 QOutEnumerator(patchIndex, volumeIndex, unknown)
888 );
889#endif
890 };
891
892 updateAlongOneCoordinateDirection(ncpX, 0);
893 updateAlongOneCoordinateDirection(ncpY, 1);
894#if Dimensions == 3
895 updateAlongOneCoordinateDirection(ncpZ, 2);
896#endif
897}
898
899
900template <class NCPFaceEnumeratorType, class QOutEnumeratorType>
902 const double* __restrict__ ncp,
903 const NCPFaceEnumeratorType& ncpEnumerator,
904 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
905 const ::tarch::la::Vector<Dimensions, double>& patchSize,
906 int patchIndex,
907 const ::tarch::la::Vector<Dimensions, int>& faceIndex,
908 int unknown,
909 double dt,
910 int normal,
911 double* __restrict__ QOut,
912 const QOutEnumeratorType& QOutEnumerator
913) {
914 const ::tarch::la::Vector<Dimensions, int> leftVolume = faceIndex;
915 ::tarch::la::Vector<Dimensions, int> rightVolume = faceIndex;
916 rightVolume(normal)++;
917
918 if (leftVolume(normal) >= 0) {
919 QOut[QOutEnumerator(
920 patchIndex,
921 leftVolume,
922 unknown
923 )] -= dt / patchSize(normal) * QOutEnumerator._numberOfDoFsPerAxisInCell
924 * ncp[ncpEnumerator(patchIndex, faceIndex, unknown)] * 0.5;
925 }
926 if (rightVolume(normal) < QOutEnumerator._numberOfDoFsPerAxisInCell) {
927 QOut[QOutEnumerator(
928 patchIndex,
929 rightVolume,
930 unknown
931 )] -= dt / patchSize(normal) * QOutEnumerator._numberOfDoFsPerAxisInCell
932 * ncp[ncpEnumerator(patchIndex, faceIndex, unknown)] * 0.5;
933 }
934}
935
936
937template <int NumberOfUnknowns, int NumberOfAuxiliaryVariables>
939 const double* __restrict__ QIn,
940 const enumerator::AoSLexicographicEnumerator& QInEnumerator,
941 const RiemannFunctor& riemannFunctor,
942 const double* __restrict__ flux,
943 const enumerator::AoSLexicographicEnumerator& fluxEnumerator,
944 const double* __restrict__ eigenvalues,
945 const enumerator::AoSLexicographicEnumerator& eigenvaluesEnumerator,
946 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
947 const ::tarch::la::Vector<Dimensions, double>& patchSize,
948 int patchIndex,
949 const ::tarch::la::Vector<Dimensions, int>& volumeIndex,
950 double t,
951 double dt,
952 int normal,
953 double* __restrict__ APDQ,
954 double* __restrict__ AMDQ,
955 const enumerator::AoSLexicographicEnumerator& DQEnumerator,
956 double* __restrict__ maxEigenvalue,
957 const enumerator::AoSLexicographicEnumerator& maxEigenvalueEnumerator
958) {
961
962 rightVolume(normal)++;
963
964 maxEigenvalue[maxEigenvalueEnumerator(patchIndex, leftVolume, 0)] = riemannFunctor(
965 QIn + QInEnumerator(patchIndex, rightVolume, 0),
966 QIn + QInEnumerator(patchIndex, leftVolume, 0),
967 flux + fluxEnumerator(patchIndex, rightVolume, 0),
968 flux + fluxEnumerator(patchIndex, leftVolume, 0),
969 eigenvalues + eigenvaluesEnumerator(patchIndex, rightVolume, 0),
970 eigenvalues + eigenvaluesEnumerator(patchIndex, leftVolume, 0),
971 getVolumeCentre(
972 patchCentre,
973 patchSize,
974 QInEnumerator._numberOfDoFsPerAxisInCell,
975 rightVolume
976 ),
977 getVolumeCentre(
978 patchCentre,
979 patchSize,
980 QInEnumerator._numberOfDoFsPerAxisInCell,
981 leftVolume
982 ),
983 getVolumeSize(patchSize, QInEnumerator._numberOfDoFsPerAxisInCell),
984 t,
985 dt,
986 normal,
987 APDQ + DQEnumerator(patchIndex, leftVolume, 0),
988 AMDQ + DQEnumerator(patchIndex, leftVolume, 0)
989 );
990}
991
992
993template <class SolverType>
995 const double* __restrict__ QIn,
996 const enumerator::AoSLexicographicEnumerator& QInEnumerator,
997 const double* __restrict__ flux,
998 const enumerator::AoSLexicographicEnumerator& fluxEnumerator,
999 const double* __restrict__ eigenvalues,
1000 const enumerator::AoSLexicographicEnumerator& eigenvaluesEnumerator,
1001 const ::tarch::la::Vector<Dimensions, double>& patchCentre,
1002 const ::tarch::la::Vector<Dimensions, double>& patchSize,
1003 int patchIndex,
1004 const ::tarch::la::Vector<Dimensions, int>& volumeIndex,
1005 double t,
1006 double dt,
1007 int normal,
1008 double* __restrict__ APDQ,
1009 double* __restrict__ AMDQ,
1010 const enumerator::AoSLexicographicEnumerator& DQEnumerator,
1011 double* __restrict__ maxEigenvalue,
1012 const enumerator::AoSLexicographicEnumerator& maxEigenvalueEnumerator
1013) {
1016
1017 rightVolume(normal)++;
1018
1019 maxEigenvalue[maxEigenvalueEnumerator(
1020 patchIndex,
1021 leftVolume,
1022 0
1023 )] = SolverType::
1024 solveRiemannProblem(
1025 QIn + QInEnumerator(patchIndex, rightVolume, 0),
1026 QIn + QInEnumerator(patchIndex, leftVolume, 0),
1027 flux + fluxEnumerator(patchIndex, rightVolume, 0),
1028 flux + fluxEnumerator(patchIndex, leftVolume, 0),
1029 eigenvalues + eigenvaluesEnumerator(patchIndex, rightVolume, 0),
1030 eigenvalues + eigenvaluesEnumerator(patchIndex, leftVolume, 0),
1031 getVolumeCentre(
1032 patchCentre,
1033 patchSize,
1034 QInEnumerator._numberOfDoFsPerAxisInCell,
1035 rightVolume
1036 ),
1037 getVolumeCentre(
1038 patchCentre,
1039 patchSize,
1040 QInEnumerator._numberOfDoFsPerAxisInCell,
1041 leftVolume
1042 ),
1043 getVolumeSize(patchSize, QInEnumerator._numberOfDoFsPerAxisInCell),
1044 t,
1045 dt,
1046 normal,
1047 APDQ + DQEnumerator(patchIndex, leftVolume, 0),
1048 AMDQ + DQEnumerator(patchIndex, leftVolume, 0),
1049 Solver::Offloadable::Yes
1050 );
1051}
1052
1053
1054template <class RiemannEnumeratorType, class QOutEnumeratorType>
1056 const double* __restrict__ leftUpdatesX,
1057 const double* __restrict__ rightUpdatesX,
1058 const double* __restrict__ belowUpdatesY,
1059 const double* __restrict__ aboveUpdatesY,
1060 const double* __restrict__ backwardUpdatesZ,
1061 const double* __restrict__ forwardUpdatesZ,
1062 const RiemannEnumeratorType& riemannEnumerator,
1063 const tarch::la::Vector<Dimensions, double>& patchCentre,
1065 int patchIndex,
1066 const tarch::la::Vector<Dimensions, int>& volumeIndex,
1067 int unknown,
1068 double dt,
1069 double* __restrict__ QOut,
1070 const QOutEnumeratorType& QOutEnumerator
1071) {
1072 auto updateAlongOneDirection =
1073 [=](
1074 const double* __restrict__ leftUpdates,
1075 const double* __restrict__ rightUpdates,
1076 const int normal
1077 ) {
1080 leftVolume(normal)--;
1081
1082 const double fluxPos = rightUpdates
1083 [riemannEnumerator(patchIndex, leftVolume, unknown)];
1084 const double fluxNeg = leftUpdates
1085 [riemannEnumerator(patchIndex, rightVolume, unknown)];
1086
1087 QOut[QOutEnumerator(
1088 patchIndex,
1089 volumeIndex,
1090 unknown
1091 )] -= dt / patchSize(normal) * QOutEnumerator._numberOfDoFsPerAxisInCell
1092 * (fluxNeg - fluxPos);
1093 };
1094
1095 updateAlongOneDirection(leftUpdatesX, rightUpdatesX, 0);
1096 updateAlongOneDirection(belowUpdatesY, aboveUpdatesY, 1);
1097#if Dimensions == 3
1098 updateAlongOneDirection(backwardUpdatesZ, forwardUpdatesZ, 2);
1099#endif
1100}
1101
1102
1104 const double* __restrict__ maxEigenvalueX,
1105 const double* __restrict__ maxEigenvalueY,
1106 const double* __restrict__ maxEigenvalueZ,
1107 enumerator::AoSLexicographicEnumerator maxEigenvalueEnumerator,
1108 int patchIndex,
1109 const tarch::la::Vector<Dimensions, int>& volumeIndex
1110) {
1111 double result = 0.0;
1112 result = std::max(
1113 result,
1114 maxEigenvalueX[maxEigenvalueEnumerator(patchIndex, volumeIndex, 0)]
1115 );
1116 result = std::max(
1117 result,
1118 maxEigenvalueY[maxEigenvalueEnumerator(patchIndex, volumeIndex, 0)]
1119 );
1120#if Dimensions == 3
1121 result = std::max(
1122 result,
1123 maxEigenvalueZ[maxEigenvalueEnumerator(patchIndex, volumeIndex, 0)]
1124 );
1125#endif
1126 return result;
1127}
#define assertion4(expr, param0, param1, param2, param3)
#define assertion8(expr, param0, param1, param2, param3, param4, param5, param6, param7)
float dt
Definition DSL_test.py:5
KeywordToAvoidDuplicateSymbolsForInlinedFunctions double maxEigenvalue(const double *const Q, int normal, const double CCZ4e, const double CCZ4ds, const double CCZ4GLMc, const double CCZ4GLMd) InlineMethod
KeywordToAvoidDuplicateSymbolsForInlinedFunctions void copySolutionAndAddSourceTerm(const double *__restrict__ QIn, const QInEnumeratorType &QInEnumerator, const SourceFunctor &sourceFunctor, const ::tarch::la::Vector< Dimensions, double > &patchCentre, const ::tarch::la::Vector< Dimensions, double > &patchSize, int patchIndex, const ::tarch::la::Vector< Dimensions, int > &volumeIndex, double t, double dt, double *__restrict__ QOut, const QOutEnumeratorType &QOutEnumerator) InlineMethod
KeywordToAvoidDuplicateSymbolsForInlinedFunctions GPUCallableInlineMethod void updateSolutionWithNonconservativeFlux(const double *__restrict__ ncpX, const double *__restrict__ ncpY, const double *__restrict__ ncpZ, const NCPFaceEnumeratorType &ncpEnumerator, const ::tarch::la::Vector< Dimensions, double > &patchCentre, const ::tarch::la::Vector< Dimensions, double > &patchSize, int patchIndex, const ::tarch::la::Vector< Dimensions, int > &volumeIndex, int unknown, double dt, double *__restrict__ QOut, const QOutEnumeratorType &QOutEnumerator) InlineMethod
KeywordToAvoidDuplicateSymbolsForInlinedFunctions double reduceMaxEigenvalue(const double *__restrict__ maxEigenvalueX, const double *__restrict__ maxEigenvalueY, const double *__restrict__ maxEigenvalueZ, enumerator::AoSLexicographicEnumerator maxEigenvalueEnumerator, int patchIndex, const tarch::la::Vector< Dimensions, int > &volumeIndex) InlineMethod
KeywordToAvoidDuplicateSymbolsForInlinedFunctions void computeRiemannSolution(const double *__restrict__ QIn, const enumerator::AoSLexicographicEnumerator &QInEnumerator, const RiemannFunctor &riemannFunctor, const double *__restrict__ flux, const enumerator::AoSLexicographicEnumerator &fluxEnumerator, const double *__restrict__ eigenvalues, const enumerator::AoSLexicographicEnumerator &eigenvaluesEnumerator, const ::tarch::la::Vector< Dimensions, double > &patchCentre, const ::tarch::la::Vector< Dimensions, double > &patchSize, int patchIndex, const ::tarch::la::Vector< Dimensions, int > &volumeIndex, double t, double dt, int normal, double *__restrict__ APDQ, double *__restrict__ AMDQ, const enumerator::AoSLexicographicEnumerator &DQEnumerator, double *__restrict__ maxEigenvalue, const enumerator::AoSLexicographicEnumerator &maxEigenvalueEnumerator) InlineMethod
KeywordToAvoidDuplicateSymbolsForInlinedFunctions void computeEigenvalues(const double *__restrict__ QIn, const QInEnumeratorType &QInEnumerator, const EigenvaluesFunctor &eigenvaluesFunctor, const ::tarch::la::Vector< Dimensions, double > &patchCentre, const ::tarch::la::Vector< Dimensions, double > &patchSize, int patchIndex, const ::tarch::la::Vector< Dimensions, int > &volumeIndex, double t, double dt, int normal, double *__restrict__ eigenvalues, const EigenvaluesEnumeratorType &eigenvaluesEnumerator) InlineMethod
KeywordToAvoidDuplicateSymbolsForInlinedFunctions void updateSolutionWithRiemannSolution(const double *__restrict__ leftUpdatesX, const double *__restrict__ rightUpdatesX, const double *__restrict__ belowUpdatesY, const double *__restrict__ aboveUpdatesY, const double *__restrict__ backwardUpdatesZ, const double *__restrict__ forwardUpdatesZ, const RiemannEnumeratorType &riemannEnumerator, const tarch::la::Vector< Dimensions, double > &patchCentre, const tarch::la::Vector< Dimensions, double > &patchSize, int patchIndex, const tarch::la::Vector< Dimensions, int > &volumeIndex, int unknown, double dt, double *__restrict__ QOut, const QOutEnumeratorType &QOutEnumerator) InlineMethod
KeywordToAvoidDuplicateSymbolsForInlinedFunctions void computeFlux(const double *__restrict__ QIn, const QInEnumeratorType &QInEnumerator, const FluxFunctor &fluxFunctor, const ::tarch::la::Vector< Dimensions, double > &patchCentre, const ::tarch::la::Vector< Dimensions, double > &patchSize, int patchIndex, const ::tarch::la::Vector< Dimensions, int > &volumeIndex, double t, double dt, int normal, double *__restrict__ flux, const FluxEnumeratorType &fluxEnumerator) InlineMethod
KeywordToAvoidDuplicateSymbolsForInlinedFunctions void computeNonconservativeFlux(const double *__restrict__ QIn, const QInEnumeratorType &QInEnumerator, const NonconservativeProductFunctor &ncpFunctor, const ::tarch::la::Vector< Dimensions, double > &patchCentre, const ::tarch::la::Vector< Dimensions, double > &patchSize, int patchIndex, const ::tarch::la::Vector< Dimensions, int > &volumeIndex, double t, double dt, int normal, double *__restrict__ ncp, const NCPFaceEnumeratorType &ncpEnumerator) InlineMethod
auto volumeIndex(Args... args)
Definition VolumeIndex.h:54
ncp
Definition swe.py:172
Simple vector class.
Definition Vector.h:134