Peano
Loading...
Searching...
No Matches
kernels.py
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
3import jinja2
4
5from enum import Enum
6
7from exahype2.solvers.PDETerms import PDETerms
8import exahype2
9
10from exahype2.solvers.LagrangeBasisWithDiagonalMassMatrix import GaussLegendreBasis
11
12
13class SolverVariant(Enum):
14 WithVirtualFunctions = 0
15 Stateless = 1
16 Accelerator = 2
17
18
19def create_abstract_solver_declarations(
20 flux_implementation,
21 ncp_implementation,
22 eigenvalues_implementation,
23 source_term_implementation,
24 point_source_implementation,
25 pde_terms_without_state,
26):
27 Template = jinja2.Template(
28 """
29 public:
30 {% if EIGENVALUES_IMPLEMENTATION=="<none>" or EIGENVALUES_IMPLEMENTATION=="<empty>" %}
31 #error Eigenvalue implementation cannot be none or empty
32 {% endif %}
33
34 {% if EIGENVALUES_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
35 /**
36 * Depending on the implementation, this variant might be slow as it
37 * lacks an inline define. Also, if you don't want to use ipo aggressively,
38 * it might be clever to put the implementation into the header.
39 */
40 #if defined(OpenMPGPUOffloading)
41 #pragma omp declare target
42 #endif
43 static double maxEigenvalue(
44 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
45 const tarch::la::Vector<Dimensions,double>& faceCentre,
46 double t,
47 double dt,
48 int normal,
49 Offloadable
50 );
51 #if defined(OpenMPGPUOffloading)
52 #pragma omp end declare target
53 #endif
54 {% endif %}
55
56 /**
57 * Determine max eigenvalue over Jacobian in a given point with solution values
58 * (states) Q. All parameters are in.
59 *
60 * @return Max eigenvalue. Result has to be positive, so we are actually speaking
61 * about the maximum absolute eigenvalue.
62 */
63 virtual double maxEigenvalue(
64 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
65 const tarch::la::Vector<Dimensions,double>& faceCentre,
66 double t,
67 double dt,
68 int normal
69 ) {% if EIGENVALUES_IMPLEMENTATION=="<user-defined>" %}= 0{% else %} final{% endif %};
70
71 {% if FLUX_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
72 /**
73 * Depending on the implementation, this variant might be slow as it
74 * lacks an inline define. Also, if you don't want to use ipo aggressively,
75 * it might be clever to put the implementation into the header.
76 */
77 #if defined(OpenMPGPUOffloading)
78 #pragma omp declare target
79 #endif
80 static void flux(
81 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
82 const tarch::la::Vector<Dimensions,double>& faceCentre,
83 double t,
84 double dt,
85 int normal,
86 double * __restrict__ F, // F[{{NUMBER_OF_UNKNOWNS}}]
87 Offloadable
88 );
89 #if defined(OpenMPGPUOffloading)
90 #pragma omp end declare target
91 #endif
92 {% endif %}
93
94 {% if FLUX_IMPLEMENTATION!="<none>" %}
95 virtual void flux(
96 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
97 const tarch::la::Vector<Dimensions,double>& faceCentre,
98 double t,
99 double dt,
100 int normal,
101 double * __restrict__ F // F[{{NUMBER_OF_UNKNOWNS}}]
102 ) {% if FLUX_IMPLEMENTATION=="<user-defined>" %}=0{% else %} final {% endif %};
103 {% endif %}
104
105 {% if NCP_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
106 /**
107 * Depending on the implementation, this variant might be slow as it
108 * lacks an inline define. Also, if you don't want to use ipo aggressively,
109 * it might be clever to put the implementation into the header.
110 */
111 #if defined(OpenMPGPUOffloading)
112 #pragma omp declare target
113 #endif
114 static void nonconservativeProduct(
115 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
116 const double * __restrict__ deltaQ, // deltaQ[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
117 const tarch::la::Vector<Dimensions,double>& faceCentre,
118 double t,
119 double dt,
120 int normal,
121 double * __restrict__ BTimesDeltaQ, // BTimesDeltaQ[{{NUMBER_OF_UNKNOWNS}}]
122 Offloadable
123 );
124 #if defined(OpenMPGPUOffloading)
125 #pragma omp end declare target
126 #endif
127 {% endif %}
128
129 {% if NCP_IMPLEMENTATION!="<none>" %}
130 virtual void nonconservativeProduct(
131 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
132 const double * __restrict__ deltaQ, // [{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}]
133 const tarch::la::Vector<Dimensions,double>& faceCentre,
134 double t,
135 double dt,
136 int normal,
137 double * __restrict__ BTimesDeltaQ // BTimesDeltaQ[{{NUMBER_OF_UNKNOWNS}}]
138 ) {% if NCP_IMPLEMENTATION=="<user-defined>" %}=0{% endif %};
139 {% endif %}
140
141 {% if SOURCE_TERM_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
142 /**
143 * Depending on the implementation, this variant might be slow as it
144 * lacks an inline define. Also, if you don't want to use ipo aggressively,
145 * it might be clever to put the implementation into the header.
146 */
147 #if defined(OpenMPGPUOffloading)
148 #pragma omp declare target
149 #endif
150 static void sourceTerm(
151 const double * __restrict__ Q,
152 const tarch::la::Vector<Dimensions,double>& volumeCentre,
153 double t,
154 double dt,
155 double * __restrict__ S,
156 Offloadable
157 );
158 #if defined(OpenMPGPUOffloading)
159 #pragma omp end declare target
160 #endif
161 {% endif %}
162
163 {% if SOURCE_TERM_IMPLEMENTATION!="<none>" %}
164 virtual void sourceTerm(
165 const double * __restrict__ Q,
166 const tarch::la::Vector<Dimensions,double>& volumeCentre,
167 double t,
168 double dt,
169 double * __restrict__ S
170 ) {% if SOURCE_TERM_IMPLEMENTATION=="<user-defined>" %}= 0{% else %} final {% endif %};
171 {% endif %}
172
173 {% if POINT_SOURCE_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
174 /**
175 * Depending on the implementation, this variant might be slow as it
176 * lacks an inline define. Also, if you don't want to use ipo aggressively,
177 * it might be clever to put the implementation into the header.
178 */
179 #if defined(OpenMPGPUOffloading)
180 #pragma omp declare target
181 #endif
182 void pointSource(
183 const double * __restrict__ Q,
184 const tarch::la::Vector<Dimensions,double>& x,
185 double t,
186 double dt,
187 double * __restrict__ S,
188 Offloadable
189 );
190 #if defined(OpenMPGPUOffloading)
191 #pragma omp end declare target
192 #endif
193 {% endif %}
194
195 {% if POINT_SOURCE_IMPLEMENTATION!="<none>" %}
196 virtual void pointSource(
197 const double * __restrict__ Q,
198 const tarch::la::Vector<Dimensions,double>& x,
199 double t,
200 double dt,
201 double * __restrict__ S
202 ) {% if SOURCE_TERM_IMPLEMENTATION=="<user-defined>" %}= 0{% else %} final {% endif %};
203 {% endif %}
204""",
205 undefined=jinja2.DebugUndefined,
206 )
207
208 d = {}
209 d["FLUX_IMPLEMENTATION"] = flux_implementation
210 d["NCP_IMPLEMENTATION"] = ncp_implementation
211 d["EIGENVALUES_IMPLEMENTATION"] = eigenvalues_implementation
212 d["SOURCE_TERM_IMPLEMENTATION"] = source_term_implementation
213 d["POINT_SOURCE_IMPLEMENTATION"] = point_source_implementation
214 d["STATELESS_PDE_TERMS"] = pde_terms_without_state
215 return Template.render(**d)
216
217
218def create_abstract_solver_definitions(
219 flux_implementation,
220 ncp_implementation,
221 eigenvalues_implementation,
222 source_term_implementation,
223 point_source_implementation,
224 pde_terms_without_state,
225):
226 Template = jinja2.Template(
227 """
228{% if EIGENVALUES_IMPLEMENTATION!="<user-defined>" and EIGENVALUES_IMPLEMENTATION!="<none>" %}
229double {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::maxEigenvalue(
230 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
231 const tarch::la::Vector<Dimensions,double>& faceCentre,
232 double t,
233 double dt,
234 int normal
235) {
236 {{EIGENVALUES_IMPLEMENTATION}}
237}
238{% endif %}
239
240{% if FLUX_IMPLEMENTATION!="<none>" and FLUX_IMPLEMENTATION!="<user-defined>" %}
241void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::flux(
242 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
243 const tarch::la::Vector<Dimensions,double>& faceCentre,
244 double t,
245 double dt,
246 int normal,
247 double * __restrict__ F // F[{{NUMBER_OF_UNKNOWNS}}]
248) {
249 {{FLUX_IMPLEMENTATION}}
250}
251{% endif %}
252
253{% if NCP_IMPLEMENTATION!="<none>" and NCP_IMPLEMENTATION!="<user-defined>" %}
254void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::nonconservativeProduct(
255 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
256 const double * __restrict__ deltaQ, // deltaQ[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}]
257 const tarch::la::Vector<Dimensions,double>& faceCentre,
258 double t,
259 double dt,
260 int normal,
261 double * __restrict__ BTimesDeltaQ // BTimesDeltaQ[{{NUMBER_OF_UNKNOWNS}}]
262) {
263 {{NCP_IMPLEMENTATION}}
264}
265{% endif %}
266
267{% if SOURCE_TERM_IMPLEMENTATION!="<user-defined>" and SOURCE_TERM_IMPLEMENTATION!="<none>" %}
268void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::sourceTerm(
269 const double * __restrict__ Q,
270 const tarch::la::Vector<Dimensions,double>& volumeCentre,
271 double t,
272 double dt,
273 double * __restrict__ S
274) {
275 {% if SOURCE_TERM_IMPLEMENTATION!="<empty>" %}
276 {{SOURCE_TERM_IMPLEMENTATION}}
277 {% else %}
278 std::fill_n(S,{{NUMBER_OF_UNKNOWNS}},0.0);
279 {% endif %}
280}
281{% endif %}
282
283{% if POINT_SOURCE_IMPLEMENTATION!="<user-defined>" and POINT_SOURCE_IMPLEMENTATION!="<none>" %}
284#error Point sources not yet implemented
285void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::pointSource(
286 const double * __restrict__ Q,
287 const tarch::la::Vector<Dimensions,double>& x,
288 double t,
289 double dt,
290 double * __restrict__ S
291) {
292 {% if POINT_SOURCE_IMPLEMENTATION!="<empty>" %}
293 {{POINT_SOURCE_IMPLEMENTATION}}
294 {% else %}
295 std::fill_n(S,{{NUMBER_OF_UNKNOWNS}},0.0);
296 {% endif %}
297}
298{% endif %}
299
300{% if EIGENVALUES_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
301#if defined(OpenMPGPUOffloading)
302#pragma omp declare target
303#endif
304double {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::maxEigenvalue(
305 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
306 const tarch::la::Vector<Dimensions,double>& faceCentre,
307 double t,
308 double dt,
309 int normal,
310 Offloadable
311) {
312 {{EIGENVALUES_IMPLEMENTATION}};
313}
314#if defined(OpenMPGPUOffloading)
315#pragma omp end declare target
316#endif
317{% endif %}
318
319{% if FLUX_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
320#if defined(OpenMPGPUOffloading)
321#pragma omp declare target
322#endif
323void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::flux(
324 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
325 const tarch::la::Vector<Dimensions,double>& faceCentre,
326 double t,
327 double dt,
328 int normal,
329 double * __restrict__ F, // F[{{NUMBER_OF_UNKNOWNS}}]
330 Offloadable
331) {
332 {% if FLUX_IMPLEMENTATION=="<none>" %}
333 abort();
334 {% else %}
335 {{FLUX_IMPLEMENTATION}}
336 {% endif %}
337}
338#if defined(OpenMPGPUOffloading)
339#pragma omp end declare target
340#endif
341{% endif %}
342
343{% if NCP_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
344#if defined(OpenMPGPUOffloading)
345#pragma omp declare target
346#endif
347void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::nonconservativeProduct(
348 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
349 const double * __restrict__ deltaQ, // deltaQ[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}]
350 const tarch::la::Vector<Dimensions,double>& faceCentre,
351 double t,
352 double dt,
353 int normal,
354 double * __restrict__ BTimesDeltaQ, // BTimesDeltaQ[{{NUMBER_OF_UNKNOWNS}}]
355 Offloadable
356) {
357 {% if NCP_IMPLEMENTATION=="<none>" %}
358 abort();
359 {% else %}
360 {{NCP_IMPLEMENTATION}}
361 {% endif %}
362}
363#if defined(OpenMPGPUOffloading)
364#pragma omp end declare target
365#endif
366{% endif %}
367
368{% if SOURCE_TERM_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
369#if defined(OpenMPGPUOffloading)
370#pragma omp declare target
371#endif
372void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::sourceTerm(
373 const double * __restrict__ Q,
374 const tarch::la::Vector<Dimensions,double>& volumeCentre,
375 double t,
376 double dt,
377 double * __restrict__ S,
378 Offloadable
379) {
380 {% if SOURCE_TERM_IMPLEMENTATION=="<none>" %}
381 abort();
382 {% else %}
383 {{SOURCE_TERM_IMPLEMENTATION}}
384 {% endif %}
385}
386#if defined(OpenMPGPUOffloading)
387#pragma omp end declare target
388#endif
389{% endif %}
390
391{% if POINT_SOURCE_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
392#if defined(OpenMPGPUOffloading)
393#pragma omp declare target
394#endif
395void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::pointSource(
396 const double * __restrict__ Q,
397 const tarch::la::Vector<Dimensions,double>& x,
398 double t,
399 double dt,
400 double * __restrict__ S,
401 Offloadable
402) {
403 assertionMsg( false, "point sources not implemented yet" );
404 {% if POINT_SOURCE_IMPLEMENTATION=="<none>" %}
405 abort();
406 {% else %}
407 {{POINT_SOURCE_IMPLEMENTATION}}
408 {% endif %}
409}
410#if defined(OpenMPGPUOffloading)
411#pragma omp end declare target
412#endif
413{% endif %}
414""",
415 undefined=jinja2.DebugUndefined,
416 )
417
418 d = {}
419 d["FLUX_IMPLEMENTATION"] = flux_implementation
420 d["NCP_IMPLEMENTATION"] = ncp_implementation
421 d["EIGENVALUES_IMPLEMENTATION"] = eigenvalues_implementation
422 d["SOURCE_TERM_IMPLEMENTATION"] = source_term_implementation
423 d["POINT_SOURCE_IMPLEMENTATION"] = point_source_implementation
424 d["STATELESS_PDE_TERMS"] = pde_terms_without_state
425 return Template.render(**d)
426
427
428def create_solver_declarations(
429 flux_implementation,
430 ncp_implementation,
431 eigenvalues_implementation,
432 source_term_implementation,
433 point_source_implementation,
434 pde_terms_without_state,
435):
436 Template = jinja2.Template(
437 """
438 public:
439 {% if SOURCE_TERM_IMPLEMENTATION=="<user-defined>" %}
440 virtual void sourceTerm(
441 const double * __restrict__ Q,
442 const tarch::la::Vector<Dimensions,double>& volumeCentre,
443 double t,
444 double dt,
445 double * __restrict__ S
446 ) override;
447 {% endif %}
448
449 {% if POINT_SOURCE_IMPLEMENTATION=="<user-defined>" %}
450 #error point sources not defined yet
451 virtual void pointSource(
452 const double * __restrict__ Q,
453 const tarch::la::Vector<Dimensions,double>& x,
454 double t,
455 double dt,
456 double * __restrict__ S
457 ) override;
458 {% endif %}
459
460 {% if EIGENVALUES_IMPLEMENTATION=="<user-defined>" %}
461 /**
462 * Determine max eigenvalue over Jacobian in a given point with solution values
463 * (states) Q. All parameters are in.
464 *
465 * @return Max eigenvalue. Result has to be positive, so we are actually speaking
466 * about the maximum absolute eigenvalue.
467 */
468 virtual double maxEigenvalue(
469 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
470 const tarch::la::Vector<Dimensions,double>& faceCentre,
471 double t,
472 double dt,
473 int normal
474 ) override;
475 {% endif %}
476
477 {% if FLUX_IMPLEMENTATION=="<user-defined>" %}
478 virtual void flux(
479 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
480 const tarch::la::Vector<Dimensions,double>& faceCentre,
481 double t,
482 double dt,
483 int normal,
484 double * __restrict__ F // F[{{NUMBER_OF_UNKNOWNS}}]
485 ) override;
486 {% endif %}
487
488 {% if NCP_IMPLEMENTATION=="<user-defined>" %}
489 virtual void nonconservativeProduct(
490 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
491 const double * __restrict__ deltaQ, // deltaQ[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}]
492 const tarch::la::Vector<Dimensions,double>& faceCentre,
493 double t,
494 double dt,
495 int normal,
496 double * __restrict__ BTimesDeltaQ // BTimesDeltaQ[{{NUMBER_OF_UNKNOWNS}}]
497 ) override;
498 {% endif %}
499
500 {% if STATELESS_PDE_TERMS and SOURCE_TERM_IMPLEMENTATION=="<user-defined>" %}
501 /**
502 * To obtain the best performance, I recommend to man inline command to
503 * this signature and to copy the implementation into the header. So it would
504 * read
505 *
506 * static inline void sourceTerm( ... ) {
507 * code here
508 * }
509 *
510 * The GPU offloading requires static functions. As we cannot overload the
511 * original (virtual) function with a static alternative, we do the
512 * TBB trick and overload by adding an additional enum. It has no semantics
513 * but helps the compiler to distinguish the different function variants.
514 */
515 #if defined(OpenMPGPUOffloading)
516 #pragma omp declare target
517 #endif
518 static void sourceTerm(
519 const double * __restrict__ Q,
520 const tarch::la::Vector<Dimensions,double>& volumeCentre,
521 const tarch::la::Vector<Dimensions,double>& volumeH,
522 double t,
523 double dt,
524 double * __restrict__ S,
525 Offloadable
526 );
527 #if defined(OpenMPGPUOffloading)
528 #pragma omp end declare target
529 #endif
530 {% endif %}
531
532 {% if STATELESS_PDE_TERMS and POINT_SOURCE_IMPLEMENTATION=="<user-defined>" %}
533 /**
534 * To obtain the best performance, I recommend to man inline command to
535 * this signature and to copy the implementation into the header. So it would
536 * read
537 *
538 * static inline void pointSource( ... ) {
539 * code here
540 * }
541 *
542 * The GPU offloading requires static functions. As we cannot overload the
543 * original (virtual) function with a static alternative, we do the
544 * TBB trick and overload by adding an additional enum. It has no semantics
545 * but helps the compiler to distinguish the different function variants.
546 */
547 #error point sources not implemented yet
548 #if defined(OpenMPGPUOffloading)
549 #pragma omp declare target
550 #endif
551 static void pointSource(
552 const double * __restrict__ Q,
553 const tarch::la::Vector<Dimensions,double>& volumeCentre,
554 const tarch::la::Vector<Dimensions,double>& volumeH,
555 double t,
556 double dt,
557 double * __restrict__ S,
558 Offloadable
559 );
560 #if defined(OpenMPGPUOffloading)
561 #pragma omp end declare target
562 #endif
563 {% endif %}
564
565 {% if EIGENVALUES_IMPLEMENTATION=="<user-defined>" and STATELESS_PDE_TERMS %}
566 /**
567 * To obtain the best performance, I recommend to man inline command to
568 * this signature and to copy the implementation into the header. So it would
569 * read
570 *
571 * static inline double maxEigenvalue( ... ) {
572 * code here
573 * }
574 *
575 * The GPU offloading requires static functions. As we cannot overload the
576 * original (virtual) function with a static alternative, we do the
577 * TBB trick and overload by adding an additional enum. It has no semantics
578 * but helps the compiler to distinguish the different function variants.
579 */
580 #if defined(OpenMPGPUOffloading)
581 #pragma omp declare target
582 #endif
583 static double maxEigenvalue(
584 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
585 const tarch::la::Vector<Dimensions,double>& faceCentre,
586 double t,
587 double dt,
588 int normal,
589 Offloadable
590 );
591 #if defined(OpenMPGPUOffloading)
592 #pragma omp end declare target
593 #endif
594 {% endif %}
595
596 {% if STATELESS_PDE_TERMS and FLUX_IMPLEMENTATION=="<user-defined>" %}
597 /**
598 * To obtain the best performance, I recommend to man inline command to
599 * this signature and to copy the implementation into the header. So it would
600 * read
601 *
602 * static inline void flux( ... ) {
603 * code here
604 * }
605 *
606 * The GPU offloading requires static functions. As we cannot overload the
607 * original (virtual) function with a static alternative, we do the
608 * TBB trick and overload by adding an additional enum. It has no semantics
609 * but helps the compiler to distinguish the different function variants.
610 */
611 #if defined(OpenMPGPUOffloading)
612 #pragma omp declare target
613 #endif
614 static void flux(
615 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
616 const tarch::la::Vector<Dimensions,double>& faceCentre,
617 double t,
618 double dt,
619 int normal,
620 double * __restrict__ F, // F[{{NUMBER_OF_UNKNOWNS}}]
621 Offloadable
622 );
623 #if defined(OpenMPGPUOffloading)
624 #pragma omp end declare target
625 #endif
626 {% endif %}
627
628 {% if STATELESS_PDE_TERMS and NCP_IMPLEMENTATION=="<user-defined>" %}
629 /**
630 * To obtain the best performance, I recommend to man inline command to
631 * this signature and to copy the implementation into the header. So it would
632 * read
633 *
634 * static inline void nonconservativeProduct( ... ) {
635 * code here
636 * }
637 *
638 * The GPU offloading requires static functions. As we cannot overload the
639 * original (virtual) function with a static alternative, we do the
640 * TBB trick and overload by adding an additional enum. It has no semantics
641 * but helps the compiler to distinguish the different function variants.
642 */
643 #if defined(OpenMPGPUOffloading)
644 #pragma omp declare target
645 #endif
646 static void nonconservativeProduct(
647 const double * __restrict__ Q, // Q[5+0],
648 const double * __restrict__ deltaQ, // deltaQ[5+0]
649 const tarch::la::Vector<Dimensions,double>& x,
650 double t,
651 double dt,
652 int normal,
653 double * __restrict__ BTimesDeltaQ, // BTimesDeltaQ[5]
654 Offloadable
655 );
656 #if defined(OpenMPGPUOffloading)
657 #pragma omp end declare target
658 #endif
659 {% endif %}
660""",
661 undefined=jinja2.DebugUndefined,
662 )
663 d = {}
664 d["FLUX_IMPLEMENTATION"] = flux_implementation
665 d["NCP_IMPLEMENTATION"] = ncp_implementation
666 d["EIGENVALUES_IMPLEMENTATION"] = eigenvalues_implementation
667 d["SOURCE_TERM_IMPLEMENTATION"] = source_term_implementation
668 d["POINT_SOURCE_IMPLEMENTATION"] = point_source_implementation
669 d["STATELESS_PDE_TERMS"] = pde_terms_without_state
670 return Template.render(**d)
671
672
673def create_solver_definitions(
674 flux_implementation,
675 ncp_implementation,
676 eigenvalues_implementation,
677 source_term_implementation,
678 point_source_implementation,
679 pde_terms_without_state,
680):
681 Template = jinja2.Template(
682 """
683{% if EIGENVALUES_IMPLEMENTATION=="<user-defined>" %}
684double {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::maxEigenvalue(
685 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
686 const tarch::la::Vector<Dimensions,double>& faceCentre,
687 double t,
688 double dt,
689 int normal
690) {
691 logTraceInWith4Arguments( "maxEigenvalue(...)", x, t, dt, normal );
692 // @todo implement
693 logTraceOut( "maxEigenvalue(...)" );
694}
695{% endif %}
696
697{% if FLUX_IMPLEMENTATION=="<user-defined>" %}
698void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::flux(
699 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
700 const tarch::la::Vector<Dimensions,double>& faceCentre,
701 double t,
702 double dt,
703 int normal,
704 double * __restrict__ F // F[{{NUMBER_OF_UNKNOWNS}}]
705) {
706 logTraceInWith4Arguments( "flux(...)", faceCentre, t, dt, normal );
707 // @todo implement
708 logTraceOut( "flux(...)" );
709}
710{% endif %}
711
712{% if NCP_IMPLEMENTATION=="<user-defined>" %}
713void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::nonconservativeProduct(
714 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
715 const double * __restrict__ deltaQ, // deltaQ[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}]
716 const tarch::la::Vector<Dimensions,double>& faceCentre,
717 double t,
718 double dt,
719 int normal,
720 double * __restrict__ BTimesDeltaQ // BTimesDeltaQ[{{NUMBER_OF_UNKNOWNS}}]
721) {
722 logTraceInWith4Arguments( "nonconservativeProduct(...)", faceCentre, t, dt, normal );
723 // @todo implement
724 logTraceOut( "nonconservativeProduct(...)" );
725}
726{% endif %}
727
728{% if SOURCE_TERM_IMPLEMENTATION=="<user-defined>" %}
729void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::sourceTerm(
730 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}]
731 const tarch::la::Vector<Dimensions,double>& volumeCentre,
732 double t,
733 double dt,
734 double * __restrict__ S // S[{{NUMBER_OF_UNKNOWNS}}]
735) {
736 logTraceInWith3Arguments( "sourceTerm(...)", volumeCentre, t, dt );
737
738 // @todo implement and ensure that all entries of S are properly set
739 for (int i=0; i<NumberOfUnknowns; i++) {
740 S[i] = 0.0;
741 }
742
743 logTraceOut( "sourceTerm(...)" );
744}
745{% endif %}
746
747{% if EIGENVALUES_IMPLEMENTATION=="<user-defined>" and STATELESS_PDE_TERMS %}
748 #if defined(OpenMPGPUOffloading)
749 #pragma omp declare target
750 #endif
751double {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::maxEigenvalue(
752 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
753 const tarch::la::Vector<Dimensions,double>& faceCentre,
754 double t,
755 double dt,
756 int normal,
757 Offloadable
758) {
759 // @todo implement
760}
761 #if defined(OpenMPGPUOffloading)
762 #pragma omp end declare target
763 #endif
764{% endif %}
765
766{% if FLUX_IMPLEMENTATION=="<user-defined>" and STATELESS_PDE_TERMS %}
767 #if defined(OpenMPGPUOffloading)
768 #pragma omp declare target
769 #endif
770void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::flux(
771 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
772 const tarch::la::Vector<Dimensions,double>& faceCentre,
773 double t,
774 double dt,
775 int normal,
776 double * __restrict__ F,
777 Offloadable
778) {
779 // @todo implement
780}
781 #if defined(OpenMPGPUOffloading)
782 #pragma omp end declare target
783 #endif
784{% endif %}
785
786{% if NCP_IMPLEMENTATION=="<user-defined>" and STATELESS_PDE_TERMS %}
787 #if defined(OpenMPGPUOffloading)
788 #pragma omp declare target
789 #endif
790void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::nonconservativeProduct(
791 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
792 const double * __restrict__ deltaQ, // deltaQ[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}]
793 const tarch::la::Vector<Dimensions,double>& faceCentre,
794 double t,
795 double dt,
796 int normal,
797 double * __restrict__ BTimesDeltaQ,
798 Offloadable
799) {
800 // @todo implement
801}
802 #if defined(OpenMPGPUOffloading)
803 #pragma omp end declare target
804 #endif
805{% endif %}
806
807{% if SOURCE_TERM_IMPLEMENTATION=="<user-defined>" and STATELESS_PDE_TERMS %}
808 #if defined(OpenMPGPUOffloading)
809 #pragma omp declare target
810 #endif
811void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::sourceTerm(
812 const double * __restrict__ Q,
813 const tarch::la::Vector<Dimensions,double>& volumeCentre,
814 double t,
815 double dt,
816 double * __restrict__ S,
817 Offloadable
818) {
819 // @todo implement
820}
821 #if defined(OpenMPGPUOffloading)
822 #pragma omp end declare target
823 #endif
824{% endif %}
825
826{% if POINT_SOURCE_IMPLEMENTATION=="<user-defined>" and STATELESS_PDE_TERMS %}
827 #if defined(OpenMPGPUOffloading)
828 #pragma omp declare target
829 #endif
830#error point sources not implemented yet
831void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::sourceTerm(
832 const double * __restrict__ Q,
833 const tarch::la::Vector<Dimensions,double>& volumeCentre,
834 double t,
835 double dt,
836 double * __restrict__ S,
837 Offloadable
838) {
839 // @todo implement
840}
841 #if defined(OpenMPGPUOffloading)
842 #pragma omp end declare target
843 #endif
844{% endif %}
845""",
846 undefined=jinja2.DebugUndefined,
847 )
848 d = {}
849 d["FLUX_IMPLEMENTATION"] = flux_implementation
850 d["NCP_IMPLEMENTATION"] = ncp_implementation
851 d["EIGENVALUES_IMPLEMENTATION"] = eigenvalues_implementation
852 d["SOURCE_TERM_IMPLEMENTATION"] = source_term_implementation
853 d["POINT_SOURCE_IMPLEMENTATION"] = point_source_implementation
854 d["STATELESS_PDE_TERMS"] = pde_terms_without_state
855 return Template.render(**d)
856
857
858def get_face_merge_implementation(patch_overlap):
859 assert False, "I don't think this one is used anymore. @todo Remove"
860
861 d = {
862 "OVERLAP": patch_overlap.dim[0]
863 / 2, # assumed to always be two, but present here anyway
864 "DOFS_PER_AXIS": patch_overlap.dim[1],
865 "UNKNOWNS": patch_overlap.no_of_unknowns,
866 }
867 template = """
868
869 #if Dimensions==2
870 constexpr int nodesPerFace = {DOFS_PER_AXIS};
871 #else
872 constexpr int nodesPerFace = {DOFS_PER_AXIS}*{DOFS_PER_AXIS};
873 #endif
874
875 constexpr int strideQ = {UNKNOWNS};
876
877 const int faceDimension = marker.getSelectedFaceNumber()%Dimensions;
878 //if the outer normal is positive, the normal points to the right so the face is on the right
879 const int faceLR = (marker.outerNormal()[faceDimension]>0.0 ? 1 : 0);
880
881 for(int i=0; i<nodesPerFace; i++){{
882 for(int j=0; j<strideQ ; j++){{
883 value[(2*i+faceLR)*strideQ+j] = neighbour.value[(2*i+faceLR)*strideQ+j];
884 }}
885 }}
886
887"""
888
889 return template.format(**d)
890
891
892def create_start_time_step_implementation_for_fixed_time_stepping(
893 normalised_time_step_size,
894):
895 """
896 We roll over all reduced data after the last time step, and we
897 plot status info in the first step.
898 """
899 return (
900 """
901 if (
902 tarch::mpi::Rank::getInstance().isGlobalMaster()
903 and
904 _maxCellH>0.0
905 and
906 isFirstGridSweepOfTimeStep()
907 ) {
908 logInfo( "step()", "Solver {{SOLVER_NAME}}:" );
909 logInfo( "step()", "t = " << _minTimeStampThisTimeStep );
910 logInfo( "step()", "dt = " << getTimeStepSize() );
911 logInfo( "step()", "h_{min} = " << _minCellHThisTimeStep );
912 logInfo( "step()", "h_{max} = " << _maxCellHThisTimeStep );
913 }
914 _timeStepSize = """
915 + str(normalised_time_step_size)
916 + """;
917"""
918 )
919
920
921def create_volumetric_solver_call(polynomial_basis, solver_variant):
922 """
923 A set of default volumetric kernel calls. You might want to use different
924 solver calls in your code depending on the context.
925
926 ## Difference to Finite Volume (FV) implementation
927
928 In the FV code base, I need the implementation routines as argument: As we work
929 with FV, the generic FV class does not know which implementations are around. If
930 you work with a domain-specific Rusanov solver, e.g., then there is no such
931 thing as an ncp.
932
933 For the DG schemes, things are different: Every DG solver in ExaHyPE 2 in principle
934 supports the ncp. So the base Python class already instantiates the corresponding
935 dictionary entries, and we can use them straightaway.
936
937 ## Arguments
938
939 solver_variant: SolverVariant
940 Picks different implementation variants
941 """
942 if (
943 isinstance(polynomial_basis, GaussLegendreBasis)
944 and solver_variant == SolverVariant.WithVirtualFunctions
945 ):
946 return """cellIntegral_patchwise_in_situ_GaussLegendre_functors(
947 cellData,
948 {{ORDER}},
949 {{NUMBER_OF_UNKNOWNS}},
950 {{NUMBER_OF_AUXILIARY_VARIABLES}},
951 [&](
952 const double * __restrict__ Q,
953 const tarch::la::Vector<Dimensions,double>& x,
954 double t,
955 double dt,
956 int normal,
957 double * __restrict__ F
958 )->void {
959 {% if FLUX_IMPLEMENTATION!="<none>" %}
960 repositories::{{SOLVER_INSTANCE}}.flux(Q,x,t,dt,normal,F);
961 {% endif %}
962 }, //flux
963 [&](
964 const double * __restrict__ Q,
965 const double * __restrict__ deltaQ,
966 const tarch::la::Vector<Dimensions,double>& faceCentre,
967 double t,
968 double dt,
969 int normal,
970 double * __restrict__ F
971 ) -> void {
972 {% if NCP_IMPLEMENTATION!="<none>" %}
973 repositories::{{SOLVER_INSTANCE}}.nonconservativeProduct(Q, deltaQ, faceCentre, t, dt, normal, F);
974 {% endif %}
975 }, //ncp
976 [&](
977 const double * __restrict__ Q,
978 const tarch::la::Vector<Dimensions,double>& x,
979 double t,
980 double dt,
981 double * __restrict__ S
982 ) -> void {
983 {% if SOURCE_TERM_IMPLEMENTATION!="<none>" %}
984 repositories::{{SOLVER_INSTANCE}}.sourceTerm(Q,x,t,dt,S);
985 {% endif %}
986 }, //source
987 [&](
988 const double * __restrict__ Q,
989 const tarch::la::Vector<Dimensions,double>& cellCentre,
990 const tarch::la::Vector<Dimensions,double>& cellH,
991 double t,
992 double dt
993 ) -> std::vector<::exahype2::dg::PointSource> {
994 {% if POINT_SOURCES_IMPLEMENTATION!="<none>" %}
995 return repositories::{{SOLVER_INSTANCE}}.pointSources(Q,cellCentre,cellH,t,dt);
996 {% else %}
997 return std::vector<::exahype2::dg::PointSource>();
998 {% endif %}
999 }, //source
1000 repositories::{{SOLVER_INSTANCE}}.QuadraturePoints1d,
1001 repositories::{{SOLVER_INSTANCE}}.MassMatrixDiagonal1d,
1002 repositories::{{SOLVER_INSTANCE}}.StiffnessOperator1d,
1003 repositories::{{SOLVER_INSTANCE}}.DerivativeOperator1d,
1004 {{ "true" if FLUX_IMPLEMENTATION!="<none>" else "false" }}, //useFlux
1005 {{ "true" if NCP_IMPLEMENTATION!="<none>" else "false" }}, //useNCP
1006 {{ "true" if SOURCE_TERM_IMPLEMENTATION!="<none>" else "false" }}, //useSource
1007 {{ "true" if POINT_SOURCES_IMPLEMENTATION!="<none>" else "false" }} //usePointSource
1008 );
1009"""
1010 elif (
1011 isinstance(polynomial_basis, GaussLegendreBasis)
1012 and solver_variant == SolverVariant.Stateless
1013 ):
1014 return """cellIntegral_patchwise_in_situ_GaussLegendre<{{SOLVER_NAME}},{{ORDER}},{{NUMBER_OF_UNKNOWNS}},{{NUMBER_OF_AUXILIARY_VARIABLES}}>(
1015 cellData,
1016 {{ "true" if FLUX_IMPLEMENTATION!="<none>" else "false" }}, //useFlux
1017 {{ "true" if NCP_IMPLEMENTATION!="<none>" else "false" }}, //useNCP
1018 {{ "true" if SOURCE_TERM_IMPLEMENTATION!="<none>" else "false" }}, //useSource
1019 {{ "true" if POINT_SOURCES_IMPLEMENTATION!="<none>" else "false" }} //usePointSource
1020 );
1021"""
1022 elif (
1023 isinstance(polynomial_basis, GaussLegendreBasis)
1024 and solver_variant == SolverVariant.Accelerator
1025 ):
1026 return """cellIntegral_patchwise_in_situ_GaussLegendre<{{SOLVER_NAME}},{{ORDER}},{{NUMBER_OF_UNKNOWNS}},{{NUMBER_OF_AUXILIARY_VARIABLES}}>(
1027 targetDevice,
1028 cellData,
1029 {{ "true" if FLUX_IMPLEMENTATION!="<none>" else "false" }}, //useFlux
1030 {{ "true" if NCP_IMPLEMENTATION!="<none>" else "false" }}, //useNCP
1031 {{ "true" if SOURCE_TERM_IMPLEMENTATION!="<none>" else "false" }}, //useSource
1032 {{ "true" if POINT_SOURCES_IMPLEMENTATION!="<none>" else "false" }} //usePointSource
1033 );
1034"""
1035 else:
1036 return "#solver variant not implemented yet"
1037
1038
1039def create_add_solver_contributions_call(polynomial_basis):
1040 if isinstance(polynomial_basis, GaussLegendreBasis):
1041 return """integrateOverRiemannSolutionsAndAddToVolume_GaussLegendre(
1042 #if Dimensions==2
1043 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(0).value, //left
1044 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(2).value, //right
1045 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(1).value, //down
1046 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(3).value, //up
1047 #else
1048 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(0).value, //left
1049 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(3).value, //right
1050 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(1).value, //down
1051 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(4).value, //up
1052 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(2).value, //front
1053 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(5).value, //back
1054 #endif
1055 {{ORDER}},
1056 {{NUMBER_OF_UNKNOWNS}},
1057 {{NUMBER_OF_AUXILIARY_VARIABLES}},
1058 marker.h(),
1059 repositories::{{SOLVER_INSTANCE}}.BasisFunctionValuesLeft1d,
1060 repositories::{{SOLVER_INSTANCE}}.MassMatrixDiagonal1d,
1061 dQdt
1062 );
1063"""
1064 else:
1065 return "#error Not supported yet"
1066
1067
1068def create_multiply_with_inverted_mass_matrix_call(polynomial_basis):
1069 if isinstance(polynomial_basis, GaussLegendreBasis):
1070 return """multiplyWithInvertedMassMatrix_GaussLegendre(
1071 cellData,
1072 {{ORDER}},
1073 {{NUMBER_OF_UNKNOWNS}},
1074 {{NUMBER_OF_AUXILIARY_VARIABLES}},
1075 repositories::{{SOLVER_INSTANCE}}.MassMatrixDiagonal1d
1076 );
1077"""
1078 else:
1079 return "#error Not supported yet"
exahype2.solvers.rkdg.kernels.create_abstract_solver_definitions
create_abstract_solver_definitions(flux_implementation, ncp_implementation, eigenvalues_implementation, source_term_implementation, point_source_implementation, pde_terms_without_state)
Definition kernels.py:225
exahype2.solvers.rkdg.kernels.create_volumetric_solver_call
create_volumetric_solver_call(polynomial_basis, solver_variant)
A set of default volumetric kernel calls.
Definition kernels.py:921
exahype2.solvers.rkdg.kernels.create_multiply_with_inverted_mass_matrix_call
create_multiply_with_inverted_mass_matrix_call(polynomial_basis)
Definition kernels.py:1068
exahype2.solvers.rkdg.kernels.create_solver_declarations
create_solver_declarations(flux_implementation, ncp_implementation, eigenvalues_implementation, source_term_implementation, point_source_implementation, pde_terms_without_state)
Definition kernels.py:435
exahype2.solvers.rkdg.kernels.create_abstract_solver_declarations
create_abstract_solver_declarations(flux_implementation, ncp_implementation, eigenvalues_implementation, source_term_implementation, point_source_implementation, pde_terms_without_state)
Definition kernels.py:26
exahype2.solvers.rkdg.kernels.create_solver_definitions
create_solver_definitions(flux_implementation, ncp_implementation, eigenvalues_implementation, source_term_implementation, point_source_implementation, pde_terms_without_state)
Definition kernels.py:680
exahype2.solvers.rkdg.kernels.create_start_time_step_implementation_for_fixed_time_stepping
create_start_time_step_implementation_for_fixed_time_stepping(normalised_time_step_size)
We roll over all reduced data after the last time step, and we plot status info in the first step.
Definition kernels.py:894
exahype2.solvers.rkdg.kernels.create_add_solver_contributions_call
create_add_solver_contributions_call(polynomial_basis)
Definition kernels.py:1039
exahype2.solvers.rkdg.kernels.get_face_merge_implementation
get_face_merge_implementation(patch_overlap)
Definition kernels.py:858