Peano
Loading...
Searching...
No Matches
CellCenteredFiniteDifferences.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
12import jinja2
13
14from abc import abstractmethod
15from enum import Enum
16
17import exahype2
19
20from exahype2.solvers.PDETerms import PDETerms
21
22from exahype2.solvers.ButcherTableau import RungeKutta_steps
23from exahype2.solvers.Storage import Storage
24
25
26
28 """
29 Abstract solver for patch-based finite diffences
30
31 All solvers in ExaHyPE are cell-centered discretisations.
32
33
34 ## Adaptive mesh refinement
35
36 We have, at the moment, no hard-coded AMR operator set available unless
37 you work with an overlap of one. In this case, you find operators in the
38 toolbox. Very few pre-manufactured operators there are ready to go for
39 higher overlaps (the injection operator is an example). In general, you
40 will have to inject your own transfer operators.
41
42 It depends on the flavour that you want to use for your interpolation and
43 restriction. A simple interpolation for an overlap of three and a patch_size
44 of five would be
45
46 self._interpolation = "tensor_product< " + self._name + ">"
47 self.add_solver_constants( "static constexpr double NormalInterpolationMatrix1d[] = {0.0, 1.0, 0.0};" )
48 self.add_solver_constants( " ""static constexpr double TangentialInterpolationMatrix1d[] = {
49 1.0,0.0,0.0,0.0,0.0,
50 1.0,0.0,0.0,0.0,0.0,
51 1.0,0.0,0.0,0.0,0.0,
52 0.0,1.0,0.0,0.0,0.0,
53 0.0,1.0,0.0,0.0,0.0,
54
55 0.0,1.0,0.0,0.0,0.0,
56 0.0,0.0,1.0,0.0,0.0,
57 0.0,0.0,1.0,0.0,0.0,
58 0.0,0.0,1.0,0.0,0.0,
59 0.0,0.0,0.0,1.0,0.0,
60
61 0.0,0.0,0.0,1.0,0.0,
62 0.0,0.0,0.0,1.0,0.0,
63 0.0,0.0,0.0,0.0,1.0,
64 0.0,0.0,0.0,0.0,1.0,
65 0.0,0.0,0.0,0.0,1.0
66 };" "" )
67 """
69 self,
70 name,
71 patch_size,
72 overlap,
73 rk_order,
74 unknowns,
75 auxiliary_variables,
76 min_meshcell_h,
77 max_meshcell_h,
78 plot_grid_properties,
79 kernel_namespace,
80 baseline_action_set_descend_invocation_order=0,
81 ):
82 """
83 name: string
84 A unique name for the solver. This one will be used for all generated
85 classes. Also the C++ object instance later on will incorporate this
86 name.
87
88 patch_size: int
89 Size of the patch in one dimension. All stuff here's dimension-generic.
90
91 overlap: int
92 That's the size of the halo layer which is half of the overlap with a
93 neighbour. A value of 1 means that a patch_size x patch_size patch in
94 2d is surrounded by one additional cell layer. The overlap has to be
95 bigger or equal to one. It has to be smaller or equal to patch_size.
96
97 unknowns: int
98 Number of unknowns per Finite Volume voxel.
99
100 auxiliary_variables: int
101 Number of auxiliary variables per Finite Volume voxel. Eventually, both
102 unknowns and auxiliary_variables are merged into one big vector if we
103 work with AoS. But the solver has to be able to distinguish them, as
104 only the unknowns are subject to a hyperbolic formulation.
105
106 min_meshcell_h: double
107 This size refers to the individual Finite Volume.
108
109 max_meshcell_h: double
110 This size refers to the individual Finite Volume.
111
112 plot_grid_properties: Boolean
113 Clarifies whether a dump of the data should be enriched with grid info
114 (such as enclave status flags), too.
115 """
116 self._name = name
117
118 self._min_meshcell_h = min_meshcell_h
119 self._max_meshcell_h = max_meshcell_h
120 self._plot_grid_properties = plot_grid_properties
121
122 self._patch_size = patch_size
123 self._overlap = overlap
124 self._rk_order = rk_order
125
126 """
127 self._unknowns and self._auxiliary_variables respectively hold the number of unknowns and
128 auxiliary variables in the equation to be computed. Unknowns are variables that change over
129 time whereas auxiliary variables can be space-dependent but don't vary over time.
130 These can be specified either as simple ints or by a dictionary
131 (e.g.) unknowns = {'a': 1, 'b': 1, 'c': 3}
132 in which the user specifies the multiplicity of the variable (the velocity has one component
133 per dimension for example.)
134 If they are specified by a dictionary then the code generates a "VariableShortcuts" file which
135 allows the user to specifiy a variable by name and automatically maps this to the right position
136 in an array for better legibility. Otherwise they must manually remember the position of each
137 variable manually.
138
139 use_var_shortcut is used to know whether or not the user passed their variables via a dict
140 variable_names and variable_pos are used internally to remember the names and respective
141 positions of variables if set by a dictionary.
142 """
143 self._use_var_shortcut = False
145 self._variable_pos = [0]
146
147 if type(unknowns) is dict:
148 self._unknowns = sum(unknowns.values())
149 self._variable_names += list(unknowns.keys())
150 for var in list(unknowns.values()):
151 self._variable_pos.append(var + self._variable_pos[-1])
152 self._use_var_shortcut = True
153 elif type(unknowns) is int:
154 self._unknowns = unknowns
155 else:
156 raise Exception(
157 "not a valid type for parameter unknowns, needs to be int or dictionary"
158 )
159
160 if type(auxiliary_variables) is dict:
161 self._auxiliary_variables = sum(auxiliary_variables.values())
162 self._variable_names += list(auxiliary_variables.keys())
163 for var in list(auxiliary_variables.values()):
164 self._variable_pos.append(var + self._variable_pos[-1])
165 self._use_var_shortcut = True
166 elif type(auxiliary_variables) is int:
167 self._auxiliary_variables = auxiliary_variables
168 else:
169 raise Exception(
170 "not a valid type for parameter auxiliary_variables, needs to be int or dictionary"
171 )
172
176
177 self._kernel_namespace = kernel_namespace
178
179 if min_meshcell_h > max_meshcell_h:
180 print(
181 "Error: min_meshcell_h ("
182 + str(min_meshcell_h)
183 + ") is bigger than max_meshcell_h ("
184 + str(max_meshcell_h)
185 + ")"
186 )
187
189 peano4.toolbox.blockstructured.ReconstructedArrayMemoryLocation.CallStack
190 )
191
193
195
197
200 self._action_set_AMR = None
207 None
208 )
211
212 self._compute_time_step_size = "#error Not yet defined"
213 self._compute_new_time_step_size = "#error Not yet defined"
216
217 self._compute_kernel_call = "#error Not yet defined"
218
224
228
232
233 self._interpolation = "#error Not yet defined"
234 self._restriction = "#error Not yet defined"
235
236 self._baseline_action_set_descend_invocation_order = baseline_action_set_descend_invocation_order
237
238 self.switch_storage_scheme(Storage.SmartPointers, Storage.SmartPointers)
239
240
241 def __str__(self):
242 result = (
243 """
244Name: """
245 + self._name
246 + """
247Type: """
248 + self.__class__.__name__
249 + """
250Patch size: """
251 + str(self._patch_size)
252 + """
253Runge-Kutta order: """
254 + str(self._rk_order)
255 + """
256Overlap: """
257 + str(self._overlap)
258 + """
259Unknowns: """
260 + str(self._unknowns)
261 + """
262Auxiliary variables: """
263 + str(self._auxiliary_variables)
264 + """
265h_meshcell_min: """
266 + str(self._min_meshcell_h)
267 + """
268h_meshcell_max: """
269 + str(self._max_meshcell_h)
270 + """
271"""
272 )
273 return result
274
275 __repr__ = __str__
276
278 coarsest_tree_level = 0
279 while (
280 domain_size * 3 ** (-coarsest_tree_level) / self._patch_size
281 > self._max_meshcell_h
282 ):
283 coarsest_tree_level += 1
284 return coarsest_tree_level
285
287 finest_tree_level = 0
288 while (
289 domain_size * 3 ** (-finest_tree_level) / self._patch_size
290 > self._min_meshcell_h
291 ):
292 finest_tree_level += 1
293 return finest_tree_level
294
295 # self._overlap = overlap
296 # self._unknowns = unknowns
297 # self._auxiliary_variables = auxiliary_variables
298
299 def get_coarsest_number_of_patches(self, domain_size):
300 return 3 ** self.get_min_number_of_spacetree_levels(domain_size)
301
302 def get_finest_number_of_patches(self, domain_size):
303 return 3 ** self.get_max_number_of_spacetree_levels(domain_size)
304
306 return self.get_coarsest_number_of_patches(domain_size) * self._patch_size
307
309 return self.get_finest_number_of_patches(domain_size) * self._patch_size
310
312 return domain_size / self.get_coarsest_number_of_compute_grid_cells(domain_size)
313
314 def get_finest_compute_grid_cell_size(self, domain_size):
315 return domain_size / self.get_finest_number_of_compute_grid_cells(domain_size)
316
317 def create_readme_descriptor(self, domain_offset, domain_size):
318 return (
319 """
320### ExaHyPE 2 solver
321
322"""
323 + str(self)
324 + """
325
326Real type of this solver: """
327 + str(type(self))
328 + """
329
330We assume that you use a domain size of (0,"""
331 + str(domain_size)
332 + """)^d. Peano 4 will cut this domain equidistantly
333and recursively into three parts along each coordinate axis. This yields a spacetree.
334
335The spacetree will at least have """
336 + str(self.get_min_number_of_spacetree_levels(domain_size))
337 + """ levels.
338The spacetree will at most have """
339 + str(self.get_max_number_of_spacetree_levels(domain_size))
340 + """ levels.
341
342The spacetree will thus span at least """
343 + str(self.get_coarsest_number_of_patches(domain_size))
344 + """ octants/cells (see clarification below) per coordinate axis. This means a """
345 + str(self.get_coarsest_number_of_patches(domain_size))
346 + """^d grid of octants.
347The spacetree will thus span at most """
348 + str(self.get_finest_number_of_patches(domain_size))
349 + """ octants/cells (see clarification below) per coordinate axis. This means a """
350 + str(self.get_finest_number_of_patches(domain_size))
351 + """^d grid of octants.
352
353ExaHyPE 2 embeds """
354 + str(self._patch_size)
355 + """^d patches of compute cells into the finest tree level.
356In the text above, we refer to the elements of this level of the tree as octants.
357The octants are squares/cubes and many papers refer to them as cells, but they are not
358the actual compute data structure. The compute data structure is the cells that
359are embedded into these finest level spacetree cells. We therefore prefer the
360term octant for the latter, whereas we use the term (compute) cell for the
361entities that are actually used for the computations, i.e. hold degrees of
362freedom, and are actually visible within Paraview, e.g.
363
364The coarsest possible mesh will consist of """
365 + str(self.get_coarsest_number_of_compute_grid_cells(domain_size))
366 + """ compute cells per coordinate axis.
367The finest possible mesh will consist of """
368 + str(self.get_finest_number_of_compute_grid_cells(domain_size))
369 + """ compute cells per coordinate axis.
370
371The coarsest mesh width of """
372 + str(self.get_coarsest_compute_grid_cell_size(domain_size))
373 + """ is thus just smaller than the maximum mesh size """
374 + str(self._max_meshcell_h)
375 + """.
376The finest mesh width of """
377 + str(self.get_finest_compute_grid_cell_size(domain_size))
378 + """ is thus just smaller than the minimum mesh size """
379 + str(self._min_meshcell_h)
380 + """.
381
382"""
383 )
384
385 @property
387 """
388
389 Add further includes to this property, if your action sets require some additional
390 routines from other header files.
391
392 """
393 return self._user_action_set_includes
394
395 @property
397 """
398
399 Add further includes to this property, if your solver requires some additional
400 routines from other header files.
401
402 """
403 return self._user_solver_includes
404
406 """!
407
408 Return number of steps required to realise the Runge-Kutta scheme
409
410 Delegate to ButcherTableau.RungeKutta_steps, which tells us for a given
411 polynomial order _rk_order how many Runge Kutta steps we have to
412 employ.
413
414 """
415 return RungeKutta_steps(self._rk_order)
416
418 """
419
420 Add further includes to this property, if your action sets require some additional
421 routines from other header files.
422
423 """
424 self._user_action_set_includes += value
425
426 def add_user_solver_includes(self, value):
427 """
428
429 Add further includes to this property, if your solver requires some additional
430 routines from other header files.
431
432 """
433 self._user_solver_includes += value
434
435 @abstractmethod
437 """
438
439 Recall in subclasses if you wanna change the number of unknowns
440 or auxiliary variables. See class description's subsection on
441 data flow.
442
443 :: Call order and ownership
444
445 This operation can be called multiple times. However, only the very
446 last call matters. All previous calls are wiped out.
447
448 If you have a hierarchy of solvers, every create_data_structure()
449 should first(!) call its parent version. This way, you always ensure
450 that all data are in place before you continue to alter the more
451 specialised versions. So it is (logically) a top-down (general to
452 specialised) run through all create_data_structure() variants
453 within the inheritance tree.
454
455
456 :: Arguments
457
458 _patch: Patch (NxNxN)
459 Actual patch data. We use Finite Volumes, so this is
460 always the current snapshot, i.e. the valid data at one point.
461
462 _patch_overlap_old, _patch_overlap_new: Patch (2xNxN)
463 This is a copy/excerpt from the two adjacent finite volume
464 snapshots plus the old data as backup. If I want to implement
465 local timestepping, I don't have to backup the whole patch
466 (see _patch), but I need a backup of the face data to be able
467 to interpolate in time.
468
469 _patch_overlap_update: Patch (2xNxN)
470 This is hte new update. After the time step, I roll this
471 information over into _patch_overlap_new, while I backup the
472 previous _patch_overlap_new into _patch_overlap_old. If I
473 worked with regular meshes only, I would not need this update
474 field and could work directly with _patch_overlap_new. However,
475 AMR requires me to accumulate data within new while I need
476 the new and old data temporarily. Therefore, I employ this
477 accumulation/roll-over data which usually is not stored
478 persistently.
479
480 """
482 (self._patch_size, self._patch_size, self._patch_size),
484 self._unknown_identifier(),
485 )
487 (
489 self._patch_size,
490 self._patch_size,
491 ),
492 self._unknowns,
493 self._unknown_identifier() + "RhsEstimates",
494 )
496 (2 * self._overlap, self._patch_size, self._patch_size),
498 self._unknown_identifier() + "Old",
499 )
501 (2 * self._overlap, self._patch_size, self._patch_size),
503 self._unknown_identifier() + "New",
504 )
506 (2 * self._overlap, self._patch_size, self._patch_size),
508 self._unknown_identifier() + "Update",
509 )
510
511 if self._cell_data_storage == Storage.CallStack:
513 self._patch, "double"
514 )
516 self._patch_estimates, "double"
517 )
518 elif self._cell_data_storage == Storage.Heap:
520 self._patch, "double"
521 )
523 self._patch_estimates, "double"
524 )
525 elif self._cell_data_storage == Storage.SmartPointers:
527 self._patch, "double"
528 )
529 self._patch_estimates.generator = (
531 self._patch_estimates, "double"
532 )
533 )
534
535 if self._face_data_storage == Storage.CallStack:
537 self._patch_overlap_old, "double"
538 )
540 self._patch_overlap_new, "double"
541 )
543 self._patch_overlap_update, "double"
544 )
545 elif self._face_data_storage == Storage.Heap:
546 self._patch_overlap_old.generator = (
548 self._patch_overlap_old, "double"
549 )
550 )
551 self._patch_overlap_new.generator = (
553 self._patch_overlap_new, "double"
554 )
555 )
556 self._patch_overlap_update.generator = (
558 self._patch_overlap_update, "double"
559 )
560 )
561 elif self._face_data_storage == Storage.SmartPointers:
562 self._patch_overlap_old.generator = (
564 self._patch_overlap_old, "double"
565 )
566 )
567 self._patch_overlap_new.generator = (
569 self._patch_overlap_new, "double"
570 )
571 )
572 self._patch_overlap_update.generator = (
574 self._patch_overlap_update, "double"
575 )
576 )
577 else:
578 assert False, "storage variant {} not supported".format(
580 )
581
582 self._patch_overlap_old.generator.merge_method_definition = (
583 peano4.toolbox.blockstructured.get_face_merge_implementation(
585 )
586 )
587 self._patch_overlap_new.generator.merge_method_definition = (
588 peano4.toolbox.blockstructured.get_face_merge_implementation(
590 )
591 )
592 self._patch.generator.merge_method_definition = (
593 peano4.toolbox.blockstructured.get_cell_merge_implementation(self._patch)
594 )
595
596 self._patch_overlap_old.generator.includes += """
597#include "peano4/utils/Loop.h"
598#include "repositories/SolverRepository.h"
599"""
600 self._patch_overlap_new.generator.includes += """
601#include "peano4/utils/Loop.h"
602#include "repositories/SolverRepository.h"
603"""
604
605 self._patch.generator.load_store_compute_flag = "::peano4::grid::constructLoadStoreComputeFlag({},{},{})".format(
609 )
610
611 self._patch_estimates.generator.load_store_compute_flag = "::peano4::grid::constructLoadStoreComputeFlag({},{},{})".format(
615 )
616
617 self._patch_overlap_old.generator.send_condition = "false"
618 self._patch_overlap_old.generator.receive_and_merge_condition = "false"
619
620 self._patch_overlap_new.generator.send_condition = "false"
621 self._patch_overlap_new.generator.receive_and_merge_condition = "false"
622
623 self._patch_overlap_update.generator.send_condition = "false"
624 self._patch_overlap_update.generator.receive_and_merge_condition = "false"
625
626 self._patch_overlap_old.generator.load_store_compute_flag = "::peano4::grid::constructLoadStoreComputeFlag({},{},{})".format(
630 )
631
632 self._patch_overlap_new.generator.load_store_compute_flag = "::peano4::grid::constructLoadStoreComputeFlag({},{},{})".format(
636 )
637
638 self._patch_overlap_update.generator.load_store_compute_flag = "::peano4::grid::constructLoadStoreComputeFlag({},{},{})".format(
640 "false",
641 "false"
642 )
643
644 self._patch.generator.includes += """
645#include "../repositories/SolverRepository.h"
646"""
647 self._patch_estimates.generator.includes += """
648#include "../repositories/SolverRepository.h"
649"""
650 self._patch_overlap_old.generator.includes += """
651#include "../repositories/SolverRepository.h"
652"""
653 self._patch_overlap_new.generator.includes += """
654#include "../repositories/SolverRepository.h"
655"""
656 self._patch_overlap_update.generator.includes += """
657#include "../repositories/SolverRepository.h"
658"""
659
660 self._cell_label = exahype2.grid.create_cell_label(self._name)
661 self._face_label = exahype2.grid.create_face_label(self._name)
662
663 @abstractmethod
665 """!
666
667 Create required action sets
668
669 Overwrite in subclasses if you wanna create different
670 action sets.
671
672 ## Call order and ownership
673
674 This operation can be called multiple times. However, only the very
675 last call matters. All previous calls are wiped out.
676
677 If you have a hierarchy of solvers, every create_data_structure()
678 should first(!) call its parent version. This way, you always ensure
679 that all data are in place before you continue to alter the more
680 specialised versions. So it is (logically) a top-down (general to
681 specialised) run through all create_data_structure() variants
682 within the inheritance tree.
683
684 :: Recreation vs backup (state discussion)
685
686 We faced some issues with action sets that should not be
687 overwritten. For example, the postprocessing should not be overwritten
688 as users might want to set it and then later on reset the number of
689 unknowns, e.g. In this case, you would loose your postprocessing if
690 create_action_sets() recreated them. So I decided to make an exception
691 here: the postprocessing step is never overwritten by the standard
692 classes.
693
694 There are further action sets which have a state, which users might
695 want to alter. The most prominent one is the AMR routines, where users
696 often alter the interpolation and restriction scheme. Here, things are
697 tricky: If we keep the default one, the action set would not pick up
698 if you changed the number of unknowns, e.g. However, if we recreated
699 the action set, we'd miss out on any changed interpolation/restriction
700 scheme. Therefore, I have to hold the interpolation and restriction
701 scheme seperatae.
702
703 ## Injecting your own guards
704
705 If you inject your own guards, you should combine them with a storage
706 predicate, i.e. _store_cell_data_default_guard() and
707 _load_cell_data_default_guard(). The action sets themselves will not
708 combine the guard with further boolean expressions. Also, you have to
709 study carefully if a predicate accepts a unique guard or a set of
710 guards.
711
712
713 ## Priorities
714
715 All action sets are given the right (default) priorities in this step.
716 You can alter them in subclasses, but it might be more appropriate to
717 set the priorities of your own action sets relative to the existing
718 ones using self._baseline_action_set_descend_invocation_order.
719
720 """
723 self, self._store_cell_data_default_guard(), "true"
724 )
725 )
728 self, self._store_cell_data_default_guard(), "false"
729 )
730 )
732 solver=self,
734 build_up_new_refinement_instructions=True,
735 implement_previous_refinement_instructions=True,
736 called_by_grid_construction=False,
737 )
740 solver=self,
742 build_up_new_refinement_instructions=False,
743 implement_previous_refinement_instructions=True,
744 called_by_grid_construction=False,
745 )
746 )
749 solver=self,
751 build_up_new_refinement_instructions=True,
752 implement_previous_refinement_instructions=True,
753 called_by_grid_construction=True,
754 )
755 )
759 )
760 )
763 )
767 )
768 )
773 False,
776 )
777 )
779 solver=self,
780 interpolation=self._interpolation,
781 restriction=self._restriction,
782 )
783
786 solver=self, guard=self._store_cell_data_default_guard()
787 )
788 )
789
790 # This one is different: it is the hook-in point for other solvers, so we create it
791 # only once if it does not exist
792 if self._action_set_postprocess_solution == None:
795 )
796 if self._action_set_preprocess_solution == None:
799 )
800
803
804 self._action_set_initial_conditions.descend_invocation_order = (
806 )
807 self._action_set_initial_conditions_for_grid_construction.descend_invocation_order = (
809 )
812 )
813 self._action_set_copy_new_faces_onto_old_faces.descend_invocation_order = (
815 )
816 self._action_set_roll_over_update_of_faces.descend_invocation_order = (
818 )
819 self._action_set_update_face_label.descend_invocation_order = (
821 )
822 self._action_set_update_cell_label.descend_invocation_order = (
824 )
825 self._action_set_handle_boundary.descend_invocation_order = (
827 )
828 self._action_set_compute_final_linear_combination.descend_invocation_order = (
830 )
831 self._action_set_project_patch_onto_faces.descend_invocation_order = (
833 )
834 self._action_set_AMR.descend_invocation_order = self._baseline_action_set_descend_invocation_order + 6
835 self._action_set_postprocess_solution.descend_invocation_order = (
837 ) # will be in touch last
838 self._action_set_preprocess_solution.descend_invocation_order = (
840 )
841
843
844
846 return (
847 "not marker.willBeRefined() "
848 + "and repositories::"
850 + ".getSolverState()!="
851 + self._name
852 + "::SolverState::GridConstruction"
853 )
854
856 return (
857 "not marker.willBeRefined() "
858 + "and repositories::"
860 + ".getSolverState()!="
861 + self._name
862 + "::SolverState::GridConstruction"
863 )
864
866 """!
867
868 Extend the guard via ands only. Never use an or, as subclasses might
869 extend it as well, and they will append further ends.
870
871 """
872 return (
873 "not marker.willBeRefined() "
874 + "and repositories::"
876 + ".getSolverState()!="
877 + self._name
878 + "::SolverState::GridConstruction"
879 )
880
882 """!
883
884 Extend the guard via ands only. Never use an or, as subclasses might
885 extend it as well, and they will append further ends.
886
887 """
888 return (
889 "not marker.hasBeenRefined() "
890 + "and repositories::"
892 + ".getSolverState()!="
893 + self._name
894 + "::SolverState::GridConstruction "
895 + "and repositories::"
897 + ".getSolverState()!="
898 + self._name
899 + "::SolverState::GridInitialisation"
900 )
901
903 """!
904
905 Extend the guard via ands only. Never use an or, as subclasses might
906 extend it as well, and they will append further ends.
907
908 """
909 return (
910 "not marker.willBeRefined() "
911 + "and repositories::"
913 + ".getSolverState()!="
914 + self._name
915 + "::SolverState::GridConstruction"
916 )
917
919 """!
920
921 Extend the guard via ands only. Never use an or, as subclasses might
922 extend it as well, and they will append further ends.
923
924 """
925 return (
926 "not marker.hasBeenRefined() "
927 + "and repositories::"
929 + ".getSolverState()!="
930 + self._name
931 + "::SolverState::GridConstruction "
932 + "and repositories::"
934 + ".getSolverState()!="
935 + self._name
936 + "::SolverState::GridInitialisation"
937 )
938
940 return self._name + "Q"
941
943 return "instanceOf" + self._name
944
945 def add_to_Peano4_datamodel(self, datamodel, verbose):
946 """!
947
948 Add all required data to the Peano4 project's datamodel
949 so it is properly built up
950
951 """
952 if verbose:
953 print("Patch data")
954 print("----------")
955 print(str(self._patch))
956 print("Patch overlap data")
957 print("----------")
958 print(str(self._patch_overlap_old))
959 print("Patch overlap data")
960 print("----------")
961 print(str(self._patch_overlap_new))
962 datamodel.add_cell(self._cell_label)
963 datamodel.add_cell(self._patch)
964 datamodel.add_cell(self._patch_estimates)
965 datamodel.add_face(self._patch_overlap_old)
966 datamodel.add_face(self._patch_overlap_new)
967 datamodel.add_face(self._patch_overlap_update)
968 datamodel.add_face(self._face_label)
969
971 """
972 Tell Peano what data to move around
973
974 Inform Peano4 step which data are to be moved around via the
975 use_cell and use_face commands. This operation is generic from
976 ExaHyPE's point of view, i.e. I use it for all grid sweep types.
977
978 """
979 step.use_cell(self._patch)
980 step.use_cell(self._patch_estimates)
981 step.use_cell(self._cell_label)
982 step.use_face(self._patch_overlap_old)
983 step.use_face(self._patch_overlap_new)
984 step.use_face(self._patch_overlap_update)
985 step.use_face(self._face_label)
986
988 return """
989#include "tarch/la/Vector.h"
990
991#include "peano4/utils/Globals.h"
992#include "peano4/utils/Loop.h"
993
994#include "repositories/SolverRepository.h"
995"""
996
998 """!
999
1000 Add your actions to init grid
1001
1002 The AMR stuff has to be the very first thing. Actually, the AMR routines'
1003 interpolation doesn't play any role here. But the restriction indeed is
1004 very important, as we have to get the face data for BCs et al. The action
1005 set order is inverted while we ascend within the tree again. Therefore, we
1006 add the AMR action set first which means it will be called last when we go
1007 from fine to coarse levels within the tree.
1008
1009 The projection onto the faces is a postprocessing step. This is different
1010 to DG, where we need face data for the current time step's solution. Here,
1011 we ensure that all halos are valid for the subsequent time step again.
1012
1013 ## Ordering
1014
1015 The order of the action sets is preserved throughout the steps down within
1016 the tree hierarchy. It is inverted throughout the backrolling.
1017
1018 This is what we want to achieve:
1019
1020 - Project solution onto faces. This happens in touchCellLastTime(). See
1021 exahype2.solvers.rkfd.actionsets.ProjectPatchOntoFaces for comments.
1022 The project will end up in QUpdate.
1023 - Roll updates over on the faces from QUpdate into Q_new. This is done
1024 by RollOverUpdateFace, which requires in turn that the getUpdated()
1025 flag is set. As the roll-over plugs into touchFaceLastTime(), it will
1026 always be called after the projection, since the latter is a cell
1027 operation.
1028 - Copy new face data Q_new into old face data Q_old, as this is the
1029 initial sweep, i.e. the old face data otherwise might hold rubbish.
1030 This is a backup operation realised through the action set
1031 BackupPatchOverlap. This one plugs into touchFaceLastTime() too.
1032 Therefore, it is important that its priority is smaller than the one
1033 of the roll-over, since we invert the call order for the touch-last
1034 events.
1035 - Restrict the data to the coarser level if we are on a hanging face.
1036
1037
1038 """
1039 assert (
1040 self._action_set_copy_new_faces_onto_old_faces.descend_invocation_order
1041 < self._action_set_roll_over_update_of_faces.descend_invocation_order
1042 )
1043
1044 step.add_action_set(self._action_set_preprocess_solution)
1045 step.add_action_set(
1047 )
1048 step.add_action_set(self._action_set_copy_new_faces_onto_old_faces)
1049 step.add_action_set(self._action_set_roll_over_update_of_faces)
1050 step.add_action_set(self._action_set_initial_conditions)
1051 step.add_action_set(self._action_set_project_patch_onto_faces)
1052 step.add_action_set(self._action_set_update_face_label)
1053 step.add_action_set(self._action_set_update_cell_label)
1054 step.add_action_set(self._action_set_AMR)
1055 step.add_action_set(self._action_set_postprocess_solution)
1056
1057 def add_actions_to_create_grid(self, step, evaluate_refinement_criterion):
1058 """
1059
1060 The boundary information is set only once. It is therefore important that
1061 we ues the face label and initialise it properly.
1062
1063 """
1065 step.add_action_set(self._action_set_update_face_label)
1066 step.add_action_set(self._action_set_update_cell_label)
1067 if evaluate_refinement_criterion:
1068 step.add_action_set(self._action_set_AMR_throughout_grid_construction)
1069 else:
1070 step.add_action_set(self._action_set_AMR_commit_without_further_analysis)
1071
1072 @property
1074 self,
1075 ):
1076 return self._plot_description
1077
1078 @plot_description.setter
1079 def plot_description(self, description):
1080 """
1081
1082 Use this one to set a description within the output patch file that tells
1083 the vis solver what the semantics of the entries are. Typicallly, I use
1084 a comma-separated list here.
1085
1086 """
1087 self._plot_description = description
1088
1089 def add_actions_to_plot_solution(self, step, output_path):
1090 """!
1091
1092 Add action sets to plotting grid sweep
1093
1094 Consult the discussion in add_actions_to_init_grid() around the order
1095 of the individual action sets.
1096
1097
1098 ## Adaptive mesh refinement
1099
1100 It is important that we have the coupling/dynamic AMR part in here, as
1101 there might be pending AMR refinement requests that now are realised.
1102 For the same reason, we need the update of the face label and the update
1103 of the cell label in here: The AMR might just propagate over into the
1104 plotting, i.e. we might create new grid entities throughout the plot.
1105 These entities (faces and cells) have to be initialised properly.
1106 Otherwise, their un-initialised data will propagate through to the next
1107 time step.
1108
1109 To make the restriction work, we have to project the solutions onto the
1110 faces.
1111
1112
1113 ## Roll over face data
1114
1115 Das ist ein Kaese, weil das nur einspringt, wenn project wahr ist
1116
1117
1118 """
1119 d = {}
1122
1123 step.add_action_set(
1125 )
1126 # step.add_action_set( self._action_set_roll_over_update_of_faces )
1127 step.add_action_set(self._action_set_update_face_label)
1128 step.add_action_set(self._action_set_update_cell_label)
1129 # step.add_action_set( self._action_set_project_patch_onto_faces )
1130
1132 filename=output_path + "solution-" + self._name,
1133 patch=self._patch,
1134 dataset_name=self._unknown_identifier(),
1135 description=self._plot_description,
1136 guard="repositories::plotFilter.plotPatch(marker) and "
1138 additional_includes="""
1139#include "exahype2/PlotFilter.h"
1140#include "../repositories/SolverRepository.h"
1141""",
1142 precision="PlotterPrecision",
1143 time_stamp_evaluation="0.5*(repositories::getMinTimeStamp()+repositories::getMaxTimeStamp())",
1144 select_dofs=self.select_dofs_to_print,
1145 )
1146 plot_patches_action_set.descend_invocation_order = self._baseline_action_set_descend_invocation_order
1147 step.add_action_set(plot_patches_action_set)
1148
1149 if self._plot_grid_properties:
1150 plot_grid_action_set = peano4.toolbox.PlotGridInPeanoBlockFormat(
1151 filename=output_path + "grid-" + self._name,
1152 cell_unknown=None,
1153 guard="repositories::plotFilter.plotPatch(marker) and "
1155 additional_includes="""
1156#include "exahype2/PlotFilter.h"
1157#include "../repositories/SolverRepository.h"
1158""",
1159 )
1160 plot_grid_action_set.descend_invocation_order = self._baseline_action_set_descend_invocation_order
1161 step.add_action_set(plot_grid_action_set)
1162
1163 pass
1164
1166 """
1167
1168 AMR
1169
1170 It is important that we do the inter-grid transfer operators before we
1171 apply the boundary conditions.
1172
1173 """
1174 d = {}
1177
1178 step.add_action_set(self._action_set_preprocess_solution)
1179 step.add_action_set(
1181 )
1182 step.add_action_set(self._action_set_roll_over_update_of_faces)
1183 step.add_action_set(self._action_set_update_face_label)
1184 step.add_action_set(self._action_set_update_cell_label)
1185 step.add_action_set(self._action_set_handle_boundary)
1186 step.add_action_set(self._action_set_update_cell)
1187 step.add_action_set(self._action_set_project_patch_onto_faces)
1188 step.add_action_set(self._action_set_compute_final_linear_combination)
1189 step.add_action_set(self._action_set_AMR)
1190 step.add_action_set(self._action_set_postprocess_solution)
1191
1192 @abstractmethod
1196 def add_implementation_files_to_project(self, namespace, output, dimensions, subdirectory=""):
1197 """
1198
1199 The ExaHyPE2 project will call this operation when it sets
1200 up the overall environment.
1201
1202 This routine is typically not invoked by a user.
1203
1204 output: peano4.output.Output
1205
1206 """
1207 templatefile_prefix = os.path.dirname(os.path.realpath(__file__)) + "/"
1208
1209 if self._solver_template_file_class_name is None:
1210 templatefile_prefix += self.__class__.__name__
1211 else:
1212 templatefile_prefix += self._solver_template_file_class_name
1213
1214 if(subdirectory):
1215 subdirectory += "/"
1216
1217 abstractHeaderDictionary = {}
1218 implementationDictionary = {}
1219 self._init_dictionary_with_default_parameters(abstractHeaderDictionary)
1220 self._init_dictionary_with_default_parameters(implementationDictionary)
1221 self.add_entries_to_text_replacement_dictionary(abstractHeaderDictionary)
1222 self.add_entries_to_text_replacement_dictionary(implementationDictionary)
1223
1224 generated_abstract_header_file = (
1226 templatefile_prefix + "Abstract.template.h",
1227 templatefile_prefix + "Abstract.template.cpp",
1228 "Abstract" + self._name,
1229 namespace,
1230 subdirectory + ".",
1231 abstractHeaderDictionary,
1232 True,
1233 True,
1234 )
1235 )
1236 generated_solver_files = (
1238 templatefile_prefix + ".template.h",
1239 templatefile_prefix + ".template.cpp",
1240 self._name,
1241 namespace,
1242 subdirectory + ".",
1243 implementationDictionary,
1244 False,
1245 True,
1246 )
1247 )
1248
1249 output.add(generated_abstract_header_file)
1250 output.add(generated_solver_files)
1251 output.makefile.add_h_file(subdirectory + "Abstract" + self._name + ".h", generated=True)
1252 output.makefile.add_h_file(subdirectory + self._name + ".h", generated=True)
1253 output.makefile.add_cpp_file(subdirectory + "Abstract" + self._name + ".cpp", generated=True)
1254 output.makefile.add_cpp_file(subdirectory + self._name + ".cpp", generated=True)
1255
1256 if self._use_var_shortcut:
1257 generated_shortcut_file = peano4.output.Jinja2TemplatedHeaderFile(
1258 os.path.dirname(os.path.realpath(__file__))
1259 + "/"
1260 + "../VariableShortcuts.template.h",
1261 "VariableShortcuts",
1262 namespace,
1263 subdirectory + ".",
1264 implementationDictionary,
1265 True,
1266 True,
1267 )
1268 output.add(generated_shortcut_file)
1269 output.makefile.add_h_file(subdirectory + "VariableShortcuts.h", generated=True)
1270
1271 def set_solver_constants(self, datastring):
1272 self._solver_constants = datastring
1273
1274 def add_solver_constants(self, datastring):
1275 self._solver_constants += datastring
1276
1278 """
1279 This one is called by all algorithmic steps before I invoke
1280 add_entries_to_text_replacement_dictionary().
1281
1282 See the remarks on set_postprocess_updated_patch_kernel to understand why
1283 we have to apply the (partially befilled) dictionary to create a new entry
1284 for this very dictionary.
1285 """
1286 d["NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS"] = self._patch.dim[0]
1287 d["OVERLAP"] = self._overlap
1288 d["HALO_SIZE"] = int(self._patch_overlap_old.dim[0] / 2)
1289 d["RK_ORDER"] = self._rk_order
1290 d["RK_STEPS"] = self.number_of_Runge_Kutta_steps()
1291 d["SOLVER_INSTANCE"] = self.get_name_of_global_instance()
1292 d["SOLVER_NAME"] = self._name
1293 d["UNKNOWN_IDENTIFIER"] = self._unknown_identifier()
1294 d["NUMBER_OF_UNKNOWNS"] = self._unknowns
1295 d["NUMBER_OF_AUXILIARY_VARIABLES"] = self._auxiliary_variables
1296 d["SOLVER_NUMBER"] = 22
1297
1298 d["ASSERTION_WITH_1_ARGUMENTS"] = "nonCriticalAssertion1"
1299 d["ASSERTION_WITH_2_ARGUMENTS"] = "nonCriticalAssertion2"
1300 d["ASSERTION_WITH_3_ARGUMENTS"] = "nonCriticalAssertion3"
1301 d["ASSERTION_WITH_4_ARGUMENTS"] = "nonCriticalAssertion4"
1302 d["ASSERTION_WITH_5_ARGUMENTS"] = "nonCriticalAssertion5"
1303 d["ASSERTION_WITH_6_ARGUMENTS"] = "nonCriticalAssertion6"
1304
1305 if self._min_meshcell_h > self._max_meshcell_h:
1306 raise Exception("min/max h are inconsistent")
1307 d["MAX_GRID_CELL_H"] = self._max_meshcell_h
1308 d["MIN_GRID_CELL_H"] = self._min_meshcell_h
1309
1310 d["SOLVER_CONSTANTS"] = self._solver_constants
1311
1312 d["SOLVER_INCLUDES"] = self.user_solver_includes
1313
1314 d[
1315 "BOUNDARY_CONDITIONS_IMPLEMENTATION"
1317 d[
1318 "REFINEMENT_CRITERION_IMPLEMENTATION"
1320 d["INITIAL_CONDITIONS_IMPLEMENTATION"] = self._initial_conditions_implementation
1321
1322 d["COMPUTE_KERNEL_CALL"] = jinja2.Template(
1323 self._compute_kernel_call, undefined=jinja2.DebugUndefined
1324 ).render(**d)
1325
1326 d["ABSTRACT_SOLVER_USER_DECLARATIONS"] = jinja2.Template(
1327 self._abstract_solver_user_declarations, undefined=jinja2.DebugUndefined
1328 ).render(**d)
1329 d["ABSTRACT_SOLVER_USER_DEFINITIONS"] = jinja2.Template(
1330 self._abstract_solver_user_definitions, undefined=jinja2.DebugUndefined
1331 ).render(**d)
1332 d["SOLVER_USER_DECLARATIONS"] = jinja2.Template(
1333 self._solver_user_declarations, undefined=jinja2.DebugUndefined
1334 ).render(**d)
1335 d["SOLVER_USER_DEFINITIONS"] = jinja2.Template(
1336 self._solver_user_definitions, undefined=jinja2.DebugUndefined
1337 ).render(**d)
1338 d["START_TIME_STEP_IMPLEMENTATION"] = jinja2.Template(
1339 self._start_time_step_implementation, undefined=jinja2.DebugUndefined
1340 ).render(**d)
1341 d["FINISH_TIME_STEP_IMPLEMENTATION"] = jinja2.Template(
1342 self._finish_time_step_implementation, undefined=jinja2.DebugUndefined
1343 ).render(**d)
1344 d["CONSTRUCTOR_IMPLEMENTATION"] = jinja2.Template(
1345 self._constructor_implementation, undefined=jinja2.DebugUndefined
1346 ).render(**d)
1347
1348 d["PREPROCESS_RECONSTRUCTED_PATCH"] = jinja2.Template(
1349 self._preprocess_reconstructed_patch, undefined=jinja2.DebugUndefined
1350 ).render(**d)
1351 d["POSTPROCESS_UPDATED_PATCH"] = jinja2.Template(
1352 self._postprocess_updated_patch, undefined=jinja2.DebugUndefined
1353 ).render(**d)
1354
1355 d["COMPUTE_TIME_STEP_SIZE"] = jinja2.Template(
1356 self._compute_time_step_size, undefined=jinja2.DebugUndefined
1357 ).render(**d)
1358 d["COMPUTE_NEW_TIME_STEP_SIZE"] = jinja2.Template(
1359 self._compute_new_time_step_size, undefined=jinja2.DebugUndefined
1360 ).render(**d)
1361 d["COMPUTE_MAX_EIGENVALUE"] = self._compute_eigenvalue
1362
1363 d["KERNEL_NAMESPACE"] = self._kernel_namespace
1364
1365 d["USE_VARIABLE_SHORTCUT"] = self._use_var_shortcut
1366 d["VARIABLE_NAMES"] = self._variable_names
1367 d["VARIABLE_POSITIONS"] = self._variable_pos
1368
1369 @property
1370 def unknowns(self):
1371 return self._unknowns
1372
1373 @property
1374 def patch_size(self):
1375 return self._patch_size
1376
1377 @property
1379 return self._auxiliary_variables
1380
1381 @patch_size.setter
1382 def patch_size(self, value):
1383 self._patch_size = value
1385 self.create_action_sets()
1386
1387 @unknowns.setter
1388 def unknowns(self, value):
1389 self._unknowns = value
1391 self.create_action_sets()
1392
1393 @auxiliary_variables.setter
1394 def auxiliary_variables(self, value):
1395 self._auxiliary_variables = value
1397 self.create_action_sets()
1398
1399 @property
1403 @preprocess_reconstructed_patch.setter
1405 """!
1406
1407 Please consult exahype2.solvers.fv.FV.preprocess_reconstructed_patch() for
1408 a documentation on this routine.
1409
1410 """
1413 self.create_action_sets()
1414
1415 @property
1416 def name(self):
1417 return self._name
1418
1419 @property
1423 @postprocess_updated_patch.setter
1424 def postprocess_updated_patch(self, kernel):
1425 """
1426
1427 Define a postprocessing routine over the data
1428
1429 The postprocessing kernel often looks similar to the following code:
1430
1431 {
1432 int index = 0;
1433 dfor( volume, {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}} ) {
1434 enforceCCZ4constraints( newQ+index );
1435 index += {{NUMBER_OF_UNKNOWNS}} + {{NUMBER_OF_AUXILIARY_VARIABLES}};
1436 }
1437 }
1438
1439
1440 Within this kernel, you have at least the following variables available:
1441
1442 - newQ This is a pointer to the whole data structure (one large
1443 array).
1444 The patch is not supplemented by a halo layer.
1445 - oldQWithHalo This is a pointer to the data snapshot before the
1446 actual update. This data is combined with the halo layer, i.e. if you
1447 work with 7x7 patches and a halo of 2, the pointer points to a 11x11
1448 patch.
1449 - marker
1450
1451 Furthermore, you can use all the symbols (via Jinja2 syntax) from
1452 _init_dictionary_with_default_parameters().
1453
1454 kernel: String
1455 C++ code that holds the postprocessing kernel
1456
1457 """
1458 self._postprocess_updated_patch = kernel
1460 self.create_action_sets()
1461
1462 @property
1463 def overlap(self):
1464 return self._overlap
1465
1466 @overlap.setter
1467 def overlap(self, value):
1468 if value < 1:
1469 raise Exception(
1470 "Halo (overlap) size has to be bigger than zero but was {}".format(
1471 value
1472 )
1473 )
1474 self._overlap = value
1476 self.create_action_sets()
1477
1478 @property
1479 def interpolation(self):
1480 return self._interpolation
1481
1482 @interpolation.setter
1483 def interpolation(self, value):
1484 """
1485
1486 Set the interpolation scheme. If you rely on a built-in operation, then this
1487 call is all you have to do. Some ExaHyPE solvers however require each solver
1488 to provide special matrices/operators for some interpolation/restriction
1489 variants. If this is the case, you still have to add these matrices manually
1490 to your solver.
1491
1492 """
1493 self._interpolation = value
1494
1496 self.create_action_sets()
1497
1498 @property
1499 def restriction(self):
1500 """
1501
1502 Set the restriction scheme. If you rely on a built-in operation, then this
1503 call is all you have to do. Some ExaHyPE solvers however require each solver
1504 to provide special matrices/operators for some interpolation/restriction
1505 variants. If this is the case, you still have to add these matrices manually
1506 to your solver.
1507
1508 """
1509 return self._restriction
1510
1511 @restriction.setter
1512 def restriction(self, value):
1513 self._restriction = value
1514
1516 self.create_action_sets()
1517
1519 self,
1520 cell_data_storage: Storage,
1521 face_data_storage: Storage,
1522 ):
1523 """
1524
1525 By default, we hold all data on the call stacks. You can explicitly switch
1526 to storage on the heap via smart pointers.
1527
1528 @see create_data_structures()
1529
1530 """
1531 assert isinstance(cell_data_storage, Storage)
1532 assert isinstance(face_data_storage, Storage)
1533
1534 self._cell_data_storage = cell_data_storage
1535 self._face_data_storage = face_data_storage
1536
1538 self.create_action_sets()
Update the cell label within a sweep.
Definition CellLabel.py:9
add_actions_to_create_grid(self, step, evaluate_refinement_criterion)
The boundary information is set only once.
add_actions_to_plot_solution(self, step, output_path)
Add action sets to plotting grid sweep.
add_user_solver_includes(self, value)
Add further includes to this property, if your solver requires some additional routines from other he...
user_action_set_includes(self)
Add further includes to this property, if your action sets require some additional routines from othe...
add_implementation_files_to_project(self, namespace, output, dimensions, subdirectory="")
The ExaHyPE2 project will call this operation when it sets up the overall environment.
add_to_Peano4_datamodel(self, datamodel, verbose)
Add all required data to the Peano4 project's datamodel so it is properly built up.
switch_storage_scheme(self, Storage cell_data_storage, Storage face_data_storage)
By default, we hold all data on the call stacks.
__init__(self, name, patch_size, overlap, rk_order, unknowns, auxiliary_variables, min_meshcell_h, max_meshcell_h, plot_grid_properties, kernel_namespace, baseline_action_set_descend_invocation_order=0)
name: string A unique name for the solver.
number_of_Runge_Kutta_steps(self)
Return number of steps required to realise the Runge-Kutta scheme.
user_solver_includes(self)
Add further includes to this property, if your solver requires some additional routines from other he...
_init_dictionary_with_default_parameters(self, d)
This one is called by all algorithmic steps before I invoke add_entries_to_text_replacement_dictionar...
create_data_structures(self)
Recall in subclasses if you wanna change the number of unknowns or auxiliary variables.
add_user_action_set_includes(self, value)
Add further includes to this property, if your action sets require some additional routines from othe...
The handling of (dynamically) adaptive meshes for finite differences.
Definition DynamicAMR.py:11
The global periodic boundary conditions are set in the Constants.h.
PostprocessSolution differs from other action sets, as I only create it once.
PreprocessSolution differs from other action sets, as I only create it once.
Project patch data onto faces, so the faces hold valid data which can feed into the subsequent Runge-...
Very simple converter which maps the patch 1:1 onto a double array.
This class plugs into the first access to a face and copies the data over.