Peano
Loading...
Searching...
No Matches
ConstrainedPoissonEquationForMarkerOnCells.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
4
5import peano4
6import peano4.datamodel
11
13
14import jinja2
15
16from abc import abstractmethod
17
18import exahype2
19
20import dastgen2
22
23
25 """
26 Use Finite Differences Poisson equation system solver to guide cell labels
27
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
31 solve
32
33 @f$ - \Delta u = f @f$
34
35 for a fixed right-hand side f with a Jacobi solver. The solution is
36 subject to the constraints
37
38 @f$ u_{\text{min}} \leq u \leq u_{\text{max}} @f$
39
40 and Dirichlet boundary conditions.
41
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$:
46
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
50 solution exists here.
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
53 solution exists here.
54 - @f$ U_{\text{max}} < u @f$: Solver B is active. Solver A is not active.
55
56
57
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
61 change.
62
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
67 are disabled.
68
69
70 ## Usage
71
72
73 To use the marker within your project, add a new solver to your project:
74
75 marker_solver = exahype2.solvers.elliptic.ConstrainedPoissonEquationOnCells( name="ConstrainedPoissonEquationOnCells" )
76 project.add_solver( marker_solver )
77
78
79 """
80
81
82 def __init__(self,
83 name,
84 relaxation_coefficient = 0.6,
85 right_hand_side = -25.0,
86 min_value = 0.0,
87 max_value = 1.0,
88 min_threshold = 0.4,
89 max_threshold = 0.6,
90 plot = False):
91 """
92
93 name: string
94 A unique name for the solver. This one will be used for all generated
95 classes.
96
97 plot: Boolean
98 Clarifies whether a dump of the data should be enriched with grid info
99 (such as enclave status flags), too.
100
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.
106
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.
110
111 """
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"
116
117 self._name = name
118 self._relaxation_coefficient = relaxation_coefficient
119 self._right_hand_side = right_hand_side
120 self._min_value = min_value
121 self._max_value = max_value
122 self._min_threshold = min_threshold
123 self._max_threshold = max_threshold
124 self._plot = plot
125
127
129 self.create_action_sets()
130
131
132 def __str__(self):
133 return """
134Name: {}
135Type: {}
136Relaxation coefficient: {}
137""".format(
138 self._name,
139 self.__class__.__name__,
141)
142
143
144 __repr__ = __str__
145
146
147 def create_readme_descriptor(self, domain_offset, domain_size):
148 return """
149### ExaHyPE 2 solver
150
151""" + str(self)
152
153
155 """
156
157 Create cell and face data structures required. For the face data,
158 we rely on the multigrid toolbox routines.
159
160 """
161 Variants = ["AnotB","AdeterminesB","BdeterminesA","BnotA"]
162
164 self._cell_data_model.data.add_attribute(dastgen2.attributes.Boolean("Invariant") )
165 self._cell_data_model.data.add_attribute(dastgen2.attributes.Double("Value") )
166 self._cell_data_model.data.add_attribute(dastgen2.attributes.Enumeration( name="Marker", variants=Variants ) )
167
168 self._face_data_model = peano4.toolbox.multigrid.cellbased.ScalarJacobiWithRediscretisation.construct_face_helper_data(self._cell_data_model)
169
170
171 @abstractmethod
173 """
174
175 Overwrite in subclasses if you wanna create different
176 action sets.
177
178 ## Call order and ownership
179
180 This operation can be called multiple times. However, only the very
181 last call matters. All previous calls are wiped out.
182
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.
189
190 ## Postprocessing algorithm
191
192 After we have updated the cell, i.e. computed the Poisson equation update,
193 we launch a series of postprocessing steps:
194
195 - We impose the value constraints, i.e. ensure the solution is within min
196 and max.
197
198
199 """
201 cell_data = self._cell_data_model,
202 face_data = self._face_data_model,
203 relaxation_coefficient = self._relaxation_coefficient,
204 rhs_expr = self._right_hand_side,
205 unknown_setter = "setValue",
206 unknown_getter = "getValue",
207 guard = "repositories::isLastGridSweepOfTimeStep() and not fineGridCellAlphaPoissonMarker.getInvariant()"
208 )
209
210 self._action_set_smoother.additional_includes = self._get_default_includes()
211 self._action_set_smoother.d[ "SOLVER_INSTANCE" ] = self.get_name_of_global_instance()
212
213 if self._postprocess_updated_cell!=None:
214 self._action_set_smoother._Template_TouchCellFirstTime_UpdateSolution += self._postprocess_updated_cell
215
216 self._action_set_smoother._Template_TouchCellFirstTime_UpdateSolution += """
217 repositories::{{SOLVER_INSTANCE}}.computeNewMarkerState(
218 fineGridCell{{CELL_NAME}},
219 fineGridFaces{{FACE_NAME}}
220 );
221
222 repositories::{{SOLVER_INSTANCE}}.constrainValue( fineGridCell{{CELL_NAME}} );
223 """
224
225
227 return self._name +"PoissonMarker"
228
229
231 return "instanceOf" + self._name
232
233
234 def add_to_Peano4_datamodel( self, datamodel, verbose ):
235 """
236
237 Add all required data to the Peano4 project's datamodel
238 so it is properly built up
239
240 """
241 if verbose:
242 print( "Constrained Poisson Equation" )
243 print( "----------" )
244 print( str(self._name) )
245 datamodel.add_cell(self._cell_data_model)
246 datamodel.add_face(self._face_data_model)
247
248
250 """
251 Tell Peano what data to move around
252
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.
256
257 """
258 step.use_cell(self._cell_data_model)
259 step.use_face(self._face_data_model)
260
261
263 return """
264#include "tarch/la/Vector.h"
265
266#include "peano4/utils/Globals.h"
267#include "peano4/utils/Loop.h"
268
269#include "repositories/SolverRepository.h"
270
271#include "Constants.h"
272"""
273
274
276 """
277
278
279 """
280 step.add_action_set( self._action_set_smoother )
281 pass
282
283
284 def add_actions_to_create_grid(self, step, evaluate_refinement_criterion):
285 """
286
287 It is important to plug into the creation, as this is the place where
288 we get the creational events.
289
290 """
291 step.add_action_set( self._action_set_smoother )
292 pass
293
294
295 def add_actions_to_plot_solution(self, step, output_path):
296 """
297
298 We plot if and only if we are asked to do so
299
300 """
301 step.add_action_set( self._action_set_smoother )
302 if self._plot:
304 filename = output_path + "marker-" + self._name,
305 cell_unknown = self._cell_data_model,
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"
312"""
313 )
314 plotter_action_set.Template_TouchCellFirstTime = """
315 int cellIndex = _writer->plotPatch(
316 marker.x()-0.5 * marker.h(),
317 marker.h()
318 );
319
320 double data[] = {
321 fineGridCell{{CELL_UNKNOWN_NAME}}.getValue(),
322 static_cast<double>(fineGridCell{{CELL_UNKNOWN_NAME}}.getMarker()),
323 static_cast<double>(fineGridCell{{CELL_UNKNOWN_NAME}}.getInvariant()),
324 };
325
326 _dataWriter->plotCell( cellIndex, data );
327"""
328 step.add_action_set( plotter_action_set )
329
330
331
333 """
334
335 AMR
336
337 It is important that we do the inter-grid transfer operators before we
338 apply the boundary conditions.
339
340 """
341 step.add_action_set( self._action_set_smoother )
342 pass
343
344
345 @abstractmethod
347 pass
348
349
350 def add_implementation_files_to_project(self,namespace,output, subdirectory=""):
351 """
352
353 The ExaHyPE2 project will call this operation when it sets
354 up the overall environment.
355
356 This routine is typically not invoked by a user.
357
358 output: peano4.output.Output
359
360 """
361 templatefile_prefix = os.path.dirname( os.path.realpath(__file__) ) + "/ConstrainedPoissonEquationForMarkerOnCells"
362
363 if(subdirectory):
364 subdirectory += "/"
365
366 implementationDictionary = {}
367 self.add_entries_to_text_replacement_dictionary(implementationDictionary)
368
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
375
376 implementationDictionary[ "MARKER" ] = self._cell_data_model.name
377
379 templatefile_prefix + ".template.h",
380 templatefile_prefix + ".template.cpp",
381 self._name,
382 namespace,
383 subdirectory + ".",
384 implementationDictionary,
385 True,
386 True)
387
388 output.add( generated_solver_files )
389 output.makefile.add_cpp_file( subdirectory + self._name + ".cpp", generated=True )
390
391
392 @property
393 def name(self):
394 return self._name
395
396
397 @property
399 return self._postprocess_updated_cell
400
401
402 @postprocess_updated_cell.setter
403 def postprocess_updated_cell(self, kernel):
404 """
405
406 Define a postprocessing routine over the data
407
408 Typically, you do something similar to
409
410 if ( not repositories::{{SOLVER_INSTANCE}}.patchCanUseStatelessPDETerms(marker.x(), marker.h(), timeStamp, timeStepSize) ) {
411 fineGridCell{{CELL_NAME}}.setValue(1.0);
412 }
413
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.
417
418
419 ## Attributes
420
421 kernel: String
422 C++ code that holds the postprocessing kernel
423
424 """
425 self._postprocess_updated_cell = kernel
427 self.create_action_sets()
cardinality is a string (you can thus use symbols as well as long as they will be defined at compile ...
Definition Double.py:6
Wrapper around C++ enumerations which is not a datatype supported natively by MPI.
Definition Enumeration.py:6
add_to_Peano4_datamodel(self, datamodel, verbose)
Add all required data to the Peano4 project's datamodel so it is properly built up.
add_actions_to_create_grid(self, step, evaluate_refinement_criterion)
It is important to plug into the creation, as this is the place where we get the creational events.
add_implementation_files_to_project(self, namespace, output, subdirectory="")
The ExaHyPE2 project will call this operation when it sets up the overall environment.
__init__(self, name, relaxation_coefficient=0.6, right_hand_side=-25.0, min_value=0.0, max_value=1.0, min_threshold=0.4, max_threshold=0.6, plot=False)
name: string A unique name for the solver.
Default superclass for any data model in Peano which is stored within the grid.
Definition DaStGen2.py:47