Peano
Loading...
Searching...
No Matches
StaticLimiting.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
4from abc import abstractmethod
5
6import jinja2
7
8import dastgen2
9import peano4.output
10import exahype2.grid
12
13from exahype2.solvers.PDETerms import PDETerms
14from exahype2.solvers.aderdg import Polynomials
15from exahype2.solvers.aderdg.Architectures import Architectures
16
17from .actionsets import (SaveNewCellData, SpreadLimiterStatus,
18 CopyAndConvertPatch, VerifyTroubledness)
19
20class StaticLimiting(object):
21 # Initialisation stores the key information about the solver and initialises a number of variables for later usage
23 self,
24 name,
25 regular_solver,
26 limiting_solver,
27 limiting_criterion_implementation = PDETerms.User_Defined_Implementation
28 ):
29 self._name = name
30
31 # Respectively solver to be used in stable regions and solver to be used in unstable regions
32 self._regular_solver = regular_solver
33 self._limiter_solver = limiting_solver
34
35 self._halo_size = int(limiting_solver._patch_overlap_old.dim[0] / 2)
36
37 """
38 These will be used by children of this class to fill them with the declaration and definition of the
39 above mandatory and optional functions.
40 """
43
46
47 """
48 Used by child classes, will contain code to be executed at the start and end respectively of the
49 time step. Typically specifies some timestep variable and advances the overall timestamp of the
50 entire solver.
51 """
54
56
57 self._physical_admissibility_criterion = limiting_criterion_implementation
58 if ( self._physical_admissibility_criterion==PDETerms.None_Implementation
59 or self._physical_admissibility_criterion==PDETerms.Empty_Implementation):
60 raise Exception("Limiting criterion cannot be none")
61
62 # create initial structures and action_sets, these are re-generated essentially anytime something about the solver is modified
65
66 """
67 Creates all the structures which are attached on the Peano grid either in a cell or on the faces.
68 These are defined as Peano patches which have given sizes to contain persistent data as well as
69 conditions for when they should be merged, sent, received, stored, loaded, etc.
70 """
71
72 @abstractmethod
74 """
75 This will serve as a face marker to communicate with neighbouring cells whether
76 a given cell is troubled, is neighbouring a troubled cell, or is fine.
77 Cells that are troubled or neighbouring a troubled cell must rollback their solution
78 to the previous one and perform one timestep using a more robust FV scheme.
79 Therefore these and their neighbours must convert the DG representation of their
80 previous solution to an FV representation and project it to their faces so that these
81 can be exchanged with their neighbours.
82 """
83 Variants = ["REGULAR","REGULAR_TO_LIMITER","LIMITER_TO_REGULAR","TROUBLED"]
84 self._cell_label = exahype2.grid.create_cell_label(self._name)
85 self._face_label = exahype2.grid.create_face_label(self._name)
86 self._cell_label.data.add_attribute(dastgen2.attributes.Enumeration( name="Troubled_Marker", variants=Variants ) )
87 self._face_label.data.add_attribute(dastgen2.attributes.Enumeration( name="Troubled_Marker", variants=Variants ) )
88
89 self._face_label.peano4_mpi_and_storage_aspect.merge_implementation += "\n _Troubled_Marker = std::max(_Troubled_Marker, neighbour.getTroubled_Marker());"
90
91 compute_time_step_size = """
92 double timeStepSize = std::min( repositories::""" + self._regular_solver.get_name_of_global_instance() + """.getAdmissibleTimeStepSize(),
93 repositories::""" + self._limiter_solver.get_name_of_global_instance() + """.getAdmissibleTimeStepSize());
94"""
95
96 self._limiter_solver._compute_time_step_size = compute_time_step_size
97 self._regular_solver._compute_time_step_size = compute_time_step_size
98
99 """
100 Creates the action sets, these are actions tied to the grid that are executed at given points of
101 grid generation, construction or steps of the ADER-DG solver. Peano handles grid structures that
102 contain data as well as operations that are executed on that data. These are the latter.
103 The guard here are the conditions for the execution of the operation, as in they are only executed
104 if the guard evaluates to true.
105 """
106
107 @abstractmethod
110 guard = ("not marker.willBeRefined() and not marker.hasBeenRefined()"
111 + " and repositories::" + self.get_name_of_global_instance() + ".getSolverState()==" + self._name + "::SolverState::TimeStepping"
112 + " and " + "repositories::" + self.get_name_of_global_instance() + ".isFirstGridSweepOfTimeStep()"),
113 transmit_cell_timestamps=True
114 )
115
117 guard = ("not marker.willBeRefined() and not marker.hasBeenRefined() and" +
118 " repositories::" + self.get_name_of_global_instance() + ".getSolverState()==" + self._name + "::SolverState::GridInitialisation"
119 ))
120
122 guard = ("not marker.willBeRefined() and not marker.hasBeenRefined() and" +
123 " repositories::" + self.get_name_of_global_instance() + ".getSolverState()==" + self._name + "::SolverState::LimiterStatusSpreadingToNeighbours"
124 ))
125
127 regularToLimiterGuard = ("not marker.willBeRefined() and not marker.hasBeenRefined()"
128 + " and repositories::" + self.get_name_of_global_instance() + ".getSolverState()==" + self._name + "::SolverState::TimeStepping"
129 + " and fineGridCell" + self._name + "CellLabel.getTroubled_Marker()<=celldata::" + self._name + "CellLabel::Troubled_Marker::REGULAR_TO_LIMITER"
130 + " and " + "repositories::" + self.get_name_of_global_instance() + ".isFirstGridSweepOfTimeStep()"),
131 limiterToRegularGuard = ( "not marker.willBeRefined() and not marker.hasBeenRefined()"
132 + " and " + "repositories::" + self.get_name_of_global_instance() + ".isLastGridSweepOfTimeStep()"
133 + " and fineGridCell" + self._name + "CellLabel.getTroubled_Marker()==celldata::" + self._name + "CellLabel::Troubled_Marker::LIMITER_TO_REGULAR")
134 )
135
136 """
137 Saving cell needs to happen before the cell is used for any computations
138 checking troubledness needs to happen in the initialization after the regular solver
139 has computed its initial condition.
140 spreading limiter status happens once in its own dedicated step after initialization
141 but before either solvers are allowed to start computing, and once in the first
142 timestep before any other operations are performed
143 copying and converting the data from regular to limiter and from limiter to regular
144 happens before they are used in the first grid traversal, so either at the very
145 beginning of the first grid traversal or at the end of the second, so that
146 regular-to-limiter and limiter-to-regular cells both have access to either version
147 and can therefore project both variants of their faces
148 """
149 self._action_set_save_new_cell_data.descend_invocation_order = 0
150 self._action_set_check_troubledness.descend_invocation_order = 1000
151 self._action_set_spread_limiter_status.descend_invocation_order = 0
152 self._action_set_copy_convert_and_project_data.descend_invocation_order = 0
153
154 self._regular_solver._action_set_prediction.guard += (
155 " and fineGridCell" + self._name + "CellLabel.getTroubled_Marker()<=celldata::" + self._name + "CellLabel::Troubled_Marker::LIMITER_TO_REGULAR"
156 )
157
158 self._regular_solver._action_set_correction.guard += (
159 " and fineGridCell" + self._name + "CellLabel.getTroubled_Marker()<=celldata::" + self._name + "CellLabel::Troubled_Marker::REGULAR_TO_LIMITER"
160 )
161
162 self._limiter_solver._action_set_update_cell.guard += (
163 " and " + "repositories::" + self.get_name_of_global_instance() + ".isLastGridSweepOfTimeStep()"
164 + " and fineGridCell" + self._name + "CellLabel.getTroubled_Marker()>=celldata::" + self._name + "CellLabel::Troubled_Marker::LIMITER_TO_REGULAR"
165 )
166
167 def create_readme_descriptor(self, domain_offset, domain_size):
168 return (
169 """ ExaHyPE 2 Limiting solver"""
170 )
171
172 def get_user_action_set_includes(self):
173 return self.user_action_set_includes
174
175 def get_user_solver_includes(self):
176 return self.user_solver_includes
177
178 def add_user_action_set_includes(self, value):
179 """
180 Add further includes to this property, if your action sets require some additional
181 routines from other header files.
182 """
183 self.user_action_set_includes += value
184
185 def _store_cell_data_default_guard(self):
186 return (
187 "not marker.willBeRefined() "
188 + "and repositories::"
189 + self.get_name_of_global_instance()
190 + ".getSolverState()!="
191 + self._name
192 + "::SolverState::GridConstruction"
193 )
194
195 def _provide_cell_data_to_compute_kernels_default_guard(self):
196 return (
197 "not marker.willBeRefined() "
198 + "and repositories::"
199 + self.get_name_of_global_instance()
200 + ".getSolverState()!="
201 + self._name
202 + "::SolverState::GridConstruction"
203 )
204
205 def _provide_face_data_to_compute_kernels_default_guard(self):
206 return (
207 "not marker.willBeRefined() "
208 + "and repositories::"
209 + self.get_name_of_global_instance()
210 + ".getSolverState()!="
211 + self._name
212 + "::SolverState::GridConstruction"
213 )
214
215 def _load_cell_data_default_guard(self):
216 return (
217 "not marker.hasBeenRefined() "
218 + "and repositories::"
219 + self.get_name_of_global_instance()
220 + ".getSolverState()!="
221 + self._name
222 + "::SolverState::GridConstruction "
223 + "and repositories::"
224 + self.get_name_of_global_instance()
225 + ".getSolverState()!="
226 + self._name
227 + "::SolverState::GridInitialisation"
228 )
229
230 def _store_face_data_default_guard(self):
231 return (
232 "not marker.willBeRefined() "
233 + "and repositories::"
234 + self.get_name_of_global_instance()
235 + ".getSolverState()!="
236 + self._name
237 + "::SolverState::GridConstruction"
238 )
239
240 def _load_face_data_default_guard(self):
241 return (
242 "not marker.hasBeenRefined() "
243 + "and repositories::"
244 + self.get_name_of_global_instance()
245 + ".getSolverState()!="
246 + self._name
247 + "::SolverState::GridConstruction "
248 + "and repositories::"
249 + self.get_name_of_global_instance()
250 + ".getSolverState()!="
251 + self._name
252 + "::SolverState::GridInitialisation"
253 )
254
255 def _store_boundary_data_default_guard(self):
256 return (
257 self._store_face_data_default_guard()
258 + " and not repositories::"
259 + self.get_name_of_global_instance()
260 + ".PeriodicBC[marker.getSelectedFaceNumber()%Dimensions]"
261 + " and not marker.hasBeenRefined() and fineGridFace"
262 + self._name
263 + "FaceLabel.getBoundary()"
264 )
265
266 def _clear_face_data_default_guard(self):
267 return (
268 "not marker.willBeRefined()"
269 + "and repositories::"
270 + self.get_name_of_global_instance()
271 + ".getSolverState()!="
272 + self._name
273 + "::SolverState::GridConstruction "
274 )
275
276 def _interpolate_face_data_default_guard(self):
277 return (
278 "repositories::"
279 + self.get_name_of_global_instance()
280 + ".getSolverState()!="
281 + self._name
282 + "::SolverState::GridInitialisation"
283 )
284
285 def _restrict_face_data_default_guard(self):
286 return "true"
287
288 def _unknown_identifier(self):
289 return self._name + "Q"
290
291 def get_name_of_global_instance(self):
292 return "instanceOf" + self._name
293
294 def _get_default_includes(self):
295 return """
296#include "tarch/la/Vector.h"
297
298#include "peano4/utils/Globals.h"
299#include "peano4/utils/Loop.h"
300
301#include "repositories/SolverRepository.h"
302"""
303
304 """
305 The following functions will be called by the Peano4 project,
306 they tell Peano which of our generated data it should attach
307 to each grid element and use.
308 """
309
310 def __str__(self):
311 result = (
312 """
313Name: """ + self._name + """
314Type: """ + self.__class__.__name__
315 )
316 return result
317
318 def add_to_Peano4_datamodel(self, datamodel, verbose):
319 datamodel.add_cell(self._cell_label)
320 datamodel.add_face(self._face_label)
321
322 def add_use_data_statements_to_Peano4_solver_step(self, step):
323 step.use_cell(self._cell_label)
324 step.use_face(self._face_label)
325
326 def add_actions_to_init_grid(self, step):
327 step.add_action_set(self._action_set_check_troubledness)
328 pass
329
330 def add_actions_to_create_grid(self, step, evaluate_refinement_criterion):
331 pass
332
333 def add_actions_to_plot_solution(self, step, output_path):
334 pass
335
336 # These actions are executed during each time step
337 def add_actions_to_perform_time_step(self, step):
338 step.add_action_set(self._action_set_save_new_cell_data)
339 step.add_action_set(self._action_set_spread_limiter_status)
340 step.add_action_set(self._action_set_copy_convert_and_project_data)
341
342 def set_plot_description(self, description):
343 self.plot_description = description
344
345 def _generate_kernels(self, namespace, output, subdirectory, dimensions):
346 d = {}
347 self._init_dictionary_with_default_parameters(d)
348 self.add_entries_to_text_replacement_dictionary(d)
349 d["DIMENSIONS"] = dimensions
350 d["NUMBER_OF_DOF_3D"] = 1 if dimensions == 2 else d["NUMBER_OF_DOF"]
351 d["NUMBER_OF_DOF_LIMITER_3D"] = 1 if dimensions == 2 else d["NUMBER_OF_DOF_LIMITER"]
352 d["GHOST_LAYER_WIDTH_3D"] = 0 if dimensions == 2 else d["GHOST_LAYER_WIDTH"]
353 d["FULL_QUALIFIED_SOLVER_NAME"] ="::".join(namespace) + "::" + self._name
354
355 kernels_namespace = namespace + ["kernels", self._name]
356 kernels_output_path = subdirectory + "kernels/" + self._name
357 kernels_templates_prefix = os.path.dirname(os.path.realpath(__file__)) + "/kernels/"
358
359 if not os.path.exists(kernels_output_path):
360 os.makedirs(kernels_output_path)
361
362 exahype2.solvers.limiting.kernels.Quadrature(d).generate_kernels(
363 kernels_namespace, output, subdirectory, kernels_templates_prefix, kernels_output_path
364 )
365
366 exahype2.solvers.limiting.kernels.Limiter(d).generate_kernels(
367 kernels_namespace, output, subdirectory, kernels_templates_prefix, kernels_output_path
368 )
369
370 """
371 This tells peano to add the solver files to the project. It generates any files that are not
372 specifically cell or face data or grid actions such as the actionsets.
373 Currently it generates solver implementation and header files.
374 """
375 def add_implementation_files_to_project(self, namespace, output, dimensions, subdirectory=""):
376 """
377 The ExaHyPE2 project will call this operation when it sets
378 up the overall environment.
379
380 This routine is typically not invoked by a user.
381
382 output: peano4.output.Output
383 """
384
385 if(subdirectory):
386 subdirectory += "/"
387
388 HeaderDictionary = {}
389 self._init_dictionary_with_default_parameters(HeaderDictionary)
390
391 generated_abstract_files = (
392 peano4.output.Jinja2TemplatedHeaderImplementationFilePair(
393 os.path.dirname(os.path.realpath(__file__)) + "/StaticLimitingAbstract.template.h",
394 os.path.dirname(os.path.realpath(__file__)) + "/StaticLimitingAbstract.template.cpp",
395 "Abstract" + self._name,
396 namespace,
397 subdirectory + ".",
398 HeaderDictionary,
399 True,
400 True
401 )
402 )
403
404 generated_solver_files = (
405 peano4.output.Jinja2TemplatedHeaderImplementationFilePair(
406 os.path.dirname(os.path.realpath(__file__)) + "/StaticLimiting.template.h",
407 os.path.dirname(os.path.realpath(__file__)) + "/StaticLimiting.template.cpp",
408 self._name,
409 namespace,
410 subdirectory + ".",
411 HeaderDictionary,
412 False,
413 True
414 )
415 )
416
417 output.add(generated_abstract_files)
418 output.makefile.add_cpp_file(subdirectory + "Abstract" + self._name + ".cpp", generated=True)
419 output.makefile.add_h_file(subdirectory + "Abstract" + self._name + ".h", generated=True)
420
421 output.add(generated_solver_files)
422 output.makefile.add_cpp_file(subdirectory + self._name + ".cpp", generated=True)
423 output.makefile.add_h_file(subdirectory + self._name + ".h", generated=True)
424
425 self._generate_kernels(namespace, output, subdirectory, dimensions)
426
427 @abstractmethod
428 def add_entries_to_text_replacement_dictionary(self, d):
429 pass
430
431 """
432 Generates a dictionary of "words" that will later be used in various templates to fill these out by adding
433 information from the solver or which was specified by the user.
434 """
435 def _init_dictionary_with_default_parameters(self, d):
436 d["SOLVER_INSTANCE"] = self.get_name_of_global_instance()
437 d["REGULAR_SOLVER_INSTANCE"] = self._regular_solver.get_name_of_global_instance()
438 d["LIMITER_SOLVER_INSTANCE"] = self._limiter_solver.get_name_of_global_instance()
439
440 d["SOLVER_NAME"] = self._name
441 d["REGULAR_SOLVER_NAME"] = self._regular_solver._name
442 d["LIMITER_SOLVER_NAME"] = self._limiter_solver._name
443
444 d["UNKNOWN_IDENTIFIER"] = self._unknown_identifier()
445 d["REGULAR_SOLVER_UNKNOWN_IDENTIFIER"] = self._regular_solver._unknown_identifier()
446 d["LIMITER_SOLVER_UNKNOWN_IDENTIFIER"] = self._limiter_solver._unknown_identifier()
447
448 d["ALIGNMENT"] = Architectures.get_alignment(self._regular_solver._architecture)
449 d["SIMD_WIDTH"] = Architectures.get_simd_width(self._regular_solver._architecture)
450 d["ORDER"] = self._regular_solver._order
451 d["NUMBER_OF_DOF"] = self._regular_solver._order + 1
452 d["NUMBER_OF_DOF_PADDED"] = Architectures.get_size_with_padding(self._regular_solver._architecture, self._regular_solver._order + 1)
453 d["NUMBER_OF_DOF_LIMITER"] = self._limiter_solver.patch_size if self._limiter_solver.patch_size >= 0 else 2 * self._regular_solver._order + 1
454 d["NUMBER_OF_DOF_LIMITER_PADDED"] = Architectures.get_size_with_padding(self._regular_solver._architecture, d["NUMBER_OF_DOF_LIMITER"])
455 d["NUMBER_OF_UNKNOWNS"] = self._regular_solver._unknowns
456 d["NUMBER_OF_UNKNOWNS_PADDED"] = Architectures.get_size_with_padding(self._regular_solver._architecture, self._regular_solver._unknowns)
457 d["NUMBER_OF_AUXILIARY_VARIABLES"] = self._regular_solver._auxiliary_variables
458 d["NUMBER_OF_AUXILIARY_VARIABLES_PADDED"] = Architectures.get_size_with_padding(self._regular_solver._architecture, self._regular_solver._auxiliary_variables)
459 d["NUMBER_OF_DATA"] = self._regular_solver._unknowns + self._regular_solver._auxiliary_variables
460 d["NUMBER_OF_DATA_PADDED"] = Architectures.get_size_with_padding(self._regular_solver._architecture, self._regular_solver._unknowns + self._regular_solver._auxiliary_variables)
461 d["NUMBER_OF_OBSERVABLES"] = 0
462 d["GHOST_LAYER_WIDTH"] = self._halo_size
463
464 d["NUMBER_OF_DOFS_PER_CELL_2D"] = (
465 (self._regular_solver.unknowns + self._regular_solver.auxiliary_variables)
466 * (self._regular_solver.order+1) * (self._regular_solver.order+1) )
467 d["NUMBER_OF_DOFS_PER_CELL_3D"] = (
468 (self._regular_solver.unknowns + self._regular_solver.auxiliary_variables)
469 * (self._regular_solver.order+1) * (self._regular_solver.order+1) * (self._regular_solver.order+1) )
470
471 d["USE_GAUSS_LOBATTO"] = ( "true" if self._regular_solver._polynomials is Polynomials.Gauss_Lobatto else "false" )
472 d["POLYNOMIAL_TYPE"] = ( "lobatto" if self._regular_solver._polynomials is Polynomials.Gauss_Lobatto else "legendre" )
473
474 d["SOLVER_INCLUDES"] = self.get_user_solver_includes()
475
476 d["SOLVER_USER_DECLARATIONS"] = jinja2.Template(
477 self._solver_user_declarations, undefined=jinja2.DebugUndefined
478 ).render(**d)
479 d["SOLVER_USER_DEFINITIONS"] = jinja2.Template(
480 self._solver_user_definitions, undefined=jinja2.DebugUndefined
481 ).render(**d)
482
483 d["START_TIME_STEP_IMPLEMENTATION"] = jinja2.Template(
484 self._start_time_step_implementation, undefined=jinja2.DebugUndefined
485 ).render(**d)
486 d["FINISH_TIME_STEP_IMPLEMENTATION"] = jinja2.Template(
487 self._finish_time_step_implementation, undefined=jinja2.DebugUndefined
488 ).render(**d)
489
490 d["CONSTRUCTOR_IMPLEMENTATION"] = jinja2.Template(
491 self._constructor_implementation, undefined=jinja2.DebugUndefined
492 ).render(**d)
493
494 d["ADMISSIBILITY_IMPLEMENTATION"] = self._physical_admissibility_criterion
495
496 d["QUADRATURE_POINTS"] = ", ".join(self._regular_solver._basis.quadrature_points())
497
498 d["USE_LIBXSMM"] = self._regular_solver._use_libxsmm
499 d["ARCHITECTURE"] = self._regular_solver._architecture
Wrapper around C++ enumerations which is not a datatype supported natively by MPI.
Definition Enumeration.py:6
create_data_structures(self)
This will serve as a face marker to communicate with neighbouring cells whether a given cell is troub...
create_readme_descriptor(self, domain_offset, domain_size)
__init__(self, name, regular_solver, limiting_solver, limiting_criterion_implementation=PDETerms.User_Defined_Implementation)