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