Towards a new vision of mesh adaptation methods and its impact on the simulation of PDEs

Loïc Gouarin

Marc Massot

The present work is the result of a team work involving

• Thomas Bellotti (CR CNRS - EM2C - Fédé Maths CS)
• Josselin Massot (IR École polytechnique, CMAP)
• Pierre Matalon (IR École polytechnique, CMAP)
• Laurent Séries (IR École polytechnique, CMAP)
• Christian Tenaud (DR CNRS - EM2C - Fédé Maths CS)

site web de l’équipe

https://github.com/hpc-maths/

Context

Burgers equation - small hat problem

\[ \partial_t u + \partial_x \left ( f(u) \right ) = 0, \quad t \geq 0, \quad x \in \mathbb{R}, \qquad f(u) = \dfrac{u^2}{2}, \]

Consider the Cauchy problem with initial cond.

\[ u^0(x)=\left\{% \begin{array}{cc} \hfill 0, & x \in ]- \infty,-1]\cup [1, + \infty[,\\ x+1, & x \in ]-1,0], \hfill \\ 1-x, & x \in[0,1[. \hfill \\ \end{array} \right. \]

• Shock formation at time \(T^* = 1\)
• Leading to irreversible solution
• RH condition governs shock dynamics

Burgers equation - small hat problem

\[ \partial_t u + \partial_x \left ( f(u) \right ) = 0, \quad t \geq 0, \quad x \in \mathbb{R}, \qquad f(u) = \dfrac{u^2}{2}, \]

Consider the Cauchy problem with initial cond.

\[ u^0(x)=\left\{% \begin{array}{cc} \hfill 0, & x \in ]- \infty,-1]\cup [1, + \infty[,\\ x+1, & x \in ]-1,0], \hfill \\ 1-x, & x \in[0,1[. \hfill \\ \end{array} \right. \]

• Shock formation at time \(T^* = 1\)
• Leading to irreversible solution
• RH condition governs shock dynamics

Burgers equation - small hat problem

\[ \partial_t u + \partial_x \left ( f(u) \right ) = 0, \quad t \geq 0, \quad x \in \mathbb{R}, \qquad f(u) = \dfrac{u^2}{2}, \]

Consider the Cauchy problem with initial cond.

\[ u^0(x)=\left\{% \begin{array}{cc} \hfill 0, & x \in ]- \infty,-1]\cup [1, + \infty[,\\ x+1, & x \in ]-1,0], \hfill \\ 1-x, & x \in[0,1[. \hfill \\ \end{array} \right. \]

• Shock location \(\varphi(t)=\sqrt{2(1+t)}-1\)
• Propagation speed shock \(\sigma(t)={1}/{\sqrt{2(1+t)}}\)
• Shock amplitude \([u]=\sqrt{{2}/{(1+t)}}\)

The key issue is to get a analytical solution of the problem involving both reversible and irreversible dynamics. Notebook MAP412 pour intégration Godunov et convergence https://jupyter_map412.gitlab.labos.polytechnique.fr/jb_2024_2025/content/chapter/11_edp_03/burgers.html

Adaptive Multiresolution

• Minimum level \(\underline{\ell}\) and maximum level \(\bar{\ell}\).
• Cells: \[ C_{\ell, k}:=\prod_{\alpha=1}^d\left[2^{-\ell} k_\alpha, 2^{-\ell}\left(k_\alpha+1\right)\right] \]
• Finest step: \(\Delta x=2^{-\bar{\ell}}\).
• Level-wise step: \(\Delta x_{\ell}:=2^{-\ell}=2^{\Delta \ell} \Delta x\).

Just introducing the notations and the notion of level difference with Delta l

Wavelets

Decomposition of the solution on a wavelet basis [Daubechies, ’88], [Mallat, ’89] to measure its local regularity. “Practical” approach by [Harten, ’95], [Cohen et al., ’03].

Projection operator

Prediction operator at order \(2 \gamma+1\)

\[ {\hat f}_{\ell+1,2 k}={f}_{\ell, k}+\sum_{\sigma=1}^\gamma \psi_\sigma\left({f}_{\ell, k+\sigma}-{f}_{\ell, k-\sigma}\right) \]

Details are regularity indicator \[ {\mathrm{d}}_{\ell, {k}}:={f}_{\ell, {k}}-{\hat{f}}_{\ell, {k}} \]

Let \(f \in W^{\nu, \infty}\) (neigh. of \(C_{\ell, k}\) ), then \[ \left|{\mathrm{d}}_{\ell, k}\right| \lesssim 2^{-\ell \min (\nu, 2 \gamma+1)}|f|_{W^{\min (\nu, 2 \gamma+1), \infty}} \]

Fast wavelet transform:

means at the finest level can be recast as means at the coarsest level + details \[ \begin{array}{rlr} {f}_{\overline{\ell}} & \Longleftrightarrow & \left({f}_{\underline{\ell}}, {{d}}_{\underline{\ell} +1}, \ldots, {d}_{\bar{\ell}}\right)\\ \end{array} \]

Mesh coarsening (static)

Local regularity of the solution allows to select areas to coarsen

\[ {{f}}_{\bar{\ell}} \rightarrow \left({f}_{\underline{\ell}}, {\mathbf{d}}_{\underline{\ell}+1}, \ldots, {\mathbf{d}}_{\bar{\ell}}\right) \rightarrow \left({f}_{\underline{\ell}}, {\tilde{\mathbf{d}}}_{\underline{\ell}+1}, \ldots, \tilde{{\mathbf{d}}}_{\bar{\ell}}\right) \rightarrow {\tilde{{f}}}_{\bar{\ell}} \] \[ \tilde{{\mathrm{d}}}_{\ell, k}= \begin{cases}0, & \text { if } \left|{\mathbf{d}}_{\ell, k}\right| \leq \epsilon_{\ell}=2^{-d \Delta \ell} \epsilon, \quad \rightarrow \quad\left\|{\mathbf{f}}_{\bar{\ell}}-\tilde{{\mathbf{f}}}_{\bar{\ell}}\right\|_{\ell^p} \lesssim \epsilon \\ {\mathrm{d}}_{\ell, k}, & \text { otherwise} \end{cases} \]

Set a small (below \(\epsilon_{\ell}\)) detail to zero \(\equiv\) erase the cell \(C_{\ell, k}\) from the structure

Examples

Equation

\[ f(x) = exp(-50x^2) \; \text{for} \; x\in[-1, 1] \]

min level	1
max level	12
ε	10-3
compression rate	96.29%
error	0.00078

Examples

Equation

\[ f(x) = \left\{ \begin{array}{l} 1 - |2x| \; \text{if} \; -0.5 < x < 0.5,\\ 0 \; \text{elsewhere} \end{array} \right. \]

min level	1
max level	12
ε	10-3
compression rate	98.49%
error	0

Examples

Equation

\[ f(x) = 1 - \sqrt{\left| sin \left( \frac{\pi}{2} x \right) \right|} \; \text{for} \; x\in[-1, 1] \]

min level	1
max level	12
ε	10-3
compression rate	96.29%
error	0.00053

Examples

Equation

\[ f(x) = \tanh(50 |x|) - 1 \; \text{for} \; x\in[-1, 1] \]

min level	1
max level	12
ε	10-3
compression rate	97.46%
error	0.002

Time evolution of PDEs

• Finite volumes with global time step \(\Delta t = \Lambda(\Delta x)\)
• Use dynamic mesh refinement

Mesh updated using “old” information at time \(t\) to accommodate the one at time \(t + \Delta t\)

• Propagation of information : add security cells
• Formation of singularities : (regularity index: \(\nu =0\), \(\mu = \min(\nu,2\gamma +1)\)) refine if \[ \left|{\mathbf{d}}_{\ell, k}\right| \geq \epsilon_{\ell}\,2^{d+\mu} \]

Finite volumes / conservation / order

Flux evaluation at interfaces between levels

Using the prediction operator allows to evaluate fluxes at the same level

Finite volume method

We use a Godunov flux for the small hat problem

\(\epsilon = 1e-2\), \(\underline{\ell} = 2\), \(\bar{\ell} = 12\)

Adaptive mesh refinement software

Mesh adaptation

If we look at all the open source software specializing in dynamic mesh adaptation, there are two main families: - patch-based, which is a hierarchical representation of the mesh: layers are placed on top of layers - cell-based, which is a flat representation of the mesh

Each has its advantages and disadvantages: - patch-based has rectangular zones for tiling and optimizing caches. But it generally requires more cells than necessary. - cell-based requires far fewer cells but requires a tree-like structure, which means that you lose the good memory locality you had with patch-based. We use space filling curves such as Morton or Hilbert to find an acceptable locality.

If we look in a little more detail at the functions offered by these software packages, we can group them into 4 main families.

There are two types of data structure, as described above: a list of blocks or a tree.

There are two main adaptation criteria: one is based on a heuristic criterion and depends on the physical problem you are looking at. This could be a gradient, for example. The other is based on a wavelet decomposition that allows you to adapt the mesh without knowing anything about the physical problem, as we have just shown with Haar wavelets and multiresolution.

As you advance in time, you can choose different time steps depending on the resolution of the grid: this is called subcycling. Otherwise, you take the same time step everywhere and, in general, it is the finest grid that guides its value to satisfy a CFL.

Finally, given that the mesh is dynamic, the load balancing must be reviewed regularly during the calculation so that, in a parallel context, all the processes have more or less the same workload. There are two types of method: one based on the space filling curve, where you cut out chunks of the same size following this curve; the other is based on solving a diffusion equation on the workload of the processes.

Open source software

Name	Data structure	Adaptation criteria	Time scheme	Load balancing
AMReX	block	heuristic	global/local	SFC
Dendro	tree	wavelet	global	SFC
Dyablo	tree	heuristic	global	SFC
Peano	tree	-	-	SFC
P4est	tree	-	-	SFC
samurai	interval	heuristic/wavelet	RK/splitting/IMEX time-space/code coupling	SFC/diffusion algorithm

samurai: create a unified framework for testing a whole range
of mesh adaptation methods with the latest generation of numerical schemes.

Here’s an overview of the software we think is interesting to look at today.

Peano and p4est don’t have any adaptation or time-scheme criteria, as they are software programs specializing solely in mesh management.

samurai

Design principles

• Compress the mesh according to the level-wise spatial connectivity along each Cartesian axis.
• Achieve fast look-up for a cell into the structure, especially for parents and neighbors.
• Maximize the memory contiguity of the stored data to allow for caching and vectorization.
• Facilitate inter-level operations which are common in many numerical techniques.

An overview of the data structure

[ start , end [ @ offset

We illustrate the samurai data structure, starting with a 2D problem to see all the components. The data structure is fully recursive and can therefore handle Nd Cartesian meshes.

The main element of samurai is an interval modeled by its beginning and non-included end. There’s also an offset parameter that will serve us in two different ways. We’ll come back to this later.

An overview of the data structure

• Level 0
• y: 0 → [0, 4[
• y: 1 → [0, 1[, [3, 4[
• y: 2 → [0, 1[, [3, 4[
• y: 3 → [0, 3[
• Level 1
• y: 2 → [2, 6[
• y: 3 → [2, 6[
• y: 4 → [2, 4[, [5, 6[
• y: 5 → [2, 6[
• y: 6 → [6, 8[
• y: 7 → [6, 7[
• Level 2
• y: 8 → [8, 10[
• y: 9 → [8, 10[
• y: 14 → [14, 16[
• y: 15 → [14, 16[

cell list

• Level 0
• x: [0, 4[, [0, 1[, [3, 4[, [0, 1[, [3, 4[, [0, 3[
• y: [0, 4[ @0
• y-offset: [0, 1, 3, 5, 6]
• Level 1
• x: [2, 6[, [2, 6[, [2, 4[, [5, 6[, [2, 6[, [6, 8[, [6, 7[
• y: [2, 8[@-2
• y-offset = [0, 1, 2, 4, 5, 6, 7]
• Level 2
• x: [8, 10[, [8, 10[, [14, 16[, [14, 16[
• y: [8, 10[@-8, [14, 16[@-12
• y-offset: [0, 1, 2, 3, 4]

cell array

The first thing we’re going to do is go through each level and each y-component, and look at the set of contiguous cells along the x-axis to construct intervals.

An overview of the data structure

• Level 0
• x: [0, 4[ @ 0, [0, 1[ @ 4, [3, 4[ @ 2, [0, 1[ @ 6, [3, 4[ @ 4, [0, 3[ @ 4
• y: [0, 4[ @ 0
• y-offset: [0, 1, 3, 5, 6]
• Level 1
• x: [2, 6[ @ 9, [2, 6[ @ 13, [2, 4[ @ 17, [5, 6[ @ 16, [2, 6[ @ 20, [6, 8[ @ 20, [6, 7[ @ 22
• y: [2, 8[ @ -2
• y-offset : [0, 1, 2, 4, 5, 6, 7]
• Level 2
• x: [8, 10[ @ 21, [8, 10[ @ 23, [14, 16[ @ 19, [14, 16[ @ 21
• y: [8, 10[ @ -8, [14, 16[ @ -12
• y-offset: [0, 1, 2, 3, 4]

Now that we have this representation, we’ll also try compressing along the y-axis.

The y-offset table lets us know how many x-intervals each y has. This is equivalent to the CSR format, where you have the same thing for rows, where you have to count the number of non-zero elements and accumulate them.

The offset is used here to find the right range of x-intervals in the 1d array of x-intervals.

Mesh constraints

• A refined cell is split into 2 in 1d, 4 in 2d and 8 in 3d equal parts.
• At a given resolution level, the size of the cells is equal.
• The size of the cells is defined by the resolution level. \[\Delta x = s2^{-level}\]
• A cell is represented by integer coordinates given its location. \[center = \Delta x (indices + 0.5)\]
• The adapted mesh is generally graded.

Identify the different types of cells

Identify the different types of cells

Identify the different types of cells

Identify the different types of cells

Algebra of sets

⋃

=

We need to add a set algebra that allows us to switch from one to the other easily.

Algebra of sets

\

=

Algebra of sets

The search of an admissible set is recursive. The algorithm starts from the last dimension (y in 2d, z in 3d,…).

The available operators in samurai are for now

• the intersection of sets,
• the union of sets,
• the difference between two sets,
• the translation of a set,
• the extension of a set.

The search for an admissible set is completely recursive. We’ll start by searching in the largest dimension and if we don’t find anything we’ll stop.

If an admissible interval is found, we look at it in the lower dimension, which means we manipulate very few intervals because we’re only looking at a subset.

Algebra of sets: Usage to MRA

initial mesh

Let’s take a closer look at how we use intervals and set algebra in mutliresolution. You could do the same with AMR. The final sets just won’t be exactly the same.

We start with this mesh.

Algebra of sets: Usage to MRA

initial mesh

we can classify cells by levels.

Algebra of sets: Usage to MRA

ghost cells used by the numerical scheme

As we said earlier, we’ll probably need ghosts to integrate our numerical scheme in space. We’ll therefore need to add cells on either side of the leaves, depending on the maximum size of the stencil.

Algebra of sets: Usage to MRA

cells needed to compute the details (all tree)

We have shown that for multiresolution, we need to calculate details. In the original version, you have the whole tree and are able to calculate details on all levels. It’s a completely recursive algorithm.

However, if we do this, we once again lose a compressed view of the mesh, adapted with more cells than necessary. We therefore propose to make this algorithm iterative and modify the mesh until it stops moving.

Algebra of sets: Usage to MRA

cells needed to compute the details (iterative way)

So we start with each leaf and add the cells needed to make the prediction for that leaf at the bottom level. This obviously depends on the order of the prediction, but we’ll assume here that we have a stencil of 1.

The idea is then to calculate the details, adapt the mesh accordingly and iterate by calculating the details again on the new mesh until it stops moving.

Algebra of sets: Usage to MRA

cells needed to compute the details (iterative way)

But we’ve got a problem: let’s imagine that these two cells are to be coarsened.

Algebra of sets: Usage to MRA

cells needed to compute the details (iterative way)

You’ll have this new mesh. You calculate the details and realize that you’ll have to refine this new cell after all. This happens because we’re in an iterative process, so we don’t have a global view of the details, but rather a local one.

We could say that it’s no big deal, but in the end we’ve put the cells back the way we started, except that we’ve added them by making a prediction. So we’ve degraded the initial values. We really don’t want to do that.

How can we correct this problem?

Algebra of sets: Usage to MRA

cells needed to compute the details (iterative way)

All we need to do is calculate the details on two successive levels to avoid this problem, so we add these cells.

Algebra of sets: Usage to MRA

cells needed to compute the details (iterative way)

In this example, a problem persists. The cell at the very bottom needs to be calculated, and since it has fine cells above it, it needs to be computed using the projection operator. However, a cell above is missing. So we need to add it.

Algebra of sets: Usage to MRA

cells needed to compute the details (iterative way)

You may be thinking that in the end we’ve added as many cells as the whole tree, and you’d be right. But these different mesh layers are rather complicated to illustrate in 1d because the mesh is quite small. This is no longer for N-D case, and what we’ve shown you in this simplistic example makes perfect sense.

Meshes description with samurai

Given a mesh with all the leaves called cells

• ghosts for the numerical scheme

samurai::for_each_interval(mesh[cells], [&](std::size_t level, auto& i)
{
    mesh[cells_and_ghosts][level].add_interval({i.start - stencil_size, i.end + stencil_size});
})

Meshes description with samurai

Given a mesh with all the leaves called cells

• ghosts for the detail computation

for(std::size_t level = min_level + 1; level <= max_level; ++level)
{
    auto subset = difference(mesh[cells_and_ghosts][level], union()[level]);

    subset([&](std::size_t level, auto& i)
    {
        auto i_above = i >> 1;
        mesh[pred][level - 1].add_interval({i_above.start - prediction_size, i_above.end + prediction_size});
    })
}

for(std::size_t level = min_level + 2; level <= max_level; ++level)
{
    samurai::for_each_interval(mesh[cells][level], [&](std::size_t level, auto& i)
    {
        auto i_above = i >> 2;
        mesh[pred][level - 2].add_interval({i_above.start - prediction_size, i_above.end + prediction_size});
    })
}

Meshes description with samurai

Given a mesh with all the leaves called cells

• ghosts where we apply the projection

for(std::size_t level = min_level; level < max_level; ++level)
{
    auto subset = intersection(mesh[reference][level], union()[level]).on(level+1);
    subset([&](std::size_t level, auto& i)
    {
        mesh[reference][level + 1].add_interval(i);
        mesh[proj][level].add_interval(i >> 1);
    })
}

Compression rates

Compression rates

Level	Num. of cells	p4est	samurai (leaves)	samurai (all)	ratio
\(9\)	66379	2.57 Mb	33.68 Kb	121 Kb	21.24
\(10\)	263767	10.25 Mb	66.64 Kb	236.8 Kb	43.28
\(11\)	1051747	40.96 Mb	132.36 Kb	467.24 Kb	87.66
\(12\)	4200559	163.75 Mb	263.6 Kb	927 Kb	176.64
\(13\)	16789627	654.86 Mb	525.9 Kb	1.85 Mb	353.98
\(14\)	67133575	2.61 Gb	1.05 Mb	3.68 Mb	709.24

We can see here that the mesh is very compressed compared to a version using a tree structure.

The total size of the mesh is multiplied by just under 4, as we build different meshes to handle ghost cells and prediction cells.

What’s important to note is that in a parallel context, it doesn’t cost much for each process to have the meshes of its neighboring domains.

Other features

• Loop algorithms over the levels and the cells
• Simplified access operator
• Helper classes to construct complex meshes
• Helper classes to construct schemes for explicit and implicit usage
• Helper classes to construct N-D operators and expressions using xtensor
• HDF5 support

Numerical Analysis and Modified Equations

Linear scalar transport equation

In this work, we are concerned with the numerical solution of the Cauchy problem associated with the linear scalar conservation law

\[ \partial_t u(t, x)+V \partial_x u(t, x)=0, \quad(t, x) \in \mathbb{R}^{+} \times \mathbb{R} \]

where \(V\) is the transport velocity, taken \(V>0\) without loss of generality. we consider 1 d problems. The extension to \(2\mathrm{~d}\) / \(3\mathrm{~d}\) problems is straightforward and usually done by tensorization [Bellotti2022] and yields analogous conclusions.

The discrete volumes are \[ C_{\ell, k}:=\left[2^{-\ell} k, 2^{-\ell}(k+1)\right], \quad k \in \{ 0,2^{\ell}-1 \}, \] for any \(\ell \in \{ \underline{\ell}, \bar{\ell} \}\). The measure of each cell at level \(\ell\) is \(\Delta x_{\ell}:=2^{-\ell}\) and we shall indicate \(\Delta x:=\Delta x_{\bar{\ell}}\). The cell centers are \(x_{\ell, k}:=\) \(2^{-\ell}(k+1 / 2)\). Finally, we shall indicate \(\Delta \ell:=\bar{\ell}-\ell\), hence \(\Delta x_{\ell}=2^{\Delta \ell} \Delta x\).

Finite Volume scheme

Finite Volume scheme at the finest level of resolution \(\bar{\ell}\) for any cell of indices \(\bar{k} \in \{ 0,2^{\bar{\ell}}-1 \}\). Explicit schemes read:

\[ \mathrm{v}_{\bar{\ell}, \bar{k}}^{n+1}=\mathrm{v}_{\bar{\ell}, \bar{k}}^n-\frac{\Delta t}{\Delta x}\left(\Phi\left(\mathrm{v}_{\bar{\ell}, \bar{k}+1 / 2}^n\right)-\Phi\left(\mathrm{v}_{\bar{\ell}, \bar{k}-1 / 2}^n\right)\right) \]

where we utilize the same linear numerical flux for the left and the right flux (conservativity) \[ \Phi\left(\mathbf{v}_{\bar{\ell}, \bar{k}-1 / 2}\right):=V \sum_{\alpha=\underline{\alpha}}^{\bar{\alpha}} \phi_\alpha \mathbf{v}_{\bar{\ell}, \bar{k}+\alpha}, \quad \Phi\left(\mathbf{v}_{\bar{\ell}, \bar{k}+1 / 2}\right):=V \sum_{\alpha=\underline{\alpha}}^{\bar{\alpha}} \phi_\alpha v_{\bar{\ell}, \bar{k}+1+\alpha} \]

Modified equations

[Carpentier et al 97] or Cauchy-Kowalewski procedure [Harten et al 87]

\[ \partial_t u\left(t^n, x_{\bar{\ell}, \bar{k}}\right)+V \partial_x u\left(t^n, x_{\bar{\ell}, \bar{k}}\right)=\sum_{h=2}^{+\infty} \Delta x^{h-1} \sigma_h \partial_x^h u\left(t^n, x_{\bar{\ell}, \bar{k}}\right) \]

• Upwind scheme \[ \partial_t u+V \partial_x u=\frac{\Delta x V}{2}(1-\lambda V) \partial_{x x} u+O\left(\Delta x^2\right) \]
• Lax-Wendroff scheme \[ \partial_t u+V \partial_x u=-\frac{\Delta x^2 V}{6}\left(1-\lambda^2 V^2\right) \partial_x^3 u+O\left(\Delta x^3\right) \]
• OSMP-3 scheme \[ \partial_t u+V \partial_x u=\frac{\Delta x^3 V}{24}\left(-\lambda^3 V^2+2 \lambda^2 V^2+\lambda V-2\right) \partial_x^4 u+O\left(\Delta x^4\right) \]

How to include MRA I

We introduce the reconstruction operator \(\hat{s}\) instead of \(s\) on the cells \(\left(\bar{\ell}, 2^{\Delta \ell} k+\delta\right)\) for any \(\delta \in \mathbb{Z}\) at the finest level

• \(\hat{s}=s\) : exact local flux reconstruction [Cohen et al. 2003].
• \(\hat{s}=0\) but \(s>0\), direct evaluation or naive evaluation [Hovhannisyan et al 2010].

\[ \mathbf{w}_{\bar{\ell}, \bar{k}}^{n+1}=\mathbf{w}_{\bar{\ell}, \bar{k}}^n-\frac{\Delta t}{\Delta x}\left(\Phi\left(\hat{\hat{\mathbf{w}}}_{\bar{\ell}, \bar{k}+1 / 2}^n\right)-\Phi\left(\hat{\hat{\mathbf{w}}}_{\bar{\ell}, \bar{k}-1 / 2}^n\right)\right) \] \[ \Phi\left(\hat{\hat{\mathbf{w}}}_{\bar{\ell}, \bar{k}-1 / 2}\right):=V \sum_{\alpha=\underline{\alpha}}^{\bar{\alpha}} \phi_\alpha \hat{\hat{\mathbf{w}}}_{\bar{\ell}, \bar{k}+\alpha} \]

How to include MRA II

Let now \((\ell, k) \in S\left(\tilde{\Lambda}^{n+1}\right)\), taking the projection yields the multiresolution scheme

\[ \mathbf{w}_{\ell, k}^{n+1}=\mathbf{w}_{\ell, k}^n-\frac{\Delta t}{\Delta x_{\ell}}\left(\Phi\left(\hat{\hat{\mathbf{w}}}_{\bar{\ell}, 2^{\Delta \ell}(k+1)+1 / 2}^n\right)-\Phi\left(\hat{\hat{\mathbf{w}}}_{\bar{\ell}, 2^{\Delta \ell} k-1 / 2}^n\right)\right) \]

\[ \Phi\left(\hat{\hat{\mathbf{w}}}_{\bar{\ell}, 2^{\Delta \ell} k-1 / 2}\right):=V \sum_{\alpha=\underline{\alpha}}^{\bar{\alpha}} \phi_\alpha \hat{\hat{\mathbf{w}}}_{\bar{\ell}, 2^{\Delta \ell} k+\alpha} \]

Some information is loss because of the averaging procedure: two different schemes we can consider for the computation of the modified equations.

Theorem

The local truncation error of the reference Finite Volume scheme and the one of the adaptive Finite Volume scheme are the same up to order \(2\hat{s}+1\) included.

Modified equations including MRA

This result establishes at which order the modified equations of the reference scheme are perturbed by the introduction of the adaptive scheme. However, it does not characterize the terms in the modified equations above order \(2\hat{s}+1\) in \(\Delta x\) (symbolic computations).

• Upwind scheme \[ \begin{array}{lr} \partial_t u+V \partial_x u=\frac{\Delta x V}{2}\left(2^{\Delta \ell}-\lambda V\right) \partial_{x x} u+O\left(\Delta x^2\right), & \text { for } \hat{s}=0 \\ \partial_t u+V \partial_x u=\frac{\Delta x V}{2}(1-\lambda V) \partial_{x x} u-\frac{\Delta x^2 V}{6}\left(1-\lambda^2 V^2\right) \partial_x^3 u+& \\ \ \ \ \ \ \ \ \ \ \ \frac{\Delta x^3 V}{24}\left(-3 \Delta \ell 2^{2 \Delta \ell}+2^{2 \Delta \ell}-\lambda^3 V^3\right) \partial_x^4 u+O\left(\Delta x^4\right), & \text { for } \hat{s}=1 \end{array} \]
• Lax-Wendroff scheme \[ \begin{array}{lr} \partial_t u+V \partial_x u=\frac{\Delta x \lambda V^2}{2}\left(2^{\Delta \ell}-1\right) \partial_{x x} u+O\left(\Delta x^2\right), & \text { for } \hat{s}=0 \\ \partial_t u+V \partial_x u=-\frac{\Delta x^2 V}{6}\left(1-\lambda^2 V^2\right) \partial_x^3 u+&\\ \ \ \ \ \ \ \ \ \ \ \frac{\Delta x^3 \lambda V^2}{24}\left(-3 \Delta \ell 2^{2 \Delta \ell}+2^{2 \Delta \ell}-\lambda^2 V^2\right) \partial_x^4 u+O\left(\Delta x^4\right), & \text { for } \hat{s}=1 \end{array} \]

Theoretical results on the global error

Theorem 2

Assume that

• The reference scheme satisfies the restricted stability condition \(\|E\| \leq 1\)
• The Harten-like scheme satisfies the restricted stability condition \(\left\|\bar{E}_{\Lambda}\right\| \leq 1\) for any \(\Lambda\).

Then, for smooth solution, in the limit \(\Delta x \rightarrow 0\) (i.e. \(\bar{\ell} \rightarrow+\infty\) ) and for \(\Delta \underline{\ell}=\bar{\ell}-\underline{\ell}\) kept fixed, we have the error estimate

\[ \left\|\mathbf{v}_{\bar{\ell}}^n-\mathbf{w}_{\bar{\ell}}^n\right\| \leq C_{t r} t^n \Delta x^{2 \hat{s}+1}+C_{m r} \frac{t^n}{\lambda \Delta x} \epsilon \]

where \(C_{t r}=C_{t r}\left(\bar{\ell}-\underline{\ell},\left(\phi_\alpha\right)_\alpha, \lambda, \hat{s}, V\right)\) and \(C_{m r}=C_{m r}\left(\bar{\ell}-\underline{\ell},\left(\phi_\alpha\right)_\alpha, \lambda, \hat{s}, s, V\right)\). \[ \left\|\mathbf{u}_{\bar{\ell}}^n-\mathbf{w}_{\bar{\ell}}^n\right\| \leq C_{r e f} t^n \Delta x^\theta+C_{t r} t^n \Delta x^{2 \hat{s}+1}+C_{m r} \frac{t^n}{\lambda \Delta x} \epsilon \]

Let us start by discussing the assumption that we have placed in the statement of Theorem 2: - The restricted stability condition \(\|E\| \leq 1\) could be replaced by a milder condition \(\|E\| \leq 1+C \Delta t\) for some constant \(C \geq 0\), see (69) in [7] and (A2) in [15]. This would not change the result. The technical assumption \(\left\|\bar{E}_{\Lambda}\right\| \leq 1\) is harder to relax and also difficult to check in practice. - The fact of considering smooth solutions comes from the fact that we want to apply the analysis of the modified equations to obtain the convergence rates, in the spirit of the Lax theorem [1]. For the same reason, we take \(\Delta x \rightarrow 0\) (or \(\bar{\ell} \rightarrow+\infty\) ). - The distance between maximum and minimum level \(\Delta \underline{\ell}=\bar{\ell}-\underline{\ell}\) has to be fixed, because otherwise the constant \(C_{\text {tr }}\) potentially explodes and dominates \(\Delta x^{2 \hat{s}+1}\) when \(\Delta x \rightarrow 0\). This would prevent us from comparing orders. Moreover, this is also reasonable from the standpoint of actual computations, where we refine the mesh to achieve convergence (or nearly so) keeping the number of different available grid levels fixed. Still, we shall also perform numerical demonstration without fixing \(\Delta \underline{\ell}=\bar{\ell}-\underline{\ell}\) to show that the modified equations that we have previously developed provide important information on the behavior of \(\left\|\mathbf{v}_{\bar{\ell}}^n-\mathbf{w}_{\bar{\ell}}^n\right\|\).

Comments on theorem

• The error estimate contains three contributions: the discretization error of the reference scheme, the perturbation error between the reference and the adaptive scheme, and the thresholding error coming from the multiresolution
• The constant \(C_{\mathrm{tr}}\) generally grows exponentially with \(\bar{\ell}-\underline{\ell}\), sometimes also involving linear terms, i.e. \(\hat{s}=1\). We have the following cases:
  • \(\theta<2 \hat{s}+1\). The error of the reference scheme dominates the perturbation introduced by the adaptive scheme \(\left\|\mathbf{u}_{\bar{\ell}}^N-\mathbf{w}_{\bar{\ell}}^N\right\| \leq C_{\mathrm{ref}} T \Delta x^\theta+C_{\mathrm{mr}} \frac{T}{\lambda \Delta x} \epsilon\). A thresholding error of the same order as the reference error \(\epsilon \sim \Delta x^{\theta+1}\).
  • \(\theta=2 \hat{s}+1\). The error of the reference scheme and the perturbation order are comparable (first example!) We have \(\left\|\mathbf{u}_{\bar{\ell}}^N-\mathbf{w}_{\bar{\ell}}^N\right\| \leq\left(C_{\mathrm{ref}}+C_{\mathrm{tr}}\right) T \Delta x^\theta+C_{\mathrm{mr}} \frac{T}{\lambda \Delta x} \epsilon\).
  • \(\theta>2 \hat{s}+1\). The perturbation introduced by the adaptive scheme dominates the error of the reference scheme. Therefore, multiresolution introduces a large perturbation that yields a different convergence rate. We have \(\left\|\mathbf{u}_{\bar{\ell}}^N-\mathbf{w}_{\bar{\ell}}^N\right\| \leq C_{\mathrm{tr}} T \Delta x^{2 \hat{s}+1}+C_{\mathrm{mr}} \frac{T}{\lambda \Delta x} \epsilon\), thus \(\epsilon \sim \Delta x^{2 \hat{s}+2}\) (AMR !)

Assume that for the choice of \(\bar{\ell}-\underline{\ell}\) at hand, we have \(C_{\mathrm{tr}} \sim C_{\mathrm{ref}}\), then w

Burgers results for upwind scheme

How to compute fluxes at the finest level

As we have seen, calculating the flux at the finest level is crucial to obtaining good error estimates. But how does this work in practice?

Let’s take a simple example.

How to compute fluxes at the finest level

Let’s imagine that our finest level is 4 levels away from where we are.

How to compute fluxes at the finest level

Now let’s imagine that we want to calculate the flow using this small cell on the left. We use an operator with s = 1, so we need 3 cells on the level below to calculate it. These 3 cells need 4 cells on the level below to be calculated. We can repeat the process indefinitely and we’ll still need 4 cells on the bottom level.

We can determine the coefficients of our linear application to go from our level to any cell on a higher level. There’s no need for an intermediate level.

How to compute fluxes at the finest level

We can do the same for the right-hand cell. We can see that the cells used are not exactly the same.

Now, if we wanted to reconstruct all the cells of the finest level, we’d only need these 5 cells, with different coefficients depending on the cell we’re trying to reconstruct.

MRA \(\epsilon = 1e-2\), \(\underline{\ell} = 1\), \(\bar{\ell} = 12\) without portion

MRA \(\epsilon = 1e-3\), \(\underline{\ell} = 1\), \(\bar{\ell} = 12\) without portion

MRA \(\epsilon = 1e-4\), \(\underline{\ell} = 1\), \(\bar{\ell} = 12\) without portion

MRA \(\epsilon = 1e-2\), \(\underline{\ell} = 1\), \(\bar{\ell} = 12\) with portion

MRA \(\epsilon = 1e-3\), \(\underline{\ell} = 1\), \(\bar{\ell} = 12\) with portion

MRA \(\epsilon = 1e-4\), \(\underline{\ell} = 1\), \(\bar{\ell} = 12\) with portion

AMR with the modulus of the gradient \(<1e-2\) as refinement criteria without portion

AMR with the modulus of the gradient \(<1e-2\) as refinement criteria with portion

AMR with the modulus of the second derivative \(<10\) as refinement criteria without portion

AMR with the modulus of the second derivative \(<100\) as refinement criteria without portion

AMR with the modulus of the second derivative \(<10\) as refinement criteria with portion

AMR with the modulus of the second derivative \(<100\) as refinement criteria with portion

References

Courses

• Master “AMS” (Analyse Modélisation Simulation) Course MSX2TA Méthodes Numériques Avancées et Calcul Haute performance pour la simulation de phénomènes complexes (L. Séries, M. M.)

Theses

• T. Bellotti, Numerical analysis of lattice Boltzmann schemes : from fundamental issues to efficient and accurate adaptive methods, PhD Thesis, Institut Polytechnique de Paris, EDMH (2023) https://theses.hal.science/tel-04266822

Publications

• T. Bellotti , L. G. , B. Graille , M. M.. High accuracy analysis of adaptive multiresolution-based lattice Boltzmann schemes via the equivalent equations, SMAI Journal of Computational Mathematics, Volume 8 (2022), pp. 161-199
• T. Bellotti, J. Massot, M. M., L. Séries, C. Tenaud, Modified equation and error analysis of adaptive multiresolution Finite Volume schemes, in preparation - to be submitted (2025)

A dual Splitting/IMEX strategy for PDEs

A dual Splitting/IMEX strategy for stiff PDEs

• A strategy has been designed (PhD M. Duarte) relying on time-adaptive operator splitting with dynamically adapted mesh (multiresolution) and error control:
  • optimal computational cost and parallelization properties (large splitting time steps)
• We aim at resolving stiff PDEs with samurai and ponio libraries with the same computation favorable properties:
  • local implicitation of the source term
  • explicit diffusion integration without von Neumann stability limit (ROCK)
  • high-order in space and time integration of convection, including shocks
  • strong acceleration through adaptation in space and time
• Stumbling block: what to do when reaction and diffusion coupled at smallest time scale (complex chem. - ignition)

A dual Splitting/IMEX strategy for stiff PDEs

• Keep reaction and diffusion coupled (IMEX)!
• No splitting time step limitation due to the coupling at small scale (complex chemistry and detailed transport in flames for example)
• Inspired from PIROCK strategy of G. Vilmart and A. Abdulle (ponio library : J. Massot)
• Adaptation in space and time with error control
• Project with CEA on hydrogen risk - DNS of hydrogen combustion with detailed transport and complex chemistry - Collaboration P.A. Masset and L. Lecointre

A dual Splitting/IMEX strategy for stiff PDEs

Ignition and Combustion (Verwer & Hundsdorfer)

Specific for RKC methods and ROCK

\[ u_t = d \Delta u + \frac{R}{\alpha \delta} (1 + \alpha - u) e^{\delta( 1 - 1/u )} \quad \text{in } [0,1]\times[0,1] \]

A dual Splitting/IMEX strategy for stiff PDEs

Ignition and Combustion (Verwer & Hundsdorfer)

Specific for RKC methods and ROCK

\[ u_t = d \Delta u + \frac{R}{\alpha \delta} (1 + \alpha - u) e^{\delta( 1 - 1/u )} \quad \text{in } [0,1]\times[0,1] \]

	PIROCK	ROCK2	ROCK4
Tolerance used	1e-5	1e-6	1e-5
L2 norm of error	1.84e-4	2.27e-4	2.69e-4
Elapsed time (s)	7.97	1.90	0.42
Nb. time steps	889	157	57
Nb. time steps rejected	0	0	0
Nb. function evaluation	-	6367	1734
Nb. function evaluation (diffusion)	8304	-	-
Nb. function evaluation (reaction)	360067655	-	-

between \(t=0.29\) and \(t=0.32\).

IMEX reasonable overhead compared to fully explicit methods (ROCK2 - ROCK4) even in this very difficult configuration

• Error control
• Good BC treatment
• No splitting errors

A dual Splitting/IMEX strategy for stiff PDEs

Belousov-Zhabotinsky (very stiff source - 3 eq) \[ \left\{ \begin{aligned} \partial_t a - D_a \, \Delta a &= \frac{1}{\mu} ( -qa - ab + fc) \\ \partial_t b - D_b \, \Delta b &= \frac{1}{\varepsilon} ( qa - ab + b\,(1-b)) \\ \partial_t c - D_c \, \Delta c &= b - c \end{aligned} \right. \]

• Error to the reference quasi-exact solution is second order in time but not of the same origin (splitting error vs. IMEX error) - but still error control
• Larger time step can be taken with IMEX while keeping a proper solution (no disastrous splitting errors - wrong wave speed)
• When optimal large splitting time step is taken, IMEX as efficient as splitting, whereas it is advantageous for smaller time steps as well as larger time steps
• No boundary condition problems
• Same computational good properties

A dual Splitting/IMEX strategy for stiff PDEs

Belousov-Zhabotinsky (very stiff source - 3 eq) \[ \left\{ \begin{aligned} \partial_t a - D_a \, \Delta a &= \frac{1}{\mu} ( -qa - ab + fc) \\ \partial_t b - D_b \, \Delta b &= \frac{1}{\varepsilon} ( qa - ab + b\,(1-b)) \\ \partial_t c - D_c \, \Delta c &= b - c \end{aligned} \right. \]

• Error to the reference quasi-exact solution is second order in time but not of the same origin (splitting error vs. IMEX error) - but still error control
• Larger time step can be taken with IMEX while keeping a proper solution (no disastrous splitting errors - wrong wave speed)
• When optimal large splitting time step is taken, IMEX as efficient as splitting, whereas it is advantageous for smaller time steps as well as larger time steps
• No boundary condition problems
• Same computational good properties

	PIROCK (tol=5e-3)	STRANG (dt=2/128)
L2 norm of error (variable b)	3.31e-2	1.14e-1
L2 norm of error (variable c)	2.92e-3	1.37e-2
Elapsed time (s)	6.91	6.70
Nb. time steps.	624	128

 

	ImEx (tol=1e-3)	STRANG (dt=2/512)
L2 norm of error (variable b)	4.94e-3	1.55e-3
L2 norm of error (variable c)	4.31e-4	1.41e-4
Elapsed time (s)	13.33	17.20
Nb. time steps.	1528	512

 

	ImEx (tol=1e-4)	STRANG (dt=2/16384)
L2 norm of error (variable b)	7.72e-6	1.81e-5
L2 norm of error (variable c)	6.97e-7	1.57e-6
Elapsed time (s)	165.12	427.81
Nb. time steps.	31728	16384

A dual Splitting/IMEX strategy for stiff PDEs

Belousov-Zhabotinsky (very stiff source - 3 eq) \[ \left\{ \begin{aligned} \partial_t a - D_a \, \Delta a &= \frac{1}{\mu} ( -qa - ab + fc) \\ \partial_t b - D_b \, \Delta b &= \frac{1}{\varepsilon} ( qa - ab + b\,(1-b)) \\ \partial_t c - D_c \, \Delta c &= b - c \end{aligned} \right. \]

• Error to the reference quasi-exact solution is second order in time but not of the same origin (splitting error vs. IMEX error) - but still error control
• Larger time step can be taken with IMEX while keeping a proper solution (no disastrous splitting errors - wrong wave speed)
• When optimal large splitting time step is taken, IMEX as efficient as splitting, whereas it is advantageous for smaller time steps as well as larger time steps
• No boundary condition problems
• Same computational good properties

Simulation with ponio / samurai

\(\epsilon = 1e-3\), \(\underline{\ell} = 2\), \(\bar{\ell} = 10\)

A dual Splitting/IMEX strategy for stiff PDEs

Ignition of a diffusion flame at 4 instants (cold fuel, hot oxydizer - [Candel & Thevenin 1995])

\[ \left\{ \begin{array}{l} \partial _{t_\star} Z + v_{x_\star} \partial _{x_\star} Z + v_{y_\star} \partial _{y_\star} Z - \left(\partial ^2_{x_\star} Z + \partial ^2_{y_\star} Z \right)= 0, \\[2.5ex] \partial _{t_\star} \theta + v_{x_\star} \partial _{x_\star} \theta + v_{y_\star} \partial _{y_\star} \theta - \left(\partial ^2_{x_\star} \theta + \partial ^2_{y_\star} \theta \right) = F(Z,\theta), \end{array} \right. \]

with

\[ F(Z,\theta)= D_a \, \phi \chi Y_{O,0} \left[ \frac{1-Z}{\phi \tau} + \frac{1}{\chi}(Z-\theta) \right] \left[ Z + \frac{\tau}{\chi}(Z-\theta) \right] \mathrm{e}^{\left(- \tau_a/(1+\tau \theta)\right)} \]

Quasi-exact solution calculated with Strang’s method using :

• a splitting time step of \(10^{-8}\), and
• the RADAUIIA method for the reaction, with absolute and relative tolerances of \(10^{-12}\).
• the ROCK4 method for diffusion, with absolute and relative tolerances of \(10^{-12}\).
• the RK3 method associated with a WENO5 spatial discretisation for convection

A dual Splitting/IMEX strategy for stiff PDEs

Evolution of the variable \(T = \theta (T_{F,0}-T_{O,0}) + T_{O,0}\) at \(t=0\), \(5\times10^{-5}\), \(1\times10^{-5}\), \(1.5\times10^{-4}\)

A dual Splitting/IMEX strategy for stiff PDEs

Strang (with RADAUIIA and ROCK4 tolerance at \(10^{-12}\)) :

	\(dt=1 \times 10^{-5}\)	\(dt=1 \times 10^{-6}\)	\(dt=1 \times 10^{-7}\)
\(|E_z|\)	\(1.55513306\times 10^{-4}\)	\(2.81342860\times10^{-6}\)	\(2.88674760\times10^{-8}\)
\(|E_{\theta}|\)	\(2.69304858\times 10^{-3}\)	\(3.74837111\times10^{-5}\)	\(3.82828537\times10^{-7}\)
Temps RD (s)	\(25+213\)	\(207+287\)	\(2060+520\)
Temps C (s)	\(13316\)	\(13471\)	\(14921\)

Dual Splitting/ImEx (ImEx tolerance set at \(10^{-4}\) - should be decreased) :

	\(dt=1 \times 10^{-5}\)	\(dt=1 \times 10^{-6}\)	\(dt=1 \times 10^{-7}\)
\(|E_z|\)	\(1.57299998\times 10^{-5}\)	\(2.71930151\times10^{-6}\)	\(2.87698144\times10^{-8}\)
\(|E_{\theta}|\)	\(1.93940574\times 10^{-3}\)	\(4.15846000\times10^{-6}\)	\(5.35535514\times10^{-6}\)
Temps RD (s)	\(28\)	\(123\)	\(1202\)
Temps C (s)	\(11221\)	\(11361\)	\(12575\)

A dual Splitting/IMEX strategy for stiff PDEs

OSMP scheme leads to roughly an order of magnitude acceleration compared to WENO5/RK3

Basis for the CEA collaboration - hydrogen risk (Full compressible Fourier-Navier-Stokes)

References

Theses

• M. Duarte, Adaptive numerical methods in time and space for the simulation of multi-scale reaction fronts, PhD Thesis, École Centrale Paris, (2011) https://theses.hal.science/tel-00667857
• L. Lecointre, Hydrogen flame acceleration in non-uniform mixtures, PhD Thesis, Université Paris Saclay, (2022) - https://theses.hal.science/tel-03879925

Publications

• M. Duarte , M. M., S. Descombes, C. Tenaud, T. Dumont, New resolution strategy for multi-scale reaction waves using time operator splitting, space adaptive multiresolution and dedicated high order implicit/explicit time integrators, SIAM SISC, vol. 34, No. 1 (2012) pp.76-104
• M. Duarte, S. Descombes, C. Tenaud, S. Candel, M. M., Time-space adaptive numerical methods for the simulation of combustion fronts, Combustion and Flame, vol. 160, No. 6 (2013) pp.1083-1101
• J. Massot, L. Séries, L. G., P. Matalon, C. Tenaud, M. M., A splitting/ImEx strategy for stiff PDEs with time adaptation and error control, invited contribution to Comptes Rendus Mécanique (2025)

Software

• J. Massot, M. M., L. Series, ponio (2024) https://hal.science/hal-04710549v1
• T. Bellotti, L. G., M. M., P. Matalon, samurai (2023) https://hal.science/hal-04545389v1

Roadmap

Here’s the samurai roadmap for the coming months and probably years.

We now have the possibility of changing the containers in which our fields are stored. We did this so that we could start using Kokkos. So we’d be very interested in working with the people at Cexa on this.

We have received several funding packages that have enabled us to hire two research engineers.

As part of numpex, we also have an engineer position to work on the I/O part of AMR methods in collaboration with Dyablo developers.

Scientific Collaborations

• Lattice Boltzmann methods and multiresolution - Thomas Bellotti (EM2C/CNRS/CS) and Benjamin Graille (LMO/Université Paris-Saclay)
• Plasma discharges and electric propulsion - Alejandro Alvarez (LPP/École polytechnique) and Louis Reboul (ONERA)
• DNS of lithium-ion batteries based on high-resolution 3D images of porous electrode microstructures - Ali Asad (TotalEnergies) and Laurent François (ONERA)
• Sharp interface method for low Mach two-phase flows - Nicolas Grenier (LISN/Université Paris-Saclay) and Christian Tenaud (EM2C/CNRS/CS)
• Low-Mach reactive flows - Christian Tenaud (EM2C/CNRS/CS)
• Interfacial flow simulation - Giuseppe Orlando and Ward Haegeman (CMAP/Ecole polytechnique), Samuel Kokh (CEA/MdlS), Joël Dupays and Clément Le Touze (ONERA), Marica Pelanti (ENSTA/IP Paris), Khaled Saleh (Aix-Marseille Université), Jean-Marc Hérard (EDF)
• Mathematical modeling and simulation of non-equilibrium plasmas for the prediction of electric propulsion - Teddy Pichard and Zoubaïr Tazakkati (CMAP/École polytechnique)
• Simulation analysis on the Hydrogen risk - Luc Lecointre, Pierre-Alexandre Masset, Etienne Studer (CEA) and Christian Tenaud (EM2C/CNRS/CS)