# Hybrid Kinetic/Fluid numerical method for the Vlasov-Poisson-BGK equation in the diffusive scaling

Tino Laidin, Thomas Rey

► To cite this version:

Tino Laidin, Thomas Rey. Hybrid Kinetic/Fluid numerical method for the Vlasov-Poisson-BGK equation in the diffusive scaling. FVCA 10 - 2023 - International Conference on Finite Volumes for Complex Applications X, Oct 2023, Strasbourg, France. pp.229-237, 10.1007/978-3-031-40860-1_24 . hal-04023578v2

HAL Id: hal-04023578

https://hal.science/hal-04023578v2

Submitted on 5 Apr 2023

**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>

Distributed under a Creative Commons Attribution 4.0 International License
# Hybrid Kinetic/Fluid numerical method for the Vlasov-Poisson-BGK equation in the diffusive scaling

Tino Laidin¹ and Thomas Rey¹

Univ. Lille, CNRS, Inria, UMR 8524 - Laboratoire Paul Painlevé, F-59000 Lille, France

**Abstract.** This short note presents an extension of the hybrid, model-adaptation method introduced in [T. Laidin, arXiv 2202.03696, 2022] for linear collisional kinetic equations in a diffusive scaling to the nonlinear mean-field Vlasov-Poisson-BGK model. The aim of the approach is to reduce the computational cost by taking advantage of the lower dimensionality of the asymptotic model while reducing the overall error. It relies on two criteria motivated by a perturbative approach to obtain a dynamic domain adaptation. The performance of the method and the conservation of mass are illustrated through numerical examples.

**Keywords:** Vlasov-Poisson-BGK equations; Diffusion scaling; Asymptotic preserving scheme; Micro-macro decomposition; Hybrid solver

## 1 The model

In this work, we are interested in the kinetic description of a system of particles interacting via both mean-field electromagnetic interaction and collisions. Such system can be modeled using the Vlasov-Poisson-BGK equation. The unknown is the probability distribution function $f = f(t, x, v) \in \mathbb{R}^+$ solution to:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial}{\partial t} f^{\epsilon} + \frac{v_x}{\epsilon} \partial_x f^{\epsilon} + \frac{E^{\epsilon}}{\epsilon} \partial_{v_x} f^{\epsilon} = \frac{1}{\epsilon^2} \mathcal{Q}(f^{\epsilon}), \\ f(0, x, v) = f_0(x, v), \end{array} \right. \qquad (1) $$

where $(t, x, v) \in \mathbb{R}^+ \times [0, x_*] \times \mathbb{R}^{d_v}$ with periodic boundary conditions in the physical space and $v = (v_x, v_y, v_z)$. The long range interactions are modeled via the self-consistent electrical field $E = E(t, x)$ solution to the Poisson equation

$$ \partial_x E^{\epsilon} = \rho^{\epsilon} - \bar{\rho} \quad \text{with} \quad \bar{\rho} = \iint f_0 \, dx \, dv. \qquad (2) $$

The short-range interactions between particles are taken into account through the linear BGK-like relaxation collision operator

$$ \mathcal{Q}(f^{\epsilon})(t, x, v) = \rho^{\epsilon}(t, x) \mathcal{M}(v) - f^{\epsilon}(t, x, v), \quad \forall (t, x, v) \in \mathbb{R}^{+} \times [0, x_*] \times \mathbb{R}^{d_v}, $$

where the global Maxwellian and local density are respectively defined as

$$ \mathcal{M}(v) = \frac{e^{-|v|^2/2}}{(2\pi)^{d_v/2}}, \quad \rho^{\epsilon}(t, x) = \int f^{\epsilon}(t, x, v) \, dv =: \langle f^{\epsilon}(t, x, \cdot) \rangle. $$

In (1), the scaling parameter $\epsilon$ is related to the Knudsen number: the ratio between the mean free path of the particles and the length scale of observation. It is now
well known [4] that the limit case $\epsilon = 0$ is described by a drift-diffusion equation on the density $\rho$: when $\epsilon \to 0$, the distribution function $f^\epsilon$ converges towards a Maxwellian distribution $\rho_M$ whose density $\rho$ is solution to

$$ \begin{cases} \partial_t \rho - \partial_x J = 0, & J = \partial_x \rho - E\rho, \\ \rho(0,x) = \rho_0(x), & \forall x \in [0, x_*]. \end{cases} \qquad (P) $$

The approximation of solutions to (1)-(2) can be computationally expensive because of the high dimensionality of the kinetic phase space. Nevertheless, using the more accurate kinetic description of the system may not be necessary in the whole computing domain because of the local validity of the fluid description (the system being in a state of thermodynamic equilibrium). This latter is less accurate in describing the kinetic zones, but computationally less costly. The aim of this work is therefore to design a hybrid kinetic/fluid scheme with an automatic domain adaptation method. It relies on a robust numerical scheme for the kinetic equation, on relevant criteria to carefully determine fluid and kinetic regions and on a smart implementation.

## 2 Macroscopic models

The aim of this section is to derive a higher-order macroscopic model from which we deduce a macroscopic coupling criterion. It generalizes the approach presented in [5]. Let us now introduce the truncated Chapman-Enskog expansion of the distribution function $f^\epsilon$ at order $K \in \mathbb{N}^*$:

$$ f^{\epsilon}(t, x, v) = \rho^{\epsilon}(t, x) \mathcal{M}(v) + \sum_{k=1}^{K} \epsilon^k h^{(k)}(t, x, v). \qquad (3) $$

By inserting (3) into the original equation (1), one can identify powers of epsilon to obtain

$$ k = 0: \qquad h^{(1)} = -\mathcal{T}(\rho^{\epsilon} \mathcal{M}), \qquad (4a) $$

$$ k = 1: \qquad h^{(2)} = -\partial_t (\rho^{\epsilon} \mathcal{M}) - \mathcal{T}(h^{(1)}), \qquad (4b) $$

$$ 2 \le k \le K-1: \qquad h^{(k+1)} = -\partial_t h^{(k-1)} - \mathcal{T}(h^{(k)}), \qquad (4c) $$

where $\mathcal{T}f = v_x \partial_x f + E \partial_{v_x} f$ is the transport operator. To obtain a hierarchy of macroscopic models, one considers different truncation orders $K$, then plugs (3) into (1) and integrates in velocity. The order $K=1$ allows us to (formally) recover the asymptotic model (P).

Let us now present the idea behind the computations in a $1D/1D$ setting for the case $K=3$. Note that the same method can be applied up to the full $3D/3D$ setting. Integrating in velocity then yields

$$ \partial_t \rho^{\epsilon} + \partial_x \langle v h^{(1)} \rangle + \epsilon^2 \partial_x \langle v h^{(3)} \rangle = \mathcal{O}(\epsilon^4). \qquad (5) $$

The set of equations (4) allows us to compute the functions $h^{(k)}$, $k=1,2,3$. Using the identity $\partial_v \mathcal{M} = -v \cdot \mathcal{M}$, one has

$$ h^{(1)} = -v \cdot \mathcal{M} J^{\epsilon}, \quad J^{\epsilon} = \partial_x \rho^{\epsilon} - E^{\epsilon} \rho^{\epsilon}. $$
$$h^{(2)} = -\mathcal{M} \partial_t \rho^\epsilon + v^2 \mathcal{M} \partial_x J^\epsilon + (1-v^2) \mathcal{M} E^\epsilon J^\epsilon,$$

and

$$h^{(3)} = v \mathcal{M} J^\epsilon + v \mathcal{M} \partial_x (\partial_t \rho^\epsilon) - v^3 \mathcal{M} \partial_{xx} J^\epsilon - (v-v^3) \mathcal{M} \partial_x (E^\epsilon J) \\ - v \mathcal{M} E^\epsilon \partial_t \rho^\epsilon - (2v-v^3) \mathcal{M} E^\epsilon \partial_x J^\epsilon - (v^3-3v) (E^\epsilon)^2 J^\epsilon.$$

It remains to compute the quantities $\partial_x \langle vh^{(1)} \rangle$ and $\partial_x \langle vh^{(3)} \rangle$. By the definition of $\mathcal{M}$, one can explicitly compute its second and fourth moments: $m_2 = 1$ and $m_4 = 3$. It yields

$$\partial_x \langle vh^{(1)} \rangle = \partial_x J^\epsilon.$$

Moreover, to avoid some approximation of mixed derivatives, we observe from (5) that

$$\partial_t \rho^\epsilon = \partial_x J^\epsilon + \mathcal{O}(\epsilon^2), \qquad (6)$$

and we deduce that

$$\partial_t J^\epsilon = \partial_{xx} J^\epsilon - E^\epsilon \partial_x J^\epsilon - \rho^\epsilon \partial_t E^\epsilon + \mathcal{O}(\epsilon^2). \qquad (7)$$

Replacing the time derivatives by their approximations (6) and (7) yields

$$\partial_x \langle vh^{(3)} \rangle = - \partial_{xxx} J^\epsilon + \partial_x \rho^\epsilon \partial_t E^\epsilon - \rho^\epsilon \partial_t (\partial_x E^\epsilon) + 3 \partial_x E^\epsilon \partial_x J^\epsilon \\ + E^\epsilon \partial_{xx} J^\epsilon + 2 J^\epsilon \partial_x (\partial_x E) + \mathcal{O}(\epsilon^2).$$

Finally, using the Poisson equation (2) and rearranging the terms, we obtain a higher order macroscopic model:

$$\partial_t \rho^\epsilon - \partial_x J^\epsilon = -\epsilon^2 \mathcal{R} + \mathcal{O}(\epsilon^4), \qquad (8)$$

where the remainder $\mathcal{R}^\epsilon$ is given by

$$\mathcal{R}^\epsilon = - \partial_{xxx} J^\epsilon + E \partial_{xx} J^\epsilon + (2\rho^\epsilon - 3\bar{\rho}) \partial_x J^\epsilon + 2J^\epsilon \partial_x \rho^\epsilon - \partial_x \rho^\epsilon \partial_t E^\epsilon. \qquad (9)$$

Furthermore, we emphasize the fact that this remainder term does not depend on the velocity variable but on the macroscopic quantities $\rho$ and $E$. It quantifies very accurately the deviation from the thermodynamical equilibrium.

## 3 Numerical scheme

In this part, we briefly recall the derivation of the micro-macro model for (1) introduced in [2]. Then, we present a micro-macro finite volume scheme that enjoys the property of being Asymptotic Preserving (AP), which is a crucial point of the hybrid method we present.

Let us decompose the distribution $f^\epsilon$ as follows:

$$f^\epsilon = \rho^\epsilon \mathcal{M} + g^\epsilon. \qquad (10)$$

We introduce the orthogonal projector $\Pi$ in $L^2(dx dv \mathcal{M}^{-1})$ on $Ker(\mathcal{D})$ defined for all $f \in L^2(dx dv \mathcal{M}^{-1})$ by:

$$\Pi f = \langle f \rangle \mathcal{M}.$$
The so-called micro equation is obtained by plugging (10) into (1) and applying $(I - \Pi)$. Moreover, plugging (10) into (1) and applying the projection $\Pi$ yields the macro equation. The micro-macro model is then given by:

$$ \partial_t g^\varepsilon + \frac{1}{\varepsilon}(\mathcal{T}g^\varepsilon - \partial_x \langle v_x g^\varepsilon \rangle \mathcal{M} + v_x \mathcal{M} J^\varepsilon) = \frac{-1}{\varepsilon^2} g^\varepsilon, \quad (\text{Micro}) $$

$$ \partial_t \rho^\varepsilon + \frac{1}{\varepsilon} \partial_x \langle v_x g^\varepsilon \rangle = 0. \quad (\text{Macro}) $$

Let us now present the discretization of the (Micro)-(Macro) model. We shall adopt a finite volume approach to discretize the phase space and present only the 1D/3D case that we shall consider in our numerical simulations.

**The mesh.** The velocity domain is restricted to the bounded symmetric cube $[-v_*, v_*]^3$. We consider a Cartesian mesh of the phase space composed of $N_v = 2L$ velocity cells in each direction arranged symmetrically around $v = 0$. Let $\mathcal{J} = \{-L+1, \dots, L\}$ and $j = (j_x, j_y, j_z) \in \mathcal{J}^3$ be a multi-index. The cells of the velocity mesh are denoted by $\mathcal{V}_j$ for $j \in \mathcal{J}^3$. Each cell $\mathcal{V}_j$ has a constant volume $\Delta v^3$ and midpoint $v_j$.

The physical domain, a torus $\mathbb{T}$ of length $x_*$, is discretized into $N_x$ primal cells

$$ \mathcal{X}_i = (x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}), \quad i \in \mathcal{J} = \mathbb{Z}/N_x\mathbb{Z}, $$

of constant length $\Delta x$ and centers $x_i$. We also define dual cells $\mathcal{X}_{i+\frac{1}{2}} = (x_i, x_{i+1})$ for $i \in \mathcal{J}$, of constant length $\Delta x$ and centers $x_{i+\frac{1}{2}}$. This defines respectively the so-called primal and staggered meshes as

$$ K_{ij} = \mathcal{X}_i \times \mathcal{V}_j \text{ and } K_{i+\frac{1}{2},j} = \mathcal{X}_{i+\frac{1}{2}} \times \mathcal{V}_j, \quad \forall (i,j) \in \mathcal{J} \times \mathcal{J}^3. $$

Finally, we set a time step $\Delta t > 0$ and define $t^n = n\Delta t$ for $n \in \mathbb{N}$. The Maxwellian is discretized as a product of 1-dimensional Gaussians in such a way that it satisfies discrete counterparts of its continuous properties, namely parity, positivity, and unit mass.

The transport terms in (Micro) are approximated using a first order upwind scheme and the time derivative is dealt with using a first order exponential time integrator. Regarding the macro equation, we take advantage of approximating the perturbation on the dual cells in position for the space derivative. The usual choice [1] to obtain an AP scheme is to implicit the stiff term of (Macro). The equation on the electrical field is solved using a centered finite difference scheme. The numerical scheme reads as follows:

**Proposition 1** Let $n \in \mathbb{N}$. Let $(g_{i+\frac{1}{2},j}^{\varepsilon,n})_{ij}$ and $(\rho_i^{\varepsilon,n})_i$ be given by the following micro-macro finite volume scheme:

$$ g_{i+\frac{1}{2},j}^{\varepsilon,n+1} = g_{i+\frac{1}{2},j}^{\varepsilon,n} e^{-\Delta t/\varepsilon^2} - \varepsilon(1 - e^{-\Delta t/\varepsilon^2}) \left( \frac{T_{i+\frac{1}{2},j}^{\varepsilon,n}}{\Delta x \Delta v^3} + \xi_j \mathcal{M}_j J_{i+\frac{1}{2}}^{\varepsilon,n} \right), $$

$$ \rho_i^{\varepsilon,n+1} = \rho_i^{\varepsilon,n} - \frac{\Delta t}{\varepsilon \Delta x} \left( \langle \xi_j g_{i+\frac{1}{2}}^{\varepsilon,n+1} \rangle_{\Delta} - \langle \xi_j g_{i-\frac{1}{2}}^{\varepsilon,n+1} \rangle_{\Delta} \right), $$

$$ E_{i+\frac{1}{2}}^{\varepsilon,n} - E_{i-\frac{1}{2}}^{\varepsilon,n} = (\rho_i^{\varepsilon,n} - \bar{\rho}) \Delta x, $$
where $\xi_j = \xi_{(j_x,j_y,j_z)} = v_{j_x} \forall j \in \mathcal{J}^3$ and $T_{i+\frac{1}{2},j}^{\epsilon,n}$ is the discretization of the transport terms. Assuming some uniform bounds in $\epsilon$ on $\rho^\epsilon$ and for a fixed mesh size $\Delta x$, $\Delta v > 0$, the scheme enjoys the AP property in the diffusion limit. This property does not depend on the initial data, and the associated limit scheme reads

$$ \rho_i^{n+1} = \rho_i^n + m_2^{\Delta v} \frac{\Delta t}{\Delta x} (J_{i+\frac{1}{2}}^n - J_{i-\frac{1}{2}}^n), \\ E_{i+\frac{1}{2}}^n - E_{i-\frac{1}{2}}^n = (\rho_i^n - \bar{\rho})\Delta x, $$

with the limit flux

$$ J_{i+\frac{1}{2}}^n = \frac{\rho_{i+1}^n - \rho_i^n}{\Delta x} - E_{i+\frac{1}{2}} \rho_{i+\frac{1}{2}}^n, $$

where $m_2^{\Delta v}$ is given by $\sum_{l \in \mathcal{J}} v_l^2 M_l \Delta v_l$.

# 4 Hybrid method

Let us now introduce the main contribution of this work, namely the hybrid method between kinetic and fluid schemes. It consists in a coupled solver that is faster than a full kinetic one to solve (1) while still being accurate. Following [3,5] we first construct a hybrid kinetic/fluid solver with a dynamic domain adaptation method and present its implementation.

## 4.1 Coupling criteria

The idea of the dynamic domain adaptation method is twofold. On the one hand, the subdomains must accurately describe the state of the solution. In particular, the fluid model has to be used only where the solution is near the local velocity equilibrium. On the other hand, the method has to be dynamic, in the sense that the subdomains are adapted at each time step. To determine in which domain each cell lies, we use two criteria based on the higher order macroscopic model (8) introduced in Section 2 and the norm of the perturbation $g^\epsilon = f^\epsilon - \rho^\epsilon \mathcal{M}$. Let $\mathcal{R}_i^{\epsilon,n}$ be a discretization of the remainder (9) at time $t^n$ in cell $\mathcal{X}_i$.

The coupling procedure will unfold as follows: $\mathcal{R}_i^{\epsilon,n}$ is computed using both the kinetic density $\rho^\epsilon$ in kinetic cells and the fluid density $\rho$ in fluid cells. Cell changes happen when one of the following situations occurs, depending on a coupling threshold investigated in the previous work [5]:

- If $\mathcal{R}_i^{\epsilon,n}$ is small (w.r.t. a fixed threshold), the solution (8) is close to the limit model (P) and one must use it;

- If $g^\epsilon$ is small (w.r.t. another fixed threshold), the solution is close to the local equilibrium in velocity, and one must consider the fluid description.

If any of those two conditions are not met then, depending on its current state, the cell either stays or becomes kinetic.

**Remark 1** Note that in a kinetic cell, the criterion on the norm of $g^\epsilon$ is mandatory. Indeed, the remainder $\mathcal{R}_i^{\epsilon,n}$ could be small because of small gradients, but the perturbation large. In this situation, one does not want to change from kinetic to fluid. As an example, one could take a distribution function at equilibrium in position and far from the Maxwellian in velocity (see [3] for details).
## 4.2 Implementation

The way the method is implemented is crucial. Indeed, the main goal of the method is to avoid the update of the perturbation which is the most computationally expensive part of the code. Therefore, from an implementation point of view, $g^\varepsilon$ is not updated in fluid regions and is set to 0 only when needed. In particular, it occurs only when a fluid cell becomes kinetic and when saving data. The interface conditions between kinetic and fluid cells are treated in the same way as in [5].

**Remark 2** Note that we want to start the resolution with the approach containing the full information on the system. Hence, the domain is initialized as fully kinetic. Moreover, let us emphasize that in practice the kinetic fluxes are explicitly computed. Therefore, the hybrid setting is an explicit method.

# 5 Numerical results

Let us now present some numerical simulations with our approach.

We shall start with the validation of the AP property of the micro-macro Vlasov-Poisson-BGK solver with exponential integrator presented in Section 3. It is a combination of the methods presented in [6,2] that have never been implemented in any work, to the best of our knowledge. Figure 1 presents the time evolution of the $L^2$ norm of the energy of solutions to equations (1–2) with different values of the relaxation parameter $\varepsilon$, for the seminal weak Landau damping initial data from [2]. We observe the convergence with respect to $\varepsilon$ of the energy. The oscillations due to the Vlasov-Poisson transport term $\mathcal{T}$ occur only in the kinetic regime, when $\varepsilon$ is large. They are then damped by the linear BGK term for smaller values of $\varepsilon$, where exponential decay of the electric field occurs.

<figure><img src="image_2.png" /><figcaption>Fig. 1: Fully kinetic scheme: Landau damping. Time evolution of $||E(t)||_2$ for different values of $\varepsilon$.</figcaption></figure>

We now turn our attention to the hybrid method. Figure 2 assesses the validity of this new method by computing the time evolution of the density of a smooth
solution. The initial condition is a Maxwellian distribution $ρ(\mathcal{M})$, whose initial
space dependent density is given by $ρ(x) = 1 + 0.05 \cos(2x)$ for $x_* = 2\pi$ with a
fixed $\varepsilon = 0.1$. Although the solution is far from the fluid description because of
this mild value of $\varepsilon$, we observe an almost perfect agreement between the fully
kinetic and the hybrid kinetic-fluid solvers. We also observe the back-and-forth
phenomenon between kinetic and fluid cells, resulting in large time in a full fluid
(albeit accurate with the kinetic equation) solver for the global equilibrium. The
speed-up tables are similar to those presented in [5], where factors up to 400
have been observed between the hybrid and the fully kinetic solver.

<figure><img src="image_3.png" /><figcaption>Fig. 2: **Comparison of the solvers.** Time evolution of the space densities for $\varepsilon = 0.1$ with a smooth initial data.</figcaption></figure>

Then this domain adaptation is investigated in Figure 3. One can observe the precise domain adaptation during time. We notice in particular the quick vanishing of the kinetic cells in favor of the fluid ones (and hence a computational speedup). This adaptation phenomenon can bring mass variation as noticed in [5], but we observe that it remains very close to the machine precision.

Finally, we investigate again the case of the weak Landau damping using the new hybrid solver. We observe in Figure 4 that this new method is able to accurately capture the oscillations induced by the transport operator $\mathcal{T}$ in short time. Nevertheless, these oscillations are destroyed by the switch to a full fluid solver, which relax exponentially because of its relaxation structure.

**Acknowledgements.** The authors were partially funded by Labex CEMPI (ANR-11-LABX-0007-01) and the MSCA DN-2022 program DATAHYKING.

References

1. Bennoune, M., Lemou, M., Mieussens, L.: Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics. J. Comput. Phys. **227** (2008)
<figure><img src="image_4.png" /><figcaption>Fig. 3: **Hybrid scheme.** Left. Time evolution of the state of the cells. Right. Evolution of the mass variation (orange crosses) and mass of $g^e$ (solid blue line) for $\varepsilon = 0.1$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4: **Landau damping revisited.** Time Evolution of $||E(t)||_2$ for $\varepsilon = 0.5$.</figcaption></figure>

2. Crouseilles, N., 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), 441–477 (2011)

3. Filbet, F., Rey, T.: A hierarchy of hybrid numerical methods for multi-scale kinetic equations. SIAM J. Sci. Comput. **37**(3), A1218–A1247 (2015)

4. Goudon, T., Poupaud, F.: Approximation by Homogenization and Diffusion of Kinetic Equations. Comm. Part. Diff. Eq. (2001)

5. Laidin, T.: Hybrid Kinetic/Fluid numerical method for the Vlasov-BGK equation in the diffusive scaling (2022). ArXiv preprint 2202.03696

6. Lemou, M.: Relaxed micro–macro schemes for kinetic equations. C. R. Acad. Sci. Paris, Ser. I **348**(7-8), 455–460 (2010)