 |
Peano
|
|
Peano
|
Some applications of ExaHyPE have to add additional mesh sweeps.
This can be due to added functionality or as they want to couple the ExaHyPE solver with another code. The latter approach is the most flexible use case. We therefore discuss this one here, as it covers the "simpler" ambition to add only an additional mesh sweep to an existing ExaHyPE solver. Therefore, our dicussions follows the Peano's generic coupling description. Please read through this discussion before continuing with this section (the description actually refers back to this page at one point).
The way how to add steps is relatively mechanical:
In this example, we plan to manually interact with the step, i.e. to insert C++ code manually into the events once the code is generated. Therefore, we pass in a True here.
This routine can also be invoked earlier. You just have to ensure that all the solvers are already properly added to the ExaHyPE project. See the method documentation.
You still have to trigger the actual mesh traversal. For this, you have to alter the main function. ExaHyPE generates a dummy main in you application's root directory which works for most plain ExaHyPE codes. However, it will never overwrite the main if you have already created a bespoke one.
Before you start manipulating the main, please have a look into the generated file repositories/StepRepository.h. The enum Steps should already contain an entry for your additional algorithmic step once you have followed the instructions from the previous section. That is, Peano knows that the there are all the ExaHyPE mesh traversals and then there's a new, additional one.
We basically let all the logic when to plot and how to pick time steps in place, but before the solver wants to do the very first grid sweep of a new time step, we postpone this step and instead run our new additional step (which we creatively named AdditionalStep in this example). Nothing stops you to extend this logic such that multiple additional steps are executed in a row. You can also obviously plug into the stage just after the solver initialisation:
The statement
runs the additional step directly after the initialisation sweep and then after each time step. If you write
the initialisation will be followed by a first time step (typically with time step size 0 as we have to fiddle out the admissible time step size). After this first one, it injects an additional sweep and from hereon one after each further time step.
The only thing you should not do is a statement similar to
The predicate isFirstGridSweepOfTimeStep() becomes true right after the first beginIteration() and stays on this value. See its documentation. That means that it is true after a primary sweep (if you have an enclave solver) or the first trial in a Runge-Kutta scheme.
As all solvers with a CFL evaluation will run an initial time step with time step size zero, plugging into the mesh traversal after the first time step has completed is sufficient. See remark above.
The important tiny detail here is the call to suspendSolversForOneGridSweep(), which will invoke suspendSolversForOneGridSweep() on each and every ExaHyPE solver. This tells them that we have a "spare" mesh traversal - logically similar to a plot - where they don't have to compute anything.
This modificated logic squeezes the additional steps in-between any two complete time steps. If a solver requires multiple mesh sweeps to realise its update, this sequence of sweeps will not be interrupted. While the realisation as sketched so far works on a single core, it very likely will crash immediately once you use domain decomposition: All ExaHyPE solvers make quite distinct decisions whether to send out data in any single grid sweep or not. Any decision has to be matched with the corresponind receive logic: If we send out data in iteration n, these data have to be reveived in iteration n+1.
Rather than modifying this complex logic, we observe that ExaHyPE has already defined one type of mesh traversal which can always squeeze into any two time steps: the plotting.
In the example above, we tell the new action set that we want to manually add new behaviour, as we set the flag
As a consequence, we will obtain a new class within our project's actionset directory. You can directly add behaviour into this class. Please note that such a class is never overwritten by subsequent Python API invocations.
Adding behaviour manually is only one way to inject behaviour. You might as well set
In this case, no new class in actionsets is generated. You will have to add additional behaviour yourself as detailed below.
Some benchmarks want to enable the additional traversals for some solver choices, while they don't need it for others. In this case, we have identified some best practices. You start from a vanilla main as generated by the Python API, and then you make modifications to this main. Main files are not overwritten by the API if it already finds and existing one, so once your modified main is added to the repository, you can be sure every checkout sees this one even though glue code is to be regenerated on a new system.
To switch the feature no/off, I recommend that you make your Python script (the simulation driver) trigger this line
whenever you need additional mesh traversals. Within the main, you can then protected the code regions of interest with
An example for the realisation of this usage pattern is benchmarks/exahhype2/ccz4/gauge-wave. Study RKDG Solver (Gauge Wave) for a discussion of their particular algorithmic steps.
The canonical way to add new functionality directly via Python follows always the same pattern:
A classic example for a sophisticated variant of the the last step is the action set implemented within applications/exahype2/ccz4/ComputeFirstDerivatives.py. It extends peano4.toolbox.blockstructured.ReconstructPatchAndApplyFunctor for a Finite Difference scheme. This one in turn inherits from peano4.solversteps.ActionSet. Once defined, we again add this to the step via add_action_set().
You can write your own action sets to it (see above), and/or you can use existing action sets within your additional mesh sweep. Each additional mesh sweep has access to all ExaHyPE solver data. Nothing stops you from using and manipulating these data in your additional mesh sweeps.
For this, I recommend that you study your solver objects' interna. Each of them creates a whole suite of solver-specific action sets which are then added to the individual steps via add_actions_to_perform_time_step(). Each solver also creates its own distinct set of data structures and ties them to vertices, faces and cells.
A good example is the ccz4.py script and the action sets within ComputeFirstDerivatives.py which can be found in applications/exahype2/ccz4. Here, we define an additional action set which manipulates the "solution" after each time step or Runge-Kutta step, respectively. It manipulates the solver's core data structure, and is used in combination with various pre-defined action sets such as exahype2.solvers.rkfd.actionsets.ProjectPatchOntoFaces.
Do not coarsen or erase the mesh in the additional grid sweeps. Some ExaHyPE solvers are really picky when it comes to when they can refine or not. A refinement in-between two Runge-Kutta steps for example is not particularly helpful. If you need the additional steps to trigger mesh modifications, use the exahype2::RefinementControl objects and add them to exahype2::RefinementControlService. The service serves as a middle layer in-between ExaHyPE and Peano.
By the time we have configured the new additional steps, you are working with a plain Peano project. You consequently have quite some freedom. One of the things you can do is to add arbitrary data to the mesh vertices, face and cell which have nothing to do with the core ExaHyPE solvers in use. You can augment the whole setup with your own data.
While the routine init_new_user_defined_algorithmic_step() ensures that additional, user-defined algorithmic steps "know" which data ExaHyPE associates with each grid entity (face, mesh, vertex), this does not hold the other way round:
If you want to add further data to the mesh, you first create this data. After that, you add the corresponding use_face(), use_vertex() or use_cell() calls to your script for your new mesh traversal. At this point, you also have to notify all of ExaHyPE's mesh traversals that we have the additional data. They also traverse the mesh and if they are not aware of some user data attached to it, this might lead to data inconsistencies. For this, you have to iterate over all of the steps associated with the Peano project and add the usage information there as well.
Unfortunately, it becomes a little bit nasty here: By the time you add the new data, the ExaHyPE project has already lowered its data representation into Peano. That is, you cannot work with the ExaHyPE data representation anymore (with solvers and their actions), but you have to alter Peano's data representation: Basically, you have to iterate over all mesh sweeps within the Peano project, and ensure that they know about the additional data that you inject into the system.
To realise the second order formulation, you have to invoke the Python script with the -so switch. If this switch is set, we add an additional mesh traversal to the arising Peano project, which we call after each and every Runge-Kutta step, as well as after the very last step.
We now need a main file that can both handle the additional mesh sweep as well as the version without an additional sweep. We don't want to write a main twice. Therefore, we either add
or we don't. This will yield a define pragma in Constants.h (generated). We query this pragma in the code when we decide which step we pick
and also within the switch statement which runs the additional step:
A naive way to implement the activation of the additional mesh sweep is the following code snippet:
This variant kicks off the manual reconstruction prior to the first mesh sweep of the Runge-Kutta scheme. This seems to work and to do the job in our setup. However, it is not enough for higher-order convergence, where we would need a preprocessing step prior to each intermediate step. We ignore this fact here.
Once we have computed the second derivatives into auxiliary variables, we have to project those guys onto the boundary. For this, we hijack the projection of the Runge-Kutta scheme. The documentation of exahype2.solvers.rkfd.actionsets.ProjectPatchOntoFaces requires K+1 predicates for a Runge-Kutta scheme of Kth order. The first K entries correspond to the Runge-Kutta scheme. The last one allows you to do something different for the initialisation. We disable all K-1 entries, and only set the very last option to true. This is clearly a hack, but if you study the generated code, it seems to work.
1.9.1