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