16from abc
import abstractmethod
26 Use Finite Differences Poisson equation system solver to guide cell labels
28 The most popular application of this solver is to guide limiters:
29 Some codes need to know where to solve their PDE with an alternative
30 more robust (though likely less effective) solver. For this, they
33 @f$ - \Delta u = f @f$
35 for a fixed right-hand side f with a Jacobi solver. The solution is
36 subject to the constraints
38 @f$ u_{\text{min}} \leq u \leq u_{\text{max}} @f$
40 and Dirichlet boundary conditions.
42 Typically, the solvers than set the value of $u=1$ in regions where
43 they need the limiter. The $u$ value will fade away from these regions,
44 and we can implement the following heuristics with
45 @f$ u_{\text{min}} < U_{\text{min}} < U_{\text{max}} < u_{\text{max}}@f$:
47 - @f$ u<U_{\text{min}} @f$: Solver A is active. Solver B is not active.
48 - @f$ U_{\text{min}} \leq u < \frac{1}{2} ( U_{\text{min}} + U_{\text{max}} )@f$: Solver
49 A is active and sets the value of solver B which is not active but its
51 - @f$ \frac{1}{2} ( U_{\text{min}} + U_{\text{max}} ) \leq u \leq U_{\text{max}} @f$: Solver
52 B is active and sets the value of solver A which is not active but its
54 - @f$ U_{\text{max}} < u @f$: Solver B is active. Solver A is not active.
58 As we use a Jacobi iterative solver, using this marker within ExaHyPE is
59 equivalent to a heat equation solver which becomes stationary eventually if
60 the areas where @f$ u=u_{\text{max}} @f$ is manually set/enforced do not
63 The solver works with a constant right-hand side and is defined over two
64 attributes per cell: There's the actual value (u in the example above)
65 and a marker which switches the Jacobi smoothing on and off per cell.
66 In cells where we fix the value, we toggle the marker so Jacobi updates
73 To use the marker within your project, add a new solver to your project:
75 marker_solver = exahype2.solvers.elliptic.ConstrainedPoissonEquationOnCells( name="ConstrainedPoissonEquationOnCells" )
76 project.add_solver( marker_solver )
84 relaxation_coefficient = 0.6,
85 right_hand_side = -25.0,
94 A unique name for the solver. This one will be used for all generated
98 Clarifies whether a dump of the data should be enriched with grid info
99 (such as enclave status flags), too.
101 relaxation_coefficient: Double (0< value <1)
102 How quickly should the solution be relaxed against real solution. This
103 is the damping parameter in the Jacobi update scheme and thus has to be
104 bigger than 0 and smaller than 1. If the coefficient approaches 1, you
105 might see oscillations.
107 right_hand_side: Double (negative)
108 Basically describes how fast the solution is pulled down to 0 once we
109 go away from the AMR transition area.
112 assert min_value<max_value,
"arguments invalid"
113 assert max_threshold<=max_value,
"arguments invalid"
114 assert min_threshold>=min_value,
"arguments invalid"
115 assert min_threshold<max_threshold,
"arguments invalid"
136Relaxation coefficient: {}
139 self.__class__.__name__,
157 Create cell and face data structures required. For the face data,
158 we rely on the multigrid toolbox routines.
161 Variants = [
"AnotB",
"AdeterminesB",
"BdeterminesA",
"BnotA"]
175 Overwrite in subclasses if you wanna create different
178 ## Call order and ownership
180 This operation can be called multiple times. However, only the very
181 last call matters. All previous calls are wiped out.
183 If you have a hierarchy of solvers, every create_data_structure()
184 should first(!) call its parent version. This way, you always ensure
185 that all data are in place before you continue to alter the more
186 specialised versions. So it is (logically) a top-down (general to
187 specialised) run through all create_data_structure() variants
188 within the inheritance tree.
190 ## Postprocessing algorithm
192 After we have updated the cell, i.e. computed the Poisson equation update,
193 we launch a series of postprocessing steps:
195 - We impose the value constraints, i.e. ensure the solution is within min
205 unknown_setter =
"setValue",
206 unknown_getter =
"getValue",
207 guard =
"repositories::isLastGridSweepOfTimeStep() and not fineGridCellAlphaPoissonMarker.getInvariant()"
217 repositories::{{SOLVER_INSTANCE}}.computeNewMarkerState(
218 fineGridCell{{CELL_NAME}},
219 fineGridFaces{{FACE_NAME}}
222 repositories::{{SOLVER_INSTANCE}}.constrainValue( fineGridCell{{CELL_NAME}} );
227 return self.
_name +
"PoissonMarker"
231 return "instanceOf" + self.
_name
237 Add all required data to the Peano4 project's datamodel
238 so it is properly built up
242 print(
"Constrained Poisson Equation" )
243 print(
"----------" )
244 print( str(self.
_name) )
251 Tell Peano what data to move around
253 Inform Peano4 step which data are to be moved around via the
254 use_cell and use_face commands. This operation is generic from
255 ExaHyPE's point of view, i.e. I use it for all grid sweep types.
264#include "tarch/la/Vector.h"
266#include "peano4/utils/Globals.h"
267#include "peano4/utils/Loop.h"
269#include "repositories/SolverRepository.h"
271#include "Constants.h"
287 It is important to plug into the creation, as this is the place where
288 we get the creational events.
298 We plot if and only if we are asked to do so
304 filename = output_path +
"marker-" + self.
_name,
306 getter =
"defined below",
307 number_of_unknows_per_cell = 3,
308 description =
"value,marker(AnotB,AdeterminesB,BdeterminesA,BnotA),invariant",
309 time_stamp_evaluation=
"0.5*(repositories::getMinTimeStamp()+repositories::getMaxTimeStamp())",
310 additional_includes=
"""
311#include "repositories/SolverRepository.h"
314 plotter_action_set.Template_TouchCellFirstTime =
"""
315 int cellIndex = _writer->plotPatch(
316 marker.x()-0.5 * marker.h(),
321 fineGridCell{{CELL_UNKNOWN_NAME}}.getValue(),
322 static_cast<double>(fineGridCell{{CELL_UNKNOWN_NAME}}.getMarker()),
323 static_cast<double>(fineGridCell{{CELL_UNKNOWN_NAME}}.getInvariant()),
326 _dataWriter->plotCell( cellIndex, data );
328 step.add_action_set( plotter_action_set )
337 It is important that we do the inter-grid transfer operators before we
338 apply the boundary conditions.
353 The ExaHyPE2 project will call this operation when it sets
354 up the overall environment.
356 This routine is typically not invoked by a user.
358 output: peano4.output.Output
361 templatefile_prefix = os.path.dirname( os.path.realpath(__file__) ) +
"/ConstrainedPoissonEquationForMarkerOnCells"
366 implementationDictionary = {}
367 self.add_entries_to_text_replacement_dictionary(implementationDictionary)
369 implementationDictionary[
"RELAXATION_COEFFICIENT" ] = self._relaxation_coefficient
370 implementationDictionary[
"RHS" ] = self._right_hand_side
371 implementationDictionary[
"MIN_VALUE" ] = self._min_value
372 implementationDictionary[
"MAX_VALUE" ] = self._max_value
373 implementationDictionary[
"MIN_THRESHOLD" ] = self._min_threshold
374 implementationDictionary[
"MAX_THRESHOLD" ] = self._max_threshold
376 implementationDictionary[
"MARKER" ] = self._cell_data_model.name
379 templatefile_prefix +
".template.h",
380 templatefile_prefix +
".template.cpp",
384 implementationDictionary,
388 output.add( generated_solver_files )
389 output.makefile.add_cpp_file( subdirectory + self._name +
".cpp", generated=
True )
402 @postprocess_updated_cell.setter
406 Define a postprocessing routine over the data
408 Typically, you do something similar to
410 if ( not repositories::{{SOLVER_INSTANCE}}.patchCanUseStatelessPDETerms(marker.x(), marker.h(), timeStamp, timeStepSize) ) {
411 fineGridCell{{CELL_NAME}}.setValue(1.0);
414 We evaluate some criterion (here, we only look at a global function, but we
415 also have access to all the other solver data in the kernel) and then set
416 the cell value to the maximum if the criterion holds.
422 C++ code that holds the postprocessing kernel
cardinality is a string (you can thus use symbols as well as long as they will be defined at compile ...
Wrapper around C++ enumerations which is not a datatype supported natively by MPI.
Default superclass for any data model in Peano which is stored within the grid.