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 os
9import exahype2
10import exahype2.kerneldsl as DSL
11
12from exahype2.solvers.LagrangeBasisWithDiagonalMassMatrix import GaussLegendreBasis
13
14class SolverVariant(Enum):
15 WithVirtualFunctions = 0
16 Stateless = 1
17 Accelerator = 2
18
19def clear_dg_kernel_file():
20 if os.path.exists("kernels"):
21 file = open("kernels/dg.h", "w")
22 file.close()
23
24def create_cell_integral_kernel_definitions(namespace):
25 if not os.path.exists("kernels"):
26 os.makedirs("kernels")
27
28 file = open("kernels/dg.h", "a")
29 file.write("""#if Dimensions == 2
30#include "dg_cellIntegral2d.h"
31#elif Dimensions == 3
32#include "dg_cellIntegral3d.h"
33#endif
34""")
35
36 create_cell_integral_kernel_definitions_2d(namespace)
37 create_cell_integral_kernel_definitions_3d(namespace)
38
39def create_cell_integral_kernel_definitions_2d(namespace):
40 template_parameters = [DSL.SyntaxTree.Argument("DGOrder", DSL.SyntaxTree.TInteger()),
41 DSL.SyntaxTree.Argument("NumberOfUnknowns", DSL.SyntaxTree.TInteger()),
42 DSL.SyntaxTree.Argument("NumberOfAuxiliaryVariables", DSL.SyntaxTree.TInteger()),
43 DSL.SyntaxTree.Argument("EvaluateFlux", DSL.SyntaxTree.TBoolean()),
44 DSL.SyntaxTree.Argument("EvaluateNonconservativeProduct", DSL.SyntaxTree.TBoolean()),
45 DSL.SyntaxTree.Argument("EvaluateSource", DSL.SyntaxTree.TBoolean()),
46 DSL.SyntaxTree.Argument("EvaluatePointSources", DSL.SyntaxTree.TBoolean())]
47
48 functor_arguments = [DSL.SyntaxTree.Argument("flux", DSL.SyntaxTree.TCustom("const Flux&")),
49 DSL.SyntaxTree.Argument("nonconservativeProduct", DSL.SyntaxTree.TCustom("const NonConservativeProduct&")),
50 DSL.SyntaxTree.Argument("sourceTerm", DSL.SyntaxTree.TCustom("const Source&")),
51 DSL.SyntaxTree.Argument("pointSource", DSL.SyntaxTree.TCustom("const PointSources&")),
52 DSL.SyntaxTree.Argument("getQuadraturePoint", DSL.SyntaxTree.TCustom("const QuadraturePoint&")),
53 DSL.SyntaxTree.Argument("QuadraturePoints1d", DSL.SyntaxTree.TCustom("const double*")),
54 DSL.SyntaxTree.Argument("MassMatrixDiagonal1d", DSL.SyntaxTree.TCustom("const double*")),
55 DSL.SyntaxTree.Argument("StiffnessOperator1d", DSL.SyntaxTree.TCustom("const double*")),
56 DSL.SyntaxTree.Argument("DerivativeOperator1d", DSL.SyntaxTree.TCustom("const double*"))]
57
58 dg_stateless_tree = DSL.Parser().parse(exahype2.solvers.rkdg.cellIntegral2d, template_parameters, functor_arguments, ["exahype2", "dg", namespace], stateless=True)
59 dg_stateless_kernel = dg_stateless_tree.print_cpp()
60 dg_stateless_call_with_measurement = dg_stateless_tree.print_definition_with_timer()
61 dg_stateless_kernel_declaration = dg_stateless_tree.print_declaration()
62 dg_stateless_call_with_measurement_declaration = dg_stateless_tree.print_declaration_with_timer()
63
64 dg_stateless_omp_accelerator_tree = DSL.Parser().parse(exahype2.solvers.rkdg.cellIntegral2d, template_parameters, functor_arguments, ["exahype2", "dg", "omp", namespace], stateless=True, use_accelerator=True)
65 dg_stateless_omp_accelerator_kernel = dg_stateless_omp_accelerator_tree.print_omp()
66 dg_stateless_omp_accelerator_call_with_measurement = dg_stateless_omp_accelerator_tree.print_definition_with_timer()
67 dg_stateless_omp_accelerator_kernel_declaration = dg_stateless_omp_accelerator_tree.print_declaration()
68 dg_stateless_omp_accelerator_call_with_measurement_declaration = dg_stateless_omp_accelerator_tree.print_declaration_with_timer()
69
70 dg_tree = DSL.Parser().parse(exahype2.solvers.rkdg.cellIntegral2d, template_parameters, functor_arguments, ["exahype2", "dg", namespace], stateless=False)
71 dg_kernel = dg_tree.print_cpp()
72 dg_call_with_measurement = dg_tree.print_definition_with_timer()
73 dg_kernel_declaration = dg_tree.print_declaration()
74 dg_call_with_measurement_declaration = dg_tree.print_declaration_with_timer()
75
76 if not os.path.exists("kernels"):
77 os.makedirs("kernels")
78
79 file = open("kernels/dg_cellIntegral2d.h", "w")
80 file.write(f"""#pragma once
81#include "exahype2/CellData.h"
82#include "exahype2/VolumeIndex.h"
83#include "exahype2/dg/DGUtils.h"
84#include "exahype2/dg/Functors.h"
85#include "peano4/utils/Loop.h"
86#include "tarch/timing/Measurement.h"
87#include "tarch/timing/Watch.h"
88#include <fstream>
89
90{dg_stateless_kernel_declaration}
91{dg_stateless_call_with_measurement_declaration}
92
93{dg_stateless_omp_accelerator_kernel_declaration}
94{dg_stateless_omp_accelerator_call_with_measurement_declaration}
95
96{dg_kernel_declaration}
97{dg_call_with_measurement_declaration}
98
99#include "dg_cellIntegral2d.cpph"
100""")
101 file.close()
102
103 file = open("kernels/dg_cellIntegral2d.cpph", "w")
104 file.write("#include \"exahype2/dg/DGUtils.h\"\n")
105 file.write(dg_stateless_kernel)
106 file.write(dg_stateless_call_with_measurement)
107 file.write(dg_stateless_omp_accelerator_kernel)
108 file.write(dg_stateless_omp_accelerator_call_with_measurement)
109 file.write(dg_kernel)
110 file.write(dg_call_with_measurement)
111 file.close()
112
113def create_cell_integral_kernel_definitions_3d(namespace):
114 template_parameters = [DSL.SyntaxTree.Argument("DGOrder", DSL.SyntaxTree.TInteger()),
115 DSL.SyntaxTree.Argument("NumberOfUnknowns", DSL.SyntaxTree.TInteger()),
116 DSL.SyntaxTree.Argument("NumberOfAuxiliaryVariables", DSL.SyntaxTree.TInteger()),
117 DSL.SyntaxTree.Argument("EvaluateFlux", DSL.SyntaxTree.TBoolean()),
118 DSL.SyntaxTree.Argument("EvaluateNonconservativeProduct", DSL.SyntaxTree.TBoolean()),
119 DSL.SyntaxTree.Argument("EvaluateSource", DSL.SyntaxTree.TBoolean()),
120 DSL.SyntaxTree.Argument("EvaluatePointSources", DSL.SyntaxTree.TBoolean())]
121
122 functor_arguments = [DSL.SyntaxTree.Argument("flux", DSL.SyntaxTree.TCustom("const Flux&")),
123 DSL.SyntaxTree.Argument("nonconservativeProduct", DSL.SyntaxTree.TCustom("const NonConservativeProduct&")),
124 DSL.SyntaxTree.Argument("sourceTerm", DSL.SyntaxTree.TCustom("const Source&")),
125 DSL.SyntaxTree.Argument("pointSource", DSL.SyntaxTree.TCustom("const PointSources&")),
126 DSL.SyntaxTree.Argument("getQuadraturePoint", DSL.SyntaxTree.TCustom("const QuadraturePoint&")),
127 DSL.SyntaxTree.Argument("QuadraturePoints1d", DSL.SyntaxTree.TCustom("const double*")),
128 DSL.SyntaxTree.Argument("MassMatrixDiagonal1d", DSL.SyntaxTree.TCustom("const double*")),
129 DSL.SyntaxTree.Argument("StiffnessOperator1d", DSL.SyntaxTree.TCustom("const double*")),
130 DSL.SyntaxTree.Argument("DerivativeOperator1d", DSL.SyntaxTree.TCustom("const double*"))]
131
132 dg_stateless_tree = DSL.Parser().parse(exahype2.solvers.rkdg.cellIntegral3d, template_parameters, functor_arguments, ["exahype2", "dg", namespace], stateless=True)
133 dg_stateless_kernel = dg_stateless_tree.print_cpp()
134 dg_stateless_call_with_measurement = dg_stateless_tree.print_definition_with_timer()
135 dg_stateless_kernel_declaration = dg_stateless_tree.print_declaration()
136 dg_stateless_call_with_measurement_declaration = dg_stateless_tree.print_declaration_with_timer()
137
138 dg_stateless_omp_accelerator_tree = DSL.Parser().parse(exahype2.solvers.rkdg.cellIntegral3d, template_parameters, functor_arguments, ["exahype2", "dg", "omp", namespace], stateless=True, use_accelerator=True)
139 dg_stateless_omp_accelerator_kernel = dg_stateless_omp_accelerator_tree.print_cpp()
140 dg_stateless_omp_accelerator_call_with_measurement = dg_stateless_omp_accelerator_tree.print_definition_with_timer()
141 dg_stateless_omp_accelerator_kernel_declaration = dg_stateless_omp_accelerator_tree.print_declaration()
142 dg_stateless_omp_accelerator_call_with_measurement_declaration = dg_stateless_omp_accelerator_tree.print_declaration_with_timer()
143
144 dg_tree = DSL.Parser().parse(exahype2.solvers.rkdg.cellIntegral3d, template_parameters, functor_arguments, ["exahype2", "dg", namespace], stateless=False)
145 dg_kernel = dg_tree.print_cpp()
146 dg_call_with_measurement = dg_tree.print_definition_with_timer()
147 dg_kernel_declaration = dg_tree.print_declaration()
148 dg_call_with_measurement_declaration = dg_tree.print_declaration_with_timer()
149
150 if not os.path.exists("kernels"):
151 os.makedirs("kernels")
152
153 file = open("kernels/dg_cellIntegral3d.h", "w")
154 file.write(f"""#pragma once
155#include "exahype2/CellData.h"
156#include "exahype2/VolumeIndex.h"
157#include "exahype2/dg/DGUtils.h"
158#include "exahype2/dg/Functors.h"
159#include "peano4/utils/Loop.h"
160#include "tarch/timing/Measurement.h"
161#include "tarch/timing/Watch.h"
162#include <fstream>
163
164{dg_stateless_kernel_declaration}
165{dg_stateless_call_with_measurement_declaration}
166
167{dg_stateless_omp_accelerator_kernel_declaration}
168{dg_stateless_omp_accelerator_call_with_measurement_declaration}
169
170{dg_kernel_declaration}
171{dg_call_with_measurement_declaration}
172
173#include "dg_cellIntegral3d.cpph"
174""")
175 file.close()
176
177 file = open("kernels/dg_cellIntegral3d.cpph", "w")
178 file.write("#include \"exahype2/dg/DGUtils.h\"\n")
179 file.write(dg_stateless_kernel)
180 file.write(dg_stateless_call_with_measurement)
181 file.write(dg_stateless_omp_accelerator_kernel)
182 file.write(dg_stateless_omp_accelerator_call_with_measurement)
183 file.write(dg_kernel)
184 file.write(dg_call_with_measurement)
185 file.close()
186
187def create_abstract_solver_declarations(
188 flux_implementation,
189 ncp_implementation,
190 eigenvalues_implementation,
191 source_term_implementation,
192 point_source_implementation,
193 pde_terms_without_state,
194):
195 Template = jinja2.Template(
196 """
197 public:
198 {% if EIGENVALUES_IMPLEMENTATION=="<none>" or EIGENVALUES_IMPLEMENTATION=="<empty>" %}
199 #error Eigenvalue implementation cannot be none or empty
200 {% endif %}
201
202 {% if EIGENVALUES_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
203 /**
204 * Depending on the implementation, this variant might be slow as it
205 * lacks an inline define. Also, if you don't want to use ipo aggressively,
206 * it might be clever to put the implementation into the header.
207 */
208 #if defined(GPUOffloadingOMP)
209 #pragma omp declare target
210 #endif
211 static void maxEigenvalue(
212 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
213 const tarch::la::Vector<Dimensions,double>& faceCentre,
214 double t,
215 double dt,
216 int normal,
217 double* maxEigenvalue,
218 Offloadable
219 );
220 #if defined(GPUOffloadingOMP)
221 #pragma omp end declare target
222 #endif
223 {% endif %}
224
225 /**
226 * Determine max eigenvalue over Jacobian in a given point with solution values
227 * (states) Q. All parameters are in.
228 *
229 * @return Max eigenvalue. Result has to be positive, so we are actually speaking
230 * about the maximum absolute eigenvalue.
231 */
232 virtual void maxEigenvalue(
233 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
234 const tarch::la::Vector<Dimensions,double>& faceCentre,
235 double t,
236 double dt,
237 int normal,
238 double* maxEigenvalue
239 ) {% if EIGENVALUES_IMPLEMENTATION=="<user-defined>" %}= 0{% else %} final{% endif %};
240
241 {% if FLUX_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
242 /**
243 * Depending on the implementation, this variant might be slow as it
244 * lacks an inline define. Also, if you don't want to use ipo aggressively,
245 * it might be clever to put the implementation into the header.
246 */
247 #if defined(GPUOffloadingOMP)
248 #pragma omp declare target
249 #endif
250 static void flux(
251 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
252 const tarch::la::Vector<Dimensions,double>& faceCentre,
253 double t,
254 double dt,
255 int normal,
256 double * __restrict__ F, // F[{{NUMBER_OF_UNKNOWNS}}]
257 Offloadable
258 );
259 #if defined(GPUOffloadingOMP)
260 #pragma omp end declare target
261 #endif
262 {% endif %}
263
264 {% if FLUX_IMPLEMENTATION!="<none>" %}
265 virtual void flux(
266 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
267 const tarch::la::Vector<Dimensions,double>& faceCentre,
268 double t,
269 double dt,
270 int normal,
271 double * __restrict__ F // F[{{NUMBER_OF_UNKNOWNS}}]
272 ) {% if FLUX_IMPLEMENTATION=="<user-defined>" %}=0{% else %} final {% endif %};
273 {% endif %}
274
275 {% if NCP_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
276 /**
277 * Depending on the implementation, this variant might be slow as it
278 * lacks an inline define. Also, if you don't want to use ipo aggressively,
279 * it might be clever to put the implementation into the header.
280 */
281 #if defined(GPUOffloadingOMP)
282 #pragma omp declare target
283 #endif
284 static void nonconservativeProduct(
285 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
286 const double * __restrict__ deltaQ, // deltaQ[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
287 const tarch::la::Vector<Dimensions,double>& faceCentre,
288 double t,
289 double dt,
290 int normal,
291 double * __restrict__ BTimesDeltaQ, // BTimesDeltaQ[{{NUMBER_OF_UNKNOWNS}}]
292 Offloadable
293 );
294 #if defined(GPUOffloadingOMP)
295 #pragma omp end declare target
296 #endif
297 {% endif %}
298
299 {% if NCP_IMPLEMENTATION!="<none>" %}
300 virtual void nonconservativeProduct(
301 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
302 const double * __restrict__ deltaQ, // [{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}]
303 const tarch::la::Vector<Dimensions,double>& faceCentre,
304 double t,
305 double dt,
306 int normal,
307 double * __restrict__ BTimesDeltaQ // BTimesDeltaQ[{{NUMBER_OF_UNKNOWNS}}]
308 ) {% if NCP_IMPLEMENTATION=="<user-defined>" %}=0{% endif %};
309 {% endif %}
310
311 {% if SOURCE_TERM_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
312 /**
313 * Depending on the implementation, this variant might be slow as it
314 * lacks an inline define. Also, if you don't want to use ipo aggressively,
315 * it might be clever to put the implementation into the header.
316 */
317 #if defined(GPUOffloadingOMP)
318 #pragma omp declare target
319 #endif
320 static void sourceTerm(
321 const double * __restrict__ Q,
322 const tarch::la::Vector<Dimensions,double>& volumeCentre,
323 double t,
324 double dt,
325 double * __restrict__ S,
326 Offloadable
327 );
328 #if defined(GPUOffloadingOMP)
329 #pragma omp end declare target
330 #endif
331 {% endif %}
332
333 {% if SOURCE_TERM_IMPLEMENTATION!="<none>" %}
334 virtual void sourceTerm(
335 const double * __restrict__ Q,
336 const tarch::la::Vector<Dimensions,double>& volumeCentre,
337 double t,
338 double dt,
339 double * __restrict__ S
340 ) {% if SOURCE_TERM_IMPLEMENTATION=="<user-defined>" %}= 0{% else %} final {% endif %};
341 {% endif %}
342
343 {% if POINT_SOURCE_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
344 /**
345 * Depending on the implementation, this variant might be slow as it
346 * lacks an inline define. Also, if you don't want to use ipo aggressively,
347 * it might be clever to put the implementation into the header.
348 */
349 #if defined(GPUOffloadingOMP)
350 #pragma omp declare target
351 #endif
352 void pointSource(
353 const double * __restrict__ Q,
354 const tarch::la::Vector<Dimensions,double>& x,
355 double t,
356 double dt,
357 double * __restrict__ S,
358 Offloadable
359 );
360 #if defined(GPUOffloadingOMP)
361 #pragma omp end declare target
362 #endif
363 {% endif %}
364
365 {% if POINT_SOURCE_IMPLEMENTATION!="<none>" %}
366 virtual void pointSource(
367 const double * __restrict__ Q,
368 const tarch::la::Vector<Dimensions,double>& x,
369 double t,
370 double dt,
371 double * __restrict__ S
372 ) {% if SOURCE_TERM_IMPLEMENTATION=="<user-defined>" %}= 0{% else %} final {% endif %};
373 {% endif %}
374
375 #if defined(GPUOffloadingOMP)
376 #pragma omp declare target
377 #endif
378 static void getQuadraturePoint(
379 const tarch::la::Vector<Dimensions,int>& index,
380 const tarch::la::Vector<Dimensions,double>& cellCentre,
381 const tarch::la::Vector<Dimensions,double>& cellSize,
382 int polynomialOrder,
383 const double* __restrict__ quadraturePoints,
384 tarch::la::Vector<Dimensions,double>* __restrict__ quadraturePositions,
385 Offloadable
386 );
387 #if defined(GPUOffloadingOMP)
388 #pragma omp end declare target
389 #endif
390
391 #if defined(GPUOffloadingOMP)
392 #pragma omp declare target
393 #endif
394 static void getQuadraturePointInFace(
395 const tarch::la::Vector<Dimensions,int>& index,
396 const tarch::la::Vector<Dimensions,double>& faceCentre,
397 const tarch::la::Vector<Dimensions,double>& cellSize,
398 int polynomialOrder,
399 int faceOrientation,
400 const double* __restrict__ quadraturePoints,
401 tarch::la::Vector<Dimensions,double>* __restrict__ quadraturePositions,
402 Offloadable
403 );
404 #if defined(GPUOffloadingOMP)
405 #pragma omp end declare target
406 #endif
407
408 void getQuadraturePoint(
409 const tarch::la::Vector<Dimensions,int>& index,
410 const tarch::la::Vector<Dimensions,double>& cellCentre,
411 const tarch::la::Vector<Dimensions,double>& cellSize,
412 int polynomialOrder,
413 const double* __restrict__ quadraturePoints,
414 tarch::la::Vector<Dimensions,double>* __restrict__ quadraturePositions
415 );
416
417 void getQuadraturePointInFace(
418 const tarch::la::Vector<Dimensions,int>& index,
419 const tarch::la::Vector<Dimensions,double>& faceCentre,
420 const tarch::la::Vector<Dimensions,double>& cellSize,
421 int polynomialOrder,
422 int faceOrientation,
423 const double* __restrict__ quadraturePoints,
424 tarch::la::Vector<Dimensions,double>* __restrict__ quadraturePositions
425 );
426""",
427 undefined=jinja2.DebugUndefined,
428 )
429
430 d = {}
431 d["FLUX_IMPLEMENTATION"] = flux_implementation
432 d["NCP_IMPLEMENTATION"] = ncp_implementation
433 d["EIGENVALUES_IMPLEMENTATION"] = eigenvalues_implementation
434 d["SOURCE_TERM_IMPLEMENTATION"] = source_term_implementation
435 d["POINT_SOURCE_IMPLEMENTATION"] = point_source_implementation
436 d["STATELESS_PDE_TERMS"] = pde_terms_without_state
437 return Template.render(**d)
438
439
440def create_abstract_solver_definitions(
441 flux_implementation,
442 ncp_implementation,
443 eigenvalues_implementation,
444 source_term_implementation,
445 point_source_implementation,
446 pde_terms_without_state,
447):
448 Template = jinja2.Template(
449 """
450{% if EIGENVALUES_IMPLEMENTATION!="<user-defined>" and EIGENVALUES_IMPLEMENTATION!="<none>" %}
451void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::maxEigenvalue(
452 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
453 const tarch::la::Vector<Dimensions,double>& faceCentre,
454 double t,
455 double dt,
456 int normal,
457 double* maxEigenvalue
458) {
459 {{EIGENVALUES_IMPLEMENTATION}}
460}
461{% endif %}
462
463{% if FLUX_IMPLEMENTATION!="<none>" and FLUX_IMPLEMENTATION!="<user-defined>" %}
464void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::flux(
465 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
466 const tarch::la::Vector<Dimensions,double>& faceCentre,
467 double t,
468 double dt,
469 int normal,
470 double * __restrict__ F // F[{{NUMBER_OF_UNKNOWNS}}]
471) {
472 {{FLUX_IMPLEMENTATION}}
473}
474{% endif %}
475
476{% if NCP_IMPLEMENTATION!="<none>" and NCP_IMPLEMENTATION!="<user-defined>" %}
477void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::nonconservativeProduct(
478 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
479 const double * __restrict__ deltaQ, // deltaQ[{{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__ BTimesDeltaQ // BTimesDeltaQ[{{NUMBER_OF_UNKNOWNS}}]
485) {
486 {{NCP_IMPLEMENTATION}}
487}
488{% endif %}
489
490{% if SOURCE_TERM_IMPLEMENTATION!="<user-defined>" and SOURCE_TERM_IMPLEMENTATION!="<none>" %}
491void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::sourceTerm(
492 const double * __restrict__ Q,
493 const tarch::la::Vector<Dimensions,double>& volumeCentre,
494 double t,
495 double dt,
496 double * __restrict__ S
497) {
498 {% if SOURCE_TERM_IMPLEMENTATION!="<empty>" %}
499 {{SOURCE_TERM_IMPLEMENTATION}}
500 {% else %}
501 std::fill_n(S,{{NUMBER_OF_UNKNOWNS}},0.0);
502 {% endif %}
503}
504{% endif %}
505
506{% if POINT_SOURCE_IMPLEMENTATION!="<user-defined>" and POINT_SOURCE_IMPLEMENTATION!="<none>" %}
507#error Point sources not yet implemented
508void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::pointSource(
509 const double * __restrict__ Q,
510 const tarch::la::Vector<Dimensions,double>& x,
511 double t,
512 double dt,
513 double * __restrict__ S
514) {
515 {% if POINT_SOURCE_IMPLEMENTATION!="<empty>" %}
516 {{POINT_SOURCE_IMPLEMENTATION}}
517 {% else %}
518 std::fill_n(S,{{NUMBER_OF_UNKNOWNS}},0.0);
519 {% endif %}
520}
521{% endif %}
522
523{% if EIGENVALUES_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
524#if defined(GPUOffloadingOMP)
525#pragma omp declare target
526#endif
527void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::maxEigenvalue(
528 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
529 const tarch::la::Vector<Dimensions,double>& faceCentre,
530 double t,
531 double dt,
532 int normal,
533 double* maxEigenvalue,
534 Offloadable
535) {
536 {{EIGENVALUES_IMPLEMENTATION}};
537}
538#if defined(GPUOffloadingOMP)
539#pragma omp end declare target
540#endif
541{% endif %}
542
543{% if FLUX_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
544#if defined(GPUOffloadingOMP)
545#pragma omp declare target
546#endif
547void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::flux(
548 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
549 const tarch::la::Vector<Dimensions,double>& faceCentre,
550 double t,
551 double dt,
552 int normal,
553 double * __restrict__ F, // F[{{NUMBER_OF_UNKNOWNS}}]
554 Offloadable
555) {
556 {% if FLUX_IMPLEMENTATION=="<none>" %}
557 abort();
558 {% else %}
559 {{FLUX_IMPLEMENTATION}}
560 {% endif %}
561}
562#if defined(GPUOffloadingOMP)
563#pragma omp end declare target
564#endif
565{% endif %}
566
567{% if NCP_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
568#if defined(GPUOffloadingOMP)
569#pragma omp declare target
570#endif
571void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::nonconservativeProduct(
572 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
573 const double * __restrict__ deltaQ, // deltaQ[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}]
574 const tarch::la::Vector<Dimensions,double>& faceCentre,
575 double t,
576 double dt,
577 int normal,
578 double * __restrict__ BTimesDeltaQ, // BTimesDeltaQ[{{NUMBER_OF_UNKNOWNS}}]
579 Offloadable
580) {
581 {% if NCP_IMPLEMENTATION=="<none>" %}
582 abort();
583 {% else %}
584 {{NCP_IMPLEMENTATION}}
585 {% endif %}
586}
587#if defined(GPUOffloadingOMP)
588#pragma omp end declare target
589#endif
590{% endif %}
591
592{% if SOURCE_TERM_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
593#if defined(GPUOffloadingOMP)
594#pragma omp declare target
595#endif
596void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::sourceTerm(
597 const double * __restrict__ Q,
598 const tarch::la::Vector<Dimensions,double>& volumeCentre,
599 double t,
600 double dt,
601 double * __restrict__ S,
602 Offloadable
603) {
604 {% if SOURCE_TERM_IMPLEMENTATION=="<none>" %}
605 abort();
606 {% else %}
607 {{SOURCE_TERM_IMPLEMENTATION}}
608 {% endif %}
609}
610#if defined(GPUOffloadingOMP)
611#pragma omp end declare target
612#endif
613{% endif %}
614
615{% if POINT_SOURCE_IMPLEMENTATION!="<user-defined>" and STATELESS_PDE_TERMS %}
616#if defined(GPUOffloadingOMP)
617#pragma omp declare target
618#endif
619void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::pointSource(
620 const double * __restrict__ Q,
621 const tarch::la::Vector<Dimensions,double>& x,
622 double t,
623 double dt,
624 double * __restrict__ S,
625 Offloadable
626) {
627 assertionMsg( false, "point sources not implemented yet" );
628 {% if POINT_SOURCE_IMPLEMENTATION=="<none>" %}
629 abort();
630 {% else %}
631 {{POINT_SOURCE_IMPLEMENTATION}}
632 {% endif %}
633}
634#if defined(GPUOffloadingOMP)
635#pragma omp end declare target
636#endif
637{% endif %}
638
639void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::getQuadraturePoint(
640 const tarch::la::Vector<Dimensions,int>& index,
641 const tarch::la::Vector<Dimensions,double>& cellCentre,
642 const tarch::la::Vector<Dimensions,double>& cellSize,
643 int polynomialOrder,
644 const double* __restrict__ quadraturePoints,
645 tarch::la::Vector<Dimensions,double>* __restrict__ quadraturePositions,
646 Offloadable
647){
648 quadraturePositions[0] = ::exahype2::dg::getQuadraturePoint(cellCentre, cellSize, index, polynomialOrder, quadraturePoints);
649}
650
651void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::getQuadraturePointInFace(
652 const tarch::la::Vector<Dimensions,int>& index,
653 const tarch::la::Vector<Dimensions,double>& faceCentre,
654 const tarch::la::Vector<Dimensions,double>& cellSize,
655 int polynomialOrder,
656 int faceOrientation,
657 const double* __restrict__ quadraturePoints,
658 tarch::la::Vector<Dimensions,double>* __restrict__ quadraturePositions,
659 Offloadable
660){
661 ::exahype2::dg::getQuadraturePointInFace(faceCentre, cellSize, index, polynomialOrder, faceOrientation, quadraturePoints, quadraturePositions);
662}
663
664void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::getQuadraturePoint(
665 const tarch::la::Vector<Dimensions,int>& index,
666 const tarch::la::Vector<Dimensions,double>& cellCentre,
667 const tarch::la::Vector<Dimensions,double>& cellSize,
668 int polynomialOrder,
669 const double* __restrict__ quadraturePoints,
670 tarch::la::Vector<Dimensions,double>* __restrict__ quadraturePositions
671){
672 quadraturePositions[0] = ::exahype2::dg::getQuadraturePoint(cellCentre, cellSize, index, polynomialOrder, quadraturePoints);
673}
674
675void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::getQuadraturePointInFace(
676 const tarch::la::Vector<Dimensions,int>& index,
677 const tarch::la::Vector<Dimensions,double>& faceCentre,
678 const tarch::la::Vector<Dimensions,double>& cellSize,
679 int polynomialOrder,
680 int faceOrientation,
681 const double* __restrict__ quadraturePoints,
682 tarch::la::Vector<Dimensions,double>* __restrict__ quadraturePositions
683){
684 ::exahype2::dg::getQuadraturePointInFace(faceCentre, cellSize, index, polynomialOrder, faceOrientation, quadraturePoints, quadraturePositions);
685}
686""",
687 undefined=jinja2.DebugUndefined,
688 )
689
690 d = {}
691 d["FLUX_IMPLEMENTATION"] = flux_implementation
692 d["NCP_IMPLEMENTATION"] = ncp_implementation
693 d["EIGENVALUES_IMPLEMENTATION"] = eigenvalues_implementation
694 d["SOURCE_TERM_IMPLEMENTATION"] = source_term_implementation
695 d["POINT_SOURCE_IMPLEMENTATION"] = point_source_implementation
696 d["STATELESS_PDE_TERMS"] = pde_terms_without_state
697 return Template.render(**d)
698
699
700def create_solver_declarations(
701 flux_implementation,
702 ncp_implementation,
703 eigenvalues_implementation,
704 source_term_implementation,
705 point_source_implementation,
706 pde_terms_without_state,
707):
708 Template = jinja2.Template(
709 """
710 public:
711 {% if SOURCE_TERM_IMPLEMENTATION=="<user-defined>" %}
712 virtual void sourceTerm(
713 const double * __restrict__ Q,
714 const tarch::la::Vector<Dimensions,double>& volumeCentre,
715 double t,
716 double dt,
717 double * __restrict__ S
718 ) override;
719 {% endif %}
720
721 {% if POINT_SOURCE_IMPLEMENTATION=="<user-defined>" %}
722 #error point sources not defined yet
723 virtual void pointSource(
724 const double * __restrict__ Q,
725 const tarch::la::Vector<Dimensions,double>& x,
726 double t,
727 double dt,
728 double * __restrict__ S
729 ) override;
730 {% endif %}
731
732 {% if EIGENVALUES_IMPLEMENTATION=="<user-defined>" %}
733 /**
734 * Determine max eigenvalue over Jacobian in a given point with solution values
735 * (states) Q. All parameters are in.
736 *
737 * @return Max eigenvalue. Result has to be positive, so we are actually speaking
738 * about the maximum absolute eigenvalue.
739 */
740 virtual void maxEigenvalue(
741 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
742 const tarch::la::Vector<Dimensions,double>& faceCentre,
743 double t,
744 double dt,
745 int normal,
746 double* maxEigenvalue
747 ) override;
748 {% endif %}
749
750 {% if FLUX_IMPLEMENTATION=="<user-defined>" %}
751 virtual void flux(
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 double * __restrict__ F // F[{{NUMBER_OF_UNKNOWNS}}]
758 ) override;
759 {% endif %}
760
761 {% if NCP_IMPLEMENTATION=="<user-defined>" %}
762 virtual void nonconservativeProduct(
763 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
764 const double * __restrict__ deltaQ, // deltaQ[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}]
765 const tarch::la::Vector<Dimensions,double>& faceCentre,
766 double t,
767 double dt,
768 int normal,
769 double * __restrict__ BTimesDeltaQ // BTimesDeltaQ[{{NUMBER_OF_UNKNOWNS}}]
770 ) override;
771 {% endif %}
772
773 {% if STATELESS_PDE_TERMS and SOURCE_TERM_IMPLEMENTATION=="<user-defined>" %}
774 /**
775 * To obtain the best performance, I recommend to man inline command to
776 * this signature and to copy the implementation into the header. So it would
777 * read
778 *
779 * static inline void sourceTerm( ... ) {
780 * code here
781 * }
782 *
783 * The GPU offloading requires static functions. As we cannot overload the
784 * original (virtual) function with a static alternative, we do the
785 * TBB trick and overload by adding an additional enum. It has no semantics
786 * but helps the compiler to distinguish the different function variants.
787 */
788 #if defined(GPUOffloadingOMP)
789 #pragma omp declare target
790 #endif
791 static void sourceTerm(
792 const double * __restrict__ Q,
793 const tarch::la::Vector<Dimensions,double>& volumeCentre,
794 const tarch::la::Vector<Dimensions,double>& volumeH,
795 double t,
796 double dt,
797 double * __restrict__ S,
798 Offloadable
799 );
800 #if defined(GPUOffloadingOMP)
801 #pragma omp end declare target
802 #endif
803 {% endif %}
804
805 {% if STATELESS_PDE_TERMS and POINT_SOURCE_IMPLEMENTATION=="<user-defined>" %}
806 /**
807 * To obtain the best performance, I recommend to man inline command to
808 * this signature and to copy the implementation into the header. So it would
809 * read
810 *
811 * static inline void pointSource( ... ) {
812 * code here
813 * }
814 *
815 * The GPU offloading requires static functions. As we cannot overload the
816 * original (virtual) function with a static alternative, we do the
817 * TBB trick and overload by adding an additional enum. It has no semantics
818 * but helps the compiler to distinguish the different function variants.
819 */
820 #error point sources not implemented yet
821 #if defined(GPUOffloadingOMP)
822 #pragma omp declare target
823 #endif
824 static void pointSource(
825 const double * __restrict__ Q,
826 const tarch::la::Vector<Dimensions,double>& volumeCentre,
827 const tarch::la::Vector<Dimensions,double>& volumeH,
828 double t,
829 double dt,
830 double * __restrict__ S,
831 Offloadable
832 );
833 #if defined(GPUOffloadingOMP)
834 #pragma omp end declare target
835 #endif
836 {% endif %}
837
838 {% if EIGENVALUES_IMPLEMENTATION=="<user-defined>" and STATELESS_PDE_TERMS %}
839 /**
840 * To obtain the best performance, I recommend to man inline command to
841 * this signature and to copy the implementation into the header. So it would
842 * read
843 *
844 * static inline void maxEigenvalue( ... ) {
845 * code here
846 * }
847 *
848 * The GPU offloading requires static functions. As we cannot overload the
849 * original (virtual) function with a static alternative, we do the
850 * TBB trick and overload by adding an additional enum. It has no semantics
851 * but helps the compiler to distinguish the different function variants.
852 */
853 #if defined(GPUOffloadingOMP)
854 #pragma omp declare target
855 #endif
856 static void maxEigenvalue(
857 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
858 const tarch::la::Vector<Dimensions,double>& faceCentre,
859 double t,
860 double dt,
861 int normal,
862 double* maxEigenvalue,
863 Offloadable
864 );
865 #if defined(GPUOffloadingOMP)
866 #pragma omp end declare target
867 #endif
868 {% endif %}
869
870 {% if STATELESS_PDE_TERMS and FLUX_IMPLEMENTATION=="<user-defined>" %}
871 /**
872 * To obtain the best performance, I recommend to man inline command to
873 * this signature and to copy the implementation into the header. So it would
874 * read
875 *
876 * static inline void flux( ... ) {
877 * code here
878 * }
879 *
880 * The GPU offloading requires static functions. As we cannot overload the
881 * original (virtual) function with a static alternative, we do the
882 * TBB trick and overload by adding an additional enum. It has no semantics
883 * but helps the compiler to distinguish the different function variants.
884 */
885 #if defined(GPUOffloadingOMP)
886 #pragma omp declare target
887 #endif
888 static void flux(
889 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
890 const tarch::la::Vector<Dimensions,double>& faceCentre,
891 double t,
892 double dt,
893 int normal,
894 double * __restrict__ F, // F[{{NUMBER_OF_UNKNOWNS}}]
895 Offloadable
896 );
897 #if defined(GPUOffloadingOMP)
898 #pragma omp end declare target
899 #endif
900 {% endif %}
901
902 {% if STATELESS_PDE_TERMS and NCP_IMPLEMENTATION=="<user-defined>" %}
903 /**
904 * To obtain the best performance, I recommend to man inline command to
905 * this signature and to copy the implementation into the header. So it would
906 * read
907 *
908 * static inline void nonconservativeProduct( ... ) {
909 * code here
910 * }
911 *
912 * The GPU offloading requires static functions. As we cannot overload the
913 * original (virtual) function with a static alternative, we do the
914 * TBB trick and overload by adding an additional enum. It has no semantics
915 * but helps the compiler to distinguish the different function variants.
916 */
917 #if defined(GPUOffloadingOMP)
918 #pragma omp declare target
919 #endif
920 static void nonconservativeProduct(
921 const double * __restrict__ Q, // Q[5+0],
922 const double * __restrict__ deltaQ, // deltaQ[5+0]
923 const tarch::la::Vector<Dimensions,double>& x,
924 double t,
925 double dt,
926 int normal,
927 double * __restrict__ BTimesDeltaQ, // BTimesDeltaQ[5]
928 Offloadable
929 );
930 #if defined(GPUOffloadingOMP)
931 #pragma omp end declare target
932 #endif
933 {% endif %}
934""",
935 undefined=jinja2.DebugUndefined,
936 )
937 d = {}
938 d["FLUX_IMPLEMENTATION"] = flux_implementation
939 d["NCP_IMPLEMENTATION"] = ncp_implementation
940 d["EIGENVALUES_IMPLEMENTATION"] = eigenvalues_implementation
941 d["SOURCE_TERM_IMPLEMENTATION"] = source_term_implementation
942 d["POINT_SOURCE_IMPLEMENTATION"] = point_source_implementation
943 d["STATELESS_PDE_TERMS"] = pde_terms_without_state
944 return Template.render(**d)
945
946
947def create_solver_definitions(
948 flux_implementation,
949 ncp_implementation,
950 eigenvalues_implementation,
951 source_term_implementation,
952 point_source_implementation,
953 pde_terms_without_state,
954):
955 Template = jinja2.Template(
956 """
957{% if EIGENVALUES_IMPLEMENTATION=="<user-defined>" %}
958void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::maxEigenvalue(
959 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
960 const tarch::la::Vector<Dimensions,double>& faceCentre,
961 double t,
962 double dt,
963 int normal,
964 double* maxEigenvalue
965) {
966 logTraceInWith4Arguments( "maxEigenvalue(...)", x, t, dt, normal );
967 // @todo implement
968 logTraceOut( "maxEigenvalue(...)" );
969}
970{% endif %}
971
972{% if FLUX_IMPLEMENTATION=="<user-defined>" %}
973void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::flux(
974 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
975 const tarch::la::Vector<Dimensions,double>& faceCentre,
976 double t,
977 double dt,
978 int normal,
979 double * __restrict__ F // F[{{NUMBER_OF_UNKNOWNS}}]
980) {
981 logTraceInWith4Arguments( "flux(...)", faceCentre, t, dt, normal );
982 // @todo implement
983 logTraceOut( "flux(...)" );
984}
985{% endif %}
986
987{% if NCP_IMPLEMENTATION=="<user-defined>" %}
988void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::nonconservativeProduct(
989 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
990 const double * __restrict__ deltaQ, // deltaQ[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}]
991 const tarch::la::Vector<Dimensions,double>& faceCentre,
992 double t,
993 double dt,
994 int normal,
995 double * __restrict__ BTimesDeltaQ // BTimesDeltaQ[{{NUMBER_OF_UNKNOWNS}}]
996) {
997 logTraceInWith4Arguments( "nonconservativeProduct(...)", faceCentre, t, dt, normal );
998 // @todo implement
999 logTraceOut( "nonconservativeProduct(...)" );
1000}
1001{% endif %}
1002
1003{% if SOURCE_TERM_IMPLEMENTATION=="<user-defined>" %}
1004void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::sourceTerm(
1005 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}]
1006 const tarch::la::Vector<Dimensions,double>& volumeCentre,
1007 double t,
1008 double dt,
1009 double * __restrict__ S // S[{{NUMBER_OF_UNKNOWNS}}]
1010) {
1011 logTraceInWith3Arguments( "sourceTerm(...)", volumeCentre, t, dt );
1012
1013 // @todo implement and ensure that all entries of S are properly set
1014 for (int i=0; i<NumberOfUnknowns; i++) {
1015 S[i] = 0.0;
1016 }
1017
1018 logTraceOut( "sourceTerm(...)" );
1019}
1020{% endif %}
1021
1022{% if EIGENVALUES_IMPLEMENTATION=="<user-defined>" and STATELESS_PDE_TERMS %}
1023 #if defined(GPUOffloadingOMP)
1024 #pragma omp declare target
1025 #endif
1026void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::maxEigenvalue(
1027 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
1028 const tarch::la::Vector<Dimensions,double>& faceCentre,
1029 double t,
1030 double dt,
1031 int normal,
1032 double* maxEigenvalue,
1033 Offloadable
1034) {
1035 // @todo implement
1036}
1037 #if defined(GPUOffloadingOMP)
1038 #pragma omp end declare target
1039 #endif
1040{% endif %}
1041
1042{% if FLUX_IMPLEMENTATION=="<user-defined>" and STATELESS_PDE_TERMS %}
1043 #if defined(GPUOffloadingOMP)
1044 #pragma omp declare target
1045 #endif
1046void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::flux(
1047 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
1048 const tarch::la::Vector<Dimensions,double>& faceCentre,
1049 double t,
1050 double dt,
1051 int normal,
1052 double * __restrict__ F,
1053 Offloadable
1054) {
1055 // @todo implement
1056}
1057 #if defined(GPUOffloadingOMP)
1058 #pragma omp end declare target
1059 #endif
1060{% endif %}
1061
1062{% if NCP_IMPLEMENTATION=="<user-defined>" and STATELESS_PDE_TERMS %}
1063 #if defined(GPUOffloadingOMP)
1064 #pragma omp declare target
1065 #endif
1066void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::nonconservativeProduct(
1067 const double * __restrict__ Q, // Q[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}],
1068 const double * __restrict__ deltaQ, // deltaQ[{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}]
1069 const tarch::la::Vector<Dimensions,double>& faceCentre,
1070 double t,
1071 double dt,
1072 int normal,
1073 double * __restrict__ BTimesDeltaQ,
1074 Offloadable
1075) {
1076 // @todo implement
1077}
1078 #if defined(GPUOffloadingOMP)
1079 #pragma omp end declare target
1080 #endif
1081{% endif %}
1082
1083{% if SOURCE_TERM_IMPLEMENTATION=="<user-defined>" and STATELESS_PDE_TERMS %}
1084 #if defined(GPUOffloadingOMP)
1085 #pragma omp declare target
1086 #endif
1087void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::sourceTerm(
1088 const double * __restrict__ Q,
1089 const tarch::la::Vector<Dimensions,double>& volumeCentre,
1090 double t,
1091 double dt,
1092 double * __restrict__ S,
1093 Offloadable
1094) {
1095 // @todo implement
1096}
1097 #if defined(GPUOffloadingOMP)
1098 #pragma omp end declare target
1099 #endif
1100{% endif %}
1101
1102{% if POINT_SOURCE_IMPLEMENTATION=="<user-defined>" and STATELESS_PDE_TERMS %}
1103 #if defined(GPUOffloadingOMP)
1104 #pragma omp declare target
1105 #endif
1106#error point sources not implemented yet
1107void {{FULL_QUALIFIED_NAMESPACE}}::{{CLASSNAME}}::sourceTerm(
1108 const double * __restrict__ Q,
1109 const tarch::la::Vector<Dimensions,double>& volumeCentre,
1110 double t,
1111 double dt,
1112 double * __restrict__ S,
1113 Offloadable
1114) {
1115 // @todo implement
1116}
1117 #if defined(GPUOffloadingOMP)
1118 #pragma omp end declare target
1119 #endif
1120{% endif %}
1121""",
1122 undefined=jinja2.DebugUndefined,
1123 )
1124 d = {}
1125 d["FLUX_IMPLEMENTATION"] = flux_implementation
1126 d["NCP_IMPLEMENTATION"] = ncp_implementation
1127 d["EIGENVALUES_IMPLEMENTATION"] = eigenvalues_implementation
1128 d["SOURCE_TERM_IMPLEMENTATION"] = source_term_implementation
1129 d["POINT_SOURCE_IMPLEMENTATION"] = point_source_implementation
1130 d["STATELESS_PDE_TERMS"] = pde_terms_without_state
1131 return Template.render(**d)
1132
1133
1134def get_face_merge_implementation(patch_overlap):
1135 assert False, "I don't think this one is used anymore. @todo Remove"
1136
1137 d = {
1138 "OVERLAP": patch_overlap.dim[0]
1139 / 2, # assumed to always be two, but present here anyway
1140 "DOFS_PER_AXIS": patch_overlap.dim[1],
1141 "UNKNOWNS": patch_overlap.no_of_unknowns,
1142 }
1143 template = """
1144
1145 #if Dimensions==2
1146 constexpr int nodesPerFace = {DOFS_PER_AXIS};
1147 #else
1148 constexpr int nodesPerFace = {DOFS_PER_AXIS}*{DOFS_PER_AXIS};
1149 #endif
1150
1151 constexpr int strideQ = {UNKNOWNS};
1152
1153 const int faceDimension = marker.getSelectedFaceNumber()%Dimensions;
1154 //if the outer normal is positive, the normal points to the right so the face is on the right
1155 const int faceLR = (marker.outerNormal()[faceDimension]>0.0 ? 1 : 0);
1156
1157 for(int i=0; i<nodesPerFace; i++){{
1158 for(int j=0; j<strideQ ; j++){{
1159 value[(2*i+faceLR)*strideQ+j] = neighbour.value[(2*i+faceLR)*strideQ+j];
1160 }}
1161 }}
1162
1163"""
1164
1165 return template.format(**d)
1166
1167
1168def create_start_time_step_implementation_for_fixed_time_stepping(
1169 normalised_time_step_size,
1170):
1171 """
1172 We roll over all reduced data after the last time step, and we
1173 plot status info in the first step.
1174 """
1175 return (
1176 """
1177 if (
1178 tarch::mpi::Rank::getInstance().isGlobalMaster()
1179 and
1180 _maxCellH>0.0
1181 and
1182 isFirstGridSweepOfTimeStep()
1183 ) {
1184 logInfo( "step()", "Solver {{SOLVER_NAME}}:" );
1185 logInfo( "step()", "t = " << _minTimeStampThisTimeStep );
1186 logInfo( "step()", "dt = " << getTimeStepSize() );
1187 logInfo( "step()", "h_{min} = " << _minCellHThisTimeStep );
1188 logInfo( "step()", "h_{max} = " << _maxCellHThisTimeStep );
1189 }
1190 _timeStepSize = """
1191 + str(normalised_time_step_size)
1192 + """;
1193"""
1194 )
1195
1196def create_volumetric_solver_call_dsl(polynomial_basis, solver_variant):
1197 """
1198 A set of default volumetric kernel calls. You might want to use different
1199 solver calls in your code depending on the context.
1200
1201 ## Difference to Finite Volume (FV) implementation
1202
1203 In the FV code base, I need the implementation routines as argument: As we work
1204 with FV, the generic FV class does not know which implementations are around. If
1205 you work with a domain-specific dg solver, e.g., then there is no such
1206 thing as an ncp.
1207
1208 For the DG schemes, things are different: Every DG solver in ExaHyPE 2 in principle
1209 supports the ncp. So the base Python class already instantiates the corresponding
1210 dictionary entries, and we can use them straightaway.
1211
1212 ## Arguments
1213
1214 solver_variant: SolverVariant
1215 Picks different implementation variants
1216 """
1217 if (
1218 isinstance(polynomial_basis, GaussLegendreBasis)
1219 and solver_variant == SolverVariant.WithVirtualFunctions
1220 ):
1221 return """cellIntegral<{{ORDER}},{{NUMBER_OF_UNKNOWNS}},{{NUMBER_OF_AUXILIARY_VARIABLES}},
1222 {{ "true" if FLUX_IMPLEMENTATION!="<none>" else "false" }}, //useFlux
1223 {{ "true" if NCP_IMPLEMENTATION!="<none>" else "false" }}, //useNCP
1224 {{ "true" if SOURCE_TERM_IMPLEMENTATION!="<none>" else "false" }}, //useSource
1225 {{ "true" if POINT_SOURCES_IMPLEMENTATION!="<none>" else "false" }} //usePointSource
1226 >(
1227 cellData,
1228 [&](
1229 const double * __restrict__ Q,
1230 const tarch::la::Vector<Dimensions,double>& x,
1231 double t,
1232 double dt,
1233 int normal,
1234 double * __restrict__ F
1235 )->void {
1236 {% if FLUX_IMPLEMENTATION!="<none>" %}
1237 repositories::{{SOLVER_INSTANCE}}.flux(Q,x,t,dt,normal,F);
1238 {% endif %}
1239 }, //flux
1240 [&](
1241 const double * __restrict__ Q,
1242 const double * __restrict__ deltaQ,
1243 const tarch::la::Vector<Dimensions,double>& faceCentre,
1244 double t,
1245 double dt,
1246 int normal,
1247 double * __restrict__ F
1248 ) -> void {
1249 {% if NCP_IMPLEMENTATION!="<none>" %}
1250 repositories::{{SOLVER_INSTANCE}}.nonconservativeProduct(Q, deltaQ, faceCentre, t, dt, normal, F);
1251 {% endif %}
1252 }, //ncp
1253 [&](
1254 const double * __restrict__ Q,
1255 const tarch::la::Vector<Dimensions,double>& x,
1256 double t,
1257 double dt,
1258 double * __restrict__ S
1259 ) -> void {
1260 {% if SOURCE_TERM_IMPLEMENTATION!="<none>" %}
1261 repositories::{{SOLVER_INSTANCE}}.sourceTerm(Q,x,t,dt,S);
1262 {% endif %}
1263 }, //source
1264 [&](
1265 const double * __restrict__ Q,
1266 const tarch::la::Vector<Dimensions,double>& cellCentre,
1267 const tarch::la::Vector<Dimensions,double>& cellH,
1268 double t,
1269 double dt
1270 ) -> std::vector<::exahype2::dg::PointSource> {
1271 {% if POINT_SOURCES_IMPLEMENTATION!="<none>" %}
1272 return repositories::{{SOLVER_INSTANCE}}.pointSources(Q,cellCentre,cellH,t,dt);
1273 {% else %}
1274 return std::vector<::exahype2::dg::PointSource>();
1275 {% endif %}
1276 }, //source
1277 [&](
1278 const tarch::la::Vector<Dimensions,int>& index,
1279 const tarch::la::Vector<Dimensions,double>& cellCentre,
1280 const tarch::la::Vector<Dimensions,double>& cellSize,
1281 int polynomialOrder,
1282 const double* __restrict__ quadraturePoints,
1283 tarch::la::Vector<Dimensions,double>* __restrict__ quadraturePositions
1284 ) -> std::vector<::exahype2::dg::PointSource> {
1285 repositories::{{SOLVER_INSTANCE}}.getQuadraturePoint(index, cellCentre, cellSize, polynomialOrder, quadraturePoints, quadraturePositions);
1286 }, //quadraturePoints
1287 repositories::{{SOLVER_INSTANCE}}.QuadraturePoints1d,
1288 repositories::{{SOLVER_INSTANCE}}.MassMatrixDiagonal1d,
1289 repositories::{{SOLVER_INSTANCE}}.StiffnessOperator1d,
1290 repositories::{{SOLVER_INSTANCE}}.DerivativeOperator1d
1291 );
1292"""
1293 elif (
1294 isinstance(polynomial_basis, GaussLegendreBasis)
1295 and solver_variant == SolverVariant.Stateless
1296 ):
1297 return """cellIntegralStateless<{{SOLVER_NAME}},{{ORDER}},{{NUMBER_OF_UNKNOWNS}},{{NUMBER_OF_AUXILIARY_VARIABLES}},
1298 {{ "true" if FLUX_IMPLEMENTATION!="<none>" else "false" }}, //useFlux
1299 {{ "true" if NCP_IMPLEMENTATION!="<none>" else "false" }}, //useNCP
1300 {{ "true" if SOURCE_TERM_IMPLEMENTATION!="<none>" else "false" }}, //useSource
1301 {{ "true" if POINT_SOURCES_IMPLEMENTATION!="<none>" else "false" }} //usePointSource
1302 >(
1303 cellData
1304 );
1305"""
1306 elif (
1307 isinstance(polynomial_basis, GaussLegendreBasis)
1308 and solver_variant == SolverVariant.Accelerator
1309 ):
1310 return """cellIntegralStateless<{{SOLVER_NAME}},{{ORDER}},{{NUMBER_OF_UNKNOWNS}},{{NUMBER_OF_AUXILIARY_VARIABLES}},
1311 {{ "true" if FLUX_IMPLEMENTATION!="<none>" else "false" }}, //useFlux
1312 {{ "true" if NCP_IMPLEMENTATION!="<none>" else "false" }}, //useNCP
1313 {{ "true" if SOURCE_TERM_IMPLEMENTATION!="<none>" else "false" }}, //useSource
1314 {{ "true" if POINT_SOURCES_IMPLEMENTATION!="<none>" else "false" }} //usePointSource
1315 >(
1316 gpuCellData.targetDevice,
1317 gpuCellData
1318 );
1319"""
1320 else:
1321 return "#solver variant not implemented yet"
1322
1323def create_volumetric_solver_call(polynomial_basis, solver_variant):
1324 """
1325 A set of default volumetric kernel calls. You might want to use different
1326 solver calls in your code depending on the context.
1327
1328 ## Difference to Finite Volume (FV) implementation
1329
1330 In the FV code base, I need the implementation routines as argument: As we work
1331 with FV, the generic FV class does not know which implementations are around. If
1332 you work with a domain-specific dg solver, e.g., then there is no such
1333 thing as an ncp.
1334
1335 For the DG schemes, things are different: Every DG solver in ExaHyPE 2 in principle
1336 supports the ncp. So the base Python class already instantiates the corresponding
1337 dictionary entries, and we can use them straightaway.
1338
1339 ## Arguments
1340
1341 solver_variant: SolverVariant
1342 Picks different implementation variants
1343 """
1344 if (
1345 isinstance(polynomial_basis, GaussLegendreBasis)
1346 and solver_variant == SolverVariant.WithVirtualFunctions
1347 ):
1348 return """cellIntegral_patchwise_in_situ_GaussLegendre_functors(
1349 cellData,
1350 {{ORDER}},
1351 {{NUMBER_OF_UNKNOWNS}},
1352 {{NUMBER_OF_AUXILIARY_VARIABLES}},
1353 [&](
1354 const double * __restrict__ Q,
1355 const tarch::la::Vector<Dimensions,double>& x,
1356 double t,
1357 double dt,
1358 int normal,
1359 double * __restrict__ F
1360 )->void {
1361 {% if FLUX_IMPLEMENTATION!="<none>" %}
1362 repositories::{{SOLVER_INSTANCE}}.flux(Q,x,t,dt,normal,F);
1363 {% endif %}
1364 }, //flux
1365 [&](
1366 const double * __restrict__ Q,
1367 const double * __restrict__ deltaQ,
1368 const tarch::la::Vector<Dimensions,double>& faceCentre,
1369 double t,
1370 double dt,
1371 int normal,
1372 double * __restrict__ F
1373 ) -> void {
1374 {% if NCP_IMPLEMENTATION!="<none>" %}
1375 repositories::{{SOLVER_INSTANCE}}.nonconservativeProduct(Q, deltaQ, faceCentre, t, dt, normal, F);
1376 {% endif %}
1377 }, //ncp
1378 [&](
1379 const double * __restrict__ Q,
1380 const tarch::la::Vector<Dimensions,double>& x,
1381 double t,
1382 double dt,
1383 double * __restrict__ S
1384 ) -> void {
1385 {% if SOURCE_TERM_IMPLEMENTATION!="<none>" %}
1386 repositories::{{SOLVER_INSTANCE}}.sourceTerm(Q,x,t,dt,S);
1387 {% endif %}
1388 }, //source
1389 [&](
1390 const double * __restrict__ Q,
1391 const tarch::la::Vector<Dimensions,double>& cellCentre,
1392 const tarch::la::Vector<Dimensions,double>& cellH,
1393 double t,
1394 double dt
1395 ) -> std::vector<::exahype2::dg::PointSource> {
1396 {% if POINT_SOURCES_IMPLEMENTATION!="<none>" %}
1397 return repositories::{{SOLVER_INSTANCE}}.pointSources(Q,cellCentre,cellH,t,dt);
1398 {% else %}
1399 return std::vector<::exahype2::dg::PointSource>();
1400 {% endif %}
1401 }, //source
1402 repositories::{{SOLVER_INSTANCE}}.QuadraturePoints1d,
1403 repositories::{{SOLVER_INSTANCE}}.MassMatrixDiagonal1d,
1404 repositories::{{SOLVER_INSTANCE}}.StiffnessOperator1d,
1405 repositories::{{SOLVER_INSTANCE}}.DerivativeOperator1d,
1406 {{ "true" if FLUX_IMPLEMENTATION!="<none>" else "false" }}, //useFlux
1407 {{ "true" if NCP_IMPLEMENTATION!="<none>" else "false" }}, //useNCP
1408 {{ "true" if SOURCE_TERM_IMPLEMENTATION!="<none>" else "false" }}, //useSource
1409 {{ "true" if POINT_SOURCES_IMPLEMENTATION!="<none>" else "false" }} //usePointSource
1410 );
1411"""
1412 elif (
1413 isinstance(polynomial_basis, GaussLegendreBasis)
1414 and solver_variant == SolverVariant.Stateless
1415 ):
1416 return """cellIntegral_patchwise_in_situ_GaussLegendre<{{SOLVER_NAME}},{{ORDER}},{{NUMBER_OF_UNKNOWNS}},{{NUMBER_OF_AUXILIARY_VARIABLES}}>(
1417 cellData,
1418 {{ "true" if FLUX_IMPLEMENTATION!="<none>" else "false" }}, //useFlux
1419 {{ "true" if NCP_IMPLEMENTATION!="<none>" else "false" }}, //useNCP
1420 {{ "true" if SOURCE_TERM_IMPLEMENTATION!="<none>" else "false" }}, //useSource
1421 {{ "true" if POINT_SOURCES_IMPLEMENTATION!="<none>" else "false" }} //usePointSource
1422 );
1423"""
1424 elif (
1425 isinstance(polynomial_basis, GaussLegendreBasis)
1426 and solver_variant == SolverVariant.Accelerator
1427 ):
1428 return """cellIntegral_patchwise_in_situ_GaussLegendre<{{SOLVER_NAME}},{{ORDER}},{{NUMBER_OF_UNKNOWNS}},{{NUMBER_OF_AUXILIARY_VARIABLES}}>(
1429 targetDevice,
1430 cellData,
1431 {{ "true" if FLUX_IMPLEMENTATION!="<none>" else "false" }}, //useFlux
1432 {{ "true" if NCP_IMPLEMENTATION!="<none>" else "false" }}, //useNCP
1433 {{ "true" if SOURCE_TERM_IMPLEMENTATION!="<none>" else "false" }}, //useSource
1434 {{ "true" if POINT_SOURCES_IMPLEMENTATION!="<none>" else "false" }} //usePointSource
1435 );
1436"""
1437 else:
1438 return "#solver variant not implemented yet"
1439
1440
1441def create_add_solver_contributions_call(polynomial_basis):
1442 if isinstance(polynomial_basis, GaussLegendreBasis):
1443 return """integrateOverRiemannSolutionsAndAddToVolume_GaussLegendre(
1444 #if Dimensions==2
1445 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(0).value, //left
1446 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(2).value, //right
1447 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(1).value, //down
1448 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(3).value, //up
1449 #else
1450 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(0).value, //left
1451 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(3).value, //right
1452 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(1).value, //down
1453 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(4).value, //up
1454 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(2).value, //front
1455 fineGridFaces{{UNKNOWN_IDENTIFIER}}RiemannSolution(5).value, //back
1456 #endif
1457 {{ORDER}},
1458 {{NUMBER_OF_UNKNOWNS}},
1459 {{NUMBER_OF_AUXILIARY_VARIABLES}},
1460 marker.h(),
1461 repositories::{{SOLVER_INSTANCE}}.BasisFunctionValuesLeft1d,
1462 repositories::{{SOLVER_INSTANCE}}.MassMatrixDiagonal1d,
1463 dQdt
1464 );
1465"""
1466 else:
1467 return "#error Not supported yet"
1468
1469
1470def create_multiply_with_inverted_mass_matrix_call(polynomial_basis):
1471 if isinstance(polynomial_basis, GaussLegendreBasis):
1472 return """multiplyWithInvertedMassMatrix_GaussLegendre(
1473 cellData,
1474 {{ORDER}},
1475 {{NUMBER_OF_UNKNOWNS}},
1476 {{NUMBER_OF_AUXILIARY_VARIABLES}},
1477 repositories::{{SOLVER_INSTANCE}}.MassMatrixDiagonal1d
1478 );
1479"""
1480 else:
1481 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:447
exahype2.solvers.rkdg.kernels.create_volumetric_solver_call_dsl
create_volumetric_solver_call_dsl(polynomial_basis, solver_variant)
A set of default volumetric kernel calls.
Definition kernels.py:1196
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:1323
exahype2.solvers.rkdg.kernels.create_cell_integral_kernel_definitions_3d
create_cell_integral_kernel_definitions_3d(namespace)
Definition kernels.py:113
exahype2.solvers.rkdg.kernels.create_cell_integral_kernel_definitions
create_cell_integral_kernel_definitions(namespace)
Definition kernels.py:24
exahype2.solvers.rkdg.kernels.create_multiply_with_inverted_mass_matrix_call
create_multiply_with_inverted_mass_matrix_call(polynomial_basis)
Definition kernels.py:1470
exahype2.solvers.rkdg.kernels.create_cell_integral_kernel_definitions_2d
create_cell_integral_kernel_definitions_2d(namespace)
Definition kernels.py:39
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:707
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:194
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:954
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:1170
exahype2.solvers.rkdg.kernels.create_add_solver_contributions_call
create_add_solver_contributions_call(polynomial_basis)
Definition kernels.py:1441
exahype2.solvers.rkdg.kernels.get_face_merge_implementation
get_face_merge_implementation(patch_overlap)
Definition kernels.py:1134