Peano
Loading...
Searching...
No Matches
RungeKuttaDG.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 os
4
5import peano4
6import peano4.datamodel
12
13from peano4.solversteps.ActionSet import ActionSet
14
15from exahype2.solvers.ButcherTableau import RungeKutta_steps
16
17from exahype2.solvers.PDETerms import PDETerms
18
19import jinja2
20import math
21
22from abc import abstractmethod
23from enum import Enum
24
25import exahype2
27
29 FaceProjections,
30)
32 compute_number_of_face_projection_quantities,
33)
34
35from exahype2.solvers.Storage import Storage
36
37
38class RungeKuttaDG(object):
39 """!
40 An abstract class for any RKDG solver of any order
41
42 This is the base class of all Runge-Kutta DG solvers. It defines what kind of action
43 sets do exist, i.e., what in principle can be done while we run through the grid. It
44 also provides some kind of very basic infrastructure, i.e. what the name of a solver
45 is, what data is to be held per face or cell, or where in the multiscale mesh
46 we actually have to hold data.
47
48 The RungeKuttaDG class cannot/should not be instantiated. Its children actually
49 decide which steps are required or invoked in which order.
50
51 They do so by setting the guards: Each data structure has a guard which controls
52 if it is to be stored. Each action set has a guard that says if it should be
53 called in this particular mesh traversal.
54 If you want to redefine/reset these guards,
55 you have to redefine create_data_structures() and reset the guards after
56 you have called the superclass' operation. I recommend that you refrain from
57 defining completely new guards. Use the predefined guards instead and
58 refine them by adding more and and or clauses.
59
60 If you want to study the mathematics routines used by all the Python-generated
61 code, read the Discontinous Galerkin page in dg.h.
62
63
64 ## Data structures
65
66 Each cell hosts a DG polynomial of order _dg_order. Technically, we store this polynomial
67 in a regular array (blockstructured data) with (_dg_order+1)^d entries. On top of this
68 current time step, we need the shots into the future. for RK(_rk_order), we need _rk_order
69 of these snapshots. They are called estimates.
70
71 On the boundary, I store both the left and right solution projections. Per RK step, we
72 project the linear combination onto the face. This is in line with the action set
73 ComputeFinalLinearCombination which simply copies over the old time step
74 into the linear combination field in the very first Runge-Kutta step.
75
76 Our DG implementation is not tied to one particular Riemann solver, even though most
77 users will likely use Rusanov. Therefore, it is not clear which data is to be projected
78 onto the faces. For plain Rusanov with only a flux function, the solutions are sufficient.
79 Other Riemann solvers or Rusanov's ncp terms however require the derivatives. So I give
80 the code the opportunity to work with different data cardinalities.
81
82 More details are provided in create_data_structures().
83
84
85 ## Control flow between this class and subclasses
86
87 There are three key routines: the constructor, create_data_structures() and
88 create_action_sets(). The constructor sets/memorises some solver variables (such as the
89 name) and then invokes the other two routines.
90
91 create_data_structures() establishes all the data structures tied to the
92 grid entities. If you want to alter the configuration of data tied to grid
93 entities, you should redefine create_data_structures(). However, any
94 subclass still should call Runge Kutta's create_data_structures() before that.
95 This will ensure that the baseline configuration of all data is in place.
96 After that, you can modify the properties.
97
98 create_action_sets() establishes the action sets, i.e., activities that are to
99 be triggered whenever you run a time step, you plot, you initialise the grid.
100 Same here: If you want to alter the configuration of the code, call
101 create_action_sets() of the superclass first, before you do fancy stuff. This
102 way, all the important action sets are in place.
103
104
105 ## General control flow and realisation of solver steps
106
107 A proper numerical time step is mapped onto a cascade of time step sweeps. Per sweep,
108 we realise the following steps (so step in a computational sense), though their exact
109 orchestration depends on the realisation chosen, i.e. some might merge some
110 subcalculations into one grid sweep while others distribute a calculation over multiple
111 sweeps.
112
113 - We combine the data of the current time step with all the shots into the future. These
114 shots are held in one large array. This is not super efficient, as we could build up
115 these arrays one by one, but I just use one large scratchpad even for the first line
116 in the Butcher tableau. See create_data_structures() for details. So all the estimates
117 and the current solution are combined according to the Butcher tableau.
118
119 It is up to the particular subclass to configure the guard of the linear combination
120 such that it is only computed when it is actually needed.
121
122 - Next, we project this linear combination onto the faces. At the same time, we compute
123 a new estimate for the volume integral. The time step size for this new estimate
124 again depends on the time step size from the Butcher tableau. The two action sets
125 both are volumetric, i.e. run per cell, but we can run them in parallel. The face
126 projection has to complete immediately (we need the data on the faces), while the
127 volumetric computation, in principle, doesn't have to finish prior to the next grid
128 sweep, as long as the linear combination remains persistent. It is only a temporary
129 thing, so will be destroyed immediately once we leave a cell.
130
131
132 ## Mandatory step of subclasses
133
134 There are different nuances of subclasses/algorithmic realisations, i.e., few solvers inherit
135 directly from RungeKuttaDG. However, even those subclasses still require you do provide some
136 additional information.
137
138 These are the fields you have to set:
139
140 _compute_time_step_size A simple string which defines a double timeStepSize. This snippet will
141 be used by both the Riemann solver and the volumetric integration. You don't have to care
142 about the Runge-Kutta scaling of time step sizes (or shifts), as this is all done by the
143 source code, but you have to feed an instruction into the solver how to determine a time step
144 size for a cell.
145
146 self._compute_new_time_step_size = "const double newTimeStepSize = timeStepSize;"
147
148
149 ## Semantics of FaceLabel
150
151 - The update time stamp is the one we update in each and every
152 Runge-Kutta step.
153
154 - The new time stamp is the valid one, i.e., the one that
155 corresponds to the face's EstimateProjection after a time step.
156
157 - The old time stamp is the one that corresponds to the (old)
158 face projection.
159
160 - We set the updated flag if and only if we have finished the
161 last Runge-Kutta step.
162 """
164 self,
165 name,
166 rk_order,
167 polynomial_basis,
168 face_projections: FaceProjections,
169 unknowns,
170 auxiliary_variables,
171 min_cell_h,
172 max_cell_h,
173 plot_grid_properties,
174 pde_terms_without_state: bool,
175 baseline_action_set_descend_invocation_order=0,
176 ):
177 """
178 Constructor of the Runge-Kutta Discontinuous Galerkin solver
179
180 :: Arguments
181
182 name: string
183 A unique name for the solver. This one will be used for all generated
184 classes. Also the C++ object instance later on will incorporate this
185 name.
186
187 dg_order: int
188 Order of the Discontinuous Galerkin approximation.
189
190 rk_order: int
191 Runge-Kutta order, i.e., time stepping order.
192
193 polynomials: Polynomials
194 Picks which kind of polynomial to use.
195
196 face_projections: Integer
197 How many projections do we need per face. If you are only
198 interested in the left and right solution along a face, pass in 1. If you
199 want to have the fluxes or derivatives as well, pass in 2. The latter version
200 is the type of data structure you need for Rusanov, e.g.
201
202 There are a couple of pre-defined values for this guy.
203
204 unknowns: int
205 Number of unknowns per Finite Volume voxel.
206
207 auxiliary_variables: int
208 Number of auxiliary variables per Finite Volume voxel. Eventually, both
209 unknowns and auxiliary_variables are merged into one big vector if we
210 work with AoS. But the solver has to be able to distinguish them, as
211 only the unknowns are subject to a hyperbolic formulation.
212
213 min_cell_h: double
214 This size refers to the individual discretisation cell.
215
216 max_cell_h: double
217 This size refers to the individual discretisation cell.
218
219 plot_grid_properties: Boolean
220 Clarifies whether a dump of the data should be enriched with grid info
221 (such as enclave status flags), too.
222
223
224 :: Attributes
225
226 All the constructor parameters are stored in object attributes. There's a few
227 more which are of interest for subclasses:
228
229 _solver_template_file_class_name: String
230 This can be used by subclasses to pick which template is to be used to create
231 the abstract solver and the actual solver blueprint.
232
233
234 _number_of_derivatives_projected_onto_face: Int
235
236 _solver_states: [String]
237 Returns a list of strings of the individual solver states. The code will
238 automatically supplement this list with the following states:
239
240 GridConstruction,
241 GridInitialisation,
242 Plotting,
243 PlottingAfterGridInitialisation
244
245 So you may implicitly assume that these guys do exist always.
246 """
247 assert rk_order >= 1, "Runge-Kutta order has to be greater or equal to one"
248
249 self._name = name
250
251 self._rk_order = rk_order
252 self._basis = polynomial_basis
253 self._face_projections = face_projections
254
255 self._volumetric_compute_kernel_call = "#error Please set the solver property _volumetric_compute_kernel_call in the Python class"
256 self._volumetric_compute_kernel_call_stateless = "#error Please set the solver property _volumetric_compute_kernel_call_stateless in the Python class"
258 "#error Please set the solver property _Riemann_compute_kernel_call in the Python class"
259 )
261 "#error Please set the solver property _Riemann_compute_kernel_call_stateless in the Python class"
262 )
263 self._pde_terms_without_state = pde_terms_without_state
264 self._add_solver_contributions_call = "#error Please set the solver property _add_solver_contributions_call in the Python class"
265 self._multiply_with_inverted_mass_matrix_call = "#error Please set the solver property _multiply_with_inverted_mass_matrix_call in the Python class"
266 self._kernel_namespace = "#error Please set the solver property _kernel_namespace in the Python class"
267
268 """
269 self._unknowns and self._auxiliary_variables respectively hold the number of unknowns and
270 auxiliary variables in the equation to be computed. Unknowns are variables that change over
271 time whereas auxiliary variables can be space-dependent but don't vary over time.
272 These can be specified either as simple ints or by a dictionary
273 (e.g.) unknowns = {'a': 1, 'b': 1, 'c': 3}
274 in which the user specifies the multiplicity of the variable (the velocity has one component
275 per dimension for example.)
276 If they are specified by a dictionary then the code generates a "VariableShortcuts" file which
277 allows the user to specify a variable by name and automatically maps this to the right position
278 in an array for better legibility. Otherwise they must manually remember the position of each
279 variable manually.
280
281 use_var_shortcut is used to know whether or not the user passed their variables via a dict
282 variable_names and variable_pos are used internally to remember the names and respective
283 positions of variables if set by a dictionary.
284 """
285 self._use_var_shortcut = False
287 self._variable_pos = [0]
288
289 if type(unknowns) is dict:
290 self._unknowns = sum(unknowns.values())
291 self._variable_names += list(unknowns.keys())
292 for var in list(unknowns.values()):
293 self._variable_pos.append(var + self._variable_pos[-1])
294 self._use_var_shortcut = True
295 elif type(unknowns) is int:
296 self._unknowns = unknowns
297 else:
298 raise Exception(
299 "Not a valid type for parameter unknowns, needs to be int or dictionary."
300 )
301
302 if type(auxiliary_variables) is dict:
303 self._auxiliary_variables = sum(auxiliary_variables.values())
304 self._variable_names += list(auxiliary_variables.keys())
305 for var in list(auxiliary_variables.values()):
306 self._variable_pos.append(var + self._variable_pos[-1])
307 self._use_var_shortcut = True
308 elif type(auxiliary_variables) is int:
309 self._auxiliary_variables = auxiliary_variables
310 else:
311 raise Exception(
312 "Not a valid type for parameter auxiliary_variables, needs to be int or dictionary."
313 )
314
315 self._min_cell_h = min_cell_h
316 self._max_cell_h = max_cell_h
317 self._plot_grid_properties = plot_grid_properties
318
322
323 if min_cell_h > max_cell_h:
324 raise Exception(
325 "min_cell_h ("
326 + str(min_cell_h)
327 + ") is bigger than max_cell_h ("
328 + str(max_cell_h)
329 + ")"
330 )
331
333
336
339
341
344
348
354
355 self._boundary_conditions_implementation = PDETerms.User_Defined_Implementation
356 self._refinement_criterion_implementation = PDETerms.Empty_Implementation
357 self._initial_conditions_implementation = PDETerms.User_Defined_Implementation
358
359 self._compute_time_step_size = "#error Not yet defined"
360
365
368
370
371 self._cell_data_storage = Storage.CallStack
372 self._face_data_storage = Storage.CallStack
373
374 self._baseline_action_set_descend_invocation_order = baseline_action_set_descend_invocation_order
375 self.switch_storage_scheme(Storage.SmartPointers, Storage.SmartPointers)
376
377
378 def __str__(self):
379 result = (
380 """
381Name: """
382 + self._name
383 + """
384Type: """
385 + self.__class__.__module__
386 + """
387Stateless PDE terms: """
388 + str(self._pde_terms_without_state)
389 + """
390RK/DG order: """
391 + str(self._rk_order)
392 + "/"
393 + str(self._basis.order)
394 + """
395Unknowns: """
396 + str(self._unknowns)
397 + """
398Auxiliary variables: """
399 + str(self._auxiliary_variables)
400 + """
401h_cell_min: """
402 + str(self._min_cell_h)
403 + """
404h_cell_max: """
405 + str(self._max_cell_h)
406 + """
407Initial conditions: """
409 + """
410Boundary conditions: """
412 + """
413Refinement criterion: """
415 + """
416Eigenvalues: """
418 + """
419Flux: """
420 + str(self._flux_implementation)
421 + """
422NCP: """
423 + str(self._ncp_implementation)
424 + """
425Source term: """
427 + """
428Point source: """
430 + """
431"""
432 )
433 return result
434
435 __repr__ = __str__
436
437
439 coarsest_tree_level = 0
440 while domain_size * 3 ** (-coarsest_tree_level) > self._max_cell_h:
441 coarsest_tree_level += 1
442 return coarsest_tree_level
443
444
446 finest_tree_level = 0
447 while domain_size * 3 ** (-finest_tree_level) > self._min_cell_h:
448 finest_tree_level += 1
449 return finest_tree_level
450
451
452 def get_coarsest_number_of_cells(self, domain_size):
453 return 3 ** self.get_min_number_of_spacetree_levels(domain_size)
454
455
456 def get_finest_number_of_cells(self, domain_size):
457 return 3 ** self.get_max_number_of_spacetree_levels(domain_size)
458
459
460 def create_readme_descriptor(self, domain_offset, domain_size):
461 return (
462 """
463### ExaHyPE 2 solver
464"""
465 + str(self)
466 + """
467
468Real type of this solver: """
469 + str(type(self))
470 + """
471
472We assume that you use a domain size of (0,"""
473 + str(domain_size)
474 + """)^d. Peano 4 will cut this domain equidistantly
475and recursively into three parts along each coordinate axis. This yields a spacetree.
476
477The spacetree will at least have """
478 + str(self.get_min_number_of_spacetree_levels(domain_size))
479 + """ levels.
480The spacetree will at most have """
481 + str(self.get_max_number_of_spacetree_levels(domain_size))
482 + """ levels.
483
484The spacetree will thus span at least """
485 + str(self.get_coarsest_number_of_cells(domain_size))
486 + """ cells per coordinate axis.
487The spacetree will thus span at most """
488 + str(self.get_finest_number_of_cells(domain_size))
489 + """ cells per coordinate axis.
490
491We use RK("""
492 + str(self._rk_order)
493 + """) for the time stepping.
494
495"""
496 )
497
498
499 @property
501 """
502 Add further includes to this property, if your action sets require some additional
503 routines from other header files.
504 """
505 return self._user_action_set_includes
506
507
508 @property
510 """
511 Add further includes to this property, if your solver requires some additional
512 routines from other header files.
513 """
514 return self._user_solver_includes
515
516
518 """
519 Add further includes to this property, if your action sets require some additional
520 routines from other header files.
521 """
522 self._user_action_set_includes += value
523
524
525 def add_user_solver_includes(self, value):
526 """
527 Add further includes to this property, if your solver requires some additional
528 routines from other header files.
529 """
530 self._user_solver_includes += value
531
532
533 @abstractmethod
535 """
536 Recall in subclasses if you wanna change the number of unknowns
537 or auxiliary variables. See class description's subsection on
538 data flow.
539
540 :: Call order and ownership
541
542 This operation can be called multiple times. However, only the very
543 last call matters. All previous calls are wiped out.
544
545 If you have a hierarchy of solvers, every create_data_structure()
546 should first(!) call its parent version. This way, you always ensure
547 that all data are in place before you continue to alter the more
548 specialised versions. So it is (logically) a top-down (general to
549 specialised) run through all create_data_structure() variants
550 within the inheritance tree.
551
552 :: Solver fields built up
553
554 _time_step: Patch ( (_dg_order+1)x(_dg_order+1)x(_dg_order+1) )
555 Actual polynomial data of the DG patch in the current time step.
556
557 _rhs_estimates: Patch ( (_dg_order+1)x(_dg_order+1)x(_dg_order+1)xno of RK steps )
558 These are the Runge-Kutta extrapolations. It is important that the
559 rk_order is mixed into the first argument, as the triple is
560 automatically truncated to a tuple if the code is translated with 2d.
561 Please note that the term rhs refers to the right-hand side of an ODE
562
563 @f$ \partial Q = F(Q) @f$
564
565 and thus does not include auxiliary variables.
566
567 _linear_combination_of_estimates: Patch ( (_dg_order+1)x(_dg_order+1)x(_dg_order+1) )
568 Non-persistent helper data structure.
569
570 _current_time_step_projection: Patch ( 2xface_projectionsx(_dg_order+1)x(_dg_order+1) )
571 This is the projection of the polynomial data of the current guess onto
572 the faces. It is important that the first dimension is only a 2. This
573 way, we can (logically) handle the projected data again as an overlap
574 of two patches with the width 1.
575
576 _estimate_projection: Patch ( Nxface_projectionsx(_dg_order+1)x(_dg_order+1) )
577 The N is usually 2. Usually this estimate projection holds the data
578 from the latest Runge-Kutta step. This implies that it never ever holds
579 the projection of the solution unless after the very first Runge-Kutta
580 step. After that, it is always a linear combination of estimates.
581
582 To avoid this plain behaviour, I make the last and final linear
583 combination project the solution again onto the faces, i.e. after
584 the very last step of the Runge-Kutta scheme, there will be a valid
585 representation of the solution _estimate_projection. Consult the
586 documentation of exahype2.solvers.rkdg.actionsets.ComputeFinalLinearCombination
587 for some further information. Also consult the semantics of the FaceLabel
588 in this context.
589
590 _old_solution_projection: Same as _estimate_projection
591 This is a backup of the final solution stored in _estimate_projection.
592 It correlates to
593
594 _face_label: FaceLabel
595 See class description. General information such as "is boundary".
596
597 _cell_label: CellLabel
598 See class description. General information such as "is enclave".
599
600 Per default, I do always store the projections onto faces, and I
601 always keep the actual time step data and the cell projections. It is
602 up to the subclasses to alter this storage behaviour.
603
604 By default, I do not exchange any face data in-between two grid sweeps.
605 You will have to overwrite the behaviour in subclasses. Please note that
606 a face is communicated between two tree/domain partitions if and only if
607 it is stored persistently and the exchange flag is set, too.
608 """
610 (
611 self._basis.dofs_per_axis,
612 self._basis.dofs_per_axis,
613 self._basis.dofs_per_axis,
614 ),
616 self._unknown_identifier(),
617 )
619 (
620 self.number_of_Runge_Kutta_steps() * (self._basis.dofs_per_axis),
621 self._basis.dofs_per_axis,
622 self._basis.dofs_per_axis,
623 ),
624 self._unknowns,
625 self._unknown_identifier() + "RhsEstimates",
626 )
628 (
629 self._basis.dofs_per_axis,
630 self._basis.dofs_per_axis,
631 self._basis.dofs_per_axis,
632 ),
634 self._unknown_identifier() + "LinearCombination",
635 )
637 (
638 compute_number_of_face_projection_quantities(self._face_projections)
639 * 2,
640 self._basis.dofs_per_axis,
641 self._basis.dofs_per_axis,
642 ),
644 self._unknown_identifier() + "EstimateProjection",
645 )
647 (
648 compute_number_of_face_projection_quantities(self._face_projections)
649 * 2,
650 self._basis.dofs_per_axis,
651 self._basis.dofs_per_axis,
652 ),
654 self._unknown_identifier() + "OldSolutionProjection",
655 )
657 (2, self._basis.dofs_per_axis, self._basis.dofs_per_axis),
658 self._unknowns,
659 self._unknown_identifier() + "RiemannSolution",
660 )
661
662 if self._cell_data_storage == Storage.CallStack:
664 self._current_time_step, "double"
665 )
667 self._rhs_estimates, "double"
668 )
669 self._linear_combination_of_estimates.generator = (
672 )
673 )
674 elif self._cell_data_storage == Storage.Heap:
675 self._current_time_step.generator = (
677 self._current_time_step, "double"
678 )
679 )
681 self._rhs_estimates, "double"
682 )
683 self._linear_combination_of_estimates.generator = (
686 )
687 )
688 elif self._cell_data_storage == Storage.SmartPointers:
689 self._current_time_step.generator = (
691 self._current_time_step, "double"
692 )
693 )
694 self._rhs_estimates.generator = (
696 self._rhs_estimates, "double"
697 )
698 )
699 self._linear_combination_of_estimates.generator = (
702 )
703 )
704
705 if self._face_data_storage == Storage.CallStack:
707 self._estimate_projection, "double"
708 )
709 self._old_solution_projection.generator = (
711 self._old_solution_projection, "double"
712 )
713 )
715 self._Riemann_solution, "double"
716 )
717 elif self._face_data_storage == Storage.Heap:
718 self._estimate_projection.generator = (
720 self._estimate_projection, "double"
721 )
722 )
723 self._old_solution_projection.generator = (
725 self._old_solution_projection, "double"
726 )
727 )
728 self._Riemann_solution.generator = (
730 self._Riemann_solution, "double"
731 )
732 )
733 elif self._face_data_storage == Storage.SmartPointers:
734 self._estimate_projection.generator = (
736 self._estimate_projection, "double"
737 )
738 )
739 self._old_solution_projection.generator = (
741 self._old_solution_projection, "double"
742 )
743 )
744 self._Riemann_solution.generator = (
746 self._Riemann_solution, "double"
747 )
748 )
749 else:
750 assert False, "Storage variant {} not supported".format(
752 )
753
754 self._estimate_projection.generator.load_store_compute_flag = "::peano4::grid::constructLoadStoreComputeFlag({},{},{})".format(
758 )
759
760 self._estimate_projection.generator.includes += """
761#include "peano4/utils/Loop.h"
762#include "repositories/SolverRepository.h"
763"""
764 self._estimate_projection.generator.merge_method_definition = (
765 peano4.toolbox.blockstructured.get_face_merge_implementation(
767 )
768 )
769
770 self._old_solution_projection.generator.load_store_compute_flag = "::peano4::grid::constructLoadStoreComputeFlag({},{},{})".format(
774 )
775 self._old_solution_projection.generator.includes += """
776#include "peano4/utils/Loop.h"
777#include "repositories/SolverRepository.h"
778"""
779
780 self._current_time_step.generator.load_store_compute_flag = "::peano4::grid::constructLoadStoreComputeFlag({},{},{})".format(
784 )
785 self._current_time_step.generator.includes += """
786#include "peano4/utils/Loop.h"
787#include "repositories/SolverRepository.h"
788"""
789 self._current_time_step.generator.merge_method_definition = (
790 peano4.toolbox.blockstructured.get_cell_merge_implementation(
792 )
793 )
794
795 self._rhs_estimates.generator.load_store_compute_flag = "::peano4::grid::constructLoadStoreComputeFlag({},{},{})".format(
799 )
800 self._rhs_estimates.generator.includes += """
801#include "peano4/utils/Loop.h"
802#include "repositories/SolverRepository.h"
803"""
804
805 self._linear_combination_of_estimates.generator.load_store_compute_flag = "::peano4::grid::constructLoadStoreComputeFlag({},{},{})".format(
807 "false",
808 "false",
809 )
810 self._linear_combination_of_estimates.generator.includes += """
811#include "peano4/utils/Loop.h"
812#include "repositories/SolverRepository.h"
813"""
814 self._estimate_projection.generator.send_condition = "false"
815 self._estimate_projection.generator.receive_and_merge_condition = "false"
816
817 self._old_solution_projection.generator.send_condition = "false"
818 self._old_solution_projection.generator.receive_and_merge_condition = "false"
819
820 self._Riemann_solution.generator.load_store_compute_flag = "::peano4::grid::constructLoadStoreComputeFlag({},{},{})".format(
824 )
825 self._Riemann_solution.generator.includes += """
826#include "repositories/SolverRepository.h"
827"""
828
829 self._Riemann_solution.generator.send_condition = "false"
830 self._Riemann_solution.generator.receive_and_merge_condition = "false"
831
832 self._current_time_step.generator.includes += """
833#include "../repositories/SolverRepository.h"
834"""
835 self._rhs_estimates.generator.includes += """
836#include "../repositories/SolverRepository.h"
837"""
838 self._estimate_projection.generator.includes += """
839#include "../repositories/SolverRepository.h"
840"""
841 self._Riemann_solution.generator.includes += """
842#include "../repositories/SolverRepository.h"
843"""
844
845 self._cell_label = exahype2.grid.create_cell_label(self._name)
846 self._face_label = exahype2.grid.create_face_label(self._name)
847
848
853 @abstractmethod
855 """
856 Overwrite in subclasses if you wanna create different
857 action sets.
858
859 :: Call order and ownership
860
861 This operation can be called multiple times over the construction of
862 a solver. However, only the very last call matters. All previous calls
863 are wiped out.
864
865 If you have a hierarchy of solvers, every create_data_structure()
866 should first(!) call its parent version. This way, you always ensure
867 that all data are in place before you continue to alter the more
868 specialised versions. So it is (logically) a top-down (general to
869 specialised) run through all create_data_structure() variants
870 within the inheritance tree.
871 """
874 self,
876 "true",
877 )
878 )
881 self,
883 "false",
884 )
885 )
887 solver=self,
889 + " and (repositories::"
891 + ".isLastGridSweepOfTimeStep() or repositories::"
893 + ".getSolverState()=="
894 + self._name
895 + "::SolverState::GridInitialisation)",
896 build_up_new_refinement_instructions=True,
897 implement_previous_refinement_instructions=True,
898 called_by_grid_construction=False,
899 )
902 solver=self,
904 build_up_new_refinement_instructions=False,
905 implement_previous_refinement_instructions=True,
906 called_by_grid_construction=False,
907 )
908 )
911 solver=self,
913 build_up_new_refinement_instructions=True,
914 implement_previous_refinement_instructions=True,
915 called_by_grid_construction=True,
916 )
917 )
920 )
921 # @todo TW The following two action sets should run in parallel I guess
923 solver=self, face_projections=self._face_projections
924 )
927 )
931 )
932 )
935 )
938 )
940 solver=self,
942 face_projections=self._face_projections,
943 )
945 solver=self, face_projections=self._face_projections
946 )
947
948 if self._action_set_postprocess_solution == None:
951 )
952 if self._action_set_preprocess_solution == None:
955 )
956
959
960 self._action_set_update_face_label.descend_invocation_order = (
962 )
963 self._action_set_update_cell_label.descend_invocation_order = (
965 )
966 self._action_set_linear_combination_of_estimates.descend_invocation_order = (
968 )
971 )
972 self._action_set_project_linear_combination_onto_faces.descend_invocation_order = (
974 ) # depends on linear combination of estimates
975 self._action_set_solve_volume_integral.descend_invocation_order = (
977 ) # the two 3s could run in parallel
978 self._action_set_handle_boundary.descend_invocation_order = (
980 )
981 self._action_set_solve_Riemann_problem.descend_invocation_order = (
983 ) # boundary has to be known
984 self._action_set_add_volume_and_face_solution.descend_invocation_order = (
986 ) # can also run in parallel, as it never occurs at same time
989 ) # depends on other guys
990 self._action_set_AMR.descend_invocation_order = self._baseline_action_set_descend_invocation_order + 5
992 self._action_set_AMR_commit_without_further_analysis.descend_invocation_order = (
994 )
995
996 self._action_set_initial_conditions.descend_invocation_order = (
998 )
999 self._action_set_initial_conditions_for_grid_construction.descend_invocation_order = (
1001 )
1002
1003 self._action_set_preprocess_solution.descend_invocation_order = (
1005 ) # touch cell first time, i.e. enter cell
1006 self._action_set_postprocess_solution.descend_invocation_order = (
1008 ) # touch cell last time, i.e. leave cell
1009
1010
1012 return (
1013 "not marker.willBeRefined() "
1014 + "and repositories::"
1016 + ".getSolverState()!="
1017 + self._name
1018 + "::SolverState::GridConstruction"
1019 )
1020
1021
1023 return (
1024 "not marker.willBeRefined() "
1025 + "and repositories::"
1027 + ".getSolverState()!="
1028 + self._name
1029 + "::SolverState::GridConstruction"
1030 )
1031
1032
1034 """!
1035 Extend the guard via ands only. Never use an or, as subclasses might
1036 extend it as well, and they will append further ends.
1037 """
1038 return (
1039 "not marker.willBeRefined() "
1040 + "and repositories::"
1042 + ".getSolverState()!="
1043 + self._name
1044 + "::SolverState::GridConstruction"
1045 )
1046
1047
1049 """!
1050 Extend the guard via ands only. Never use an or, as subclasses might
1051 extend it as well, and they will append further ends.
1052 """
1053 return (
1054 "not marker.hasBeenRefined() "
1055 + "and repositories::"
1057 + ".getSolverState()!="
1058 + self._name
1059 + "::SolverState::GridConstruction "
1060 + "and repositories::"
1062 + ".getSolverState()!="
1063 + self._name
1064 + "::SolverState::GridInitialisation"
1065 )
1066
1067
1069 """!
1070 Extend the guard via ands only. Never use an or, as subclasses might
1071 extend it as well, and they will append further ends.
1072 """
1073 return (
1074 "not marker.willBeRefined() "
1075 + "and repositories::"
1077 + ".getSolverState()!="
1078 + self._name
1079 + "::SolverState::GridConstruction"
1080 )
1081
1082
1084 """!
1085 Extend the guard via ands only. Never use an or, as subclasses might
1086 extend it as well, and they will append further ends.
1087 """
1088 return (
1089 "not marker.hasBeenRefined() "
1090 + "and repositories::"
1092 + ".getSolverState()!="
1093 + self._name
1094 + "::SolverState::GridConstruction "
1095 + "and repositories::"
1097 + ".getSolverState()!="
1098 + self._name
1099 + "::SolverState::GridInitialisation"
1100 )
1101
1102
1104 return self._name + "Q"
1105
1106
1108 return "instanceOf" + self._name
1109
1110
1111 def add_to_Peano4_datamodel(self, datamodel, verbose):
1112 """
1113 Add all required data to the Peano4 project's datamodel
1114 so it is properly built up
1115 """
1116 if verbose:
1117 print("Polynomial basis")
1118 print("----------")
1119 print(str(self._basis))
1120 print("Face projections")
1121 print("----------")
1122 print(str(self._face_projections))
1123 datamodel.add_cell(self._cell_label)
1124 datamodel.add_cell(self._current_time_step)
1125 datamodel.add_cell(self._rhs_estimates)
1126 datamodel.add_cell(self._linear_combination_of_estimates)
1127 datamodel.add_face(self._face_label)
1128 datamodel.add_face(self._estimate_projection)
1129 datamodel.add_face(self._old_solution_projection)
1130 datamodel.add_face(self._Riemann_solution)
1131
1132
1134 """
1135 Tell Peano what data to move around
1136
1137 Inform Peano4 step which data are to be moved around via the
1138 use_cell and use_face commands. This operation is generic from
1139 ExaHyPE's point of view, i.e. I use it for all grid sweep types.
1140 """
1141 step.use_cell(self._cell_label)
1142 step.use_cell(self._current_time_step)
1143 step.use_cell(self._rhs_estimates)
1144 step.use_cell(self._linear_combination_of_estimates)
1145 step.use_face(self._face_label)
1146 step.use_face(self._estimate_projection)
1147 step.use_face(self._old_solution_projection)
1148 step.use_face(self._Riemann_solution)
1149
1150
1152 return """
1153#include "tarch/la/Vector.h"
1154
1155#include "peano4/utils/Globals.h"
1156#include "peano4/utils/Loop.h"
1157
1158#include "repositories/SolverRepository.h"
1159"""
1160
1162 """!
1163 Add the action sets to the grid initialisation
1164
1165 The AMR stuff has to be the very first thing. Actually, the AMR routines'
1166 interpolation doesn't play any role here. But the restriction indeed is
1167 very important, as we have to get the face data for BCs et al. The action
1168 set order is inverted while we ascend within the tree again. Therefore, we
1169 add the AMR action set first which means it will be called last when we go
1170 from fine to coarse levels within the tree.
1171
1172 ## Ordering
1173
1174 The order of the action sets is preserved throughout the steps down within
1175 the tree hierarchy. It is inverted throughout the backrolling.
1176
1177 This is what we want to achieve:
1178
1179 - Restrict the data to the coarser level if we are on a hanging face.
1180 """
1181 step.add_action_set(self._action_set_preprocess_solution)
1182 step.add_action_set(self._action_set_initial_conditions)
1183 step.add_action_set(self._action_set_update_face_label)
1184 step.add_action_set(self._action_set_update_cell_label)
1185 step.add_action_set(self._action_set_AMR)
1186 step.add_action_set(self._action_set_postprocess_solution)
1187
1188
1189 def add_actions_to_create_grid(self, step, evaluate_refinement_criterion):
1190 """
1191 @todo:
1192 The boundary information is set only once. It is therefore important that
1193 we ues the face label and initialise it properly.
1194 """
1196 step.add_action_set(self._action_set_update_face_label)
1197 step.add_action_set(self._action_set_update_cell_label)
1198 if evaluate_refinement_criterion:
1199 step.add_action_set(self._action_set_AMR_throughout_grid_construction)
1200 else:
1201 step.add_action_set(self._action_set_AMR_commit_without_further_analysis)
1202
1203
1204 def set_plot_description(self, description):
1205 """
1206 Use this one to set a description within the output patch file that tells
1207 the vis solver what the semantics of the entries are. Typicallly, I use
1208 a comma-separated list here.
1209 """
1210 self.plot_description = description
1211
1212
1213 def add_actions_to_plot_solution(self, step, output_path):
1214 """!
1215 Dump snapshot of solution
1216
1217 Consult the discussion in add_actions_to_init_grid() around the order
1218 of the individual action sets.
1219
1220 It is important that we have the coupling/dynamic AMR part in here, as
1221 there might be pending AMR refinement requests that now are realised.
1222 For the same reason, we need the update of the face label and the update
1223 of the cell label in here: The AMR might just propagate over into the
1224 plotting, i.e. we might create new grid entities throughout the plot.
1225 These entities (faces and cells) have to be initialised properly.
1226 Otherwise, their un-initialised data will propagate through to the next
1227 time step.
1228
1229 To make the restriction work, we have to project the solutions onto the
1230 faces.
1231 """
1232 d = {}
1235
1236 # step.add_action_set( self._action_set_couple_resolution_transitions_and_handle_dynamic_mesh_refinement )
1237 # step.add_action_set( self._action_set_roll_over_update_of_faces )
1238 step.add_action_set(self._action_set_update_face_label)
1239 step.add_action_set(self._action_set_update_cell_label)
1240
1241 mapping = []
1242 for z in self._basis.quadrature_points:
1243 for y in self._basis.quadrature_points:
1244 for x in self._basis.quadrature_points:
1245 mapping.append((x, y, z))
1246
1248 filename=output_path + "solution-" + self._name,
1249 patch=self._current_time_step,
1250 dataset_name=self._unknown_identifier(),
1251 description=self.plot_description,
1252 mapping=mapping,
1253 guard="repositories::plotFilter.plotPatch(marker) and "
1255 additional_includes="""
1256#include "exahype2/PlotFilter.h"
1257#include "../repositories/SolverRepository.h"
1258""",
1259 precision="PlotterPrecision",
1260 time_stamp_evaluation="repositories::getMinTimeStamp()",
1261 plot_cell_data=False,
1262 select_dofs=self.select_dofs_to_print,
1263 )
1264 plot_patches_action_set.descend_invocation_order = self._baseline_action_set_descend_invocation_order
1265 step.add_action_set(plot_patches_action_set)
1266
1267 if self._plot_grid_properties:
1268 plot_grid_properties_action_set = peano4.toolbox.PlotGridInPeanoBlockFormat(
1269 filename=output_path + "grid-" + self._name,
1270 cell_unknown=None,
1271 guard="repositories::plotFilter.plotPatch(marker) and "
1273 additional_includes="""
1274#include "exahype2/PlotFilter.h"
1275#include "../repositories/SolverRepository.h"
1276""",
1277 plot_cell_data=False,
1278 )
1279 plot_grid_properties_action_set.descend_invocation_order = (
1281 )
1282 step.add_action_set(plot_grid_properties_action_set)
1283
1284 pass
1285
1286
1288 """
1289 The tricky part here is that we have to get the order right.
1290
1291 - Update of the labels: This can be done as very first step, as it
1292 might feed into follow-up steps.
1293 - Linear combination: This has to be the first step of a Runge-Kutta
1294 scheme. Without this first step, there's absolutely nothing we can
1295 do properly.
1296 - Couple resolutions and AMR: Once we have the linear combination,
1297 we can couple different resolutions.
1298 - Project linear combination onto faces: This should be done as soon
1299 as possible but after we've determined the linear combination. After
1300 all, we need the solution representation on the faces asap.
1301 - Solve volume integral: can run in parallel to the projection onto
1302 the faces.
1303 - Handle boundary: doesn't really matter when we insert it, as it plugs
1304 into the first face load, while all the other stuff so far are volumetric
1305 operations.
1306 It is important that we do the inter-grid transfer operators before we
1307 apply the boundary conditions.
1308 - Solve Riemann problem: the constraint here is that it has to come after
1309 the boundary handling.
1310 - Add volume and face solution: last
1311
1312 If we use enclave tasking, we have to be careful when we insert the merger.
1313 See SeparateSweepsWithEnclaveTasking.add_actions_to_perform_time_step() for
1314 a discussion of the details.
1315 """
1316 d = {}
1319
1320 step.add_action_set(self._action_set_preprocess_solution)
1321 step.add_action_set(self._action_set_update_face_label)
1322 step.add_action_set(self._action_set_update_cell_label)
1323 step.add_action_set(self._action_set_linear_combination_of_estimates)
1324 step.add_action_set(
1326 )
1328 step.add_action_set(self._action_set_solve_volume_integral)
1329 step.add_action_set(self._action_set_handle_boundary)
1330 step.add_action_set(self._action_set_solve_Riemann_problem)
1331 step.add_action_set(self._action_set_add_volume_and_face_solution)
1332 step.add_action_set(
1334 )
1335 step.add_action_set(self._action_set_AMR)
1336 step.add_action_set(self._action_set_postprocess_solution)
1337
1338
1339 @abstractmethod
1341 pass
1342
1343
1344 def add_implementation_files_to_project(self, namespace, output, dimensions, subdirectory=""):
1345 """
1346 The ExaHyPE project will call this operation when it sets
1347 up the overall environment.
1348
1349 This routine is typically not invoked by a user.
1350
1351 output: peano4.output.Output
1352 """
1353 templatefile_prefix = os.path.dirname(os.path.realpath(__file__)) + "/"
1354
1355 if self._solver_template_file_class_name is None:
1356 templatefile_prefix += self.__class__.__name__
1357 else:
1358 templatefile_prefix += self._solver_template_file_class_name
1359
1360 if(subdirectory):
1361 subdirectory += "/"
1362
1363 abstractHeaderDictionary = {}
1364 implementationDictionary = {}
1365 self._init_dictionary_with_default_parameters(abstractHeaderDictionary)
1366 self._init_dictionary_with_default_parameters(implementationDictionary)
1367 self.add_entries_to_text_replacement_dictionary(abstractHeaderDictionary)
1368 self.add_entries_to_text_replacement_dictionary(implementationDictionary)
1369
1370 generated_abstract_header_file = (
1372 templatefile_prefix + "Abstract.template.h",
1373 templatefile_prefix + "Abstract.template.cpp",
1374 "Abstract" + self._name,
1375 namespace,
1376 subdirectory + ".",
1377 abstractHeaderDictionary,
1378 True,
1379 True,
1380 )
1381 )
1382 generated_solver_files = (
1384 templatefile_prefix + ".template.h",
1385 templatefile_prefix + ".template.cpp",
1386 self._name,
1387 namespace,
1388 subdirectory + ".",
1389 implementationDictionary,
1390 False,
1391 True,
1392 )
1393 )
1394
1395 output.add(generated_abstract_header_file)
1396 output.add(generated_solver_files)
1397 output.makefile.add_h_file(subdirectory + "Abstract" + self._name + ".h", generated=True)
1398 output.makefile.add_h_file(subdirectory + self._name + ".h", generated=True)
1399 output.makefile.add_cpp_file(subdirectory + "Abstract" + self._name + ".cpp", generated=True)
1400 output.makefile.add_cpp_file(subdirectory + self._name + ".cpp", generated=True)
1401
1402 if self._use_var_shortcut:
1403 generated_shortcut_file = peano4.output.Jinja2TemplatedHeaderFile(
1404 os.path.dirname(os.path.realpath(__file__))
1405 + "/"
1406 + "../VariableShortcuts.template.h",
1407 "VariableShortcuts",
1408 namespace,
1409 subdirectory + ".",
1410 implementationDictionary,
1411 True,
1412 True,
1413 )
1414 output.add(generated_shortcut_file)
1415 output.makefile.add_h_file(subdirectory + "VariableShortcuts.h", generated=True)
1416
1417
1418 def set_solver_constants(self, datastring):
1419 self._solver_constants = datastring
1420
1421
1422 def add_solver_constants(self, datastring):
1423 self._solver_constants += datastring
1424
1425
1427 """
1428 This one is called by all algorithmic steps before I invoke
1429 add_entries_to_text_replacement_dictionary().
1430
1431 See the remarks on set_postprocess_updated_cell_kernel to understand why
1432 we have to apply the (partially befilled) dictionary to create a new entry
1433 for this very dictionary.
1434 """
1435 d["DG_ORDER"] = self._basis.order
1436 d["RK_ORDER"] = self._rk_order
1437 d["RK_STEPS"] = self.number_of_Runge_Kutta_steps()
1438
1439 d["SOLVER_INSTANCE"] = self.get_name_of_global_instance()
1440 d["SOLVER_NAME"] = self._name
1441 d["UNKNOWN_IDENTIFIER"] = self._unknown_identifier()
1442 d["NUMBER_OF_UNKNOWNS"] = self._unknowns
1443 d["NUMBER_OF_AUXILIARY_VARIABLES"] = self._auxiliary_variables
1444
1445 d["NUMBER_OF_DOFS_PER_CELL_2D"] = (self._basis.dofs_per_axis) * (
1446 self._basis.dofs_per_axis
1447 )
1448 d["NUMBER_OF_DOFS_PER_CELL_3D"] = (
1449 (self._basis.dofs_per_axis)
1450 * (self._basis.dofs_per_axis)
1451 * (self._basis.dofs_per_axis)
1452 )
1453
1454 d["NUMBER_OF_DOFS_PER_FACE_2D"] = self._basis.dofs_per_axis
1455 d["NUMBER_OF_DOFS_PER_FACE_3D"] = (self._basis.dofs_per_axis) * (
1456 self._basis.dofs_per_axis
1457 )
1458
1459 d["ASSERTION_WITH_1_ARGUMENTS"] = "nonCriticalAssertion1"
1460 d["ASSERTION_WITH_2_ARGUMENTS"] = "nonCriticalAssertion2"
1461 d["ASSERTION_WITH_3_ARGUMENTS"] = "nonCriticalAssertion3"
1462 d["ASSERTION_WITH_4_ARGUMENTS"] = "nonCriticalAssertion4"
1463 d["ASSERTION_WITH_5_ARGUMENTS"] = "nonCriticalAssertion5"
1464 d["ASSERTION_WITH_6_ARGUMENTS"] = "nonCriticalAssertion6"
1465
1466 if self._min_cell_h > self._max_cell_h:
1467 raise Exception("min/max h are inconsistent")
1468 d["MAX_CELL_H"] = self._max_cell_h
1469 d["MIN_CELL_H"] = self._min_cell_h
1470
1471 d["SOLVER_CONSTANTS"] = self._solver_constants
1472
1473 d["SOLVER_INCLUDES"] = self.user_solver_includes
1474
1475 # d[ "SOURCE_TERM_CALL"] = jinja2.Template(self._source_term_call, undefined=jinja2.DebugUndefined).render( **d )
1476 # d[ "RIEMANN_SOLVER_CALL"] = jinja2.Template(self._Riemann_solver_call, undefined=jinja2.DebugUndefined).render( **d )
1477 d["PREPROCESS_RECONSTRUCTED_CELL"] = jinja2.Template(
1478 self._preprocess_cell, undefined=jinja2.DebugUndefined
1479 ).render(**d)
1480 d["POSTPROCESS_UPDATED_CELL_AFTER_RUNGE_KUTTA_STEP"] = jinja2.Template(
1482 undefined=jinja2.DebugUndefined,
1483 ).render(**d)
1484 d["POSTPROCESS_UPDATED_CELL_AFTER_FINAL_LINEAR_COMBINATION"] = jinja2.Template(
1486 undefined=jinja2.DebugUndefined,
1487 ).render(**d)
1488 d[
1489 "BOUNDARY_CONDITIONS_IMPLEMENTATION"
1491 d[
1492 "REFINEMENT_CRITERION_IMPLEMENTATION"
1494 d["INITIAL_CONDITIONS_IMPLEMENTATION"] = self._initial_conditions_implementation
1495 d["ABSTRACT_SOLVER_USER_DECLARATIONS"] = jinja2.Template(
1496 self._abstract_solver_user_declarations, undefined=jinja2.DebugUndefined
1497 ).render(**d)
1498 d["ABSTRACT_SOLVER_USER_DEFINITIONS"] = jinja2.Template(
1499 self._abstract_solver_user_definitions, undefined=jinja2.DebugUndefined
1500 ).render(**d)
1501 d["SOLVER_USER_DECLARATIONS"] = jinja2.Template(
1502 self._solver_user_declarations, undefined=jinja2.DebugUndefined
1503 ).render(**d)
1504 d["SOLVER_USER_DEFINITIONS"] = jinja2.Template(
1505 self._solver_user_definitions, undefined=jinja2.DebugUndefined
1506 ).render(**d)
1507 d["START_TIME_STEP_IMPLEMENTATION"] = jinja2.Template(
1508 self._start_time_step_implementation, undefined=jinja2.DebugUndefined
1509 ).render(**d)
1510 d["FINISH_TIME_STEP_IMPLEMENTATION"] = jinja2.Template(
1511 self._finish_time_step_implementation, undefined=jinja2.DebugUndefined
1512 ).render(**d)
1513 d["CONSTRUCTOR_IMPLEMENTATION"] = jinja2.Template(
1514 self._constructor_implementation, undefined=jinja2.DebugUndefined
1515 ).render(**d)
1516 d["COMPUTE_TIME_STEP_SIZE"] = jinja2.Template(
1517 self._compute_time_step_size, undefined=jinja2.DebugUndefined
1518 ).render(**d)
1519 d["COMPUTE_NEW_TIME_STEP_SIZE"] = jinja2.Template(
1520 self._compute_new_time_step_size, undefined=jinja2.DebugUndefined
1521 ).render(**d)
1522
1523 d["FLUX_IMPLEMENTATION"] = self._flux_implementation
1524 d["NCP_IMPLEMENTATION"] = self._ncp_implementation
1525 d["SOURCE_TERM_IMPLEMENTATION"] = self._source_term_implementation
1526 d["POINT_SOURCES_IMPLEMENTATION"] = self._point_sources_implementation
1527
1528 d["COMPUTE_MAX_EIGENVALUE"] = self._compute_eigenvalue
1529
1530 d["STATELESS_PDE_TERMS"] = self._pde_terms_without_state
1531 d["KERNEL_NAMESPACE"] = self._kernel_namespace
1532
1533 d["USE_VARIABLE_SHORTCUT"] = self._use_var_shortcut
1534 d["VARIABLE_NAMES"] = self._variable_names
1535 d["VARIABLE_POSITIONS"] = self._variable_pos
1536
1537 self._basis.init_dictionary_with_default_parameters(d, False)
1538
1539 # Has to come last, as we already use the other entries
1540 d["VOLUMETRIC_COMPUTE_KERNEL_CALL"] = jinja2.Template(
1541 self._volumetric_compute_kernel_call, undefined=jinja2.DebugUndefined
1542 ).render(**d)
1543 d["VOLUMETRIC_COMPUTE_KERNEL_CALL_STATELESS"] = jinja2.Template(
1544 self._volumetric_compute_kernel_call_stateless, undefined=jinja2.DebugUndefined
1545 ).render(**d)
1546 d["RIEMANN_COMPUTE_KERNEL_CALL"] = jinja2.Template(
1547 self._Riemann_compute_kernel_call, undefined=jinja2.DebugUndefined
1548 ).render(**d)
1549 d["RIEMANN_COMPUTE_KERNEL_CALL_STATELESS"] = jinja2.Template(
1550 self._Riemann_compute_kernel_call_stateless, undefined=jinja2.DebugUndefined
1551 ).render(**d)
1552 d["ADD_SOLVER_CONTRIBUTIONS_CALL"] = jinja2.Template(
1553 self._add_solver_contributions_call, undefined=jinja2.DebugUndefined
1554 ).render(**d)
1555 d["MULTIPLY_WITH_INVERTED_MASS_MATRIX_CALL"] = jinja2.Template(
1557 undefined=jinja2.DebugUndefined,
1558 ).render(**d)
1559
1560
1561 @property
1562 def unknowns(self):
1563 return self._unknowns
1564
1565
1566 @property
1568 return self._auxiliary_variables
1569
1570
1571 @unknowns.setter
1572 def unknowns(self, value):
1573 self._unknowns = value
1575 self.create_action_sets()
1576
1577
1578 @auxiliary_variables.setter
1579 def auxiliary_variables(self, value):
1580 self._auxiliary_variables = value
1582 self.create_action_sets()
1583
1584
1586 return RungeKutta_steps(self._rk_order)
1587
1588
1590 self,
1591 flux=None,
1592 ncp=None,
1593 eigenvalues=None,
1594 boundary_conditions=None,
1595 refinement_criterion=None,
1596 initial_conditions=None,
1597 source_term=None,
1598 point_source=None,
1599 additional_action_set_includes="",
1600 additional_user_includes="",
1601 ):
1602 """
1603 If you pass in User_Defined, then the generator will create C++ stubs
1604 that you have to befill manually. If you pass in None_Implementation, it
1605 will create nop, i.e., no implementation or defaults. Any other string
1606 is copied 1:1 into the implementation. If you pass in None, then the
1607 set value so far won't be overwritten.
1608 """
1609 if boundary_conditions is not None:
1610 self._boundary_conditions_implementation = boundary_conditions
1611 if refinement_criterion is not None:
1612 self._refinement_criterion_implementation = refinement_criterion
1613 if initial_conditions is not None:
1614 self._initial_conditions_implementation = initial_conditions
1615
1616 if refinement_criterion == exahype2.solvers.PDETerms.None_Implementation:
1617 assert False, "refinement criterion cannot be none"
1618 if initial_conditions == exahype2.solvers.PDETerms.None_Implementation:
1619 assert False, "initial conditions cannot be none"
1620
1621 if flux is not None:
1622 self._flux_implementation = flux
1623 if ncp is not None:
1624 self._ncp_implementation = ncp
1625 if eigenvalues is not None:
1626 self._eigenvalues_implementation = eigenvalues
1627 if source_term is not None:
1628 self._source_term_implementation = source_term
1629 if point_source is not None:
1630 self._point_sources_implementation = point_source
1631
1632 self.add_user_action_set_includes(additional_action_set_includes)
1633 self.add_user_solver_includes(additional_user_includes)
1634
1636 self.create_action_sets()
1637
1638
1639 @property
1644 @property
1649 @postprocess_updated_cell_after_Runge_Kutta_step.setter
1651 """
1652 Define a postprocessing routine over the data
1653
1654 This routine allows you to update the patch data immediately after the
1655 patch update. The routine provides the whole patch, and therefore you
1656 can write stuff similar to
1657
1658 my_solver.postprocess_updated_cell = " ""
1659 {
1660 constexpr int itmax = ({{DG_ORDER}} + 1) * ({{DG_ORDER}} + 1) * ({{DG_ORDER}} + 1);
1661 int index = 0;
1662 for (int i = 0; i < itmax; i++) {
1663 applications::exahype2::ccz4::enforceCCZ4constraints(dQdt + index);
1664 index += {{NUMBER_OF_UNKNOWNS}};
1665 }
1666 }
1667
1668
1669 ## Difference to Finite Volume solvers
1670
1671 Different to the Finite Volume solvers, you don't have the auxiliary
1672 variables in this routine. You also don't have the new solution, but the
1673 time derivative instead.
1674
1675
1676 ## Available constants
1677
1678 This list is not complete. You might want to consult the generated code to
1679 spot more variables that are on offer. Whenever I use the brackets below,
1680 these are symbolic constants which will be befilled with the proper constants
1681 when the postprocessing code snippet is inserted into the generated code.
1682 The other variables are directly available, i.e., no text replacement is done
1683 here.
1684
1685 - {{DG_ORDER}} Spatial discretisation order.
1686 - timeStamp
1687 - timeStepSize
1688 - dQdt
1689
1690 ## Runge-Kutta steps
1691
1692 This postprocessing is called after each and every Runge-Kutta step.
1693 """
1696 self.create_action_sets()
1697
1698
1699 @postprocess_updated_cell_after_final_linear_combination.setter
1701 """
1702 Define a postprocessing routine over the data
1703
1704 This routine allows you to update the patch data immediately after the
1705 patch update. The routine provides the whole patch, and therefore you
1706 can write stuff similar to
1707
1708 my_solver.postprocess_updated_cell = " ""
1709 {
1710 constexpr int itmax = ({{DG_ORDER}} +1 ) * ({{DG_ORDER}} + 1) * ({{DG_ORDER}} + 1);
1711 int index = 0;
1712 for (int i = 0; i < itmax; i++) {
1713 applications::exahype2::ccz4::enforceCCZ4constraints(dQdt + index);
1714 index += {{NUMBER_OF_UNKNOWNS}};
1715 }
1716 }
1717
1718
1719 ## Difference to Finite Volume solvers
1720
1721 Different to the Finite Volume solvers, you don't have the auxiliary
1722 variables in this routine. You also don't have the new solution, but the
1723 time derivative instead.
1724
1725
1726 ## Available constants
1727
1728 This list is not complete. You might want to consult the generated code to
1729 spot more variables that are on offer. Whenever I use the brackets below,
1730 these are symbolic constants which will be befilled with the proper constants
1731 when the postprocessing code snippet is inserted into the generated code.
1732 The other variables are directly available, i.e., no text replacement is done
1733 here.
1734
1735 - {{DG_ORDER}} Spatial discretisation order.
1736
1737
1738 ## Runge-Kutta steps
1739
1740 This postprocessing is called after each and every Runge-Kutta step.
1741 """
1744 self.create_action_sets()
1745
1746
1748 self,
1749 cell_data_storage: Storage,
1750 face_data_storage: Storage,
1751 ):
1752 """
1753 By default, we hold all data on the call stacks. You can explicitly switch
1754 to storage on the heap via smart pointers.
1755
1756 @see create_data_structures()
1757 """
1758 assert isinstance(cell_data_storage, Storage)
1759 assert isinstance(face_data_storage, Storage)
1760
1761 self._cell_data_storage = cell_data_storage
1762 self._face_data_storage = face_data_storage
1763
1765 self.create_action_sets()
Update the cell label within a sweep.
Definition CellLabel.py:9
An abstract class for any RKDG solver of any order.
_load_face_data_default_guard(self)
Extend the guard via ands only.
create_readme_descriptor(self, domain_offset, domain_size)
_store_cell_data_default_guard(self)
Extend the guard via ands only.
_load_cell_data_default_guard(self)
Extend the guard via ands only.
add_to_Peano4_datamodel(self, datamodel, verbose)
Add all required data to the Peano4 project's datamodel so it is properly built up.
add_actions_to_init_grid(self, step)
Add the action sets to the grid initialisation.
add_implementation_files_to_project(self, namespace, output, dimensions, subdirectory="")
The ExaHyPE project will call this operation when it sets up the overall environment.
_store_face_data_default_guard(self)
Extend the guard via ands only.
add_actions_to_plot_solution(self, step, output_path)
Dump snapshot of solution.
add_actions_to_create_grid(self, step, evaluate_refinement_criterion)
create_data_structures(self)
Recall in subclasses if you wanna change the number of unknowns or auxiliary variables.
create_action_sets(self)
Overwrite in subclasses if you wanna create different action sets.
set_implementation(self, flux=None, ncp=None, eigenvalues=None, boundary_conditions=None, refinement_criterion=None, initial_conditions=None, source_term=None, point_source=None, additional_action_set_includes="", additional_user_includes="")
If you pass in User_Defined, then the generator will create C++ stubs that you have to befill manuall...
user_action_set_includes(self)
Add further includes to this property, if your action sets require some additional routines from othe...
user_solver_includes(self)
Add further includes to this property, if your solver requires some additional routines from other he...
switch_storage_scheme(self, Storage cell_data_storage, Storage face_data_storage)
By default, we hold all data on the call stacks.
_init_dictionary_with_default_parameters(self, d)
This one is called by all algorithmic steps before I invoke add_entries_to_text_replacement_dictionar...
add_user_action_set_includes(self, value)
Add further includes to this property, if your action sets require some additional routines from othe...
__init__(self, name, rk_order, polynomial_basis, FaceProjections face_projections, unknowns, auxiliary_variables, min_cell_h, max_cell_h, plot_grid_properties, bool pde_terms_without_state, baseline_action_set_descend_invocation_order=0)
Constructor of the Runge-Kutta Discontinuous Galerkin solver.
add_actions_to_perform_time_step(self, step)
The tricky part here is that we have to get the order right.
add_use_data_statements_to_Peano4_solver_step(self, step)
Tell Peano what data to move around.
set_plot_description(self, description)
Use this one to set a description within the output patch file that tells the vis solver what the sem...
add_user_solver_includes(self, value)
Add further includes to this property, if your solver requires some additional routines from other he...
The linear combination of the Runge Kutta trials has to be projected onto the faces,...
Action set determining the final solution using a linear combination of previous Runge-Kutta guesses.
Default DG interpolation/restriction action set.
Definition DynamicAMR.py:13
The linear combination of the Runge Kutta trials has to be projected onto the faces,...
Computes a linear combination of all of the estimates according to the Butcher Tableau.
PostprocessSolution differs from other action sets, as I only create it once.
PreprocessSolution differs from other action sets, as I only create it once.
Solve the actual Riemann problem and add stuff back to the solution.
Very simple converter which maps the patch 1:1 onto a double array.
Action set (reactions to events)
Definition ActionSet.py:6