Stéphane Popinet
email: s.popinet@niwa.co.nz

Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries

Keywords: Incompressible flow, adaptive mesh refinement, approximate projection method, complex geometry 1991 MSC: 65M06, 65N06, 65N55, 76D05, 76-04

An adaptive mesh projection method for the time-dependent incompressible Euler equations is presented. The domain is spatially discretised using quad/octrees and a multilevel Poisson solver is used to obtain the pressure. Complex solid boundaries are represented using a volume-of-fluid approach. Second-order convergence in space and time is demonstrated on regular, statically and dynamically refined grids. The quad/octree discretisation proves to be very flexible and allows accurate and efficient tracking of flow features. The source code of the method implementation is freely available.

Introduction

Efficient techniques for the numerical simulation of low Mach number flows have a large range of applications: from fundamental fluid mechanics studies such as turbulence or interfacial flows, to engineering and environmental problems. For time-dependent flows, the finite speed of propagation of sound waves can lead to strong restrictions on the maximum value of the timestep. While filtering techniques can be applied to try to lift this constraint, a better approach is to assume that the fluid considered is strictly incompressible. This introduces an elliptic problem for the pressure which expresses the instantaneous propagation of pressure information throughout the entire domain. In practice, this leads to the fundamental change from a spatially explicit to a spatially implicit problem.

Projection methods and multigrid solvers have proved an efficient combination to solve this type of problem [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Peyret | Computational Methods for Fluid Flow[END_REF][START_REF] Brandt | Guide to multigrid development, Multigrid Methods[END_REF][START_REF] Wesseling | An introduction to Multigrid Methods[END_REF]. More recently, these techniques have been extended through the use of higher-order, unconditionally stable advection schemes [START_REF] Bell | A second-order projection method for the incompressible Navier-Stokes equations[END_REF][START_REF] Bell | A second-order projection method for variable density flows[END_REF].

Another characteristic of fluid flows is the very wide range of spatial scales often encountered: shocks in compressible flows, interfaces between immiscible liquids, turbulence intermittency, boundary layers and vorticity generation near solid boundaries are just a few examples. Consequently, in recent years a number of researchers have investigated the use of adaptive mesh refinement, where the spatial discretisation is adjusted to follow the scale and temporal evolution of flow structures [START_REF] Berger | Adaptive mesh refinement for shock hydrodynamics[END_REF][START_REF] Coirier | An adaptively-refined, Cartesian, cell-based scheme for the Euler and Navier-Stokes equations[END_REF][START_REF] Howell | An adaptive mesh projection method for viscous incompressible flow[END_REF].

For compressible flows, two main approaches have been developed: the hierarchical structured grid approach of Berger and Oliger (Adaptive Mesh Refinement, AMR) [START_REF] Berger | Adaptive mesh refinement for shock hydrodynamics[END_REF] and quad/octree based discretisations [START_REF] Coirier | An adaptively-refined, Cartesian, cell-based scheme for the Euler and Navier-Stokes equations[END_REF][START_REF] Khokhlov | Fully threaded tree algorithms for adaptive refinement fluid dynamics simulations[END_REF]. The AMR framework uses classical algorithms on regular Cartesian grids of different resolutions arranged hierarchically. The only modification necessary is to allow coupling between grids at different levels through the boundary conditions. Quad/octree discretisations, on the other hand, deal with various levels of refinement locally through the use of finite-difference operators adapted to work at fine/coarse cell boundaries.

The AMR framework has been extended to incompressible flows by Minion [START_REF] Minion | A projection method for locally refined grids[END_REF], Almgren et al. [START_REF] Almgren | A numerical method for the incompressible Navier-Stokes equations based on an approximate projection[END_REF] and Howell and Bell [START_REF] Howell | An adaptive mesh projection method for viscous incompressible flow[END_REF] but we are not aware of any quad/octree implementation of adaptive mesh refinement for incompressible flows. The natural hierarchical nature of tree-based discretisations is well suited for multigrid implementations. Moreover, we believe that the flexibility and simplicity of mesh refinement and coarsening of quad/octrees can be a significant advantage when dealing with complex solid boundaries or evolving interfacial flows.

Complex solid boundaries are usually represented using boundary-following structured curvilinear grids or unstructured grids. While boundary conditions can be easily and accurately applied on such grids, grid generation can be a difficult and time consuming process. In recent years, "Cartesian grids" [START_REF] Berger | A rotated difference scheme for Cartesian grids in complex geometries[END_REF][START_REF] Quirk | An alternative to unstructured grids for computing gas dynamics flows around arbitrarily complex two-dimensional bodies[END_REF][START_REF] Almgren | A Cartesian grid projection method for the incompressible Euler equations in complex geometries[END_REF][START_REF] Ye | An accurate cartesian grid method for viscous incompressible flows with complex immersed boundaries[END_REF][START_REF] Calhoun | Solving the advection-diffusion equation in irregular geometries[END_REF] and "immersed boundary" [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF][START_REF] Saiki | Numerical simulations of a cylinder in a uniform flow: application of a virtual boundary method[END_REF][START_REF] Fadlun | Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations[END_REF] techniques have known a regain of interest because they greatly simplify the grid generation process. This flexibility comes at the cost of a more complex treatment of boundary conditions at solid boundaries.

In this light, we present a numerical method for solving the incompressible Euler equations, combining a quad/octree discretisation, a projection method and a multilevel Poisson solver. Advection terms are discretised using the robust second-order upwind scheme of Bell, Colella and Glaz [START_REF] Bell | A second-order projection method for the incompressible Navier-Stokes equations[END_REF] and complex solid boundaries are treated through a Cartesian volume-of-fluid approach. On a uniform grid without solid boundaries, the approach presented reduces to the approximate projection method described by Martin [START_REF] Martin | An adaptive cell-centered projection method for the incompressible Euler equations[END_REF][START_REF] Martin | A cell-centered adaptive projection method for the incompressible euler equations[END_REF]. Solid boundaries are treated using a combination of a Poisson solver similar to the one studied by Johansen and Colella [START_REF] Johansen | A cartesian grid embedded boundary method for poisson's equation on irregular domains[END_REF][START_REF] Day | Embedded boundary algorithms for solving the Poisson equation on complex domains[END_REF] and of a cell-merging technique for the advection scheme [START_REF] Quirk | An alternative to unstructured grids for computing gas dynamics flows around arbitrarily complex two-dimensional bodies[END_REF]. In contrast to classical AMR strategies, adaptive refinement is performed at the fractional timestep. While we restrict this description to two-dimensional flows for clarity, the extension to three dimensions is straightforward: the source code of the three-dimensional parallel implementation [START_REF] Popinet | The Gerris Flow Solver[END_REF] can be freely accessed, redistributed and modified under the terms of the Free Software Foundation General Public License.

Spatial discretisation

The domain is spatially discretised using square (cubic in 3D) finite volumes organised hierarchically as a quadtree (octree in 3D) [START_REF] Samet | Applications of Spatial Data Structures[END_REF]. This type of discretisation has been used and studied extensively for image processing and computer graphics applications [START_REF] Gonzalez | Digital Image Processing[END_REF][START_REF] Samet | Applications of Spatial Data Structures[END_REF] and more recently applied to the solution of the Euler equations for compressible flows [START_REF] Coirier | An adaptively-refined, Cartesian, cell-based scheme for the Euler and Navier-Stokes equations[END_REF][START_REF] Khokhlov | Fully threaded tree algorithms for adaptive refinement fluid dynamics simulations[END_REF]. An example of spatial discretisation and the corresponding tree representation is given in figure 1. In what follows we will refer to each finite volume as a cell. The length of a cell edge is denoted by h. Each cell may be the parent of up to four children (eight in 3D). The root cell is the base of the tree and a leaf cell is a cell without any child. The level of a cell is defined by starting from zero for the root cell and by adding one every time a group of four descendant children is added. Each cell C has a direct neighbour at the same level in each direction d (four in 2D, six in 3D), noted N d . Each of these neighbours is accessed through a face of the cell, noted C d . In order to handle embedded solid boundaries, we also define mixed cells which are cut by a solid boundary.

To simplify the calculations required at the cell boundaries, we add the constraints illustrated in figure 2: (a) the levels of direct neighbouring cells cannot differ by more than one. (b) the levels of diagonally neighbouring cells can not differ by more than one. (c) All the cells directly neighbouring a mixed cell must be at the same level.

While not fundamentally necessary, these constraints greatly simplify the gradient and flux calculations presented in this article. Constraints (a) and (b) have little impact on the flexibility of the discretisation (they only impose gradual refinement by increments of two). Constraint (c) is more restrictive as it forces all the cells cut by the interface to be at the same level (i.e. the whole solid boundary must be described at the same resolution). It is also important to note that a major restriction of the quad/octree structure is that it imposes a locally spatially isotropic refinement. This can be an issue in highly non-isotropic flows (i.e. boundary layers, large scale atmospheric flows etc. . . ). A limited solution is to use a rectangle instead of a square as root cell, thus resulting in a fixed refinement ratio between the corresponding spatial directions. A more general (and complicated) approach would be to use the "variable quadtree" approach of Berger et al. [START_REF] Berger | Aspects (and aspect ratios) of Cartesian mesh methods[END_REF].

In practice, the choice of a data structure to represent the tree is conditioned by the following requirements: At present, we use the fully-threaded tree structure presented by Khokhlov [START_REF] Khokhlov | Fully threaded tree algorithms for adaptive refinement fluid dynamics simulations[END_REF] which allows (a) and (b) to be performed in O(1) operations (versus O(log N ) for a standard pointer-based structure). Operations (c) are performed in O(N log N ) using the standard pointer-based tree description (N is the number of cells traversed). Other modern quad/octree representations might be as good or better (in particular, the linear quadtree encoding of Balmelli et al. [START_REF] Balmelli | Quadtrees for embedded surface visualization: constraints and efficient data structures[END_REF] is noteworthy).

The primitive variables of the Euler equations (velocity U and pressure p) are all defined at the centre of the cells. In mixed cells, the solid boundary is defined through a volume-of-fluid type approach. Specifically, we define:

• the volume fraction a as the ratio of the volume occupied by the fluid to the total volume of the cell, • the surface fraction in direction d, s d as the ratio of the area of face C d occupied by the fluid to the total area of the face.

This solid boundary description assumes that the geometries represented do not possess features with spatial scales smaller than the mesh size. In particular, sharp angles or thin bodies cannot be represented correctly. This can be an issue for some applications, but more importantly, as argued by Day et al. [START_REF] Day | Embedded boundary algorithms for solving the Poisson equation on complex domains[END_REF], it will restrict the efficiency of the multigrid solver.

Computing the volume and area fractions can be expressed in terms of boolean operations (intersection, union, difference) between curves (in 2D) or volumes (in 3D). This is a difficult problem to solve in a robust manner (due to the limited precision of arithmetic operations in computers). Because of their numerous practical applications, robust geometrical operations have attracted considerable attention from the computational geometry community in recent years [START_REF] Edelsbrunner | Simulation of simplicity: A technique to cope with degenerate cases in geometric algorithms[END_REF][START_REF] Clarkson | Safe and effective determinant evaluation[END_REF][START_REF] Milenkovic | Robust polygon modeling[END_REF][START_REF] Shewchuk | Adaptive precision floating-point arithmetic and fast robust geometric predicates[END_REF]. Drawing from these results, we use the boolean operations implemented in the GTS Library [START_REF] Popinet | The GNU Triangulated Surface Library[END_REF] based on an approach similar to that presented by Aftosmis et al. [START_REF] Aftosmis | Robust and efficient Cartesian mesh generation for component-based geometry[END_REF].

Temporal discretisation

We consider a constant density, incompressible and inviscid fluid. Given a velocity field

U(x, y, t) = (u(x, y, t), v(x, y, t)),
and a pressure field p = p(x, y, t) defined at location (x, y) and time t, on some domain Ω with a solid wall boundary ∂Ω, the incompressible Euler evolution equations for U are

U t = -uU x -vU y -∇p, ∇ • U = 0.
The boundary condition for the velocity at solid wall boundaries is the no-flow condition U(x, y, t)

• n = 0 for (x, y) ∈ ∂Ω,
where n is the outward unit vector on ∂Ω.

We use a classical fractional-step projection method [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Peyret | Computational Methods for Fluid Flow[END_REF][START_REF] Brown | Accurate projection methods for the incompressible Navier-Stokes equations[END_REF]. At any given timestep n, we assume that the velocity at time n, U n and the fractional step pressure p n-1/2 are known at cell centres. In a first step, a provisional value U is computed using

U -U n ∆t = -A n+1/2 , (1)
where A n+1/2 is an approximation to the advection term [(U • ∇)U] n+1/2 . The new velocity U n+1 is then computed by applying an approximate projection operator to U which also yields the fractional step pressure p n+1/2 .

Poisson equation

The projection method relies on the Hodge decomposition of the velocity field as

U = U + ∇φ, (2)
where

∇ • U = 0 in Ω and U • n = 0 on ∂Ω. (3)
Taking the divergence of (2) yields the Poisson equation

∇ 2 φ = ∇ • U , (4)
while the normal component of (3) yields the boundary condition

∂φ ∂n = U • n on ∂Ω.
The divergence-free velocity field is then defined as

U = U -∇φ,
where φ is obtained as the solution of the Poisson problem (4). This defines the projection of the velocity U onto the space of divergence-free velocity fields.

In the context of the approximate projection method we are using here, the discrete formulation of the projection operator will depend on where the velocity field is discretised relative to the pressure field. We will use both an exact projection for face-centred advection velocities and an approximate projection for the final projection of the cellcentred velocities. The detail of these two projections does not influence the general description of the Poisson solver.

Relaxation operator

In practice, the spatially discretised Poisson problem results in a linear system of equations with the pressure at cell centres as unknowns:

L(φ) = ∇ • U . ( 5
)
where L is a discretisation of the Laplacian. This system can be solved through iterative methods (Jacobi, Gauss-Seidel) using a relaxation operator. If we consider a discretisation cell C of boundary ∂C, using the divergence theorem, the integration of (4) yields

C C C N d N d N d φ φ φ φ d φ d ∇ d φ ∇ d φ ∇ d φ (a) (b) (c)
∂C ∇φ • n = C ∇ • U , ( 6
)
where n is the outward unit normal of ∂C. In the case of a cubic discretisation cell, the discrete equivalent of ( 6) can be written as

d s d ∇ d φ = ha∇ • U , (7)
where d is the direction, s d the surface fraction in direction d and a the fluid volume fraction of the cell. Johansen and Colella [START_REF] Johansen | A cartesian grid embedded boundary method for poisson's equation on irregular domains[END_REF] have shown that this discretisation is second-order accurate if the right-hand side is defined at the geometric centre of the partial cell and the gradient at the geometric centre of the partial faces. Expressing the gradient at the geometric centre of the partial face requires interpolation of the full-face-centered gradients. While this is relatively simple on a regular Cartesian grid, this is more difficult within the adaptive framework we are using. Consequently we have chosen to use the full-face-centered gradient even in mixed cells. The following description thus applies to both full and mixed cells.

To construct the relaxation operator, we assume that the face gradient can be expressed as a linear function of the pressure at the centre of the cell

∇ d φ = α d φ + β d ,
where the α are constants and the β are linear functions of the values of the pressure in the adjacent discretisation cells.

In practice, three cases must be considered for the construction of the gradient operator (figure 3). If the neighbour of the cell in direction d, N d is at the same level and is a leaf cell, the gradient is simply ∇ d φ = (φ dφ)/h where φ d is the value of φ at the centre of N d . Using the notation above: 4. Second-order interpolation used for the gradient calculation at fine/coarse cell boundaries.

α d = -1/h and β d = φ d /h.

PSfrag replacements

N d N ⊥d N ⊥d N d φ d φ d φ 7 φ φ ⊥d φ ⊥d φ 3 φ 4 φ 6 ∇ d φ C Fig.
Figure 4 illustrates the case where N d is at a lower level (case 3.b). In order to maintain the second-order accuracy of the gradient calculation, it is necessary to use a threepoint interpolation procedure. The gradient ∇ d φ is computed by fitting a parabola through points φ 6 , φ and either φ 7 or φ d . By construction, N d is at the same level as C. If N d is a leaf cell, ∇ d φ can be expressed as

h∇ d φ = - φ 3 - φ d 5 + 8 15 φ 6 , (8)
where the value of the pressure at the centre of N d , φ d has been used. If N d is not a leaf cell, an interpolated value for the pressure φ 7 is constructed by averaging the values of its children closest to C (indicated by • in figure 4). The gradient is then given by

h∇ d φ = - 2 9 φ - 8 27 φ 7 + 14 27 φ 6 . (9)
The pressure φ 6 must itself be interpolated from φ d and from the values in the neighbouring cells in directions perpendicular to d. Due to the corner refinement constraint (figure 2.b), these cells ( N ⊥d and N ⊥d ) are guaranteed to be at the same level as N d .

The values φ 3 and φ 4 are derived using the same averaging procedure if N ⊥d and N ⊥d are not leaf cells. This leads to the following four cases:

φ 6 =                    15 16 φ d -3 32 φ ⊥d + 5 32 φ ⊥d if N ⊥d and N ⊥d are leaf cells 5 6 φ d -1 14 φ ⊥d + 5 21 φ 3 if N ⊥d is a leaf cell φ d -1 7 φ 4 + 1 7 φ ⊥d if N ⊥d is a leaf cell 8 9 φ d -1 9 φ 4 + 2 9 φ 3 otherwise (10)
The gradient ∇ d φ can still be expressed as a linear function of φ. The corresponding values of α d and β d can be calculated by using ( 8), ( 9) and [START_REF] Khokhlov | Fully threaded tree algorithms for adaptive refinement fluid dynamics simulations[END_REF].

In the third case, N d is at the same level but is not a leaf cell (figure 3.c). The gradient is simply constructed as minus the average of the gradients constructed from the children cells of N d closest to C (indicated by • in figure 3.c). These gradients are in turn computed using the interpolation technique described above (case 3.b). This approach ensures that the pressure gradient fluxes across coarse/fine boundaries are consistent. The extension to three dimensions is straightforward.

Once the α and β coefficients have been computed for each cell face of the domain, using (7) a relaxation operator can be defined as

R(φ, ∇ • U ) : φ ← ha∇ • U -d s d β d d s d α d . (11)
In the case where all the cells are on the same level and there are no solid boundaries (regular Cartesian grid), the operator reduces to the classical stencil

R(φ, ∇ • U ) : φ ← d φ d -h 2 ∇ • U n ,
where n is the number of directions (4 in 2D, 6 in 3D).

This operator, together with the interpolation procedure described above, has several desirable properties. It is second-order accurate in space at coarse/fine cell boundaries and uses a consistent flux estimation. In the case of cells cut by solid boundaries, the flux calculation is only first-order accurate in space, however.

Boundary conditions

Cells on the boundary of the domain or mixed cells may not have neighbours in all directions. If values for the pressure φ d are required in one of these directions, either by the gradient operator ∇ d φ or by the interpolation formula [START_REF] Khokhlov | Fully threaded tree algorithms for adaptive refinement fluid dynamics simulations[END_REF], they are set as equal to φ (the pressure at the centre of the cell considered). For cells entirely contained within the fluid, this is equivalent to a classical second-order accurate implementation of Neumann boundary conditions for the pressure.

M 0 M 1 M 2 M 3 M 4 M 5

Multilevel acceleration

The point relaxation defined by R can be accelerated using a multigrid technique [START_REF] Brandt | Guide to multigrid development, Multigrid Methods[END_REF][START_REF] Wesseling | An introduction to Multigrid Methods[END_REF]. When using quad/octrees, different choices are possible for the construction of the multilevel hierarchy. We have chosen to define a multilevel M l of depth l as the set of cells C which satisfy either of the conditions:

• level of C is equal to l,
• C is a leaf cell of level smaller than l.

An example of such a hierarchy is given in figure 5. This is probably not the best possible hierarchy for multigrid acceleration, in the sense that not all cells get coarser when moving from one level to the next. It is relatively easy to manually generate a possibly better hierarchy such as illustrated in figure 6. However, the systematic generation of such optimised hierarchies involves a set of rules substantially more complicated than the two conditions given above. In practice, if the simple rules are used, the traversal of the cells belonging to M l is straightforward to implement when using a pointer-based quad/octree structure.

Using this multilevel hierarchy, we apply a classical multigrid "V-cycle" using the correction form of the linear system [START_REF] Bell | A second-order projection method for the incompressible Navier-Stokes equations[END_REF].

L(φ + δφ) = ∇ • U ⇐⇒ L(δφ) = R with R = ∇ • U -L(φ).
The residual R is first computed on all the cells of the deepest level M L as

R L = ∇ • U - 1 ha d s d ∇ d φ. M 0 M 1 M 2 M 3 M 4 M 5
Fig. 6. Example of optimised multilevel hierarchy.

The residual is then transfered recursively on all the coarser levels as a volume weighted average

R l = i ah 2 R l+1 i ah 2
, where i designates the summation over all the children of the cell considered. The value of the pressure correction δφ is then computed exactly on the coarsest level. This value is used as the initial guess on the next finer level. Straight injection is used i.e. the initial guess δφ in each cell of M l is set as the value of δφ in its parent cell. The relaxation operator R is then applied a few times (using Jacobi iterations) and the resulting solution is used as initial guess on the next finer level. This is repeated recursively down to level L where the resulting correction is applied to φ. The whole Vcycle is repeated until the residual on the finest level is suitably small. This algorithm can be summarised as:

Compute R L on M L while aR L ∞ >
for l = L -1 to 0 Compute R l using weighted average of R l+1 end for Apply relaxation operator R(δφ, R 0 ) to M 0 down to convergence for l = 1 to L Get initial guess for δφ in cells at level l using straight injection from level l -1 Apply r times relaxations

R(δφ, R l ) to M l end for Correct φ on M L using δφ Compute R L on M L end while
It is important to note that, when applied to level M l , the relaxation operator should not use any cell of level larger than l (on which the solution for δφ is not yet defined). More specifically, when computing the gradient operator as described in the previous section, all the cells at level l must be considered as leaf cells even if they have children at level l + 1. This multigrid algorithm also differs from a classical implementation where a prerelaxation is applied before transferring the residual onto the coarser level [START_REF] Brandt | Guide to multigrid development, Multigrid Methods[END_REF]. In a classical multigrid the solution computed at each level is thus a correction to the correction at a deeper level. Such a scheme is difficult to implement on the multilevel quadtree hierarchy illustrated in figures 5 and 6 because, depending on the way the refined patches are laid out, it would require the storage of multiple corrections for the cells used as boundary conditions for refined patches. The scheme we propose solves this problem by dealing on all levels only with the correction to the pressure on the finest level. Of course, the convergence rate of such a "half" V-cycle is less than the convergence rate of the classical version, but tests have shown that the increased speed of such a simplified V-cycle more than compensate for the decrease in convergence rate.

In the following, we generally stop the V-cycle iterations when the maximum volumeweighted residual aR L ∞ is smaller than 10 -3 and we apply r = 4 iterations of the relaxation operator at each level.

Numerical validation

We are interested in two main properties of the multilevel Poisson solver: the speed of convergence for each V-cycle iteration and the spatial order of the method as the grid is refined. Given the way the relaxation operator is constructed, the method is expected to be globally second-order accurate on both regular and refined grids. If solid boundaries are used, the method should be first-order accurate near the solid boundaries and second-order accurate elsewhere.

We define the volume-weighted norm of a variable e as

ae p = i |e i | p a i h 2 i a i h 2 , (12)
where i designates the summation over all the leaf cells of the domain. An ∞-norm, ae ∞ , is the maximum over all the leaf cells of the absolute value of e. Knowing two solutions defined on domains of maximum refinement L 1 and L 2 , the rate of convergence in a given norm p can be estimated as

O p = log e 1 p e 2 p (L 2 -L 1 ) log 2 . ( 13
)
The convergence rate, O p = n, indicates n th -order accuracy, i.e. the leading term in the truncation error scales as O(h n ). A first test illustrates convergence on a regular Cartesian grid for a smooth pressure solution. We consider a square domain of size unity centred on the origin, with Neumann boundary conditions on all sides. The divergence is set in each cell as

∇ • U (x, y) = -π 2 (k 2 + l 2 ) sin(πkx) sin(πly), (14)
with k = l = 3. The exact solution of the Poisson equation with this source term is

φ(x, y) = sin(πkx) sin(πly) + κ, (15)
where κ is an arbitrary constant. The initial guess for the pressure is a constant field. Seven levels of refinement are used which results in a Cartesian discretisation of 2 7 × 2 7 = 128 × 128. We apply ten iterations of the V-cycle with r = 4 iterations of the relaxation operator at each level. Figure 7 illustrates the evolution of the maximum norm of the residual. A reduction factor (ratio of the residuals before and after the V-cycle) of about 25 per V-cycle is obtained.

To estimate the order of the solver, we solved the same problem on regular grids of increasing resolution. For each grid size, the norm of the error on the solution is calculated using the computed solution and the exact solution given by [START_REF] Almgren | A Cartesian grid projection method for the incompressible Euler equations in complex geometries[END_REF], where κ is taken as the average value of the computed pressure over the entire domain. Figure 8 illustrates the evolution of the error as a function of the depth of refinement L (i.e. a regular Cartesian grid of size 2 L × 2 L ). The order of convergence is computed as indicated above. As expected for this simple problem, the method shows second-order convergence in all norms.

For the moment, only the classical stencil on regular meshes has been used. In order to test the accuracy of the gradient operator in the case of coarse/fine mesh boundaries, we use the following test. A domain is first discretised with L -2 levels of refinement. Two more levels are then added only in the cells contained within a circle centred on the origin and of radius 1/4. The resulting discretisation for L = 6 is illustrated in figure 9. The same simple problem is then solved on this mesh. Figure 10 gives the convergence rate of the residual for a mesh with L = 7. The residual reduction factor is about 15 per V-cycle. The order of the solver for the same problem is illustrated in figure 11. Close to second-order convergence in all norms is obtained which confirms that the gradient operator described previously is second-order accurate at coarse/fine mesh boundaries.

In order to test the ability of the method in presence of solid boundaries, we set up a series of tests with a variety of solid geometries. The corresponding solutions of the Poisson equation are illustrated in figure 12. All problems use the source term defined by ( 14) where x and y are the coordinates of the geometric centre of the cell considered [START_REF] Johansen | A cartesian grid embedded boundary method for poisson's equation on irregular domains[END_REF] . A circular solid boundary centred on the origin and of radius 1/4 is used for (a). A star-shaped solid boundary defined in polar coordinates as r(θ) = 0.237 + 0.079 cos(6θ), is used in problem (b) and an ellipse centred on the origin measuring 3 4 × 5 8 in problem (c). All problems use Neumann conditions on all boundaries. Figure 13 illustrates the convergence speed for the three problems with L = 7 levels of refinement. The "star" problem (b) is notably more difficult to solve with an average reduction factor of only five per V-cycle. This is due to the limitation of the volume-of-fluid representation of the solid boundaries. As mentioned earlier, the features of the solid boundaries are only represented accurately if their spatial scale is comparable to the mesh size. For the "star" problem, while the geometry is represented correctly on the finest level, it is not well represented on all the coarser levels used by the multigrid procedure. Cases (a) and (b) do not have this problem because the smallest spatial scales of the solid boundaries (circle and ellipse) are comparable to the domain size.

The evolution of the error with resolution and the associated convergence order is given in figure 14. As the exact solution of the problem is not known analytically, Richardson extrapolation is used. That is, the error for a given level of refinement L is computed by taking the solution at level L + 1 as reference. A combination of solid boundaries and refinement is tested using a discretisation with L -2 levels of refinement on the whole domain plus two levels added only in cells cut by the solid boundary (a discretisation example is given in figure 15 for problem 12.b and L = 6). Figures 16 and17 illustrate the convergence speed and the order of the method using this discretisation.

The convergence is close to second-order (asymptotically in L) for all norms in all cases. The second-order convergence of the maximum error ae ∞ is surprising as the discretisation of the pressure gradient fluxes is only first-order accurate near solid boundaries (as described in section 4). This first-order error in the pressure gradient fluxes should lead to an O(1) truncation error of the Laplacian operator. Johansen and Colella [START_REF] Johansen | A cartesian grid embedded boundary method for poisson's equation on irregular domains[END_REF] have demonstrated that a scheme with an O(1) truncation error will lead to an O(h) error on the solution for the pressure, in contradiction to the O(h 2 ) convergence we obtain here.

To try to clarify this issue we present truncation and solution errors for the test case used in [START_REF] Johansen | A cartesian grid embedded boundary method for poisson's equation on irregular domains[END_REF][START_REF] Day | Embedded boundary algorithms for solving the Poisson equation on complex domains[END_REF]. The embedded boundary is defined by the curve, r(θ) = 0.30 + 0.15 cos 6θ.

The divergence is set in each full cell as

∇ • U (r, θ) = 7r 2 cos 3θ.
The exaction solution for this system is φ(r, θ) = r 4 cos 3θ. A mesh similar to figure 15 is used, with two levels of refinement added near the embedded boundary. In mixed cells, in order to be able to use Neumann boundary conditions at the solid surface while retaining the exact solution, the flux of the gradient of the exact solution through the boundary is subtracted from the divergence, giving

∇ • U (r, θ) = 7r 2 cos 3θ - s(n x ∇ x φ + n y ∇ y φ) ah ,
where s is the length of the embedded boundary contained within the cell, n x and n y are the components of the outward-pointing unit normal to the solid boundary. The gradients of the exact solution are defined as

∇ x φ = 4x 4 -3x 2 y 2 -3y 4 r , (16)
∇ y φ = xy (5x 2 + 9y 2 ) r , (17)
where x and y are the coordinates of the center of mass of the piece of embedded boundary contained within the cell. The results are summarized in figure 18. Figure 18.(a) gives the error norms and corresponding orders of convergence of the computed solution as functions of the level of refinement L. Figure 18.(b) illustrates the volume-weighted truncation error of the numerical Laplacian L defined in section 4.1. As expected the max-norm of the truncation error of the numerical Laplacian is O(1) due to the O(h) error in the pressure gradient fluxes in mixed cells, while the orders of the 1-and 2-norm are close to one and one-half respectively. However, while one would expect only first order convergence of the max-norm of the error on the solution, second-order convergence in all norms is obtained as illustrated in figure 18.(a). This confirms the results obtained for the previous tests and implies that second-order convergence in all norms can be obtained for practical problems even if the truncation error on the Laplacian is O(1). The discrepancy between our results and the theoretical study of Johansen and Colella could be explained if second-order converging errors in the bulk of the flow were always larger than first-order converging errors in mixed cells for all the tests we performed. This seems unlikely but if this were the case, it would be necessary to find more stringent test cases than used in this study or in [START_REF] Johansen | A cartesian grid embedded boundary method for poisson's equation on irregular domains[END_REF][START_REF] Day | Embedded boundary algorithms for solving the Poisson equation on complex domains[END_REF]. Further work in this direction would be useful.

Finally, figure 19 shows how the solver scales with problem size. Problem 12.a was solved on successively finer grids and the average residual reduction factor was com- puted as

aR 0 ∞ aR n ∞ 1 n
, where R i is the residual after i V-cycle have been applied and n is the total number of V-cycles (10 in this test). The residual reduction factor decreases approximately linearly with resolution level L. The computational cost of solving a problem with 2 2L = N 2 degrees of freedom (in 2D) thus scales as O(N 2 log N ) as expected from a multigrid scheme.

Advection term

We use a conservative formulation for the evaluation of the advection term. Given a cell C of boundary ∂C, using the divergence theorem and the non-divergence of the velocity field, the finite volume advection term A n+1/2 of (1) can be computed as

C A n+1/2 = C [(U • ∇)U] n+1/2 = C [∇ • (UU)] n+1/2 = ∂C U n+1/2 (U n+1/2 • n),
where n is the outward unit normal of ∂C. In the case of our cubic discretisation cell this can be written

ahA n+1/2 = d s d U n+1/2 d u n+1/2 d , ( 18
)
where

U n+1/2 d
is the velocity at the centre of the face in direction d at time n + 1/2 and u n+1/2 d is the normal component of the velocity at the centre of the face in direction d at time n + 1/2. In order to compute these time-and face-centred values, we use a Godunov procedure [START_REF] Bell | A second-order projection method for the incompressible Navier-Stokes equations[END_REF] i.e. the leading terms of a Taylor series of the velocity of the form

U n+1/2 d = U n + h 2 ∂ d U n + ∆t 2 ∂ t U n + O h 2 , ∆t 2 ,
where ∂ d designates the spatial derivative in direction d. Using the Euler equations, the temporal derivative can be replaced by spatial derivatives yielding

U n+1/2 d = U n + h 2 - ∆t 2 v n d ∂ d U n - ∆t 2 v n ⊥d ∂ ⊥d U n - ∆t 2 ∇p n ,
where ⊥d is the direction perpendicular to d in 2D (in 3D the sum over the two perpendicular directions) and v d is the velocity component in direction d at the centre of the cell. Given a cell face, two values of U n+1/2 d can be constructed, one for each cell sharing this face. In the original Godunov method for compressible fluids an unique value is constructed from these two values by solving a Riemann problem. In the incompressible case, simple upwinding is sufficient [START_REF] Bell | A second-order projection method for the incompressible Navier-Stokes equations[END_REF].

Following [START_REF] Martin | An adaptive cell-centered projection method for the incompressible Euler equations[END_REF][START_REF] Martin | A cell-centered adaptive projection method for the incompressible euler equations[END_REF] we use a simplified upwind scheme of the form

U n+1/2 d (C) = U n + h 2 min 1 -v n d ∆t h , 1 ∂ d U n - ∆t 2 v n ⊥d ∂ ⊥d U n ( 19
)
where ∂ ⊥d U n is the upwinded derivative in direction ⊥d

∂ ⊥d U n =      ∇ ⊥d U n if v n ⊥d < 0 ∇ ⊥d U n if v n ⊥d > 0 ( 20
)
∇ ⊥d is computed as in section 4.1 and ⊥d is the direction opposite to ⊥d. The cellcentred derivative ∂ d U n is computed by fitting a parabola through the centre of C and of its neighbours in directions d and d. If the neighbours are on different levels, an interpolation or averaging procedure similar to that presented in section 4.1 is used. In the case of neighbouring cells at the same level, this procedure reduces to the classical second-order accurate centred difference scheme. We also do not use any slope limiters on the derivatives as we do not expect strong discontinuities in the velocity field for incompressible flows. Slope limiters can easily be added in this scheme if necessary.

Given the time-and face-centred values U

U n+1/2 d (C) = U n+1/2 d (N d ) =              U n+1/2 d (C) if u n d > 0 U n+1/2 d (N d ) if u n d < 0 1 2 ( U n+1/2 d (C) + U n+1/2 d (N d )) if u n d = 0 ( 21
)
where N d is the neighbour of C in direction d. If N d is at a lower level (figure 20) and u n d ≤ 0, the value upwinded from N d at the centre of the face (marked by •) In order to compute the advection term using [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF], we first need to construct the faceand time-centred normal velocities u n+1/2 d

PSfrag replacements C N d N ⊥d U d U d U d U U u d C d
. If we want the method to be conservative, these normal velocities have to be discretely divergence-free. In a first step, normal velocities are constructed for both sides of each cell face using [START_REF] Saiki | Numerical simulations of a cylinder in a uniform flow: application of a virtual boundary method[END_REF] and [START_REF] Fadlun | Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations[END_REF] where v n d and v n ⊥d are the corresponding components of the centred velocity U n . The upwind state u d is then selected for each face using [START_REF] Martin | An adaptive cell-centered projection method for the incompressible Euler equations[END_REF] where u n d is obtained by linear interpolation of the relevant component of the centred velocities U n (C) and U n (N d ). To make this set of normal velocities divergence-free we then apply a projection step by solving

L(φ) = ∇ • u , (22)
where ∇ • u is the finite-volume divergence of the normal velocity field, expressed for each cell as

∇ • u = 1 ah d s d u d . (23)
By correcting u with the pressure solution, we obtain a set of face-and time-centred normal velocities

u n+1/2 d = u d -∇ d φ. (24)
When correcting the normal velocities, we also calculate a cell-centred value for the pressure gradient by simple averaging of face gradients

∇ d φ = ∇ d φ -∇ dφ 2 . (25)
To compute the advection term A n+1/2 , we first need to re-predict the face-centred velocities U n+1/2 d , this time using u n+1/2 d rather than averages from cell-centred values. Again, in a first step, normal velocities are constructed for both sides of each cell face using [START_REF] Saiki | Numerical simulations of a cylinder in a uniform flow: application of a virtual boundary method[END_REF] and [START_REF] Fadlun | Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations[END_REF] where we now take

v n d = (u n+1/2 d -u n+1/2 d
)/2. A unique value U d is then selected for each face using [START_REF] Martin | An adaptive cell-centered projection method for the incompressible Euler equations[END_REF]. A face-centred pressure gradient ∇φ is then computed by linear interpolation from the average cell-centred values ∇ φ(C) and ∇ φ(N d ) (or its children). The predicted value is then obtained as

U n+1/2 d = U d -∇φ.
Note that we could re-use the face-and time-centred normal velocities u n+1/2 d as predicted values (the tangential component would still need to be recalculated), however, we have found this approach to be unstable for flow around sharp angles. The spatial filtering of the pressure gradient provided by the averaging procedure seems to be necessary to ensure stability in this particular case.

Small-cell problem

To obtain the provisional cell-centred velocity field U using (1), it is necessary to divide the finite volume advection term [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF] by the volume of the cell (ah 2 ) to get A n+1/2 . This leads to the classical CFL stability condition

U ∆t ah ≤ 1,
which expresses the condition that a cell should not "overflow" during a given timestep. In the general case, the fluid fraction a can be arbitrarily small with a corresponding restrictive condition on the maximum timestep. This is traditionally referred to as the "small-cell problem". A number of approaches exist to work around this problem: cell merging [START_REF] Noh | Cel: a time-dependent, two-space-dimensional, coupled Eulerian-Lagrangian code[END_REF][START_REF] Quirk | An alternative to unstructured grids for computing gas dynamics flows around arbitrarily complex two-dimensional bodies[END_REF], redistribution [START_REF] Almgren | A Cartesian grid projection method for the incompressible Euler equations in complex geometries[END_REF] or special difference schemes [START_REF] Berger | A rotated difference scheme for Cartesian grids in complex geometries[END_REF]. We have chosen to use a simple cell-merging technique similar to that presented by Quirk [START_REF] Quirk | An alternative to unstructured grids for computing gas dynamics flows around arbitrarily complex two-dimensional bodies[END_REF]. At initialisation time, after the volume and area fractions have been computed, all small cells are assigned a pointer to their biggest neighbour B (as measured by the fluid fraction a). For a given timestep, the advection term ah 2 A n+1/2 is then computed in all cells as described above. To compute the advection update to the velocity, small cells are first grouped with adjacent mixed or full cells using the following recursive algorithm:

Group (G, C) if C does not already belong to any group then Add C to G if C is a small cell then Group (G, B(C)) end if for each direction d if N d is a small cell then Group (G, N d ) end if end for end if
where G is the resulting group of cells. For each group, the weighted averaged update is computed as

A n+1/2 G = G ah 2 A n+1/2 G ah 2 .
Each cell in the group then receives a fraction of the update proportional to its volume

A n+1/2 ≈ ah 2 G ah 2 A n+1/2 G .
This is equivalent to using a "virtual" cell formed by all the cells in the group. The CFL stability now depends on the total volume of the group of cells. In practice, choosing to define small cells as cells for which a < 1/2 ensured stability in all the cases we tested.

Approximate projection

While it is easy to formulate an exact projection operator for MAC (staggered, facebased) discretisation of the velocity field, it is difficult to do the same for a cell-centred discretisation. This is due to the spatial decoupling of the stencils used for the relaxation operator. This can cause numerical instabilities in the pressure field and makes efficient implementation of multigrid techniques difficult [START_REF] Lai | A projection method for reacting flows in the zero Mach number limit[END_REF][START_REF] Howell | An adaptive mesh projection method for viscous incompressible flow[END_REF]. Attempts to couple neighbouring pressure cells through asymmetric operators have been unsuccessful [START_REF] Strikwerda | Finite difference methods for the Stokes and Navier-Stokes equations[END_REF].

Drawing from these conclusions, Almgren, Bell and Szymczak [START_REF] Almgren | A numerical method for the incompressible Navier-Stokes equations based on an approximate projection[END_REF] dropped the requirement of exact discrete non-divergence of the projected cell-centred velocity field and proposed to use an approximate Laplacian operator well-behaved with respect to spatial coupling. Following Lai [START_REF] Lai | A projection method for reacting flows in the zero Mach number limit[END_REF], Minion [START_REF] Minion | A projection method for locally refined grids[END_REF] and Martin [START_REF] Martin | An adaptive cell-centered projection method for the incompressible Euler equations[END_REF] we use an approximate projection based on face-centred interpolation of the cell-centred velocity field. In a first step, face-centred normal components of the velocity are constructed by interpolation of the cell-centred provisional velocity U . This normal (MAC) velocity field is then projected using the exact projection operator (following steps [START_REF] Martin | A cell-centered adaptive projection method for the incompressible euler equations[END_REF] to [START_REF] Day | Embedded boundary algorithms for solving the Poisson equation on complex domains[END_REF]) and average cell-centred pressure gradients are constructed (using ( 25)). These pressure gradients are then used to correct U to obtain the approximately divergence-free velocity field U n+1 .

A detailed study of the stability of the approximate projection can be found in [START_REF] Rider | Approximate projection methods for incompressible flows: Implementation, variants and robustness[END_REF][START_REF] Almgren | Approximate projection methods: Part I. Inviscid analysis[END_REF]. The use of pressure filters was found to be necessary in some cases (long, quasistationary simulations) to avoid a gradual build-up of non-divergence-free velocity modes. We do not use pressure filters in the current version of the code but did not encounter any noticeable numerical instabilities for the various tests we performed.

It is also important to note that even if the resulting cell-centred velocity field is not exactly divergence-free, the face-centred normal advection field u n+1/2 d is exactly discretely divergence-free, so that the advection scheme is exactly conservative. This is particularly important for the treatment of variable density flows.

7 Adaptive mesh refinement Using a tree-based discretisation, it is relatively simple to implement a fully flexible adaptive refinement strategy.

In a first step, all the leaf cells which satisfy a given criterion are refined (as well as their neighbours when necessary, in order to respect the constraints described in figure 2). This step could be repeated recursively but we generally assume that the flow is evolving slowly (compared to the frequency of adaptation) so that only one pass is necessary.

In a second step, we consider the parent cells of all the leaf cells (i.e. the immediately coarser discretisation). All of these cells which do not satisfy the refinement criterion are coarsened (i.e. become leaf cells).

The values of the cell-centred variables for newly created or coarsened cells must be initialised. For newly coarsened cells, it is consistent to compute these values as the volume weighted average of the values of their (defunct) children, so that quantities such as momentum are preserved exactly. For newly created cells, the solution is less obvious. In particular, it is desirable that momentum and vorticity are locally preserved. Unfortunately, this is not simple to achieve in practice. We have chosen a simple linear interpolation procedure using the parent cell value and its gradients. Given a newly created cell C with parent cell P, the new cell-centred value v(C) is obtained as

v(C) = v(P) + ∆ x ∇ x v(P) + ∆ y ∇ y v(P),
where (∆ x , ∆ y ) are the coordinates of the centre of C relative to the centre of P. This formula guarantees local conservation of momentum but tends to introduce numerical noise in the vorticity field. A better choice may be higher-order interpolants such as bicubic interpolation.

On the new discretisation, there is no guarantee that the velocity field is divergencefree anymore. A projection step is then needed. To avoid the cost of an extra projection step when adapting the grid, we perform the grid refinement at the fractional timestep, using the provisional velocity field U , just before the approximate projection is applied.

Various choices are possible for the refinement criterion. An attractive option would be to use Richardson extrapolation to obtain a numerical approximation of the truncation error of the whole scheme [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF][START_REF] Martin | An adaptive cell-centered projection method for the incompressible Euler equations[END_REF]. For the moment, we use a simple criterion based on the norm of the local vorticity vector. Specifically, a cell is refined whenever

h ∇ × U max U > τ,
where max U is evaluated over the entire domain. The threshold value τ can be interpreted as the maximum acceptable angular deviation (caused by the local vorticity) of a particle travelling at speed max U across the cell.

The computational cost of this algorithm is small compared to the cost of the Poisson solver. It can be applied at every timestep with a negligible overall penalty (less than 5% of the total cost).

Numerical results

Following Minion [START_REF] Minion | A projection method for locally refined grids[END_REF] and Almgren et al. [START_REF] Almgren | A conservative adaptive projection method for the variable density incompressible Navier-Stokes equations[END_REF], we present two convergence tests illustrating the second-order accuracy of our method for flows without solid boundaries.

The first problem uses a square unit domain with periodic boundary conditions in both directions. The initial conditions are taken as u(x, y) = 1 -2 cos(2πx) sin(2πy), v(x, y) = 1 + 2 sin(2πx) cos(2πy).

The exact solution of the Euler equations for these initial conditions is u(x, y, t) = 1 -2 cos(2π(xt)) sin(2π(yt)), v(x, y, t) = 1 + 2 sin(2π(xt)) cos(2π(yt)), p(x, y, t) =cos(4π(xt))cos(4π(yt)).

As in [START_REF] Almgren | A conservative adaptive projection method for the variable density incompressible Navier-Stokes equations[END_REF] nine runs are performed on grids with L = 5, 6 and 7 levels of refinement (labelled "uniform") and with one (labelled r = 1) or two (labelled r = 2) additional 1 Errors and convergence orders in the x-component of the velocity for a simple periodic problem. levels added only within the square defined by the points (-0.25, -0.25) and (0, 0). The length of the run for each case is 0.5, the CFL number is 0.75. For each run both the L 2 and L ∞ norms of the error in the x-component of the velocity is computed using [START_REF] Almgren | A numerical method for the incompressible Navier-Stokes equations based on an approximate projection[END_REF] for both the whole domain (labelled "domain") and the refined region only (labelled "patch"). Table 1 gives the errors and order of convergence obtained.

Patch L 2 L = 5 O 2 L = 6 O 2 L = 7 r =
L = 5 O ∞ L = 6 O ∞ L = 7 r = 1 1.73e-2 1.82 4.89e-3 1.91 1.30e-3 r = 2 1.58e-2 1.41 5.96e-3 1.84 1.66e-3 Domain L ∞ L = 5 O ∞ L = 6 O ∞ L = 7 Uniform 1.
Close to second-order convergence is obtained (asymptotically in L) for the L 2 and L ∞ norms on both uniform and refined domains. The values obtained are comparable to that in [START_REF] Minion | A projection method for locally refined grids[END_REF][START_REF] Almgren | A conservative adaptive projection method for the variable density incompressible Navier-Stokes equations[END_REF]. The error in the refined patch is comparable to the error at the resolution of the base grid. This is expected, given the arbitrary placement of the refined patch, the error is controlled mainly by the surrounding coarse cells.

The second test is the four-way vortex merging problem of Almgren et al. [START_REF] Almgren | A conservative adaptive projection method for the variable density incompressible Navier-Stokes equations[END_REF]. It demonstrates the convergence of the method when refinement is placed appropriately.

Four vortices are placed in the unit-square, centred at (0, 0), (0.09, 0), (-0.045, 0.045 √ 3)

and (-0.045, -0.045 √ 3) and of strengths -150, 50, 50, 50 respectively. The profile of each vortex centred around (x i , y i ) is 1 + tanh(100(0.03r i ))

2 ,

where r i = (xx i ) 2 + (yy i ) 2 . To initialise the velocity field, we use this vorticity as the source term in the Poisson equation for the streamfunction ψ

∇ 2 ψ = ∇ × U .
Each component of the velocity field is then calculated from the streamfunction. Noflow boundary conditions are used on the four sides of the domain and the simulations are ran to t = 0.25 using a CFL of 0.9.

Five different discretisations are used, each time with up to L levels of refinement: a uniform grid, a grid using static refinement in concentric circles of decreasing radius and the dynamic adaptive refinement described in section 7. The "circle" grid is constructed by starting from a uniform grid with four levels of refinement and by successively adding one level to all the cells contained within circles centred on the origin and of radii:

• For the dynamically refined grid, the vorticity-based criterion is applied at every timestep with a threshold τ = 4 × 10 -3 . As we do not have an analytical solution for this problem, Richardson extrapolation is used.

Figure 21 illustrates the evolution of the vorticity and of the adaptively refined grid for L = 8. The most refined level closely follows the three outer vortices as they orbit the central one. Far from the vortices, a very coarse mesh is used (l = 3). One may note a few isolated patches of refinement scattered at the periphery of the outer vortices. They are due to the numerical noise added to the vorticity by the interpolation procedure necessary to fill in velocity values for newly created cells. As mentioned in section 7, this could be improved by using higher-order interpolants. This numerical noise is small enough that it does not compromise the convergence properties of the adaptive method (as shown below).

Table 2 summarises the results obtained for the first twelve calculations. For fine enough grids close to second-order convergence is obtained for both norms and for the three discretisations used. The norms of the error on the various grids are also comparable for a given resolution.

L = 6 O ∞ L = 7 O ∞ L = 8 O ∞ L = 9
Uniform 4.46e-1 1.25 and for two levels of refinement. A PC-compatible Pentium 350 MHz machine was used. The total number of leaf cells advanced for the whole calculation is given as well as the corresponding average speed. For L = 8, a speedup of about nine is obtained when using the statically refined "circle" grid, and thirteen when using the adaptive technique. Both the "circle" and adaptive discretisations are notably slower (per cell) than the uniform discretisation. This is mainly due to the interpolations necessary to compute the pressure gradient at coarse/fine cell boundaries when solving the Poisson equation (figure 4). It is also interesting to note that the CPU times obtained are very close to those reported by Almgren et al. [START_REF] Almgren | A conservative adaptive projection method for the variable density incompressible Navier-Stokes equations[END_REF] for the same problem (keeping in mind that they used a Cartesian AMR technique on a DEC Alpha computer and solved a viscous flow).

To demonstrate the convergence properties of the method in the presence of solid boundaries, we use a test case initially presented by Almgren et al. [START_REF] Almgren | A Cartesian grid projection method for the incompressible Euler equations in complex geometries[END_REF]. A diverging t = 0.05 t = 0.15 t = 0.25 channel is constructed in a 4×1 domain by restricting the fluid flow through the curves Wind is coming at a right angle from the right of the vessel. The stream ribbons and cross-section at sea level are coloured according to the norm of the velocity. The spatial resolution is about 50 centimetres near the ship and is adapted dynamically (using the vorticity criterion) down to a minimum scale of one meter elsewhere in the flow. The resulting mesh is composed of about 350,000 leaf cells in established regime. Figure 24 shows a vertical and horizontal cross-section through the adapted octree mesh for the same timestep. We are in the process of comparing these results to experimental measurements, which will be the subject of a future publication.

Conclusion

The combination of a quad/octree discretisation, an approximate projection method, a multigrid Poisson solver and a volume-of-fluid embedded description of solid boundaries proves to be a feasible and efficient technique for the numerical solution of the time-dependent incompressible Euler equations. This approach differs from the classical Cartesian AMR technique [START_REF] Berger | Adaptive mesh refinement for shock hydrodynamics[END_REF][START_REF] Almgren | A numerical method for the incompressible Navier-Stokes equations based on an approximate projection[END_REF][START_REF] Howell | An adaptive mesh projection method for viscous incompressible flow[END_REF] by treating the connection between levels of refinement at the cell operator level rather than through boundary conditions between refined patches. These operators can be designed to be spatially second-order accurate and to use consistent (conservative) flux estimations at coarse/fine boundaries. This fine-grained description allows almost full flexibility in the placement and shape of refined regions. Moreover, the refinement and coarsening process is naturally implemented by the quad/octree structure and does not need specialised algorithms for grid generation [START_REF] Berger | An algorithm for point clustering and grid generation[END_REF]. The mesh adaptation can thus be performed for every timestep with minimum overhead.

The price to pay for this flexibility is the loss of the array-based, cache-and accessefficient structured grids which are at the core of the Cartesian AMR technique. While more thorough investigation would be necessary, we show that similar performances to AMR can be achieved using our quad/octree approach. Moreover, we believe that in the case of small and complicated structures (such as interfaces between fluids or shocks) the flexibility of this approach can more than compensate for this overhead (given that a Cartesian AMR technique would require a large number of refined patches to cover the small structures, leading to substantial overheads in boundary conditions and most probably to the loss of cache-efficiency).

Future developments include extension to the incompressible variable-density Navier-Stokes equations and interfacial flows, using VOF [START_REF] Gueyffier | Volume of fluid interface tracking with smoothed surface stress methods for three-dimensional flows[END_REF] and marker techniques [START_REF] Popinet | A front tracking algorithm for the accurate representation of surface tension[END_REF]. Using sub-cycling in time on different levels of refinement [START_REF] Almgren | A conservative adaptive projection method for the variable density incompressible Navier-Stokes equations[END_REF] would also be a useful extension of the algorithm presented.

Finally, by providing an open source version of the code which can be freely redistributed and modified [START_REF] Popinet | The Gerris Flow Solver[END_REF], we hope to encourage research and collaboration in this field.

Fig. 1 .

 1 Fig. 1. Example of quadtree discretisation and corresponding tree representation.

Fig. 2 .

 2 Fig.2. Additional constraints on the quadtree discretisation. The refinement necessary to conform to the given constraint is indicated by the dotted lines.

( a )

 a for any given cell, efficient access to neighbouring cells. (b) for any given cell, efficient access to cell level and spatial coordinates. (c) efficient traversal of:• all leaf cells,• all cells at a given level, • all mixed cells.

Fig. 3 .

 3 Fig. 3. Three cases for face-centered gradient calculation. (a) Cells at the same level. (b) Fine-coarse boundary. (c) Coarse-fine boundary.

Fig. 5 .

 5 Fig. 5. Example of simple multilevel hierarchy.

Fig. 7 .

 7 Fig. 7. Speed of convergence of the Poisson solver for a simple problem, L = 7. (a) Evolution of the residual. (b) Reduction factor.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Order of convergence of the Poisson solver for a simple problem. (a) Evolution of the error and (b) order of convergence as functions of resolution.

Fig. 10 .

 10 Fig. 10. Speed of convergence of the Poisson solver for a simple problem discretised using a mesh similar to figure 9 with L = 7. (a) Evolution of the residual. (b) Reduction factor.

Fig. 11 .

 11 Fig. 11. Order of convergence of the Poisson solver for a simple problem discretised using a mesh similar to figure 9. (a) Evolution of the error and (b) order of convergence as functions of resolution.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Contour plots of the solution of Poisson problems with solid boundaries.

Fig. 14 .

 14 Fig. 14. Evolution of the error and associated convergence order for Poisson problems with solid boundaries.

Fig. 15 .Fig. 16 .

 1516 Fig. 15. Boundary-refined mesh for problem 12.b, L = 6.

Fig. 17 .

 17 Fig. 17. Evolution of the error and associated convergence order for Poisson problems with refined solid boundaries.

Fig. 18 .

 18 Fig. 18. Evolution of the error and associated convergence order for the Neumann Poisson problem of [23,24] using locally refined solid boundaries. (a) Error on the solution. (b) Volume-weighted truncation error on the Laplacian of the exact solution.

Fig. 19 .

 19 Fig. 19. Average residual reduction factor for problem 12.a as a function of resolution L.

  ) we then choose the upwind state

Fig. 20 .

 20 Fig. 20. Upwinding in the case of neighbouring cells at different levels. Linear interpolation is used to derive the value on the right side of C d .

  L = 6: 0.25, 0.15 • L = 7: 0.25, 0.2, 0.15 • L = 8: 0.25, 0.2, 0.175, 0.15 • L = 9: 0.25, 0.2, 0.175, 0.1625, 0.15 • L = 10: 0.25, 0.225, 0.2, 0.175, 0.1625, 0.15

Fig. 21 .

 21 Fig. 21. Contour plots of vorticity (left) and adaptive grids used (right) for the four-way vortex merging calculation. The lines on the pictures in the right column represent the boundaries between levels of refinement (with a maximum of L = 8 levels).

Fig. 22 .

 22 Fig. 22. Contour plot of the error on the x-component of the velocity estimated for a solution with L = 6 levels of refinement.practical engineering-type problem. The air flow around the vessel RV Tangaroa of the National Institute of Water and Atmospheric Research has been simulated by solving the 3D time-dependent incompressible Euler equations around a CAD model. Figure23is a snapshot in time of the developed turbulent flow. Wind is coming at a right angle from the right of the vessel. The stream ribbons and cross-section at sea level are coloured according to the norm of the velocity. The spatial resolution is about 50 centimetres near the ship and is adapted dynamically (using the vorticity criterion) down to a minimum scale of one meter elsewhere in the flow. The resulting mesh is composed of about 350,000 leaf cells in established regime. Figure24shows a vertical and horizontal cross-section through the adapted octree mesh for the same timestep. We are in the process of comparing these results to experimental measurements, which will be the subject of a future publication.

Fig. 23 .

 23 Fig. 23. Airflow around RV Tangaroa. The stream ribbons and cross-section at sea level are coloured according to the norm of the velocity.

Fig. 24 .

 24 Fig. 24. Adaptive mesh. The horizontal and vertical cross-sections illustrate the three-dimensional adaptive octree.

Table 3

 3 gives the CPU time and the size of the problems solved for all three grids

	Domain				L 2
		L = 6	O 2	L = 7	O 2	L = 8	O 2	L = 9
	Uniform 2.61e-2 1.31 1.05e-2 1.97 2.68e-3 2.11 6.19e-4
	Circle	2.61e-2 1.33 1.04e-2 1.96 2.68e-3 1.98 6.81e-4
	Adaptive 2.66e-2 1.35 1.04e-2 2.07 2.47e-3 2.05 5.96e-4
	Domain				L ∞

  Timings for uniform, circle and adaptive grids for the four-way vortex merging problem.

			1.87e-1 1.95 4.84e-2 1.85 1.34e-2
	Circle	4.49e-1 1.27 1.86e-1 1.94 4.85e-2 1.87 1.33e-2
	Adaptive 4.45e-1 1.26 1.86e-1 1.94 4.86e-2 1.83 1.37e-2
	Table 2
	Errors and convergence orders in the x-component of the velocity for the four-way vortex
	merging problem.
			CPU Time	Cells advanced
			Total (s) µs/cell	Number
	Uniform L = 8	1486	167	8,912,896
	Circle	L = 8	166	222	746,368
	Adaptive L = 8	117	286	409,632
	Uniform L = 9	13,034	166	78,643,200
	Circle	L = 9	1024	183	5,608,960
	Adaptive L = 9	764	326	2,342,200
	Table 3

Acknowledgments

This work was funded by the New Zealand Foundation for Research, Science and Technology and by the Marsden Fund of the Royal Society of New Zealand.

y top and y bot , defined as

and y top = 1y bot , with y 1 = 0.2 and y 2 = 10 -6 . Neumann boundary conditions for the pressure are set at the inlet (x = 0) and at the solid boundaries. A fixed unity inflow velocity is set at the inlet and simple outflow boundary conditions at the outlet (the pressure and the gradients of all the components of the velocity are set to zero at x = 4). The simulations are ran to t = 1 using a CFL of 0.8. Three simulations are performed on uniform grids with L = 5, 6 and 7 levels of refinement. Table 4 and 5 show the errors and convergence rates obtained. As in [START_REF] Almgren | A Cartesian grid projection method for the incompressible Euler equations in complex geometries[END_REF] we calculate errors both on the full domain ("All cells") and on the part of the domain covered by cells at level 5 entirely contained within the fluid ("Full level 5 cells"). Columns labelled "5-6" give the error computed on the mesh with 5 levels of refinement using the mesh with 6 levels of refinement as reference (and similarly for columns labelled "6-7"). For both components of the velocity close to first-order convergence is obtained for the L ∞ norm and close to second-order convergence for the L 1 norm, as expected from a solution globally second-order accurate but first-order accurate at the boundaries. Figure 22 confirms that the error is concentrated near solid boundaries. The maximum error in either component is small (less than one percent of the magnitude of the velocity).

Finally, we present an application of the three-dimensional version of the code to a