Stéphane Abide
email: stephane.abide@univ-perp.fr

Belkacem Zeghmati

Multigrid defect correction and fourth-order compact scheme for Poisson's equation

This paper presents an analysis of a multigrid defect correction to solve a fourth-order compact scheme discretization of the Poisson's equation. We focus on the formulation, which arises in the velocity/pressure decoupling methods encountered in computational fluid dynamics. Especially, the Poisson's equation results of the divergence/gradient formulation and Neumann boundary conditions are prescribed. The convergence rate of a multigrid defect correction is investigated by means of an eigenvalues analysis of the iteration matrix. The stability and the mesh-independency are demonstrated. An improvement of the convergence rate is suggested by introducing the damped Jacobi and Incomplete Lower Upper smoothers. Based on an eigenvalues analysis, the optimal damping parameter is proposed for each smoother. Numerical experiments confirm the findings of this analysis for periodic domain and uniform meshes which are the working assumptions. Further numerical investigations allow us to extend the results of the eigenvalues analysis to Neumann boundary conditions and non-uniform meshes. The Hodge-Helmholtz decomposition of a vector field is carried out to illustrate the computational efficiency, especially by making comparisons with a second-order discretization of the Poisson's equation solved with a state of art of algebraic multigrid method.

Introduction

Since the early days of computing, many research works have been published on the development of numerical solutions of Poisson's equation. This elliptic partial differential equation has many applications in various areas of applied physics such as acoustic, electromagnetism, heat and mass transfer, etc. For instance, the Poisson solver is an essential part for solving incompressible fluid flows. In this case, the computational cost saving still remains of primary importance, that has led to an important amount of papers in this research area. The combination of a discrete spatial approximation and a linear solver is the key point of such methods of solutions for elliptic partial differential equations. For instance, the progress in multigrid methods has allowed us to achieve a good scalability for finite-element, finite-difference or finite-volume discretizations. In order to reduce the computational cost, yet another popular way is to discretize the equations with high-order schemes like the spectral methods. High-order methods achieve a better accuracy than lower-order discretizations with identical grid size. The compact finite difference schemes belong to this category of high-order schemes. They have the particular feature to combine the accuracy with linear computational complexity. For this reason, the compact finite difference schemes have a good scaling potential. However, the combination of such discretizations with an efficient linear solver for elliptic boundary value problems still remains an open issue. In this way, the main objective of the present work is to propose an iterative method of solving the Poisson equation discretized based on a fourth-order compact finite difference schemes. A special care is paid to the singular Poisson equation arising from incompressible fluid flow solvers. Here the discrete formulation results from the combination of the discrete divergence with the discrete gradient. The Neumann boundary conditions are also considered.

The compact finite difference schemes achieve a high-order truncation error while keeping the smallest computational stencil. One of pioneer way to obtain those schemes consists in approximating the partial derivative equation with a secondorder finite difference and in replacing the higher derivatives in the truncation error by the lower-order derivatives of equation itself. The well-known Mehrstellen discretization of the Poisson equation, proposed by Collatz [START_REF] Collatz | The Numerical Treatment of Differential Equations[END_REF], is among the first published papers relating to this procedure. A large amount of works follows this trend [START_REF] Spotz | High-Order compact finite difference schemes for computational mechanics[END_REF][START_REF] Spotz | A high-order compact formulation for the 3D Poisson equation[END_REF][START_REF] Spotz | Formulation and experiments with high-order compact schemes for nonuniform grids[END_REF][START_REF] Gupta | High accuracy multigrid solution of the 3D convection diffusion equation[END_REF][START_REF] Ge | High accuracy iterative solution of convection diffusion equation with boundary layers on nonuniform grids[END_REF][START_REF] Ge | A transformation-free HOC scheme and multigrid method for solving the 3D Poisson equation on nonuniform grids[END_REF][START_REF] Wang | Sixth order compact scheme combined with multigrid method and extrapolation technique for 2D Poisson equation[END_REF][START_REF] Gupta | Comparison of second-and fourth-Order discretizations for multigrid Poisson solvers[END_REF][START_REF] Pardhanani | A stable multigrid strategy for convection-diffusion using high order compact discretization[END_REF]. For instance, Spotz and co-workers proposed high-order compact schemes (HOCS) formulation for solving numerically the two-dimensional convection/diffusion equation [START_REF] Spotz | High-Order compact finite difference schemes for computational mechanics[END_REF], the three-dimensional Poisson equation [START_REF] Spotz | A high-order compact formulation for the 3D Poisson equation[END_REF], with extension to non-uniform grids [START_REF] Spotz | Formulation and experiments with high-order compact schemes for nonuniform grids[END_REF]. The contributions to HOCS of Zhang and co-workers concern also elliptic boundary value problems [START_REF] Gupta | High accuracy multigrid solution of the 3D convection diffusion equation[END_REF][START_REF] Ge | High accuracy iterative solution of convection diffusion equation with boundary layers on nonuniform grids[END_REF][START_REF] Ge | A transformation-free HOC scheme and multigrid method for solving the 3D Poisson equation on nonuniform grids[END_REF][START_REF] Wang | Sixth order compact scheme combined with multigrid method and extrapolation technique for 2D Poisson equation[END_REF]. A threedimensional HOCS of the convection/diffusion equation was proposed [START_REF] Gupta | High accuracy multigrid solution of the 3D convection diffusion equation[END_REF]. Furthermore, an analytical mesh transformation was developed to deal with boundary layer type problems [START_REF] Ge | High accuracy iterative solution of convection diffusion equation with boundary layers on nonuniform grids[END_REF], and a transformation-free mesh of three-dimensional Poisson equation [START_REF] Ge | A transformation-free HOC scheme and multigrid method for solving the 3D Poisson equation on nonuniform grids[END_REF] was recently proposed. A six-order accuracy was achieved with the compact scheme presented by Wang and Zhang [START_REF] Wang | Sixth order compact scheme combined with multigrid method and extrapolation technique for 2D Poisson equation[END_REF]. These works illustrated a first set of methods to design accurate multi-dimensional schemes for elliptic boundary value problems. The common feature is that they lead to a sparse linear system, which has to be solved efficiently. Iterative methods for the solution of linear systems of equations are usually used. The Krylov projection type methods [START_REF] Spotz | High-Order compact finite difference schemes for computational mechanics[END_REF] or multigrid methods [START_REF] Pardhanani | A stable multigrid strategy for convection-diffusion using high order compact discretization[END_REF][START_REF] Gupta | Comparison of second-and fourth-Order discretizations for multigrid Poisson solvers[END_REF] have been successfully developed. The recent geometric multigrid procedure designed by Wang and Zhang [START_REF] Wang | Sixth order compact scheme combined with multigrid method and extrapolation technique for 2D Poisson equation[END_REF] has a convergence factor independent of the mesh size. This illustrates how compact finite difference schemes could combine accuracy with a linear computational complexity. In the framework of computational fluid dynamics, these schemes have been applied to the ψω formulation of the incompressible Navier-Stokes equations [START_REF] Zhang | Numerical simulation of 2D square driven cavity using fourth-order compact finite difference schemes[END_REF][START_REF] Spotz | High-Order compact scheme for the steady stream-function vorticity equations[END_REF]. It should be mentioned that in the computational fluid dynamics the Poisson's equation generally arises from the combination of gradient and divergence operators, especially with the class of projection methods [START_REF] Temam | Une méthode d'approximation des solutions des équations Navier Stokes[END_REF][START_REF] Chorin | Numerical solution of the Navier Stokes equations[END_REF].

Yet another popular way to derive the compact schemes is based on Taylor series expansions. Here linear combinations on the discrete derivative (interpolation) and the function involved are determined by fulfilling requirements on accuracy and/or resolution. The core of the method was detailed in the seminal paper of Lele [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF], and still remains a starting point for many other high-order schemes [START_REF] Fu | Modified weighted compact scheme with global weights for shock capturing[END_REF][START_REF] Chu | A three-point sixth-order staggered combined compact difference scheme[END_REF][START_REF] Zhou | A new family of high-order compact upwind difference schemes with good spectral resolution[END_REF]. Their derivations are based on a one-dimensional computational stencil. The multidimensional extension can be simply done by means of a tensor product. However, in this case, the resulting linear system of equations is not sparse which introduces several difficulties to get a solution efficiently, especially for Poisson equation arising in the incompressible fluid flow solvers [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF][START_REF] Pozrikidis | A note on the regularization of the discrete PoissonNeumann problem[END_REF]. Nevertheless, several direct and iterative methods of solutions have been used [START_REF] Vedy | A high-order finite difference method for incompressible fluid turbulence simulations[END_REF][START_REF] Abide | A 2D compact fourth-order projection decomposition method[END_REF][START_REF] Schiestel | A Hermitian-Fourier Numerical method for solving the incompressible Navier-Stokes equations[END_REF][START_REF] Knikker | Study of a staggered fourth-order compact scheme for unsteady incompressible viscous flows[END_REF]. For instance, the diagonalization method, which was first used in combination with spectral methods [START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF], still remains retained as a method of solutions for linear systems arising from compact schemes discretization [START_REF] Abide | A 2D compact fourth-order projection decomposition method[END_REF][START_REF] Tyliszczak | A high-order compact difference algorithm for half-staggered grids for laminar and turbulent incompressible flows[END_REF]. This direct method has an unfavorable computational complexity, and seems not to be suitable for large scale computing. An another way based on FFT leads to efficient numerical solvers, but it introduces some limitations on the number of stretching mesh directions [START_REF] Zhuang | A high-order fast direct solver for singular Poisson equations[END_REF][START_REF] Laizet | A numerical strategy to combine high-order schemes, complex geometry and parallel computing for high resolution {DNS} of fractal generated turbulence[END_REF]. The early work of Schiestel and Viazzo [START_REF] Schiestel | A Hermitian-Fourier Numerical method for solving the incompressible Navier-Stokes equations[END_REF] introduces an iterative procedure between the momentum and the pressure correction equations. Their method is based on a variant of the defect correction method [START_REF]Defect Corrections and Multigrid Iterations[END_REF]. Knikker [START_REF] Knikker | Study of a staggered fourth-order compact scheme for unsteady incompressible viscous flows[END_REF] suggested a defect correction procedure as method of solutions for the Poisson's equation discretized with a fourth-order compact finite difference scheme. The principle consists in solving a lower order accurate discretization of the Poisson's equation and in correcting the residual, which is computed using a higher-order discretization [START_REF]Defect Corrections and Multigrid Iterations[END_REF]. Introducing an under-relaxation coefficient, this iterative method has demonstrated the efficiency for the expected application, but no details concerning the convergence rate were specified. Such an analysis was given by Brüger et al. [START_REF] Brüger | High order accurate solution of the incompressible NavierStokes equations[END_REF] for an iterative Krylov method preconditioned with an Incomplete Lower Upper (ILU) factorization. Since the matrix was not available explicitly, the ILU factorization concerned a second-order finite differences approximation of the Poisson equation. The analysis of the convergence factor showed a condition number which behaves as o(h -1 ), h being the grid space.

The above review presents results of several works relating the efficiency of the Poisson solvers based on high-order compact schemes combined with multigrid. In these works, the Poisson equation did not arise from the combination of the discrete divergence and the gradient operators. For this reason, these methods are not applicable to the Poisson equation occurring in the incompressible fluid flow solvers. In this context, only few works considered such discretizations with an iterative linear solver, also there is a little amount of information on the convergence factor. In the present work, we are going to give some insights into defect correction applied to the Poisson equation with a particular focus on the convergence factor analysis. Yet another aim of this work is to show how the convergence factor can be improved by using a smoothing step. Section 1 deals with the spatial discretization and the analysis of the convergence factor of the defect correction method. Then several numerical experiments highlight the accuracy and the influence of the smoother operator on the convergence factor. The computational cost issues are also addressed with the Hodge-Helmholtz decomposition of a vector field, which mimics the pressure equation in incompressible fluid flow solvers.

Numerical methods

Fourth-order compact finite difference scheme

The two-dimensional Poisson equation with the Neumann boundary conditions is considered, viz.,

-  ∂ x ∂ x + ∂ y ∂ y  φ = f (x, y), (x, y) ∈ Ω, ∂ n φ = g(x, y), (x, y) ∈ ∂Ω, ( 1
)
where Ω is a bounded domain and ∂Ω is the boundary of Ω. The solution φ, the forcing term f and the boundary values function g are assumed to be sufficiently smooth. This equation is defined on a square domain Ω = [0, 1] × [0, 1] and a uniform mesh size h x = h y = h = n -1 is considered. The function φ in Eq. ( 1) is defined on the nodes (x i-1/2 , y j-1/2 ) with x i-1/2 = (i -1/2)h and y j-1/2 = (j -1/2)h, where 1 ≤ i, j ≤ n. This corresponds to the classical layout of pressure like variables on a staggered grid. According to the projection methods [START_REF] Temam | Une méthode d'approximation des solutions des équations Navier Stokes[END_REF][START_REF] Chorin | Numerical solution of the Navier Stokes equations[END_REF], the Poisson equation is built as the composition of the divergence and gradient operators. The operators are defined on the staggered grid. This implies that the gradient and the divergence are evaluated at the cell faces and at the cell centers, respectively. In particular, the components of the gradient operator are evaluated at the nodes (x i , y j-1/2 ) and (x i-1/2 , y j ), whereas the discrete divergence is computed at the cell centers (x i-1/2 , y j-1/2 ). The notation δ cf x(y) refers to the cartesian components of the gradient operator and δ fc

x(y) refers to the derivatives involved in the divergence operator. Thus the discrete Poisson problem, expressed with tensorial notations, reads

-  δ fc x δ cf x ⊗ I + I ⊗ δ fc y δ cf y  φ = f (2)
with the following Neumann boundary conditions

 δ cf x φ   0,j-1/2 = g 0,j-1/2 δ cf x φ   n x ,j-1/2 = g n x ,j-1/2 ,  δ cf y φ   i-1/2,0 = g i-1/2,0 δ cf y φ   i-1/2,n y = g i-1/2,n y . (3)
Each discrete derivative δ cf x(y) and δ fc x(y) is discretized with the fourth-order compact finite-difference. The cell-to-face derivative for the inner nodes is given by 1 24

φ ′ i-1 + 11 12 φ ′ i + 1 24 φ ′ i+1 = 1 h (φ i+1/2 -φ i-1/2 ) + o  h 4  , (4)
and the face-to-cell corresponds to the shift version of Eq. ( 4)

1 24 φ ′ i-3/2 + 11 12 φ ′ i-1/2 + 1 24 φ ′ i+1/2 = 1 h (φ i -φ i-1 ) + o  h 4  . (5)
The boundary nodes discretization of the face-to-cell derivative is

φ ′ -1/2 + 23φ ′ 1/2 = 1 h (26φ 0 -25φ 1 + φ 2 ) + o  h 3  . (6)
The cell-to-face derivative boundary relation is given by

φ ′ 0 = 1 24h (-23φ -1/2 + 21φ 1/2 + 3φ 3/2 -φ 5/2 ) + o  h 3  . (7)
Using Eqs. ( 4)-( 7), each discrete derivative may be rewritten in the matrix form, viz.,

M cf x(y) φ ′ = D cf x(y) φ and M fc x(y) φ ′ = D fc x(y) φ. (8)
With this notation, M cf x(y) and M fc x(y) are tridiagonal non singular matrices. The discrete fourth-order formulation of the Poisson equation is then

L ′ φ = -   M fc x  -1 D fc x  M cf x  -1 D cf x ⊗ I + I ⊗  M fc y  -1 D fc y  M cf y  -1 D cf y  φ = f . (9)
The finite difference scheme for the boundary conditions Eq. ( 3) is

g 0,j-1/2 = 1 24h (-23φ -1/2,j-1/2 + 21φ 1/2,j-1/2 + 3φ 3/2,j-1/2 -φ 5/2,j-1/2 ), g n x ,j-1/2 = - 1 24h (-23φ n x +1/2,j-1/2 + 21φ n x -1/2,j-1/2 + 3φ n x -3/2,j-1/2 -φ n x -5/2,j-1/2 ), g i-1/2,0 = 1 24h (-23φ i-1/2,-1/2 + 21φ i-1/2,1/2 + 3φ i-1/2,3/2 -φ i-1/2,5/2 ), g i-1/2,n y = - 1 24h (-23φ i-1/2,n y +1/2 + 21φ i-1/2,n y -1/2 + 3φ i-1/2,n y -3/2 -φ i-1/2,n y -5/2 ). (10)
Non-uniform meshes are considered via the introduction of an analytic mesh mapping. If x(X ) is a non singular mesh mapping, with X referring to the uniform mesh coordinates, a derivative is evaluated as δ x = x -1 X δ X . For further details see Ref. [START_REF] Knikker | Study of a staggered fourth-order compact scheme for unsteady incompressible viscous flows[END_REF]. It should be noted that the discrete operator L ′ is a dense matrix due to the inverses of M cf x(y) and M fc x(y) . This makes difficult to develop efficient methods to solve Eq. ( 9), and specific strategies have to be drawn. The direct method of diagonalization may be used to solve Eqs. ( 9) and [START_REF] Pardhanani | A stable multigrid strategy for convection-diffusion using high order compact discretization[END_REF] [START_REF] Vedy | A high-order finite difference method for incompressible fluid turbulence simulations[END_REF], or by its multidomain extension [START_REF] Abide | A 2D compact fourth-order projection decomposition method[END_REF]. Few iterative methods of solutions are available [START_REF] Brüger | High order accurate solution of the incompressible NavierStokes equations[END_REF][START_REF] Knikker | Study of a staggered fourth-order compact scheme for unsteady incompressible viscous flows[END_REF]. This way is investigated in the next section by considering an analysis of the multigrid defect correction method applied to Eqs. ( 9) and (10).

Convergence analysis of the defect correction

The multigrid defect correction solves the problem L ′ φ = c by means of an auxiliary problem Lφ = c for which an efficient solution is available [START_REF]Defect Corrections and Multigrid Iterations[END_REF][START_REF] Trottenberg | Multigrid[END_REF]. In our formulation, L ′ corresponds to the fourth-order discrete formulation described by Eqs. ( 9) and [START_REF] Pardhanani | A stable multigrid strategy for convection-diffusion using high order compact discretization[END_REF], which can be represented by a dense matrix, being evaluated with a linear time complexity. The auxiliary discrete problem L refers to the standard five-points stencil resulting in the second-order finite differences discretization of the Poisson equation Eq. ( 1). The defect correction method is based on the following iterative procedure

Lφ (k+1) = Lφ (k) -  L ′ φ (k) -f  , 1 ≤ k ≤ n. ( 11
)
The initial guess function φ (1) is arbitrarily chosen. At a first sight, each defect iteration requires the solution of the auxiliary linear system L. However, as indicated in [START_REF] Trottenberg | Multigrid[END_REF], one iteration of a multigrid solver is sufficient to ensure the convergence of the iterative procedure. This has been confirmed by our experience. The analysis of the convergence factor of this iterative process is based on the spectral radius of the iteration matrix G std = I -L -1 L ′ . The convergence holds for a spectral radius of G strictly lower than 1, and the lowest spectral radius is suitable to obtain a better residual reduction rate, named convergence factor. The analytical calculation of the spectral radius becomes rapidly a cumbersome task. This is addressed by considering the following eigenfunctions [START_REF]Defect Corrections and Multigrid Iterations[END_REF] ϕ(x, y) = e iθ x x/h e iθ y y/h π ≤ θ x , θ y ≤ π .

If the periodic boundary conditions are assumed, the replacing (θ x , θ y ) by (k x π h, k y π h) is equivalent to consider a discrete spectrum, and the rigorous Fourier analysis theory is recovered. A direct computation of the eigenvalues of the five points second order and the fourth order discretization gives, accordingly,

Lϕ = 4 h 2  1 - 1 2 cos θ x - 1 2 cos θ y  ϕ = 4 h 2 f 2th (θ x , θ y ), (13)
and

L ′ ϕ = 4 h 2  sin 2 (θ x /2)  11 12 + 1 12 cos θ x  2 + sin 2  θ y /2   11 12 + 1 12 cos θ y  2  ϕ = 4 h 2 f 4th (θ x , θ y ). (14)
The eigenvalues λ sf of the iteration matrix G sf = I -L -1 L ′ are deduced from the relations

λ sf (θ x , θ y ) = 1 -f 4th (θ x , θ y )/f 2th (θ x , θ y ). (15)
Fig. 1 represents

  λ sf   as a function of θ x and θ y .
This contour plot includes ten isovalues ranging from 0.05 to 0.40, which shows that the highest eigenvalues are located on the nodes (θ x , θ y ) = (0, ±π ) and (θ x , θ y ) = (±π , 0). So the following inequality is readily deduced:

  λ sf   ≤ |λ(0, ±π )| = |λ(±π , 0)| = 11 25 . ( 16
)
Since the eigenvalues are strictly lower than 1, the defected correction method is a convergent iterative procedure. Moreover, the convergence factor q sf = 11/25 is independent of the grid parameter h. This feature is important to design an efficient method of solution because this is a prerequisite for the linear time algorithm. According to Fig. 1, the convergence factor is strongly governed by the high frequencies: the highest eigenvalues are located close to the domain's boundaries. A classical technique to improve the convergence factor consists in adding to the defect correction a smoothing sweep, which can be viewed as a preconditioner. In this case, the iterative procedure Eq. ( 11) becomes

φ = φ (k) -T  L ′ φ (k) -c  , (17)
Lφ (k+1) = Lφ -  L ′ φ -c  . ( 18
)
The influence of the smoothing sweep T onto the convergence factor is examined for two cases: the damped Jacobi and ILU factorization. The relaxation operator modifies the iteration matrix as G = (I -L -1 L ′ )(I -TL). Its spectral radius governs the convergence factor of the whole procedure. Here, the Jacobi relaxation [START_REF]Defect Corrections and Multigrid Iterations[END_REF] is considered. It is based on the diagonal of the matrix L arising from the second-order discretization: ωh 2 /4. Indeed, the eigenvalues λ jac are given by

λ jac (θ x , θ y , ω) = λ sf (θ x , θ y )  1 -ωf 4th (θ x , θ y )  , (19)
where ω is the relaxation factor. In a similar way, the formal eigenfunctions of the ILU smoothing operator are given by [START_REF] Wesseling | An Introduction to Multigrid Methods[END_REF] 4

h 2  1 - 1 2 cos θ x - 1 2 cos θ y + 1 4 + 2 √ 2 cos(θ x + θ y )  = 4 h 2 f ilu (θ x , θ y ), (20)
and, the eigenvalues of the iteration matrix are readily computed as

λ ilu (θ x , θ y , ω) = λ sf (θ x , θ y )  1 -ω f 4th (θ x , θ y ) f ilu (θ x , θ y )  . (21)
Conversely, the calculation of the spectral radius of the iteration matrix is defined as

ρ(ω) = max -π<θ x ,θ y <π   λ(θ x , θ y , ω)   , (22)
where the former notation holds for the damped Jacobi and ILU smoother. The spectral radius is a function of the relaxation factor ω, and might admit an optimal value which minimizes it. In this work, a numerical study is carried out to evaluate the optimal relaxation factor for the both smoothers. The convergence factor Eq. ( 22) is plotted in Fig. 2 for each smoother.

The figure clearly indicates the existence of a minimum which corresponds to the optimal relaxation factor. The spectral radius of the Jacobi smoother is minimal for ω jac = 0.4763, yielding a convergence factor q jac = 0.1635 which is lower than the smoother-free method. The optimal relaxation factor of the ILU smoothing operator is ω ilu = 0.6751, and the associated convergence factor is q ilu = 0.0611. To complete the description of the convergence properties, the contour plot of the eigenvalues for the optimal relaxation factor is represented in Fig. 3.

As the smoother-free method, these plots show that the highest eigenvalues are located in the high frequency range. It should be noted that the introduction of the smoothing operator improves the reduction factor of the defect correction. Thus, the ILU smoother gives a factor of improvement of 8 in comparison with the smoother-free method. We remind that a uniform periodic mesh is assumed in the former analysis, and that the knowledge of the behavior with the Neumann boundary conditions and non-uniform meshes is of major importance in practical applications. This aspect and the implementation verification are addressed in the next section by means of various numerical experiments.

Implementation details

The Neumann boundary conditions lead to a singular problem: the solutions of Poisson's equation are defined up to a constant. There are several methods to consider this special feature. In this work, it is addressed by imposing the Dirichlet condition on a unique boundary node in the auxiliary problem L. Next, both the discrete approximations L and L ′ are scaled by h 2 , and the stopping criterion is the L 2 -norm of the residual. The auxiliary problem L is not solved at each iteration of the defect correction. More particularly, only one iteration of a multigrid method is sufficient to get the convergence of the defect correction. Here, the Algebraic Multigrid method (AMG) provided by the HYPRE library is used. The ILU decomposition smoother, is also provided by the HYPRE library. The fill-in level set to zero (ILU(0)) corresponds to the above mentioned analysis. In the numerical experiments, the fill-in level dependency on the convergence factor is investigated, in this case associated notation is ILU(k).

Numerical experiments

In this section, several numerical experiments are carried out to validate the proposed multigrid defect correction. The fourth-order accuracy of the Poisson solver is first demonstrated. Then, the influence of the smoothing operator and the relaxation parameter on the convergence factor are outlined. In the last section, the Hodge-Helmholtz decomposition is considered to illustrate how the proposed defected correction is able to deal efficiently with the Poisson equation arising from simulations of incompressible fluid flow.

Spatial accuracy

The multigrid defect correction introduces a second-order finite difference discretization for the auxiliary problem, whereas the defect is evaluated via a higher accurate scheme. When the defect is converged, the high-order accurate approximation is expected. This important feature is assessed by computing the order of accuracy p from the numerical errors computed with several mesh sizes. The following set of four problems is considered [START_REF] Zhuang | A high-order fast direct solver for singular Poisson equations[END_REF]:

1. φ (x, y) = (xy) 3.5 [1 -cos (xy)] (five times differentiable), 2. φ (x, y) = x 4.5 + y 4.5 (four times differentiable), 3. φ (x, y) = (x + y) 2.5 sin (x) (three times differentiable), 4. φ (x, y) = (x + y) 2.5 (two times differentiable).

The order of differentiability of these solutions is between 2 and 5. The computations are performed on a unit square domain, and both uniform and non-uniform meshes are considered. The non-uniform mesh is obtained from the following analytical mapping x(X ) = X + β sin [(8X + 1)π ] , [START_REF] Schiestel | A Hermitian-Fourier Numerical method for solving the incompressible Navier-Stokes equations[END_REF] where X stands for the uniform mesh of size

h = n -1 = n -1 x = n -1 y
and where the parameter β is determined to prescribe the smallest space step h min = 0.6h. The metrics x X = x ′ (X) and y Y = y ′ (Y ) are determined from the analytical mesh transform Eq. ( 23). The order of accuracy is computed from the L2-error denoted ϵ n which is the root mean square of the numerical error on a grid of size n. The order of accuracy p is computed as p = log 10  ϵ n 1 /ϵ n 2  / log 10 (n 1 /n 2 ) for the problems (1)-( 4). Section 2 details three variants of the multigrid defect correction: the Jacobi and ILU smoothers and the smoother-free one. Preliminary computations show that the three methods give identical solutions with the machine accuracy precision. Thus, for the sake of clarity, only the numerical errors computed with the smoother-free version are presented. Tables 1 and2 present the numerical error as a function of the mesh size for each problem. Table 1 refers to the error computed on the uniform mesh. The numerical error shows that the order of accuracy achieves the expected fourth-order only with the solution [START_REF] Collatz | The Numerical Treatment of Differential Equations[END_REF].

As noted by Zhuang and Sun [START_REF] Zhuang | A high-order fast direct solver for singular Poisson equations[END_REF], the order decreases accordingly to the order of differentiability of the solution in the other cases. Table 2 lists the numerical errors on the non-uniform meshes.

As for the previous test with uniform mesh, the use of non-uniform meshes and sufficiently fine meshes leads to a similar conclusion: the scheme is fourth-order accurate except for the problem 4. These findings demonstrate that the present defect multigrid method retains the fourth-order accuracy on both uniform and non-uniform meshes. The next section investigates the convergence factor of the three different versions of the defect correction method.

Convergence factor of the multigrid defect correction

Section 2 details the multigrid defect correction and the benefit on the convergence factor in introducing a smoother. The theoretical convergence factor was derived from a uniform mesh and a periodic domain assumptions. The existence of an optimal relaxation factor ω for the Jacobi and the ILU smoothers has been demonstrated. In this section, several numerical experiments are first carried out to validate the theoretical predictions, and then to assess their validity out of the range of validity inherent in postulated assumptions. Thus, the above discussion is based on the convergence properties which include the convergence factor q, and the defect reductions d k at iterations k. The defect reduction is computed as suggested by Trottenberg et al. [START_REF] Trottenberg | Multigrid[END_REF]. The solution u = 0 of the homogeneous Poisson's problem f = 0 is considered, and the initial guess of the multigrid defect correction is initialized with a random field of zero mean value. The empirical convergence factor is computed from the defects at the iterations m 0 = 10 and m 1 = 20, with the relation

q = m 1 -m 0  d m 1 /d m 0 . (24)
The optimal relaxation factor is computed from the minimum of the curve representing the convergence factor q as a function of the relaxation factor ω. Once, the optimal relaxation factor is determined for the Jacobi and the ILU(0) smoothers, their associated defect d k versus the iteration numbers are considered for the two grid sizes 128×128 and 1024×1024. The defect of the smoother-free version ω = 0 is also considered on these two grids. This allows us to evaluate the mesh-independency of the present method. From the relaxation factor and the defect data, we are able to assess the convergence features of the present multigrid defect correction. Fig. 4 presents the convergence features, on a periodic domain with uniform mesh, of the multigrid defect correction for smoother-free, Jacobi and ILU(0) smoothers. It should be noted from Fig. 4(a) that the convergence factor q with respect to the relaxation parameter ω is in good agreement with the predictions of Section 2. The computed convergence factors are close to the analytical curves (full and dashed line) in the whole range of the relaxation factor. This allows us to validate the value of the optimal relaxation factor and tour numerical implementation. Next, Fig. 4(b) outlines the mesh-independency of the defect reduction for the three versions. Indeed, the defects reduction rate is identical regardless of the mesh size. It can be stated that the ILU(0) has the best defect reduction factor, as demonstrated in the former section. Fig. 5 presents the convergence features when the periodic boundary conditions are substituted by Neumann ones. So, the assumptions made for the analytical derivation of the convergence factor are not valid. However, Fig. 5(a) shows that the computed convergence factors fit with theoretical ones for the Jacobi and the ILU(0) smoothers. Fig. 5(b) also shows the mesh-independency of the three implemented versions. From these remarks, the Neumann boundary conditions seem to be effect-less on the convergence features.

Fig. 6 presents the convergence properties for a periodic domain with a non-uniform mesh. In this case this is the mesh which does not fulfill the initial assumptions. This times, the convergence factor at the optimal relaxation factor is altered for both the Jacobi and the ILU(0) smoothers. In fact, the convergence factor is, respectively, about 0.3 and 0.15 for the Jacobi and ILU(0) smoothers in contrast to 0.18 and 0.08 as previously observed. Only the smoother-free version (ω = 0) conserves a convergence close to 0.44. The defect convergence history for the Jacobi and ILU(0) smoothers are also affected: to reach the stopping criteria, few iterations are needed. This feature involves a slight loss of the mesh-independency.

In the same way, Fig. 7 shows the influence on the convergence factor of the Neumann boundary conditions associated with a non-uniform mesh. In this case all the initial assumptions are avoid. The curves follow the same trend observed in the periodic/non-uniform configuration. The convergence factor at the optimal relaxation factor is degraded and a loss of the mesh-independency of the defect reduction is noted. Thus it can be concluded that for the singular Poisson problem the convergence factor is strongly influenced by a non-uniform mesh.

Here, the influence of the fill-in parameter for the ILU smoother is investigated.

The convergence properties as a function of the fill-in parameter are detailed in Fig. 8 when Neumann boundary conditions and non-uniform mesh are considered. An important improvement is noted when the fill-in parameter increases. In Fig. 8(a), ILU(2), ILU(4) and ILU [START_REF] Pardhanani | A stable multigrid strategy for convection-diffusion using high order compact discretization[END_REF] give a better convergence factor than the ILU(0). Moreover, the optimal relaxation factor ω ilu opt = 0.6751 predicted by theory is recovered by using ILU [START_REF] Spotz | High-Order compact finite difference schemes for computational mechanics[END_REF]. In comparison with ILU(2), ILU(4) and ILU [START_REF] Pardhanani | A stable multigrid strategy for convection-diffusion using high order compact discretization[END_REF] do not lead to a more efficient scheme within the meaning of convergence factor. Regardless of the computational cost issue, the following conclusions can be drawn. First, the numerical experiments show that predicted optimal relaxation factors are in agreement with those computed, at least with non-uniform meshes. If non-uniform meshes are considered the optimal relaxation factors are altered, but can be recovered if ILU(4) is considered for the smoothing step. From the theory and the numerical experiments, there follows a clear conclusion: the defect correction associated with the ILU(4) smoother is the best combination in the meaning of the convergence factor. The complementary side of the computational cost issue is addressed in the next section. The computational time of the present multigrid defect correction is analyzed within a framework closed to the computational fluid dynamics: the Hodge-Helmholtz decomposition.

Hodge-Helmholtz decomposition of a vector field

A vector field can be decomposed as the sum of a divergence free vector field and a conservative vector field as follows

u ⋆ = u + ∇φ, (25)
where ∇ • u = 0 and φ stands for a scalar field. In computational fluid dynamics of the incompressible flows, a such decomposition is a key point to obtain an approximate velocity fulfilling divergence-free constraint with the machine accuracy. Taking the divergence of the decomposition Eq. ( 25), the following Poisson equation is formulated

∇ • ∇φ = ∇ • u ⋆ . (26)
Here, homogeneous Neumann boundary conditions are considered. The computational efficiency of the multigrid defect correction is evaluated with the decomposition of the following vector velocity field u ⋆ = (cos(2π x) sin(2π y) + π sin(4π x), sin(2π x) cos(2π y) + π sin(4π y)) .

The Hodge-Helmholtz decomposition is given by  φ = - (cos(4πx) + cos(4π y)) u = (cos(2π x) sin(2π y), sin(2π x) cos(2π y)) .

(

) 28
The Neumann boundary conditions and the non-uniform mesh Eq. ( 23) are also considered. The stopping criterion is set to 10 -10 . In order to investigate the performance of the multigrid defect correction, the absolute numerical error ϵ of the pressure φ, is plotted as a function of the CPU-time in Fig. 9(a).

This figure provides two major informations: the order of accuracy of the scheme and the dependency of the computational time with the grid size. For instance, if a scheme of order p (ϵ ∝ h p ) has a linear complexity algorithm (t cpu ∝ n 2 ∝ 1/h 2 ), then the numerical error varies with the computational time according to ϵ ∝ t -p/2 cpu . Fig. 9(a) illustrates this feature for the smoother-free, Jacobi, ILU(4) and ILU(10) smoothers. First, independently on the smoother, slopes range in 1.8-1.9. This confirms the fourth-order accuracy of the scheme and the linear complexity algorithm of this method of solutions. This figure also illustrates that the most efficient smoother is the ILU(4) one. For instance, an error of the order of 10 -7 is reached with the ILU(4) smoother with more than twice less computational resource than the smoother-free version. It should be also noted that the ILU [START_REF] Pardhanani | A stable multigrid strategy for convection-diffusion using high order compact discretization[END_REF] smoother needs a slightly more computational resource than the ILU(4) one for similar results for the largest mesh. This is in agreement with the conclusion of the previous section. Finally, comparisons between the ILU(4) smoother and a second order discretization of Poisson equation solved by the algebraic multigrid solver of HYPRE are provided in Fig. 9(b). The numerical error as a function of the computational time is plotted. First, the slope associated with the second order is about -1 as expected. Then, for this test case, the plot clearly shows the superiority of the fourth-order defect correction with the ILU(4) smoother associated with a multigrid Poisson solver. Indeed, the error level of the fourth-order is never achieved by the second-order discretization, at least with a reasonable cpu time.

Conclusion

In this work, a multigrid defect correction method is investigated to solve a fourth-order compact scheme approximation of the singular Poisson equation. An eigenvalues analysis of the iteration matrix is developed, and provides several insights on the convergence factor and the mesh-independency of the spectral radius. It is demonstrated that the introduction of a smoothing sweep improves the convergence factor of the iterative method, the latter being dependent on a relaxation factor. A procedure is proposed to determine the optimal value of this relaxation factor both for the Jacobi and ILU smoothers. The better convergence factor has been found for the ILU smoother: a factor of improvement of 8 is determined with the ILU smoother in comparison with the smoother-free version. This analysis is based on the assumption of a periodic domain and uniform mesh. Several numerical experiments were performed. First, the fourth-order accuracy is demonstrated on the uniform and non-uniform meshes with the Neumann boundary conditions. Then, the experiments on the uniform meshes corroborate the existence of an optimal relaxation factor for each smoother and the mesh-independency of the convergence factor. It has been found that a non-uniform mesh alterates the optimal convergence factor and the meshindependency if the Jacobi and the ILU(0) smoothers are used. Furthermore, the smoother ILU(4) allows us to recover the theoretical convergence factor and the mesh independency. Finally, the computational performance is assessed by means of the Hodge-Helmholtz decomposition test case. It clearly shows that within the meaning of computational effort, the fourth-order accurate scheme is better than second order, thanks to the proposed efficient multigrid defect correction.

Fig. 1 .

 1 Fig. 1. Eigenvalues plot of the iteration matrix for the smoother-free defect correction.

Fig. 2 .

 2 Fig.2. Convergence factor q versus the relaxation factor ω.

Fig. 3 .

 3 Fig. 3. Eigenvalues plot of the iteration matrix for the Jacobi (a) and the ILU (b) smoothers.

Fig. 4 .

 4 Fig. 4. Convergence factor versus relaxation factor (a), defect reduction versus the iteration number (b) with periodic boundary conditions and uniform meshes.

Fig. 5 .

 5 Fig. 5. Convergence factor versus relaxation factor (a), defect reduction versus the iteration number (b) with Neumann boundary conditions and uniform meshes.

Fig. 6 .

 6 Fig. 6. Convergence factor versus relaxation factor (a), defect reduction versus the iteration number (b) with periodic boundary conditions and nonuniform meshes.

Fig. 7 .

 7 Fig. 7. Convergence factor versus relaxation factor (a), defect reduction versus the iteration number (b) with periodic boundary conditions and nonuniform mesh.

Fig. 8 .

 8 Fig. 8. Influence of the fill-in parameter; convergence factor versus relaxation factor (a), defect reduction versus the iteration number (b).

Fig. 9 .

 9 Fig. 9. Numerical error on φ versus cpu time, for different smoother (a), comparison with the second-order discretization (b).

Table 1

 1 Numerical L2-error and order of accuracy with a uniform mesh.

	N	Problem 1		Problem 2		Problem 3		Problem 4
		ϵ	p	ϵ	p	ϵ	p	ϵ	p
	128	2.63e-09	0.00	2.42e-08	0.00	4.77e-09	0.00	1.40e-07	0.00
	256	1.60e-10	4.04	2.04e-09	3.57	3.89e-10	3.62	2.48e-08	2.50
	512	9.85e-12	4.02	1.74e-10	3.55	3.25e-11	3.58	4.40e-09	2.50
	1024	6.12e-13	4.01	1.50e-11	3.54	2.75e-12	3.56	7.78e-10	2.50
	2048	2.12e-14	4.85	1.29e-12	3.53	2.16e-13	3.67	1.38e-10	2.50

Table 2

 2 Numerical L2-error and order of accuracy with a non-uniform mesh.

	N	Problem 1		Problem 2		Problem 3		Problem 4
		ϵ	p	ϵ	p	ϵ	p	ϵ	p
	128	6.51e-06	0.00	1.49e-05	0.00	1.28e-06	0.00	6.32e-07	0.00
	256	3.99e-07	4.03	8.74e-07	4.09	3.95e-08	5.02	1.47e-07	2.10
	512	2.44e-08	4.03	5.20e-08	4.07	1.93e-09	4.35	3.82e-08	1.95
	1024	1.51e-09	4.02	3.14e-09	4.05	1.30e-10	3.89	8.68e-09	2.14
	2048	9.33e-11	4.01	1.92e-10	4.04	8.73e-12	3.90	1.91e-09	2.18


Florent Goujon

Aziz Ghoufi

Patrice Malfreyt
email: patrice.malfreyt@uca.fr

Size-effects on the surface tension near the critical point : Monte Carlo simulations of the Lennard-Jones fluid

Keywords: Surface tension, critical point, Monte Carlo simulations, slab-based tail method

We report Monte Carlo (MC) simulations of the Lennard-Jones (LJ) fluid at the liquid-vapor interface in the critical region. A slab-based tail method is associated with the MC simulations to approach as close as possible the critical point (T * = 0.98T * C ). We investigate then the impact of system-sizes on the surface tension and coexisting densities by considering very large box dimensions for which the surface tension is independent of system-sizes at low temperatures.

Introduction

Most of the phenomena in surface science adhesion, wetting, lubrication [START_REF] Mendonca | Interactions and ordering of ionic liquids at a metal surface[END_REF][START_REF] Mendonca | Nonequilibrium molecular simulations of new ionic lubricants at metallic surfaces: Prediction of the friction[END_REF], nucleation [START_REF] Auer | Suppression of crystal nucleation in polydisperse colloids due to increase of the surface free energy[END_REF][START_REF] Prestipino | Systematic improvement of classical nucleation theory[END_REF] involve the combination of liquid-liquid, liquid-vapor, and liquid-solid interfaces. An understanding and control of the surface tension is at the heart of many important industrial and practical processes. The molecular simulations have been extensively used this last decade [START_REF] Ghoufi | Computer modelling of the surface tension of the gas-liquid and liquid-liquid interface[END_REF] to reproduce the interfacial tension at different thermodynamic conditions and to describe the interface region in terms of specific arrangements [START_REF] Ferrando | A transferable force field to predict phase equilibria and surface tension of ethers and glycol ethers[END_REF][START_REF] Mendonca | Nonequilibrium molecular simulations of new ionic lubricants at metallic surfaces: Prediction of the friction[END_REF].

Since 1974 and so far, the calculation of the surface tension from two-phase simulations has been an area of active research. This long-standing activity is explained by the fact that a number of factors such as the surface area (See Figures 1a and1b) [START_REF] Orea | Oscillatory surface tension due to finite-size effects[END_REF][START_REF] Gonzalez-Melchor | Stress anisotropy induced by periodic boundary conditions[END_REF][START_REF] Errington | Calculation of surface tension via area sampling[END_REF][START_REF] Biscay | Calculation of the surface tension from monte carlo simulations: Does the model impact on the finite-size effects?[END_REF][START_REF] Werth | The inflence of the liquid slab thickness on the planar vapor-liquid interfacial tension[END_REF][START_REF] Goujon | Controlling the long-range corrections in atomistic monte carlo simulations of two-phase systems[END_REF], the cutoff radius and the corresponding long-range corrections [START_REF] Salomons | Atomistic simulation of liquid-vapour coexistence:binary mixtures[END_REF][START_REF] Blokhuis | Tail corrections to the surface tension of a Lennard-Jones liquid-vapour interface[END_REF][START_REF] Guo | Long range corrections to thermodynamic properties of inhomogeneous systems with planar interfaces[END_REF][START_REF] Goujon | Direct monte carlo simulations of the equilibrium properties of n-pentane liquid-vapor interface[END_REF][START_REF] Janecek | Long range corrections in inhomogeneous simulations[END_REF][START_REF] Ibergay | Molecular simulations of the n-alkane liquid-vapor interface: Interfacial properties and their long range corrections[END_REF][START_REF] Shen | Comparative study of the effect of tail corrections on surface tension determined by molecular simulation[END_REF][START_REF] Míguez | Influence of the long-range corrections on the interfacial properties of molecular models using Monte Carlo simulation[END_REF], the truncation of the potential and the force as well as the method used for the calculation of the surface tension [START_REF] Kirkwood | The statistical mechanical theory of surface tension[END_REF][START_REF] Irving | The statistical mechanical theory of transport processes .IV. The equations of hydrodynamics[END_REF][START_REF] Buff | Some considerations of surface tension[END_REF][START_REF] Buff | Spherical interface. ii. molecular theory[END_REF][START_REF] Gloor | Test-area simulation method for the direct determination of the interfacial tension of systems with continuous or discontinuous potentials[END_REF][START_REF] Ghoufi | Multiple histogram reweighting method for the surface tension calculation[END_REF][START_REF] Laird | Determination of the solid-liquid interfacial free energy along a coexistence line by gibbs-cahn integration[END_REF] make the simulated properties dependent on initial conditions. All these issues have been addressed in a recent review [START_REF] Ghoufi | Computer modelling of the surface tension of the gas-liquid and liquid-liquid interface[END_REF]. A number of recommendations have been proposed to resolve these issues leading now to an accurate calculation of the surface tension for reduced temperatures less than 0.85 T * C . For higher temperatures, the situation is far from being identical. Indeed, the cutoff-dependence on the surface tension and coexisting densities is still stronger [START_REF] Goujon | Can we approach the gas-liquid critical point using slab simulations of two coexisting phases ?[END_REF]. Recently, we have shown that it is possible to approach the gas-liquid critical point of the Lennard-Jones fluid by performing two-phase simulations of a slab geometry. These simulations have used a slab-based tail [START_REF] Janecek | Long range corrections in inhomogeneous simulations[END_REF] methodology with large cutoff values. The main conclusions of this work are that Monte Carlo simulations with the addition of long-range corrections to the energy during the course of the simulations using a reduced cutoff of r * c = 6.4 reproduce quantitatively the dependencies of the densities and surface tensions on the temperature close to the critical point. In the case of the modeling of methane, the real value of the cutoff corresponds to r c = 6.4σ = 24 Å which is twice that used in most atomistic simulations [START_REF] Ferrando | A transferable force field to predict phase equilibria and surface tension of ethers and glycol ethers[END_REF][START_REF] Ghoufi | Multiple histogram reweighting method for the surface tension calculation[END_REF].

Whereas a reduced cutoff of r * c = 3.0 is large enough to avoid any cutoff-dependence of the surface tension with a slab-based tail method for temperatures far from the critical point (see Ref.12), an increase to r * c = 6.4 is required to provide reliable surface tensions and coexisting densities in the critical region [START_REF] Goujon | Can we approach the gas-liquid critical point using slab simulations of two coexisting phases ?[END_REF]. A key question arises now : for this large cutoff value of r * c = 6.4, how do the size-effects affect the values of surface tension and coexisting densities of the critical region ? The size-effects will be evaluated trough the variations of the interfacial area A = L x × L y and the longitudinal dimension L z where z represents the direction normal to the interface. It has been established, at relatively low temperatures, that the size-effects led to oscillatory behaviors of the surface tension with the surface area (Figure 1a) and the L * z dimension (Figure 1b). To address this issue, we propose here to investigate LJ liquid-vapor systems close to the critical point with large surface areas and longitudinal dimensions. Figures 1a and1b show the values of surface areas (L * x1 , L * x2 , L * x3 ) and longitudinal dimensions (L * z1 , L * z2 , L * z3 ) considered here. The details of the simulation in terms of box sizes and number of LJ particles are given in Table 1.

Methodology

All the properties reported here refer to the Lennard-Jones (LJ) particle and are expressed in reduced LJ units where T * = k B T/ , γ * = γσ 2 / , L * = L/σ, p * = pσ 3 / , ρ * = ρσ 3 where T * , γ * , L * , p * , ρ * represent the reduced temperature, surface tension, length, pressure and density, respectively. MC simulations were performed in the constant-NVT ensemble. All simulations were organized in cycles. Each cycle consisted of N translations. The equilibration phase was composed of 500 000 cycles and the production phase of 2 × 10 6 cycles up to 5×10 10 translation moves. The maximum displacement which was adjusted during the equilibrium phase to give an acceptance ratio of 0.4, was found to be close to 0.23 at T * = 1.26 and to 0.26 at T * = 1.28 in reduced LJ units. The thermodynamic and mechanical properties were calculated every 20 cycles requiring the storage of 100 000 configurations. The statistical errors for these properties were estimated using 5 superblocks averages of 20 000 configurations. The MC calculations were carried out at two reduced temperatures T * = 1.26 and 1.28 for a cutoff radius r * c = 6.4. In the Janeček approach, the long-range correction to the truncated potential is added to the truncated u ST (r i j ) potential. Considering a system of N atoms, the total configurational energy U TOT is

U TOT = N-1 i=1 N j=i+1 u ST (r i j ) + U LRC (1)
where r i j is the pair separation distance and u ST (r i j ) is the spherically-truncated Lennard-Jones potential defined by

u ST (r i j ) = u LJ (r i j ) r i j < r c 0 r i j ≥ r c . ( 2
)
u LJ is the Lennard-Jones (LJ) potential defined by and σ corresponding to energy and size parameters, respectively. The simulation box is divided into N s slabs of width ∆z = 0.13 (in LJ units). Each slab, which is parallel to the interface, has a volume V s = L x L y ∆z and z k defines the centre of the kth slab. In Equation ( 1), U LRC is defined as

U LRC = 1 2 N s k=1 u lrc (z k ) (3)
where the long-range correction energy of the slab k is

u lrc (z k ) = ρ(z k )V s N s j=1 ρ(z j )w(|z j -z k |)∆z (4)
where the sum is over all the slabs in the box. ρ(z k ) defines the density number of the slab k. The contribution w(ξ) = w(|z jz k |) is calculated by assuming a uniform distribution of atoms in the slab. For the Lennard-Jones potential, the function w(ξ) is

w(ξ) =                    4π σ 2        1 5 σ r c 10 - 1 2 σ r c 4        for ξ ≤ r c 4π σ 2        1 5 σ ξ 10 - 1 2 σ ξ 4        for ξ > r c (5)
Note that Equation (3) defines the long-range correction energy as a sum over slabs whereas the original paper [START_REF] Janecek | Long range corrections in inhomogeneous simulations[END_REF] considers the correction in terms of a molecule i at z i . Finally there is an improvement to the original version of Janeček's method proposed by MacDowell and Blas [START_REF] Macdowell | Surface tension of fully flexible Lennard-Jones chains: role of long-range corrections[END_REF], which avoids the decomposition of the long-range correction of the energy into local contributions and increases the speed of the calculation. The long-range corrections to the energy are applied at each MC move although the density profile used to calculated the long-range corrections is updated every 10 MC cycles [START_REF] Goujon | Controlling the long-range corrections in atomistic monte carlo simulations of two-phase systems[END_REF]. The intrinsic part of the surface tension, γ I , calculated using the Irving and Kirkwood [START_REF] Kirkwood | The statistical mechanical theory of surface tension[END_REF][START_REF] Irving | The statistical mechanical theory of transport processes .IV. The equations of hydrodynamics[END_REF][START_REF] Rowlinson | Molecular Theory of Capillarity[END_REF][START_REF] Walton | The pressure tensor at the planar surface of a liquid[END_REF][START_REF] Walton | The pressure tensor at the planar surface of a liquid[END_REF][START_REF] Martinez-Ruiz | Effect of dispersive long-range corrections to the pressure tensor. the vapour-liquid interfacial properties of the lennard-jones system revisited[END_REF]] definition, uses the components p N (z) and p T (z) of the pressure tensor as a function of z

γ I = L z /2 -L z /2 p N (z) -p T (z) dz. (6)
Equation ( 6) is a mechanical definition of γ based upon the force acting across a unit area in the z-plane for one interface. Due to the truncation of the potential, the long range corrections to the surface tension are calculated with the Janeček approach as

γ lrc (z k ) = V s A p N,lrc (z k ) -p T,lrc (z k ) (7)
The normal and tangential components of the long-range corrections of the pressure tensor are given by

p N,lrc (z k ) = ρ(z k ) N s j=1 ρ(z j )π zz (|z j -z k |)∆z (8)
and

p T,lrc (z k ) = 1 2 p xx,lrc (z k ) + p yy,lrc (z k ) = ρ(z k ) N s j=1 ρ(z j ) 1 2 π xx (|z j -z k |) + π yy (|z j -z k |) ∆z (9)
where the operational expressions of the functions π xx (ξ), π yy (ξ) and π zz (ξ) can be found in Ref.12.

Table 1. Reduced surface tensions of the LJ fluid calculated at different surface areas A and L * z dimensions where N represents the total number of particles in the box. The different contributions of the surface tension, i.e, intrinsic, long range and total parts are given in LJ reduced units. < γ > is calculated using the Irving-Kirkwood definition. The subscripts give the accuracy of the last decimal(s), i.e., 0.104 37 means 0.104 ± 0.037. The surface area A 1 corresponds to 12.9 × 12.9 and the L z1 dimension to 200.9 in reduced units.

L * x × L * y L * z N Box γ * I γ * LRC γ * T * =

Results and discussions

Figure 2 shows the profiles of the difference between the normal and tangential pressure components at T * = 1.28 for both surface areas. First, we observe that the profiles of p * N (z * )p * T (z * ) and γ * (z * ) are in line with an expected behavior of liquid-vapor interface at mechanical equilibrium. However, the peaks at the interface are much less marked at T * = 1.28. The profiles of γ * (z * ) shows that the bulk phases do not provide any contribution to the surface tension and lead to flat portions on the profiles. Only the interface contributes to this property indicating that it is stable and well-located. In other words, it is possible to simulate explicit interfaces of LJ fluid very close to the critical point (0.98 T * C ). T * C was estimated to 1.31 by GCMC-FSC simulations [START_REF] Potoff | Surface tension of the three-dimensional Lennard-Jones fluid from histogram reweighting Monte Carlo simualtions[END_REF]. We also observe that a larger system-size (A 3 L z1 ) leads to a different surface tension. The difference of surface tension between these two surface areas comes from the contribution of the interfaces, with peaks of slightly larger magnitudes for A 1 L z1 . The fluctuations in the bulk phases are identical between the two systems-sizes. As result, the investigation of the impact of the system-sizes on the surface tension is relevant in the critical region at these two temperatures.

We now focus on the heart of this work, namely the effects of the system-sizes close to the critical point. Figure 3 and Table 1 show the dependencies of the surface tension on the interfacial area (L * x = √ (A)) and on the longitudinal dimension (L * z ) at both temperatures T * = 1.26 and T * = 1.28. We also report the impact of the box dimensions on the liquid and vapor densities. First, we check that the surface tensions and coexisting densities respect the temperature dependence even if the variation of temperature (∆T * = 0.02) is weak. The decreases of the surface tension and liquid density and the increase of the vapor density are always observed despite the system-size. Second, we are able to provide very small values of surface tensions (smaller than 0.027 in reduced units) with a slab-based tail method. The long range correction to the surface tension contributes to about 15 % of the total value. This tail contribution can never be neglected even at high temperatures. With intensive MC calculations, we obtain very small values of fluctuations in the surface tension making the comparison of the interfacial properties with different system-sizes relevant. Third, for lower temperatures (T * < 0.85T * C ), the maximum variation of the surface tension in the first oscillations of γ * = f (L *

x ) reaches 30% with the LJ truncated potential and decreases to 15% with the slab-based tail methodology used here (see Figure 1). When the surface tension does not change any more with L *

x (from L * x > 11), the values of γ are identical within 0.8%. As shown in part b) of Figure 1, these same variations in γ * = f (L * z ) are about 20% with the Janeček method and the surface tension is independent on L * z from L * z = 60. The dimensions of the simulation boxes simulated here (see Figure 1 and Table 1) are significantly larger than those recommended to avoid any size-effects at smaller temperatures. Nevertheless, the size-effects are still noticeable for these very large box dimensions close to the critical point. Actually, parts a) and c) of Figure 3 show variations of surface tension in the order of 35% with the interfacial area (Figure 3a) and the longitudinal dimension (Figure 3c). It means that it is not possible to reach a surface tension independent on the system-size close to the critical region within the range of system-sizes investigated here. As expected as we approach T C , the size-effects become more important. We can draw the same conclusions with the coexisting densities. The liquid density supports variations in the order of 20% whereas the vapor density changes by 45% with L *

x and L * z . This is an interesting result since it has been shown [START_REF] Chen | Area dependence of the surface-tension of a Lennard-Jones fluid from molecular-dynamics simulations[END_REF][START_REF] Orea | Oscillatory surface tension due to finite-size effects[END_REF] in the past that the coexisting densities were not dependent on system-sizes for temperatures lower than 0.90T * C . The density of Figure 4 confirm the dependencies of the coexisting densities on system-sizes. For the same interfacial area A 1 , the density profiles show a clear dependence on L * z . The density profiles can be fitted to a hyperbolic tangent function even if the agreement is less satisfactory in the interfacial region. The density profiles in the bulk regions are less flat than at lower temperatures and the interface is more diffuse. This is the price to be paid when we approach the critical point. In addition, a weak variation (∆T * = 0.2) of temperature is clearly observable on the profiles. Anyway, these profiles are rather well defined to allow an estimate of the coexisting densities. These density profiles show that the interfacial thickness increases significantly as we approach T C . It is well-known that the interfacial width coincides with the bulk correlation length [START_REF] Stillinger | Capillary waves and the inherent density profile for the liquid-vapor interface[END_REF][START_REF] Russier | Statistical field theory for liquid vapor interface[END_REF][START_REF] Mayoral | Hyperscaling relationship between the interfacial tension of liquids and their correlation length near the critical point[END_REF] and thus diverges at T * C as d * ∼ (T * C -T * ) ν . We have shown in a recent paper [START_REF] Goujon | Can we approach the gas-liquid critical point using slab simulations of two coexisting phases ?[END_REF] that this divergence of the interfacial thickness starts at T * = 1.28 and that the slope ν obtained by Monte Carlo simulations was equal to -0.52 ± 0.02 whereas the mean-field theory [START_REF] Rowlinson | Molecular Theory of Capillarity[END_REF] predicts ν = -1/2. In Figure 4b, we have represented the interfacial thickness d * as a function of T * C -T * . The values were taken from a recent work [START_REF] Goujon | Can we approach the gas-liquid critical point using slab simulations of two coexisting phases ?[END_REF] for which the dimensions of the system match are a bit smaller than those of A 1 L z1 . The values d * calculated here are reported on this figure for different box dimensions. We observe a good agreement with the value of d * calculated in A 1 L z1 and strong deviations when the longitudinal box dimension L * z is changing. Actually, the strong dependence of the interfacial thickness on L * z , illustrated in part a) of this figure, is confirmed on this curve. We cannot conclude on the scaling exponent since additional simulations must be carried out with the box dimensions used in this paper. To complete this comparison, we have represented in Figure 4c the surface tension γ * as a function of T * C -T * . Our MC simulations [START_REF] Goujon | Can we approach the gas-liquid critical point using slab simulations of two coexisting phases ?[END_REF] predict a slope µ of 1.12 ± 0.0.2 in γ * ∼ (T * C -T * ) whereas µ = 1.26 in the mean-field theory [START_REF] Rowlinson | Molecular Theory of Capillarity[END_REF]. We observe here that changing L z impacts strongly on the values of γ * when approaching T * C .

Summary

The combination of a slab-based tail method with Monte Carlo simulations has been used to investigate the dependencies of the coexistence properties close to the critical point at T * = 0.98T * C . It has been established at low temperatures (T * < 0.85T * C ) that the surface tension shows an oscillatory behavior with the interfacial area for interfacial areas smaller than (L *

x × L * y < 11 × 11). In this work, we took the route of using large box dimensions for investigating the critical region. The main conclusion is that the surface tension and the coexisting densities are system-sizes dependent even with very large box dimensions. The surface tension is dependent on both the interfacial area and the longitudinal dimension. The variation of the surface tension as a function of L *

x or L * z can reach 35%. The liquid and vapor density significantly change with the box dimensions. None of these dependencies can be attributed to large statistical fluctuations and unstable interfaces. All these factors have been controlled by doing intensive Monte Carlo simulations with a slab-based tail method. We also report a strong dependence of the surface tension and the interfacial thickness on the longitudinal dimension when approaching T * C . 1. We add for discussion the fits of these density profiles to a hyperbolic tangent function [START_REF] Goujon | Can we approach the gas-liquid critical point using slab simulations of two coexisting phases ?[END_REF]; Log-log plots of the b) interfacial width d * and c) reduced surface tension γ * versus T * C -T * where T * C was taken to 1.31 [START_REF] Potoff | Surface tension of the three-dimensional Lennard-Jones fluid from histogram reweighting Monte Carlo simualtions[END_REF]. The dotted lines represent the least-square fits [START_REF] Goujon | Can we approach the gas-liquid critical point using slab simulations of two coexisting phases ?[END_REF]. The interfacial widths and surface tensions calculated in this work are represented for different box dimensions as indicated in the legend.

• We investigate the size-effects on the surface tension in the critical region • We approach the critical as close as possible with a slab-based tail Monte Carlo method. • The surface tension and coexisting densities are system-sizes dependent in the critical region.

Figure 1 .

 1 Figure 1. Surface tensions [10] calculated with the truncated LJ potential and the LJ potential modified with the Guo and Lu slab-based tail methodology at T * = 0.8 for a) different interfacial areas and b) L * z dimensions. The dotted vertical lines represent the values of interfacial area and of longitudinal dimensions L * z investigated here.

)Figure 2 .

 2 Figure 2. Profiles of the difference p * N (z * )p * T (z * ) along the direction z * normal to the interface at T * = 1.28. The integral of these profiles are represented by dotted lines on the right axis.

Figure 3

 3 Figure 3. a) and c) Reduced surface tensions and b) and d) coexisting densities as a function of the interfacial area (L * x = √ (A)) and the longitudinal dimension L * z at two reduced temperatures T * = 1.26 and 1.28. In parts b) and d), the liquid and vapor densities are represented by red points and blue points on the right and left axes, respectively.

Figure 4

 4 Figure 4. a) Density profiles of the LJ fluid at T * = 1.26 and T * = 1.28 for different system-sizes as indicated in Table1. We add for discussion the fits of these density profiles to a hyperbolic tangent function[START_REF] Goujon | Can we approach the gas-liquid critical point using slab simulations of two coexisting phases ?[END_REF]; Log-log plots of the b) interfacial width d * and c) reduced surface tension γ * versus T * C -T * where T * C was taken to 1.31[START_REF] Potoff | Surface tension of the three-dimensional Lennard-Jones fluid from histogram reweighting Monte Carlo simualtions[END_REF]. The dotted lines represent the least-square fits[START_REF] Goujon | Can we approach the gas-liquid critical point using slab simulations of two coexisting phases ?[END_REF]. The interfacial widths and surface tensions calculated in this work are represented for different box dimensions as indicated in the legend.


Julie Golliot
email: julie.golliot@ensc.fr

Michèle Timsit
email: m.timsit@ramsaygds.fr

Cathy Herrera
email: c.herrera@ramsaygds.fr

Elodie Fontugne Clinique

Rgds Provence-Bourbonne

Alexandre Abellard
email: alexandre.abellard@univ-tln.fr

Michel Durampart
email: michel.durampart@univ-tln.fr

Specification, Use and Impact of the Persuasive Serious Game S'TIM in a Rehabilitation Process for Patients with Dysexecutive Syndrome

Keywords: Social and professional topics → Assistive technologies, People with disabilities, • Applied computing → Computerassisted instruction, Serious-Game,Rehabilitation,Dysexecutive Syndrome,Tactile Table,Cognitive Impact,Organisational Change

The aim of this multidisciplinary study is to specify and develop a Serious-Game (SG) to immerge patients with a dysexecutive syndrome in a virtual word. With an elaborate scenario and various challenges, the SG we named S'TIM is used on a robotised and easy-to-use touch table. The high stakes for patients are firstly to break anosognosia and intrinsically motivate patients to implicate them in their rehabilitation; Secondly to enable them to reach a sufficient metacognition level to develop their own strategy and select the most relevant in each context. Finally, to facilitate these strategies transfer in daily-life. Changes in organisations will also be observed.

INTRODUCTION

New technologies are more and more developed and used by people of all ages. In a context of cognitive rehabilitation, most of the tools are traditional with activities on papers or table games. Improvements are observed but the transfer in daily-life is very limited and patients don't recover their autonomy. The ecological validity of conventional tests (with only one task mostly initiated by therapists) seems to be an issue, since they don't take care of the patients' environment. This artificial situation doesn't reproduce complexity of daily life activities and patients don't use their executive functions in order to elaborate strategies and adapt their behaviour. Ecological practices have been elaborated but they are hard to use (lack of therapists, legislation). Therefore, rehabilitation specialists are looking for a new tool to help the patients to recover their capacities and daily-life.

CONTEXT

In this PhD Project (preparation stage), we focus on patients with dysexecutive syndrome which can be observed after a cerebrovascular accident, a head injury with frontal lobe damage or degenerative diseases, amongst other things. Executive functions are implicated in cognitive processing requiring coordination of several sub-processes to achieve a particular goal [START_REF] Godefroy | Syndromes frontaux et dysexécutifs[END_REF] when reflex actions are insufficient. Executive troubles impair capacities of planning, strategies elaboration, updating, shifting, inhibition and attention. There are cognitive, emotional and behavioural consequences. Patients lose their autonomy and are unable to recognize their disorders (anosognosia) which limits their implication in rehabilitation.

MATERIAL AND METHOD

The aim of this research program is to specify and develop a persuasive serious-game to immerge patients in a virtual world with an elaborate scenario with various challenges. We base our research on the self-determination theory highlight three basic needs: competence, relatedness and mostly autonomy [START_REF] Ryan | Overview of self-determination theory: An organismic dialectal perspective[END_REF]. According to this psychological theory, patient intrinsic motivation is a central element to facilitate socio-cognitive and socioaffective changes, since motivation is higher when people do not act for rewards. The device has to be designed in order to make the most of seriousgames characteristics. Interactivity, feedbacks, playful challenges, based on rules and with a clear goal [START_REF] Alvarez | Introduction au Serious Game[END_REF]. As it uses persuasive communication according to Teresa de la Hera Conde-Pumpido PhD Reports DH'18, April 23-26, 2018, Lyon, France

RESULTS

The chosen device is a robotised touch table of 46 inches which allows all degrees of movements (low-high, horizontal-vertical or any another degree). It is easy-to-use by people of all age, with or without wheelchair, sight troubles, and so on. Other applications are already used on the device and therapists haven't noticed difficulties with any patients. Each detail of the serious-game and its use was discussed with all the experts listed above to take into account the various needs. Furthermore, this device is designed for patients and therapists ; thus, their implication since the beginning of the process reduces objections and rejections. With therapists' care, patients should be implicated, acquire pertinent strategies, be intrinsically motivated (self-determination theory), improve their self-respect, their self-efficiency and be more empathic. Collected data will enable qualitative and quantitative analyses. The serious-game development is made on Unity3D, a game development platform. A first prototype using the main decided mechanisms has been developped. Patients and therapists' feedbacks give tools to improve the future version (Fig. 1). With therapists' care, S'TIM will transfer patients at the heart of their rehabilitation.

LIMITS AND PERSPECTIVES

This technology presents some limits. Patients could reject it or have worse results than expected if they have mental difficulties. Moreover, only some therapists were implicated in the development, so others could reject the device if the accompanying is insufficient. For patients, three stakes are expected and analyzed in clinical studies. First, to break the anosognosia and intrinsically motivate patients to implicate them in their rehabilitation. Second, to enable them to reach a sufficient metacognition level to develop their own strategy and select the most relevant in each context. Finally, to facilitate these strategies transfer in daily-life.

Standardized psychometric tests will be conducted before the study, after the two-months rehabilitation training and three months later in order to validate the serious-game pertinence (GREFEX set for cognitive evaluation, HAD scale for thymic criterion and IDSC scale for the behavior). Therapists will keep their fundamental role. They will not have to give some background but to contextualize knowledge to favour transfer (didactic action). Patients should be more implicated, intrinsically motivated, participate more [START_REF] Groleau | L'articulation technologies -organisations: des pistes pour une approche communicationnelle[END_REF] and be less dependent.

Standardized situations with the SG will allow the study of quantitative data and to objective changes. Finally, articulation between goals, actor system and technic dimension will create organisational change [START_REF] Barroy | Les dispositifs numériques de formation: un enjeu de changements organisationnels[END_REF]. The introduction of this technology has to be guided to avoid a situation with more disorders than solutions and optimize the device use.

CONCLUSION

This new, innovating and multidisciplinary study presents high stakes for the future of patients with dysexecutive syndrome. Theories of information and communication sciences, educational sciences, cognitive sciences, psychology and neuropsychology are implicated in the design and the use of this new technology by rehabilitation centers already established. With therapists' implication, S'TIM will increase patient autonomy and transfer them at the heart of their rehabilitation and give them a chance to go back to their daily life.

Figure 1 :

 1 Figure 1: Screenshot of S'TIM.

View publication stats View publication stats