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

HAL
open science

# An Unstructured Forward-Backward Lagrangian Scheme for Transport Problems

Martin Campos Pinto

► To cite this version:

Martin Campos Pinto. An Unstructured Forward-Backward Lagrangian Scheme for Transport Problems. Finite Volumes for Complex Applications VIII - Hyperbolic, Elliptic and Parabolic Problems, 2017, Lille, France. pp.265-273, 10.1007/978-3-319-57394-6_29. hal-04041196

HAL Id: hal-04041196

https://hal.science/hal-04041196

Submitted on 23 Mar 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.

Public Domain
# An Unstructured Forward-Backward Lagrangian Scheme for Transport Problems

Martin Campos Pinto

**Abstract** In a recent work [3], a novel method has been proposed and analyzed to reconstruct accurate backward transport flows, based on point markers pushed forward. When used in conjunction with a reliable particle code, this approach provides a simple tool to improve the accuracy of density approximations. In this article we report on an extension of the method to unstructured set of markers, which are generic in particle codes. The resulting approximations have the same order of convergence, both in the a priori estimates and in the numerical simulations.

**Key words:** Particle methods, Semi-Lagrangian methods, Transport equations, A priori error estimates, Remapped particle methods

MSC (2010): 76M28, 35F10, 65M12

## 1 Introduction

Consider a transport equation

$$ \partial_t f(t,x) + u(t,x) \cdot \nabla f(t,x) = 0, \quad t \in [0,T], \quad x \in \mathbb{R}^d \quad (1) $$

associated with an initial data $f^0: \mathbb{R}^d \to \mathbb{R}$ and a velocity field $u: [0,T] \times \mathbb{R}^d \to \mathbb{R}^d$. If $u$ is smooth, e.g. $L^\infty(0,T; W^{1;\infty}(\mathbb{R}^d))$ [13], we can define characteristic trajectories $X(t) = X(t;s,x)$ solutions to the ODEs $X'(t) = u(t,X(t))$, $X(s) = x$ on $[0,T]$, for all $x \in \mathbb{R}^d$ and $s \in [0,T]$. The corresponding flow $F_{s,t}: x \mapsto X(t)$ is then invertible and satisfies $(F_{s,t})^{-1} = F_{t,s}$. In particular, the transported density reads

<footnote>Martin Campos Pinto</footnote>

<footnote>CNRS, Sorbonne Universités, UPMC Univ Paris 06, UMR 7598, Laboratoire Jacques-Louis Lions, 4, place Jussieu 75005, Paris, France</footnote>

<footnote>e-mail: campos@ljll.math.upmc.fr</footnote>
$$f(t,x) = f^0((F_{0,t})^{-1}(x)) \quad \text{for } t \in [0,T], x \in \mathbb{R}^d. \quad (2)$$

In general $u$ depends on $f$ through some self-consistent coupling. Following [3] we assume that we are given an accurate particle solver that (i) pushes forward arbitrary sets of markers along the characteristic flow, and (ii) computes reliable approximations to the velocity field $u$, at least for moderate simulation times $T' < T$. Using this solver as a black box we may then consider that $u$ is given and that we can apply the exact forward flow $F_{ex}^{n,n+1}$ on small time steps of size $\Delta t$, on a given set of point markers.

The method that we describe below follows a series of works [10, 5, 8, 1, 2] where accurate approximations of transported densities are obtained through enhanced representations of the transport flow. In its most recent version [3], it implements the fundamental idea that local descriptions of the characteristic flow can be computed using a rather inexpensive method and then exploited to accurately reconstruct the transported density [6]. In practice the method studied in [3] combines key tools from the usual forward and backward Lagrangian methods. It consists of

* pushing forward given markers along the characteristic trajectories, like in a standard particle method, and

* representing the density on a grid at given time steps, like in a backward semi-Lagrangian method.

The crux of the method is then to use the markers pushed forward to approximate the backward flow between two time steps. The approximated density is then transported as in a standard Backward Semi-Lagrangian (BSL) method [14]. As explained in [2, 3], the strength of this approach over a standard particle method with smooth remappings (interpolations) is a lower diffusivity and higher convergence rate, and compared to the BSL method it has the advantage of avoiding a backward time integration of the trajectories. Owing to its hybrid nature we call it a *Forward-Backward Lagrangian* method.

# 2 The Forward-Backward Lagrangian method

## 2.1 Backward flow reconstruction

The method relies on local approximations of the backward flow that are valid close to the marker positions $x_k^n$. For simplicity we restrict our presentation to first order flow approximations. Following [2] we define

$$B_{h,k}^{0,n} : x \mapsto x_k^0 + D_k^n (x - x_k^n) \quad (3)$$

with $D_k^n$ a $d \times d$ matrix approximating the backward Jacobian $J_{B_{ex}^{0,n}}(x_k^n)$ of the backward flow, that can be computed from the position of the neighboring markers.
To compute a global approximation to the backward flow we subdivide the computational domain $\Omega$ in cartesian cells of size $h$, with centers denoted $\xi_i = ih, i \in \mathbb{Z}^d$, to avoid a confusion with the particle positions. To any $\xi_i$ we associate a nearby marker, e.g. the closest one,

$$k^*(n, i) := \argmin_{k \in \mathbb{Z}^d} \|x_k^n - \xi_i\|_{\infty}$$

and its corresponding backward flow (3). The global approximation to $B_{\text{ex}}^{0,n}$ is then obtained by smoothly patching these local approximations. Given a partition of unity $\sum_{i \in \mathbb{Z}^d} S(x - i) = 1$ involving a compactly supported, non-negative shape function $S$ (e.g., a B-spline), we set

$$B_h^{0,n}(x) := \sum_{i \in \mathbb{Z}^d} B_{h,k^*(n,i)}^{0,n}(x) S\left(\frac{x - \xi_i}{h}\right). \quad (4)$$

## 2.2 Remapped FBL method

Used in conjunction with a standard particle code, the above technique can be used in several ways to derive numerical schemes that improve the accuracy of the particle approximations.

For instance, the density can be approximated at any time step $n$ by using the approximated flow (4) in the Lagrangian formula

$$f_h^{n, \text{fbl}}(x) := f^0(B_h^{0,n}(x)), \quad x \in \mathbb{R}^d. \quad (5)$$

Assuming that the underlying particle code pushes the markers along accurate trajectories, this reconstruction will be accurate as long as the associated characteristic flow remains smooth.

In many cases however, the regularity of the flow deteriorates over time and so does the accuracy of its approximations. To reduce this effect a simple method then consists of restarting the transport problem from time to time, namely before the approximated flow becomes too inaccurate. In the literature these restart time steps are often called *remappings*, and they essentially consist of re-initializing both the approximated density and the flow markers. After a restart indeed we must solve a new transport problem, where the characteristic flow has been reset to the identity mapping of $\mathbb{R}^d$. This comes at a price, which is the approximation error on the transported density. Formally the method reads as follows.

1. The two ingredients of the method are initialized: the positions of the markers $x_k^0, k \in \mathbb{Z}^d$, are computed, and the initial density $f^0$ is approximated on some grid of size $h$. This grid is a priori independent of the markers, and many methods can be used. B-spline interpolations or quasi-interpolations are simple and efficient, see e.g. [15, 2]. We denote the corresponding approximation by $f_h^0 = A_h f^0$.
2. Letting $m_0 = 0 < m_1 < m_2 < \dots$, denote the initial and subsequent remapping steps, on the r-th remapping cycle, $r = 0, 1, \dots$, we do:

a. For $n = m_r, \dots, m_{r+1} - 1$, push all the markers $x_k^{n+1} = F^n(x_k^n)$, $k \in \mathbb{Z}^d$.

b. Define the FBL approximation $f_h^{m_{r+1}, \text{fbl}} := f^{m_r} \circ B_h^{m_r, m_{r+1}}$ to $f(m_{r+1}\Delta t)$.

c. Compute a new approximated density $f_h^{m_{r+1}} := A_h f_h^{m_{r+1}, \text{fbl}}$ for the next cycle.

d. Re-initialize the markers to prepare the local flow approximations (3) between the present remapping time $m_{r+1}\Delta t$ and the future times $n\Delta t$.

To determine the method it thus remains to specify how the matrices $D_{h,k}^n$ involved in the local flows (3) are computed from the markers positions, and how the latter must be initialized (and re-initialized) so that these matrices approximate well the Jacobian matrices of the backward flow.

# 3 Flow reconstructions with structured or unstructured markers

## 3.1 A method using structured markers

In the structured version proposed and studied in [3] following [2], the flow markers are initialized on a cartesian grid,

$$x_k^0 = hk, \quad k \in \mathbb{Z}^d. \qquad (6)$$

After pushing forward these markers over $n$ time steps, one computes the deformation matrix $D_k^n$ approximating the Jacobian matrix of the backward flow at the particle position $x_k^n$, namely

$$J_{B_{\text{ex}}^{0,n}}(x_k^n) = (\partial_j(B_{\text{ex}}^{0,n})_i(x_k^n))_{1 \le i,j \le d},$$

as follows. First one approximates the derivatives of the forward flow $F_{\text{ex}}^{0,n}$ with finite differences involving the current particle positions $x_k^n = F_{\text{ex}}^{0,n}(x_k^0)$. With a centered formula we define

$$J_k^n := \left( \frac{(x_{k+e_j}^n - x_{k-e_j}^n)_i}{2h} \right)_{1 \le i,j \le d} \approx J_{F_{\text{ex}}^{0,n}}(x_k^0) \qquad (7)$$

and using the relation $J_{B_{\text{ex}}^{0,n}}(x_k^n)J_{F_{\text{ex}}^{0,n}}(x_k^0) = I_d$ which follows by differentiating the identity $x = B_{\text{ex}}^{0,n}(F_{\text{ex}}^{0,n}(x))$ at $x_k^0$, we approximate $J_{B_{\text{ex}}^{0,n}}(x_k^n)$ with

$$D_k^n := (J_k^n)^{-1}. \qquad (8)$$
Re-initializing the markers on some time step $n$ then simply consists by using a new set of markers located on the cartesian grid (6). We note that by doing this one forgets the previous marker positions $x_k^n$.

## 3.2 A new method using unstructured markers

Our extension of the above centered finite difference formulas to unstructured set of markers relies on the following notion of admissible simplices and parallelotopes.

**Definition 1.** Let $\alpha > 0$. An ordered simplex $(x_0, \dots, x_d)$ of $\mathbb{R}^d$ is called *admissible* if the unit vectors

$$e_i = \frac{c - x_i}{\|c - x_i\|_2}, \quad i = 0, \dots, d-1, \qquad (9)$$

defined with $c = \frac{1}{2}(x_d + x_0)$ the center of the last edge, form a matrix satisfying

$$|\det(e_0, \dots, e_{d-1})| \ge \alpha. \qquad (10)$$

The associated *admissible parallelotope* is obtained by adding the vertices

$$x_{i+d} = x_0 + x_d - x_i, \quad i = 1, \dots, d-1, \qquad (11)$$

so that $(x_0, x_i, x_d, x_{i+d})$ forms a parallelogram with center $c$.

On every remapping step (including the initial step), the unstructured markers are then *prepared* in two steps. Below we consider a fixed value for $\alpha < 1$.

1. Inside every cell $C_j$ of some cartesian mesh of resolution $\mathcal{O}(h)$, determine whether there exists an admissible simplex of markers $(x_0, \dots, x_d)$. If not, insert new markers to form one (e.g. with a random algorithm).

2. Add $d$ auxiliary markers corresponding to the $d-1$ remaining vertices (11) and center $x_{2d} = \frac{1}{2}(x_0 + x_d)$ of the associated parallelotope $X_j = (x_0, \dots, x_{2d})$.

The backward flow close to a marker $x_k^n$ is then approximated as follows. Denoting by $(\bar{x}_0, \dots, \bar{x}_{2d})$ the position at time $t^n$ of the pointers $X_j = (x_0, \dots, x_{2d})$ corresponding to the cell $C_j$ containing $x_k^0$, we first approximate the derivative of the forward flow $F^{0,n}$ along each unit vector $e_i = (x_{i+d} - x_i)/\|x_{i+d} - x_i\|_2$, see (9), by

$$\delta_i := \frac{\bar{x}_{i+d} - \bar{x}_i}{\|x_{i+d} - x_i\|_2} = \frac{F^{0,n}(x_{i+d}) - F^{0,n}(x_i)}{\|x_{i+d} - x_i\|_2}, \quad i = 0, \dots, d-1. \qquad (12)$$

This is a second order formula for $J_{F_{\text{ex}}^{0,n}}(x_{2d})e_i$, where $x_{2d}$ is the center of the paralletope. The forward Jacobian matrix at $x_{2d}$ is then approximated by

$$J_k^n := (\delta_0, \dots, \delta_{d-1})E^{-1}$$

where $E$ is the matrix $(e_0, \dots, e_{d-1})$, and the backward Jacobian matrix by
$$D_k^n := (J_k^n)^{-1} = E(\delta_0, \dots, \delta_{d-1})^{-1}.$$

We observe that this matrix is actually determined by the parallelotope $X_j$ and hence only depends on $j$. Moreover, since the center of the parallelotope is also a marker $x_{k(j)}^0 \in C_j$ that has been pushed forward, $D_k^n = D_{k(j)}^n$ approximates the backward Jacobian matrix at $x_{k(j)}^n$. The local backward flow (3) associated to any $x_k^n$ that was initially (i.e. at the last remapping step) in the cell $C_j$ is then be defined as

$$B_{h,k}^{0,n} = B_{h,k(j)}^{0,n} : x \mapsto x_{k(j)}^{0} + D_{k(j)}^{n}(x - x_{k(j)}^{n}). \quad (13)$$

From the property (10) satisfied by the admissible parallelotope we can derive a priori estimates similar to the ones of the structured case [3].

**Lemma 1.** Let $x_k^n$ be a marker initially located in a cell $C_j$ associated with an admissible parallelotope $X_j$ with center $x_{k(j)}^0$. Then the approximated forward Jacobian defined above satisfies the a priori estimate

$$\|J_k^n - J_{F_{\text{ex}}}^{0,n}(x_{k(j)}^0)\|_{\infty} \le C h^q |F_{\text{ex}}^{0,n}|_{q+1}, \quad q \in \{1, 2\} \quad (14)$$

with a constant depending only on $d$ and $\alpha$. In addition, there exists $h^* > 0$ such that for all $h < h^*$, $J_k^n$ is invertible and the following a priori estimate holds

$$\|D_k^n - J_{B_{\text{ex}}^{0,n}}(x_{k(j)}^n)\|_{\infty} \le C \min_{q \in \{1,2\}} (h^q |F_{\text{ex}}^{0,n}|_{q+1}) |F_{\text{ex}}^{0,n}|_1^{2(d-1)}. \quad (15)$$

# 4 Numerical results

The efficiency of the structured method has been assessed on several transport problems in [3], see also [4] for a smooth particle approximation to Vlasov-Poisson plasmas based on similar flow reconstructions. To validate the unstructured version we use a passive transport problem of [11] which involves a swirling velocity field $u(t,x) := \cos(\pi t/T) \operatorname{curl} \phi(x)$ with $\phi(x) := -\frac{1}{\pi} \sin^2(\pi x_1) \sin^2(\pi x_2)$ and $T = 5$. The time symmetry yields a reversible problem: at $t = T/2$ the solutions reach a maximum stretching, and they revert to their initial value at $t = T$. For the initial data we consider a smooth hump centered on $\bar{x} = (0.5, 0.7)$ with approximate radius 0.2, $f^0(x) := \frac{1}{2}(1 + \operatorname{erf}(\frac{1}{3}(11 - 100\|x - \bar{x}\|_2)))$. Figure 1 shows the profile of accurate solutions at initial, half and final times. $L^2$ convergence curves achieved by several particle methods are then plotted in Figure 2. Here the particle pusher is a RK4 scheme with time step $\Delta t = T/100 = 0.05$ that has been taken small enough to have no significant effect on the final accuracy, and all the remappings are performed with cubic splines. As a reference, the top left panel shows results obtained with a standard Forward Semi-Lagrangian (FSL) scheme (i.e., a smooth particle method with periodic remappings), see e.g [12, 7, 9]. The curves in the top right panel are obtained with the structured FBL version and confirm (i) the improved accuracy of
this approach, and (ii) the need for much less remappings. In the two bottom panels we then show results corresponding to the unstructured FBL method which seem to exhibit similar accuracies. On the left the approximate flows are computed as described above, using 5 markers per cell corresponding to the vertices and centers of the admissible parallelograms. On the right a small variant is tested, where the center marker is discarded and replaced by the current marker $x_k^h$ in the backward flow (13). This represents an error of order $h$ on the initial and current positions of the center marker, but if the flow is $W^{2,\infty}$ it also corresponds to an error of order $h$ on the Jacobian matrix, hence the resulting approximated flow is again of order $h^2$ which is confirmed by the numerical convergence rate.

<figure><img src="image_2.png" /><figcaption>Fig. 1 Initial, intermediate and final profiles of the solution to the reversible test case in the text.</figcaption></figure>

# 5 Conclusion and perspectives

An unstructured version of a recent Forward-Backward particle method has been described and validated, showing second order accuracy on passive transport problems. These results represent an encouraging step towards the implementation of such reconstruction methods within standard particle codes. Further comparisons with standard advection methods should be performed to investigate the merits of this approach, and extensions to nonlinear transport problems should be addressed, including problems with discontinuous flows.

# References

1. Alard, C., Colombi, S.: A cloudy Vlasov solution. Monthly Notices of the Royal Astronomical Society **359**(1), 123–163 (2005)

2. Campos Pinto, M.: Towards smooth particle methods without smoothing. Journal of Scientific Computing (2014)

3. Campos Pinto, M., Charles, F.: From particle methods to forward-backward Lagrangian schemes. hal.archives-ouvertes.fr (2016). [hal-01385676](https://hal-ouvertes.fr/01385676)
<figure><img src="image_3.png" /><figcaption>**Fig. 2** $L^2$ convergence curves (errors at $t = T$ vs. number of particles) for the test case described in the text, using several particle methods with varying remappings periods $\Delta t_r$. Results obtained with a standard FSL scheme are shown for comparison on the top left panel and as a gray curve in the other panels. The other curves are obtained with a structured FBL method and an unstructured one using admissible simplices with parameter $\alpha = 0.5$. The lower panels show that similar results are obtained with or without using the center marker, see the text for details.</figcaption></figure>

4. Campos Pinto, M., Sonnendrücker, E., Friedman, A., Grote, D., Lund, S.: Noiseless Vlasov-Poisson simulations with linearly transformed particles. Journal of Computational Physics **275**(C), 236–256 (2014)

5. Cohen, A., Perthame, B.: Optimal Approximations of Transport Equations by Particle and Pseudoparticle Methods. SIAM Journal on Mathematical Analysis **32**(3), 616–636 (2000)

6. Colombi, S., Alard, C.: A “metric” semi-Lagrangian Vlasov-Poisson solver (2016). Submitted.

7. Cotter, C., Frank, J., Reich, S.: The remapped particle-mesh semi-Lagrangian advection scheme. Quarterly Journal of the Royal Meteorological Society **133**(622), 251–260 (2007)

8. Cottet, G.H., Koumoutsakos, P., Salihi, M.: Vortex Methods with Spatially Varying Cores. Journal of Computational Physics **162**(1), 164–185 (2000)

9. Crouseilles, N., Respaud, T., Sonnendrücker, E.: A forward semi-Lagrangian method for the numerical solution of the Vlasov equation. Computer Physics Communications **180**(10), 1730–1745 (2009)

10. Hou, T.: Convergence of a Variable Blob Vortex Method for the Euler and Navier-Stokes Equations. SIAM Journal on Numerical Analysis **27**(6), 1387–1404 (1990)

11. LeVeque, R.: High-resolution conservative algorithms for advection in incompressible flow. SIAM Journal on Numerical Analysis pp. 627–665 (1996)

12. Nair, R., Scroggs, J., Semazzi, F.: A forward-trajectory global semi-Lagrangian transport scheme. Journal of Computational Physics **190**(1), 275–294 (2003)

13. Raviart, P.A.: An analysis of particle methods. In: Numerical methods in fluid dynamics (Como, 1983), pp. 243–324. Lecture Notes in Mathematics, Berlin (1985)

14. Sonnendrücker, E., Roche, J., Bertrand, P., Ghizzo, A.: The semi-Lagrangian method for the numerical resolution of the Vlasov equation. Journal of Computational Physics **149**(2), 201–220 (1999)

15. Unser, M., Daubechies, I.: On the approximation power of convolution-based least squares versus interpolation. Signal Processing, IEEE Transactions on **45**(7), 1697–1711 (1997)