samurai/ponio: two complementary software libraries
for the efficient implementation of adaptive space/time
numerical schemes

Sébastien Dubois, Loïc Gouarin, Alexandre Hoffmann, Josselin Massot, Marc

Massot, Pierre Matalon, Laurent Series

► To cite this version:

Sébastien Dubois, Loïc Gouarin, Alexandre Hoffmann, Josselin Massot, Marc Massot, et al. · samurai/ponio: two complementary software libraries for the efficient implementation of adaptive space/time numerical schemes. Modeling Summer Visit, NASA Ames Research Center, Jul 2025, Mountain View, CA, United States. [hal-05503161](https://hal.science/hal-05503161v1)

HAL Id: hal-05503161

https://hal.science/hal-05503161v1

Submitted on 10 Feb 2026

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# samurai/ponio: two complementary software libraries for the efficient implementation of adaptive space/time numerical schemes

Sébastien Dubois<fnref n="1" />

Loïc Gouarin<fnref n="1" />

Alexandre Hoffmann<fnref n="1" />

Josselin Massot<fnref n="1" />

Marc Massot<fnref n="1" />

Pierre Matalon<fnref n="1" />

Laurent Series<fnref n="1" />

¹CMAP, CNRS, École polytechnique, Institut polytechnique de Paris

## Abstract

We present two complementary open-source C++ libraries for the numerical simulation of multiscale physical systems. **samurai** provides an efficient adaptive Cartesian meshing framework, while **ponio** offers a wide range of time-integration schemes for ordinary differential equations. Together, they enable the rapid implementation and testing of high-performance, adaptive space-time numerical methods for PDEs. Both libraries emphasize modularity, data locality, and modern C++ design. Their capabilities are illustrated through applications developed during the Modeling Summer Visit 2025, including two-phase flows, rarefied gases, and reactive flows for aerospace applications relevant to NASA.

**Keywords:** adaptive mesh, AMR, adaptive multiresolution, time-stepping schemes, HPC.

## 1 Introduction

Numerical simulation of complex physical phenomena requires the implementation of efficient numerical schemes for the discretization of both space and time. The optimal choice of numerical methods depends strongly on the problem under consideration, and the landscape of available approaches continues to evolve with the regular emergence of new algorithms. At the same time, limited computational resources and the demand for ever greater accuracy to resolve fine-scale physical features make adaptive methods essential for the efficient treatment of multiscale problems [9]. These observations highlight two complementary needs: (i) to quickly assess existing methods on complex problems such as the ones that can be found in aerospace engineering, and (ii) to easily prototype and test new numerical schemes. To address them, the pair of open-source libraries **samurai**<fnref n="1" />/ponio<fnref n="2" /> provides efficient implementations of reference numerical techniques for the fast solution of partial differential equations. **samurai** is a parallel Cartesian mesher designed for mesh adaptation, and **ponio** is an ODE solver (method of lines) that offers a wide range of time-stepping methods readily available to the user. Each library can be used independently,

and **samurai** alone enables the development of coupled space-time schemes with mesh adaptation. Both are implemented in modern C++, emphasizing flexibility and ease of integration. In the following, Section 2 presents the main features of **samurai**, Section 3 describes **ponio**, and Section 4 illustrates how both libraries have been applied in the context of the Modeling Summer Visit 2025 to a range of problems, including two-phase flows, rarefied gases, and combustion models.

## 2 The **samurai** library

Adaptive mesh refinement (AMR) [4] and multiresolution analysis (MRA) [12, 7] techniques are essential for resolving steep gradients, shocks, or localized features, while saving computational resources and memory trace in domain regions where the solution is smooth. **samurai** is a C++ library built around a data structure specifically designed to efficiently handle adaptive Cartesian meshes. The originality of **samurai**'s approach hinges on a special data structure enabling both AMR and MRA, improving data contiguity and facilitating the search for neighbours in stencil-based computations.

<footnote n="1">¹https://github.com/hpc-maths/samurai</footnote>

<footnote n="2">²https://github.com/hpc-maths/ponio</footnote>
## 2.1 The interval-based data structure

Most adaptive meshing frameworks rely on a quad-tree/octree structure, individual leaves representing a single cell (cell-based refinement) or a uniform block of cells (block and patch-based refinement). For instance, the meshing libraries `p4est` [6] and `AMREx` [18] use octrees. In contrast, the fundamental building block in `samurai` is the *interval*. An interval represents a sequence of identical, contiguous cells along a Cartesian direction, storing only the indices of the first and last cell and implicitly including all cells in between. This process, applies in the *x*-direction, enables a first step of compression. In a second step, the *y*-coordinates associated with these *x*-intervals are also grouped whenever they are contiguous, yielding a list of *y*-intervals. The association between each *y* coordinate and its corresponding list of *x*-intervals is organized in a manner reminiscent of the *compressed sparse row (CSR)* format widely used for sparse matrices. Additional axes (*z*, etc.) are handled recursively in the same way, which makes the approach dimension-independent. `samurai` is therefore not limited to 3D, but natively handles any dimension.

At each level of refinement, `samurai` maintains such a multidimensional interval structure, and the complete adaptive mesh is obtained as a hierarchy of these level-wise data structures. This approach drastically reduces memory usage and avoids the recursive pointer traversal typical of tree-based structures.

## 2.2 Field storage and stencil-based computations

When a data field is attached to the mesh, its memory layout is determined by the numbering of the underlying cells, since the stored values will follow the same order in a contiguous memory space. In classical octree-based approaches, cells are typically indexed using a space-filling curve [16] (e.g., Hilbert or Morton ordering). Such orderings improve spatial locality, which benefits cache efficiency, but they do not guarantee strict contiguity along Cartesian directions. This lack of axis-aligned contiguity can hinder vectorized computations, particularly for stencil-based numerical schemes.

In `samurai`, cells are numbered level by level and, within each level, interval by interval along the *x*-direction. As a result, data fields are stored contiguously along every *x*-interval, providing both memory alignment and cache-friendly access. This layout is especially advantageous for stencil computations: for instance, given a two-cell stencil to compute fluxes, once the stencil is anchored on two neighbouring cells (displayed either horizontally or vertically), it can slide along the *x*-direction to vectorize the computation of multiple fluxes, up until the end of the common interval where the stencil has been anchored.

Stencil-based schemes require, for each cell, quick access to its neighbors. A common strategy is to precompute and store these relationships in a correspondence table. While this allows fast lookup, it comes at the

cost of significant additional memory. In `samurai`, neighbor relations are computed on the fly, without storing explicit connectivity. The interval-based data structure provides a set of simple geometric operations—union, intersection, difference, translation, etc.—that make this possible. For example, translating an *x*-interval upward yields, in a single operation, the upper neighbors of all its cells. Only the index of the first cell in the translated interval needs to be retrieved from the mesh’s internal storage; the rest follow contiguously.

## 2.3 Mesh adaptation

The interval-based mesh hierarchy in `samurai` supports both AMR and MRA. Since the MSV projects presented here all rely on MRA, we give below a concise and deliberately simplified overview of this technique. Unlike AMR, which is a heuristic refinement strategy, the MRA approach is driven by wavelet compression. It does not require a user-defined refinement criterion. Instead, refinement and coarsening are guided by an intrinsic regularity test. The principle is as follows: projection and prediction operators are defined so that any local polynomial of degree up to a prescribed order is exactly preserved when projected onto coarser cells and then predicted back to the fine cells. This test is applied to the numerical solution in each cell. The discrepancy between the predicted value and the actual value defines the so-called *detail*. By design, the detail will be small if the solution is smooth, and will grow as large as the regularity of the solution is low. Should a cell be coarsened, the detail quantifies the information that cannot be recovered by the prediction operator. A user-defined tolerance specifies the acceptable loss of information, and cells are refined or coarsened to ensure that all details remain below this threshold. Compared to AMR, which is heuristics-driven, MRA presents the advantage of providing rigorous error control, directly tied to the actual smoothness of the solution.

## 2.4 Flux-based user interface

`samurai` provides dedicated tools for implementing finite volume schemes. In practice, the user only needs to define how a flux is computed at the interface between two cells, given a chosen stencil size around that interface. The implementation is fully dimension-independent: the same flux definition applies in 1D, 2D, or 3D. Mesh traversal is handled internally by `samurai`, which successively provides the stencil values needed for flux evaluation. As a result, the user can write schemes as if working on a uniform mesh, without having to account for level jumps. When such a jump occurs, `samurai` automatically supplies the appropriate stencil values: values predicted into ghost cells on the coarse side and values from real cells on the fine side of the interface. The following code illustrates the implementation of the Laplace operator (the flux is a gradient) using the classical second-order scheme with a two-cell stencil.
```cpp
[](
    FluxValue<cfg>& flux,
    const StencilData<cfg>& data,
    const StencilValues<cfg>& u
)
{
    const int L = 0; // left
    const int R = 1; // right
    double dx = data.cell_length;
    flux = (u[R] - u[L]) / dx;
};
```

Interfaces for both conservative and non-conservative fluxes are provided. In the conservative case, a single flux is computed: the contribution is added to the left cell and subtracted from the right one. In the non-conservative case (useful, e.g., in two-phase flow models), two distinct fluxes (incoming and outgoing) must be computed.

The application of MRA requires that all field values involved in flux computations be predicted at the finest level of resolution via the MRA prediction operator. To balance accuracy and cost, however, the prediction level can be chosen by the user. In either case, the flux function itself remains unchanged—the only difference lies in the stencil values provided by *samurai*.

# 3 The ponio library

ponio is a C++ library providing a collection of ready-to-use numerical integration schemes for ordinary differential equations (ODEs). It integrates first-order systems of the form

$$ \dot{u} = f(t, u), \quad (3.1) $$

where $t$ is the time variable, $u$ is the unknown state, and $f$ a temporal flux function.

ponio also allows users to extend its capabilities by defining custom Butcher tableaus, and provides Python tools for the automatic analysis of numerical methods.

## 3.1 Automatic method analysis and code generation

Each Runge-Kutta (RK) method is implemented directly from its Butcher tableau, which serves as the sole input. From this tableau, dedicated and optimized C++ code is automatically generated. User-defined tableaus can be imported in JSON format, triggering an analysis that computes the method's theoretical properties, such as its stability domain, order of convergence, and, when applicable, its stage order. The data of the 60+ Butcher tableaus already imported is available online<fnref n="3" />, covering a wide variety of classical and modern RK schemes. The same analysis tools are also employed for more advanced numerical schemes, such as extended-stability and IMEX methods.

## 3.2 Catalogue of methods

A wide range of Runge-Kutta (RK) methods are available, including explicit RK schemes and their associated

Lawson formulations [13], exponential RK methods, and diagonally implicit Runge-Kutta (DIRK) methods. Adaptive time-stepping is supported through embedded RK schemes, such as the Dormand-Prince family for explicit methods and ESDIRK(5,4) for implicit ones.

Several extended-stability methods are implemented: (i) with an arbitrary number of stages: the Runge-Kutta Chebyshev method RKC2 and the Runge-Kutta Legendre methods [15] of first and second orders (RKL1 and RKL2); (ii) with a dynamically adjusted number of stages: ROCK2 [2] and ROCK4 [1], available with both fixed and adaptive time steps. This kind of method is useful for diffusion problems or other problems with very large negative eigenvalues.

Additionally, ponio includes operator-splitting approaches such as Strang [8] and Lie splitting, as well as several (IMEX) schemes, offering a flexible framework for stiff and multiscale problems alike.

## 3.3 PIROCK: a dedicated scheme for reaction-diffusion problems

ponio also provides time-integration methods specifically designed for reaction-diffusion systems, among which the most prominent is the PIROCK [3] scheme.

Reaction-diffusion equations present major challenges for numerical integration, owing to the stiffness of the reaction operator and the spectrum of the diffusion operator. They are commonly treated using operator-splitting techniques, where each subproblem is integrated with a dedicated solver. Although effective, such approaches introduce a splitting error.

To overcome this limitation, the PIROCK method combines an extended-stability explicit scheme (ROCK2) for the diffusion operator—avoiding the solution of large linear systems—with a DIRK scheme for the reaction operator, which only requires solving small, independent, local systems. These two integrators are coupled within an IMEX framework, where the stabilization stages are performed independently of the reaction solve, followed by only four coupling stages.

Compared to the original PIROCK scheme, ponio extends the method with a flexible mechanism for error estimation. A naïve computation of the norm on an adaptive mesh tends to overestimate the error in the finest cells; ponio corrects this by weighting the error according to cell size and, when necessary, normalizing each field component when their magnitudes differ significantly.

## 3.4 A versatile user interface

Using C++ templates, ponio relies directly on the user's data structures to store the system state (i.e. the data type of $u$ in formula (3.1)). This can range from a single scalar (for scalar-valued problems) to any array or vector type, as long as it supports either basic linear operations or random-access. In the context of MSV 2025, *samurai* field types are typically employed as state containers. The example below illustrates how to define the flux function $f$ for the heat equation

<footnote n="3">³https://hpc-maths.github.io/ponio/</footnote>
tion $\dot{u} = \Delta u$ using *samurai*, and how to integrate it with *ponio's* RKC2 solver configured with 20 stages.

```cpp
// problem construction
auto u_ini = samurai::make_scalar_field
    <double>("u", mesh);
using Field = decltype(u_ini);
auto laplacian = samurai::
    make_laplacian_order2<Field>(
    auto f = [&](double t, auto&& u, auto& du)
    {
        du = laplacian(u);
    };

// ponio solver
auto method = ponio::runge_kutta::
    chebyshev::explicit_rkc2<20>(
    double T = 1;
    double dt = 1e-2;
auto sol_range = ponio::makeSolution_range(
    f, method, un_ini, {0., T}, dt);

// time loop
while (it sol->time < T)
{
    ++it sol;
    save(it sol->state);
}
```

# 4 Application to MSV projects

samurai and ponio have been used in MSV projects modeling two-phase flows [11], rarefied gases [5], reactive flows [10] and high-order code coupling [17]. Below, we give details on how the *samurai/ponio* frameworks have been leveraged for the first three projects.

## 4.1 Two-phase flow with scale separation

The two-phase flow model of [14] is extended to non-isothermal regimes by incorporating complete equations of state and accounting for energy exchanges during inter-scale mass transfer. This ensures total energy conservation and consistent entropy production in high-enthalpy flows such as ablation or fuel atomization.

Numerical discretizations of both models have been implemented in the *samurai* framework with Rusanov and HLLC fluxes. Both diffuse and sharp interface approaches are supported. Recent developments focus on improving numerical robustness and validating the fully compressible model through direct comparison between these two interface representations within the same *samurai* environment.

Note that the default MRA prediction operator is built as a constant linear combination of field values, without enforcing any exterior constraint. In particular, when predicting volumic fractions on the finest level, the resulting values may very well be out of admissible bounds (below 0 or above 1). To solve this issue, a custom-made prediction operator has been implemented.

## 4.2 Hybrid dense-rarefied flows

The simulation of hybrid dense-rarefied flows is addressed, where the transition between Navier-Stokes and kinetic regimes remains difficult to capture with traditional hybrid Navier-Stokes/direct simulation Monte Carlo methods. This work aims to develop a deterministic Eulerian fluid solver capable of spanning both continuous and rarefied regimes, reducing reliance on heuristic coupling between solvers.

The approach is based on the method of moments, which derives reduced fluid models from the Boltzmann equation by tracking macroscopic moments of the velocity distribution instead of individual molecular states. The main challenges include choosing suitable closure relations and ensuring consistent coupling between different moment hierarchies. A relaxation scheme is designed to smoothly connect these models across flow regions. The resulting 2D solver is implemented using the *samurai* library, leveraging its adaptive mesh refinement and error control capabilities to efficiently resolve local rarefaction gradients.

## 4.3 Multicomponent reactive flows

This work proposes a dual splitting/IMEX strategy for efficiently solving multicomponent reacting flow equations that exhibit strong stiffness and multiscale coupling in space and time. Traditional explicit schemes fail due to stability constraints, while fully implicit ones are too costly.

The approach combines the advantages of adaptive operator splitting with IMEX Runge-Kutta methods. The convection term is integrated using high-order monotonicity-preserving one-step schemes, whereas the reaction-diffusion subsystem is handled with the PIROCK method (cf. Section 3.3), offering stability and error control. This hybrid technique aims to retain the low computational cost of splitting methods while overcoming their limitations in strongly coupled regimes. It is tested on canonical reaction-diffusion systems (ignition and Belousov-Zhabotinsky problems) and will later be on combustion applications with increasing chemical complexity, including premixed and diffusion flames. The implementation relies on adaptive spatial and temporal discretization with *samurai* and *ponio* to achieve optimal accuracy at minimal cost, targeting performance comparable to or better than existing adaptive-splitting strategies.

# 5 Conclusion

The development of *samurai* and *ponio* stems from a common goal: to provide researchers with reliable, modular tools for designing adaptive numerical solvers. Both libraries emphasize clarity, composability, and reproducibility, making them valuable testbeds for new numerical ideas. Future developments will focus on tighter coupling between spatial and temporal adaptivity, extended parallelization strategies, and the inclusion of additional physical models.
# Acknowledgements

The authors want to thank the NASA hosts for their invitation, and acknowledge the financial support of the CIEDS (project OPEN-NUM-DEF) and the Fondation de l'École polytechnique through the HPC@Maths Initiative.

# References

[1] A. Abdulle. Fourth order Chebyshev methods with recurrence relation. *SIAM Journal on Scientific Computing*, 23(6):2041–2054, 2002.

[2] A. Abdulle and A. A. Medovikov. Second order Chebyshev methods based on orthogonal polynomials. *Numerische Mathematik*, 2001.

[3] A. Abdulle and G. Vilmart. PIROCK: A swiss-knife partitioned implicit–explicit orthogonal Runge–Kutta Chebyshev integrator for stiff diffusion–advection–reaction problems with or without noise. *Journal of Computational Physics*, 242:869–888, 2013.

[4] M. J. Berger and P. Colella. Local adaptive mesh refinement for shock hydrodynamics. *Journal of Computational Physics*, 82(1):64–84, 1989.

[5] P. Bernigaud, A. Clout, D. Biasone, S. Dubois, A. Hoffmann, and M. Massot. Simulation of hybrid dense-rarefied flows relying on relaxation schemes and a hierarchy of moment methods. Technical report, NASA, 2025.

[6] C. Burstedde, L. C. Wilcox, and O. Ghattas. p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees. *SIAM Journal on Scientific Computing*, 33(3):1103–1133, 2011.

[7] A. Cohen, S. Kaber, S. Müller, and M. Postel. Fully adaptive multiresolution finite volume schemes for conservation laws. *Mathematics of Computation*, 72(241):183–225, 2003.

[8] S. Descombes, M. Duarte, and M. Massot. Operator splitting methods with error estimator and adaptive time-stepping. Application to the simulation of combustion phenomena. In R. Glowinski, S. Osher, and W. Yin, editors, *Splitting Methods in Communication, Imaging, Science, and Engineering*, pages 1–13. Springer International Publishing, 2015.

[9] M. P. Duarte. *Adaptive numerical methods in time and space for the simulation of multi-scale reaction fronts*. Theses, Ecole Centrale Paris, Dec. 2011.

[10] L. Gouarin, J. Massot, M. Massot, P. Matalon, L. Series, and C. Tenaud. A dual splitting/ImEx strategy for multicomponent reacting flows with space/time adaptation and error control. Technical report, NASA, 2025.

[11] W. Haegeman, G. Orlando, S. Kokh, M. Massot, L. Gouarin, P. Matalon, and S. Dubois. A two-phase flow model with scale separation for the separated-to-disperse phase transition. Technical report, NASA, 2025.

[12] A. Harten. Multiresolution algorithms for the numerical solution of hyperbolic conservation laws. *Communications on Pure and Applied Mathematics*, 48(12):1305–1342, 1995.

[13] J. D. Lawson. Generalized Runge-Kutta processes for stable systems with large Lipschitz constants. *SIAM Journal on Numerical Analysis*, 4(3):372–380, 1967.

[14] A. Loison, S. Kokh, T. Pichard, and M. Massot. A unified two-scale gas–liquid multi-fluid model with capillarity and interface regularization through a mass transfer between scales. *International Journal of Multiphase Flow*, 177:104857, 2024.

[15] C. D. Meyer, D. S. Balsara, and T. D. Aslam. A stabilized Runge–Kutta–Legendre method for explicit super-time-stepping of parabolic and mixed equations. *Journal of Computational Physics*, 257:594–626, 2014.

[16] H. Sagan. *Space-Filling Curves*. Springer New York, 1994.

[17] A. E. Simon, L. François, A. Asad, E. Quémerais, and M. Massot. Development of a high-order space-time coupling HPC library for multiphysics applications. Technical report, NASA, 2025.

[18] W. Zhang, A. Myers, K. Gott, A. Almgren, and J. Bell. AMReX: Block-structured adaptive mesh refinement for multiphysics applications. *The International Journal of High Performance Computing Applications*, 35(6):508–526, 2021.