Erwan Deriaz

High-order Adaptive Mesh Refinement Multigrid Poisson Solver in any Dimension

Keywords: Poisson solver, Adaptive Mesh Refinement, Multigrid, Compact Finite Difference, High order scheme Mathematics Subject Classification (2000) 65N50, 65N55, 65Y02

A numerical method to solve the d-dimensional Poisson equation with 2pth-order accuracy for arbitrary d and p integers is proposed in the Adaptive Mesh Refinement framework. This method relies on vertex-centered mesh refinement and interpolation, on compact finite difference schemes and on multigrid algorithm. It is compared to other existing methods. And, in extensive numerical experiments, a sixth-order version of it in dimensions two to six and a tenth-order version in dimension three are tested.

Introduction

Solving the Poisson Equation in dimension d:

∆ u = v, with ∆ = d ∑ i=1 ∂ 2 i (1.1)
plays a fundamental role in the Physics modelization and simulation. Amongst many of its applications, this ubiquitous equation appears: E. Deriaz Institut Jean Lamour, Matériaux-Métallurgie-Nanosciences-Plasmas-Surfaces UMR 7198 -CNRS -Université de Lorraine, Campus Artem, 2 allée André Guinier -BP 50840, 54011 Nancy Cedex Tel.: +33(0)372742557 E-mail: erwan.deriaz@univ-lorraine.fr in the computation of potential forces which derive from gradients of potentials of the kind a = ∇ϕ with ∆ ϕ = ±ρ where ρ represents a density of mass or of charges and ϕ a gravitational field (plus case) or an electric field (minus case) for instances, to numerically solve the Incompressible Navier-Stokes Equation where the pressure is obtained through ∆ P = -ρdiv ((u • ∇)u), to model dispersive effects as in the Heat Equation ∂ t u = ∆ u which needs to be solved with implicit numerical methods of the kind (Id -ν∆ )u n+1 = u n with ν > 0, related to Poisson Equation.

Some of the solutions to these equations develop singularities in space and/or in time.

In this case the Adaptive Mesh Refinement (AMR) methods whose principle consists in automatically adapting the grid to the numerical solution, give access to all the scales, from the largest to the smallest with discrepancies of few orders of magnitudes as in the cosmological simulations [START_REF] Teyssier | Cosmological hydrodynamics with adaptive mesh refinement. A new high resolution code called RAMSES[END_REF]. Compared to the uniform grid methods, these methods cut down the memory needs of the simulations, and cut down also their computational times when their compression rates drop below 10% [START_REF] Cottet | Semi-Lagrangian particle methods for high-dimensional Vlasov-Poisson systems[END_REF][START_REF] Deriaz | Six-Dimensional Adaptive Simulation of the Vlasov Equations Using a Hierarchical Basis[END_REF]. They allow to treat boundary layers [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF], sudden disruptions [START_REF] Sarto | A multigrid AMR algorithm for the study of magnetic reconnection[END_REF] and to reduce the cost of large simulations in cosmology [START_REF] Teyssier | Cosmological hydrodynamics with adaptive mesh refinement. A new high resolution code called RAMSES[END_REF], in ab initio simulations [START_REF] Genovese | Efficient and accurate three-dimensional Poisson solver for surface problems[END_REF] and in geological simulations [START_REF] Kronbichler | High accuracy mantle convection simulation through modern numerical methods[END_REF][START_REF] Isaac | Recursive Algorithms for Distributed Forests of Octrees[END_REF].

Nevertheless these new possibilities offered by the AMRs come at the cost of more complexity in the implementation because of the versatility of the data structure, and more numerical difficulties because of the non-uniformity of the grid. In particular, as the discontinuities of the space step trigger discontinuities in the differentials of the error functions, composing discrete differential operators decreases the final order of accuracy -this effect does not appear in uniform grids where the final order of accuracy of the composed differential operator is usually equal to the minimum order of its differential components. Hence, in AMRs, the need for composed discrete operators [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF] or high-order schemes [START_REF] Barad | A fourth order accurate adaptive mesh refinement method for Poisson's equation[END_REF][START_REF] Sarto | A multigrid AMR algorithm for the study of magnetic reconnection[END_REF]. The parallelization of AMRs, in particular for distributed memory, constitutes a real challenge overcome only by a few number of codes [START_REF] Teyssier | Cosmological hydrodynamics with adaptive mesh refinement. A new high resolution code called RAMSES[END_REF][START_REF] Brix | Adaptive multi-resolution Methods: Practical issues on Data Structures, Implementation and Parallelization[END_REF][START_REF] Isaac | Recursive Algorithms for Distributed Forests of Octrees[END_REF].

In uniform Cartesian grids with periodic boundary conditions the Fast Fourier Transform (FFT) with its complexity in O(N pt log(N pt )) -N pt denotes the number of points-and its spectral accuracy beats all concurrent numerical methods. In the 80's, the multigrid methods [START_REF] Hackbusch | Multi-grid methods and applications[END_REF] with their complexity in O(N pt log(N pt )) (the log(N pt ) factor stands for the number of iterations necessary to reach the accuracy corresponding to the increase of the number of points N pt ) opened the door to efficient numerical methods suited to non periodic boundaries and immersed boundaries. Their principle (to separate scales to apply Gauss Seidel iterations) inspired the preconditioning of powerful Linear Solvers establishing the algebraic multigrid methods. These are blind to the underlying grid structure and can be used in the adaptive grids as do [START_REF] Kronbichler | High accuracy mantle convection simulation through modern numerical methods[END_REF][START_REF] Duarte | A numerical strategy to discretize and solve the Poisson equation on dynamically adapted multiresolution grids for time-dependent streamer discharge simulations[END_REF] for instance. In the 90's, the Fast Multipole Method [START_REF] Mckenney | A Fast Poisson Solver for Complex Geometries[END_REF][START_REF] Greengard | A Direct Adaptive Poisson Solver of Arbitrary Order Accuracy[END_REF][START_REF] Beatson | A short course on fast multipole methods[END_REF] based on the integral solution of the Poisson Equation and on the properties of its Green kernel, appeared as a concurrent method efficiently handling the adaptive context and the presence of boundaries. This research area continues to develop actively nowadays [START_REF] Hejlesen | A high order solver for the unbounded Poisson equation[END_REF][START_REF] Askham | An adaptive fast multipole accelerated Poisson solver for complex geometries[END_REF].

In the following, we present and investigate an original geometric multigrid method in the AMR framework to solve the Poisson Equation in any dimension at a high order of accuracy. This multigrid algorithm is implemented inside a fully-threaded tree structure programmed in C, but it should also work in block-structured AMR implementations like [START_REF] Barad | A fourth order accurate adaptive mesh refinement method for Poisson's equation[END_REF][START_REF] Miniati | Block structured adaptive mesh and time refinement for hybrid, hyperbolic+N-body systems[END_REF]. Compared to other state-of-the-art Poisson solvers [START_REF] Gholami | Multigrid? A comparative Study of State-Of-the-Art Poisson Solvers for Uniform and Nonuniform Grids in the Unit Cube[END_REF], geometric multigrid Poisson solvers yield very good results regarding the efficiency and the accuracy, similar to Fast Multipole Methods and better than algebraic multigrid methods. The method presented here relies on finite-difference compact schemes introduced by Collatz in [START_REF] Collatz | The numerical treatment of differential equations[END_REF] and more precisely on the original high-order schemes presented in [START_REF] Deriaz | Compact finite difference schemes of arbitrary order for the Poisson equation in arbitrary dimensions[END_REF]. It generalizes to any order and any dimension a fourth-order Poisson solver used to run a two-dimensional magnetic reconnection problem [START_REF] Sarto | A multigrid AMR algorithm for the study of magnetic reconnection[END_REF].

The course of the present paper goes as follows: Section 2 underpins the difference between vertex-based and cell-based AMRs, introduces the use of flags to attribute roles to the points, details the interpolation process and the refinement criterion; Section 3 recalls the original high-order compact finite-difference Poisson schemes from [START_REF] Deriaz | Compact finite difference schemes of arbitrary order for the Poisson equation in arbitrary dimensions[END_REF], unfolds the algorithm of the AMR Poisson solver and studies its stopping criterion, its cost and its convergence; Section 4 is concerned with extended numerical experiments, it tests the convergence of the AMR Poisson solver, its accuracy, its computational efficiency, its OpenMP parallelization, the gain obtained from initial guess and finally it statistically validates the theoretical refinement criterion.

Vertex-based AMR

Most of the AMR schemes rely on finite volume discretization. It is the most natural way to subdivide a computation cell [START_REF] Duarte | A numerical strategy to discretize and solve the Poisson equation on dynamically adapted multiresolution grids for time-dependent streamer discharge simulations[END_REF]. It genuinely conserves the volume and allows to use all the powerful finite volume tools: fluxes, slope limiters, ENO/WENO adaptive stencils.

On the other side, finite difference based AMRs reach high order accuracies more easily [START_REF] Sarto | A multigrid AMR algorithm for the study of magnetic reconnection[END_REF]. It necessitates more complexity in the implementation but reward comes with its compactness and ability to reach high dimensions [START_REF] Deriaz | Six-Dimensional Adaptive Simulation of the Vlasov Equations Using a Hierarchical Basis[END_REF].

Differences between cell-based and vertex-based AMRs

In a d-dimensional space, a cell-based AMR refines a cell at level j by cutting it into 2 d finer cells at level j + 1. Given a cell C j,k ⊂ R d : The two-dimensional case is illustrated by Fig. 2.1.

C j,k = d ∏ i=1 2 -j k i , 2 -j (k i + 1) (2.1) it splits into ε∈{0,1} d C j+1,2k+ε = ε∈{0,1} d d ∏ i=1 2 -j-1 (2k i + ε i ), 2 -j-1 (2k i + ε i + 1) . (2.2) 2k 1 level j + 1 2k 1 + 1 k 1 + 1 level j 2k 2 + 2 2k 2 + 1 2k 1 + 2 k 2 + 1 k 1 k 2 2k 2
A vertex-based AMR refines a point at level j by introducing the 2 d -1 middle points located at its right (by convention) at level j + 1:

2 -j k = 2 -j k i 1≤i≤d spawns the points ε∈{0,1} d 2 -j-1 (2k + ε) . (2.3)
One point at level j becomes 2 d points at level j + 1 and one of these is exactly at the same place, when ε = 0. Hence this vertex-centered structure eases the interpolation process. Fig. 2.2 compares this spawning process with the cell-based splitting.

2k 1 level j + 1 2k 1 + 1 level j 2k 2 + 1 k 1 k 2 2k 2

From interpolability criterion to activation flag

A first difficulty is that the left/right symmetry is broken. The cell-based splitting naturally conserves the left/right symmetry, while the vertex-based spawning breaks it.

A solution to this problem is to symmetrize the new mesh by activating only the middle new vertices which have their two neighboring older vertices. Hence it is necessary to put flags indicating whether a vertex is active or not.

This consideration generalizes to any d-dimensional refinement: at level j, given ε ∈ {0, 1} d the vertex 2 -j (2k+ε) is activated if and only if the vertices 2 -j (2k+2ε ) such that ε ∈ {0, 1} d with ε i ≤ ε i ∀i ∈ [1, d], are activated.

In an approximation perspective, at level j, the older vertices upload their corresponding values at level j -1. Then the new vertices are interpolated from their direct neighbors with a 2nd-order accuracy.

Conversely this 2nd-order interpolation process can be seen as an activation criterion on the new points. A vertex is activated if and only if it can be interpolated through this process.

A similar criterion with 2pth-order interpolability permits to generalize this 2ndorder interpolability criterion and to increase the order of accuracy of the interpolation in the AMR to any 2pth-order with p ≥ 1.

Remark 2.1 In the following numerical experiments, we will need as high as 12thorder interpolation to test the 10th-order accurate AMR three-dimensional Poisson solver.

Vertex typing

Some of the interpolable vertices are located outside any continuous active domain. For instance the fourth order two-dimensional case of Fig. 2.3 allows us to distinguish three kinds of vertices: -the vertices inside the computational domain, -the vertices remaining to the boundary, -and the vertices which are only necessary for the interpolation. Typing these vertices gives us the opportunity to do the right operations at the right places. This typing relies on the fact that for all k node of the level j + 1, the least interpolable vertex 2 -j-1 (2k i + 1) i is active if and only if all its neighbors are active. As a result if it is active (interpolable) then it is located inside the continuous domain of interpolable vertices.

Remark 2.2

The boundary can be made thicker with two points large boundaries, depending on the stencil used for the Poisson solver.

In the C program, all these flags are stored in unsigned long long variables (i.e. on 64 = 2 6 bits) inside nodes, allowing the dimensionality to reach d = 6.

Interpolet

Wavelet theory provides the most comprehensive frame to translate the interpolation process into the basis element of a multiresolution analysis [START_REF] Deslauriers | Symmetric iterative interpolation processes[END_REF][START_REF] Kahane | Fourier Series and Wavelets[END_REF]. We look for an interpolatory compactly supported function x → ϕ p (x) verifying:

ϕ p (0) = 1, ϕ p (k) = 0 ∀k ∈ Z *
that induces a 2pth-order interpolation of the points in 1 2 Z by the points in Z i.e.

Ih 2 ( f ) = ∑ k∈Z ϕ p -k + 1 2 f (kh) = f h 2 + O(h 2p ) ∀ f ∈ C 2p .
Then the function ϕ p satisfies the self-similarity relation:

ϕ p (x) = ϕ p (2x) + ∑ ∈Z c (p) 2 +1 (ϕ p (2x -(2 + 1)) with c (p)
2 +1 the coefficients associated to a 2p-th order interpolation at the point 0 by the points 1 + 2Z. The coefficients associated to the shortest compact supports [-2p + 1, 2p -1] are displayed in table 2.1. And the corresponding basis functions ϕ p (x) are plotted in figure 2.4. These functions are called Interpolets [START_REF] Deslauriers | Symmetric iterative interpolation processes[END_REF].

According to [START_REF] Kahane | Fourier Series and Wavelets[END_REF][START_REF] Deriaz | Six-Dimensional Adaptive Simulation of the Vlasov Equations Using a Hierarchical Basis[END_REF] the coefficients c

(p) 2 +1 are symmetric c (p) -2 -1 = c (p)
2 +1 , non zero only for < p and given by c (p)

2 +1 = p 2 4p-3 2p -1 p -1 (-1) 2 + 1 2p -1 p + . (2.4) For fixed, c (p)
2 +1 is monotonous in p and tends to 2 π

(-1)
2 +1 which is the filter associated to the sinc function also called the Shannon wavelet [START_REF] Kahane | Fourier Series and Wavelets[END_REF][START_REF] Deriaz | Orthogonal Helmholtz decomposition in arbitrary dimension using divergence-free and curl-free wavelets[END_REF], hence

∀x ∈ R lim p→∞ ϕ p (x) = sinc(x) = sin(πx) πx . (2.5)
The sinc function is associated to the Fourier interpolation. 2 +1 for 2nd to 8th orders. In dimension d, the multidimensional interpolet is obtained by tensor product:

2 + 1 -7 -5 -3 -1 1 3 5 7 c (1) 2 +1 1 2 1 2 c (2) 2 +1 -1 16
Φ(x) = d ∏ i=1 ϕ p (x i ).
In the following, we denote by I 2p the 2pth-order interpolation at mid-point using the coefficients c

(p) 2 +1 from Table 2.1.

Refinement criterion

In order to single out the areas of interest which need grid refinement to approximate a given function f , we expand f into a wavelet basis or a hierarchical basis (ϕ λ ) λ with λ = ( j, k):

f (x) = ∑ j≥0 ∑ k∈[0,2 j -1] d d λ ( f )ϕ(2 j x -k) = ∑ λ d λ ( f )ϕ λ (x). (2.6)
If |d λ | > ε j a threshold which depends on the level j, then the node N λ containing ϕ λ is tagged as area of interest. This area extends horizontally to the neighbors of the node N λ , vertically upward to all their ancestors and vertically downward to all its children.

A particular care should be given to the choice of the refinement criterion ε j since it conditions the efficiency of the Adaptive Refinement. If the numerical scheme consistency i.e. the numerical accuracy behaves like O h p f (p) then a best expansion basis (ϕ λ ) should satisfy d λ = O h p f (p) for any smooth f function. The norm equivalence [START_REF] Cohen | Maximal spaces with given rate of convergence for thresholding algorithms[END_REF]:

f B s q (L q ) = O   (2 s j 2 -d j q d λ ) k∈[0,2 j -1] d q j≥0 q   (2.7)
points out the exact criterion minimizing the L q -norm of the sth-derivative of f :

|d λ | ≥ 2 -s+ d q j ε 0 (2.8)
with ε 0 > 0 a constant in agreement with the desired number of discretization points. For instance, in our case we want to solve the Poisson equation ∆ u = v and to minimize the error on the result u controlling the AMR thanks to the wavelet expansion of v = ∑ λ d λ ϕ λ . As we can see u as a s = -2nd derivative of v we should take

ε j = 2 2+ d q j
ε 0 to minimize the L q error on u.

If we just have a wavelet expansion

v = ∑ λ d λ ϕ λ such that d λ = O h p f (p ) with p = p we can consider in first approximation that d λ ( f (p-p ) ) = O 2 (p-p ) j d λ ( f ) . Hence as d λ ( f (p-p ) ) = O h p f (p) , the criterion 2 (p-p ) j |d λ ( f )| ≥ 2 -s+ d q j ε 0 i.e. |d λ | ≥ 2 -s+p -p+ d q j
ε 0 yields a good result if the function f does not have too many disparities regarding its local scales, i.e. f (k) (x) ∼ ξ k (x) f (x) with the local frequency ξ (x) bounded on the whole domain. Otherwise, the criterion should just stick to the solid estimate (2.7) and follow (2.8) or something in between.

The coarse-graining induced by a high-order large compact-support wavelet ϕ may constitute a serious flaw and make a lower-order smaller compact support wavelet expansion more optimal in practice.

Summary

To summarize this vertex-based AMR construction, it goes along the following the steps:

1. Spot the areas of interest thanks to a wavelet expansion. 2. Allocate nodes (with 2 d vertices each) in the tree structure so that a minimal active area encompasses the areas of interest. 3. Given these possibly active vertices and depending on the interpolations needed further, activate the vertices. 4. Attribute types to the vertices: inside the computational domain, on the boundary or just necessary for the interpolation process.

Remark 2.3 All the typed vertices are active. The others are deactivated and may be removed from the memory as it may represent a significant load of memory in the high-dimensional cases.

3 High order AMR Poisson solver

In the following, we present in detail a generalization to any dimension and any order of the two-dimensional fourth-order AMR Poisson solver introduced in [START_REF] Sarto | A multigrid AMR algorithm for the study of magnetic reconnection[END_REF] for the adaptive simulation of the magnetic reconnection. This algorithm necessitates three main ingredients:

vertex-based AMR presented in the previous section, multigrid methods [START_REF] Hackbusch | Multi-grid methods and applications[END_REF][START_REF] Schaffer | Higher Order Multi-Grid Methods[END_REF], high-order compact finite difference schemes for Poisson equations [START_REF] Collatz | The numerical treatment of differential equations[END_REF][START_REF] Spotz | High-order compact finite difference methods[END_REF][START_REF] Iserles | A First Course in the Numerical Analysis of Differential Equations[END_REF][START_REF] Deriaz | Compact finite difference schemes of arbitrary order for the Poisson equation in arbitrary dimensions[END_REF].

Compact finite difference schemes for Poisson equation

These elegant methods were introduced by [START_REF] Albrecht | Taylors-Entwicklungen und finite Ausdrücke für ∆ u und ∆ ∆ u[END_REF][START_REF] Collatz | The numerical treatment of differential equations[END_REF]. The reference [START_REF] Deriaz | Compact finite difference schemes of arbitrary order for the Poisson equation in arbitrary dimensions[END_REF] describes a method to construct arbitrary order compact finite difference schemes for Poisson equation in arbitrary dimension. In particular it derives sixth-order any-dimensional schemes and a tenth-order three-dimensional scheme satisfying a Discrete Maximum Principle condition which implies the stability of these schemes and the convergence of the associated Gauss-Seidel iterations used in the multigrid methods. We call stencil a discrete function A :

Z d → R d , k → a k .
From this function we form an approximation A of the Laplacian ∆ . For any function u sufficiently regular at 0,

lim h→0 1 h 2 ∑ k∈Z d a k u(kh) = ∆ u(0) = d ∑ i=1 ∂ 2 i u(0).
By extension the A stencil may also denote the discrete operator

A A : R d (Z d ) → R d (Z d ) , (u k ) k∈Z d → 1 h 2 ∑ k ∈Z d a k u k+k k∈Z d .
To define a compact finite difference scheme we need two stencils A and B with A an approximation of the Laplace operator and B an approximation of the Identity. The compact scheme is 2pth-order if for all discretized function u h = (u(kh)) k with u C 2p+2 -regular, 1

h 2 Au h = B(∆ u) h + O(h 2p ). (3.1)
If a 0 < 0 and a k ≥ 0 for all k = 0, then the compact scheme (A, B) verifies the Discrete Maximum Principle. This property helps constructing accurate and stable schemes.

Relying on the one dimensional case, it is possible to show that for any dimension there exist such schemes at arbitrary orders. With the growing order of the compact scheme, the radii r A and r B of the stencils A and B, defined as r A = max a k =0 k ∞ become larger. Nevertheless these radii remain much smaller than in the non compact finite difference case. For instance, the 6th-order compact scheme has radius r A = 1 while we have r A = 3 for the non compact scheme.

Multigrid Poisson algorithm

Let us consider embedded continuous subsets of [0, 1] d

Ω J max ⊂ • • • ⊂ Ω j+1 ⊂ Ω j ⊂ • • • ⊂ Ω 0 = [0, 1] d
with periodic boundary condition. The set Ω j defines a discrete set of grid points

ω j = Ω j ∩ 2 -j Z d .
We denote ω c j the complementary of

ω j inside [0, 1] d ∩ 2 -j Z d .
Then, inside ω j we distinguish the interior ω j of ω j and its border ∂ ω j such that dist( ω j , ω c j ) = max(r A , r B ) and ∂ ω j = ω j \ ω j .

Doing so, we obtain boundaries ∂ ω j thick enough to provide for the boundary conditions of the compact scheme (3.1) on each level j. We remind that at level j we have the space step h = 2 -j . We also stress that in the algorithm, from the computer memory point of view,

x λ = x λ if λ = λ even if 2 -j k = 2 -j k .
With all these elements established, the AMR multigrid algorithm at arbitrary order and dimension presented and tested in this manuscript to solve ∆ u = v works as follows:

1. The initial arguments given to the routine are:

(a) a typed adaptive grid (ω j ) j∈[0,J max ] , (b) the right hand side discrete function v, (c) an initial guess u 0 , taken equal to zero if there is no prior estimate of the solution u (d) and a relative error tolerance tol > 0.

2. On the initialization step we make:

(a) u = u 0 , (b) v = vv , so that v = 0, (c) interpolation of the boundaries: for j = 1 to J max , for

x λ ∈ ∂ ω j do v(x λ ) = I 2p (ω j-1 )v u(x λ ) = I 2p+2 (ω j-1 )u.
(d) We form the right hand member Bv:

for j = 1 to J max , for x λ ∈ ω j do f (x λ ) = B v(x λ ) = ∑ k b k v(x j,k+k ).
(e) We start with the residue res = tol × v ∞ .

Beginning of the loop

1 while res ∞ ≥ tol × v ∞ do (a) for j = J max to 1 do i. for x λ ∈ ω j do res(x λ ) = f (x λ ) - 1 h 2 A u(x λ ) tmp(x λ ) = h 2 a 0 res(x λ ).
ii. Execute twice the Gauss-Seidel iteration:

for x λ ∈ ω j do tmp(x λ ) = tmp(x λ ) + h 2 a 0 res(x λ ) - 1 h 2 A tmp(x λ ) .
iii. Update res and u at the level j for x λ ∈ ω j do

res(x λ ) = res(x λ ) - 1 h 2 A tmp(x λ ) u(x λ ) = u(x λ ) + tmp(x λ ).
iv. Pass the residual and u to the level j -1 with the tensorial stencil 1 4 1 2

1 4 d for x j-1,k ∈ ω j-1 such that x j,2k ∈ ω j do res(x j-1,k ) = 1 4 1 2 1 4 d res(x j,2k ) u(x j-1,k ) = 1 4 1 2 1 4 d u(x j,2k
).

(b) S 0 = res(x 0,0 ) for j = 0 to J max , for

x λ ∈ ω j do f (x λ ) = f (x λ ) -S 0 .
(c) Prediction-correction step for j = 1 to J max , for

x λ ∈ ω j do u(x λ ) = u(x λ ) + I 2p+2 (ω j-1 )tmp(x λ ) tmp(x λ ) = tmp(x λ ) + I 2p+2 (ω j-1 )tmp(x λ ).
4. At the end of the iterative process, the result is available in u.

Computational cost

As one can see, this special case of geometric multigrid method naturally extends the original multigrid algorithm [START_REF] Hackbusch | Multi-grid methods and applications[END_REF] to a vertex-based AMR. In particular considering points instead of volume cells as most AMR papers do [START_REF] Duarte | A numerical strategy to discretize and solve the Poisson equation on dynamically adapted multiresolution grids for time-dependent streamer discharge simulations[END_REF][START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF][START_REF] Teyssier | Cosmological hydrodynamics with adaptive mesh refinement. A new high resolution code called RAMSES[END_REF] makes it easier to stay close to the original finite-difference multigrid design. The computational costs are similar to the uniform grid case although it includes the overhead induced by the more complex AMR data structure, cf Fig. 4.6. The observed minimum tolerance after which the residual stops decreasing lies in between ε and 10 4 ε where ε = 10 -16 stands for the computer error in double precision, cf Fig. 4.5, depending on the accuracy of the approximation (i.e. ∼ h min ). The greater is the accuracy, the larger will be the minimum bound since a good approximation of it is given by

2d u ∞ v ∞ h -2 min ε or more precisely sup x∈[0,1] d 2d |u(x)| v ∞ h(x) -2 ε .
Depending on the number of dimensions and on the compact scheme used, this tolerance is reached in about fifteen iterations. The most costly operations are the application of A four times for each iteration, once for computing the residual (step 3(a)i.), twice in the Gauss-Seidel method (step 3(a)ii.) and once when updating the residual (step 3(a)iii.), and the computation of the right hand member f i.e. applying B once (step 2(d)). The cost of applying a stencil depends on the number of its non-zero elements. Hence the total cost C of the Poisson solver depends on the cardinality of the stencil A mainly and of the stencil B in a lesser extent and can be written:

C = (#B + N it (4#A + 2 × 3 + 2(2p + 2))) N pt (3.2)
where #B and #A stand for the numbers of non zero elements of the stencils B and A, N it the number of iterations (of the loop 3.), N pt the number of points and 2p the order of the method. We neglected the small additional costs dependent on the size of the boundaries including them in the fuzzy definition of N pt . The cardinalities #A and #B directly depend on the order p and on the dimensionality d. They should roughly evolve as O(p d-1 ) for #A and O(p d ) for #B although in the non compact case they have the values #A = 2d p + 1 and #B = 1.

For instance, the 3-dimensional 6th-order HOC stencil [START_REF] Spotz | High-order compact finite difference methods[END_REF] verifies #A = 27 (or 16 if we take advantage of the tensorial structure of the stencil), #B = 25 and only needs 12 iterations to converge to the computer rounding error so its cost is given by C = 1585N pt . Of course it is possible to decrease this cost by fixing a larger error tolerance and taking fewer iterations. For instance converging to 10 -6 only takes 5 iterations then the cost is given by C = 675N pt .

This compares advantageously to the costs given for the pth-order Fast Multipole Method [START_REF] Beatson | A short course on fast multipole methods[END_REF] whose optimal implementation in 3 dimensions yields

C = 200N pt p + 3.5N pt p 2 .
Which means C = 1326N pt for the 6th-order p = 6 case. This small computation confirms the findings from [START_REF] Gholami | Multigrid? A comparative Study of State-Of-the-Art Poisson Solvers for Uniform and Nonuniform Grids in the Unit Cube[END_REF] which states that the two methods have equivalent costs.

Convergence of the multigrid iterations

As this AMR multigrid algorithm works very similarly to the original one [START_REF] Hackbusch | Multi-grid methods and applications[END_REF], it should be possible to prove its convergence similarly. If we adopted the technique presented in [START_REF] Teyssier | Cosmological hydrodynamics with adaptive mesh refinement. A new high resolution code called RAMSES[END_REF] which solves the Poisson equation recursively on the ω j domains for j decreasing then the proof of convergence is exactly the same as in the uniform case. Here there is an additional correction from the finer levels to the coarser ones because the residual is locally computed on the finest level. This influences the rate of convergence.

For instance if r A = 2 and the first boundary goes through the pivotal points (2k i ) i which straightly get their values from level j, and not through the points in the middle of the cells (2k i + 1) i which are the most dependent on interpolation (see Fig. 3.1 to visualize the situation) then the convergence rate shifts from 1/2 to 1/8 which represents a non-negligible gain of efficiency.

In order to reduce the number of iterations N it in (3.2), the initial guess u 0 chosen as close as possible to u in the H 2 -norm allows a dramatic cut of the computational costs. In the simulations with time evolution it is possible to first extrapolate the solution u from the previous time steps in order to provide an accurate initial guess.

The H 2 -norm criterion comes from the computation of the residual as res = v -∆ u 0 which is the pivotal operation of the Gauss Seidel iterative algorithm. Hence a noisy initial guess, even close to the solution in the L ∞ -norm, has little interest.

This specificity of the geometric multigrid methods is a strong incitement to use them in simulations. Let us notice that the algebraic multigrid methods can benefit from similar initial guess when they rely on iterative solvers [START_REF] Kronbichler | High accuracy mantle convection simulation through modern numerical methods[END_REF].

Numerical experiments

High order compact schemes (A,B) tested

The finite difference compact schemes used to test the present high order AMR Poisson solver were obtained in [START_REF] Deriaz | Compact finite difference schemes of arbitrary order for the Poisson equation in arbitrary dimensions[END_REF]. All of them satisfy a condition implying the Discrete Maximum Principle and favorable to the fast convergence of the Gauss-Seidel iterations. This condition states: a 0 < 0 as small as possible and all the other a k coefficients positive if k = 0.

Tenth-order three-dimensional scheme

This scheme was designed to be tenth-order and to have few points in the stencil A although it does not realizes the minimum.

The A stencil is composed of 59 points: the central point 0, six symmetric points of type (1, 0, 0), twelve of type (1, 1, 0), eight of type (1, 1, 1), 24 of type (2, 1, 1) and eight of type (2, 2, 2), so it has a radius r A = 2. The B stencil has 113 points and a radius r B = 4. All their localizations and their values are indicated in [START_REF] Deriaz | Compact finite difference schemes of arbitrary order for the Poisson equation in arbitrary dimensions[END_REF].

ω j-1 ∂ ω j ω j ω j-1
∂ ω j ω j Fig. 3.1 Two ways of fixing the boundary with respect to the interior domain ω j : setting the limit just after the middle points (on the left) or just after the pivotal points (on the right). The first case is more efficient with respect to memory management when r A is even, the second case converges much faster.

Sixth-order any-dimensional schemes

This group of schemes generalizes the HOC (High Order Compact) stencil [START_REF] Albrecht | Taylors-Entwicklungen und finite Ausdrücke für ∆ u und ∆ ∆ u[END_REF][START_REF] Spotz | High-order compact finite difference methods[END_REF] to any dimensions. Their stencils have tensorial structures so their complexities only increase linearly with the number of dimensions d.

For α, p, β , γ real numbers, the A stencil is of the form:

A = α[p 1 p] d + γ | γ -β -γ | γ (4.1)
namely

a 0 = α + β , a (1,0,...,0) = α p + γ, and a (1, . . . , 1 k times ,0,...,0) = α p k for k ≥ 2.
To reach the sixth-order accuracy, computations presented in [START_REF] Deriaz | Compact finite difference schemes of arbitrary order for the Poisson equation in arbitrary dimensions[END_REF] lead to:

p = 1 3 , α = 3 2 3 5 d-2 , β = - 25 + 2d 6 and γ = 1 6 .
The B stencil is built similarly: The radii of A and B are r A = 1 and r B = 2. Applying the stencil A has a 1 + 5d complexity instead of the usual O(d 3 ) minimum complexity of compact none-tensorial type stencils, and of the 1 + 6d complexity of the non-compact stencil.

B = ω[q 1 q] d + b 01 | 0 | b 01 -0 -λ -0 -b 01 | 0 | b 01 namely b 0 = ω + λ , b (1, . . . , 1 k times ,0,...,0) = ω q k for k ≥ 1,

Gaussian test case

In the periodic domains [-0.5, 0.5] d we consider the Gaussian function:

u(x) = F(r) = exp - r 2 σ 2 + b (4.2)
with r = |x| and the parameter σ small enough so that exp(-0.5 2 /σ 2 ) ε the computer rounding error. For lower dimensions we take σ = 0.005 as in [START_REF] Duarte | A numerical strategy to discretize and solve the Poisson equation on dynamically adapted multiresolution grids for time-dependent streamer discharge simulations[END_REF] so exp(-0.5 2 /σ 2 ) = exp(-100 2 ) ≈ 10 -4343 . The number b is chosen such that

[-0.5,0.5] d u(x) dx = 0.
With periodic boundary conditions, the solution u of the equation ∆ u = v is determined modulo an additive constant. The non uniformity of the grid makes the condition u(x) dx = 0 difficult to reach in practice. Hence we choose a b minimizing the L ∞ error from the discrete solution to the exact one. The argument function v, right hand member of the equation ∆ u = v, is given by:

v(x) = 4|x| 2 σ 4 - 2d σ 2 exp - |x| 2 σ 2 + c
with c a free real constant such that v(x) dx = 0 in the discrete form.

Verification of the order of convergence and of the accuracy

To assert the accuracies of the schemes, we begin with a test where starting from a coarse non uniform grid we half its mesh length everywhere once, then twice, and so on, computing the resulting L ∞ , L 2 and L 1 errors. That way the orders of convergence appear clearly in two to six dimensions. All the following numerical experiments were run on a 16 cores, 32 threads, 32GB RAM Dell Precision T5610 computer. 4.4 show the convergence tests of the 6th-order and 10th-order schemes respectively. These can be compared to the same experiment in [START_REF] Duarte | A numerical strategy to discretize and solve the Poisson equation on dynamically adapted multiresolution grids for time-dependent streamer discharge simulations[END_REF] where a 2nd-order algebraic multigrid Poisson solver induces a 10 -6 L 2 -error in the finest 3D case. Given the fact that the L 2 -norm of the exact solution is approximately 5 × 10 -4 , one can measure the improvement of accuracy that the present geometric multigrid method provides.

Computational efficiency

The cost estimate Eq. (3.2) is verified quite accurately knowing that each operation (rawly speaking accessing an element of the grid) costs 10 -8 s. The OpenMP parallelization of the C code proves useful up to eight threads.

In Fig. 4.5 we visualize the convergence of the iterative process and see the number of necessary loop iterations in the Multigrid Poisson algorithm, paragraph 3.2. For instance, no more than 10 iterations are necessary to reach a residual equal 10 -11 ∆ u ∞ criterion ending the iterative process, with the 6th-order 2D scheme.

As expected theoretically due to the tensorial structure of the stencil, and as shown in Fig. 4.6, the cost depends bilinearly on the number of points and on the number of dimensions. Above four dimensions, the shape of the adaptive grid plays a prominent role as illustrated by the two six-dimensional cases: in the most costly one, the interior points only represent 25% of the total number of points, while in the other one they represent 95% of it. Regarding the 10th-order 3-dimensional scheme, its computer times amount only to one and a half those of the 6th-order scheme.

For this experiment Fig. 4.6, the minimum number of interior points was taken equal to 15e+3 in 2D, 3e+5 in 3D, 1.36e+6 in 4D, 2.78e+6 in 5D and 17e+6 in 6D. These figures are multiplied by a number approximately equal to 2 d to plot the second curve, except for the 5D and 6D cases where the factors were taken equal to 10 and 3.3 respectively.

From now on we skip the 6-dimensional experiments because of their enormous memory needs.

Based on three runs of the Poisson algorithm with the residual inferior to 10 -10 ∆ u ∞ ending criterion, Fig. Compared to p4est [START_REF] Isaac | Recursive Algorithms for Distributed Forests of Octrees[END_REF], these relatively poor parallelization results show that there is room for improvement.

Initial guess

Starting from an initial guess u G close to solution u, in the sense ∆ u G -∆ u small compared to ∆ u , allows to decrease the number of iterations necessary to reach convergence. To test this concept and prove its efficiency, we take

u G (x) = u(x) + ε cos(ξ x 1 ) so ∆ u G (x) = ∆ u(x) -ξ 2 ε cos(ξ x 1 ), (4.3)
and we look at the number of iterations needed to reach a residual equal to 10 -10 ∆ u ∞ . Starting without any initial guess i.e. u G = 0, the algorithm needs exactly 10 iterations with the 6th-order and 10th-order stencils to reach this ending condition. In order to discard the approximation error, we replace in the definition of u G Eq. (4.3) the exact solution u by the converged numerical solution of the discrete problem.

For the experiment Fig. 4.8 we used the 3-dimensional Gaussian with σ = 0.005 as u function, and for both cases, the sixth-order and the tenth-order schemes, the same adaptive grid with 1, 139, 121 interior points distributed on scales going from 32 3 to 4096 3 uniform grids equivalent. The computational time comes out approximately the same (approximately 6.5s for 10 iterations with 8 threads) for both cases, in accordance with Fig. 4.7. The sixth-order scheme ends up with a 3.2 × 10 -9 L ∞error and a 1.1 × 10 -9 L 2 -error while the tenth-order one ends up with a 9 × 10 -12 L ∞ -error and a 3.7 × 10 -13 L 2 -error.

Oscillating Gaussian test case

The Gaussian test case contains a limited oscillating behavior which does not allow to explore the possibilities of the adaptive refinement. Hence we introduce the Oscillating Gaussian test case.

In dimension 2 ≤ d ≤ 6 we consider the following radial function:

F ε (r) = e -A r 2 cos B r 2 + ε α (4.4)
with A = 150, ε > 0 a parameter which creates oscillations at an approximative distance r = √ ε from 0. For N ∈ N * a given number we take α = -log 4N log 4ε and B = π 1/α 4•2 1/α so that the function F ε (r) = e -A r 2 cos (X(r))

with X(0) = 2Nπ and X( 1 2 ) = π 2 (1 + 4ε) α ∼ π 2 contains N oscillations from the center of the domain to its boundary.

Function analysis

Function r → F ε (r) is C ∞ on (-0.5, 0.5) and it is sufficiently small at 0.5 (∼ 10 -20 ) so its C 1 discontinuity at the periodic boundary remains numerically benign. Suppose u

(x) = F ε (r) with r 2 = ∑ d i=1 x 2 i then, if F ε (r) = G(r 2 ), v(x) = ∆ u(x) = 2 d G (r 2 ) + 4 r 2 G (r 2 ).
This provides the right hand member for the numerical experiments.

We infer the maximal frequency by studying the variations of the function:

g(r) = B r 2 + ε α . (4.5)
We obtain g (r) = -2αrB α (r 2 + ε) -α-1 maximal for r = r M = ε 2 α+1 . At this distance to 0, g (r M ) ∼ 2πNα √ ε , the maximal frequency reached by u, provides the smallest scale for the numerical scheme.

Validation of the refinement criterion

The diversity of behaviors displayed by this test function enables to test the refinement criterion. Let us challenge the non-linear approximation theoretical predictions and try to catch the best criteria.

As shown in Fig. 4.3 and 4.4, the present Poisson solver has orders p = 6 or p = 10 so the L q -error ε L q = O(h -p ). For an adaptive grid in dimension d the local mesh size h(x) = 2 -j(x) varies as a discontinuous function of the total number of points N pt . It is more or less given by j(x) = log K v,α (x)N In the following we look for an optimal α minimizing the compacity factor K v,α i.e. driving to the most efficient non-linear approximation. To assert the study presented in part 2.5, we ran a large number of two-dimensional instances of the sixthorder adaptive Poisson solver with a criterion based on the sixth-order interpolation of the right-hand member v. In Fig. 4.11 we plot the total number of optimal cases when α = k is optimal for various ε values in F, Eq. (4.4), and various numbers of points N pt . It more or less looks like a Gaussian distribution centered on α = 2 which is the theoretical value obtained in part 2.5. 4.4) and the number of points N pt we retained the optimal α providing the smallest value for K v,α in Eq. (4.6). The theoretically optimal value is α = 2.

Conclusion

The present AMR geometric multigrid Poisson solver proves computationally efficient compared to state-of-the-art modern Poisson solvers and fitted for any order of accuracy and any dimension.

It relies on vertex-centered Cartesian non-uniform grids, although there is no reason why it should not work with cell-based AMR. With cell-based AMR, its rate of convergence which conditions its efficiency should be expected between 0.125, the best vertex-centered case and 0.5, the worse one, cf Its use of the 10th-order compact scheme developed in [START_REF] Deriaz | Compact finite difference schemes of arbitrary order for the Poisson equation in arbitrary dimensions[END_REF] demonstrates its ability to increase its order of accuracy. Maybe such an increase could be interesting for delicate simulations such as ab initio molecular dynamics [START_REF] Genovese | Efficient and accurate three-dimensional Poisson solver for surface problems[END_REF]. Despite all the experimental efforts, it was possible to demonstrate only a statistical optimal refinement criterion in agreement with the non-linear approximation theory. Thanks to its possibility to start from an initial guess, this geometric solver finds a very efficient use in non stationary equations, such as Navier-Stokes, Schrödinger or Vlasov-Poisson equations. In these time dependent equations an initial guess is provided by temporal extrapolation.

There are quite a few tracks to extend this work:

the 6th-order scheme combined with a 4th-order immersed boundary method [START_REF] Hosseinverdi | An efficient, high-order method for solving Poisson equation for immersed boundaries: Combination of compact difference and multiscale multigrid methods[END_REF] would result in a sixth-order AMR immersed boundary scheme, the production of an any-dimensional high-order MPI-parallelized code sounds feasible in a near future, there is also the possibility to directly compute the gradient of the potential from the density function as it is done in Fast Multipole Methods [START_REF] Askham | An adaptive fast multipole accelerated Poisson solver for complex geometries[END_REF].

Fig. 2 . 1

 21 Fig. 2.1 Splitting of a two-dimensional cell into four cells. The numbers indicate their positions in their respective levels.

Fig. 2 . 2

 22 Fig. 2.2 Two-dimensional spawning of the new three vertices from the initial one. The numbers indicate their numbering in their respective levels.

Fig. 2 . 3

 23 Fig. 2.3 Instance of fourth-order consistent distribution of vertices with their types. Only one among the horizontal couple of interpolatory points or the vertical one is necessary.

Fig. 2 . 4

 24 Fig. 2.4 Interpolets of 2nd, 4th, 6th and 8th orders tending to the sinc function.

  and b (2,0,...,0) = b 01 ,

Fig. 4 . 1

 41 Fig. 4.1 Gaussian with σ = 0.5 6 much larger than σ = 0.005.

Fig. 4 . 2

 42 Fig. 4.2 First two instances of the 3-dimensional grid with maximum refinement J = 9 corresponding to a 512 3 uniform grid (left) and J = 10, a 1024 3 grid (right) used to test the accuracy of the 6th-order scheme. Different colors correspond to different point concentrations.

Fig. 4 . 3

 43 Fig. 4.3 Order of convergence for the 6th-order scheme in various dimensions. The dashed line represents the 6th-order slope.

Figure 4 .

 4 Figure 4.2 shows two successive instances of grid. The grid points are represented by dots whose colors change when the mesh length halves. Figures 4.3and 4.4 show the convergence tests of the 6th-order and 10th-order schemes respectively. These can be compared to the same experiment in[START_REF] Duarte | A numerical strategy to discretize and solve the Poisson equation on dynamically adapted multiresolution grids for time-dependent streamer discharge simulations[END_REF] where a 2nd-order algebraic multigrid Poisson solver induces a 10 -6 L 2 -error in the finest 3D case. Given the fact that the L 2 -norm of the exact solution is approximately 5 × 10 -4 , one can measure the improvement of accuracy that the present geometric multigrid method provides.

  Figure 4.2 shows two successive instances of grid. The grid points are represented by dots whose colors change when the mesh length halves. Figures 4.3and 4.4 show the convergence tests of the 6th-order and 10th-order schemes respectively. These can be compared to the same experiment in[START_REF] Duarte | A numerical strategy to discretize and solve the Poisson equation on dynamically adapted multiresolution grids for time-dependent streamer discharge simulations[END_REF] where a 2nd-order algebraic multigrid Poisson solver induces a 10 -6 L 2 -error in the finest 3D case. Given the fact that the L 2 -norm of the exact solution is approximately 5 × 10 -4 , one can measure the improvement of accuracy that the present geometric multigrid method provides.

Fig. 4 . 4

 44 Fig. 4.4 Convergence for the 3-dimensional 10th-order scheme in various norms. The dashed line represents the 10th-order slope.

  4.7 represents the average time of execution per interior point when applying an OpenMP task parallelization. Increasing the number of threads un-

Fig. 4 . 5

 45 Fig. 4.5 Convergence process of the algorithm 3.2. The first iterations achieve the largest reductions of the residual.

Fig. 4 . 6

 46 Fig. 4.6 Time per interior point to perform ten iterations of the Multigrid Poisson algorithm when varying the number of points and the dimension.

  til 8 to 11 allows to decrease the time of execution. The tenth-order scheme scales better than the tensorial sixth-order scheme because at each Gauss-Seidel iteration this latter needs to visit a node three times with lighter operations when the first one only necessitates one visit with heavier operations. So memory access is more concentrated in the tenth-order case.The observed cost from Fig.4.7 for the tenth order scheme, 3.210 -5 × 10 8 = 3, 200 operations, is in agreement with the theoretical cost C = 3, 305 obtained from Eq. (3.2) with #B = 113, #A = 59, N it = 12 and p = 10.

Fig. 4 . 7

 47 Fig. 4.7 Time per interior point to converge the Multigrid Poisson algorithm to ∆ u n -∆ u < 10 -10 ∆ u when varying the number of threads. This tests the strong scaling. The dashed line represents an ideal scaling.

Fig. 4 . 8

 48 Fig. 4.8 Number of iterations needed to reach a residual equal to 10 -10 ∆ u ∞ ( ∆ u ∞ = 2.4 × 10 5 ) and stop the iterative process, starting from the initial guess Eq. (4.3). We generate a new point by multiplying ξ by 3. At ξ 2 ε = 10 -5 we shift ξ from 2187 × 2π back to 2π and ε from 10 -8 to 4.8 × 10 -2 .

Fig. 4 . 9

 49 Fig.4.9 Oscillating Gaussian with N = 3 oscillations centered at 0, (4ε) 0.418 /2 and (4ε) 0.279 /2.

- 1

 1 /d ptwhere K v,α (x) is a local constant depending on the right-hand-side function v and on the refinement criterion α.On the whole domain and in average we would expect a relation of the kindε L q = K v,α N -p/d pt . (4.6)We verify this relation experimentally in Fig.4.10. The wavelet coefficients d λ represent the local error of the 4th-order interpolation of the right-hand member v so it is unclear which α parameter should be optimal. The experimental results verify the general trend of Eq. (4.6) although it also shows a slightly chaotic behavior.

Fig. 4 .

 4 Fig. 4.10 Experimental L q -error as a function of the number of points for the case ε = 10 -4 and N = 3 in Eq. (4.4). The theoretical slopes ε L q = KN -p/d pt ( p d = 6 2 , 10 3 , 6 3 from left to right) are represented by dashed lines.

Fig. 4 .

 4 Fig. 4.11 Test of the optimality of the refinement criterion α. Varying ε in F ε of Eq. (4.4) and the number of points N pt we retained the optimal α providing the smallest value for K v,α in Eq. (4.6). The theoretically optimal value is α = 2.

Fig. 3 . 1

 31 Part 3.4.

Table 2 . 1

 21 Interpolation coefficients c

	(p)

The algorithm uses two intermediate variables: res which stores the residual and tmp which stores the increments of u provided by the Gauss-Seidel algorithm.

Acknowledgements The author is thankful to Marc Massot and Christian Tenaud for instructive discussions about algebraic and geometric multigrids.