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

HAL
open science

# A random choice scheme for scalar advection

Thierry Gallouët, Olivier Hurisse, Samuel Kokh

► To cite this version:

Thierry Gallouët, Olivier Hurisse, Samuel Kokh. A random choice scheme for scalar advection. International Journal for Numerical Methods in Fluids, In press, 10.1002/fld.5218. hal-04053798

HAL Id: hal-04053798

https://hal.science/hal-04053798

Submitted on 28 Jul 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
# A Random Choice Scheme for Scalar Advection

Thierry Gallouët, Olivier Hurisse, Samuel Kokh

## Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Extension of the GRU scheme to the one-dimensional advection problem</td><td>3</td></tr><tr><td>2.1</td><td>Description of the numerical scheme</td><td>3</td></tr><tr><td>2.1.1</td><td>The prediction step</td><td>3</td></tr><tr><td>2.1.2</td><td>The projection step or Glimm Random Update</td><td>4</td></tr><tr><td>2.2</td><td>Properties of the numerical scheme</td><td>4</td></tr><tr><td>2.2.1</td><td>Transitional probabilities</td><td>4</td></tr><tr><td>2.2.2</td><td>Statistical consistency of the scheme</td><td>5</td></tr><tr><td>2.2.3</td><td>Maximum principle and preservation of the monotonicity</td><td>5</td></tr><tr><td>2.2.4</td><td>Behavior at the local extrema and on locally uniform sequences</td><td>6</td></tr><tr><td>2.2.5</td><td>A word on the convergence of the scheme</td><td>9</td></tr><tr><td>3</td><td>Numerical tests for assessing the behavior of the scheme</td><td>9</td></tr><tr><td>3.1</td><td>Approximating solutions for the linear advection problem</td><td>10</td></tr><tr><td>3.2</td><td>Approximating solutions with non-uniform velocity fields</td><td>13</td></tr><tr><td>3.3</td><td>Extension for two-dimensional problems</td><td>15</td></tr><tr><td>4</td><td>Discussion on the choice of the low-discrepancy sequences</td><td>21</td></tr><tr><td>4.1</td><td>Some quasi/pseudo-random sequences and simple metrics for assessing their quality</td><td>21</td></tr><tr><td>4.2</td><td>Comparison of the results for several families of sequences</td><td>24</td></tr><tr><td>4.3</td><td>A short analysis of the results for the Halton-Van de Corput sequences</td><td>26</td></tr><tr><td>5</td><td>Conclusion</td><td>27</td></tr><tr><td>6</td><td>Appendices</td><td>27</td></tr><tr><td>6.1</td><td>A more general result for the statistical consistency of the GRU step</td><td>27</td></tr><tr><td>6.2</td><td>A class of analytical solutions for the 1D convection problem with non-uniform velocity and density</td><td>28</td></tr><tr><td>6.3</td><td>Computing a Halton-Van der Corput sequence</td><td>29</td></tr></table>

## Abstract

This paper is dedicated to a numerical method based on a random choice as proposed in Glimm's scheme. It is applied to the problem of advection of a scalar quantity. The numerical scheme proposed here relies on a fractional step approach for which: the first step is performed using any classical finite-volume scheme, and the second step is a cell-wise update. This second step is a projection based on a random choice. The resulting scheme possesses a very low level of numerical diffusion. In order to assess the capabilities of this approach, several test cases have been investigated including convergence studies with respect to the mesh-size. The algorithm performs very well on one-dimensional and multi-dimensional problems. This algorithm is very easy to implement even for multi-processor computations.

## 1 Introduction

Scalar advection is a problem that has been widely investigated in the field of numerical simulation. Two important mathematical properties of scalar advection (whithout diffusion operator), and that are preserved at a numerical level for industrial applications, is that both local extrema and sharp fronts should be preserved. Unfortunately, for most classical numerical methods, numerical diffusion tends to smear the scalar profiles and local extrema may be altered. A large amount of methods have been proposed using low- or high-order algorithms in order to simulate such advection problems. The numerical diffusion for high-order methods is in general much lower than for first-order methods. However, one technical difficulty for the former approaches is to maintain sharp fronts in weak solutions while limiting spurious oscillations that may appear at this discontinuities. Avoiding such patterns in high-order methods requires sophisticated techniques, this
difficulty being even more important when dealing with unstructured meshes for multi-dimensional domains. Among the first-order methods with a low level of numerical diffusion and a good accuracy for weak solutions, the Glimm's scheme proposed in [8] may be put forward with regards to the present work.

In Glimm's scheme, the technique for updating cell values is completely different from that which is classically used for other finite-volume techniques. The new cell values are indeed not obtained by assembling fluxes at the boundaries of each cells, but by picking one value among the possible ones (in the sense of the system) in a random manner. This results in a very low numerical diffusion and it enhances the scheme to have a convergence rate of 1 even for weak solutions. At our knowledge, this is the only first-order finite-volume scheme that allows to approximate solutions of non-linear systems with a convergence rate of 1 even for weak solutions. This is a particularly interesting property. Unfortunately, the Glimm's scheme is restricted to 1D problems.

Recently, a scheme has been proposed on the basis of Glimm's idea for simulating front propagation [10, 7]. Even if it was restricted to the specific class of the problem of advection of a Heaviside function, it benefited from the advantage to be able to handle 2D or 3D problems. The projection step nicknamed GRU [10] (for Glimm Random Update) is based on the idea of Glimm [8], but also on the analysis of the Upwind scheme proposed in [4]. The GRU step is used as a projection step which is performed at each time-iteration, just after the prediction of the advection through classical finite-volume methods (for instance the prediction step can be done using the Upwind scheme). This projection step is based on a random choice. One of the main advantage of this technique is that the fronts defined by the Heaviside functions remain perfectly sharp in the sense that no intermediate value is created.

A theoretical study of convergence of this algorithm has been proposed in [7] for a planar front. The convergence in probability with order 1 has been proved, as for the Glimm's scheme. This result can also be connected to the results of [4]. In work [4], the authors present a probabilistic analysis of the classical finite-volume Upwind scheme. The approximated solutions of the latter are built as the expectation of a stochastic process, where the characteristics are stochastic and follow a Markov chain. It arises from this analysis that the effective order of 1/2 of the Upwind scheme can be associated with the fluctuations of the stochastic characteristics around the average characteristic, and to the expectation operator applied to recover the Upwind scheme. The basic idea of the scheme proposed in [10, 7] is thus that, following Glimm's idea, selecting only one sample of this stochastic process avoid to apply an expectation operator. Doing so, the effective convergence rate is higher than 1/2 and accuracy of the approximated solutions are significantly improved.

In the present paper, the GRU step is adapted to the more general problem of the advection of a scalar quantity. We intend here to perform some tests in the framework of the compressible flows. For that purpose, we consider in this paper the following one-dimensional system in conservative form:

$$ \left\{ \begin{array}{l} \partial_t (\rho) + \partial_x (\rho U) = 0, \\ \partial_t (\rho\phi) + \partial_x (\rho U\phi) = 0, \end{array} \right. \qquad (1) $$

where the initial condition for the scalar quantity $\phi(t, x)$ is: $\forall x$, $\phi(t = 0, x) = \phi_0(x)$. System (1) involves a non-negative density $\rho(t, x) > 0$ and a velocity $U(t, x)$ which are assumed to be bounded. Obviously, since system (1) is only based on two equations and three unknowns, $\rho$ or $U$ have to be specified. In the sequel, we focus on two classes of solutions. The first one relies on the linear advection problem for which the velocity field and the density are constant and uniform, $U(t, x) = U_0 > 0$ and $\rho(t, x) = \rho_0$ for all $(t, x)$. Secondly, the specific class of solutions described in appendix 6.2 are considered. The latter possesses non-uniform but smooth profiles for the density and for the velocity.

The paper is organized as follows. After a presentation of the scheme in section 2.1, some of its properties are studied in section 2.2. One of the most interesting point is that no intermediate values are created by the scheme. Indeed, at any iteration the cell values are chosen in the existing cell values, so that they already existed in the approximated initial solution. For simple 1D configuration (i.e. with uniform meshes, and for uniform density and velocity), the scheme corresponds to the classical Glimm's scheme [8]. But on the contrary to the latter, the scheme proposed here extends very easily to the multi-dimensional setting. For assessing the behavior of this scheme, several tests are performed in section 3 using the GRU step associated with the Upwind scheme. Comparisons between this scheme and the sole Upwind scheme (without the GRU step) are proposed on the basis of analytical solutions and convergence curves. An extension of the scheme to the 2D (or 3D) advection problem is also tested on unstructured meshes with very satisfactory results, as in [10] for the transport of Heaviside functions. At last, section 4 provides some elements of comparison of the numerical results when considering different pseudo- or quasi-random sequences. This last section is a try to assess the behavior of the GRU scheme with regards to the theoretical results obtained in [11] for the Glimm's scheme. The numerical tests of section 4 seem to show that the use of pseudo-random sequences with very low discrepancy is crucial, as pointed out in a theoretical manner in [11].
## 2 Extension of the GRU scheme to the one-dimensional advection problem

The scheme proposed here is a direct extension of the method introduced in [10]. The latter is based on a two-step algorithm involving: a first step accounting for the convection terms, and for which a wide choice of numerical scheme are possible; followed by a projection step based on a random choice. The new feature in this method obviously entirely relies on the projection step which has been denoted in [10] by the acronym GRU, for Glimm Random Update. In [10], the method has been applied to the front advection and it has been shown that for simple 1D cases, it is equivalent to the well-known Glimm's scheme originally proposed in [8] and studied for instance in [1, 2, 3]. In this section, it is shown that it can be easily extended to the advection of any scalar function.

It should be emphasized that, even if Upwind-like schemes are considered in the rest of the paper, many other schemes may be used for the prediction step, see [10].

### 2.1 Description of the numerical scheme

As mentioned above, the method studied here is based on a prediction step, followed by a projection step. For the prediction step, we restrict here to an Upwind scheme described in section 2.1.1. Moreover, we consider system (1) where both the velocity and the density are supposed to be analytically known.

In the following, we use a finite-volume approach: at any time $t^n$ and for all cells $i$, the value of $\phi$ is approximated by $\phi_i^n$. The value of the scalar quantity $\phi_i^n$ is first updated into $\phi_i^{n+1,*}$ by the prediction step described in section 2.1.1. Afterwards, $\phi_i^{n+1}$ is obtained from the predicted value $\phi_i^{n+1,*}$ through the GRU projection step described in section 2.1.2. For the 1D test cases we consider the domain $[0, 1]$ split into cells of length $\Delta x_i$. Cell $i$ thus corresponds to the interval $[x_{i-1/2}, x_{i+1/2}]$, with $x_{i+1/2} = x_{i-1/2} + \Delta x_i$. The center of cell $i$ is $x_i = x_{i-1/2} + \Delta x_i/2$. The time-step at iteration $n$ is denoted by $\Delta t^n$.

#### 2.1.1 The prediction step

In this work, we use a specific scheme based on the classical Upwind scheme for which the analytical solution for $\rho$ and $U$ is supposed to be known. The properties of this scheme have been investigated in [10]. The scheme reads:

$$ \left\{ \begin{aligned} \rho_i^{n+1,*} \phi_i^{n+1,*} &= \rho_i^n \phi_i^n - \frac{\Delta t^n}{\Delta x_i} (Q_{i+1/2}^n \phi_{i+1/2}^n - Q_{i-1/2}^n \phi_{i-1/2}^n), \\ \rho_i^{n+1,*} &= \rho_i^n - \frac{\Delta t^n}{\Delta x_i} (Q_{i+1/2}^n - Q_{i-1/2}^n), \end{aligned} \right. \qquad (2) $$

with the mass fluxes: $Q_{i+1/2}^n = \rho_{i+1/2}^n U_{i+1/2}^n$. In our particular case, the density and the velocity are given by explicit formula, for instance by a solution as those presented in appendix 6.2. We thus choose here: $\rho_{i+1/2}^n = \rho(t^n, x_{i+1/2})$, $U_{i+1/2}^n = U(t^n, x_{i+1/2})$ and $\rho_i^n = \rho(t^n, x_i)$. Moreover, it should be noted that in general $\rho_i^{n+1,*}$ may be different from $\rho_i^{n+1} = \rho(t^n + \Delta t^n, x_i)^1$. The interfacial value $\phi_{i+1/2}^n$ is then:

$$ \phi_{i+1/2}^n = \begin{cases} \phi_i^n, & \text{if } U_{i+1/2}^n > 0 \\ \phi_{i+1}^n, & \text{otherwise.} \end{cases} \qquad (3) $$

This scheme is based on the classical Upwind scheme but it embeds the analytical solutions for $\rho$ and $U$. It is named Upwind in the following but it should not be mistaken with the classical Upwind scheme, except for the particular cases of constant and uniform density and velocity where both schemes coincide.

Some properties of this scheme are given in [10] when considering uniform meshes, i.e. for $\Delta x_i = \Delta x$, and constant time steps $\Delta t^n = \Delta t$. In particular it is proven that it is monotonicity preserving for $\phi$ providing that the time-step $\Delta t > 0$ fulfills the CFL constraint:

$$ \Delta t < \frac{\Delta x}{2 \max_i \left( |U_{i-1/2}^n| \right)}. \qquad (4) $$

It should be noted that, in order to get this property, it is mandatory for the scheme (2) to be such that $\rho_i^{n+1,*} > 0$. For the particular class of solutions that are proposed in appendix 6.2 and that are used for the numerical tests of section 3, it was shown in [10] that a constraint on the time-step may arise. When considering the particular solutions used here, this additional constraint is less restrictive that the CFL constraint (4). Hence, for all the test cases proposed in this work, we have $\rho_i^{n+1,*} > 0$ and prediction step (3) with uniform meshes and a constant time step is thus monotonicity preserving. This is an essential property for the GRU projection step.

<footnote n="1">¹Indeed, thanks to a Taylor expansion, one can easily obtain that: $\rho_i^{n+1,*} - \rho_i^{n+1} = o(\Delta t)$. See the dedicated appendix in [10] for details.</footnote>
## 2.1.2 The projection step or Glimm Random Update

The GRU step consists in projecting the values $\phi_i^{n+1,*}$ obtained after the prediction step. For that purpose, we consider a quasi-random number $\omega^n$ that follows a uniform distribution in $[0, 1]$ and denoted by $\mathcal{U}(0, 1)$. The projection step then simply reads:

$$ \phi_i^{n+1} = \begin{cases} \phi_{i,m}^n, & \text{if } \phi_i^{n+1,*} < \omega_i^n, \\ \phi_{i,M}^n, & \text{otherwise;} \end{cases} \qquad (5) $$

where $\phi_{i,m}^n$ and $\phi_{i,M}^n$ are respectively the local minimum and the local maximum of the cell values at time $t^n$ when considering the set of the upwind cells with respect to the mass fluxes $Q_{i-1/2}^n$ and $Q_{i+1/2}^n$:

$$ \phi_{i,m}^n = \min_{j \in I_{i,up}^n} (\phi_j^n), \qquad (6) $$

$$ \phi_{i,M}^n = \max_{j \in I_{i,up}^n} (\phi_j^n), \qquad (7) $$

$$ I_{i,up}^n = \{i\} \cup \{(i-1)\delta_{(Q_{i-1/2}^n > 0)}\} \cup \{(i+1)\delta_{(Q_{i+1/2}^n < 0)}\}. \qquad (8) $$

It should be noted that the set $I_{i,up}^n$ contains at least the local cell $i$, it is thus never empty. When $Q_{i+1/2}^n \le 0$, the index $i+1$ is added to the set $I_{i,up}^n$, and when $Q_{i-1/2}^n > 0$ the index $i-1$ is added to the set $I_{i,up}^n$. The number $\omega_i^n$ is a renormalization of the number $\omega^n$ over $[\phi_{i,m}^n, \phi_{i,M}^n]$:

$$ \omega_i^n = \phi_{i,m}^n + \omega^n(\phi_{i,M}^n - \phi_{i,m}^n), \qquad (9) $$

so that $\omega_i^n$ follows the uniform distribution $\mathcal{U}(\phi_{i,m}^n, \phi_{i,M}^n)$. An important point to be quoted here is that the same $\omega^n$ is used for all the cells. This is a cornerstone of the algorithm, as noticed in [1, 2] for the Glimm's scheme. Moreover, the proof of convergence proposed in [7] for the advection of indicator functions also strongly relies on this feature.

In a practical point of view, $\omega^n$ is chosen in low-discrepancy sequences with values in $[0, 1]$. They are well-suited for the Glimm's scheme [1, 11, 3] as well as for the GRU step [10]. In the present work, Halton-Van der Corput sequences are used, except in section 4 where some tests have been performed with other sequences.

## 2.2 Properties of the numerical scheme

In this section we investigate some properties of the whole scheme described in the previous section. As the projection step involves a random choice, these properties are often expressed using probability of occurrence. Nevertheless, the simplicity of the scheme allows to use very simple probabilistic tools, which facilitate the understanding of the behavior of the scheme.

For the sake of simplicity, we assume that the velocity field and the density are such that $\partial_x (\rho U) = 0$, and that $U(t, x) > 0$. The first equation of system (1) then implies that the density is constant: $\partial_t (\rho) = 0$. The mass flow rate denoted by $Q$ is thus uniform with $Q(t) = \rho(x)U(t, x) > 0$. Prediction step (2) can then be simplified since the discrete mass fluxes are $Q_{i+1/2}^n = Q(t^n)$ for all interface $i + 1/2$. Moreover, the approximated density is constant in each cell $\rho_i^{n+1,*} = \rho_i^n = \rho_i^0$, and the update of the approximated value for the scalar $\phi$ reads:

$$ \phi_i^{n+1,*} = \phi_i^n - \beta_i^n (\phi_i^n - \phi_{i-1}^n) = (1 - \beta_i^n) \phi_i^n + \beta_i^n \phi_{i-1}^n, \qquad (10) $$

where we have set:

$$ \beta_i^n = \frac{Q(t^n) \Delta t^n}{\rho_i^0 \Delta x_i}. $$

It is an important to quote that $\beta_i^n$ does not depend on $\phi_i^n$ here, it is thus not stochastic.

Provided that $\beta_i^n < 1$ for all cells $i$, it can easily be shown that scheme (10) is monotonicity preserving and TVD [9]. The time step $\Delta t^n$ at iteration $n$ is chosen according to the constraint: $\max_i(\beta_i^n) < 1$, so that we have $\beta_i^n \in [0, 1]$. Therefore, it straightforwardly arises from prediction step (10) that $\phi_i^{n+1,*} \in [\phi_{i,m}^n, \phi_{i,M}^n]$.

### 2.2.1 Transitional probabilities

Before pursuing this section, let us exhibit two probabilities that will be used throughout section 2.2. Considering prediction step (10), update (5)-(9) gives the following transitional probabilities on $\phi_i^{n+1}$ when knowing the sequence $(\phi_j^n)_j$ at iteration $n$:

$$ \mathcal{P}(\phi_i^{n+1} = \phi_{i,m}^n | (\phi_i^n)_i) = \alpha_i^{n,*}, \quad \text{and} \quad \mathcal{P}(\phi_i^{n+1} = \phi_{i,M}^n | (\phi_i^n)_i) = 1 - \alpha_i^{n,*}, \qquad (11) $$
with

$$ \alpha_{i}^{n,*} = \begin{cases} 1 & \text{if } \phi_{i,m}^{n} = \phi_{i,M}^{n}, \\ \left( \frac{\phi_{i,M}^{n} - \phi_{i}^{n+1,*}}{\phi_{i,M}^{n} - \phi_{i,m}^{n}} \right) & \text{otherwise.} \end{cases} \qquad (12) $$

When $\phi_i^n \neq \phi_{i-1}^n$, we thus have $\alpha_i^{n,*} = 1 - \beta_i^n$ if $\phi_i^n < \phi_{i-1}^n$ and $\alpha_i^{n,*} = \beta_i^n$ if $\phi_i^n > \phi_{i-1}^n$. It should be noted that by definition we have $\alpha_i^{n,*} \in [0, 1]$, this is a consequence of the fact that $\phi_i^{n+1,*}$ and $\omega_i^n$ both belong to the same interval $[\phi_{i,m}^n, \phi_{i,M}^n]$. The property of preservation of the monotonicity by the prediction step is thus essential here. Moreover, the case $\phi_{i,m}^n = \phi_{i,M}^n$ is singular. Indeed, it implies that $\phi_i^n = \phi_{i-1}^{n-1}$ and that $\phi_i^{n+1,*}$ can be equivalently projected onto $\phi_{i,m}^n$ or $\phi_{i,M}^n$, i.e. onto $\phi_i^n$ or $\phi_{i-1}^n$. For this peculiar situation, we chose here $\alpha_i^{n,*} = 1$, but $\alpha_i^{n,*} = 0$ is a strictly equivalent choice. Furthermore, the projection step written using formulae (5)-(9) naturally degenerates in the case where $\phi_{i,m}^n = \phi_{i,M}^n$.

### 2.2.2 Statistical consistency of the scheme

On the basis of the transitional probabilities (11), the transitional expectation of $\phi_i^{n+1}$ knowing $(\phi_j^n)_j$ can be written:

$$ \begin{aligned} E(\phi_i^{n+1} | (\phi_i^n)_i) &= \mathcal{P}(\phi_i^{n+1} = \phi_{i,m}^n | (\phi_i^n)_i) \quad \phi_{i,m}^n + \mathcal{P}(\phi_i^{n+1} = \phi_{i,M}^n | (\phi_i^n)_i) \quad \phi_{i,M}^n, \\ E(\phi_i^{n+1} | (\phi_i^n)_i) &= \alpha_i^{n,*} \phi_{i,m}^n + (1 - \alpha_i^{n,*}) \phi_{i,M}^n, \\ E(\phi_i^{n+1} | (\phi_i^n)_i) &= \phi_i^{n+1,*}. \end{aligned} $$

Thanks to prediction step (10), one can easily obtain that:

$$ E(\phi_i^{n+1} | (\phi_i^n)_i) = \phi_i^n - \beta_i^n (\phi_i^n - \phi_{i-1}^n). \qquad (13) $$

By integrating relation (13) over all the realizations of $(\phi_j^n)_j$, and since $\beta_i^n$ is not stochastic, we obtain:

$$ E(\phi_i^{n+1}) = E(\phi_i^n) - \beta_i^n (E(\phi_i^n) - E(\phi_{i-1}^n)). \qquad (14) $$

Relation (14) exactly corresponds to prediction step (10) applied to the expectation of $(\phi_j^n)_j$. It expresses the statistical consistency of the GRU step with the prediction step; or in other word that no statistical bias with respect to the prediction step is introduced during the projection step. Moreover, it can be deduced from (14) that:

(i) statistically the whole scheme preserves the monotonicity,

(ii) statistically a maximum principle holds,

(iii) statistically the scheme is conservative.

These three properties are directly inherited from the prediction step which is performed in the sequel through the Upwind scheme. Counterparts of properties (i) and (iii) do not hold for a single realization, i.e. for a single sequence $(\omega^n)_n$, while it will be shown that a maximum principle holds for each realization (i.e. the counterpart of point (ii)).

**Remark.** An additional result for the statistical consistency of the algorithm may be found in appendix 6.1. In particular, it appears that the use of the uniform distribution in the GRU step is a key point.

**Remark.** In a more general setting where $\rho$ and/or $U$ depend on $\phi$, some bias could arise in the statistical consistency relation. Indeed, expectations and integrations above should be done on products that contain dependent terms. This difficulty could arise from the coefficient $\beta_i^n$ that could depend on $\phi$ and that could become stochastic. Or, for instance, when considering the classical partial-mass $\rho(\Phi) \times \Phi$ in two-phase flow models, we would have:

$$ E(\rho(\phi) \phi) = E(\rho(\phi)) \times E(\phi) + E(\rho(\phi)' \phi'), $$

where $\Psi' = \Psi - E(\Psi)$ denotes the fluctuation of a quantity $\Psi$ with respect to the expectation $E(\cdot)$. The additional correlation term $E(\rho(\phi)'\phi')$ is in general not null and the statistical consistency could be expressed in a more complex form than (14). This point may have important consequences when attempting to extend the present approach for instance to multiphase flows.

### 2.2.3 Maximum principle and preservation of the monotonicity

It has been shown in the previous section that statistical consistency relation (14) and the form of prediction step (10) allow to deduce that the scheme statistically preserves the monotonicity and that a maximum principle statistically holds (i.e. on $E(\phi_i^n)$). We examine here these properties for a single realization, i.e. for one particular sequence $(\omega^n)_n$.

**Maximum principle.**

First, since the prediction step preserves the monotonicity, the definition of the projection step obviously
implies that all the values $(\phi_i^{n+1})_i$ at iteration $n+1$ are equal to one of the value of the sequence $(\phi_i^n)_i$ at iteration $n$. Hence the maximum principle naturally holds for the whole scheme. In particular, this also implies that at any iteration $n \ge 0$ and for any cell $i$, the value $\phi_i^n$ necessarily belongs to the set $\{\phi_i^0, i = 1..N_c\}$, where $N_c$ is the number of cells in the mesh. Moreover, between two successive iterations the number of the different values reached by the sequence $(\phi_i^n)_i$ decreases or remains constant, but it can not increase because no new cell-value is created.

### Preservation of the monotonicity.

Let us assume that at iteration $n$ the sequence $(\phi_i^n)_i$ is decreasing, so that for all $i$ we have $\phi_i^n > \phi_{i+1}^n$. We have chosen here a uniform positive velocity field, we then easily get that the local minimum and maximum used for the projection step are: $\phi_{i,m}^n = \phi_i^n$ and $\phi_{i,M}^n = \phi_{i-1}^n$, and thus we have

$$ \phi_{i,m}^n = \phi_{i+1,M}^n. \qquad (15) $$

Moreover, whatever $\omega^n$ is, we have $\phi_i^{n+1} \ge \phi_{i,m}^n$ and $\phi_{i+1,M}^n \ge \phi_{i+1}^{n+1}$, so that:

$$ \phi_i^{n+1} \ge \phi_{i,m}^n = \phi_{i+1,M}^n \ge \phi_{i+1}^{n+1}, $$

which enables to conclude that $(\phi_i^{n+1})_i$ is non-increasing.

It is important to emphasize that a decreasing sequence may become non-increasing with a non-null probability. Let us estimate this probability by considering a decreasing sequence $(\phi_i^n)_i$. Such a situation occurs at iteration $n+1$ between cells $i$ and $i+1$ if and only if: $\phi_{i+1}^{n+1} = \phi_i^{n+1}$. Because of (15), this equality can be obtained if and only if $\phi_{i+1}^{n+1} = \phi_{i+1,M}^n$ and $\phi_i^{n+1} = \phi_{i,m}^n$, or in other words for:

$$ \left\{ \begin{array}{l} \omega_i^n \ge \phi_i^{n+1,*}, \\ \text{and} \\ \omega_{i+1}^n < \phi_{i+1}^{n+1,*}. \end{array} \right. \qquad (16) $$

Since it has been assumed $U > 0$, we get from the projection step that $\phi_{i,m}^n = \phi_i^n$ and $\phi_{i,M}^n = \phi_{i-1}^n$, and thus that:

$$ \omega_i^n = \phi_i^n + \omega^n(\phi_{i-1}^n - \phi_i^n). $$

By combining this result and prediction step (10) in relations (16), one obtain the equivalent relations:

$$ \left\{ \begin{array}{l} \omega^n \ge \beta_i^n, \\ \text{and} \\ \omega^n < \beta_{i+1}^n. \end{array} \right. \qquad (17) $$

If $\beta_i^n \ge \beta_{i+1}^n$, then (17) is not possible. On the contrary, when $\beta_i^n < \beta_{i+1}^n$, (17) can be fulfilled with a probability $\beta_{i+1}^n - \beta_i^n$ ($\omega^n$ is chosen according to a uniform distribution $U(0,1)$). Therefore, it can be concluded that for a decreasing sequence at iteration $n$, we can get $\phi_{i+1}^{n+1} = \phi_i^{n+1}$ with the probability:

$$ \mathcal{P} ( \phi_{i+1}^{n+1} = \phi_i^{n+1} | (\phi_j^n)_j \text{ decreasing} ) = \max (0, \beta_{i+1}^n - \beta_i^n). $$

Moreover, if we indeed get $\phi_{i+1}^{n+1} = \phi_i^{n+1}$ at iteration $n+1$ from a decreasing sequence $(\phi_j^n)_j$ at iteration $n$, then we necessary have $\phi_{i-1}^{n+1} \ne \phi_i^{n+1}$ and $\phi_{i+1}^{n+1} \ne \phi_{i+2}^{n+1}$. This is quite straightforward since in this case we have:

$$ \phi_i^{n+1} = \phi_{i,m}^n = \phi_i^n < \phi_{i-1,m}^n = \phi_{i-1}^n $$

and

$$ \phi_{i+1}^{n+1} = \phi_{i+1,M}^n = \phi_i^n > \phi_{i+2,M}^n = \phi_{i+1}^n. $$

In other words, starting with a decreasing sequence $(\phi_j^n)_j$ at iteration $n$, the sequence at iteration $n+1$ contains sets of at most 2 neighboring cells sharing the same value.

Obviously, the same properties hold when considering increasing sequences. These results can also be locally extended to sequences that are only locally monotonic.

### 2.2.4 Behavior at the local extrema and on locally uniform sequences

In the previous section, it has been shown that the scheme fulfills the maximum principle and that a locally strictly monotonic sequence can become locally uniform from an iteration to the next one. Let us now study two other properties. Firstly, even if the maximum principle holds, it will be seen that local extrema may be lost. Secondly, the variation of the width of a plateau in a locally uniform sequence is studied.

#### Loss of a local extremum.

Let us consider a sequence $(\phi_j^n)_j$ with a local maximum at cell $i+1$. This situation can be depicted by figure (1b) with $N=1$ and with $\phi_{i-1} = \phi_a$ and $\phi_{i+3} = \phi_c$. The local extremum of cell $i+1$ is lost between iteration $n$ and iteration $n+1$ if and only if $\phi_{i+1}^{n+1} = \phi_{i+1,m}^n$ and $\phi_{i+2}^{n+1} = \phi_{i+2,m}^n$. Considering projection
step (5) and the choice $U > 0$, we have: $\phi_{i+1,m}^n = \phi_i^n$ and $\phi_{i+1,M}^n = \phi_{i+1}^n$, together with $\phi_{i+2,m}^n = \phi_{i+2}^n$ and $\phi_{i+2,M}^n = \phi_{i+1}^n$. Hence, we obtain the following thresholds in cells $i+1$ and $i+2$ for the projection step:

$$ \left\{ \begin{array}{l} \omega_{i+1}^n = \phi_i^n + \omega^n(\phi_{i+1}^n - \phi_i^n), \\ \omega_{i+2}^n = \phi_{i+2}^n + \omega^n(\phi_{i+1}^n - \phi_{i+2}^n). \end{array} \right. $$

After some straightforward calculus, the projection step associated with prediction step (10) then gives:

$$ \left\{ \begin{array}{l} \phi_{i+1}^{n+1} = \phi_{i+1,m}^n, \\ \phi_{i+2}^{n+1} = \phi_{i+2,m}^n, \end{array} \right. \iff \left\{ \begin{array}{l} \omega_{i+1}^n \ge \phi_{i+1}^{n+1,*}, \\ \omega_{i+2}^n \ge \phi_{i+2}^{n+1,*}, \end{array} \right. \iff \left\{ \begin{array}{l} (\phi_i^n - \phi_{i+1}^n) \ge (\phi_i^n - \phi_{i+1}^n)(\beta_{i+1}^n + \omega^n), \\ \omega^n(\phi_{i+1}^n - \phi_{i+2}^n) \ge (\phi_{i+1}^n - \phi_{i+2}^n)\beta_{i+2}^n. \end{array} \right. \quad (18) $$

In the present case, we have assumed a local maxima in cell $i+1$ at iteration $n$, we then have: $\phi_i^n - \phi_{i+1}^n < 0$ and $\phi_{i+1}^n - \phi_{i+2}^n > 0$. As a consequence, the last set of relations of (18) is equivalent to:

$$ \left\{ \begin{array}{l} 1 \le (\beta_{i+1}^n + \omega^n), \\ \text{and} \\ \omega^n \ge \beta_{i+2}^n, \end{array} \right. \iff \left\{ \begin{array}{l} \omega^n \ge 1 - \beta_{i+1}^n, \\ \text{and} \\ \omega^n \ge \beta_{i+2}^n. \end{array} \right. \quad (19) $$

Due to the uniform distribution for $\omega^n$, this event occurs with a probability:

$$ \mathcal{P}(\phi_{i+1}^{n+1} = \phi_{i+1,m}^n \text{ and } \phi_{i+2}^{n+1} = \phi_{i+2,m}^n \mid \phi_{i-1}^n = \phi_i^n < \phi_{i+1}^n \text{ and } \phi_{i+1}^n > \phi_{i+2}^n = \phi_{i+3}^n) = \min(1 - \beta_{i+2}^n, \beta_{i+1}^n). \quad (20) $$

A similar result can be obtained for a local minimum. Moreover, in the particular case with uniform coefficients $\beta_i^n = \beta_0^n$, probability (20) becomes $\min(1-\beta_0^n, \beta_0^n)$. It thus reaches a maximum of $1/2$ for $\beta_0^n = 1/2$.

### Variation of the width of a plateau.

When the approximated solution contains several successive uniform plateaux, their width in terms of number of cells can vary with the time-iterations. Such patterns may obviously be part of the approximated solution, or they may have been created by the projection step from a strictly monotonic sequence as depicted in section 2.2.3. Let us estimate some probabilities on the width this uniform plateaux. Two typical situations considered here are depicted in figure 1 (with $\phi_b > \phi_c$, the cases $\phi_b < \phi_c$ have not been represented). It is assumed that $\phi_a \ne \phi_b$ and $\phi_b \ne \phi_c$. In the following, $N_b^n$ denotes the width of the plateau with value $\phi_b$ in figure 1. It corresponds to the number of successive cells that share the same value $\phi_b$ at iteration $n$. It is assumed that $N_b^n = N \ge 2$. After prediction step (10), one obtain:

$$ \left\{ \begin{array}{l} \phi_{i+1}^{n+1,*} = \beta_{i+1}^n \phi_a + (1 - \beta_{i+1}^n) \phi_b, \\ \phi_{i+N+1}^{n+1,*} = \beta_{i+N+1}^n \phi_b + (1 - \beta_{i+N+1}^n) \phi_c, \\ \phi_k^{n+1,*} = \phi_k^*, \quad \text{if } k \notin \{i+1, i+N+1\}, \end{array} \right. \quad (21) $$

and the projection thresholds read:

$$ \left\{ \begin{array}{l} \omega_{i+1}^n = m_{a,b} + \omega^n(M_{a,b} - m_{a,b}), \\ \omega_{i+N+1}^n = m_{b,c} + \omega^n(M_{b,c} - m_{b,c}), \end{array} \right. \quad (22) $$

where for the sake of readability it has been set: $m_{a,b} = \min(\phi_a, \phi_b)$, $m_{b,c} = \min(\phi_b, \phi_c)$, $M_{a,b} = \max(\phi_a, \phi_b)$ and $M_{b,c} = \max(\phi_b, \phi_c)$. From the definition of the scheme and thanks to (21) and (22), some straightforward calculus lead to:

$$ \phi_{i+1}^{n+1,*} < \omega_{i+1}^n \iff B_{i+1}^n = \frac{\phi_b - m_{a,b}}{M_{a,b} - m_{a,b}} + \beta_{i+1}^n \frac{\phi_a - \phi_b}{M_{a,b} - m_{a,b}} < \omega^n, \quad (23) $$

and

$$ \phi_{i+N+1}^{n+1,*} < \omega_{i+N+1}^n \iff B_{i+N+1}^n = \frac{\phi_c - m_{b,c}}{M_{b,c} - m_{b,c}} + \beta_{i+N+1}^n \frac{\phi_b - \phi_c}{M_{b,c} - m_{b,c}} < \omega^n. \quad (24) $$

It should be remarked that when $\phi_a > \phi_b$, the terms in the second equality of (23) are:

$$ (\phi_b - m_{a,b})/(M_{a,b} - m_{a,b}) = 0 \quad \text{and} \quad (\phi_a - \phi_b)/(M_{a,b} - m_{a,b}) = 1, $$

and thus we simply get $B_{i+1}^n = \beta_{i+1}^n$. Conversely, when $\phi_a < \phi_b$, it yields:

$$ (\phi_b - m_{a,b})/(M_{a,b} - m_{a,b}) = 1 \quad \text{and} \quad (\phi_a - \phi_b)/(M_{a,b} - m_{a,b}) = -1, $$

and hence $B_{i+1}^n = 1 - \beta_{i+1}^n$. Obviously, the same remark holds for the ratios involving $\phi_b$ and $\phi_c$ in (24) and for the definition of $B_{i+N+1}^n$ which can be simplified in: $B_{i+N+1}^n = \beta_{i+N+1}^n$ if $\phi_b > \phi_c$, and $B_{i+N+1}^n = 1 - \beta_{i+N+1}^n$ otherwise. These notations will be helpful hereafter for expressing general expressions for the probabilities.

From now, let us focus on the specific case where $\phi_a > \phi_b > \phi_c$, i.e. on the case depicted by figure 1b. After the projection step, the initial front between cells $i$ and $i+1$ (resp. between cells $i+N$ and $i+N+1$): remains between cells $i$ and $i+1$ (resp. between cells $i+N$ and $i+N+1$) when $\omega_{i+1}^n > \phi_{i+1}^{n+1,*}$ (resp. $\omega_{i+N+1}^n > \phi_{i+N+1}^{n+1,*}$) or is translated at the interface between cells $i+1$ and $i+2$ (resp. between cells $i+N+1$ and $i+N+2$) when $\omega_{i+1}^n < \phi_{i+1}^{n+1,*}$ (resp. $\omega_{i+N+1}^n < \phi_{i+N+1}^{n+1,*}$). There are thus four possible configurations
at iteration $n+1$. Among these four situations, two are such that $\mathcal{N}_b^{n+1} = \mathcal{N}_b^n$ and for the two others $\mathcal{N}_b^{n+1}$ changes while the event $|\mathcal{N}_b^{n+1} - \mathcal{N}_b^n| \ge 2$ is not possible. Therefore, the three following possibilities are obtained for $\mathcal{N}_b^{n+1}$<fnref n="2" />:

$$ \mathcal{N}_b^{n+1} = \mathcal{N}_b^n \iff \left\{ \begin{array}{ll} \omega_{i+1}^n < \phi_{i+1}^{n+1,*} & \text{and } \omega_{i+N+1}^n < \phi_{i+N+1}^{n+1,*}, \\ & \text{or} \\ \omega_{i+1}^n > \phi_{i+1}^{n+1,*} & \text{and } \omega_{i+N+1}^n > \phi_{i+N+1}^{n+1,*}, \end{array} \right. \qquad (25) $$

$$ \mathcal{N}_b^{n+1} = \mathcal{N}_b^n - 1 \iff \left( \omega_{i+1}^n < \phi_{i+1}^{n+1,*} \text{ and } \omega_{i+N+1}^n > \phi_{i+N+1}^{n+1,*} \right), \qquad (26) $$

$$ \mathcal{N}_b^{n+1} = \mathcal{N}_b^n + 1 \iff \left( \omega_{i+1}^n > \phi_{i+1}^{n+1,*} \text{ and } \omega_{i+N+1}^n < \phi_{i+N+1}^{n+1,*} \right). \qquad (27) $$

Thanks to relations (23) and (24) the probabilities for the three events above can be explicitly written using $B_{i+1}^n$ and $B_{i+N+1}^n$. The probability to maintain the number of cells in the plateau with value $\phi_b$ is then:

$$ \mathcal{P}(\mathcal{N}_b^{n+1} = \mathcal{N}_b^n | \mathcal{N}_b^n) = 1 + \min(B_{i+1}^n, B_{i+N+1}^n) - \max(B_{i+1}^n, B_{i+N+1}^n), \qquad (28) $$

and the probabilities to decrease or to increase the number of cells with value $\phi_b$ are respectively:

$$ \mathcal{P}(\mathcal{N}_b^{n+1} = \mathcal{N}_b^n - 1 | \mathcal{N}_b^n) = \max(0, B_{i+1}^n - B_{i+N+1}^n), \qquad (29) $$

and

$$ \mathcal{P}(\mathcal{N}_b^{n+1} = \mathcal{N}_b^n + 1 | \mathcal{N}_b^n) = \max(0, B_{i+N+1}^n - B_{i+1}^n). \qquad (30) $$

It is an important point to be quoted that when $B_{i+N+1}^n = B_{i+1}^n$, the number of cells in the plateau with value $\phi_b$ is maintained from iteration $n$ to $n+1$ with probability 1 (that point is developed hereafter in this section).

The expectation of the number of cells in the plateau $\phi_b$ knowing $\mathcal{N}_b^n$:

$$ \begin{aligned} E(\mathcal{N}_b^{n+1} | \mathcal{N}_b^n) = & (\mathcal{N}_b^n - 1) \mathcal{P}(\mathcal{N}_b^{n+1} = \mathcal{N}_b^n - 1 | \mathcal{N}_b^n) \\ & + \mathcal{N}_b^n \mathcal{P}(\mathcal{N}_b^{n+1} = \mathcal{N}_b^n | \mathcal{N}_b^n) \\ & + (\mathcal{N}_b^n + 1) \mathcal{P}(\mathcal{N}_b^{n+1} = \mathcal{N}_b^n + 1 | \mathcal{N}_b^n), \end{aligned} \qquad (31) $$

can be expressed by using probabilities (28), (29) and (30). It yields:

$$ E(\mathcal{N}_b^{n+1} | \mathcal{N}_b^n) = \mathcal{N}_b^n - (B_{i+1}^n - B_{i+N+1}^n). \qquad (32) $$

Since the second term on the right hand side of relation (32) is not stochastic, we can easily integrate relation (32) over all the possible $\mathcal{N}_b^n$ and we get:

$$ E(\mathcal{N}_b^{n+1}) = E(\mathcal{N}_b^n) - (B_{i+1}^n - B_{i+N+1}^n). \qquad (33) $$

In fact the second term on the right hand side of (33) (and (32)) corresponds to a compression/dilatation term that represents the variation of the number of cells in the plateau $\phi_b$ when the mesh size is not in agreement with the density profile (which has been assumed constant here). For the present case with $\phi_a > \phi_b > \phi_c$, we have:

$$ E(\mathcal{N}_b^{n+1}) = E(\mathcal{N}_b^n) \iff B_{i+N+1}^n = B_{i+1}^n \iff \beta_{i+N+1}^n = \beta_{i+1}^n \iff \rho_{i+1}^0 \Delta x_{i+1} = \rho_{i+N+1}^0 \Delta x_{i+N+1}. $$

This case corresponds to a specific situation where the mesh size is adapted to the density profile for obtaining the same mass $\rho_j^0 \Delta x_j$ in the two cells $j = i + 1$ and $j = i + N + 1$. It should be noted that in this case the probability (28) is equal to 1. When turning to a non-adapted mesh, such that $\rho_{i+1}^0 \Delta x_{i+1} > \rho_{i+N+1}^0 \Delta x_{i+N+1}$ for instance, we get:

$$ \rho_{i+1}^0 \Delta x_{i+1} > \rho_{i+N+1}^0 \Delta x_{i+N+1} \iff \beta_{i+N+1}^n < \beta_{i+1}^n \iff B_{i+N+1}^n < B_{i+1}^n \iff E(\mathcal{N}_b^{n+1}) > E(\mathcal{N}_b^n). $$

Moreover, from probabilities (28), (29) and (30), we obtain:

$$ \left\{ \begin{array}{ll} \mathcal{P}(\mathcal{N}_b^{n+1} = \mathcal{N}_b^n | \mathcal{N}_b^n) = 1 + B_{i+1}^n - B_{i+N+1}^n < 1, \\ \mathcal{P}(\mathcal{N}_b^{n+1} = \mathcal{N}_b^n - 1 | \mathcal{N}_b^n) = 0, \\ \mathcal{P}(\mathcal{N}_b^{n+1} = \mathcal{N}_b^n + 1 | \mathcal{N}_b^n) = B_{i+N+1}^n - B_{i+1}^n > 0. \end{array} \right. $$

These relations state that if cell $i+1$ contains more mass than cell $i+N+1$, then the scheme has to statistically reduce the number of cells in the plateau $\phi_b$ in order to statistically maintain a constant mass. This behavior is obviously in agreement with the behavior of the solution of system of equations (1) with the additional constraint $\partial_x(\rho U) = 0$ that has been retained all along section 2.2.

<footnote n="2">²For the other cases (where $\phi_a > \phi_b > \phi_c$ does not hold), the inequalities on the right-hand side of the events (25), (26) and (27) have simply to be changed.</footnote>
<figure><img src="image_2.png" /><figcaption>(a) Situation with a plateau of width N with a decreasing sequence.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(b) Situation with a plateau of width N that forms a local extremum.</figcaption></figure>

Figure 1: Variation of the width of a plateau, case with an uniform velocity $U_0 > 0$. The letter M denotes the local maximum $\phi_{i,M}^n$, m denotes the local minimum $\phi_{i,m}^n$, and * the value $\phi_{i}^{n+1,*}$. We have $\phi_a \neq \phi_b$ and $\phi_b \neq \phi_c$.

### 2.2.5 A word on the convergence of the scheme

One of the great advantages of the Glimm's scheme is that it has a convergence rate of 1 even for weak solutions involving discontinuities. Classical first-order schemes have a convergence rate of 1/2 for such solutions. The explanation of this difference can be found in the work [4]. In this paper, the authors present a probabilistic analysis of the classical finite-volume Upwind scheme. The approximated solutions of the latter are built as the expectation of a stochastic process, where the characteristics are stochastic and follow a Markov chain. It arises from this analysis that the effective order of 1/2 of the Upwind scheme can be associated with the fluctuations of the stochastic characteristics around the average characteristic, and to the expectation operator applied to recover the Upwind scheme. The basic idea of the scheme proposed here is thus that, following Glimm's idea, selecting only one sample of this stochastic process avoid to apply an expectation operator. Doing so, we can expect to get an effective convergence rate higher than 1/2 and we could wait for a better accuracy for the approximated solutions. This idea has also been applied in [10] for the GRU projection step on Heaviside functions. Convergence rates close to 0.85 have been retrieved together with a great improvement of the accuracy. In section 3, the behavior of the GRU projection step is assessed considering analytical test cases for scalar advection.

It has been proved in [7] that the GRU projection step associated with the Upwind scheme for the prediction step converges in probability with order 1. This proof has been obtained for solutions based on a Heaviside function forming a planar front and for a uniform velocity field. For the setting studied in this work, the proof proposed in [7] could be extended in order to obtain the convergence of the scheme when considering piecewise constant functions.

## 3 Numerical tests for assessing the behavior of the scheme

In the following, several numerical tests are performed on the basis of analytical solutions. Two velocity fields and two initial conditions are considered. In section 3.1, the case of a constant and uniform velocity field is investigated. It clearly belongs to the class of problems that have been retained in section 2.2 for studying the properties of the scheme. On the contrary, the velocity fields used in section 3.2 are not uniform and they correspond to a compression or a dilatation of the initial condition for $\phi$. These velocity fields arise from the class of analytical solutions described in appendix 6.2. The two initial conditions used for the different test cases are:

$$ \phi_0^s(x) = \begin{cases} 1 + \frac{1}{2} \sin \left( \frac{x - 0.1}{0.3 - 0.1} \right), & \text{if } 0.1 < x < 0.3, \\ 0 & \text{otherwise,} \end{cases} \qquad (34) $$

and

$$ \phi_0^c(x) = \begin{cases} 1.5, & \text{if } 0.1 < x < 0.2, \\ 0.5, & \text{if } 0.2 < x < 0.275, \\ 1, & \text{if } 0.275 < x < 0.3, \\ 0 & \text{otherwise.} \end{cases} \qquad (35) $$

These initial conditions correspond to weak solutions of system (1) since they both involves discontinuities. It should be noted that in section 3.1 the initial condition $\phi_0^c(x - 0.3)$ will be considered for performing periodic simulations. The numerical tests are performed on the interval $[0, 1]$, for uniform meshes and for a CFL number equal to 1/2.

In the projection/GRU step, $\omega^n$ should be chosen according to the uniform distribution. In practice, it has been noticed that equidistributed low-discrepancy sequences give better results [3]. In the following we consider the Halton-Van der Corput sequences as proposed in [3]. Additional tests have been performed in section 4 with other sequences. The Halton-Van der Corput sequences have already been tested for the
GRU step in [10]. The Halton-Van der Corput sequences depend on a couple of parameters ($K_1, K_2$) where $K_1$ and $K_2$ are prime numbers greater or equal to 3, and such that $K_1 > K_2$. Each choice ($K_1, K_2$) gives a different deterministic sequence of numbers in (0, 1) which is denoted by $V_{K_1,K_2}$ in the following. For the projection step, $\omega^n$ is chosen as the $n^{th}$ element of the sequence $V_{K_1,K_2}$: $\omega^n = V_{K_1,K_2}(n)$. Obviously, each couple ($K_1, K_2$) leads to a different approximated solution. More details on the computation of these sequences may be found in appendix 6.3.

In order to produce convergence curves, several couples ($K_1, K_2$) are considered. If we denote by $\mathcal{P}(n)$ the sequence containing the $n$ first prime numbers greater or equal to 3, the couples that are used in the following are all the couples of the set $\mathcal{V}(n) = \{(K_1, K_2) \in \mathcal{P}(n) \times \mathcal{P}(n); 2 < K_2 < K_1\}$. In particular, the results of the present section have been computed for $n = 8$, with $\mathcal{P}(8) = \{3, 5, 7, 11, 13, 17, 19, 23\}$, and therefore 28 different couples ($K_1, K_2$) have been used for each mesh size.

For each mesh size $\Delta x$ and for each couple ($K_1, K_2$), the $L_1$-error between the approximated solution and the analytical solution, $\varepsilon(\Delta x, K_1, K_2)$, is computed. Then, for each mesh size, the $L_1$-errors obtained for all the couple ($K_1, K_2$) are averaged:

$$ \varepsilon(\Delta x) = \frac{1}{card(\mathcal{V}(8))} \sum_{(K_1, K_2) \in \mathcal{V}(8)} \varepsilon(\Delta x, K_1, K_2). \quad (36) $$

We then get an average error estimate $\varepsilon(\Delta x)$ for each mesh size. It is an important point to be quoted here that $\varepsilon(\Delta x)$ is the average of the errors between each approximated solution and the analytical solution and not the error between the average of the approximated solutions and the analytical solution. Since the $L_1$-error operator is not linear, these two errors are different. Moreover, according to section 2.2.2 the average of the approximated solutions statistically converges to the approximated solution of the Upwind scheme. Hence, with the error of the mean approximated solution we would recover the effective order of convergence 1/2 of the classical Upwind scheme, while using the GRU step we could expect an effective order of convergence closer to 1, see the discussion of section 2.2.5.

At last, an extension of the GRU step for the two-dimensional advection problem is tested in section 3.3. These 2D results have been performed considering unstructured meshes and analytical solutions. The GRU step performs very well on these test cases. These satisfactory results could be of great interest for industrial applications.

## 3.1 Approximating solutions for the linear advection problem

We consider here a uniform velocity $U(t, x) = U_0 = 1$ and a uniform density $\rho(t, x) = 1$. The initial condition for $\phi$ is given by $\phi(0, x) = \phi_0(x - 0.3)$, where $\phi_0 = \phi_0^c$ or $\phi_0 = \phi_0^s$ is respectively defined by (35) or by (34). By translating $\phi_0$ from +0.3, the initial profile for $\phi$ is centered on $x = 0.5$. We then consider periodic boundary conditions, so that the analytical solution becomes periodic with a period equal to 1. Several final times are considered: $T_f \in \{1, 2, 5, 10, 20, 50, 100, 1000\}$, respectively corresponding to 1, 2, 5, 10, 20, 50, 100 and 1000 periods of the solution. The approximated solutions obtained for 200 cells are plotted in figure 2 and 3 considering the Halton-Van der Corput sequence $V_{5,3}$. Convergence curves for the two initial conditions and for the periods $\{1, 2, 5, 10, 20, 50, 100, 1000\}$ are reported in figure 4.

On these results, the projection step clearly improves the quality of the approximated solutions. In figure 2, it can be remarked for the test case with $\phi_0^s$ that plateaux are created at each of the three extrema. The values of these extrema have been altered which is possible thanks to the results of section 2.2.4. Moreover, for the initial condition $\phi_0^c$, it clearly appears that the plateaux are well advected by the scheme, which is also in agreement with the results of section 2.2.4. The convergence curves in figure 4a and 4b confirm that the decrease in time of the accuracy of the approximated solution is small, at least with respect to the Upwind scheme. Such a behavior is typical from anti-diffusive schemes as the downwind scheme proposed in [5].

In a quantitative point of view, we recover for the Upwind scheme the classical convergence rate of 1/2, whereas with the GRU step the effective convergence rate is around 0.85. This effective rate of convergence was already observed for Heaviside functions in [10]. In terms of accuracy on coarse meshes, the GRU step provides a gain of one order of magnitude for $\phi_0^c$ and half an order of magnitude for $\phi_0^s$. The former contains more discontinuities and is piecewise uniform, which is clearly in favor of the GRU step.
Comparisons of app. solutions for 200 cells

<figure><img src="image_4.png" /><figcaption>Figure 2: Approximated solutions for the initial condition $\phi_0^c$ for different periods $\{1, 2, 5, 10, 20, 50, 100, 1000\}$ and for a projection step using $(K_1, K_2) = (5, 3)$. The mesh contains 200 cells.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Comparisons of app. solutions for 200 cells - zoom</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 3: Approximated solutions for the initial condition $\phi_0^s$ for different periods $\{1, 2, 5, 10, 20, 50, 100, 1000\}$ and for a projection step using $(K_1, K_2) = (5, 3)$. The mesh contains 200 cells.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 4: Convergence curves for the average error $\epsilon(\Delta x)$ on the set of Halton-Van der Corput sequences $V(8)$ and for the periods $\{1, 2, 5, 10, 20, 50, 100\}$. The black plain lines represent the curves for the Upwind scheme, while the red dashed lines denote the results for the Upwind scheme with GRU projection.</figcaption></figure>

## 3.2 Approximating solutions with non-uniform velocity fields

In this section we consider the three analytical solutions for $U$ and $\rho$ proposed in appendix 6.2 (see figure 19) with the parameters and the final time $T_f$:

$$ (i) \quad U_0 = 2, R_0 = 2, B_0 = -1, T_f = 0.4, $$

$$ (ii) \quad U_0 = -2, R_0 = 2 \ln(2), B_0 = -2, T_f = 0.3, $$

$$ (iii) \quad U_0 = 2, R_0 = 2 \ln(2), B_0 = 2, T_f = 0.3. $$

The velocity and the density are thus not uniform, but the velocity remains constant in time. The first and third sets of parameters correspond to decreasing velocity with respect to $x$ and the second set to an increasing velocity. For sets of parameters (ii) and (iii), the velocity vanishes at $x = 1/2$. The initial conditions for $\phi$ are $\phi_0^c$ and $\phi_0^s$ for set of parameters (i), while for sets (ii) and (iii) the initial conditions are $x \mapsto \phi_0^c(x - 0.3)$ and $x \mapsto \phi_0^s(x - 0.3)$. Due to the divergence of the velocity field, these initial profiles are compressed (for (i) and (iii)) or dilated (for (ii)) along the streamline. It should be noted that for sets of parameters (ii) and (iii), the initial conditions overlaps the stationary point $U(x = 1/2, t) = 0$. Since the velocity is constant, the CFL constraint has been fulfilled by using a constant time-step.

Figure 5 shows the approximated solutions obtained with the Upwind scheme with and without the GRU projection step for the three sets of parameters and for the two initial conditions. These results have been obtained with a (5, 3) Halton-Van der Corput sequence. The convergence curves for the average error $\epsilon(\Delta x)$ and for each individual errors $\epsilon(\Delta x, K_1, K_2)$ with $(K_1, K_2) \in V(8)$ have been plotted in figure 6. As for the test of the previous section, the effective convergence rate for the scheme with the GRU step is close to 0.82 whereas for the Upwind scheme without projection step the convergence rate reaches its classical value 1/2 (for weak solutions involving discontinuities). Moreover, accuracy is globally improved by the projection step. Nonetheless, it should be noted that the improvements of the GRU step is not important for coarse meshes on the test case with parameters (iii) (see figure 6f). This can be explained by the fact that the initial profile is compressed and becomes thin. For coarse meshes, the number of points in the profile becomes small and the probability to lose peaks and to flatten the approximated solution becomes important (see section 2.2.4).
<figure><img src="image_8.png" /><figcaption>(a) Test cases for the set of parameter (i) and with $\Phi_0^c$ as initial condition.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) Test cases for the set of parameter (ii) and with $\Phi_0^c$ as initial condition.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(e) Test cases for the set of parameter (iii) and with $\Phi_0^c$ as initial condition.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Test cases for the set of parameter (i) and with $\Phi_0^s$ as initial condition.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(d) Test cases for the set of parameter (ii) and with $\Phi_0^s$ as initial condition.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(f) Test cases for the set of parameter (iii) and with $\Phi_0^s$ as initial condition.</figcaption></figure>

Figure 5: Comparison of the approximated solutions for different meshes and for the three sets of parameters. The GRU step has been performed here using the (5, 3) Halton-Van der Corput sequence. The black dotted lines correspond to the initial conditions, and the black dashed lines correspond to the final solutions.
<figure><img src="image_14.png" /><figcaption>Figure 6: Convergence curves for the average error $\epsilon(\Delta x)$ (red line with squares) and for $\epsilon(\Delta x, K_1, K_2)$ (turquoise plain lines) on the set of Halton-Van der Corput sequences $V(8)$.</figcaption></figure>

## 3.3 Extension for two-dimensional problems

The two-dimensional counterpart of system (1) reads:

$$ \left\{ \begin{array}{l} \partial_t (\rho) + \partial_x (\rho U_x) + \partial_y (\rho U_y) = 0, \\ \partial_t (\rho \phi) + \partial_x (\rho U_x \phi) + \partial_y (\rho U_y \phi) = 0, \end{array} \right. \qquad (37) $$

where $U_x$ and $U_y$ denote respectively the velocity components along the x and y axis. As for the one-dimensional case, the density $\rho$ and the velocity field $(U_x, U_y)$ are given. For this 2D problem, the prediction step can be straightforwardly extended from section 2.1.1. In the following, the projection step of section
2.1.2 is also extended in a naive manner by modifying definitions (6) and (7) of $\phi_{i,m}^n$ and $\phi_{i,M}^n$:

$$ \phi_{i,m}^n = \min_{j \in \mathcal{N}_i} (\phi_j^n), \qquad (38) $$

$$ \phi_{i,M}^n = \max_{j \in \mathcal{N}_i} (\phi_j^n); \qquad (39) $$

where the set $\mathcal{N}_i$ corresponds to the following neighborhood of cell $i$: $\mathcal{N}_i$ contains the cell-index $i$ of the local cell, and the indices of all the cells that share a common face with cell $i$ and such that the normal mass-flux at this face is oriented inward to cell $i$. It should be noted that $\mathcal{N}_i$ thus does not contain the cells that share only one vertex with cell $i$. It gathers the cell $i$ and all the neighboring cells that are upwind of cell $i$ in the sense of the mass flow.

In order to assess the accuracy of the scheme, numerical approximations are computed on unstructured meshes. These meshes are built on the basis of structured meshes with $N \times N$ square cells. Then, some vertices are modified in order to “twist” the cells around the center (0.5, 0.5) of the square domain $[0, 1] \times [0, 1]$. The transformation of the vertices is based on the polar coordinates centered on the center of the domain $[0, 1] \times [0, 1]$. We set:

$$ r(x, y) = \sqrt{(x - 0.5)^2 + (y - 0.5)^2}, $$

with,

$$ (r_x(x, y), r_y(x, y)) = \left( \frac{0.5 - y}{r}, \frac{x - 0.5}{r} \right) \quad \text{if } r \neq 0; \quad \text{and} \quad r_x(0.5, 0.5) = r_y(0.5, 0.5) = 0.5 $$

and the transformation $(X'_p, Y'_p)$ of coordinates $(X_p, Y_p)$ of vertex $p$ is:

$$ \left\{ \begin{array}{l} X'_{p} = X_{p} + r_{x}(X_{p}, Y_{p}) D_{0} e^{-(\delta r)^{2}}, \\ Y'_{p} = Y_{p} + r_{y}(X_{p}, Y_{p}) D_{0} e^{-(\delta r)^{2}}, \end{array} \right. \qquad (40) $$

where $D_0 = 0.08$ and $\delta r = (r(X_p, Y_p) - R_0)/\sigma$, with $R_0 = 0.2$ and $\sigma = 0.06$. Examples of such meshes are shown in figure 7.

A classical class of problems for assessing scalar transport algorithms consists in simulating the solid rotation of shapes. We consider here the solid rotation around the center (0.5, 0.5) of the computational domain $[0, 1] \times [0, 1]$. The analytical solution for the density is $\rho(t, x, y) = 1.0$ and for the velocity $U_x(t, x, y) = \theta_0(0.5 - y)$ and $U_y(t, x, y) = \theta_0(x - 0.5)$, where $\theta_0 = 1$. With these choices, the initial profile for $\phi$ rotates around the point (0.5, 0.5). We consider the final solution at time $t_{end} = 2\Pi$, so that the initial profile for $\phi$ has achieved a complete rotation, and thus we have $\phi(0, x, y) = \phi(t_{end}, x, y)$. Two initial profiles are tested. We have $\phi(0, x, y) = 0$ everywhere, except in the circle with a center located at (0.5, 0.75) and with a radius of 0.15. For defining the initial conditions, the distance to the center of this circle (0.5, 0.75) is introduced:

$$ \tilde{r}(x, y) = \sqrt{(x - 0.5)^2 + (y - 0.75)^2}. $$

The first initial condition is piecewise uniform:

$$ \phi(0, x, y) = \begin{cases} 2 & \text{if } \tilde{r}(x, y) < 0.15 \text{ and } x > 0.5 \text{ and } y > 0.75, \\ 1.5 & \text{if } \tilde{r}(x, y) < 0.15 \text{ and } x \le 0.5 \text{ and } y > 0.75, \\ 1.0 & \text{if } \tilde{r}(x, y) < 0.15 \text{ and } x \le 0.5 \text{ and } y \le 0.75, \\ 0.5 & \text{if } \tilde{r}(x, y) < 0.15 \text{ and } x > 0.5 \text{ and } y \le 0.75, \\ 0 & \text{otherwise;} \end{cases} \qquad (41) $$

while the second one is piecewise regular but contains a discontinuity curve (which corresponds to the circle $\tilde{r}(x, y) = 0.15$):

$$ \phi(0, x, y) = \begin{cases} 1.5 + \cos(\Pi \tilde{r}(x, y)) & \text{if } \tilde{r}(x, y) < 0.15 \\ 0 & \text{otherwise.} \end{cases} \qquad (42) $$

Convergence curves are plotted in figure 8 for both initial conditions, see respectively figure 8a for initial condition (41) and figure 8b for initial condition (42). Several remarks arise from these 2D results. First of all, it should be noticed that the extension of the GRU projection step to 2D case seems to perform in a satisfactory manner. Indeed, the approximated solutions converge to the exact solution with an effective rate of convergence close to 0.82 or 0.85. It should be noted that for the test case with initial condition (42), the Upwind scheme has not yet reached its asymptotic rate of convergence 1/2 on the finest meshes.

The behavior of the scheme is also illustrated by sets of figures 9-10 and by sets of figures 11-12. In the two former sets, approximated solutions have been plotted at $t_{end}$ for several meshes and for $(K_1, K_2) = (5, 3)$; and in the two latter sets the approximated solutions computed for a mesh with $400 \times 400$ cells has been plotted for several time-iterations and for $(K_1, K_2) = (13, 11)$. The sequence $(K_1, K_2) = (13, 11)$ produces a very accurate final approximated solution for the mesh $400 \times 400$, in particular it is more accurate than the
<figure><img src="image_15.png" /><figcaption>Figure 7: Examples of unstructured “twisted” meshes used for the 2D test case: 20 × 20 cells, 50 × 50 cells and 100 × 100 cells.</figcaption></figure>

one obtained with $(K_1, K_2) = (5, 3)$ (see figure 9).

When considering initial condition (41) and (42), the GRU step improves the accuracy of the approximated solutions with respect to the sole Upwind scheme. It can be observed in figure 9 (resp. 10), by comparing with figure 11 (resp. 12), that the quality of the approximated solutions strongly depends on the choice of the parameters $(K_1, K_2)$. Indeed, even if fluctuations are observed, the convergence rate is the same for all the couples of parameters. Nonetheless, the accuracy changes with the parameters. This can clearly be observed for instance in figure 8 when focusing on the finest meshes. For the sequences that have been used here, the difference between the best approximated solution and the worst one (with respect to the error $\varepsilon(\Delta x, K_1, K_2)$) can be greater than one order of magnitude.

<figure><img src="image_16.png" /><figcaption>Figure 8: Convergence curves for the 2D test cases for meshes from 20 × 20 cells to 1500 × 1500 cells.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 9: Approximated final solutions for the 2D test cases with (41) as initial condition and for $(K_1, K_2) = (5, 3)$. The first row of plots shows the approximated exact solutions (and approximated initial solutions), the second row shows the approximated solutions computed with the Upwind scheme, and the third row shows the approximated solutions computed with the GRU step. The first column corresponds to the mesh $100 \times 100$, the second one to the mesh $200 \times 200$, the third one to the mesh $400 \times 400$, the fourth one to the mesh $800 \times 800$ and the last one to the mesh $1500 \times 1500$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 10: Approximated final solutions for the 2D test cases with (42) as initial condition and for $(K_1, K_2) = (5, 3)$. The first row of plots shows the approximated exact solutions (and approximated initial solutions), the second row shows the approximated solutions computed with the Upwind scheme, and the third row shows the approximated solutions computed with the GRU step. The first column corresponds to the mesh $100 \times 100$, the second one to the mesh $200 \times 200$, the third one to the mesh $400 \times 400$, the fourth one to the mesh $800 \times 800$ and the last one to the mesh $1500 \times 1500$.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 11: Approximated solutions for the 2D test case with initial condition (41) and for $(K_1, K_2) = (13, 11)$ for the mesh $400 \times 400$. Several successive approximated solutions are plotted along one full rotation from figure (a) to (j).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 12: Approximated solutions for the 2D test case with initial condition (42) and for $(K_1, K_2) = (13, 11)$ for the mesh $400 \times 400$. Several successive approximated solutions are plotted along one full rotation from figure (a) to (j).</figcaption></figure>
# 4 Discussion on the choice of the low-discrepancy sequences

The scheme presented in section 2 is based on the uniform distribution $U(0, 1)$, whereas the numerical tests of section 3 have been performed using a family of deterministic low-discrepancy sequences. This is a classical approach used for Glimm's methods or for Monte-Carlo simulations, since the uniform distribution is a purely mathematical object that has to be approximated when turning to numerical simulations. The convergence of the Glimm's scheme has been studied in [11] when deterministic low-discrepancy sequences are used instead of the uniform distribution. It has been shown in [11] that:

* the sequences have to be equidistributed in the sense that the measure of discrepancy (see in section 4.1) has to tend to zero when the length of the sequences tend to infinity;

* the rate of convergence depends on how well the sequences are equidistributed.

In [7], a proof of convergence for a uniform distribution has been proposed for the GRU scheme, but there is no result available when pseudo-random or quasi-random sequences are used. In this section, a numerical investigation is carried in order to assess the counterpart for the GRU scheme of the theoretical result of [11] obtained for the Glimm's scheme. This numerical study is done on the basis of the test case of section 3.2 with $\Phi_s^0$ as an initial condition (see equation (34)) and by using several pseudo-random or quasi-random sequences for computing the approximated solutions of the corresponding initial condition problem.

## 4.1 Some quasi/pseudo-random sequences and simple metrics for assessing their quality

A wide range of low-discrepancy sequences and pseudo-random sequences have been proposed in the literature and it is not possible to be exhaustive here. Hence, we only consider in this section four classical low-discrepancy or pseudo-random sequences: the Halton-Van de Corput sequences (reader can refer to appendix 6.3 for details), the Hammersley sequences, the Linear Congruential Generator (nicknamed LCG in the following) and the standard pseudo-random generator of the C library (stdlib.h). It should be noted that the C pseudo-random generator is a LCG with specific choices of parameters that ensure that values of the sequences are not repeated. It will be nicknamed *rand*() in the following. Practical information on the computation of the Hammersley and Halton-Van de Corput sequences can be found for instance in [14], and the LCG method is widely studied in [12].

All these methods involve a set of parameters, different sets of parameters leading to different sequences. This point has already been noticed in section 3 for the Halton-Van de Corput sequences. In order to limit the number of computations, we restrict here to particular choices for these sets of parameters. The Hammersley sequence used in the following is the 2D sequence, $s = 2$, for $p_1 = 2$ when considering the notation of [14]. Only the first component of the sequence has been kept for our purpose. The LCG is used with the multiplier $\lambda = 1$, the modulus $m = 1$ and with the increment chosen as the golden ratio $r = (\sqrt{5} - 1)/2$. This set of parameters is known to lead to a sequence with a good quality for the LCG. At last, for the Halton-Van der Corput sequence, the classical parameters $k_1 = 5$ and $K_2 = 3$ have been retained here.

The first 100 elements of these sequences are plotted in figure 13. A first obvious remark is that for the *rand*() sequence the points visually seems less well distributed in the plot. Moreover, for the LCG and the Halton-Van der Corput sequences, some sets of points are almost aligned. Some objective metrics are then needed to estimate the quality of these sequences. We focus in this section on four criteria among many others [12, 6]. Let us assume that we have a sequence $X = (X_i)_{i=1..N}$ with $\forall i, X_i \in [0, 1]$. The average value and the variance of the sequence X are respectively defined as:

$$ \bar{X} = \frac{1}{N} \sum_{i=1}^{N} X_i \quad \text{and} \quad \sigma(X) = \frac{1}{N} \sum_{i=1}^{N} (X_i - \bar{X})^2. $$

It is recalled that the mean and the variance for the uniform distribution $U(0, 1)$ are respectively equal to $1/2$ and $1/12$. The following criteria are then defined for all these sequences.

* The mean criterion is: $\varepsilon_N^m(X) = |\bar{X} - \frac{1}{2}|$.

* The variance criterion is: $\varepsilon_N^v(X) = |\sigma(X) - \frac{1}{12}|$.

* The star-discrepancy reads:

$$ D_N^*(X) = \sup_{I \in \mathcal{I}(0,1)} \left| \frac{A(I, X)}{N} - |I| \right|, \quad (43) $$

where $\mathcal{I}(0, 1)$ corresponds to the set of all the intervals of the form $[0, b]$ with $0 \le b \le 1$. For any interval $I = [0, b]$, $|I| = b$ is its Lebesgue measure. The quantity $A(I, \bar{X})$ is the number of elements of the sequence X that belong to the interval I.

* The auto-correlation with one shift is defined if $\sigma(X) \neq 0$ as:

$$ Corr_N^1(X) = \frac{1}{\sigma(X)} \frac{\sum_{i=2}^{N} ((X_i - \bar{X}) \times (X_{i-1} - \bar{X}))}{N-1}. $$
We use here the star-discrepancy defined by (43) for the sake of simplicity. It is indeed more easy to compute that the discrepancy $D_N(X)$ for which the set of intervals $\mathcal{I}(0, 1)$ in definition (43) is replaced by the set of the intervals of the form $[a, b]$ with $0 \le a < b \le 1$. Moreover, some bounds on the discrepancy can be obtained from the star-discrepancy:

$$D_N^*(X) \le D_N(X) \le 4 D_N^*(X),$$

when considering one-dimensional sequences, see [6] for instance. A sequence is equidistributed if the discrepancy tends to zero as $N$ tends to $+\infty$.

These four criteria are plotted in figures 14 and 15 for sequences that contain up to 1000000 elements. For the Halton-Van der Corput, the Hammersley and the LCG sequences, the mean criterion and the variance criterion tend to zero when $N$ tends to $+\infty$ whereas it decreases very slowly for the sequences obtained by $rand()$. On the contrary, it should be noted that the three former sequences have almost constant auto-correlations which are respectively equal to $-0.452, -0.714$ and $-0.416$. They have thus not been plotted in figure 14. The difference between the $rand()$ sequence on the one hand and the three other sequences on the other hand is clearly exhibited by these figures. The $rand()$ sequence is a pseudo-random sequence that is built for providing a low auto-correlation. The latter clearly tends to zero when $N$ increases and it is always lower than for the three other sequences. This mimics the independence of two tries for a uniform distribution. On the contrary, the three other sequences are quasi-random sequences for which the equidistribution property is the most important feature. It can be observed in figure 15 that the $rand()$ sequence has a very large discrepancy with respect to the three other sequences. Up to $N \sim 120000$ (i.e. $\ln(120000) \sim 11.7$), the star-discrepancy diminishes with a small slope, and it remains always greater than that of the three quasi-random sequences. These three sequences behave similarly for the star-discrepancy. In particular, it is very low even for the short sequences (i.e. for small values of $N$), which is an important point (see also section 4.2).

<figure><img src="image_21.png" /><figcaption>Figure 13: First 100 elements of the sequences used in the present work.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure 14: Mean criterion, variance criterion and auto-correlation for the different sequences with respect to ln(N) (N is the length of the sequence). A ln-scale is used for the mean and variance criteria. The auto-correlation for the Halton-Van der Corput, Hammersley and LCG sequences are almost constant respectively equal to -0.452, -0.714 and -0.416. They have thus not been plotted here.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 15: Neperian logarithm of the star-discrepancy $D_N^*(X)$ for the different sequences X with respect to ln(N).</figcaption></figure>
## 4.2 Comparison of the results for several families of sequences

In the previous section, the quality of the different sequences have been examined. In this section, the influence of this quality on the simulation results with the GRU scheme is studied. For that purpose, we consider the test case of section 3.2 with a non-vanishing non-uniform velocity field (set of parameters (i)) and the function $\Phi_s^0$ defined by (34) as an initial condition. For this test case, the approximated solutions obtained by using the four sequences described in section 4.1 are compared in terms of the error $\varepsilon(\Delta x)$ defined by equation (36).

In figure 16, approximated solutions have been plotted for a mesh with 200 cells. It can be observed that the results obtained with the *rand()* sequences (for several seeds) is not as good as the results for the three other sequences. The convergence curves for the error $\varepsilon(\Delta x)$ is plotted in figure 17. It should be noted that for the *rand()* sequences and for the Halton-Van der Corput sequences, the convergence curve for a single sequence can fluctuate, as already shown in section 3. In figure 17, the convergence curves for these two families of sequences have been averaged over 28 sets of parameters in order to reduce these fluctuations. The convergence curve clearly highlight that *rand()* generator provides sequences with a poor quality for the GRU scheme. Both accuracy and convergence rate are degraded when using *rand()*: the classical upwind scheme even performs better than the GRU step associated to the upwind scheme. On the contrary, for the three other quasi-random sequences, an effective convergence rate close to 0.85 is observed together with a very good accuracy. These three sequences provide very similar accuracy.

In section 4.1, it has been remarked that the star-discrepancy for *rand()* remains high for sequences with less that ~ 120000 elements (i.e. see figure 15 with ln(120000) ~ 11.7) and is slightly improved for larger sequences. For the test case used here, due to CFL number of 0.45 and to the final time of the solution, a mesh that contains ~ 64000 cells requires ~ 120000 time-iterations and thus a sequence of ~ 120000 elements. Hence, to go beyond the limitation mentioned above for the discrepancy, the simulations with *rand()* have been performed up to 256000 cells (i.e. ln(256000) ~ 12.46). No noticeable improvement of the convergence curve is observed for these finest meshes. Indeed, even with a sequence containing more than 120000 elements, the 120000 first time-iterations are performed thanks to a sequence with a poor discrepancy. It seems that this “pollution” of the approximated solution on the first time-iterations impacts deeply the quality of the latter even afterwards. The Halton-Van der Corput, Hammersley and LCG sequences have a good star-discrepancy from the very first elements as shown by figure 15 (in figure 18 the star-discrepancy of Halton-Van der Corput sequences has been plotted on coarser meshes). This seems to be an important feature for the GRU scheme to be accurate.

<figure><img src="image_24.png" /><figcaption>Figure 16: Approximated solutions obtained for the test case with $\Phi_s^0$ as an initial condition and a non-vanishing non-uniform velocity field (set of parameters (i)).</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure 17: Convergence curves for several sequences. The results obtained with the Halton-Van de Corput and the *rand()* sequences have been averaged over several sets of parameters. The convergence curve for *rand()* has been computed up to 256000 cells.</figcaption></figure>
4.3 A short analysis of the results for the Halton-Van der Corput sequences

As mentioned in the sections above, the Halton-Van der Corput sequences exhibit fluctuations of the error
with respect to the mesh size and to the choice of the parameters involved in these sequences. In figure 18,
the results for each Halton-Van der Corput sequence in V(8) are plotted with respect to the length of the
sequence N, which is directly linked to the number of iterations and thus, through the CFL number and tend,
to the mesh size Δx. The error ε(Δx) is shown on the left and the star-discrepancy is shown on the right
for the same values of the length N. The decrease of the star-discrepancy with respect to N has a higher
slope than the decrease of ε(Δx). In that sense, the convergence of the sequence towards an equidistributed
sequence is faster than the convergence of the approximated solutions towards the exact solution when N
increases. Hence, the former may not be a limitation for the latter.

<figure><img src="image_26.png" /><figcaption>(a) Convergence curves for ε(Δx) with respect to the length of the sequence N.
(b) Convergence of the star-discrepancy with respect to the length of the sequence N.</figcaption></figure>

<caption>Figure 18: Convergence with respect to the length of the sequence *N* for the Halton-Van der Corput sequences of *V*(8).</caption>
# 5 Conclusion

In this work, an extension of the GRU step [10] has been proposed. It was first restricted to the advection of Heaviside functions and the extension tested here allows to apply this stochastic projection step to the advection of any scalar quantity. Some properties or behaviors of the GRU projection step have been studied when coupled to an Upwind-like scheme for the prediction step. In particular, it appears that this scheme has a very low level of numerical diffusion. One of the origin of this behavior lies in the fact that between two time-iterations the scheme does not create intermediate values. The approximated solution at iteration $n+1$ is built by choosing each cell-value among the cell-values of the approximated solution at iteration $n$. This particularity make the GRU scheme very close to the Glimm's scheme, even if the former is restricted to scalar advection. One major consequence and advantage of this feature is that the discontinuities in the approximated solutions are computed very accurately. On the other hand, one main drawback of the method could be that it is not conservative for a single realization (i.e. one quasi-random sequence), but it is only statistically conservative.

The numerical tests of section 3 clearly show that the GRU step largely improves the accuracy of the approximated solutions when compared to the sole Upwind scheme. As expected, a convergence rate close to 1 (the effective convergence rate is around 0.82, as in [10]) has been retrieved for piecewise regular functions involving discontinuities. This rate of convergence is higher than that of more classical first-order finite-volume scheme, which is in general around 1/2 for weak solutions. The interest of the whole scheme is meaningful if one consider that the implementation of the GRU step is very easy. It corresponds to a projection step that can be added after an existing prediction step, and that only involves a cell and its neighboring cells.

From section 4 it arises that the sequences used in a numerical point of view mimics the uniform distribution must be low-discrepancy sequences. For these quasi-random sequences the low level of discrepancy simulates the independance of the uniform distribution even if they exhibit a high level of auto-correlation. This has been numerically observed for the GRU scheme in section 4 and a theoretical result of convergence with deterministic sequences should be investigated on the basis of [7]. Such a theoretical result was obtained for the Glimm's scheme in [11].

At last, let us mention some further investigations. First, if a proof of convergence has been proposed in [7] in the case of Heaviside functions, it should be possible to extend this proof to the scalar case, at least for simple one-dimensional configurations. Furthermore, one of the feature of the GRU step with respect to the Glimm's scheme is that it can be applied for 1D, 2D or 3D domains. Actually, the GRU step extends easily to 2D and 3D domains and its implementation remains very easy even for multi-processor computing.

# Acknowledgement

Computational facilities were provided by EDF R&D and by the NEEDS-M2SIR project.

# 6 Appendices

## 6.1 A more general result for the statistical consistency of the GRU step

Let us assume that the prediction step is achieved thanks to a three-point scheme that can be written:

$$ \phi_i^{n+1,*} = G_i(\phi_{i-1}^n, \phi_i^n, \phi_{i+1}^n). $$

This prediction step is supposed to be such that $\phi_i^{n+1,*} \in [\phi_{i,m}^n, \phi_{i,M}^n]$. Let us also assume that the random number $w^n$ is chosen accordingly to a probability density function $w \in [0, 1] \mapsto f(w)$, and that the threshold projection in cell $i$ is still defined by (9). The projection step then gives the transitional probabilities (11) with:

$$ \alpha_i^{n,*} = \int_{\delta_i^n}^{1} f(x)dx = 1 - \int_{0}^{\delta_i^n} f(x)dx, \qquad (44) $$

where:

$$ \delta_i^n = \frac{\phi_i^{n+1,*} - \phi_{i,m}^n}{\phi_{i,M}^n - \phi_{i,m}^n} \quad \text{if} \quad \phi_{i,m}^n \neq \phi_{i,M}^n, \quad \text{and} \quad \delta_i^n = 1 \text{ otherwise.} $$

A necessary condition to get the statistical consistency is that the expectation of $\phi_i^{n+1}$ knowing $(\phi_j^n)_j$ has to be equal to $\phi_i^{n+1,*}$:

$$ E(\phi_i^{n+1}|(\phi_j^n)_j) = \phi_i^{n+1,*}. \qquad (45) $$

Thanks to the definition of the GRU step, we have here:

$$ E(\phi_i^{n+1}|(\phi_j^n)_j) = \alpha_i^{n,*} \phi_{i,m}^n + (1 - \alpha_i^{n,*}) \phi_{i,M}^n, \qquad (46) $$
therefore, using relations (45) and (46), the necessary condition to get the statistical consistency reads:

$$ \alpha_{i}^{n,*} \phi_{i,m}^{n} + (1 - \alpha_{i}^{n,*}) \phi_{i,M}^{n} = \phi_{i}^{n+1,*}, \quad (47) $$

Let us now assume that $\phi_{i,m}^n \neq \phi_{i,M}^n$, so that the equality (47) implies that:

$$ \alpha_{i}^{n,*} = 1 - \delta_{i}^{n}, $$

or using definition (44) for the quantity $\alpha_i^{n,*}$:

$$ \delta_{i}^{n} = \int_{0}^{\delta_{i}^{n}} f(x) dx. \quad (48) $$

A first important result is that the only density probability $f$ that fulfills (48) for all $\delta_i^n \in [0, 1]$ is the density probability of the uniform distribution $U(0, 1)$, that is $f(x) = 1 \forall x \in [0, 1]$. We thus assume from now that $f$ denotes the density function for $U(0, 1)$, so that we have:

$$ E(\phi_i^{n+1} | (\phi_j^n)_j) = \phi_i^{n+1,*} = \mathcal{G}_i(\phi_{i-1}^n, \phi_i^n, \phi_{i+1}^n) \quad (49) $$

By integrating expectation (49) over all the possible sequences $(\phi_j^n)_j$, we obtain:

$$ E(\phi_i^{n+1}) = E(\mathcal{G}_i(\phi_{i-1}^n, \phi_i^n, \phi_{i+1}^n)) $$

If the function $(x, y, z) \mapsto \mathcal{G}_i(x, y, z)$ is linear for all $i$, we then have:

$$ E(\phi_i^{n+1}) = \mathcal{G}_i(E(\phi_{i-1}^n), E(\phi_i^n), E(\phi_{i+1}^n)), \quad (50) $$

which gives the statistical consistency of the projection step with the prediction step. The following result can thus be stated if the prediction step is achieved using a three point linear scheme that preserves monotoncity: a necessary condition for the GRU step to allow to get the statistical consistency with the prediction step in the sense of relation (50) is that it has to be performed using a uniform distribution on (0, 1).

The condition stated above is necessary and not sufficient. The particular case $\phi_{i,m}^n = \phi_{i,M}^n$ is not studied here in details, but it should be noted that it implies additional constraint for the statistical consistency. In particular, the scheme $\mathcal{G}_i$ has to be defined in agreement with the GRU step. For instance, the statistical consistency can be lost when using the GRU step defined in section 2.1.2 and the Lax-Friedrichs scheme for the prediction step.

## 6.2 A class of analytical solutions for the 1D convection problem with non-uniform velocity and density

We propose here a very simple class of analytical solution for system (1) on the domain $(t, x) \in [0, T_0] \times [0, 1]$, with a bounded final time $T_0 > 0$. This class of solutions possesses non-uniform but regular profiles for the density (which is non-negative) and for the velocity, and the scalar field $\phi$ can be piecewise continuous. When considering regular density and velocity, the first equation of (1) can be written:

$$ \partial_t (\ln(\rho)) + U \partial_x (\ln(\rho)) + \partial_x (U) = 0. \quad (51) $$

Since system (1) is based on two equations and three unknowns, we can specify one among this unknown. If the density $\rho$ is chosen such that:

$$ \rho(t, x) = e^{(R_0 (x+B_0 t))}, \quad (52) $$

where $R_0$ and $B_0$ are uniform and constant, then through equation (51) the velocity must fulfill:

$$ \partial_x (U) + R_0 (U + B_0) = 0. \quad (53) $$

Equation (53) can obviously be exactly integrated on $[0, T_0] \times [0, 1]$. It yields:

$$ U(t, x) = U(t, x=0) e^{(-R_0 x)} - B_0 \left(1 - e^{(-R_0 x)}\right). \quad (54) $$

Let us assume that the velocity $U(t, x=0)$ does not depend on time: $U(t, x=0) = U_0$, so that the velocity (54) does not depend on time. If $(1 + U_0/B_0) > 0$, the velocity vanishes for a unique abscissa $\tilde{x}$ defined by:

$$ 0 = U_0 e^{(-R_0 \tilde{x})} - B_0 \left(1 - e^{(-R_0 \tilde{x})}\right) \iff \tilde{x} = \frac{1}{R_0} \ln \left(1 + \frac{U_0}{B_0}\right). \quad (55) $$

We also assume that the initial profile of the scalar field is defined by a function $\phi_0$ which is piecewise continuous:

$$ \phi(t=0, x) = \phi_0(x), \quad (56) $$

According to system (1), and since both $\rho$ and $U$ are regular, the solution for $\phi$ with initial condition (56) is advected with the velocity $U$. Let us consider a fluid particle which is located at $x_0$ at time $t=0$. Provided
that the velocity does not vanish for $t \in [0, t^*]$, this fluid particle is located at time $t^*$ at the abscissa $x^*$ given by the relation:

$$t^* = \int_{x_0}^{x^*} \frac{dx}{U(x)}. \qquad (57)$$

It should be emphasized that relation (57) holds because the velocity does not depend on time. The velocity (54) can be written in the form: $U = b + e^{(ax)}$, and the integral in (57) can thus be exactly expressed. Indeed, we have:

$$\int_{x_1}^{x_2} \frac{dx}{b + e^{(ax)}} = \frac{-1}{ab} \ln \left( \frac{1 + be^{(-ax_2)}}{1 + be^{(-ax_1)}} \right),$$

which leads through (57) to the following expression of $x_0$ as a function of $t^*$ and $x^*$:

$$x_0(t^*, x^*) = x^* + B_0 t^* + \frac{1}{R_0} \ln \left( 1 - \left(1 + U_0/B_0\right) e^{(-R_0 x^*)} \left(1 - e^{(-R_0 B_0 t^*)}\right) \right). \qquad (58)$$

Hence, if $U_0$ and $B_0$ are chosen so that the velocity never vanishes, the solution for the scalar field is:

$$\phi(t, x) = \phi_0(x_0(t, x)), \qquad (59)$$

where the position $x_0(t, x)$ is given by relation (58). When the velocity field vanishes at a point $\tilde{x}$, solution (59) defined by characteristic curve (58) has to be defined following the same way but by considering separately the domains $x < \tilde{x}$ and $x > \tilde{x}$.

In figure (19), three different analytical solutions are plotted on $[0, 1]$. These solutions are used in section (3.2) for assessing the behavior of the numerical schemes. The first two solutions on the left of figure (19) correspond to solution with a vanishing velocity at $\tilde{x} = 1/2$.

<figure><img src="image_27.png" /><figcaption>Figure 19: Three analytical solutions. Left: dilatation of the scalar profile with a vanishing velocity at $\tilde{x} = 1/2$, middle: compression of the scalar profile with a vanishing velocity at $\tilde{x} = 1/2$, right: compression of the scalar profile for a non vanishing velocity.</figcaption></figure>

## 6.3 Computing a Halton-Van der Corput sequence

In the numerical tests the random number $\omega_n$ introduced in the GRU step has been replaced by using $(K_1, K_2)$ Halton-Van der Corput sequences. These are low discrepancy sequences on $(0, 1)$. The parameters of the sequence, $K_1$ and $K_2$, are two integers relatively prime and such that $K_1 > K_2 \ge 3$. The nth element, $\omega_n$, of the sequence is computed using the following algorithm [13]:

$$\omega_n = \sum_{i=0}^{m} A_i K_1^{-(i+1)},$$
with:

$$A_i = rem(K_2a_i, K_1) \quad \text{and} \quad n = \sum_{i=0}^{m} a_i K_i,$$

and where $rem(a, b)$ is the remainder of the division of $a$ by $b$.
References

[1] A. J. CHORIN, *Random choice solution of hyperbolic systems*, Journal of Computational Physics, 22 (1976), pp. 517–533.

[2] P. COLELLA, *Analysis of the effect of operator splitting and of the sampling procedure on the accuracy of glimm's method*, Report No. LBL-8774. California Univ., Berkeley (USA). Lawrence Berkeley Lab., (1978).

[3] ——, *Glimm's method for gas dynamics*, SIAM Journal on Scientific and Statistical Computing, 3 (1982), pp. 76–110.

[4] F. DELARUE AND F. LAGOUTIÈRE, *Probabilistic analysis of the upwind scheme for transport equations*, Archive for rational mechanics and analysis, 199 (2011), pp. 229–268.

[5] B. DESPRÉS AND F. LAGOUTIÈRE, *Un schéma non linéaire anti-dissipatif pour l'équation d'advection linéaire*, Comptes Rendus de l'Académie des Sciences-Series I-Mathematics, 328 (1999), pp. 939–943.

[6] M. DRMOTA AND R. F. TICHY, *Sequences, discrepancies and applications*, Springer, 2006.

[7] T. GALLOUËT AND O. HURISSE, *Convergence of a multidimensional glimm-like scheme for the transport of fronts*, IMA Journal of Numerical Analysis, (2021).

[8] J. GLIMM, *Solutions in the large for nonlinear hyperbolic systems of equations*, Communications on pure and applied mathematics, 18 (1965), pp. 697–715.

[9] A. HARTEN, P. D. LAX, AND B. V. LEER, *On upstream differencing and godunov-type schemes for hyperbolic conservation laws*, SIAM review, 25 (1983), pp. 35–61.

[10] O. HURISSE, *On the use of glimm-like schemes for transport equations on multidimensional domain*, International Journal for Numerical Methods in Fluids, 93 (2021), pp. 1235–1268.

[11] T.-P. LIU, *The deterministic version of the glimm scheme*, Communications in Mathematical Physics, 57 (1977), pp. 135–148.

[12] H. NIEDERREITER, *Quasi-monte carlo methods and pseudo-random numbers*, Bulletin of the American mathematical society, 84 (1978), pp. 957–1041.

[13] F. TORO ELEUTERIO, *Riemann solvers and numerical methods for fluid dynamics: A practical introduction*, Springer, Berlin, (1997).

[14] T.-T. WONG, W.-S. LUK, AND P.-A. HENG, *Sampling with hammersley and halton points*, Journal of graphics tools, 2 (1997), pp. 9–24.