domain: phys spi

text:
Nour Dahmen
email: nourallah.dahmen@onera.fr

Jérôme Droniou

François Rogier

A cost-effective nonlinear extremum-preserving finite volume scheme for highly anisotropic diffusion on Cartesian grids, with application to radiation belt dynamics

Keywords: Radiation belts, Finite volume method, Anisotropic diffusion equation, Monotonicity, Discrete maximum principle

We construct a new nonlinear finite volume (FV) scheme for highly anisotropic diffusion equations, that satisfies the discrete minimum-maximum principle. The construction relies on the linearized scheme satisfying less restrictive monotonicity conditions than those of an M -matrix, based on a weakly regular matrix splitting and using the Cartesian structure of the mesh (extension to quadrilateral meshes is also possible). The resulting scheme, obtained by expressing fluxes as nonlinear combinations of linear fluxes, has a larger stencil than other nonlinear positivity preserving or minimum-maximum principle preserving schemes. Its larger "linearized" stencil, closer to the actual complete stencil (that includes unknowns appearing in the convex combination coefficients), enables a faster convergence of the Picard iterations used to compute the solution of the scheme.

Steady state dimensionless numerical tests as well as simulations of the highly anisotropic diffusion in electron radiation belts show a second order of convergence of the new scheme and confirm its computational efficiency compared to usual nonlinear FV schemes.

Introduction

Finite volume methods have been favored for a long time in engineering communities for solving diffusion equations resulting from fluid flows. Unlike finite difference methods, they are specifically designed to provide a local conservativity property to the numerical resolution, fundamental to conservation laws [START_REF] Eymard | Finite volume methods[END_REF]. However, important numerical properties related to accuracy, computational efficiency and physical relevance are also required to ensure a reliable numerical solution [START_REF] Droniou | Finite volume schemes for diffusion equations: introduction to and review of modern methods[END_REF].

One of the most sought after properties is monotonicity, which embodies (for linear schemes) the respect of the discrete minimum-maximum principle (MMP) and the preservation of the positivity of the numerical solution [START_REF] Droniou | Finite volume schemes for diffusion equations: introduction to and review of modern methods[END_REF]. Monotonicity is obtained by ensuring a monotonic discretization matrix. However, this property is hardly achieved by linear FV schemes for any mesh type or in case of highly anisotropic diffusion problems (with anisotropy ratios > 10 6 ), typically encountered in applications like reservoir engineering, radionuclide migration or radiation belt dynamics. For the last application precisely, the numerical constraints are sometimes so high that an incomplete diffusive frame is considered to prevent negative values and spurious non-physical oscillations from occurring [START_REF] Varotsou | Three-dimensional test simulations of the outer radiation belt electron dynamics including electron-chorus resonant interactions[END_REF][START_REF] Dahmen | On the modelling of highly anisotropic diffusion for electron radiation belt dynamic codes[END_REF][START_REF] Dahmen | Développement d'un nouveau coeur numérique pour le code de calcul salammbô de modélisation des ceintures de radiation terrestres[END_REF].

On quadrilateral grids, linear schemes with 9-point stencils cannot be monotonic in case of highly anisotropic diffusion tensor [START_REF] Kershaw | Differencing of the diffusion equation in Lagrangian hydrodynamic codes[END_REF][START_REF] Buet | On the non existence of monotone linear schema for some linear parabolic equations[END_REF]; recovering monotonicity requires the usage of sometimes unreasonably large stencils [START_REF] Potier | Un schéma linéaire vérifiant le principe du maximum pour des opérateurs de diffusion très anisotropes sur des maillages déformés[END_REF] (which comes with a computational cost). To preserve local stencils and the monotonicity of the scheme, it is necessary to abandon the principle of linearity of the scheme and adopt a flexible nonlinear "monotonic" construction based on a convex combination of linear upwind fluxes [START_REF] Droniou | Finite volume schemes for diffusion equations: introduction to and review of modern methods[END_REF][START_REF] Schneider | Comparison of finite-volume schemes for diffusion problems[END_REF], inspired by [START_REF] Bertolazzi | A second-order maximum principle preserving finite volume method for steady convection-diffusion problems[END_REF]. The convex combination weights are properly chosen to either ensure extremum principle preserving properties, or at least non negative numerical solutions [START_REF] Droniou | Construction and convergence study of schemes preserving the elliptic local maximum principle[END_REF][START_REF] Sheng | The finite volume scheme preserving extremum principle for diffusion equations on polygonal meshes[END_REF][START_REF] Sheng | A nonlinear convex combination in the construction of finite volume scheme satisfying maximum principle[END_REF][START_REF] Potier | Correction non linéaire et principe du maximum pour la discrétisation d'opérateurs de diffusion avec des schémas volumes finis centrés sur les mailles[END_REF][START_REF] Yuan | Monotone finite volume schemes for diffusion equations on polygonal meshes[END_REF][START_REF] Potier | A nonlinear finite volume scheme satisfying maximum and minimum principles for diffusion operators[END_REF][START_REF] Cancès | Monotone correction for generic cell-centered finite volume approximations of anisotropic diffusion equations[END_REF][START_REF] Lipnikov | Interpolation-free monotone finite volume method for diffusion equations on polygonal meshes[END_REF][START_REF] Lipnikov | Minimal stencil finite volume scheme with the discrete maximum principle[END_REF][START_REF] Lipnikov | Monotone finite volume schemes for diffusion equations on unstructured triangular and shape-regular polygonal meshes[END_REF]. This usually amounts to ensuring that a linearized version of the scheme has a so-called M -matrix structure, which satisfies the following properties:

• Irreducible with strict positive diagonal entries,

• Non negative diagonal entries (Z-matrix),

• Weak diagonal dominance (strict for at least one row).

However, turning to nonlinear schemes impose a trade-off as the positivity preserving principle and the MMP are no longer equivalent. Some issues can also arise when needing to eliminate intermediate unknowns in the construction, several options being available but not all of them leading to preserving the monotonicity of the scheme [START_REF] Sheng | A nonlinear convex combination in the construction of finite volume scheme satisfying maximum principle[END_REF][START_REF] Agelas | A nine-point finite volume scheme for the simulation of diffusion in heterogeneous media[END_REF][START_REF] Sheng | The finite volume scheme preserving extremum principle for diffusion equations on polygonal meshes[END_REF][START_REF] Sheng | Construction of nonlinear weighted method for finite volume schemes preserving maximum principle[END_REF]. Besides, nonlinear resolution introduces an additional computational cost that can quickly become substantial, especially in the context of radiation belt diffusion, as reported in recent studies [START_REF] Dahmen | On the modelling of highly anisotropic diffusion for electron radiation belt dynamic codes[END_REF].

Starting from the general quest of designing robust numerical schemes for highly anisotropic and inhomogeneous diffusion equations, and the specific considerations imposed by radiation belts dynamics that imperatively requires physically relevant solutions and optimal resolution methods, we design a new nonlinear FV full-monotonicity preserving scheme on non-uniform Cartesian grids (as used in radiation belts simulations [START_REF] Varotsou | Three-dimensional test simulations of the outer radiation belt electron dynamics including electron-chorus resonant interactions[END_REF][START_REF] Subbotin | Three-dimensional modeling of the radiation belts using the versatile electron radiation belt (verb) code[END_REF][START_REF] Su | Steerb: A three-dimensional code for storm-time evolution of electron radiation belt[END_REF]). The main innovation of our design, compared to usual nonlinear FV scheme, is the usage of convex combinations of upwind linear fluxes with a larger stencil. Precisely, the latter matches in its "linearized form" the complete nonlinear stencil accounting for the unknowns appearing in the convex combination coefficients. As a consequence, this feature leads to an important reduction of the number of nonlinear iterations required to compute the scheme's solution. The usage of a larger stencil prevents the linearized schemes to have an M -matrix. The monotonicity of our improved method is however preserved using the conditions described in [START_REF] Nordbotten | Monotonicity of control volume methods[END_REF], which are less restrictive than the conditions imposed by an M -matrix construction, and are based on a weakly regular splitting of the discretization matrix. Thanks to this construction, the new scheme preserves both positivity and the MMP, and numerically achieves a second order rate of convergence. The outline of this paper is as follows. In section 2, we describe the studied diffusion equation, the FV mesh and the construction of typical nonlinear FV schemes found in the literature.

The construction of the new scheme, called Relaxed NonLinear MultiPoint Flux Approximation (R-NLMPFA) due to its usage of relaxed monotonicity conditions, is the purpose of Section 3, in which we also recall the alternative monotonicity conditions underpinning it. In Section 4, numerical results are presented on steady state dimensionless cases, assessing the numerical properties of the new scheme compared to typical nonlinear schemes, with a special attention paid to the computational cost. In section 5, we succinctly present our application context of electron radiation belt modelling. The latter relies on a highly inhomogenous, anisotropic diffusion equation inducing huge numerical constraints on finite difference based codes. Then, the R-NLMPFA is tested over a practical case of electron radiation belt simulation to confirm its relevance. Finally, we summarize our work in section 6 and give some conclusions.

Nonlinear finite volume schemes for anisotropic diffusion equations

Diffusion equation and finite volume construction on a Cartesian mesh

We consider the following general steady state diffusion problem

   ∇ • (D∇f ) + S = 0 in Ω, f = f on ∂Ω, (1)
with Ω ⊂ R 2 a rectangular domain (extension to 3D is possible) and ∂Ω its boundary on which Dirichlet boundary conditions are imposed. D : Ω → M 2 (R) is the diffusion tensor (uniformly symmetric definite positive and bounded matrix-valued) and S ∈ L 2 (Ω) is the source term.

The finite volume framework is based on the local balance of fluxes, obtained after integrating the PDE in (1) over an elementary volume K ⊂ Ω, and introducing the flux F K,σ crossing through each edge σ ∈ ∂K, see [START_REF] Eymard | Finite volume methods[END_REF]:

K SdV = - ∂K D∇f • n K dl = σ∈∂K F K,σ (2)
with

F K,σ = - σ D∇f • n K,σ dl. (3)
Here, n K is the outer normal vector to ∂K, and n K,σ is the normal vector to the edge σ pointing outside K.

In FV methods [START_REF] Droniou | Finite volume schemes for diffusion equations: introduction to and review of modern methods[END_REF], F K,σ is approximated by a consistent expression F K,σ on a set T of disjoint mesh cells covering Ω, that is Ω = K∈T K. Each cell K ∈ T is is equipped with a point c K of coordinates (x K , y K ) (often referred to as its "center", although it does not have to be the center of mass) and associated to a set E K of edges such that ∂K = σ∈E K σ; we denote by |σ| the measure of σ. E ext is the set of exterior edges such that ∂Ω = σ∈Eext σ, S K denotes the approximation of S on the cell K, f σ is the approximation of f on the edge σ and D(K) the diffusion tensor evaluated in K (this evaluation could be at the center of the cell, or the average over the cell).

For K and L, two neighbouring cells in the mesh T and sharing the edge σ, FV methods impose by definition (see ( 3)) the conservativity property

F K,σ + F L,σ = 0. ( 4
)
The second key ingredient in the design of a FV method is to impose at the discrete level the flux balance [START_REF] Droniou | Finite volume schemes for diffusion equations: introduction to and review of modern methods[END_REF], that is: for all cell K,

σ∈E K F K,σ = S K .
As this research work is primarily related to the electron radiation belt numerical modelling application, we consider from now on a Cartesian grid case, i.e., a rectangular and conforming but possibly nonuniform mesh as represented in Figure 1. In fact, this grid type is used in the most recent application of the finite volume method in the electron radiation belt context within ONERA's Salammbô code [START_REF] Varotsou | Three-dimensional test simulations of the outer radiation belt electron dynamics including electron-chorus resonant interactions[END_REF][START_REF] Dahmen | On the modelling of highly anisotropic diffusion for electron radiation belt dynamic codes[END_REF][START_REF] Dahmen | Développement d'un nouveau coeur numérique pour le code de calcul salammbô de modélisation des ceintures de radiation terrestres[END_REF][START_REF] Bourdarie | Electron radiation belt data assimilation with an ensemble kalman filter relying on the salammbô code[END_REF]; it is also currently used in the finite-difference based codes of the radiation belts community [START_REF] Subbotin | Three-dimensional modeling of the radiation belts using the versatile electron radiation belt (verb) code[END_REF][START_REF] Su | Steerb: A three-dimensional code for storm-time evolution of electron radiation belt[END_REF][START_REF] Glauert | Three-dimensional electron radiation belt simulations using the bas radiation belt model with new diffusion models for chorus, plasmaspheric hiss, and lightning-generated whistlers[END_REF], so designing a FV method specifically adapted to such a grid can ease transition for said community. As a consequence of the choice of a Cartesian grid, we can easily adjust the cell centers so that lines joining them are parallel to the axes.

c K (x K , y K ) c L c M K c M ′ L c M ′ K c M L σ (x σ , y σ ) σ M K σ M ′ K σ M ′ L σ M L e x e y
F K,σ = µ 1 F 1 + µ 2 F 2 . (5)
The weights µ 1 and µ 2 of the convex combination are such as µ 1 , µ 2 ≥ 0 and µ 1 +µ 2 = 1; their choice is made to either ensure positivity (commonly called monotonic scheme) or the MMP (extremum principle preserving scheme).

F 1 and F 2 are associated respectively to K and L; their expressions are based on equation [START_REF] Varotsou | Three-dimensional test simulations of the outer radiation belt electron dynamics including electron-chorus resonant interactions[END_REF] and defined inside their corresponding cell. Since D∇f • n = Dn • ∇f (as D is symmetric), this yields the following approximations

F 1 = -D(K)n K,σ • (∇f ) K |σ| and F 2 = -D(L)n L,σ • (∇f ) L |σ|. (6)
Recalling that we use a Cartesian grid, consider the case n K,σ = e x and the diffusion tensor expressed as

D =   D xx D xy D xy D yy   . (7)
We aim to define (∇f ) K (respectively (∇f ) L ) according to the direction of the vector Dn K,σ (respectively Dn L,σ ) and using only the unknowns from cells neighbouring K and whose closure touches σ,

that is f K , f σ , f σ M K , f σ M ′ K , following the notations in Figure 1 (respectively f L , f σ , f M L , f M ′ L for the cell L).
The approximations of ∇f in each cell are "upwind" (thus leading to what we referred to as "upwind fluxes") in the sense that they are adapted to the directions of Dn in each cell (see [START_REF] Potier | Un schéma linéaire vérifiant le principe du maximum pour des opérateurs de diffusion très anisotropes sur des maillages déformés[END_REF])

and to the sign of D xy , with an appropriate choice of transverse edge unknowns. The expressions in (6) become:

F 1 =              -D xx (K) f σ -f K |x σ -x K | + D xy (K) f σ M K -f K |y σ M K -y K | |σ| if D xy (K) ≥ 0, -D xx (K) f σ -f K |x σ -x K | + (-D xy (K)) f σ M ′ K -f K |y σ M ′ K -y K | |σ| if D xy (K) < 0, F 2 =              --D xx (L) f σ -f L |x σ -x L | -D xy (L) f σ M L -f L |y σ M L -y L | |σ| if D xy (L) ≥ 0, --D xx (L) f σ -f L |x σ -x K | -(-D xy (L)) f σ M ′ L -f L |y σ M ′ L -y L | |σ| if D xy (L) < 0, (9)
where, for any edge σ, (x σ , y σ ) are the coordinates of the intersection of σ with the line (orthogonal to σ) joining the centers of the cells on each side of σ (see Figure 1).

The approximations of ∇f can also be recast with cell center unknowns, which yields the following alternative forms of F 1 and F 2 :

F 1 =            -D xx (K) f L -f K |x L -x K | + D xy (L) f M K -f K |y M K -y K | |σ| if D xy (K) ≥ 0, -D xx (K) f L -f K |x L -x K | + (-D xy (K)) f M ′ K -f K |y M ′ K -y K | |σ| if D xy (K) < 0, F 2 =            --D xx (L) f K -f L |x L -x K | -D xy (L) f M L -f L |y M L -y L | |σ| if D xy (L) ≥ 0, --D xx (L) f K -f L |x L -x K | -(-D xy (L)) f M ′ L -f L |y M ′ L -y L | |σ| if D xy (L) < 0. (10)
A first scheme that can be constructed using the flux-weighing approach is the NonLinear Two-Point Flux approximation scheme (NLTPFA) [START_REF] Droniou | Finite volume schemes for diffusion equations: introduction to and review of modern methods[END_REF][START_REF] Schneider | Comparison of finite-volume schemes for diffusion problems[END_REF][START_REF] Yuan | Monotone finite volume schemes for diffusion equations on polygonal meshes[END_REF][START_REF] Potier | Schéma volumes finis monotone pour des opérateurs de diffusion fortement anisotropes sur des maillages de triangles non structurés[END_REF]. This scheme achieves with the expressions ( 5) and ( 9) a two-point flux structure 4). The two-point structure is obtained through a special choice of µ 1 and µ 2 that eliminates the edge unknowns (see Appendix A). It provides a discretization matrix (once the nonlinear coefficients α K,L and β K,L are frozen)

F K,σ = α K,L (f )f K -β K,L (f )f L , (11)
with α K,L (f ), β K,L (f ) ≥ 0, and α K,L (f ) = β L,K (f ) to ensure (
that has an M -matrix structure [START_REF] Droniou | Finite volume schemes for diffusion equations: introduction to and review of modern methods[END_REF], which only ensures the positivity of the solution.

To obtain the full monotonicity, a Local Maximum Principle (LMP) structure is required [START_REF] Droniou | Construction and convergence study of schemes preserving the elliptic local maximum principle[END_REF].

In the context of nonlinear FV schemes, this is achieved with the NonLinear Multi-Point Flux Approximation scheme (NLMPFA) [START_REF] Droniou | Construction and convergence study of schemes preserving the elliptic local maximum principle[END_REF][START_REF] Droniou | Finite volume schemes for diffusion equations: introduction to and review of modern methods[END_REF][START_REF] Schneider | Comparison of finite-volume schemes for diffusion problems[END_REF][START_REF] Sheng | The finite volume scheme preserving extremum principle for diffusion equations on polygonal meshes[END_REF][START_REF] Cancès | Monotone corrections for generic cell-centered finite volume approximations of anisotropic diffusion equations[END_REF], whose fluxes are

F K,σ = N τ K,N (f K -f N ), (12)
with discretization weights τ K,N ≥ 0 (> 0 if N is the cell on the other side of σ from K), and the sum usually extended to the cells N that share an edge or vertex with K. In order to obtain this cell-centered structure for F K,σ , we use the cell-centered versions of F 1 and F 2 in (10) and we follow the construction reported in Appendix A.

From both schemes, a nonlinear algebric system is obtained that can be written

A(X)X = B ( 13
)
with A the discretization matrix of nonlinear scheme, X = (f K ) K∈T the unknown vector and B the right-hand side vector containing the source terms S K and the boundary terms f σ . The previous Algorithm 1 Picard algorithm for solving [START_REF] Sheng | A nonlinear convex combination in the construction of finite volume scheme satisfying maximum principle[END_REF] 1: Initialize X 0 ≥ 0 and s ← 0

K L M K M ′ L M ′ K M L σ σ M K σ M ′ K σ M ′ L σ M L e x e y
2: while ||X s+1 -X s ||∞ ||X s ||∞ ≥ ϵ do 3:
Construct the linearized system A(X s )X = B(X s ) (after freezing the nonlinear terms of the scheme)

4:

Solve A(X s )X s+1 = B(X s ) 5:

s ← s + 1

6: end while equation is solved using the Picard algorithm (Algorithm 1), with ϵ representing the stopping criteria Remark 1 (Stencil of linearized NLTPFA and NLMPFA). After freezing the nonlinear coefficients, both schemes have a linearized stencil of 5 points, which thus does not develop on all unknowns involved in the definitions of the linear fluxes F 1 and F 2 (which would lead to a 9-point stencil).

In fact, the NLTPFA scheme considers only K and L as unknowns in each discretization direction.

The NLMPFA scheme, on the other hand, considers only the cell center unknowns involved in F 1 , that is K, L and M K or M ′ K (see Figure 2), while F 2 and its cell center unknowns are only used to estimate the convex combination weights during the nonlinear resolution phase and do not enter explicitly in the linearized form of F K,σ (see Appendix A).

3. Nonlinear monotonic FV scheme with improved stencil

Monotonicity conditions

As previously presented, the monotonicity of nonlinear FV schemes is usually based on the Mmatrix construction argument, imposing a constrained 5-point linearized stencil of discretization on both NLTPFA and NLMPFA schemes. The new scheme that we present here is based on the weaker monotonicity conditions of [START_REF] Nordbotten | Monotonicity of control volume methods[END_REF], derived from the weakly regular splitting of the discretization matrix [START_REF] Berman | Nonnegative matrices in the mathematical sciences[END_REF].

Based on this decomposition, [START_REF] Nordbotten | Monotonicity of control volume methods[END_REF] establishes a set of conditions on the discretization matrix terms to ensure its monotonicity. Precisely, let us consider a discretization matrix of a nine-point locally conservative linear scheme, solving (1) on a two dimensional quadrilateral grid and that can be written as

9 k=1 m i,j k f k = Q i,j for all (i, j). (14)
Here, m i,j k represents the discretization weight of the 9-point stencil scheme, associated to the neighbouring cell k (as indexed in Figure 3) of the grid cell (i, j), and Q i,j is constructed from the source term and boundary conditions.

Assume a global ordering of cells starting from the lower left corner of the mesh grid (i = 1,j = 1) and running over each grid row from bottom to top. When reported on the discretization matrix A, the terms m i,j k refers to the matrix elements disposed in the row of cell (i, j) (so m i,j 1 is the diagonal element, m i,j 2 and m i,j 6 the first secondary diagonal terms of rank 1 and -1, etc.). The monotonicity conditions presented in [START_REF] Nordbotten | Monotonicity of control volume methods[END_REF] consists in the following set of 10 inequalities imposed on m i,j k :

1 (i, j)
m i,j 1 > 0, (A0) m i,j 2 < 0, (A1a) m i,j 6 < 0, (A1c)
m i,j 1 + m i,j 2 + m i,j 6 > 0, (A2) m i,j 4 < 0, (A1b)
m i,j 8 < 0, ( A1d
)
m i,j 2 m i,j-1 4 -m i,j-1 3 m i,j 1 > 0, ( A3a
)
m i,j 6 m i,j-1 4 -m i,j-1 5 m i,j 1 > 0, ( A3b
)
m i,j 2 m i,j+1 8 -m i,j+1 9 m i,j 1 > 0, ( A3c
)
m i,j 6 m i,j+1 8 -m i,j+1 7 m i,j 1 > 0. (A3d)
By interpreting the terms m i,j k as the elements of A, one can notice that:

• Condition (A0) imposes a strictly positive main diagonal entry, similar as for an M -matrix.

• Conditions (A1) impose strictly negative entries only on 4 secondary diagonal entries, unlike the M -matrix that requires negative entries on all extra diagonal terms.

• Condition (A2) imposes a dominance of the main diagonal entries only on 2 extra-diagonal terms instead of a diagonal dominance by row imposed by the M -matrix structure.

• Conditions (A3) impose conditions on cross products involving other secondary diagonals.

Our idea is to design a nonlinear scheme whose linearized version satisfies these conditions; the scheme is based on the flux weighing method ( 5), but these weaker conditions enable us to remove the 5-point stencil constraint of NLTPFA and NLMPFA on the linearized matrix.

Construction of the new scheme

The new nonlinear FV scheme is based on the convex combination of linear fluxes taken in their cell-centered form [START_REF] Bertolazzi | A second-order maximum principle preserving finite volume method for steady convection-diffusion problems[END_REF], which can be put in the following general form on a Cartesian grid

F 1 = λ 1 (f K -f L ) + ν 1 (f K -f M ), F 2 = -λ 2 (f L -f K ) -ν 2 (f L -f N ). (15)
Here, λ 1 , λ 2 , ν 1 , ν 2 > 0, M and N refer to the neighbouring cells of K and L chosen accordingly to the signs of D xy (K) and 2).

D xy (L) (that is, M = M K or M ′ K and N = M L or M ′ L , see Figure
Then, we apply a convex decomposition of the transverse flux branches in F 1 and F 2

F 1 = λ 1 (f K -f L ) + c 1 ν 1 (f K -f M ) + (1 -c 1 )ν 1 (f K -f M ) G 1 , F 2 = -λ 2 (f L -f K ) -c 2 ν 2 (f L -f N ) -(1 -c 2 )ν 2 (f L -f N ) G 2 , (16)
with 0 < c 1 , c 2 < 1. The decomposition parameters c 1 and c 2 are carefully selected (see below) to ensure that the linearized form of the scheme satisfies the monotonicity conditions (A0)-(A3). The expression of F K,σ becomes

F K,σ = (µ 1 λ 1 + µ 2 λ 2 )(f K -f L ) + c 1 µ 1 ν 1 (f K -f M ) + c 2 µ 2 ν 2 (f N -f L ) + [µ 1 G 1 -µ 2 G 2 ]. ( 17
)
We construct µ 1 , µ 2 in a usual way, to either annihilate the bracketed expression in the last equation or express it in a more appropriate way:

     µ 1 = |G 2 | |G 1 | + |G 2 | , µ 2 = |G 1 | |G 1 | + |G 2 | if |G 1 | + |G 2 | ̸ = 0, µ 1 = µ 2 = 1 2 otherwise.
This leads to the following cases:

• if G 1 G 2 ≥ 0 then G 1 |G 2 | = |G 1 |G 2 and thus µ 1 G 1 -µ 2 G 2 = 0.
The expression [START_REF] Cancès | Monotone correction for generic cell-centered finite volume approximations of anisotropic diffusion equations[END_REF] becomes

F K,σ = (µ 1 λ 1 + µ 2 λ 2 )(f K -f L ) + c 1 µ 1 ν 1 (f K -f M ) + c 2 µ 2 ν 2 (f N -f L ). ( 18
) • if G 1 G 2 < 0 then |G 1 |G 2 = -G 1 |G 2 | so µ 1 G 1 -µ 2 G 2 = 2µ 1 G 1 .
The expression [START_REF] Cancès | Monotone correction for generic cell-centered finite volume approximations of anisotropic diffusion equations[END_REF] becomes

F K,σ = (µ 1 ν 1 + µ 2 ν 2 )(f K -f L ) + (2 -c 1 )µ 1 ν 1 (f K -f M ) + c 2 µ 2 ν 2 (f N -f L ). ( 19
)
These expressions show that the flux discretization retains all the cell centered unknowns used by F 1 and F 2 , namely, f K , f L , f M and f N . This leads, after freezing the non-linearities in µ 1 and µ 2 , to a 9-point stencil scheme instead of the 5-point stencil obtained for the NLTPFA and the NLMPFA schemes.

The branches 18)-( 19) prevent the fluxes from having the form [START_REF] Sheng | The finite volume scheme preserving extremum principle for diffusion equations on polygonal meshes[END_REF], so the scheme's matrix will not be, in general, an M-matrix. The key idea is however to tune the coefficients c 1 , c 2 in these branches to ensure that this matrix respect the relaxed monotonicity conditions (A0)-(A3). These coefficients can be determined cell-by-cell, or globally using the most demanding local monotonicity conditions. We choose the latter, and we also associate a unique couple (c 1 , c 2 ) to each discretization direction. Appendix B reports the analytical developments followed to obtain such a suitable couple (c 1 , c 2 ), which only depends on the coefficients λ r , ν r in the linear fluxes [START_REF] Yuan | Monotone finite volume schemes for diffusion equations on polygonal meshes[END_REF] and not on the frozen nonlinear coefficients µ r (and thus not on the Picard iterations), see (B.10).

c 2 µ 2 ν 2 (f N -f L ) and c 2 µ 2 ν 2 (f N -f L ) in (
Due to the cell-centered discretization logic of the new scheme, similar to the NLMPFA, with "relaxed" monotonicity conditions, we refer to the new scheme as Relaxed NLMPFA or R-NLMPFA.

Positivity and extremum preserving properties

We write the scheme in the following form

S(U )U = R, (20)
where U = ((f K ) K∈T , (f σ ) σ∈Eext ) is a vector containing the unknowns of the scheme and the discretized boundary conditions, and R represents the contribution of the volume source S. In [START_REF] Lipnikov | Monotone finite volume schemes for diffusion equations on unstructured triangular and shape-regular polygonal meshes[END_REF], the nonlinearities in S(U ) correspond to those (frozen in the nonlinear resolution procedure) in µ 1 , µ 2 in (17), while the linear term U gathers contributions of f K , f L , f M , f N explicitly visible in ( 18)- [START_REF] Lipnikov | Minimal stencil finite volume scheme with the discrete maximum principle[END_REF].

On the discrete level, the positivity and minimum-maximum preserving principles can be expressed in the following forms if R ≥ 0 and

f σ ≥ 0 for all σ ∈ E ext , then f K ≥ 0 for all K ∈ T , (21)
if R = 0 then min

σ∈Eext f σ ≤ f K ≤ max σ∈Eext f σ for all K ∈ T . ( 22
)
Theorem 1 (Monotonicity of the scheme). The finite volume scheme based on the fluxes defined by (15)-( 19), whose discretization matrix respects (A0)-(A3) when its nonlinear terms are frozen, preserves at each nonlinear iteration both positivity (21) and minimum-maximum principle [START_REF] Sheng | Construction of nonlinear weighted method for finite volume schemes preserving maximum principle[END_REF], achieving consequently full monotonicity.

Proof. The proof of the positivity of the linearized scheme is straightforward due to the monotonicity of its matrix. For the extremum principle we consider R = 0. By the expressions ( 18)-( 19), at each linearization step the matrix S(U s ) satisfies S(U s )1 = 0 (1 is the vector with all components equal to 1). Thus, the vectors V = (max σ∈Eext f σ )1 -U s+1 and W = U s+1 -(min σ∈Eext f σ )1 satisfy v σ ≥ 0 for all σ ∈ E ext and w σ ≥ 0 for all σ ∈ E ext , and S(U s )V = S(U s )W = 0. Applying the positivity principle to V and W then yields min σ∈Eext f σ ≤ f K ≤ max σ∈Eext f σ .

Numerical tests for stationary problems

We test here the numerical performances of R-NLMPFA and we compare them with the results of the NLTPFA and NLMPFA schemes (partially taken from [START_REF] Dahmen | On the modelling of highly anisotropic diffusion for electron radiation belt dynamic codes[END_REF][START_REF] Dahmen | Développement d'un nouveau coeur numérique pour le code de calcul salammbô de modélisation des ceintures de radiation terrestres[END_REF]). First, we start with a numerical test involving a uniform but highly anisotropic diffusion tensor to measure the computational gain achieved by R-NLMPFA. Then we move to three numerical tests related to the positivity, the minimum principle, and the minimum-maximum principle, using an inhomogeneous highly anisotropic diffusion tensor. Finally, we measure the accuracy of R-NLMPFA with a convergence study.For all these numerical tests, we also report the choice of c 1 and c 2 used in the R-NLMPFA scheme construction, following the conditions in Appendix B.

The following notations are used in the section: N u is the number of cell center unknowns, f min and f max , the extremal values of f inside the domain, and N iter is the number of nonlinear Picard iterations to reach convergence fixed by a stopping criteria of ϵ = 10 -6 . R under = N 1 Nu is the ratio of undershoots, with N 1 the number of grid cells where f K < 0 when testing positivity, and N 1 the number of grid cells where f K < min σ∈Eext f σ when testing the minimum principle preserving property. R over = N 2 Nu is the ratio of overshoots, with N 2 the number of grid cells where f K > max σ∈Eext f σ when testing the maximum principle preserving property. To evaluate the rate of convergence, we use the discrete counterpart of the L 2 -norm:

∥f ∥ 2,D = K∈T |K|f 2 K 1/2 (23)
with |K| the area of K. The relative approximation error is calculated based on the previously defined norm:

Err 2 = ∥f -f ref ∥ 2,D ∥f ref ∥ 2,D , (24)
where f ref is the analytical solution chosen by the method of manufactured solution.

Highly anisotropic uniform diffusion tensor test

First, we consider the problem (1) in the unit square domain Ω = (0, 0.5) 2 with S = 0. We fix the following diffusion tensor:

D =   10 7 10 3 10 3 1   . ( 25
)
The latter tensor is uniform, to focus solely on its anisotropy ratio that is ≈ 1, 11 • 10 7 . We fix f (x, y) = sin(πx) sin(πy). The steady state diffusion equation is solved by the NLTPFA, NLMPFA and the R-NLMPFA schemes for three levels of mesh refinements and the Picard algorithm is initialized with unitary function (f 0 = 1). We report in Table 1 the number of iterations N iter needed by each scheme to achieve nonlinear convergence.

The R-NLMPFA scheme requires much fewer nonlinear iterations (up to 6 times fewer) than the NLMPFA, thus preserving the maximum principle at a much lower computational cost. R-NLMPFA shows a similar number of nonlinear iterations as NLTPFA, all the while achieving full monotonicity that NLTPFA lacks. If, instead of the residual

||X s+1 -X s || ∞ ||X s || ∞ , (26)
we use in Algorithm 1 the residual

||A(X s+1 )X s+1 -B s+1 || 2 ||B s || 2 , (27)
then similar computational performances are obtained (with actually fewer iterations for NLTPFA and NLMPFA); see Table 2 (c 1 and c 2 remains the same). The residual ( 27) is global and sensitive to the matrix system convergence that can be attained sometimes before the local convergence of the solution measured by [START_REF] Bourdarie | Electron radiation belt data assimilation with an ensemble kalman filter relying on the salammbô code[END_REF]. Given that this latter residual is therefore actually more demanding, we continue using it in the following tests in this section.

Positivity preserving test

We consider the problem (1) in the unit square domain Ω = (0, 1) 2 . We fix the following diffusion tensor:

D(x, y) = 1 x 2 + y 2   αx 2 + y 2 (α -1)xy (α -1)xy x 2 + αy 2   . ( 28
)
The latter tensor is nonuniform, has 1 and α as eigenvalues, and thus an anisotropy ratio of 1 α ; this ratio is fixed at 10 9 in our tests, as this is the typical magnitude of anisotropy ratios in radiation belts applications [START_REF] Dahmen | On the modelling of highly anisotropic diffusion for electron radiation belt dynamic codes[END_REF][START_REF] Dahmen | Développement d'un nouveau coeur numérique pour le code de calcul salammbô de modélisation des ceintures de radiation terrestres[END_REF]. The source S adopts the following form:

S(x, y) =    1 if (x, y) ∈ [0.25, 0.75] × [0.25, 0.75], 0 otherwise. ( 29
)
We fix homogeneous Dirichlet conditions at x = 0, y = 0, y = 1, and a no-flux boundary at x = 1 (to reflect typical real boundary conditions of radiation belts models). The results are reported in Table 3, in which we can see that all schemes preserve the positivity of the solution; R-NLMPFA however requires fewer Picard iterations N iter , with a reduction of up to 27% compared to NLTPFA and up to 67% compared to NLMPFA.

Minimum principle preserving test

The second numerical is related to the minimum preserving principle. We consider the same diffusion problem as in Section 4.2, except for a full domain Dirichlet boundary condition f = 1 on ∂Ω. The results are reported in Table 4 (c 1 and c 2 remains the same as in 4.2).

NLMPFA and the R-NLMPFA respect the minimum principle (f K ≥ 1 for all K ∈ T ) as they satisfy an extremum principle preserving property in addition to the monotonicity of their discretization matrix. This is not the case of the NLTPFA scheme that is only positivity preserving and presents minimum values under 1. The NLTPFA scheme is the most efficient on this test, but at the price of developing many undershoots. Among the two schemes that preserve the mini- mum principle, the R-NLMPFA outperforms the NLMPFA scheme with up to 36% fewer nonlinear iterations.

Minimum-maximum principle preserving test

For the third test, the domain is reduced to Ω = (0, 0.5) 2 , and we take D as in [START_REF] Potier | Schéma volumes finis monotone pour des opérateurs de diffusion fortement anisotropes sur des maillages de triangles non structurés[END_REF]. We impose f (x, y) = sin(πx) sin(πy) on the boundary at x = 0, y = 0, y = 0.5, and no flux condition at x = 0.5;

we also set S = 0. The results are synthesised in Table 5. For this test, all three schemes preserve the minimum and the maximum principles. The R-NLMPFA confirms its efficiency in terms of nonlinear convergence with lower N iter values than NLTPFA (up to 33% fewer iterations) and NLMPFA (up to 55% fewer iterations).

Convergence study

We consider the problem (1) in the domain Ω = (0, 0.5) 2 , with the diffusion tensor (28) and α = 10 -6 . We fix f ref (x, y) = sin(πx) sin(πy) and S is computed from this exact solution. The results of the convergence study of the three schemes are reported in Figure 4 and Table 6. The R-NLMPFA shows a second order of convergence, with error magnitudes similar to the NLMPFA scheme. This might be explained by their similar cell-centered approach. Again the R-NLMPFA is more competitive with up to 35% fewer iterations than the NLMPFA scheme and up to 24% fewer iterations than the NLTPFA scheme.

Application to radiation belts

Earth's electron radiation belts are structures evolving around Earth and crossing the majority of satellite orbits (see Figure 5) [START_REF] Horne | Space weather impacts on satellites and forecasting the earth's electron radiation belts with spacecast[END_REF]. Filled with highly energetic electrons trapped by Earth's magnetic field, they can present a rapidly varying and asymmetric dynamic correlated to solar activity.

The usual theoretical development adopted to describe radiation belt dynamics [START_REF] Schulz | Particle diffusion in the radiation belts[END_REF][START_REF] Roederer | Dynamics of magnetically trapped particles[END_REF] leads to a 3 dimension diffusion equation, governing the evolution of a phase space density (PSD) function

f (M eV -1 .s -3
) and expressed in the space (y, E, L * ) with

• y = sin α eq the sine of the particle's equatorial pitch angle, the angle between its velocity vector and the magnetic field at the equator (see Figure 5).

• E the particle's energy expressed in M eV .

• L * the radial distance of the magnetic field line (dipole approximation) traveled by the trapped particle at the equator (see Figure 5) [START_REF] Roederer | Dynamics of magnetically trapped particles[END_REF]. Due to the multi-scale, localized and highly asymmetric physical processes involved in radiation belt dynamics, the local diffusion frame expressed in (y, E) is spatially inhomogeneous and presents high levels of anisotropy as shown in Figure 6 for D taken at L * = 4.49. The 2D tensor anisotropy ratio is also represented in Figure 7.

∂f ∂t = 1 G ∇ • (GD∇f ) (30)
with

D =   D yy D yE D yE D EE   (31)
and G the Jacobian of the transformation from the canonical space to (y, E). Based on this diffusion model, the radiation belt study community implemented several physical codes like the pioneering ONERA's Salammbô code [START_REF] Varotsou | Three-dimensional test simulations of the outer radiation belt electron dynamics including electron-chorus resonant interactions[END_REF][START_REF] Bourdarie | Electron radiation belt data assimilation with an ensemble kalman filter relying on the salammbô code[END_REF][START_REF] Beutier | A three-dimensional analysis of the electron radiation belt by the salammbô code[END_REF], the British Antarctic Survey (BAS)

Radiation Belt model code [START_REF] Glauert | Three-dimensional electron radiation belt simulations using the bas radiation belt model with new diffusion models for chorus, plasmaspheric hiss, and lightning-generated whistlers[END_REF], the Storm-Time Evolution of Electron Radiation Belt (STEERB) code [START_REF] Su | Steerb: A three-dimensional code for storm-time evolution of electron radiation belt[END_REF] and the Versatile Electron Radiation Belt (VERB) code [START_REF] Subbotin | Three-dimensional modeling of the radiation belts using the versatile electron radiation belt (verb) code[END_REF]. Cumulative improvements in the modelling of physical processes occurring in the inner magnetosphere led to more and more accurate depictions of the radiation belt dynamics by these codes, especially during the main phases of geomagnetic events. However, their finite difference based numerical resolution started to show its limit to tackle all the induced numerical challenges [START_REF] Dahmen | On the modelling of highly anisotropic diffusion for electron radiation belt dynamic codes[END_REF][START_REF] Dahmen | Développement d'un nouveau coeur numérique pour le code de calcul salammbô de modélisation des ceintures de radiation terrestres[END_REF].

Recently, a transition to a finite volume based resolution was initiated in the Salammbô code.

This was motivated by the suitability of the FV method for diffusion equations and to take advantage of the conservativity property. In order to cope with the high levels of anisotropy and take into account the tricky cross diffusion terms D yE without numerical artefacts, the use of nonlinear finite volume schemes seemed inevitable. Thus, NLPTFA and NLMPFA schemes were adapted and tested over the tricky pitch angle and energy diffusion [START_REF] Berman | Nonnegative matrices in the mathematical sciences[END_REF], with the results reported in [START_REF] Dahmen | On the modelling of highly anisotropic diffusion for electron radiation belt dynamic codes[END_REF].

In this section, we test the R-NLMPFA in the radiation belt context and we focus on its computational cost, which is compared to NLTPFA and NLMPFA. We consider the transient diffusion problem (30) evaluated over T simu = 90000s. We adopt the 2D diffusion tensor D evaluated at the plan L * = 4.49 and showed in Figure 6. Note that the diffusion coefficients were derived from estimations of resonant particle-particle interactions wave-particle (using the WAPI code [START_REF] Sicard-Piet | Wapi: A new model for the wave particle interaction[END_REF]).

We start with the initial state f init = 10 30 and we impose the boundary conditions presented in Table 7.

The Picard algorithm is initialized with the output of a linear FV scheme, cleaned from overshoots and undershoots, and we adopt the residual [START_REF] Glauert | Three-dimensional electron radiation belt simulations using the bas radiation belt model with new diffusion models for chorus, plasmaspheric hiss, and lightning-generated whistlers[END_REF], more adapted to the high dispersion of the solution. We also fix a maximum number of nonlinear iterations N iter,max = 300, and we consider that the Picard algorithm did not converge above that number. We report in Table 8 the extremum values of f , the average number of Picard iterations per time iteration N iter,avg , the percentage of time iterations where the Picard algorithm did not converge (N iter > 300).

Table 7: Boundary conditions adopted for the transient simulation. For all the time steps ∆t we tested, the NLPTFA scheme showed the least N iter,avg values (at maximum 2 Picard iterations per time iteration), followed very closely by the R-NLMPFA (at maximum 4 Picard iterations per time iteration). Hence, the transition from a only positivity preserving scheme to a full monotonicity preserving scheme in our application does not sensibly increase the computational cost. Focusing on the two schemes that also preserve extremum values, the R-NLMPFA prevails with much lower N iter,avg values than the NLMPFA schemes. Most notably, we did not observe with R-NLMPFA any oscillating profiles in the residuals, contrary to those observed with the NLMPFA scheme (see Figure 8). This annoying behavior, frequently reported with the NLMPFA scheme, is mainly reponsible for the non-convergence of the Picard algorithm at some time iterations. The R-NLMPFA very likely avoids their occurrence due to its relaxed linearized stencil.

Boundary location Boundary condition type Explanation

Conclusion

We presented in this paper the construction of a new nonlinear finite volume scheme, R-NLMPFA, adapted to highly anisotropic diffusion equations. This scheme shares the convex combination structure adopted by typical nonlinear FV schemes and a cell-centered formulation. However, it stands out with its monotonicity preserving property obtained by a set of conditions less restrictive than those needed for an M-matrix construction. This construction leads to a full monotonicity preserving scheme, i.e., a positivity and minimum-maximum principle preserving scheme, in addition to a larger linearized stencil over Cartesian grids. Adimensionnal numerical tests involving diffusion tensors with very high anisotropy ratios showed an important reduction of the computational cost with the R-NLMPFA, especially when compared to the analogous extremum preserving scheme NLMPFA. These results were also confirmed on the numerically challenging radiation belts application. The next step consists in generalizing the construction of the R-NLMPFA to three dimensions and/or quadrilateral meshes, which seems achievable since the monotonicity conditions at the core of the construction are still valid for these types of mesh. We unify the previous expression for more clarity by introducing the parameter θ σp such that

θ σp =      c σp 1 if G σp 1 G σp 2 ≥ 0 2 -c σp 1 if G σp 1 G σp 2 < 0
Then (B.1) becomes

F K,σ i = (µ σp 1 λ σp 1 + µ σp 2 λ σp 2 )(f K -f L σp ) + θ σp µ σp 1 ν σp 1 (f K -f M σp ) + c σp 2 µ σp 2 ν σp 2 (f N σp -f L σp ). (B.2)
We define a couple (c 1 , c 2 ) on each discretization direction, i.e., a unique couple for each σ p is applied on each mesh cell. Thus, we have to determine four pairs of weighing coefficients. Imposing the conservativity property (4) links the coefficients on both sides of each edge: This reduces the number of weights to be determined to 4, namely (c σ 1 1 , c σ 1 2 , c σ 2 1 , c σ 2 2 ).

We then apply the condition (A0), that is m i,j 1 > 0. From (B.2) we find This shows, since all these coefficients are strictly positive by construction, that (A0) is always valid.

m i,j 1 = µ σ 1 1 λ σ 1 1 + µ σ 1 2 λ σ 1 2 + θ σ 1 µ σ 1 1 ν σ 1 1 + µ σ 2 1 λ σ 2 1 + µ σ 2 2 λ σ 2 2 + θ σ 2 µ σ 2 1 ν σ 2 1 + µ σ 3 1 λ σ 3 1 + µ σ 3 2 λ σ 3 2 + θ σ 3 µ σ 3 1 ν σ 3 1 + µ σ 4 1 λ σ 4 1 + µ σ
We now move to conditions (A1): m i,j 2 < 0 , m i,j 4 < 0 , m i,j 6 < 0 , m i,j 8 < 0.

The expressions of these coefficients depend on the location of M σp in the local cell numbering around the cell K (see Figure 3), which in turn depends on the sign of D xy (K). To take this into account and later on for (A3)) we define the following parameters

ψ σp 1 =      1 if D xy (K) ≥ 0 0 if D xy (K) < 0 (B.5) m i,j 1 > 0.
The expressions of m 3 , m 5 , m 7 , m 9 depend on the location of M σp and N σp which is defined by the sign of D xy on K and L σp . To take that into account, we define the following parameter:

ψ σp 2 =      1 if D xy (L σp ) ≥ 0, 0 if D xy (L σp ) < 0, (B.8)
and from (B.2) we deduce the expressions of m i,j 3 , m i,j 5 , m i,j 7 , m i,j 9 :

Figure 1 :

 1 Figure 1: Cartesian mesh notations

ThusD

  (K)n K,σ = D(K)e x = D xx (K)e x + D xy (K)e y , D(L)n L,σ = -D(L)e x = -D xx (L)e x -D xy (L)e y .

Figure 2 :

 2 Figure 2: Construction and stencil unknowns of the FK,σ flux discretization on a Cartesian grid, for to the NLTPFA scheme (in red) and for to the NLMPFA scheme (in blue: stencil of F1, in green: stencil of F2 -solid line refers to the case Dxy ≥ 0, dashed line refers to the case Dxy < 0). Note that, in NLMPFA, F2 is eliminated from the final expression of FK,σ.

Figure 3 :

 3 Figure 3: Local cell numbering for a nine-point scheme on a Cartesian grid, adapted from [25].

Figure 4 :

 4 Figure 4: Rates of convergence of NLTPFA, NLMPFA and R-NLMPFA

Figure 5 :

 5 Figure 5: Electron radiation belts representation with the L * parameter and the equatorial pitch angle αeq, in addition to the Van Allen probes orbit, adapted from [34].

Figure 6 :

 6 Figure 6: 2D mapping of the diffusion tensor terms Dyy, DEE and DyE, taken at L * = 4.49.

Figure 7 :

 7 Figure 7: Ratio of anisotropy of the diffusion tensor taken at L * = 4.49.

E

  min = 5.664 • 10 -4 M eV f = 10 30 M eV -1 s -3Low energy seeding E max = 12.561 M eV f = 0 Limit of magnetic trapping at very high energy

Figure 8 :

 8 Figure 8: Residuals at each Picard iteration step, at n × ∆t = 70000s for ∆t = 10000s

Figure B. 9 :

 9 Figure B.9: edge notations on a cell K.

4 2 λ σ 4 2 + θ σ 4 µ σ 4 1 ν σ 4 1 .(B. 4 ) 4 r=1 (µ σr 1 λ σr 1 + µ σr 2 λ σr 2 +

 21442 which can be written in the following convinient form m i,j 1 = θ σr µ σr 1 ν σr 1 ) .

Table 1 :

 1 Anisotropic uniform diffusion tensor test: comparison of NLTPFA, NLMPFA and R-NLMPFA on uniform Cartesian grids.

		N u	20 × 20	40 × 40	80 × 80
	Scheme
	NLTPFA	N iter	4	5	8
	NLMPFA	N iter	14	25	45
	R-NLMPFA N iter	4	5	7
		c 1	8.33 • 10 -6 8.33 • 10 -6 8.33 • 10 -6
		c 2	1.67 • 10 -5 1.67 • 10 -5 1.67 • 10 -5

Table 2 :

 2 Anisotropic uniform diffusion tensor test using the residual (27): comparison of NLTPFA, NLMPFA and R-NLMPFA on uniform Cartesian grids.

		N u 20 × 20 40 × 40 80 × 80
	Scheme
	NLTPFA	N iter	2	2	3
	NLMPFA	N iter	4	6	9
	R-NLMPFA N iter	2	2	2

Table 3 :

 3 Postivity preserving numerical test: comparison of the NLTPFA, NLMPFA and R-NLMPFA on uniform Cartesian grids.

		N u	20 × 20	40 × 40	80 × 80
	Scheme
	NLTPFA	f N iter	81	141	222
	NLMPFA	f min 2.82 • 10 -7 9.15 • 10 -10 2.71 • 10 -13
		N iter	153	310	433
	R-NLMPFA f min 1.57 • 10 -7	3.6 • 10 -10 8.85 • 10 -14
		N iter	68	102	193
		c 1	2.88 • 10 -4 1.64 • 10-4 8.76 • 10 -5
		c 2	5.76 • 10 -4	3.28 • 10 -4 1.75 • 10-4

min 6.44 • 10 -10 7.55 • 10 -15 2.52 • 10 -23

Table 4 :

 4 Minimum preserving test: comparison of the NLTPFA, NLMPFA and R-NLMPFA on uniform Cartesian grids.

		N u	20 × 20 40 × 40 80 × 80
	Scheme
	NLTPFA	f min	0.9983	0.9984	0.9987
		N iter	62	83	66
		R under 14.25% 14.31% 13.96%
	NLMPFA	f min	1.00	1.00	1.00
		N iter	66	107	201
	R-NLMPFA	f min	1.00	1.00	1.00
		N iter	58	93	128

Table 5 :

 5 minimum-maximum preserving test: comparison of the NLTPFA, NLMPFA and R-NLMPFA on uniform Cartesian grids.

		N u	20 × 20	40 × 40	80 × 80
	Scheme
	NLTPFA	f min 2.22 • 10 -15 4.67 • 10 -27 3.60 • 10 -44
		f max	0.9955	0.9988	0.9997
		N iter	95	181	338
	NLMPFA	f min 2.10 • 10 -9 1.63 • 10 -13 3.73 • 10 -18
		f max	0.9950	0.9988	0.9995
		N iter	126	250	606
	R-NLMPFA f min 7.30 • 10 -10 4.94 • 10 -14 6.80 • 10 -20
		f max	0.9918	0.9982	0.9995
		N iter	63	127	271
		c 1	2.88 • 10 -4	1.64 • 10 -4	8.76 • 10 -5
		c 2	5.76 • 10 -4	3.28 • 10 -4	1.75 • 10 -4

Table 6 :

 6 Convergence study: comparison of the NLTPFA, NLMPFA and R-NLMPFA on uniform Cartesian grids, Picard iterations initialised with f 0 = 1.

		N u	20 × 20	40 × 40	80 × 80
	Scheme
	NLTPFA	N iter	63	112	185
	NLMPFA	N iter	79	134	216
	R-NLMPFA N iter	66	101	140
		c 1	2.88 • 10 -4 1.64 • 10 -4 8.78 • 10 -5
		c 2	5.75 • 10 -4 3.28 • 10 -4 1.76 • 10 -4

Table 8 :

 8 Transient simulation results with the NLTPFA,NLMPFA and R-NLMPFA schemes for different ∆t, at Tsimu = 90000s.

	Scheme	∆t	100s	1000s	10000s
	NLTPFA	f min	3.96 • 10 26	4.05 • 10 26	4.97 • 10 26
		f max	1.000001191 • 10 30 1.000000448 • 10 30 9.99 • 10 29
		N iter,avg	1	1	2
	NLMPFA	f min	3.95 • 10 26	4.06 • 10 26	4.96 • 10 26
		f max	9.99 • 10 26	9.99 • 10 26	9.99 • 10 29
		N iter,avg	1.23	7.95	242.66
		% no-conv	0	0	66.66%
	R-NLMPFA	f min	4.12 • 10 26	4.34 • 10 26	5.18 • 10 26
		f max	9.99 • 10 29	9.99 • 10 29	9.99 • 10 29
		N iter,avg	1	1.01	3.11

Acknowledgments

This research was supported by CNES -The French Space Agency, ONERA -The French Aerospace Lab, and the University of Monash.

Code availability

The convergence study python code was made available and can be accessed in the following link : https://github.com/ndahmen/R-NLMPFA

Appendix A. Construction of NLTPFA and NLMPFA convex combination weights

The NLTPFA scheme construction is based on the forms (9) of F 1 and F 2 . Injecting these expressions in [START_REF] Dahmen | Développement d'un nouveau coeur numérique pour le code de calcul salammbô de modélisation des ceintures de radiation terrestres[END_REF] gives

with a 1 , a 2 > 0 containing the edge unknowns, responsible for the loss of the M -matrix structure.

The expression µ 1 a 1 -µ 2 a 2 is annihilated by imposing µ 1 a 1 -µ 2 a 2 = 0 and µ 1 + µ 2 = 1, which yields

The NLMPFA scheme construction, on the other hand, is based on the cell-centered form of F 1 and F 2 as reported in the equation [START_REF] Bertolazzi | A second-order maximum principle preserving finite volume method for steady convection-diffusion problems[END_REF], which can be framed in the following general form (with the same notation as in subsection 3.2 -so M and N are chosen accordingly to the signs of D xy (K) and D xy (L))

Injecting these expressions in (5) yields

The NLMPFA scheme aims to recast the latter flux expression in the multi-point form [START_REF] Sheng | The finite volume scheme preserving extremum principle for diffusion equations on polygonal meshes[END_REF], by discarding the transverse flux branch G ′ 2 responsible for the loss of the M -matrix structure. Thus, µ 1 , µ 2 are chosen to either annihilate the bracketed expression in the last equation, or express it in a form that is appropriate to the LMP structure:

In other words:

and from (B.2) we deduce the expressions of m i,j 2 , m i,j 4 , m i,j 6 , m i,j 8 :

Similarly as for (A0), the conditions (A1) are then valid for any choice of c 1 and c 2 , since 0 < µ 1 , µ 2 , c 1 , c 2 < 1 and ν 1 , ν 2 > 0.

Let us now consider (A2), that is, m i,j 1 + m i,j 2 + m i,j 6 > 0 which translates using (B.2) into

This leads to conditions on c σ 1 2 , c σ 2 2 , c σ 3 2 and c σ 4 2 :

Combined with (B.7), this leads to the five following inequalities to be ensured by

From here, one can assume a relation between the coefficients c σ 1 1 , c σ 1 2 , c σ 2 1 and c σ 2 2 to norrow 25 their choices from B.10. We can fix the following relations

and

Adding these new relations in B.10 gives the following conditions on c 1

2 ) min i,j+1 (λ σ 4 1 , λ σ 4 2 ) max i,j+1 (ν σ 1 2 , ν σ 4 2 )A ′ All that remains is to select the smallest c 1 over the whole mesh grid for the numerical solving.