Peano 4
Loading...
Searching...
No Matches
ProjectPatchOntoFaces.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 peano4
4import jinja2
5
6from .AbstractRKFDActionSet import AbstractRKFDActionSet
7
8from peano4.solversteps.ActionSet import ActionSet
9
10from exahype2.solvers.ButcherTableau import ButcherTableau
11
12
14 """!
15
16 Project patch data onto faces, so the faces hold valid data which can feed into the subsequent Runge-Kutta step
17
18 This class is based upon peano4.toolbox.blockstructured.ProjectPatchOntoFaces.
19 Yet, it differs fundamentally in the way the inner-most loop is realised:
20 The standard projection from the toolbox takes a solution from a patch
21 and copies parts of this solution onto the faces. This is not what we want
22 to have here: we have the solution plus a number of estimates, and we want
23 to write exactly the type of Runge-Kutta linear combination onto the face
24 that we need to set the proper halo data on the adjacent cells in the next
25 grid sweep. Therefore, this class expects the Butcher tableau and first
26 computes the Runge-Kutta linear combination. The outcome of this is then
27 projected onto the faces.
28
29 The routine is always called as one of the last things, i.e. we have a
30 valid estimate in. To some degree, we anticipate what will happening next
31 in the Runge-Kutta scheme.
32
33 Once more: Each cell hosts the p predictions created by Runge-Kutta of the
34 pth order. These are stored in the field rhs estimates. The projection
35 combines these predictions using the Butcher tableau and projects the
36 linear combination onto the face. The action set never ever uses what's
37 stored in the value array. It always computes this linear combination.
38
39 ## Auxiliary parameters
40
41 Auxiliary parameters are not subject to the PDE in ExaHyPE's jargon. They
42 typically carry material parameters or similar. We project them onto the
43 faces, but here we don't do linear combinations due to the Butcher
44 tableau, as the rhs estimates do not hold any auxiliary parameters at all.
45
46 ## Last Runge-Kutta step
47
48 The last Runge-Kutta step is slightly different. After the last
49 Runge-Kutta estimate, the solver combines all the estimates into a new
50 solution. Consult the action set ComputeFinalLinearCombination for details.
51 After this has happened, we could take the plain solution stored within the
52 cell and write this guy onto the faces.
53
54 Instead, we ignore the face that the action set ComputeFinalLinearCombination
55 would do all of this. Instead, we compute the final linear combination
56 from the rhs guesses (redundantly) and project the outcome onto the faces.
57
58 ## Invocation order
59
60 To be able to recompute the linear combination of rhs guesses from the
61 Runge-Kutta scheme, this projection has to be invoked before we determine
62 the final linear combination per cell and overwrite the solution in the
63 patch. This means its descend_invocation_order has to be higher than
64 the one from the linear combination computation, as touchCellLastTime()
65 inverts this order.
66
67 ## Write data
68
69 The data is written onto the faces' QUpdate data field. Faces in the
70 Runge-Kutta scheme hold three different Q representations:
71
72 1. The old solution from the previous time step.
73 2. The current solution.
74 3. The update solution, i.e. the one we have just projected.
75
76 It is the your job to ensure that the updated solution is later on copied
77 over into the new solution, so it is taken into account by the halo/patch
78 reconstruction.
79
80 ## Realisation
81
82 I did not find a convenient, simple way how to tweak the projection from
83 the toolbox according to this scheme. So I ended up with copying the whole
84 thing over and change the critical things myself.
85
86 ## Boundary flags
87
88 Another difference compared to the toolbox projection is that the projection also sets the isUpdated() flag and the
89 UpdatedTimeStamp(). As the projection writes to these two updated records, it is
90 important that you roll it over afterwards. This is done via the mapping
91 RollOverUpdatedFace.
92
93 It is important to study this action set in combination with DynamicAMR. In the
94 documentation of the latter I explain why we need the guard
95
96 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
97 not marker.hasBeenRefined()
98 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
99
100 if we want to support dynamic coarsening.
101
102 ## Use action set in additional helper sweeps
103
104 The action set expects users to specify one predicate per Runge-Kutta step.
105 This means there is one if that checks which Runge-Kutta step is active.
106 This one is typically combined with two other checks: Is the solver
107 currently handling an unrefined cell (we only solve on the finest level)
108 and are we in the right phase of an enclave solve. The latter is relevant
109 if and only if we work with enclave tasking.
110
111 If the action set is used by an additional grid sweep outside of the actual
112 time stepping, the checks become simpler:
113
114 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
115 project_patch_onto_faces = exahype2.solvers.rkfd.actionsets.ProjectPatchOntoFaces(my_solver)
116 project_patch_onto_faces.guards = [ my_solver._store_cell_data_default_guard() for x in range(0,my_solver.number_of_Runge_Kutta_steps()+1) ]
117 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
118
119 Here, we basically say "for each unrefined cell always for each Runge-Kutta
120 step as well as for the final linear combination". The latter one is
121 reflected by the +1.
122
123 Besides this, you also have to roll over the data:
124
125
126 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
127 roll_over_projected_faces = exahype2.solvers.rkfd.actionsets.RollOverUpdatedFace(my_solver,
128 my_solver._store_face_data_default_guard(),
129 )
130 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
131
132 All solvers at the moment send out their face data in the suspended mode,
133 so you should be fine here.
134
135
136 """
137
138
139
140 __Template_TouchCellLastTime = """
141 logTraceIn( "touchCellLastTime(...)---ProjectPatchOntoFaces" );
142 {% for PREDICATE_NO in range(0,PREDICATES|length) %}
143 if ({{PREDICATES[PREDICATE_NO]}}) {
144
145 ::exahype2::enumerator::AoSLexicographicEnumerator enumeratorWithAuxiliaryVariables( 1, {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}}, 0, {{NUMBER_OF_UNKNOWNS}}, {{NUMBER_OF_AUXILIARY_VARIABLES}});
146 ::exahype2::enumerator::AoSLexicographicEnumerator enumeratorWithoutAuxiliaryVariables( {{RK_STEPS}}, {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}}, 0, {{NUMBER_OF_UNKNOWNS}}, 0 );
147
148 // Set the variable
149 // double timeStepSize
150 {{COMPUTE_TIME_STEP_SIZE}}
151
152 for(int d=0; d<Dimensions; d++) {
153 /**
154 * d-loop over all dimensions except d. The vector k's entry d is set
155 * to 0. We start with the left/bottom face, i.e. the one closer to the
156 * coordinate system's origin.
157 */
158 dfore(k,{{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}},d,0) {
159 for (int i=0; i<{{OVERLAP}}; i++) {
160 tarch::la::Vector<Dimensions,int> patchCell = k;
161 tarch::la::Vector<Dimensions,int> overlapCell = k;
162 patchCell(d) = i;
163 overlapCell(d) = i+{{OVERLAP}};
164
165 // I could use the enumerator as well. Doesn't really matter here. One
166 // thing is tied to ExaHyPE, the other one to the blockstructured toolbox.
167 int overlapCellSerialised = toolbox::blockstructured::serialiseVoxelIndexInOverlap(overlapCell,{{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}},{{OVERLAP}},d);
168 logDebug(
169 "touchCellLastTime(....)",
170 "project " << patchCell << "->" << overlapCell <<
171 "(" << enumeratorWithAuxiliaryVariables(0,patchCell,0) << "->" << overlapCellSerialised*({{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}) << "): " <<
172 fineGridCellEulerRKFDQ.value[ enumeratorWithAuxiliaryVariables(0,patchCell,0) ]
173 );
174 for (int j=0; j<{{NUMBER_OF_UNKNOWNS}}; j++) {
175 {{FACES_ACCESSOR}}(d).value[overlapCellSerialised*({{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}})+j] =
176 {{CELL_ACCESSOR}}.value[ enumeratorWithAuxiliaryVariables(0,patchCell,j) ];
177
178 {% if PREDICATE_NO<PREDICATES|length %}
179 {% for WEIGHT_NO in range(0,BUTCHER_TABLEAU_WEIGHTS[PREDICATE_NO+1]|length) %}
180 {% if BUTCHER_TABLEAU_WEIGHTS[PREDICATE_NO+1][WEIGHT_NO]!=0 %}
181 {{FACES_ACCESSOR}}(d).value[overlapCellSerialised*({{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}})+j] +=
182 {{BUTCHER_TABLEAU_RELATIVE_TIME_STEP_SIZES[PREDICATE_NO+1]}} * timeStepSize * {{BUTCHER_TABLEAU_WEIGHTS[PREDICATE_NO+1][WEIGHT_NO]}} *
183 fineGridCell{{UNKNOWN_IDENTIFIER}}RhsEstimates.value[ enumeratorWithoutAuxiliaryVariables({{PREDICATE_NO}},patchCell,j) ];
184 {% endif %}
185 {% endfor %}
186 {% endif %}
187 }
188
189 for (int j={{NUMBER_OF_UNKNOWNS}}; j<{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}; j++) {
190 {{FACES_ACCESSOR}}(d).value[overlapCellSerialised*({{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}})+j] =
191 {{CELL_ACCESSOR}}.value[ enumeratorWithAuxiliaryVariables(0,patchCell,j) ];
192 }
193
194 patchCell(d) = i+{{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}}-{{OVERLAP}};
195 overlapCell(d) = i;
196
197 overlapCellSerialised = toolbox::blockstructured::serialiseVoxelIndexInOverlap(overlapCell,{{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}},{{OVERLAP}},d);
198 logDebug(
199 "touchCellLastTime(....)",
200 "project " << patchCell << "->" << overlapCell <<
201 "(" << enumeratorWithAuxiliaryVariables(0,patchCell,0) << "->" << overlapCellSerialised*({{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}) << "): " <<
202 fineGridCellEulerRKFDQ.value[ enumeratorWithAuxiliaryVariables(0,patchCell,0) ]
203 );
204 for (int j=0; j<{{NUMBER_OF_UNKNOWNS}}; j++) {
205 {{FACES_ACCESSOR}}(d+Dimensions).value[overlapCellSerialised*({{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}})+j] =
206 {{CELL_ACCESSOR}}.value[ enumeratorWithAuxiliaryVariables(0,patchCell,j) ];
207
208 {% if PREDICATE_NO<PREDICATES|length %}
209 {% for WEIGHT_NO in range(0,BUTCHER_TABLEAU_WEIGHTS[PREDICATE_NO+1]|length) %}
210 {% if BUTCHER_TABLEAU_WEIGHTS[PREDICATE_NO+1][WEIGHT_NO]!=0 %}
211 {{FACES_ACCESSOR}}(d+Dimensions).value[overlapCellSerialised*({{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}})+j] +=
212 {{BUTCHER_TABLEAU_RELATIVE_TIME_STEP_SIZES[PREDICATE_NO+1]}} * timeStepSize * {{BUTCHER_TABLEAU_WEIGHTS[PREDICATE_NO+1][WEIGHT_NO]}} *
213 fineGridCell{{UNKNOWN_IDENTIFIER}}RhsEstimates.value[ enumeratorWithoutAuxiliaryVariables({{PREDICATE_NO}},patchCell,j) ];
214 {% endif %}
215 {% endfor %}
216 {% endif %}
217 }
218
219 for (int j={{NUMBER_OF_UNKNOWNS}}; j<{{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}}; j++) {
220 {{FACES_ACCESSOR}}(d+Dimensions).value[overlapCellSerialised*({{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}})+j] =
221 {{CELL_ACCESSOR}}.value[ enumeratorWithAuxiliaryVariables(0,patchCell,j) ];
222 }
223 } // overlap depth loop
224 } // dfore, i..e loop over submanifold
225 } // loop over dimensions
226 } // if predicate holds, i.e. right RK step
227 {% endfor %} // all the predicates are done now
228
229 for (int d=0; d<Dimensions; d++) {
230 {{FACE_METADATA_ACCESSOR}}(d).setUpdated(1,true);
231 {{FACE_METADATA_ACCESSOR}}(d).setUpdatedTimeStamp(1,{{CELL_METADATA_ACCESSOR}}.getTimeStamp());
232 {{FACE_METADATA_ACCESSOR}}(d+Dimensions).setUpdated(0,true);
233 {{FACE_METADATA_ACCESSOR}}(d+Dimensions).setUpdatedTimeStamp(0,{{CELL_METADATA_ACCESSOR}}.getTimeStamp());
234 }
235 logTraceOut( "touchCellLastTime(...)---ProjectPatchOntoFaces" );
236"""
237
238 def __init__(self, solver):
239 AbstractRKFDActionSet.__init__(self, solver)
240 self._guards = []
242
243 @property
244 def guards(self):
245 if self._guards == []:
246 raise Exception("Guards are not initialised")
247 return self._guards
248
249 @guards.setter
250 def guards(self, new_guards):
251 if (
252 new_guards != []
253 and len(new_guards) != self._solver.number_of_Runge_Kutta_steps() + 1
254 ):
255 raise Exception(
256 "Expect one guard per Runge Kutta step and a projection throughout the initialisation. Have {} steps but got guards {}".format(
257 self._solver.number_of_Runge_Kutta_steps(), new_guards
258 )
259 )
260 self._guards = new_guards
261
263 return __name__.replace(".py", "").replace(".", "_")
264
265 def get_body_of_operation(self, operation_name):
266 result = "\n"
267 if operation_name == ActionSet.OPERATION_TOUCH_CELL_LAST_TIME:
268 d = {}
269 if self._solver._patch_overlap_update.dim[0] % 2 != 0:
270 print(
271 "Error: Patch associated to face has to have even number of cells. Otherwise, it is not a symmetric overlap."
272 )
273 assert patch_overlap.dim[0] % 2 == 0
274 if self._solver._patch.dim[0] != self._solver._patch.dim[1]:
275 print("Error: Can only handle square patches.")
276 assert patch.dim[0] == patch.dim[1]
277 if self._solver._patch_overlap_update.dim[1] != self._solver._patch.dim[0]:
278 print("Error: Patch of overlap and patch of cell have to match")
279 assert (
280 self._solver._patch_overlap_update.dim[1]
281 == self._solver._patch.dim[0]
282 )
283
284 d["PREDICATES"] = self.guardsguards
285 d["FACE_METADATA_ACCESSOR"] = (
286 "fineGridFaces" + self._solver._face_label.name
287 )
288 d["CELL_METADATA_ACCESSOR"] = (
289 "fineGridCell" "" + self._solver._cell_label.name
290 )
291 d["FACES_ACCESSOR"] = (
292 "fineGridFaces" + self._solver._patch_overlap_update.name
293 )
294 d["CELL_ACCESSOR"] = "fineGridCell" + self._solver._patch.name
295 d["BUTCHER_TABLEAU_WEIGHTS"] = self._butcher_tableau.weight_matrix()
296 d[
297 "BUTCHER_TABLEAU_RELATIVE_TIME_STEP_SIZES"
298 ] = self._butcher_tableau.time_step_sizes()
299
300 self._solver._init_dictionary_with_default_parameters(d)
301 self._solver.add_entries_to_text_replacement_dictionary(d)
302
303 result = jinja2.Template(self.__Template_TouchCellLastTime).render(**d)
304 pass
305 return result
306
307 def get_includes(self):
308 return """
309#include "peano4/utils/Loop.h"
310#include "toolbox/blockstructured/Enumeration.h"
311#include "exahype2/fd/BoundaryConditions.h"
312#include "exahype2/enumerator/enumerator.h"
313""" + AbstractRKFDActionSet.get_includes(
314 self
315 )
Project patch data onto faces, so the faces hold valid data which can feed into the subsequent Runge-...
get_action_set_name(self)
You should replicate this function in each subclass, so you get meaningful action set names (otherwis...
get_body_of_operation(self, operation_name)
Return actual C++ code snippets to be inserted into C++ code.
Action set (reactions to observations made by the grid traversal observer)
Definition ActionSet.py:6