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