 |
Peano
|
|
Peano
|
All solvers are declared within the Python file SBH.py. Each solver therein is a subclass of a CCZ4 solver, tailored exactly towards the single black hole setup. Each solver therein is also available to be used from the command line by picking a variant through --solver.
The lionshare of the compute time is burnt within the Finite Volume patches. I followed a series of optimisation steps to get this part of the code faster:
I have written a bespoke orchestration scheme for this solver with benchmarks::exahype2::ccz4::MulticoreOrchestration. This orchestration is similar to hardcoded if you use a single solver, but it can prioritise and distinguish tasks once you use the coupled variant. In this case, only Finite Volume patches for example go to a GPU, while all the FD4 stuff stays on the host.
The section below discusses various simulation-specific implementation choices.
limited Even simple initial conditions can make our higher-order Finite Differences solver fail: The higher order solution representation is not capable to resolve shocks properly (or steep gradients), and the solution therefore starts to oscillate. At the same time, experiments with Finite Volumes yield solutions which a subject to very strong numerical diffusion. The solution quality is not good. The black hole disappears.
There are two fundamental approaches to tackle this dilemma: On the one hand, we can introduce artificial damping terms, i.e. numerical diffusion, in the FD4 scheme. In the context of Finite Differences, this is usually done via Kreiss-Oliger (KO) dissipation. ExaHyPE's FD4 scheme indeed supports KO term. In the present case study, we however switch them off! Instead, we accept that both FD4 and FV have their pros and cons and try to use each solver where it unfolds its full potential.
The image above illustrates our concept: There's an FD4 solver which computes the solution in the whole domain. This solution is is illustrated via the horizontal cut. In an area around the black hole (cross-shaped here; we visualise a slice through this domain), an additional Finite Volume solver is active. In each time step, the FD4 solver sets the boundary conditions of the Finite Volume solver. Those are the dark blue cells around the Finite Volume patches. With these coupling boundary conditions and global Sommerfeld conditions for the FD4 scheme, both solvers advance in time by one time step. After that one, the Finite Volume scheme overwrites the FD4 scheme in points which are covered by both numerical schemes (restriction).
This coupling strategy follows the default coupling schema of ExaHyPE. Different to the generic coupling recipes there, we let the higher-order solver cover the whole domain. As FV only will be used in a tiny region, this overhead is acceptable. Furthermore, we do not use any fancy colouring to decide which solver is active where. We hard-code this information, as we know where our black hole is.
Mathematically, the Finite Volume solver serves as a priori limiter to the FD4 solution. It is a priori, as we do not wait for the FD4 solution to determine if we need Finite Volumes. We overwrite always. For the realisation, we try to keep as much coupling logic as possible within the FD4 scheme: The limiter knows where to compute, but all the remaining logic is realised with the FD4 solver, as this one knows what the data layout of FD4 is. Furthermore, this separation of responsibilities allows us, in theory, to use the same limiter with RKDG as well.
As we combine solvers, we have to ensure that their time step sizes match. This means that the Finite Volume patch size has to be calibrated with \( p^{2}=16 \) relative to the high order's patch size. If they had the same mesh size, the higher order schemes would take significantly smaller time step sizes. However, we want to havee the same spatial and temporal accuacy, so we need a smaller mesh size for the Finite Volume scheme. The file benchmarks/exahype2/ccz4/single-black-hole/convergence-study.py and the solvers within applications/exahype2/ccz/CCZ4Solver.py hold further details.
As the FD4 scheme works with subpatches per octree cell (octant) of size \( 3^d \) or \( 6^d \) respectively, we get Finite Volume patches around the centre which are of the size \( (3\cdot 16)^d = 48^d \) or \( (6\cdot 16)^d = 96^d \). They host up to 884736 volumes. They are large.
Although we recalibrate the patch sizes to take the different spatial discretisation orders into account, this not yet ensure that both the Finite Differences and Finite Volume solver run in-sync. The Finite Volume scheme will be subject to a higher numerical dissipation, while the Finite Differences scheme might run into numerical oscillations/instabilities. At the same time, the limiter will always run at a significantly higher resolution than its higher-order counterpart and capture the singularity. It might therefore happen that it actually yields a smaller admissibile time step size. Therefore, it remains important that both solvers agree on a common time step size after each time step, even though both solvers should come up with an admissible time step size of roughly the same order.
We try to keep all coupling and other logic within the higher order solver. Our intention is to have a generic limiter which can work with multiple solvers, once we have multiple higher order schemes. As a consequence, I introduced a new benchmarks::exahype2::ccz4::CCZ4SBH_FV::reduceAdmissibleTimeStepSize() which allows us to reset the admissible time step size. Next, each higher order solver overwrites startTimeStep(), computes the minimum over the available time step sizes, and sets this minimum as own admissible size as well as the admissible time step size of the limiter.
Our FV solver is, through various inheritance levels, a subclass of peano.exahype2.solvers.fv.EnclaveTasking which in turn inherits from peano.exahype2.solvers.fv.FV. The solver is a text-book FV scheme, with one major added functionality: We add the implementation the routines benchmarks::exahype2::ccz4::CCZ4SBH_FV::isCellOverlappingWithBHImpactArea(), benchmarks::exahype2::ccz4::CCZ4SBH_FV::areBothAdjacentCellsOverlappingWithBHImpactArea(), and their cousins. That is, we implement the actual logic where to apply the limiter within the C++ class.
Within the Python code, we follow localisation recommendations for solver coupling and overwrite the routines
such that only data within a certain radius around the black hole are actually held. An analogous scheme holds for the faces. For the implementation of these guards, we rely on the previously introduced benchmarks::exahype2::ccz4::CCZ4SBH_FV::isCellOverlappingWithBHImpactArea() and benchmarks::exahype2::ccz4::CCZ4SBH_FV::areBothAdjacentCellsOverlappingWithBHImpactArea().
Their implementation is straightforward: We store these cells which overlap with this black hole influence region. Consequently, we store only faces persistently which have a cell left and right which actually holds data. Finally, we have to create faces temporarily - similar to hanging faces - around this Finite Volume area of interest. We need those as the Finite Volume solver will, for example, project its solution onto a patch's boundaries, even though we can throw them away afterwards, i.e. we do not store them.
It is the routine peano.exahype2.solvers.fv.FV.create_action_sets() which configures all the dynamic behaviour. After calling the superclass implementation, we modifiy the two guards _action_set_update_cell.guard and _action_set_merge_enclave_task_outcome.guard to ensure that the solution is also only computed on those cells which hold meaningful data. Without those guards, the solver would invoke the compute kernels everywhere, even though the majority of the mesh holds only garbage solution which is neither stores nor loaded persistently in-between two mesh traversals.
Finally, we set the priorities of the Finite Volume solver up. By default, any enclave tasks has the same (default) priority. We know however that the Finite Volume tasks are extremely expensive - notably compared to their higher order cell countersparts. In an ideal world, all enclave tasks should be postponed into a queue throughout the grid traversal besides the Finite Volume tasks which should go straight onto a core or the GPU. To re-prioritise, we add
to the constructor.
To realise the coupling, we have to augment the FD4 solver by some additional pre- and post-processing steps:
The primary solver takes its reconstructed solution, i.e. the Q values within the patch plus its halo (of size 3 for FD) and projects this information down into the FV solver. This should happen
I initially thought that the projection could be triggered by the FD4 solver after each time step. However, this is not a good idea, as we do not have consistent face data at this point. We however need the faces to have proper halo information to be able to interpolate correctly. So it has be done before the solver kicks off.
Associating it with the FD4 solver's preprocess_reconstructed_solution would be convenient, as these data hold the halo information already, i.e. they provide access to the patch data plus its halo. However, plugging into this routine makes assumptions on the evaluation order: We have to know that the FV code has not started any calculation yet. It is more convenient to let the coupling happen within either solver's _action_set_preprocess_solution action set as a preamble of its actual time step. This way, it will definitely be kicked off before we actually run into the timestepping.
Once the fifth step has terminated, the Finite Volume solver's patches have a properly befilled QNew on the faces. This QNew usually holds the halo copy from neighbouring patches. Where these neighbouring patches do not exist in the FV scheme, we now have interpolated data from the FD4 solution.
Our scheme couples the FD4 solver to the Finite Volume scheme over the faces. However, we realise the coupling volumetric, i.e we project volumetric data from FD4 onto FV data (temporarily) and from there project the volumetric Finite Volume representation again onto the faces. This not particularly elegant, and there might be way faster realisations which omit this temporary volumetric step, but is a realisation which can be implemented with the on-board routines and the logic of it is rather clear.
The coupling from the Finite Volume scheme to FD4 in contrast is a real volumetric coupling: We overwrite the FD4's self._action_set_postprocess_solution with an instance of exahype2.solvers.rkfd.actionsets.PatchWisePostprocessSolution. This postprocessing is be executed after both the Finite Volume and the Differences scheme have updated their patch. This is ensured by the fact that postprocessing steps always plug into touchCellLastTime (cmp Peano's generic discussion of the order of events).
For the realisation, we employ the built-in routines. For example, we call toolbox::blockstructured::restrictCellIntoOverlappingCell_inject() for every cell which holds Finite Volume data, i.e. overlaps with the black hole. This is slightly different to ExaHypE's generic coupling sketch, where I recommend to have at least one whole cell overlap between any two solvers: I do not only inject data from Finite Volume patches back which are completely surrounded by further FV patches. It seemed to me to be a waste of compute power to do this. Rather make the black hole influence area smaller.
We have played around with different restriction variants:
1.9.8