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

HAL
open science

# A Parareal in time numerical method for the collisional Vlasov equation in the hyperbolic scaling

Tino Laidin, Thomas Rey

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

## To cite this version:

Tino Laidin, Thomas Rey. A Parareal in time numerical method for the collisional Vlasov equation in the hyperbolic scaling. 2025. hal-04929730

HAL Id: hal-04929730

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

Preprint submitted on 4 Feb 2025

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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A PARAREAL IN TIME NUMERICAL METHOD FOR THE
COLLISIONAL VLASOV EQUATION IN THE HYPERBOLIC SCALING

TINO LAIDIN AND THOMAS REY

ABSTRACT. We present the design of a multiscale parareal method for kinetic equations in the fluid dynamic regime. The goal is to reduce the cost of a fully kinetic simulation using a parallel in time procedure. Using the multiscale property of kinetic models, the cheap, coarse propagator consists in a fluid solver and the fine (expensive) propagation is achieved through a kinetic solver for a collisional Vlasov equation. To validate our approach, we present simulations in the 1D in space, 3D in velocity settings over a wide range of initial data and kinetic regimes, showcasing the accuracy, efficiency, and the speedup capabilities of our method.

KEYWORDS: Parareal algorithm, Vlasov-BGK equation, collisional kinetic equation, multiscale model, rarefied gas dynamics, OpenMP, MPI.

2020 MATHEMATICS SUBJECT CLASSIFICATION: 82B40, 82C70 76P05, 65M08.

# 1. INTRODUCTION

The simulation of many complex phenomena involving interacting particles can be modelled through either fluid or kinetic descriptions. However, the validity of the models varies highly depending on the application. In particular, to accurately describe rarefied gases or charged particles in a device, a fluid description given by the Euler, Navier-Stokes or drift-diffusion equations may break down. Typically, it occurs around shocks or because of boundary layers and a kinetic description is therefore necessary. Nevertheless, the main drawback of the latter approach lies in its numerical cost. Indeed, deterministic kinetic simulations suffer even more strongly than usual solvers from the curse of dimensionality because of the large size of the phase space. As a consequence, one wants to resort to kinetic simulation sparingly.

In this work, we are interested in the approximation of solutions to the following scaled collisional Vlasov equation:

$$ egin{cases} \partial_t f^\epsilon(t, x, v) + \frac{1}{\epsilon^\alpha} (v \cdot \nabla_x f^\epsilon(t, x, v) - E \cdot \nabla_v f^\epsilon(t, x, v)) = \frac{1}{\epsilon^{\alpha+1}} Q(f^\epsilon), \\ f(0, x, v) = f_0(x, v), \end{cases} \quad (1) $$

where $f(t, x, v)$ is the distribution function of particles, $t \ge 0$, $x \in \Omega_x \subset \mathbb{R}^{d_x}$, $v \in \mathbb{R}^{d_v}$. In position, the domain $\Omega_x$ is bounded, and boundary conditions will be presented later. This equation is used to model several phenomena including rarefied gases ($E = 0$) or charged particles in a tokamak, or electronic device where $E$ typically is gradient of a potential that solves a Poisson equation. The parameter $\epsilon$ is the Knudsen number, which is the ratio between the mean free path of particles and the typical length scale of observation. Depending on the value of $\alpha$, the asymptotic model as $\epsilon$ tends to 0 will either be given by a hydrodynamic system ($\alpha = 0$) [3, 25], such as Euler's equations, or a diffusive one ($\alpha = 1$) [45, 46], typically a drift-diffusion equation.

A main component of (1) is the Boltzmann-like collision operator $Q$. In the following, we will assume that it satisfies the following classical properties:

<footnote>(Tino Laidin) UNIV BREST, CNRS UMR 6205, LABORATOIRE DE MATHÉMATIQUES DE BRETAGNE AT-LANTIQUE, F-29200 BREST, FRANCE</footnote>

<footnote>(Thomas Rey) UNIVERSITÉ CÔTE D'AZUR, CNRS, LJAD, PARC VALROSE, F-06108 NICE, FRANCE</footnote>

<footnote>E-mail addresses: tino.laidin@univ-brest.fr, thomas.rey@univ-cotedazur.fr.</footnote>
* It preserves mass, momentum and energy:

$$ \int_{\mathbb{R}^3} Q(f) \begin{pmatrix} 1 \\ v \\ \frac{|v|^2}{2} \end{pmatrix} \mathrm{d}v = 0_{\mathbb{R}^5}. \qquad (2) $$

* It dissipates the Boltzmann entropy:

$$ \int_{\mathbb{R}^3} Q(f) \log(f) \mathrm{d}v \le 0. \qquad (3) $$

* Its equilibria are given by Maxwellian distributions:

$$ Q(f) = 0 \iff f = M_{\rho, u, \theta}(v) := \frac{\rho}{(2\pi\theta)^{d/2}} \exp\left(-\frac{|v-u|^2}{2\theta}\right), \quad \forall v \in \mathbb{R}^3, \qquad (4) $$

where $\rho$, $u$ and $\theta$ denote the density, mean velocity and temperature respectively. In regard to property (4) one obtains, at least formally [29], that in the limit $\varepsilon \to 0$, the distribution $f^\varepsilon$ tends towards a Maxwellian distribution whose moments are solution to a fluid model. In the hydrodynamic scaling ($\alpha = 0$), the moments solve Euler's equations for a monoatomic gas:

$$ \left\{ \begin{array}{l} \partial_t \rho + \nabla_x \cdot (\rho u) = 0, \\ \partial_t (\rho u) + \nabla_x \cdot (\rho u \otimes u) + \nabla_x (\mathbb{P}) - \rho E = 0, \\ \partial_t \mathcal{E} + \nabla_x \cdot ((\mathcal{E} + \mathbb{P})u) - \rho u \cdot E = 0, \end{array} \right. \qquad (5) $$

where $\mathbb{P} = \rho\theta I$ is the pressure and the energy $\mathcal{E}$ is related to the temperature through the relation,

$$ \mathcal{E} = \frac{1}{2} (\rho|u|^2 + d_v \rho\theta). $$

In the diffusive scaling ($\alpha = 1$), they solve a drift-diffusion equation ($u = 0, \theta = 1$):

$$ \partial_t \rho - m_2 \mathrm{div}_x J = 0, \quad \text{where } J = \nabla_x \rho - E\rho, \qquad (6) $$

and $m_2$ is the space averaged second moment. A widely used collision operator that we consider in this work, and which is simpler than the quadratic Boltzmann operator, is the so-called BGK operator [5]. It satisfies the above properties and is a relaxation operator with Maxwellian equilibrium:

$$ Q(f) = \frac{\tau}{\varepsilon}(M_{\rho,u,\theta} - f), \qquad (7) $$

where $\tau$ depends on the moments $\rho$ and $\theta$. This simple model gives the correct asymptotic models as $\varepsilon \to 0$, and an extended study using a Chapman-Enskog expansion [50, 29] shows that higher order systems can be obtained. For example, in the hydrodynamic scaling, one can obtain the Navier-Stokes equations:

$$ \left\{ \begin{array}{l} \partial_t \rho + \nabla_x \cdot (\rho u) = 0, \\ \partial_t (\rho u) + \nabla_x \cdot (\rho u \otimes u) + \nabla_x (\mathbb{P}) - \rho E = \varepsilon \nabla_x \cdot (\mu\sigma(u)), \\ \partial_t \mathcal{E} + \nabla_x \cdot ((\mathcal{E} + \mathbb{P})u) - \rho u \cdot E = \varepsilon \nabla_x \cdot (\mu\sigma(u)u + \kappa \nabla_x \theta), \end{array} \right. \qquad (8) $$

with the strain rate tensor $\sigma(u)$ given by

$$ \sigma(u) = (\nabla_x u + (\nabla_x u)^T) - \frac{2}{3}(\nabla_x u)I. $$

Several approaches have been designed to significantly reduce the computational cost of kinetic solvers. An inspiration for this work are the so-called *hybrid methods*. The idea behind these approaches is to achieve a coupling between a cheap, low dimensional, fluid model and the expensive, high dimensional, kinetic one. This type of techniques [53, 8, 18] mainly relies on some kind of domain decomposition in position or velocity. To accurately describe the solution, domain indicators that assess the validity of the cheap fluid description over the expensive kinetic one are often needed [38, 52, 10, 32]. In addition to the computation of subdomains (that can be costly depending on the physical regime considered), one must also deal with boundary conditions between different kinds of models.

The last type of methods aiming at reducing the cost of kinetic simulations that we would like to mention are the so-called Dynamical Low Rank methods, that directly deal with the matrix of the discretized system [31, 15] through some type of SVD-like decomposition.
**The parareal method.** Our goal in this work is to take a new strategy, which focuses on the time integration of the multiscale kinetic equation considered. More precisely, we design a *Parallel in Time* (PinT) method that takes advantage of parallel computing architectures. Such approximation paradigm gained a lot of popularity in the past years with the introduction of the so-called parareal method in the seminal series of work [39, 41]. It consists in parallelizing the time integration of a dynamical system using an iterative process based on a predictor-corrector framework, but has never been used in the framework of collisional kinetic equation. Let us mention for the interested reader the very recent book [23] that provides a comprehensive introduction to this topic.

Let us now present the parareal approach for the simple ODE

$$ \begin{cases} \frac{d}{dt}u(t) = f(u), & t \in [0, T], \\ u(0) = u^0. \end{cases} $$

The goal of the parareal algorithm is to approximate the solution at some fixed discrete times $T^n$,

$$ U^n \approx u(T^n), \quad n \in \{0, \dots, N_g\}, \quad N_g \in \mathbb{N}. $$

To do so, instead of using a single time integrator (or propagator), two solvers are coupled through an iterative process: the *parareal iterations*. The first solver will be denoted by $\mathcal{G}(T^n, T^{n+1}, U^n)$. It must be cheap to compute numerically, but yet convergent towards the solution to the ODE. It will provide some numerical *guesses* of the behavior of the solution on the coarse grid. The second solver, that we will denote by $\mathcal{F}(T^n, T^{n+1}, U^n)$, must be very accurate, and will be used to improve the aforementioned guesses.

While the propagator $\mathcal{G}$ is less costly, its drawback is naturally its precision. However, it will be corrected by the accurate, computationally expensive, propagator $\mathcal{F}$. The algorithm then unfolds as follows:

(1) Divide the time domain in $N_g \in \mathbb{N}$ subintervals $[T^n, T^{n+1}]$, $n \in \{0, \dots, N_g - 1\}$;

(2) Perform a first coarse guess:

$$ U^{n+1,0} = \mathcal{G}(T^n, T^{n+1}, U^{n,0}) = u^0; $$

(3) Refine the guess through the parareal iterations:

$$ U^{n+1,k+1} = \mathcal{G}(T^n, T^{n+1}, U^{n,k+1}) + \mathcal{F}(T^n, T^{n+1}, U^{n,k}) - \mathcal{G}(T^n, T^{n+1}, U^{n,k}), $$

where $k = 1, 2, \dots$ denotes the $k^{th}$ parareal iteration.

The advantage of this iterative algorithm is its ability to compute the expensive fine propagations in parallel. Therefore, as long as the number of parareal iterations remains small enough to obtain a given accuracy, one can expect a reduction of the computational cost of the time integration.

Parallel in time methods are now widely utilized and studied in various contexts. The convergence of the algorithm was investigated in [24, 21, 11], including variations of the original method. The uses are plenty, and we refer to [4] for neutron transport, [19, 49] for the Navier-Stokes equations, [22] for Hamiltonian systems and [47] for turbulence in plasmas. We also refer to a recent work on the acceleration of the method [40, 42]. It was furthermore observed that the parareal method suffers from poor convergence for hydrodynamic systems where convection dominates diffusion [20, 14, 43, 2, 16]. In particular, in the case of Burgers' equation it was noticed that for a fixed discretization, the longer the integration domain, the slower was the convergence [21].

In the context of multiscale systems, the parareal algorithm has been used in [11, 34, 26, 27, 48, 6] by leveraging reduced order models as coarse solvers. The idea is to use a simpler model, instead of a simpler time integrator, as a coarse propagator and a richer model for the fine propagations. These ideas are reminiscent to the ones introduced in the co-called Heterogeneous Multiscale Method framework from the series of works [13, 12, 1].

Expanding on these advancements, our main goal in this work is to build a multiscale parareal method for kinetic equations to perform accurate simulations of the observable macroscopic moments. A fluid model, either given by a hydrodynamic or diffusive limit of the underlying kinetic equation, is used as a coarse propagator. The fine propagation is in turn achieved by the resolution of the kinetic model. To this aim, we resort to *Asymptotic*
Preserving (AP) schemes [30, 28] that remain stable for any value of the Knudsen number $\epsilon$ in (1). While they provide an accurate approximation in the limit $\epsilon \to 0$, their main drawback is that their cost remains the one of a kinetic scheme even in a fluid regime. We will demonstrate in our numerical experiments that our new method can solve this issue.

**Plan of the paper.** This work is organized as follows. In Section 2 we present the multiscale parareal algorithm for kinetic equations and discuss some of its properties. Sections 3 is dedicated to the implementation of the method and in particular to the presentation of the solvers used. We then discuss the parallelization of the method in Section 4. Finally, we thoroughly assess the accuracy and speedup properties of the method through numerical experiments in Section 5.

## 2. MULTISCALE PARAREAL ALGORITHM

We introduce in this section the multiscale parareal algorithm. The formalism we shall employ follows the one introduced in [34] and the references therein.

### 2.1. Link between the scales.
The crucial point of a multiscale method is the link between the different scales of the model. In the context of kinetic equations a natural way to go from kinetic to fluid is to consider the projection of a distribution $f$ towards its moments $U = (\rho, u, \theta)^T$ defined in Definition 1.

**Definition 1.** For a distribution function $f = f(t, x, v) \in L^1((1 + v)^2 \, d\nu)$, we define the projection $\mathcal{P}f(t, x) = U(t, x)$ as

$$ \mathcal{P}f = U = \begin{pmatrix} \int_{\mathbb{R}^3} f \, d\nu \\ \frac{1}{\rho} \int_{\mathbb{R}^3} v f \, d\nu \\ \frac{1}{3\rho} \int_{\mathbb{R}^3} |v - u|^2 f \, d\nu \end{pmatrix} = \begin{pmatrix} \rho \\ u \\ \theta \end{pmatrix}. $$

Note that in Definition 1, we start from a distribution that contains all the information and project it onto a manifold where the moments solve a system valid only for Maxwellians. On the opposite, to lift a macroscopic data $U(t, x) \in \mathbb{R}^5$ to a distribution $f(t, x, v)$ a first idea is to consider the validity of the fluid model used. In the case of Euler's equations, we therefore chose to reconstruct a Maxwellian whose moments are given by $U$.

**Definition 2.** For a macroscopic data $U(t, x) \in \mathbb{R}^5$, we define the lifting $\mathcal{L}U(t, x, v) = f(t, x, v)$ as

$$ \mathcal{L}U = f = \frac{\rho}{(2\pi\theta)^{3/2}} \exp\left(-\frac{|v-u|^2}{2\theta}\right). $$

**Remark 1.** Note that a finer way to handle the lifting and fluid propagations is to consider the Chapman-Enskog expansion of the distribution:

$$ f^{\epsilon} = \mathcal{M}_{\rho^{\epsilon}, u^{\epsilon}, \theta^{\epsilon}} + \sum_{l=1}^{\infty} \epsilon^{l} g^{(l)}, \qquad (9) $$

where the perturbations $g^{(l)}$ can be explicitly computed and depend on $f$ only through its moments (some explicit computations for the BGK model we shall consider in this work can be found in the book [50]).

We have already mentioned that in the hydrodynamic scaling, the moments of a Maxwellian distribution are solution to Euler's equations. By considering the first order perturbation in (9), these moments can be shown to solve the Compressible Navier-Stokes equations (8) (CNS) that can be understood as a first order correction of Euler's equations. Consequently, a finer method would be to use a CNS solver as a coarse propagator. Note that even with this choice, the cost of the coarse integration remains much lower than a kinetic propagation. In addition, instead of only lifting the moments to a Maxwellian, one can consider adding the perturbations $g^{(l)}$. It allows for a finer reconstruction of the distribution and the consideration of far from equilibrium phenomena.
**Definition 3.** With the same hypothesis of Definition 2, we define the lifting

$$ \mathcal{L}^{(L)} U(t, x, v) = f(t, x, v) $$

of order L by setting

$$ f = \mathcal{L}^{(L)} U = \frac{\rho}{(2\pi\theta)^{3/2}} \exp\left(-\frac{|v-u|^2}{2\theta}\right) + \sum_{l=1}^{L} \varepsilon^l g^{(l)}. \quad (10) $$

2.2. **Semi-discretization in time.** Let us now consider the time interval $[0, T]$ and divide it in $N_g \in \mathbb{N}$ uniform subintervals $[T^n, T^{n+1}]$. We are interested in the approximation of the moment vector $U$ at the fixed discrete times $T^n$. We denote by $U^{n,k}$ this approximation at a parareal iteration $k$. Furthermore, we also need to introduce a fine discretization that will be used by the fine propagator. Let $N_f \in \mathbb{N}$ denote the number of uniform fine subintervals $[t^n, t^{n+1}]$. Typically, $N_f$ shall be greater than $N_g$. We can then define the coarse (resp. fine) time steps and times:

$$ T^n = n\Delta t_g, \quad \Delta t_g = \frac{T}{N_g}, \quad t^n = n\Delta t_f, \quad \Delta t_f = \frac{T}{N_f}. $$

Note that the output of the algorithm will therefore contain $N_g$ fixed snapshots at times $T^n$ of the moments. We briefly omit the definitions of the kinetic and fluid propagators $\mathcal{F}$ and $\mathcal{G}$ as they will be discussed below. However, it is important to mention that both the coarse and fine time steps are in practice the maximum time steps allowed within a propagator. The local time steps inside the solvers are naturally subject to a stability condition, and we define

$$ \Delta t_{\text{loc}}^{n, \mathcal{G}} = \min \{ \Delta t_{\text{stab}}^{n, \mathcal{G}}, \Delta t_g \}, \quad \Delta t_{\text{loc}}^{n, \mathcal{F}} = \min \{ \Delta t_{\text{stab}}^{n, \mathcal{F}}, \Delta t_f \}, $$

where $\Delta t_{\text{stab}}^{n, \mathcal{G}}$ and $\Delta t_{\text{stab}}^{n, \mathcal{F}}$ denote the time steps prescribed by the stability conditions of the numerical schemes used. The method can then be summarized in Algorithm 1.

<caption>**Algorithm 1** Multiscale kinetic parareal Algorithm</caption>

**Require:** $U^{0,0}$
1: **for** $n = 1, \dots, N_g$ **do**
2: &nbsp;&nbsp;&nbsp;&nbsp; $U^{n,0} \leftarrow \mathcal{G}(U^{n-1,0})$
3: **end for**
4: **while** $k \le K$ **or** error $\ge \text{tol}$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp; **for** $n = 1, \dots, N_g$ **do**
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\Delta^n = \mathcal{P}\mathcal{F}\mathcal{L}(U^{n-1,k-1}) - \mathcal{G}(U^{n-1,k-1})$
7: &nbsp;&nbsp;&nbsp;&nbsp; **end for**
8: **for** $n = 1, \dots, N_g$ **do**
9: &nbsp;&nbsp;&nbsp;&nbsp; $U^{n,k+1} = \mathcal{G}(U^{n-1,k}) + \Delta^n$
10: **end for**
11: Compute successive errors on the moments and $k \leftarrow k + 1$
12: **end while**

As a stopping criterion of the algorithm, we simply consider the error between two successive parareal iterations:

$$ \text{error} = \max_n |U^{n,k+1} - U^{n,k}|, $$

or a fixed number of iterations $K \in \mathbb{N}$.

An interesting property of Algorithm 1 is that, at parareal iteration $k$, one can in fact start the loops (lines 5 and 8) at $k$ instead of 1. Indeed, since one propagates a data from time 0 that is fixed, it is unnecessary to recompute the first 1 to $k$ iterations as they would not be modified. This can be seen through a simple ODE example by focusing on the first time intervals and parareal iterations: for $n = 1$, to update $U^{1,k+1}$, one needs to compute

$$ U^{1,k} = \mathcal{G}(T^0, T^1, U^{0,1}) + \mathcal{F}(T^0, T^1, U^{0,0}) - \mathcal{G}(T^0, T^1, U^{0,0}). $$
Since $U^{0,k} = U^0$ for all $k$, the two coarsely propagated data cancel each other, and we obtain that the first parareal iteration automatically corrects the initial guess towards the fine propagation:

$$U^{1,k} = \mathcal{F}(T^0, T^1, U^{0,0}), \quad \forall k \ge 1.$$

Consequently, when computing the next parareal iteration $k = 2$, $U^{1,k}$ is now fixed and one can make the same observation as before but now for $n = 2$:

$$U^{2,k} = \mathcal{F}(T^1, T^2, U^{1,k}) = \mathcal{F}(T^1, T^2, \mathcal{F}(T^0, T^1, U^{0,0})), \quad \forall k \ge 2.$$

By induction, at parareal iteration $k$, the $1$ to $k$ first times are therefore fixed, and don't need to be recomputed. A consequence is that after $N_g$ parareal iterations, the outcome of Algorithm 1 is the same as a fully fine simulation.

Building upon these observations, an optimized version of this algorithm is given by Algorithm 2.

<caption>**Algorithm 2 Optimized multiscale kinetic parareal Algorithm**</caption>

**Require:** $U^{0,0}$
1: **for** $n = 1, \dots, N_g$ **do**
2: &nbsp;&nbsp;&nbsp;&nbsp; $U^{n,0} \leftarrow \mathcal{G}(U^{n-1,0})$
3: **end for**
4: **while** $k \le K$ **or** error $\ge \text{tol}$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp; **for** $n = k, \dots, N_g$ **do**
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\Delta^n = \mathcal{P}\mathcal{F}\mathcal{L}(U^{n-1,k-1}) - \mathcal{G}(U^{n-1,k-1})$
7: &nbsp;&nbsp;&nbsp;&nbsp; **end for**
8: **for** $n = k, \dots, N_g$ **do**
9: &nbsp;&nbsp;&nbsp;&nbsp; $U^{n,k+1} = \mathcal{G}(U^{n-1,k}) + \Delta^n$
10: **end for**
11: Compute successive error on the moments and $k \leftarrow k + 1$
12: **end while**

$\triangleright$ First coarse guess

$\triangleright$ Parareal iterations

$\triangleright$ Compute the jumps in parallel

$\triangleright$ Sequential correction

From an optimization point of view, note also that steps 5-7 of both Algorithms 1 and 2 contain all the expensive computations. More particularly we chose to do the lifting and the projection alongside the kinetic propagation to better distribute the workload. Indeed, since the iterations of this loop are independant from one another, which is the core of the parareal method, it can be parallelized to reduce the computation time.

### 3. NUMERICAL SCHEMES

In this section, we shall detail our choice of fine and coarse propagators as well as the discrete versions of the lifting and projection operators.

From now on, let us place ourselves in the $d_x = 1$ and $d_v = 3$ setting that we shall use for our numerical experiments. The $(x, v)$-phase space is discretized in a finite volume fashion using control volumes $K_{ij} = \mathcal{X}_i \otimes \mathcal{V}_j$, for $i \in \{1, \dots, N_x\}$ and $j = (j_x, j_y, j_z) \in \{1, \dots, N_v\}^3$. The discrete unknowns are defined as

$$f_{ij}^n \approx \frac{1}{\Delta x_i \Delta v_j} \int_{K_{ij}} f(T^n, x, v) \, dx \, dv, \quad U_i^n \approx \frac{1}{\Delta x_i} \int_{\mathcal{X}_i} U(T^n, x) \, dx,$$

with $\Delta x_i$ the volume of $\mathcal{X}_i$ and $\Delta v_j$ the volume of the cube $\mathcal{V}_j$.

#### 3.1. Kinetic schemes.
The main difficulty to approximate scaled kinetic equations is the stability of the time integration with respect to the small parameters. In both hydrodynamical and diffusive scalings, AP schemes must therefore be considered to ensure tractable simulations for any value of the Knudsen number.

**Hydrodynamic scaling.** In the hydrodynamic scaling, several AP schemes have been developed, and most of them rely on an implicit-explicit (IMEX) type time discretization
[44]. In this work, we consider a very simple IMEX discretization where the stiff relaxation term is implicit. The semi-discretized scheme can be written as:

$$ \frac{f^{n+1} - f^n}{\Delta t} + [v \cdot \nabla_x f^n + E \cdot \nabla_v f^n] = \frac{\tau^{n+1}}{\varepsilon} (M_{\rho^{n+1}, u^{n+1}, \theta^{n+1}} - f^{n+1}). \quad (11) $$

Thanks to the conservation properties of the equation, the moments $\rho$, $u$ and $\theta$ at time $n+1$ used to build the Maxwellian $M^{n+1}$ can be obtained explicitly, yielding an efficient, fully explicit in time procedure [7, 17]. More advanced methods such as a micro-macro approach could also be considered here [36]. For the phase space discretization, the numerical fluxes are of upwind type for the transport in position and of Rusanov type for the transport in velocity [37, 54]. The fully discretized kinetic system now reads

$$ \frac{f_{ij}^{n+1} - f_{ij}^{n}}{\Delta t} + \left[ \mathcal{F}_{i+\frac{1}{2},j}^{n} - \mathcal{F}_{i-\frac{1}{2},j}^{n} + \mathcal{G}_{j+\frac{1}{2}}^{n} - \mathcal{G}_{j-\frac{1}{2}}^{n} \right] = \frac{\tau_{i}^{n+1}}{\varepsilon} \left( \mathcal{M}_{\rho_{i}^{n+1},u_{i}^{n+1},\theta_{i}^{n+1}} - f_{ij}^{n+1} \right), \quad (12) $$

with the numerical fluxes

$$ \mathcal{F}_{i+\frac{1}{2},j}^{n} = v^{+}f_{ij}^{n} + v^{-}f_{i+1,j}^{n}, \quad \mathcal{G}_{i,j+\frac{1}{2}}^{n} = \frac{E_i}{2}(f_{ij}^{n} + f_{i,j+1}^{n}) + \frac{E_{max}}{2}(f_{i,j+1}^{n} - f_{i,j+1}^{n}), \quad (13) $$

where $a^+ = \max\{a, 0\}$, $a^- = -\min\{a, 0\}$ and $E_{max} = \max E_i$.

It is widely known that such an approach yields a stable, yet first order discretization of the collisional kinetic equation (1). A natural improvement to this fine propagator would be to use fluxes limiters [55].

**Diffusive scaling.** In the diffusive scaling, one could adopt the micro-macro approach introduced in [36, 9] together with an exponential time integrator for the microscopic equation [35, 33]. This approach will be explored in a future work.

**3.2. Fluid schemes.** Coarse propagators aim at giving a decent guess of the behavior of the moments at a low numerical cost. We distinguish two types of schemes, depending on the scaling considered at the kinetic level.

**Systems of conservation laws.** In the framework of multiscale kinetic equations, systems of conservation laws typically arise as hydrodynamic limits of kinetic equations. They can generally be written as:

$$ \partial_t U + \nabla_x \cdot F(U) = S(U), \quad (14) $$

where $F$ is the flux function and $S$ a potential source term [37, 54]. As a coarse propagator, we rely on a simple finite volume discretization with Rusanov-type fluxes.

The Euler equations in this setting can be expressed as

$$ \frac{\partial V}{\partial t} + \frac{\partial h(V)}{\partial x} = 0, $$

where

$$ V = \begin{bmatrix} \rho \\ \rho u_x \\ \rho u_y \\ \rho u_z \\ \mathcal{E} \end{bmatrix}, \quad h(V) = \begin{bmatrix} \rho u_x \\ \rho u_x^2 + \rho \theta \\ \rho u_x u_y \\ \rho u_x u_z \\ u_x(\mathcal{E} + \rho \theta) \end{bmatrix} $$

One can use the following explicit formula to compute the coarse integrator

$$ V_{i,j}^{n+1} = V_{i,j}^{n} - \left[ H_{i+\frac{1}{2}}^{n} - H_{i-\frac{1}{2}}^{n} \right], \quad (15) $$

where the Rusanov flux $H$ is given by

$$ H_{i+\frac{1}{2}}^{n} = \frac{1}{2}(h(V_{i+1}^{n}) - h(V_{i}^{n})) - \frac{1}{2}S_{i+\frac{1}{2}}^{n}(U_{i+1}^{n} - U_{i}^{n}), $$

with

$$ S_{i+\frac{1}{2}}^{n} = \max(|u_{x,i+1}^{n}| + c_{i+1}, |u_{x,i}^{n}| + c_{i,j}), $$

and $c_i^n = \sqrt{\theta_i^n}$ is the sound speed.

**Drift-diffusion equation.** Explicit discretizations of diffusive systems usually suffer from a parabolic stability condition. However, since the cost of such a scheme is essentially
negligible compared to a full kinetic simulation, we do not consider an implicit-in-time procedure. The fluid system is therefore approximated using a central finite difference scheme that is nothing but the asymptotic limit of the micro-macro scheme mentioned above [35, 33].

**Remark 2.** As is the case for the fine propagator, it is common to consider the use of slope of flux limiters when solving systems of conservation laws. Nevertheless, we noticed in practice that additional spatial accuracy on the predictor does not accelerate the convergence of the parareal algorithm. On the contrary, a slower convergence was observed. Our interpretation is that the numerical viscosity is beneficial. This corroborates with the fact the parareal algorithm is known to work best when the solution is regular, as was already mentioned in the series of works [20, 14, 43, 2, 16].

### 3.3. Discrete lifting and projection operators.
Let us finish the complete description of our numerical propagators by describing the discrete lifting and projection operators.

**Discrete lifting operator.** A natural way to define the discrete lifting operator is to compute the pointwise discrete Maxwellian associated to the moments $U_i = (\rho_i, u_i, \theta_i)^T$, $i \in \{1, \dots, N_x\}$, namely

$$ (\mathcal{L}U)_{ij} = \frac{\rho_i}{(2\pi\theta_i)^{3/2}} \exp\left(-\frac{|v_j - u_i|^2}{2\theta_i}\right), \quad j = (j_x, j_y, j_z) \in \{1, \dots, N_v\}^3. \quad (16) $$

We could also define higher order lifting as defined in the continuous case (3) and proceed in the same way using pointwise computations and finite difference approximations of the derivatives.

**Discrete projection operator.** To project a given discrete distribution $f_{ij}$ defined on a discrete phase space towards its discrete moments $(\rho_i, u_i, \theta_i)^T$, we consider a simple first order quadrature:

$$ \rho_i = \sum_j f_{ij} \Delta v_j, \quad (\rho u)_i = \sum_j v_j f_{ij} \Delta v_j, \quad 3(\rho \theta)_i = \sum_j |v_j - u_i|^2 f_{ij} \Delta v_j. $$

Note that more accurate quadratures could be considered depending on the velocity discretization and a desired accuracy.

## 4. PARALLELIZATION OF THE METHOD

In order to achieve a significant gain in computational time, one needs to rely on an efficient parallelization. The method is currently implemented using an OpenMP paradigm to validate its concept. Its real performance should stand out when properly deployed on an heterogeneous supercomputer along with an MPI parallelization over several threads. The parallelization of this type of method on distributed memory architecture is known to be challenging and such an implementation is briefly discussed in Appendix A.

**Theoretical speed up.** In an ideal setting, i.e. by omitting the cost of communications and synchronizations between cores/threads, one can derive an estimate on the number of parareal iterations to obtain a speedup with PinT-types methods.

Let us denote by $T_{\text{Kin}}$ (resp. $T_{\text{Fluid}}$) the maximum time for the fine, kinetic (resp. coarse, fluid) solver to evolve an initial data from time $T^n$ to time $T^{n+1}$. Note that we take the maximum over all propagations because, on each sub time-interval, the stability condition may vary, impacting the optimal local time step. In addition, one needs to take into account the cost of the projection and lifting operators since they are fully non-local and can therefore be significantly expensive. We shall denote by $T_{\text{Lift}}$ and $T_{\text{Proj}}$ these costs. The ideal numerical cost of Algorithm 1 on $N_p$ threads with $N_g$ sub-time intervals is then given by the formula:

$$ T_{\text{Parareal}} = T_{\text{Fluid}} + N_g k \left( \frac{T_{\text{Lift}} + T_{\text{Proj}} + T_{\text{Kin}} + T_{\text{Fluid}}}{N_p} + T_{\text{Fluid}} \right). $$

Consequently, the ideal number of parareal iterations to be less costly than a fully kinetic simulation should satisfy:

$$ T_{\text{Parareal}} \leq N_g T_{\text{Kin}} $$
<table><caption>TABLE 1. Computing architecture.</caption><tr><td># CPUs</td><td>Intel(R) Xeon(R) Platinum 8268 CPU @ 2.90GHz x2</td></tr><tr><td># cores</td><td>48</td></tr><tr><td>RAM</td><td>750 GB</td></tr></table>

or stated differently,

$$k_{\text{opt}} \le \left[ \frac{N_g T_{\text{Kin}} - T_{\text{Fluid}}}{N_g \left( \frac{T_{\text{Lift}} + T_{\text{Proj}} + T_{\text{Kin}} + T_{\text{Fluid}}}{N_p} + T_{\text{Fluid}} \right)} \right].$$

**Practical implementation.** The performance of the method and scalability of the current implementation for this paper were made in Fortran90 (gfortran 9.4.0) with an OpenMP parallelization that will be discussed in the Section 5. These tests were ran on the architecture detailed in Table 1.

In a shared memory paradigm, the parallelization of the parallel loop in algorithm 2 is achieved using a dynamic scheduling as the number of iteration is not fixed.

## 5. NUMERICAL RESULTS

In this section, we present some numerical results for simulations of the 1D/3D Vlasov-BGK equation (1) in the hydrodynamic scaling. We consider Algorithm 2 and assess its accuracy as well as its performance through several experiments. In the following, unless specified otherwise, the phase space is uniformly discretized using 200 × 32 × 32 × 32 cells, namely 200 points in position and 32 in each velocity direction. Note that the method was only implemented in the hydrodynamic scaling with an Euler solver at the fluid level and the lifting does not consider any perturbations and then any higher order models. Extensions to the cases mentioned in Section 3 will be implemented in future works.

### 5.1. Test 1: Sod shock tube.
We first consider a Riemann problem where the initial data is given by:

$$f_0(x, v) = \mathcal{M}_{\rho(x), u(x), \theta(x)}(v), \quad x \in [0, 2], \quad v \in [-8, 8],$$

where

$$(\rho(x), u(x), \theta(x)) = \begin{cases} (1, 0, 0, 0, 1) & \text{if } x < 1, \\ (0.125, 0, 0, 0, 0.8) & \text{if } x \ge 1. \end{cases}$$

The exterior force is set to 0, and the regime is rarefied: $\epsilon = 10^{-2}$ at the kinetic level. Boundary conditions are of the absorbing type. We set a simulation interval $[0, T]$ with $T = 0.5$ and discretize it uniformly using $N_g = 200$ points for the coarse grid and $N_f = 800$ points for the fine grid. As a stopping criterion, we set a maximum of $K = 80$ iterations or when then consecutive error (2.2) is smaller than $10^{-8}$. Note that this threshold is much smaller than both the fine and coarse time increments.

We present in Figure 1 snapshots of the moments obtained with Algorithm 2 as well as the ones obtained through fully kinetic and fully fluid simulations. We can observe a very good agreement between the kinetic and parareal densities, even if the results of the fully fluid equation (that we recall is used as coarse propagator) are far from exact. Note that we also notice that higher order moments are less well captured, and this observation also holds as we look for larger times. Nevertheless, considering we used the simplest fluid model and reconstruction, the results are promising.

**Convergence.** Figure 2 illustrates the convergence of the algorithm, by plotting in loglog scale the successive errors (2.2) at each parareal iteration. A first observation is that the algorithm indeed converges well. In particular, the error decreases exponentially fast with the number of parareal iterations, which is encouraging in the sense that one may need only a few iterations to reach a desired accuracy.
<figure><img src="image_4.png" /><figcaption>FIGURE 1. Test 1 - Sod shock tube, ε = 10⁻²: Snapshots of the density (Top), x mean velocity (Middle) and Temperature (Bottom) at times Tⁿ = 0.05 (Left), 0.15 (Middle) and 0.25 (Right).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 2. Test 1 - Sod shock tube, ε = 10⁻²: Convergence of the successive errors.</figcaption></figure>

5.2. Test 2: Blast waves. We now consider a second Riemann problem consisting of two hot, energetic shock waves converging asymmetrically towards a same point, namely a blast wave-type collision. The initial data is at a local equilibrium in velocity.

$$f_0(x, v) = \mathcal{M}_{\rho(x), u(x), \theta(x)}(v), \quad x \in [0, 2], \quad v \in [-8, 8],$$
<figure><img src="image_6.png" /><figcaption>FIGURE 3. Test 2 - Blast waves, ε = 10⁻²: Snapshots of the density (Top), x mean velocity (Middle) and Temperature (Bottom) at times Tⁿ = 0.1 (Left), 0.23 (Middle) and 0.3 (Right).</figcaption></figure>

with the moments defined by:

$$ (\rho(x), u(x), \theta(x)) = \begin{cases} (1, 1, 0, 0, 2) & \text{if } x < 0.4, \\ (1, 0, 0, 0, 0.25) & \text{if } 0.4 \le x < 1.6, \\ (1, -1, 0, 0, 2) & \text{if } x \ge 1.6. \end{cases} $$

We consider the same setting as in Test 1, and we fix the number of parareal iterations to K = 10. We present in Figure 3 snapshots of the moments. Note that with only 10 iterations, the consecutive error is only of the order of 10⁻³. We observe that the solution has not yet converged towards the kinetic one, and this is even more striking for large time. Nevertheless, this behavior is expected as the parareal algorithm first corrects for early times and then propagates this correction.

In Figure 4, we investigate the main source of error by plotting the point-wise difference between the kinetic and parareal moments. According to Figure 3 we observe that the deviation between kinetic and parareal solvers is essentially localized near shocks, where the regularity of the solution is low. Again, such behavior is not surprising as it was already observed in previous works [20, 14, 43, 2, 16].

We also report in Figure 5 the convergence of the algorithm where we observe the same exponential decay with respect of the number of parareal iteration, as was noticed in Test 1.

5.3. Test 3: Exterior forces. We conclude this numerical investigation of the multiscale PinT method by considering the case of a smooth initial data along with a non-zero exterior force. Moreover, we consider this time an initial data far from the Maxwellian equilibrium
<figure><img src="image_7.png" /><figcaption>FIGURE 4. Test 2 - Blast waves, $\epsilon = 10^{-2}$: Snapshots of the pointwise difference on the density (Top), $x$ mean velocity (Middle) and Temperature (Bottom) at times $T^n = 0.1$ (Left), 0.23 (Middle) and 0.3 (Right).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 5. Test 2 - Blast waves, $\epsilon = 10^{-2}$: Convergence of the successive errors.</figcaption></figure>

given by two opposite beams with constant densities:

$$f_0(x, v) = \mathcal{M}_{\rho_1(x), u_1(x), \theta_1(x)}(v) + \mathcal{M}_{\rho_2(x), u_2(x), \theta_2(x)}(v), \quad x \in [0, 2], \quad v \in [-8, 8],$$

where the two sets of moments are given by:

$$(\rho_1(x), u_1(x), \theta_1(x)) = (1, 1, 0, 0, 1)$$

and

$$(\rho_2(x), u_2(x), \theta_2(x)) = (1, -1, 0, 0, 1).$$
<figure><img src="image_9.png" /><figcaption>FIGURE 6. Test 3 - Exterior forces: Plot of the external force $E(x)$.</figcaption></figure>

<table><caption>TABLE 2. Test 3 - Exterior forces, $\epsilon = 10^{-5}$:</caption><thead><tr><th># Threads</th><th>1</th><th>2</th><th>4</th><th>8</th><th>16</th><th>32</th></tr></thead><tbody><tr><th>Runtime (s)</th><td>278.7</td><td>153.5</td><td>84.3</td><td>48.9</td><td>30.6</td><td>22.9</td></tr></tbody></table>

We set the exterior forces as

$$E(x) = -5x^4(x - 2)^4(x - 1).$$

This exterior force is illustrated in Figure 6 and will concentrate the mass towards the center of the domain $x = 1$. A natural follow up to this test will be to consider a coupling with Poisson's equation.

The discretization is chosen as $100 \times 256 \times 16 \times 16$ cells, namely 100 points in position, 256 points in the $v_x$ directions and 32 in the $v_y$ and $v_z$ directions. The specific refinement in $v_x$ is only meant to reduce the consistency error of our kinetic solver as there now is transport in this particular direction. Note that the number of cells in the other directions is unchanged compared to the previous test cases. In addition, we assume periodic boundary conditions and a Knudsen number of $10^{-5}$.

The resulting moments are presented in Figure 7. We can observe that the mass indeed concentrates towards the center of the domain because of the exterior force. In particular, the solution being smooth, we obtain a very good agreement between kinetic and parareal moments. In addition, we also illustrate the convergence of the algorithm in the next Paragraph along with the performance of the method.

**Performance.** We conclude this test by first investigating the parallel efficiency of our implementation. The scalability is assessed by measuring the time to perform the first, most expensive parareal iteration. The results obtained on the architecture presented in Table 1 are illustrated in Table 2. The current implementation scales as around 0.737 with the number of threads which is very good considering the parallelization method used. A refinement of this implementation in MPI is discussed in Appendix A.

Finally, we report in Table 3 the successive error, the runtimes and the speedup obtained with our new method. Firstly, we observe a significant reduction of the computational time with a factor that naturally depends on the number of parareal iterations. For a reasonable error of $10^{-6}$, the observed speedup reaches 4.52 between the full, direct kinetic integration, and the new multiscale PinT one. We believe that this is particularly significant, specially when one considers our simplified $d_x = 1$ setting.

## 6. CONCLUSION

We presented in this work the first multiscale, parareal in time numerical solver for computing efficiently solution to collisional, Vlasov-type kinetic equation. We introduced our method, analysed some of its numerical properties, proposed a first implementation in a shared memory framework and conducted a thorough testing of this new method. We
<figure><img src="image_10.png" /><figcaption>FIGURE 7. Test 3 - Exterior forces, $\epsilon = 10^{-5}$: Snapshots of the density (Top), $x$ mean velocity (Middle) and Temperature (Bottom) at times $T^n = 0.1$ (Left), 0.25 (Middle) and 0.5 (Right).</figcaption></figure>

<table><caption>TABLE 3. Test 3 - Exterior forces, $\epsilon = 10^{-5}$: Performance of Algorithm 2 on 32 threads.</caption><thead><tr><td>Parareal Iterations</td><td>Final error</td><td>Runtime (s)</td><td>Kinetic Runtime</td><td>Speedup</td></tr></thead><tbody><tr><td>3</td><td>1.387E-03</td><td>71.7</td><td>672.8</td><td>9.38</td></tr><tr><td>6</td><td>7.782E-06</td><td>153.0</td><td>672.8</td><td>4.52</td></tr><tr><td>9</td><td>2.421E-08</td><td>218.4</td><td>672.8</td><td>3.08</td></tr></tbody></table>

exhibited some very serious speedup (up to 10 times faster than a full kinetic run) using a regular, not state-of-the-art computer. We finally presented a possible distributed memory version of our algorithm that we hope to test in bigger computational clusters in future works.

In addition, it is worth mentioning that we expect even larger computational gains in more complex cases such as higher dimensions in positions and Boltzmann-type collision operators. Moreover, as it was for example observed in [40], a speedup can also be more significant for large time simulations, which have not been considered in this work.

On a final note, more sophisticated versions of the parareal algorithm have been recently constructed with improvement of performance in mind and we refer to [40, 42] for more details that we also wish to consider in future works.

ACKNOWLEDGEMENT

The authors would like to warmly thank M. Bessemoulin-Chattard for the careful reading and feedback on early versions of the manuscript that happened at the end of T. Laidin's
PhD. T. Laidin has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 865711). T. Rey received funding from the European Union's Horizon Europe research and innovation program under the Marie Sklodowska-Curie Doctoral Network DataHyKing (Grant No. 101072546), and is supported by the French government, through the UniCA<sub>JEDI</sub> Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01.

## APPENDIX A. MPI IMPLEMENTATION OF THE MULTISCALE PARAREAL METHOD

This appendix is dedicated to proposing an MPI implementation of Algorithm 2. We present the procedure to be implemented and discuss the challenges it poses.

**Workload distribution.** Let us consider the parallelization of Algorithm 2 over $N_p \in \mathbb{N}$ processors with distributed memory. We recall that we aim at computing the expensive kinetic propagations in parallel to reduce the simulation time.

A first observation is that after each parareal iteration, the actual work to be done diminishes. Therefore, one must adjust the load distribution at each iteration. Let us now consider a discretization of the position variable using $N_x$ points and $N_g$ points for the coarse time discretization. From a practical point of view, the memory allocation of the jumps arrays ($\Delta^n$ in Algorithm 2) should be done locally on each processor, and only once, to avoid multiple allocation/deallocation overheads. Indeed, since they depend on the size of the spatial discretization, their memory footprint is potentially large. Therefore, the master processor shall be allocated with the full time discretization which is of size $5N_xN_g$, where 5 corresponds to the number of moments. The remaining processors are allocated their maximum workload, namely their workload at the first parareal iteration. This quantity, denoted by **chunk** in the following, can be defined using integer arithmetic to evenly distribute the time iterations.

The workload at parareal iteration $k$, the starting and ending time indices for each processor are then computed as

$$ 
\begin{aligned}
work &= N_g - k; \\
chunk &= \lfloor work / N_p \rfloor; \\
remainder &= \text{mod}(work, N_p); \\
start &= \text{rank} * chunk + \min(remainder, \text{rank}) + 1; \\
end &= (\text{rank} + 1) * chunk + \min(remainder, \text{rank} + 1); \\
myChunk &= end - start.
\end{aligned}
$$

At later iterations, when each process has less fine propagations to deal with, one can adjust using array indexing instead of dealing with multiple allocation/deallocation. Our strategy is presented in Algorithm 3 and we shall now discuss the communication aspect of our approach.

**Communications.** The main challenge to implement this strategy will be to efficiently send time-chunks of data that are of size $5N_x \times \text{myChunk}$ which are potentially large. It is however important to note that since the workload (in time) diminishes with parareal iterations, an empirical threshold for the number of processors to use should be considered to keep the amount of communications in check. Our approach is to gather all the local-in-time jumps on the master processor (step 14) and update the solution. Then, the full updated solution of size $5N_xN_g$ is broadcast to each processor (step 20) to perform its work for the next parareal iteration. Note that this strategy may not be the most efficient depending on the size of the problem and computing architecture. For example, we refer to [51] where message passing cost and optimization is discussed.

Finally, another computational difficulty will be the actual storage of such data that needs to be MPI friendly in order to avoid the use of many buffer arrays and have an optimized browsing of the memory.
<caption>**Algorithm 3** MPI implementation: multiscale kinetic parareal Algorithm 2.</caption>

**Require:** $U^{0,0}$
1: **for** $n = 1, \dots, N_g$ **do**
2: &nbsp;&nbsp;&nbsp;&nbsp; $U^{n,0} = \mathcal{G}(U^{n-1,0})$
3: **end for**
4: **while** $k \le K$ **or** error $\ge \text{tol}$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp; **work** = $N_g - k$
6: &nbsp;&nbsp;&nbsp;&nbsp; **chunk** = $[\text{work}/N_p]$
7: &nbsp;&nbsp;&nbsp;&nbsp; **remainder** = $\text{mod}(\text{work}, N_p)$
8: &nbsp;&nbsp;&nbsp;&nbsp; **start** = $\text{rank} * \text{chunk} + \min(\text{remainder}, \text{rank}) + 1$
9: &nbsp;&nbsp;&nbsp;&nbsp; **end** = $(\text{rank}+1) * \text{chunk} + \min(\text{remainder}, \text{rank}+1)$
10: **myChunk** = **end-start**
11: **for** $n = start, \dots, end$ **do**
12: &nbsp;&nbsp;&nbsp;&nbsp; $\Delta^n = \mathcal{P}\mathcal{F}\mathcal{L}(U^{n-1,k-1}) - \mathcal{G}(U^{n-1,k-1})$
13: **end for**
14: Gather chunks of jumps $\Delta^n$ to master process
15: **if** $\text{rank}$ is master **then**
16: &nbsp;&nbsp;&nbsp;&nbsp; **for** $n = k, \dots, N_g$ **do**
17: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $U^{n,k+1} = \mathcal{G}(U^{n-1,k}) + \Delta^n$
18: &nbsp;&nbsp;&nbsp;&nbsp; **end for**
19: **end if**
20: Broadcast updated solution $U^{n,k+1}$;
21: Compute successive error on the moments and
22: $k = k + 1$
23: **end while**

REFERENCES

[1] ABDULLE, A., E. W., ENGQUIST, B., AND VANDEN-EIJNDEN, E. The heterogeneous multiscale method. *Acta Numer.* **21** (2012), 1–87.

[2] BAL, G. On the convergence and the stability of the parareal algorithm to solve partial differential equations. In *Domain decomposition methods in science and engineering. XV. Selected papers of the 15th international conference on domain decomposition*, Berlin, Germany, July 21–25, 2003. Berlin: Springer, 2005, pp. 425–432.

[3] BARDOS, C. W., GOLSE, F., AND LEVERMORE, D. M. Fluid dynamic limits of kinetic equations. i. formal derivations. *J. Statist. Phys.* **63** (1991), 323–344.

[4] BAUDRON, A.-M., LAUTARD, J.-J., MADAY, Y., AND MULA, O. The parareal in time algorithm applied to the kinetic neutron diffusion equation. In *Domain Decomposition Methods in Science and Engineering XXI*. pringer International Publishing, 2014, pp. 437–445.

[5] BHATNAGAR, P. L., GROSS, E. P., AND KROOK, M. A model for collision processes in gases. I: Small amplitude processes in charged and neutral one-component systems. *Phys. Rev., II. Ser.* **94** (1954), 511–525.

[6] BOSSUYT, I., VANDEWALLE, S., AND SAMAEY, G. Monte-carlo/moments micro-macro parareal method for unimodal and bimodal scalar mckean-vlasov sdes. Preprint arXiv:2310.11365, 2023.

[7] CORON, F., AND PERTHAME, B. Numerical passage from kinetic to fluid equations. *SIAM J. Numer. Anal.* **28**, 1 (1991), 26–42.

[8] CROUSEILLES, N., DEGOND, P., AND LEMOU, M. A hybrid kinetic/fluid model for solving the gas dynamics Boltzmann-BGK equation. *J. Comput. Phys.* **199**, 2 (2004), 776–808.

[9] CROUSEILLES, N., AND LEMOU, M. An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits. *Kinet. Relat. Mod.* **4**, 2 (2011), 441–477.

[10] DEGOND, P., DIMARCO, G., AND MIEUSSENS, L. A moving interface method for dynamic kinetic-fluid coupling. *J. Comput. Phys.* **227**, 2 (2007), 1176–1208.

[11] DUARTE, M., MASSOT, M., AND DESCOMBES, S. Parareal operator splitting techniques for multi-scale reaction waves: numerical analysis and strategies. *ESAIM, Math. Model. Numer. Anal.* **45**, 5 (2011), 825–852.

[12] E, W. *Principles of multiscale modeling*. Cambridge: Cambridge University Press, 2011.
[13] E. W., AND ENGQUIST, B. The heterogeneous multiscale methods. Commun. Math. Sci. 1, 1 (2003), 87–132.

[14] EGHBAL, A., GERBER, A. G., AND AUBANEL, E. Acceleration of unsteady hydrodynamic simulations using the parareal algorithm. J. Comput. Sci. 19 (2017), 57–76.

[15] EINKEMMER, L., HU, J., AND WANG, Y. An asymptotic-preserving dynamical low-rank method for the multi-scale multi-dimensional linear transport equation. J. Comput. Phys. 439 (2021), 21. Id/No 110353.

[16] FARHAT, C., AND CHANDESRIS, M. Time-decomposed parallel time-integrators: Theory and feasibility studies for fluid, structure, and fluid-structure applications. Int. J. Numer. Methods Eng. 58, 9 (2003), 1397–1434.

[17] FILBET, F., AND JIN, S. An asymptotic preserving scheme for the ES-BGK model of the Boltzmann equation. J. Sci. Comput. 46, 2 (2011), 204–224.

[18] FILBET, F., AND REY, T. A hierarchy of hybrid numerical methods for multiscale kinetic equations. SIAM J. Sci. Comput. 37, 3 (2015), a1218–a1247.

[19] FISCHER, P. F., HECHT, F., AND MADAY, Y. A parareal in time semi-implicit approximation of the Navier-Stokes equations. In Domain decomposition methods in science and engineering. XV. Selected papers of the 15th international conference on domain decomposition, Berlin, Germany, July 21–25, 2003. Berlin: Springer, 2005, pp. 433–440.

[20] GANDER, M. J. Analysis of the parareal algorithm applied to hyperbolic problems using characteristics. Bol. Soc. Esp. Mat. Apl., SèMA 42 (2008), 21–35.

[21] GANDER, M. J., AND HAIRER, E. Nonlinear convergence analysis for the parareal algorithm. In Domain decomposition methods in science and engineering XVII. Selected papers based on the presentations at the 17th international conference on domain decomposition methods, St. Wolfgang/Strobl, Austria, July 3–7, 2006. Berlin: Springer, 2008, pp. 45–56.

[22] GANDER, M. J., AND HAIRER, E. Analysis for parareal algorithms applied to Hamiltonian differential equations. J. Comput. Appl. Math. 259 (2014), 2–13.

[23] GANDER, M. J., AND LUNET, T. Time Parallel Time Integration. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2024.

[24] GANDER, M. J., AND VANDEWALLE, S. Analysis of the parareal time-parallel time-integration method. SIAM J. Sci. Comput. 29, 2 (2007), 556–578.

[25] GOLSE, F., AND SAINT-RAYMOND, L. The Navier-Stokes limit for the Boltzmann equation. C. R. Acad. Sci., Paris, Sér. I, Math. 333, 9 (2001), 897–902.

[26] GRIGORI, L., HIRSTOAGA, S. A., NGUYEN, V.-T., AND SALOMON, J. Reduced model-based parareal simulations of oscillatory singularly perturbed ordinary differential equations. J. Comput. Phys. 436 (2021), 18. Id/No 110282.

[27] GRIGORI, L., HIRSTOAGA, S. A., AND SALOMON, J. A parareal algorithm for a highly oscillating Vlasov-Poisson system with reduced models for the coarse solving. Comput. Math. Appl. 130 (2023), 137–148.

[28] JIN, S. Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations. SIAM J. Sci. Comput. 21, 2 (1999), 441–454.

[29] JÜNGEL, A. Transport Equations for Semiconductors, vol. 773 of Lecture Notes in Physics. Springer Berlin Heidelberg, Berlin, Heidelberg, 2009.

[30] KLAR, A. Asymptotic-induced domain decomposition methods for kinetic and drift diffusion semiconductor equations. SIAM J. Sci. Comput. 19, 6 (1998), 2032–2050.

[31] KOCH, O., AND LUBICH, C. Dynamical low-rank approximation. SIAM Journal on Matrix Analysis and Applications 29, 2 (2007), 434–454.

[32] KOLOBOV, V. I., ARSLANBEKOV, R. R., ARISTOV, V. V., FROLOVA, A. A., AND ZABELOK, S. A. Unified solver for rarefied and continuum flows with adaptive mesh and algorithm refinement. J. Comput. Phys. 223, 2 (2007), 589–608.

[33] LAIDIN, T. Hybrid kinetic/fluid numerical method for the Vlasov-BGK equation in the diffusive scaling. Kinet. Relat. Mod. (2023).

[34] LEGOLL, F., LELIÈVRE, T., AND SAMAEY, G. A micro-macro parareal algorithm: application to singularly perturbed ordinary differential equations. SIAM J. Sci. Comput. 35, 4 (2013), a1951–a1986.

[35] LEMOU, M. Relaxed micro-macro schemes for kinetic equations. C. R., Math., Acad. Sci. Paris 348, 7-8 (2010), 455–460.

[36] LEMOU, M., AND MIEUSSENS, L. A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit. SIAM J. Sci. Comput. 31, 1 (2008), 334–368.

[37] LEVEQUE, R. J. Finite volume methods for hyperbolic problems. Camb. Texts Appl. Math. Cambridge: Cambridge University Press, 2002.

[38] LEVERMORE, C. D., MOROKOFF, W. J., AND NADIGA, B. T. Moment realizability and the validity of the Navier–Stokes equations for rarefied gas dynamics. Physics of Fluids 10 (1998), 3214–3226.

[39] LIONS, J.-L., MADAY, Y., AND TURINICI, G. A “parareal” in time discretization of PDE’s. C. R. Acad. Sci., Paris, Sér. I, Math. 332, 7 (2001), 661–668.

[40] MADAY, Y., AND MULA, O. An adaptive parareal algorithm. J. Comput. Appl. Math. 377 (2020), 17. Id/No 112915.

[41] MADAY, Y., AND TURINICI, G. A parareal in time procedure for the control of partial differential equations. C. R., Math., Acad. Sci. Paris 335, 4 (2002), 387–392.

[42] NGUYEN, V.-T. Acceleration techniques of the Parareal algorithm for solving some differential equations. Theses, Sorbonne Université, May 2022.
[43] NIELSEN, A. S., BRUNNER, G., AND HESTHAVEN, J. S. Communication-aware adaptive parareal with application to a nonlinear hyperbolic system of partial differential equations. *J. Comput. Phys.* **371** (2018), 483–505.

[44] PARESCHI, L., AND RUSSO, G. Implicit-explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation. *J. Sci. Comput.* **25**, 1-2 (2005), 129–155.

[45] POUPAUD, F. Diffusion approximation of the linear semiconductor Boltzmann equation: Analysis of boundary layers. *Asymptotic Anal.* **4**, 4 (1991), 293–317.

[46] POUPAUD, F., AND SOLER, J. Parabolic limit and stability of the Vlasov-Fokker-Planck system. *Math. Models Methods Appl. Sci.* **10**, 7 (2000), 1027–1045.

[47] SAMADDAR, D., NEWMAN, D. E., AND SÁNCHEZ, R. Parallelization in time of numerical simulations of fully-developed plasma turbulence using the parareal algorithm. *J. Comput. Phys.* **229**, 18 (2010), 6558–6573.

[48] SAMAEY, G., AND SLAWIG, T. A micro/macro parallel-in-time (parareal) algorithm applied to a climate model with discontinuous non-monotone coefficients and oscillatory forcing. Preprint arXiv:1806.04442, 2023.

[49] STEINER, J., RUPRECHT, D., SPECK, R., AND KRAUSE, R. Convergence of parareal for the Navier-Stokes equations depending on the Reynolds number. In *Numerical mathematics and advanced applications – ENUMATH 2013. Proceedings of ENUMATH 2013, the 10th European conference on numerical mathematics and advanced applications, Lausanne, Switzerland, August 26–30, 2013*. Cham: Springer, 2015, pp. 195–202.

[50] STRUCHTRUP, H. Macroscopic transport equations for rarefied gas flows. *Approximation methods in kinetic theory*. Interact. Mech. Math. Berlin: Springer, 2005.

[51] THAKUR, R., RABENSEIFNER, R., AND GROPP, W. Optimization of collective communication operations in mpich. *Int. J. High Perform. Comput. Appl.* **19**, 1 (2005), 49–66.

[52] TIWARI, S. Coupling of the Boltzmann and Euler equations with automatic domain decomposition. *J. Comput. Phys.* **144**, 2 (1998), 710–726.

[53] TIWARI, S., AND KLAR, A. An adaptive domain decomposition procedure for Boltzmann and Euler equations. *J. Comput. Appl. Math.* **90**, 2 (1998), 223–237.

[54] TORO, E. F. Riemann solvers and numerical methods for fluid dynamics. *A practical introduction*, 3rd ed. ed. Berlin: Springer, 2009.

[55] VAN LEER, B. Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method. *J. Comput. Phys.* **32**, 1 (1979), 101–136.