Higher-dimensional PDE Solver Meshes

Higher-dimensional Adaptive Cartesian Meshes for PDE Solvers

Almost all mainstream PDE solvers in my work are 2d or 3d (with a few exceptions). However, a lot of the problems of people using my code are somehow higher-dimensional. The prime example is solvers where we don’t know some parameters precisely. In this case, the thing most people do is to repeatedly run simulations with varying parameters to match the outcomes they are after or to get a feeling how the solution reacts to variations in some input parameters. Within optimisation algorithms, these simulations often run one after the other. They form the forward cycle. After each one, we tweak the material parameters of interest and then start again. Within a Bayesian framework, multiple of these simulations would run for different parameters in parallel to sample the input parameter distribution.

We could visualise such a traditional approach as in the picture below: Several solver instances are run “in parallel” (we can deploy them to different compute jobs, as they are not coupled). Here, the domain is 2d, but we could think of the problem as a 3d thing with the parameter space yielding an additional dimension.

The illustration assumes that all the simulations yield the same adaptivity pattern. If this is the case, we could actually run only a 2d simulation and embed the parameter space into each cell; effectively an additional 1d cell embedded into a 2d cell.

This is a standard technique in climate and weather codes – illustrated above – and it means that we get very high workload per cell, but pay the price for all the overhead such as load balancing only once. Obviously, the scheme has flaws in the context of unknown parameters: The spatial adaptivity pattern has to be the same for all parameters, it is not clear how we handle adaptivity in the parameter dimension, and we notably have to know about all parameter dimensions a priori: while we can use a 3d baseline kit and up to N material parameters, we have to know this N from the outset. A popular extension of this scheme uses an adaptive (3d) mesh for the parameter space and forms a 3d+3d space via a tensor product approach. Can we do better?

Vision

A classic spacetree (octree) algorithm embeds the d-dimensional computational domain into a cube and then refines this cube following some adaptivity criterion. Each refinement steps yields new, smaller cubes equivalent to an adaptive Cartesian mesh. A simple octree is depicted below.

We hypothesize that we can take this formalism and allow it to unfold into arbitrary dimensions. That is, the mesh starts with a d-dimensional adaptive mesh, and then each cell can either refine within these d dimensions or unfold (extend the hypercube) into an additional dimension. They can also collapse (parameter) dimensions. Neither the number of (unfolded) dimensions nor the combinations of which parameter dimensions have been covered has to be known a priori. Further to the unfolding, a mesh can also grow along any dimension: It might initially only cover [0,1] along one direction, but then we let it grow and get additional cells – for example to cover [0,2].

We assume that such an algorithm can host, for example, a wave equation. We don’t have to know how many material parameters (e.g. for a wave running through a heterogeneous medium or subject to some parameterised terms) make a difference and we notably don’t have to take them into account in the solve process/mesh right from the start. Instead, we can study, per cell if needed, if the parameters make a difference and originally start with a sole 3d simulation. As soon as we spot that the parameters make a huge difference to the outcome, we call this parameter a new dimension and unfold the mesh along this new dimension. From this point on, we continue to develop the solution. If, at one point, the parameter makes no significant difference anymore – for example as the wave fades out or leaves the domain – we collapse the domain again.

Preparatory work

We have prepared a mesh management algorithm (and software) which yields data access patterns with very strong spatial and temporal locality. The approach works for any dimension. In principle, this is already brilliant, as the unfolding needs a mesh traversal scheme that does not first run through the 3 main dimensions and then covers all the unfolded ones. Instead, we want to stay local, i.e. if we unfold, we want to handle the additional dimensions around the time we handle a cell. The published approach also has a very low memory footprint, which again suits our needs. The paper is the SISC one

Peano—A Traversal and Storage Scheme for Octree-Like Adaptive Cartesian Multiscale Grids

The software also has a quite mature API which we use, for example, for ExaHyPE or our SPH code. That is, if we had the unfolding/collapsing mesh feature, we could implement a textbook sensitivity analysis over an ExaHyPE benchmark and demonstrate that the idea works. The API got published in TOMS:
The Peano Software—Parallel, Automaton-based, Dynamically Adaptive Grid Traversals

Challenge

• Design and implement the unfolding/collapse algorithm.
• Prototype a simple benchmark on top of it which demonstrates how this approach can be used for sensitivity analysis.
• Realise domain decomposition on top of the dimension-flexible mesh. This is not trivial, as the actual data exchange along cuts in arbitrary-dimensional meshes can be challenging. Obviously, we do not only want to cut in the spatial dimensions.
• Showcase how the new algorithm can be integrated into ExaHyPE.