samurai - the present and beyond

Loïc Gouarin
Pierre Matalon

New containers

  • Add flexibility in container library dependencies
  • Transparent for field usage
  • Two available implementations : xtensor and Eigen
  • In a near future, Kokkos implementation will be provided

Algebra of sets

The algebra is not complete in samurai.

For example, it is not possible to write

auto set = translation(intersection(mesh1, mesh2), {1, 0});

A comprehensive rewrite of the system is underway to address these limitations.

With this new implementation, we can write:

auto set = translate(intersection(translate(mesh1, {-2}).on(4),
                                  translate(mesh2,  {2}).on(3)),
                     {5}).on(10);

Load balancing

The goal is to evenly divide the workload amongst processors.

Mesh parts (and associated data) must be exchanged between neighbours.

Load balancing using Morton curve

Load balancing using Hilbert curve

Loïc Strafella’s work

Load balancing


What was achieved?

  • Load balancing using different algorithms
    • Morton curve
    • Hilbert curve
    • Diffusion algorithm

What remains?

  • All the techniques previously used are based on a exchanging cells
  • Understand the differences between this view and an interval view
  • Try various strategies and find a generic criteria to achieve good load balancing

Flux computation at the finest level

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

Flux computation at the finest level

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

Flux computation at the finest level

What was achieved?

  • Implementation of a general algorithm which can reconstruct the solution at upper/lower levels
    • used for accurate flux computation
    • used for transferring a field from one mesh to another
    • used to reconstruct the solution at the finest level

What remains?

  • More complicated in higher dimensions…
  • Option to de/activate the flux computation at the finest level
  • Fix the graduation according to the prediction order

Performance optimization

2 new recruits: Sébastien Dubois and Alexandre Hoffmann

  • xtensor
  • changes in internal data structures and algorithms
  • vectorization / batch processing

Packaging

samurai is available on conda-forge and conan.

In the NumPEx context, samurai must also be installed with Spack and Guix. In the next few months, we will provide a way to install it through these package managers.

Restart

It’s coming!

… but some issues remain:

  • How to select a list of fields in a restart file?
  • How to give them back to the user?
  • Can we have a generic output for restart and plot?
  • Is it convenient to have post processing scripts?

Samurai/ponio

An 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

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

NASA Roadmap covering the five projects

Requested

  • Flux computation at the finest level
  • Restart
  • Sequential optimizations
  • Load balancing
  • Various finite volume schemes available out of the box (MUSCL, Rusanov, HLL, …)
  • Monotony-preserving prediction (MUSCL?)
  • Specific scaling according to the spatial direction
  • Distinguish between a scalar field and a vector field

Nice-to-have

  • Post-processing scripts
  • Problem greater than 3 dimensions (e.g. plasma problems)

ponio

What you can use now?

  • A large collection of explicit and diagonally implicit Runge-Kutta methods with adaptive time step
  • Splitting methods (Lie/Strang)
  • Some extended stability methods (RKC, ROCK2, ROCK4 and RKL), including adaptive time step for ROCK 2 and 4
  • ImEx strategy with PIROCK (see BZ example)

What is planned?

  • Adding adaptive time steps to PIROCK (already available in ROCK2 and ROCK4)
  • Adaptive Strang method
  • Adding projection operator between stages of RK method (for two-phase flow simulations)
  • Exploring time-space strategies for asymptotic preserving schemes, or HLL scheme (for plasma project)