samurai

Loïc Gouarin

HPC@Maths team
CMAP / CNRS - Ecole polytechnique

Patched-based representation

Advantages

• compact data structure
• rectangular zones work well for tiling and optimizing caches
• efficient to work with neighbors and the different levels

Disadvantages

• requires more cells than necessary
  • grid hierarchy: overlapping
  • larger refinement zone

Cell-based representation

Advantages

• no grids hierarchy
• keep only needed cells

Disadvantages

• more difficult to find the neighbors
• more difficult to work with the different levels
• tree-like structure
• more difficult to add ghost cells hierarchy

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.

Adaptive methods

AMR

Advantages

• based on heuristic criteria (gradient, second derivative, …)
• easy to implement

Disadvantages

• physical problem understanding is important
• add potentially more cells that needed

MRA

Advantages

• based on wavelet decomposition
• error control
• independent of the physical problem
• only needed cells

Disadvantages

• difficult to implement with the current data structures
• can be costly

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 enable caching and vectorization.
• Manage and facilitate complex operations between meshes (numerical schemes, MRA, …).

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
• x: [0, 4[, [0, 1[, [3, 4[, [0, 1[, [3, 4[, [0, 3[
• y: [0, 4[
• 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[
• y-offset = [0, 1, 2, 4, 5, 6, 7]
• Level 2
• x: [8, 10[, [8, 10[, [14, 16[, [14, 16[
• y: [8, 10[, [14, 16[
• y-offset: [0, 1, 2, 3, 4]

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.

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

\

=

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 projection 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.

Usage example: find jumps

auto jump_set = intersection(
                    translate(mesh[1], {1}).on(0),
                    mesh[0]
                );

Usage example: the projection operator

auto proj_set = intersection(mesh[level], mesh[level + 1])
          .on(level);

proj_set([&](auto i)
{
    u(level, i) = 0.5*(u(level+1, 2*i) + u(level+1, 2*i+1));
});

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.

Parallelism in samurai

CPU

• MPI implementation
• Algebra of set gives the data to exchange
• Various load balancing strategies (WIP)
  • SFC
    • Hilbert
    • Morton
  • diffusion
  • Metis
  • Scotch

GPU

• POC (3 months of work)
  • Two strategies
    • the mesh and the algebra of set keep on the CPU
    • the mesh and an adapted algebra is on the GPU
  • Identify the bottlenecks and the difficulties

Future collaboration with NVIDIA

NumPEx interactions

High-order Space-Time Coupling HPC Library

Project Overview

• Development of a high-order space-time coupling HPC library for multiphysics applications
• Numerical analysis of the multistep coupling scheme in terms of precision and stability
• Implemented in CwIPI² C++ library

Key Achievements

• Dynamic time step adaptation capability
• Demonstrated on conjugate heat transfer test case (solid-fluid)
• Uses dimensionless equations for both systems

Antoine Simon - PhD thesis funded by NumPEx in collaboration with ONERA

A. Asad, L. François, L. Gouarin, M. Massot

High-resolution methods with error control

Context

• Multiscale reaction-diffusion-convection PDEs (combustion, stiff systems)
• Low-Lewis-number mixtures (hydrogen): thermo-diffusive instabilities
• Need an a priori control of errors both in time and space

Objective

• New generation of high-resolution methods with error control
• New generation of open-source simulation codes
• Applications: hydrogen risk (CEA), plasma & solar physics, nonlinear chemistry

Approach

• ImEx (Implicit/Explicit) additive Runge-Kutta schemes
• Controlled operator coupling instead of operator splitting
• Couple ImEx integrators with adaptive multiresolution (samurai, ponio)
• Discontinuous Galerkin + multiwavelets for local \(h\)-\(p\) control

Framework

• M2 internship \(\rightarrow\) PhD thesis: Fadel Kicha
• Funded by NumPEx – EXA-MA project
• CMAP / HPC@Maths & EM2C
• L. Gouarin, J. Massot, M. Massot, L. Series, C. Tenaud

This is an M2 internship subject within the EXA-MA project, aiming at high-resolution methods with a priori error control for stiff multiscale PDEs. The idea is to move from operator splitting to fully coupled ImEx Runge-Kutta integrators, and to combine them with adaptive multiresolution (samurai, ponio) and Discontinuous Galerkin schemes for an h-p control in space.

Exa-DI interaction - GT Block-structured AMR

Project overview

• Identify open source software to compare
• Define physical problem suitable for AMR benchmarking
• Numerical scheme selection
• Establish performance metrics

Key achievements

• First tests with Euler equations
• 4 selected software
  • Dyablo
  • samurai
  • AthenaK
  • Astro

Future Work

• Define reference solutions (double mach reflection, Riemann problems,…)
• Compare accuracy and performance of the different AMR methods and codes

Exa-DOST interaction - AMR I/O

Objective

• Identify I/O requirements for AMR libraries
• Define a common I/O strategy for both AMR libraries
• Allow data compression without or with loss

Methodology

• Bibliography on wavelet compression
• Analysis of various formats for AMR data

Sylvain Joube (Post-doc - 24 months - started Nov. 2025)

J. Bigot, M. Delorme, A. Durocher, L. Gouarin, M. Stauffert

Exa-SofT interaction - Task-based parallelism

Objective

• Express samurai’s adaptive computations as a StarPU task graph
• Target asynchronous and heterogeneous (CPU/GPU) execution
• Improve load balancing across refinement levels

Methodology

• Identify the right task granularity from the interval-based data structure (per level / per sub-domain)
• Delegate scheduling and dependencies to the StarPU runtime

Involved Person

• Théo Grandsart (M2 internship - PhD funding secured)
• Nathalie Furmento
• Loic Gouarin
• Marc Massot
• Raymond Namyst
• Pierre-André Wacrenier

People involved

• 4 core developers
• users and developers communities
• Various collaborations and research programs

Real-World Applications

• Interfacial flow simulation
• Simulation analysis on the Hydrogen risk
• Plasma discharges simulations
• Lithium battery simulation
• Combustion DNS
• Code coupling simulations

Samurai

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

Hands-on samurai

https://hpc-maths.github.io/2025-hands-on-samurai