3from .CellCenteredFiniteDifferences
import CellCenteredFiniteDifferences
12 ReconstructPatchAndApplyFunctor,
21 Update one cell, i.e. compute Runge-Kutta step on it
23 This routine is very similar to the Finite Volume step. The big difference in
24 the context of Runge-Kutta is that we always have to construct a linear
25 combination of the input data and then write the new estimate for the
26 right-hand side into the large estimate field.
28 So we operate in two steps: First, we let the default block-structured code
29 reconstruct the old time step data. After that, we add the rhs estimates
30 onto this reconstructed data. The latter can be read as a correction step
31 to the reconstructed values. It is correct if and only if the haloes are
32 properly initialised with a corrected value, as we cannot correct the halo
36 # Correction terms with Runge-Kutta trials
38 We current the values with a linear combination of all of the estimates according to the
41 The linear combination is a volumetric representation which includes both
42 the degrees of freedom and the auxiliary variables. However, the auxiliary
43 variables are not developing over time. In Runge-Kutta, I have estimates
44 for the right-hand side, i.e. for the derivatives
46 @f$ \partial _t Q = F(Q) @f$
48 This is the stuff stored in RhsEstimates. It does not comprise any auxiliary
49 variables. So I have to copy the auxiliary variables from the last valid
50 time step every time I reconstruct a new guess.
56 SolveRiemannProblemsOverPatch = jinja2.Template(
58 double timeStamp = fineGridCell{{SOLVER_NAME}}CellLabel.getTimeStamp();
61 // double timeStepSize
62 {{COMPUTE_TIME_STEP_SIZE}}
64 {% for PREDICATE_NO in range(0,PREDICATES|length) %}
65 if ({{PREDICATES[PREDICATE_NO]}}) {
66 ::exahype2::enumerator::AoSLexicographicEnumerator enumeratorWithAuxiliaryVariablesOnReconstructedPatch( 1, {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}}, {{HALO_SIZE}}, {{NUMBER_OF_UNKNOWNS}}, {{NUMBER_OF_AUXILIARY_VARIABLES}});
67 ::exahype2::enumerator::AoSLexicographicEnumerator enumeratorWithoutAuxiliaryVariables( {{RK_STEPS}}, {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}}, 0, {{NUMBER_OF_UNKNOWNS}}, 0 );
68 int primaryVarsIndices[25]={ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
69 10,11,12,13,14,15,16,17,18,19,
71 int auxiliaryVarsIndices[34]={23,24,25,26,27,28,29,30,31,32,
72 33,34,35,36,37,38,39,40,41,42,
73 43,44,45,46,47,48,49,50,51,52,
76 dfor( dof, {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}} ) {
77 {% if RECONSTRUCTION_STAGES==0 %}
78 for (int unknown=0; unknown<{{NUMBER_OF_UNKNOWNS}}; unknown++) {
79 {% for WEIGHT_NO in range(0,BUTCHER_TABLEAU_WEIGHTS[PREDICATE_NO]|length) %}
80 {% if BUTCHER_TABLEAU_WEIGHTS[PREDICATE_NO][WEIGHT_NO]!=0 %}
81 oldQWithHalo[ enumeratorWithAuxiliaryVariablesOnReconstructedPatch(0,dof,unknown) ] +=
82 timeStepSize * {{BUTCHER_TABLEAU_WEIGHTS[PREDICATE_NO][WEIGHT_NO]}} *
83 fineGridCell{{UNKNOWN_IDENTIFIER}}RhsEstimates.value[ enumeratorWithoutAuxiliaryVariables({{WEIGHT_NO}},dof,unknown) ];
88 for (int unknown : primaryVarsIndices) {
89 {% for WEIGHT_NO in range(0,BUTCHER_TABLEAU_WEIGHTS[PREDICATE_NO]|length) %}
90 {% if BUTCHER_TABLEAU_WEIGHTS[PREDICATE_NO][WEIGHT_NO]!=0 %}
91 oldQWithHalo[ enumeratorWithAuxiliaryVariablesOnReconstructedPatch(0,dof,unknown) ] +=
92 timeStepSize * {{BUTCHER_TABLEAU_WEIGHTS[PREDICATE_NO][WEIGHT_NO]}} *
93 fineGridCell{{UNKNOWN_IDENTIFIER}}RhsEstimates.value[ enumeratorWithoutAuxiliaryVariables({{WEIGHT_NO}},dof,unknown) ];
97 for (int unknown : auxiliaryVarsIndices) {
98 oldQWithHalo[ enumeratorWithAuxiliaryVariablesOnReconstructedPatch(0,dof,unknown) ] = //abuse of RHS, auxiliary variables are stored there.
99 fineGridCell{{UNKNOWN_IDENTIFIER}}RhsEstimates.value[ enumeratorWithoutAuxiliaryVariables({{PREDICATE_NO}},dof,unknown) ];
105 newQ = fineGridCell{{UNKNOWN_IDENTIFIER}}RhsEstimates.value + enumeratorWithoutAuxiliaryVariables({{PREDICATE_NO}},0,0);
109 {% for PREDICATE_NO in range(0,PREDICATES_RECONSTRUCTION|length) %}
110 if ({{PREDICATES_RECONSTRUCTION[PREDICATE_NO]}}) {
111 ::exahype2::enumerator::AoSLexicographicEnumerator enumeratorWithAuxiliaryVariablesOnReconstructedPatch( 1, {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}}, {{HALO_SIZE}}, {{NUMBER_OF_UNKNOWNS}}, {{NUMBER_OF_AUXILIARY_VARIABLES}});
112 ::exahype2::enumerator::AoSLexicographicEnumerator enumeratorWithoutAuxiliaryVariables( {{RK_STEPS}}, {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}}, 0, {{NUMBER_OF_UNKNOWNS}}, 0 );
113 int primaryVarsIndices[25]={ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
114 10,11,12,13,14,15,16,17,18,19,
116 int auxiliaryVarsIndices[34]={23,24,25,26,27,28,29,30,31,32,
117 33,34,35,36,37,38,39,40,41,42,
118 43,44,45,46,47,48,49,50,51,52,
121 dfor( dof, {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}} ) {
122 for (int unknown : primaryVarsIndices) {
123 {% for WEIGHT_NO in range(0,BUTCHER_TABLEAU_WEIGHTS[PREDICATE_NO]|length) %}
124 {% if BUTCHER_TABLEAU_WEIGHTS[PREDICATE_NO][WEIGHT_NO]!=0 %}
125 oldQWithHalo[ enumeratorWithAuxiliaryVariablesOnReconstructedPatch(0,dof,unknown) ] +=
126 timeStepSize * {{BUTCHER_TABLEAU_WEIGHTS[PREDICATE_NO][WEIGHT_NO]}} *
127 fineGridCell{{UNKNOWN_IDENTIFIER}}RhsEstimates.value[ enumeratorWithoutAuxiliaryVariables({{WEIGHT_NO}},dof,unknown) ];
131 for (int unknown : auxiliaryVarsIndices) {
132 oldQWithHalo[ enumeratorWithAuxiliaryVariablesOnReconstructedPatch(0,dof,unknown) ] = 0;
133 // actually no need for this initialization as we do overwritten in reconstruction, but we keep it for safety.
138 newQ = fineGridCell{{UNKNOWN_IDENTIFIER}}RhsEstimates.value + enumeratorWithoutAuxiliaryVariables({{PREDICATE_NO}},0,0);
143 {{PREPROCESS_RECONSTRUCTED_PATCH}}
145 assertion2( tarch::la::greaterEquals( timeStamp, 0.0 ), timeStamp, timeStepSize );
146 assertion2( tarch::la::greaterEquals( timeStepSize, 0.0 ), timeStamp, timeStepSize );
148 ::exahype2::fd::validatePatch(
150 {{NUMBER_OF_UNKNOWNS}},
151 {{NUMBER_OF_AUXILIARY_VARIABLES}},
152 {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}},
153 {{HALO_SIZE}}, // halo
154 std::string(__FILE__) + "(" + std::to_string(__LINE__) + "): " + marker.toString()
155 ); // previous time step has to be valid
157 double subTimeStamp=timeStamp;
158 {% for PREDICATE_NO in range(0,PREDICATES|length) %}
159 if ({{PREDICATES[PREDICATE_NO]}}) {
160 subTimeStamp += {{BUTCHER_TABLEAU_RELATIVE_TIME_STEP_SIZES[PREDICATE_NO]}}*timeStepSize;
164 ::exahype2::CellData patchData( oldQWithHalo, marker.x(), marker.h(), subTimeStamp, timeStepSize, newQ );
166 {% if RECONSTRUCTION_STAGES==0 %}
167 ::exahype2::fd::{{KERNEL_NAMESPACE}}::{{COMPUTE_KERNEL_CALL}}
169 if ({{IS_UNDER_RECONSTRUCTION}}) {
170 ::exahype2::fd::fd4::reconstruct_first_derivatives(
172 {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}}, // int numberOfGridCellsPerPatchPerAxis,
174 {{NUMBER_OF_UNKNOWNS}},
175 {{NUMBER_OF_AUXILIARY_VARIABLES}}
177 } else if ({{IS_UNDER_RKSTEP}}) {
178 ::exahype2::fd::{{KERNEL_NAMESPACE}}::{{COMPUTE_KERNEL_CALL}}
182 ::exahype2::fd::validatePatch(
184 {{NUMBER_OF_UNKNOWNS}},
185 {{NUMBER_OF_AUXILIARY_VARIABLES}},
186 {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}},
188 std::string(__FILE__) + "(" + std::to_string(__LINE__) + "): " + marker.toString()
189 ); // outcome has to be valid
195 ReconstructPatchAndApplyFunctor.__init__(
198 patch_overlap=solver._patch_overlap_new,
199 functor_implementation=
"""
200#error please switch to your Riemann solver of choice
202 reconstructed_array_memory_location=solver._reconstructed_array_memory_location,
203 guard=
"not marker.willBeRefined() and not marker.hasBeenRefined()",
204 add_assertions_to_halo_exchange=
False,
212 +
"""CellLabel.setHasUpdated(false);
222 This is our plug-in point to alter the underlying dictionary
227 self.
_solver._init_dictionary_with_default_parameters(d)
228 self.
_solver.add_entries_to_text_replacement_dictionary(d)
231 "{} and repositories::{}.getSolverState()=={}::SolverState::RungeKuttaSubStep{}".format(
232 self.
_solver._store_cell_data_default_guard(),
233 self.
_solver.get_name_of_global_instance(),
237 for step
in range(0, self.
_solver.number_of_Runge_Kutta_steps())
239 if self.
_solver._reconstruction_stages != 0 :
240 d[
"PREDICATES_RECONSTRUCTION"] = [
241 "{} and repositories::{}.getSolverState()=={}::SolverState::ReconstructionSubStep{}".format(
242 self.
_solver._store_cell_data_default_guard(),
243 self.
_solver.get_name_of_global_instance(),
247 for step
in range(0, self.
_solver._reconstruction_stages)
250 re_condition_text=
"repositories::{}.getSolverState()=={}::SolverState::ReconstructionSubStep{}".format(
251 self.
_solver.get_name_of_global_instance(),
255 for step
in range(1,self.
_solver._reconstruction_stages):
256 re_condition_text +=
" or repositories::{}.getSolverState()=={}::SolverState::ReconstructionSubStep{}".format(
257 self.
_solver.get_name_of_global_instance(),
261 rk_condition_text=
"repositories::{}.getSolverState()=={}::SolverState::RungeKuttaSubStep{}".format(
262 self.
_solver.get_name_of_global_instance(),
266 for step
in range(1,self.
_solver.number_of_Runge_Kutta_steps()):
267 rk_condition_text +=
" or repositories::{}.getSolverState()=={}::SolverState::RungeKuttaSubStep{}".format(
268 self.
_solver.get_name_of_global_instance(),
272 d[
"IS_UNDER_RECONSTRUCTION"]=re_condition_text
273 d[
"IS_UNDER_RKSTEP"]=rk_condition_text
278 "BUTCHER_TABLEAU_RELATIVE_TIME_STEP_SIZES"
290#include "exahype2/enumerator/enumerator.h"
291#include "exahype2/fd/PatchUtils.h"
292#include "tarch/NonCriticalAssertions.h"
294 + self.
_solver._get_default_includes()
295 + self.
_solver.user_action_set_includes
301 You should replicate this function in each subclass, so you get
302 meaningful action set names (otherwise, it will be
303 AbstractRKFDActionSet0,1,2,...).
306 return __name__.replace(
".py",
"").replace(
".",
"_") +
"_UpdateCell"
311 Probably the simplest solver you could think off.
313 This particular variant works only for Runge-Kutta order greater or equal
317 :: Write your own specialisation ::
319 self._preprocess_reconstructed_patch
320 Has to hold any preprocessing, but it also has to set the doubles
321 timeStepSize and timeStamp to valid data.
323 self._postprocess_updated_patch
324 You don't have to redefine this one, but if you want to alter the
325 time step size, then this is the right place to do so. If you don't
326 alter timeStepSize, the code will automatically continue with
327 the current one subject to a preprocessing in the next step.
341 plot_grid_properties,
346 Instantiate a generic FV scheme with an overlap of 1.
351 super(SeparateSweeps, self).
__init__(
360 plot_grid_properties,
366 "Runge-Kutta order has to be at least 1 but was {}".format(rk_order)
379 Call the superclass' create_data_structures() to ensure that all the data
380 structures are in place, i.e. each cell can host a patch, that each face hosts
381 patch overlaps, and so forth. These quantities are all set to defaults. See
382 FV.create_data_structures().
384 After that, take the patch overlap (that's the data stored within the faces)
385 and ensure that these are sent and received via MPI whenever they are also
386 stored persistently. The default in FV is that no domain boundary data exchange
387 is active. Finally, ensure that the old data is only exchanged between the
388 initialisation sweep and the first first grid run-through.
393 initialisation_sweep_guard =
"""(
394 repositories::{}.getSolverState()=={}::SolverState::GridInitialisation
398 first_iteration_after_initialisation_guard =
"""(
399 repositories::{}.getSolverState()=={}::SolverState::RungeKuttaSubStep0AfterGridInitialisation or
400 repositories::{}.getSolverState()=={}::SolverState::PlottingAfterGridInitialisation
413 first_iteration_after_initialisation_guard
417 first_sweep_of_time_step_or_plotting_guard =
"""(
418 repositories::{}.getSolverState()=={}::SolverState::RungeKuttaSubStep0AfterGridInitialisation or
419 repositories::{}.getSolverState()=={}::SolverState::ReconstructionSubStep0 or
420 repositories::{}.getSolverState()=={}::SolverState::PlottingAfterGridInitialisation or
421 repositories::{}.getSolverState()=={}::SolverState::Plotting or
422 repositories::{}.getSolverState()=={}::SolverState::Suspended
431 first_sweep_of_time_step_or_plotting_guard =
"""(
432 repositories::{}.getSolverState()=={}::SolverState::RungeKuttaSubStep0AfterGridInitialisation or
433 repositories::{}.getSolverState()=={}::SolverState::RungeKuttaSubStep0 or
434 repositories::{}.getSolverState()=={}::SolverState::PlottingAfterGridInitialisation or
435 repositories::{}.getSolverState()=={}::SolverState::Plotting or
436 repositories::{}.getSolverState()=={}::SolverState::Suspended
446 last_sweep_of_time_step_or_plotting_initialisation =
"""(
447 repositories::{}.getSolverState()=={}::SolverState::GridInitialisation or
448 repositories::{}.getSolverState()=={}::SolverState::RungeKuttaSubStep{} or
449 repositories::{}.getSolverState()=={}::SolverState::PlottingAfterGridInitialisation or
450 repositories::{}.getSolverState()=={}::SolverState::Plotting or
451 repositories::{}.getSolverState()=={}::SolverState::Suspended
460 self.
_patch_estimates.generator.load_store_compute_flag =
"::peano4::grid::constructLoadStoreComputeFlag({},{},{})".format(
465 + first_sweep_of_time_step_or_plotting_guard
470 + last_sweep_of_time_step_or_plotting_initialisation
477 Call superclass routine and then reconfigure the update cell call.
478 Only the UpdateCell action set is specific to a single sweep.
480 This operation is implicity called via the superconstructor.
482 Each guard here should combine the storage predicate with the
495 "{} and repositories::{}.getSolverState()=={}::SolverState::RungeKuttaSubStep{}".format(
503 "{} and repositories::{}.getSolverState()=={}::SolverState::GridInitialisation".format(
512 "{} and repositories::{}.getSolverState()=={}::SolverState::ReconstructionSubStep{}".format(
525#include "exahype2/CellData.h"
527 + super(SeparateSweeps, self).user_action_set_includes
537 refinement_criterion,
540 additional_action_set_includes,
541 additional_user_includes,
544 If you pass in User_Defined, then the generator will create C++ stubs
545 that you have to befill manually. If you pass in None_Implementation, it
546 will create nop, i.e. no implementation or defaults. Any other string
547 is copied 1:1 into the implementation. If you pass in None, then the
548 set value so far won't be overwritten.
550 if boundary_conditions
is not None:
552 if refinement_criterion
is not None:
554 if initial_conditions
is not None:
556 if memory_location
is not None:
559 if refinement_criterion == exahype2.solvers.PDETerms.None_Implementation:
560 assert False,
"refinement criterion cannot be none"
561 if initial_conditions == exahype2.solvers.PDETerms.None_Implementation:
562 assert False,
"initial conditions cannot be none"
568 if eigenvalues
is not None:
570 if source_term
is not None:
575 == peano4.toolbox.blockstructured.ReconstructedArrayMemoryLocation.HeapThroughTarchWithoutDelete
577 == peano4.toolbox.blockstructured.ReconstructedArrayMemoryLocation.HeapWithoutDelete
579 == peano4.toolbox.blockstructured.ReconstructedArrayMemoryLocation.ManagedSharedAcceleratorDeviceMemoryThroughTarchWithoutDelete
582 "memory mode without appropriate delete chosen, i.e. this will lead to a memory leak"
592 d: Dictionary of string to string
Abstract solver for patch-based finite diffences.
_provide_cell_data_to_compute_kernels_default_guard(self)
_action_set_project_patch_onto_faces
_user_action_set_includes
_store_cell_data_default_guard(self)
Extend the guard via ands only.
_load_cell_data_default_guard(self)
Extend the guard via ands only.
get_name_of_global_instance(self)
number_of_Runge_Kutta_steps(self)
Return number of steps required to realise the Runge-Kutta scheme.
_solver_template_file_class_name
_baseline_action_set_descend_invocation_order
_boundary_conditions_implementation
_reconstructed_array_memory_location
create_action_sets(self)
Create required action sets.
_action_set_compute_final_linear_combination
_initial_conditions_implementation
_refinement_criterion_implementation
Probably the simplest solver you could think off.
add_entries_to_text_replacement_dictionary(self, d)
d: Dictionary of string to string in/out argument
_initial_conditions_implementation
__init__(self, name, patch_size, overlap, rk_order, unknowns, auxiliary_variables, min_meshcell_h, max_meshcell_h, plot_grid_properties, kernel_namespace)
Instantiate a generic FV scheme with an overlap of 1.
_reconstructed_array_memory_location
set_implementation(self, flux, ncp, source_term, eigenvalues, boundary_conditions, refinement_criterion, initial_conditions, memory_location, 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...
_eigenvalues_implementation
user_action_set_includes(self)
Add further includes to this property, if your action sets require some additional routines from othe...
_refinement_criterion_implementation
_boundary_conditions_implementation
create_data_structures(self)
Call the superclass' create_data_structures() to ensure that all the data structures are in place,...
create_action_sets(self)
Call superclass routine and then reconfigure the update cell call.
_source_term_implementation
_solver_template_file_class_name
get_action_set_name(self)
You should replicate this function in each subclass, so you get meaningful action set names (otherwis...
_add_action_set_entries_to_dictionary(self, d)
This is our plug-in point to alter the underlying dictionary.
get_includes(self)
Return include statements that you need.
_Template_TouchCellFirstTime_Preamble
SolveRiemannProblemsOverPatch
__init__(self, solver)
patch: peano4.datamodel.Patch Patch which is to be used