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 exahype2
4import exahype2.dslhype as DSL
5import os
6
7from exahype2.solvers.rkdg.actionsets.ProjectLinearCombinationOfEstimatesOntoFaces import FaceProjections
8from exahype2.solvers.rkdg.actionsets.ProjectLinearCombinationOfEstimatesOntoFaces import compute_number_of_face_projection_quantities
9
10def create_Riemann_solver_kernel_definitions():
11 if not os.path.exists("kernels"):
12 os.makedirs("kernels")
13
14 file = open("kernels/dg.h", "a")
15 file.write("""#if Dimensions == 2
16#include "dg_riemann2d.h"
17#elif Dimensions == 3
18#include "dg_riemann3d.h"
19#endif
20""")
21
22 create_Riemann_solver_kernel_definitions_2d()
23 create_Riemann_solver_kernel_definitions_3d()
24
25def create_Riemann_solver_kernel_definitions_2d():
26 template_parameters = [DSL.SyntaxTree.Argument("DGOrder", DSL.SyntaxTree.TInteger()),
27 DSL.SyntaxTree.Argument("NumberOfUnknowns", DSL.SyntaxTree.TInteger()),
28 DSL.SyntaxTree.Argument("NumberOfAuxiliaryVariables", DSL.SyntaxTree.TInteger()),
29 DSL.SyntaxTree.Argument("EvaluateFlux", DSL.SyntaxTree.TBoolean()),
30 DSL.SyntaxTree.Argument("EvaluateNonconservativeProduct", DSL.SyntaxTree.TBoolean())]
31
32 functor_arguments = [DSL.SyntaxTree.Argument("flux", DSL.SyntaxTree.TCustom("const Flux&")),
33 DSL.SyntaxTree.Argument("nonconservativeProduct", DSL.SyntaxTree.TCustom("const NonConservativeProduct&"))]
34
35 dg_riemann_x_stateless_tree = DSL.Parser().parse(exahype2.solvers.rkdg.rusanov.solveRiemannProblemX2D, template_parameters, functor_arguments, ["exahype2", "dg", "rusanov"], stateless=True, use_accelerator=False)
36 dg_riemann_x_stateless_kernel = dg_riemann_x_stateless_tree.print_cpp()
37 dg_riemann_x_stateless_call_with_measurement = dg_riemann_x_stateless_tree.print_definition_with_timer()
38 dg_riemann_x_stateless_kernel_declaration = dg_riemann_x_stateless_tree.print_declaration()
39 dg_riemann_x_stateless_call_with_measurement_declaration = dg_riemann_x_stateless_tree.print_declaration_with_timer()
40
41 dg_riemann_y_stateless_tree = DSL.Parser().parse(exahype2.solvers.rkdg.rusanov.solveRiemannProblemY2D, template_parameters, functor_arguments, ["exahype2", "dg", "rusanov"], stateless=True, use_accelerator=False)
42 dg_riemann_y_stateless_kernel = dg_riemann_y_stateless_tree.print_cpp()
43 dg_riemann_y_stateless_call_with_measurement = dg_riemann_y_stateless_tree.print_definition_with_timer()
44 dg_riemann_y_stateless_kernel_declaration = dg_riemann_y_stateless_tree.print_declaration()
45 dg_riemann_y_stateless_call_with_measurement_declaration = dg_riemann_y_stateless_tree.print_declaration_with_timer()
46
47 if not os.path.exists("kernels"):
48 os.makedirs("kernels")
49
50 file = open("kernels/dg_riemann2d.h", "w")
51 file.write(f"""#pragma once
52#include "exahype2/CellData.h"
53#include "exahype2/VolumeIndex.h"
54#include "exahype2/dg/DGUtils.h"
55#include "exahype2/dg/Functors.h"
56#include "peano4/utils/Loop.h"
57#include "tarch/timing/Measurement.h"
58#include "tarch/timing/Watch.h"
59#include <fstream>
60
61{dg_riemann_x_stateless_kernel_declaration}
62{dg_riemann_x_stateless_call_with_measurement_declaration}
63
64{dg_riemann_y_stateless_kernel_declaration}
65{dg_riemann_y_stateless_call_with_measurement_declaration}
66
67namespace exahype2::dg::rusanov {{
68template <class SolverType,int DGOrder,int NumberOfUnknowns,int NumberOfAuxiliaryVariables,bool EvaluateFlux,bool EvaluateNonconservativeProduct>
69 void solveRiemannProblemStateless(const tarch::la::Vector<Dimensions,double>& faceCentre, const tarch::la::Vector<Dimensions,double>& cellSize, double t, double dt, int faceNumber, const double* __restrict__ faceData, double* __restrict__ solution);
70}}
71
72#include "dg_riemann2d.cpph"
73""")
74 file.close()
75
76 file = open("kernels/dg_riemann2d.cpph", "w")
77 file.write(dg_riemann_x_stateless_kernel)
78 file.write(dg_riemann_x_stateless_call_with_measurement)
79 file.write(dg_riemann_y_stateless_kernel)
80 file.write(dg_riemann_y_stateless_call_with_measurement)
81 file.write("""
82namespace exahype2::dg::rusanov {
83 template <class SolverType,int DGOrder,int NumberOfUnknowns,int NumberOfAuxiliaryVariables,bool EvaluateFlux,bool EvaluateNonconservativeProduct>
84 void solveRiemannProblemStateless(const tarch::la::Vector<Dimensions,double>& faceCentre, const tarch::la::Vector<Dimensions,double>& cellSize, double t, double dt, int faceNumber,const double* __restrict__ faceData, double* __restrict__ solution){
85 switch(faceNumber % Dimensions) {
86 case 0:
87 solveRiemannProblemXStateless<
88 SolverType,
89 DGOrder,
90 NumberOfUnknowns,
91 NumberOfAuxiliaryVariables,
92 EvaluateFlux,
93 EvaluateNonconservativeProduct
94 >(
95 faceCentre,
96 cellSize,
97 t,
98 dt,
99 faceData,
100 solution
101 );
102 break;
103 case 1:
104 solveRiemannProblemYStateless<
105 SolverType,
106 DGOrder,
107 NumberOfUnknowns,
108 NumberOfAuxiliaryVariables,
109 EvaluateFlux,
110 EvaluateNonconservativeProduct
111 >(
112 faceCentre,
113 cellSize,
114 t,
115 dt,
116 faceData,
117 solution
118 );
119 break;
120 }
121 }
122}\n""")
123
124 file.close()
125
126def create_Riemann_solver_kernel_definitions_3d():
127 template_parameters = [DSL.SyntaxTree.Argument("DGOrder", DSL.SyntaxTree.TInteger()),
128 DSL.SyntaxTree.Argument("NumberOfUnknowns", DSL.SyntaxTree.TInteger()),
129 DSL.SyntaxTree.Argument("NumberOfAuxiliaryVariables", DSL.SyntaxTree.TInteger()),
130 DSL.SyntaxTree.Argument("EvaluateFlux", DSL.SyntaxTree.TBoolean()),
131 DSL.SyntaxTree.Argument("EvaluateNonconservativeProduct", DSL.SyntaxTree.TBoolean())]
132
133 functor_arguments = [DSL.SyntaxTree.Argument("flux", DSL.SyntaxTree.TCustom("const Flux&")),
134 DSL.SyntaxTree.Argument("nonconservativeProduct", DSL.SyntaxTree.TCustom("const NonConservativeProduct&"))]
135
136 dg_riemann_x_stateless_tree = DSL.Parser().parse(exahype2.solvers.rkdg.rusanov.solveRiemannProblemX3D, template_parameters, functor_arguments, ["exahype2", "dg", "rusanov"], stateless=True, use_accelerator=False)
137 dg_riemann_x_stateless_kernel = dg_riemann_x_stateless_tree.print_cpp()
138 dg_riemann_x_stateless_call_with_measurement = dg_riemann_x_stateless_tree.print_definition_with_timer()
139 dg_riemann_x_stateless_kernel_declaration = dg_riemann_x_stateless_tree.print_declaration()
140 dg_riemann_x_stateless_call_with_measurement_declaration = dg_riemann_x_stateless_tree.print_declaration_with_timer()
141
142 dg_riemann_y_stateless_tree = DSL.Parser().parse(exahype2.solvers.rkdg.rusanov.solveRiemannProblemY3D, template_parameters, functor_arguments, ["exahype2", "dg", "rusanov"], stateless=True, use_accelerator=False)
143 dg_riemann_y_stateless_kernel = dg_riemann_y_stateless_tree.print_cpp()
144 dg_riemann_y_stateless_call_with_measurement = dg_riemann_y_stateless_tree.print_definition_with_timer()
145 dg_riemann_y_stateless_kernel_declaration = dg_riemann_y_stateless_tree.print_declaration()
146 dg_riemann_y_stateless_call_with_measurement_declaration = dg_riemann_y_stateless_tree.print_declaration_with_timer()
147
148 dg_riemann_z_stateless_tree = DSL.Parser().parse(exahype2.solvers.rkdg.rusanov.solveRiemannProblemZ3D, template_parameters, functor_arguments, ["exahype2", "dg", "rusanov"], stateless=True, use_accelerator=False)
149 dg_riemann_z_stateless_kernel = dg_riemann_z_stateless_tree.print_cpp()
150 dg_riemann_z_stateless_call_with_measurement = dg_riemann_z_stateless_tree.print_definition_with_timer()
151 dg_riemann_z_stateless_kernel_declaration = dg_riemann_z_stateless_tree.print_declaration()
152 dg_riemann_z_stateless_call_with_measurement_declaration = dg_riemann_z_stateless_tree.print_declaration_with_timer()
153
154 if not os.path.exists("kernels"):
155 os.makedirs("kernels")
156
157 file = open("kernels/dg_riemann3d.h", "w")
158 file.write(f"""#pragma once
159#include "exahype2/CellData.h"
160#include "exahype2/VolumeIndex.h"
161#include "exahype2/dg/DGUtils.h"
162#include "exahype2/dg/Functors.h"
163#include "peano4/utils/Loop.h"
164#include "tarch/timing/Measurement.h"
165#include "tarch/timing/Watch.h"
166#include <fstream>
167
168{dg_riemann_x_stateless_kernel_declaration}
169{dg_riemann_x_stateless_call_with_measurement_declaration}
170
171{dg_riemann_y_stateless_kernel_declaration}
172{dg_riemann_y_stateless_call_with_measurement_declaration}
173
174{dg_riemann_z_stateless_kernel_declaration}
175{dg_riemann_z_stateless_call_with_measurement_declaration}
176
177namespace exahype2::dg::rusanov {{
178template <class SolverType,int DGOrder,int NumberOfUnknowns,int NumberOfAuxiliaryVariables,bool EvaluateFlux,bool EvaluateNonconservativeProduct>
179 void solveRiemannProblemStateless(const tarch::la::Vector<Dimensions,double>& faceCentre, const tarch::la::Vector<Dimensions,double>& cellSize, double t, double dt, int faceNumber, const double* __restrict__ faceData, double* __restrict__ solution);
180}}
181
182#include "dg_riemann3d.cpph"
183""")
184 file.close()
185
186 file = open("kernels/dg_riemann3d.cpph", "w")
187 file.write(dg_riemann_x_stateless_kernel)
188 file.write(dg_riemann_x_stateless_call_with_measurement)
189 file.write(dg_riemann_y_stateless_kernel)
190 file.write(dg_riemann_y_stateless_call_with_measurement)
191 file.write(dg_riemann_z_stateless_kernel)
192 file.write(dg_riemann_z_stateless_call_with_measurement)
193 file.write("""
194namespace exahype2::dg::rusanov {
195 template <class SolverType,int DGOrder,int NumberOfUnknowns,int NumberOfAuxiliaryVariables,bool EvaluateFlux,bool EvaluateNonconservativeProduct>
196 void solveRiemannProblemStateless(const tarch::la::Vector<Dimensions,double>& faceCentre, const tarch::la::Vector<Dimensions,double>& cellSize, double t, double dt, int faceNumber, const double* __restrict__ faceData, double* __restrict__ solution){
197 switch(faceNumber % Dimensions) {
198 case 0:
199 solveRiemannProblemXStateless<
200 SolverType,
201 DGOrder,
202 NumberOfUnknowns,
203 NumberOfAuxiliaryVariables,
204 EvaluateFlux,
205 EvaluateNonconservativeProduct
206 >(
207 faceCentre,
208 cellSize,
209 t,
210 dt,
211 faceData,
212 solution
213 );
214 break;
215 case 1:
216 solveRiemannProblemYStateless<
217 SolverType,
218 DGOrder,
219 NumberOfUnknowns,
220 NumberOfAuxiliaryVariables,
221 EvaluateFlux,
222 EvaluateNonconservativeProduct
223 >(
224 faceCentre,
225 cellSize,
226 t,
227 dt,
228 faceData,
229 solution
230 );
231 break;
232 case 2:
233 solveRiemannProblemZStateless<
234 SolverType,
235 DGOrder,
236 NumberOfUnknowns,
237 NumberOfAuxiliaryVariables,
238 EvaluateFlux,
239 EvaluateNonconservativeProduct
240 >(
241 faceCentre,
242 cellSize,
243 t,
244 dt,
245 faceData,
246 solution
247 );
248 break;
249 }
250 }
251 }\n""")
252
253 file.close()
254
255def create_Riemann_solver_call_dsl(polynomial_basis,
256 face_projections: FaceProjections):
257 if isinstance( polynomial_basis, exahype2.solvers.GaussLegendreBasis) and face_projections==FaceProjections.Solution:
258 return """solveRiemannProblemStateless<
259 {{SOLVER_NAME}},
260 {{ORDER}},
261 {{NUMBER_OF_UNKNOWNS}},
262 {{NUMBER_OF_AUXILIARY_VARIABLES}},
263 {{ "true" if FLUX_IMPLEMENTATION!="<none>" else "false" }}, //useFlux
264 {{ "true" if NCP_IMPLEMENTATION!="<none>" else "false" }} //useNCP
265 >(
266 marker.x(),
267 marker.h(),
268 timeStamp,
269 timeStepSize,
270 marker.getSelectedFaceNumber(),
271 fineGridFace{{UNKNOWN_IDENTIFIER}}EstimateProjection.value,
272 fineGridFace{{UNKNOWN_IDENTIFIER}}RiemannSolution.value
273 );
274"""
275
276def create_Riemann_solver_call(polynomial_basis,
277 face_projections: FaceProjections):
278 """
279 number_of_face_projections: Integer
280 How many quantities are to be projected onto the face. If you pass in one,
281 this means that only the left and right values are projected onto the face.
282 If you pass in two, we store the values plus the projections of the
283 derivative along the normal per face.
284 """
285
286 if isinstance( polynomial_basis, exahype2.solvers.GaussLegendreBasis) and face_projections==FaceProjections.Solution:
287 return """solveRiemannProblem_pointwise_in_situ(
288 [&](
289 const double * __restrict__ Q,
290 const tarch::la::Vector<Dimensions,double>& x,
291 double t,
292 double dt,
293 int normal,
294 double * __restrict__ F
295 )->void {
296 {% if FLUX_IMPLEMENTATION!="<none>" %}
297 repositories::{{SOLVER_INSTANCE}}.flux(Q,x,t,dt,normal,F);
298 {% endif %}
299 }, //flux
300 [&](
301 const double * __restrict__ Q,
302 const double * __restrict__ deltaQ,
303 const tarch::la::Vector<Dimensions,double>& x,
304 double t,
305 double dt,
306 int normal,
307 double * __restrict__ F
308 ) -> void {
309 {% if NCP_IMPLEMENTATION!="<none>" %}
310 repositories::{{SOLVER_INSTANCE}}.nonconservativeProduct(Q, deltaQ, x, t, dt, normal, F);
311 {% endif %}
312 }, //ncp
313 [&](
314 const double * __restrict__ Q,
315 const tarch::la::Vector<Dimensions,double>& x,
316 double t,
317 double dt,
318 int normal,
319 double* maxEigenvalue
320 ) -> void {
321 repositories::{{SOLVER_INSTANCE}}.maxEigenvalue(Q, x, t, dt, normal, maxEigenvalue);
322 }, // maxEigenvalue
323 marker.x(),
324 marker.h(),
325 timeStamp,
326 timeStepSize,
327 {{ORDER}},
328 {{NUMBER_OF_UNKNOWNS}},
329 {{NUMBER_OF_AUXILIARY_VARIABLES}},
330 marker.getSelectedFaceNumber(),
331 repositories::{{SOLVER_INSTANCE}}.QuadraturePoints1d,
332 {{ "true" if FLUX_IMPLEMENTATION!="<none>" else "false" }}, //useFlux
333 {{ "true" if NCP_IMPLEMENTATION!="<none>" else "false" }}, //useNCP
334 fineGridFace{{UNKNOWN_IDENTIFIER}}EstimateProjection.value,
335 fineGridFace{{UNKNOWN_IDENTIFIER}}RiemannSolution.value
336 );
337"""
338 elif isinstance( polynomial_basis, exahype2.solvers.GaussLegendreBasis) and number_of_face_projections==2:
339 return """solveRiemannProblem_pointwise_in_situ_with_gradient_projection(
340 [&](
341 const double * __restrict__ Q,
342 const tarch::la::Vector<Dimensions,double>& x,
343 double t,
344 double dt,
345 int normal,
346 double * __restrict__ F
347 )->void {
348 {% if FLUX_IMPLEMENTATION!="<none>" %}
349 repositories::{{SOLVER_INSTANCE}}.flux(Q,x,t,dt,normal,F);
350 {% endif %}
351 }, //flux
352 [&](
353 const double * __restrict__ Q,
354 const double * __restrict__ deltaQ,
355 const tarch::la::Vector<Dimensions,double>& x,
356 double t,
357 double dt,
358 int normal,
359 double * __restrict__ F
360 ) -> void {
361 {% if NCP_IMPLEMENTATION!="<none>" %}
362 repositories::{{SOLVER_INSTANCE}}.nonconservativeProduct(Q, deltaQ, x, t, dt, normal, F);
363 {% endif %}
364 }, //ncp
365 [&](
366 const double * __restrict__ Q,
367 const tarch::la::Vector<Dimensions,double>& x,
368 double t,
369 double dt,
370 int normal,
371 double* maxEigenvalue
372 ) -> void {
373 repositories::{{SOLVER_INSTANCE}}.maxEigenvalue(Q, x, t, dt, normal, maxEigenvalue);
374 }, // maxEigenvalue
375 marker.x(),
376 marker.h(),
377 timeStamp,
378 timeStepSize,
379 {{ORDER}},
380 {{NUMBER_OF_UNKNOWNS}},
381 {{NUMBER_OF_AUXILIARY_VARIABLES}},
382 marker.getSelectedFaceNumber(),
383 repositories::{{SOLVER_INSTANCE}}.QuadraturePoints1d,
384 {{ "true" if FLUX_IMPLEMENTATION!="<none>" else "false" }}, //useFlux
385 {{ "true" if NCP_IMPLEMENTATION!="<none>" else "false" }}, //useNCP
386 fineGridFace{{UNKNOWN_IDENTIFIER}}EstimateProjection.value,
387 fineGridFace{{UNKNOWN_IDENTIFIER}}RiemannSolution.value
388 );
389"""
390 else:
391 assert False, "not implemented"
392 return "#not implemented"
exahype2.solvers.LagrangeBasisWithDiagonalMassMatrix.GaussLegendreBasis
The Gauss-Legendre Basis is by construction the only basis which yields diagonal mass matrices.
Definition LagrangeBasisWithDiagonalMassMatrix.py:610
exahype2.solvers.rkdg.rusanov.kernels.create_Riemann_solver_call
create_Riemann_solver_call(polynomial_basis, FaceProjections face_projections)
number_of_face_projections: Integer How many quantities are to be projected onto the face.
Definition kernels.py:277
exahype2.solvers.rkdg.rusanov.kernels.create_Riemann_solver_call_dsl
create_Riemann_solver_call_dsl(polynomial_basis, FaceProjections face_projections)
Definition kernels.py:256