Peano 4
Loading...
Searching...
No Matches
ComputeFirstDerivatives.py
Go to the documentation of this file.
2# New action set for second order CCZ4 formualtion
3#
4import peano4
5import jinja2
6
7from exahype2.solvers.ButcherTableau import ButcherTableau
8from exahype2.solvers.rkfd.OneSweepPerRungeKuttaStep import ReconstructLinearCombination
9
10
12 """!
13
14 Explicit reconstruction of derivatives for Finite Volumes
15
16 """
17 def __init__(self,
18 solver
19 ):
20 super( ComputeFirstDerivativesFD4, self ).__init__(patch=solver._patch,
21 patch_overlap=solver._patch_overlap_new,
22 functor_implementation="""
23::exahype2::CellData patchData(
24 oldQWithHalo,
25 marker.x(), marker.h(),
26 0.0, // timeStamp,
27 0.0, // timeStepSize,
28 newQ
29);
30::exahype2::fd::fd4::reconstruct_first_derivatives(patchData,"""+str(solver._patch.dim[0])+","+str(int(solver._patch_overlap_new.dim[0]/2))+","+str(solver._patch.no_of_unknowns)+""",0);
31""",
32 reconstructed_array_memory_location=solver._reconstructed_array_memory_location,
33 guard="not marker.willBeRefined() and not marker.hasBeenRefined()",
34 add_assertions_to_halo_exchange=False,
35 )
36 pass
37
38
40 return False
41
42
43 def get_includes(self):
44 return super( ComputeFirstDerivativesFV, self ).get_includes() + """
45 #include "exahype2/CellData.h"
46 #include "tarch/NonCriticalAssertions.h"
47 """
48
49
51 return __name__.replace(".py", "").replace(".", "_")
52
53
54
56 """!
57
58 Explicit reconstruction of first derivatives for FD4 discretisation
59
60 FD4 is implemented as Runge-Kutta scheme. So we first have to recompute
61 the current solution guess according to the Butcher tableau. Then we feed
62 this reconstruction into the derivative computation. In theory, this is all
63 simple, but in practice things require more work.
64
65 We assume that developers follow @ref page_exahype_coupling "ExaHyPE's generic recommendation how to add additional traversals".
66 Yet, this is only half of the story. If they follow the vanilla blueprint
67 and plug into a step after the last time step
68
69 ~~~~~~~~~~~~~~~~~~~~~~~~~~
70 if (repositories::isLastGridSweepOfTimeStep()
71 and
72 repositories::StepRepository::toStepEnum( peano4::parallel::Node::getInstance().getCurrentProgramStep() ) != repositories::StepRepository::Steps::AdditionalMeshTraversal
73 ) {
74 peano4::parallel::Node::getInstance().setNextProgramStep(
75 repositories::StepRepository::toProgramStep( repositories::StepRepository::Steps::AdditionalMeshTraversal )
76 );
77 continueToSolve = true;
78 }
79 ~~~~~~~~~~~~~~~~~~~~~~~~~~
80
81 then the reconstruction is only done after the last and final Runge-Kutta
82 step. Alternatively, users might plug into the solver after each
83 Runge-Kutta step. In principle, we might say "so what", as all the
84 right-hand sides after the final step are there, so we can reconstruct
85 the final solution. However, there are two pitfalls:
86
87 1. A solver might not hold the rhs estimates persistently in-between two
88 time steps.
89 2. The solver's state always is toggled in startTimeStep() on the solver
90 instance's class. Therefore, an additional grid sweep after the last
91 Runge-Kutta sweep cannot be distinguished from an additional sweep
92 after the whole time step. They are the same.
93
94 Consequently, we need special treatment for the very last Runge-Kutta step:
95
96 - If we have an additional grid sweep after an intermediate RK step, we use
97 the right-hand side estimate of @f$ \partial _t Q = F(Q) @f$ to store the
98 outcome. Furthermore, we scale the derivative with 1/dt, as the linear
99 combination of RK will later on multiple this estimate with dt again.
100 - If we are in the very last grid sweep, we work with the current estimate.
101
102 The actual derviative calculation is "outsourced" to the routine
103
104 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
105 ::exahype2::CellData patchData(
106 oldQWithHalo,
107 marker.x(), marker.h(),
108 0.0, // timeStamp => not used here
109 0.0, // timeStepSize => not used here
110 newQ
111 );
112 ::exahype2::fd::fd4::reconstruct_first_derivatives(
113 patchData,
114 {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}}, // int numberOfGridCellsPerPatchPerAxis,
115 {{HALO_SIZE}},
116 {{NUMBER_OF_UNKNOWNS}},
117 {{NUMBER_OF_AUXILIARY_VARIABLES}}
118 );
119 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
120
121 Wrapping up oldQWithHalo and newQ in a CellData object is a slight
122 overkill (we could have done with only passing the two pointers, as we
123 don't use the time step size or similar anyway), but we follow the
124 general ExaHyPE pattern here - and in theory could rely on someone else
125 later to deploy this to GPUs, e.g.
126
127 Right after the reconstruction, we have to project the outcome back again
128 onto the faces. Here, we again have to take into account that we work
129 with a Runge-Kutta scheme. The cell hosts N blocks of cell data for the
130 N shots of Runge-Kutta. The reconstruction writes into the nth block
131 according to the current step. Consequently, the projection also has to
132 work with the nth block. exahype2.solvers.rkfd.actionsets.ProjectPatchOntoFaces
133 provides all the logic for this, so we assume that users add this one to
134 their script. As the present action set plugs into touchCellFirstTime()
135 and the projection onto the faces plugs into touchCellLastTime(), the
136 order is always correct no matter how you prioritise between the different
137 action sets.
138
139 Further to the projection, we also have to roll over details. Again,
140 exahype2.solvers.rkfd.actionsets.ProjectPatchOnto has more documentation
141 on this.
142
143 """
144 def __init__(self,
145 solver,
146 is_enclave_solver
147 ):
148 super( ComputeFirstDerivativesFD4RK, self ).__init__(patch=solver._patch,
149 patch_overlap=solver._patch_overlap_new,
150 functor_implementation="""
151#error Will be inserted later down in this Python file
152""",
153 reconstructed_array_memory_location=solver._reconstructed_array_memory_location,
154 guard="not marker.willBeRefined() and not marker.hasBeenRefined()",
155 add_assertions_to_halo_exchange=False,
156 )
157 self._solver = solver
159 self._use_enclave_solver = is_enclave_solver
160 pass
161
162
163 ComputeDerivativesOverCell = jinja2.Template( """
164 double timeStamp = fineGridCell{{SOLVER_NAME}}CellLabel.getTimeStamp();
165
166 // Set the variable
167 // double timeStepSize
168 {{COMPUTE_TIME_STEP_SIZE}}
169
170 {% for PREDICATE_NO in range(0,PREDICATES|length) %}
171 if ({{PREDICATES[PREDICATE_NO]}}) {
172 ::exahype2::enumerator::AoSLexicographicEnumerator enumeratorWithAuxiliaryVariablesOnReconstructedPatch( 1, {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}}, {{HALO_SIZE}}, {{NUMBER_OF_UNKNOWNS}}, {{NUMBER_OF_AUXILIARY_VARIABLES}});
173 ::exahype2::enumerator::AoSLexicographicEnumerator enumeratorWithoutAuxiliaryVariables( {{RK_STEPS}}, {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}}, 0, {{NUMBER_OF_UNKNOWNS}}, 0 );
174
175 dfor( dof, {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}} ) {
176 for (int unknown=0; unknown<{{NUMBER_OF_UNKNOWNS}}; unknown++) {
177 {% for WEIGHT_NO in range(0,BUTCHER_TABLEAU_WEIGHTS[PREDICATE_NO]|length) %}
178 {% if BUTCHER_TABLEAU_WEIGHTS[PREDICATE_NO][WEIGHT_NO]!=0 %}
179 oldQWithHalo[ enumeratorWithAuxiliaryVariablesOnReconstructedPatch(0,dof,unknown) ] +=
180 {{BUTCHER_TABLEAU_RELATIVE_TIME_STEP_SIZES[PREDICATE_NO]}} * timeStepSize * {{BUTCHER_TABLEAU_WEIGHTS[PREDICATE_NO][WEIGHT_NO]}} *
181 fineGridCell{{UNKNOWN_IDENTIFIER}}RhsEstimates.value[ enumeratorWithoutAuxiliaryVariables({{PREDICATE_NO-1}},dof,unknown) ];
182 {% endif %}
183 {% endfor %}
184 }
185 }
186
187 newQ = fineGridCell{{UNKNOWN_IDENTIFIER}}RhsEstimates.value + enumeratorWithoutAuxiliaryVariables({{PREDICATE_NO}},0,0);
188 }
189 {% endfor %}
190
191 assertion2( tarch::la::greaterEquals( timeStamp, 0.0 ), timeStamp, timeStepSize );
192 assertion2( tarch::la::greaterEquals( timeStepSize, 0.0 ), timeStamp, timeStepSize );
193
194 ::exahype2::fd::validatePatch(
195 oldQWithHalo,
196 {{NUMBER_OF_UNKNOWNS}},
197 {{NUMBER_OF_AUXILIARY_VARIABLES}},
198 {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}},
199 {{HALO_SIZE}}, // halo
200 std::string(__FILE__) + "(" + std::to_string(__LINE__) + "): " + marker.toString()
201 ); // previous time step has to be valid
202 """ + """
203 // take into account that linear combination has already been computed,
204 // and rhs even might not be held persistently
205 if ( repositories::isLastGridSweepOfTimeStep() ) {
206 newQ = fineGridCell{{UNKNOWN_IDENTIFIER}}.value;
207 }
208
209 ::exahype2::CellData patchData(
210 oldQWithHalo,
211 marker.x(), marker.h(),
212 0.0, // timeStamp => not used here
213 0.0, // timeStepSize => not used here
214 newQ
215 );
216 ::exahype2::fd::fd4::reconstruct_first_derivatives(
217 patchData,
218 {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}}, // int numberOfGridCellsPerPatchPerAxis,
219 {{HALO_SIZE}},
220 {{NUMBER_OF_UNKNOWNS}},
221 {{NUMBER_OF_AUXILIARY_VARIABLES}}
222 );
223
224 if ( not repositories::isLastGridSweepOfTimeStep() ) {
225 ::exahype2::enumerator::AoSLexicographicEnumerator enumerator( 1, {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}}, 0, {{NUMBER_OF_UNKNOWNS}}, 0 );
226 dfor( dof, {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}} ) {
227 for (int unknown=0; unknown<{{NUMBER_OF_UNKNOWNS}}; unknown++) {
228 newQ[ enumerator(0,dof,unknown) ] *= 1.0 / timeStepSize;
229 }
230 }
231 }
232 """)
233
234
235# """+str(solver._patch.dim[0])+","+str(int(solver._patch_overlap_new.dim[0]/2))+","+str(solver._patch.no_of_unknowns)+"""
236
237
239 return False
240
241
242 def get_includes(self):
243 return super( ComputeFirstDerivativesFD4RK, self ).get_includes() + """
244 #include "exahype2/CellData.h"
245 #include "exahype2/fd/fd4/FD4.h"
246 #include "tarch/NonCriticalAssertions.h"
247 """
248
250 """!
251
252 This is our plug-in point to alter the underlying dictionary
253
254 """
255 super(ComputeFirstDerivativesFD4RK, self)._add_action_set_entries_to_dictionary(d)
256
257 self._solver._init_dictionary_with_default_parameters(d)
258 self._solver.add_entries_to_text_replacement_dictionary(d)
259
260 d["BUTCHER_TABLEAU_WEIGHTS"] = self._butcher_tableau.weight_matrix()
261 d["BUTCHER_TABLEAU_RELATIVE_TIME_STEP_SIZES"] = self._butcher_tableau.time_step_sizes()
262
263 if self._use_enclave_solver:
264 d["PREDICATES"] = self._solver._primary_sweeps_of_Runge_Kutta_step_on_cell
265 else:
266 d["PREDICATES"] = [
267 "{} and repositories::{}.getSolverState()=={}::SolverState::RungeKuttaSubStep{}".format(
268 self._solver._store_cell_data_default_guard(),
269 self._solver.get_name_of_global_instance(),
270 self._solver._name,
271 step,
272 )
273 for step in range(0, self._solver.number_of_Runge_Kutta_steps())
274 ]
275
276 # Has to come after we've set the predicates, as we use these
277 # fields in here already
278 d["CELL_FUNCTOR_IMPLEMENTATION"] = self.ComputeDerivativesOverCell.render(
279 **d
280 )
281
282
284 return __name__.replace(".py", "").replace(".", "_")
285
Explicit reconstruction of first derivatives for FD4 discretisation.
get_includes(self)
Return include statements that you need.
__init__(self, solver, is_enclave_solver)
patch: peano4.datamodel.Patch Patch which is to be used
_add_action_set_entries_to_dictionary(self, d)
This is our plug-in point to alter the underlying dictionary.
user_should_modify_template(self)
Is the user allowed to modify the output.
Explicit reconstruction of derivatives for Finite Volumes.
get_includes(self)
Return include statements that you need.
__init__(self, solver)
patch: peano4.datamodel.Patch Patch which is to be used
get_action_set_name(self)
Return unique action set name.
user_should_modify_template(self)
Is the user allowed to modify the output.