domain: math phys info info.info-ds

text:
Victorien Menier

Adrien Loseille

Frédéric Alauzet

Multigrid Strategies Coupled with Anisotropic Mesh Adaptation

This paper is an attempt to combine a full multigrid (FMG) algorithm with the classical anisotropic mesh adaptation loop. The FMG algorithm aims at improving the convergence rate of the flow solver, while mesh adaptation tends to improve the tradeoff between CPU time and accuracy of the solution. Coupling these two methods consists, at each stage of the adaptation loop, in running a multigrid simulation using the previously adapted coarser meshes. The main benefits of this coupling are that (i) the convergence of mesh adaptation is improved thanks to the properties stated by the FMG theory, and (ii) anisotropic coarsening is automatically handled by mesh adaptation. First, the implementation and the validation of the multigrid procedures are described, and then the coupling is detailed. Results are presented in 2D and 3D for various flow prescriptions, including a real-life configuration.

Introduction

Multigrid methods have been widely used for algebraic problems since their original development over thirty years ago. [START_REF] Brandt | Algebraic multigrid (AMG) for sparse matrix equations[END_REF][START_REF] Brezina | Adaptive algebraic multigrid[END_REF] Interest in these methods has since become even greater, thanks to their ability to efficiently solve problems arising from partial differential equations. For instance, during a CFD simulation using an implicit time integration, a linear system is solved at each iteration of the flow solver. The Newton method used for solving this linear system can be dramatically accelerated using multigrid algorithms, thus improving the convergence rate of the whole simulation. [START_REF] Trottenberg | Multigrid[END_REF][START_REF] Mccormick | Multilevel Adaptive Methods for Partial Differential Equations[END_REF][START_REF] Sorensen | A multigrid accelerated time-accurate inviscid compressible fluid flow solution algorithm employing mesh movement and local remeshing[END_REF] Multigrid simulations require a sequence of meshes of different resolutions, and methods of generating these meshes can be classified into three main categories. The simplest manner to generate coarser meshes is to build a hierarchical set of embedded meshes, which presents serious limitations, one of which is that bad quality elements are created during the process, which badly impacts the numerical solution. [START_REF] Carré | Multigrid Strategies for CFD Problems on Non-Structured Meshes[END_REF] Another method is the volume agglomeration technique, which consists in agglomerating the finite volume cells of the dual mesh. [START_REF] Francescatto | A Semi-Coarsening Strategy for Unstructured MG with Agglomeration[END_REF][START_REF] Koobus | Unstructured volume-agglomeration MG : solution of the Poisson equation[END_REF] In this paper, we use a third approach: the generation of non-nested unstructured coarse meshes. [START_REF] Guillard | Node-Nested Multi-Grid With Delaunay Coarsening[END_REF][START_REF] Mavriplis | Adaptive meshing techniques for viscous flow calculations on mixed element unstructured meshes[END_REF][START_REF] Carré | Multigrid Strategies for CFD Problems on Non-Structured Meshes[END_REF][START_REF] Mller | Anisotropic adaptation and multigrid for hybrid grids[END_REF] The metric corresponding to the initial mesh is scaled, and used as an input by an adaptive remesher. Note that stiff problems (shock waves, boundary layers, etc.) can cause a breakdown in efficiency of multigrid methods, [START_REF] Mavriplis | Multigrid Strategies for Viscous Flow Solvers on Anisotropic Unstructured Meshes[END_REF] due to high anisotropy. This is faced using specific anisotropic coarsening strategies. [START_REF] Brandt | Multigrid Techniques[END_REF][START_REF] Francescatto | A Semi-Coarsening Strategy for Unstructured MG with Agglomeration[END_REF][START_REF] Mavriplis | Multigrid Strategies for Viscous Flow Solvers on Anisotropic Unstructured Meshes[END_REF][START_REF] Mavriplis | Adaptive meshing techniques for viscous flow calculations on mixed element unstructured meshes[END_REF][START_REF] Mesri | Automatic coarsening of three dimensional anisotropic unstructured meshes for multigrid applications[END_REF] The full multigrid (FMG) algorithm [START_REF] Carré | Multigrid Strategies for CFD Problems on Non-Structured Meshes[END_REF][START_REF] Trottenberg | Multigrid[END_REF] consists in several stages. At each stage, a multigrid simulation is performed, the output solution is linearly interpolated from one stage to the next and used as a restart solution. The complexity of the finest mesh used is increased at each stage. This algorithm presents similarities with the classical mesh adaptation loop, which has proved its efficiency for improving the accuracy of the solution while decreasing the computational wall clock time. [START_REF] Castro-Díaz | Anisotropic Unstructured Mesh Adaptation for Flow Simulations[END_REF][START_REF] Frey | Anisotropic mesh adaptation for CFD computations[END_REF][START_REF] Gruau | 3D tetrahedral, unstructured and anisotropic mesh generation with adaptation to natural and multidomain metric[END_REF][START_REF] Li | 3D anisotropic mesh adaptation by mesh modification[END_REF][START_REF] Alauzet | High Order Sonic Boom Modeling by Adaptive Methods[END_REF] This paper presents a coupling of the mesh adaptation loop and the FMG algorithm, both of which can benefit from it. Multigrid algorithms use non-adapted meshes (except in the boundary layer regions), which are not optimal in terms of sizes and directions, for capturing physical phenomena which are anisotropic or located in small areas of the computational domain. Moreover, anisotropic coarsening can automatically be handled using mesh adaptation. According to the FMG theory, [START_REF] Carré | On the Application of FMG to Variational Approximation of Flow Problems[END_REF] multigrid schemes can improve the convergence of the solution in the mesh adaptation loop. Indeed, it states that, if the solution is fully converged at stage 1, then it is sufficient to converge the solution by one order of magnitude at stages 2, 3, etc. in order to achieve the global convergence on the finest mesh. This is of main interest for the adaptive process, since it states that it is not always necessary to fully converge the solution at intermediary steps of the mesh adaptation loop. Thus, it provides clues on how to reduce the number of solver iterations at intermediary steps, without affecting the accuracy of the final solution.

The paper is organized as follows. In Section I, the implementation of the non-adapted multigrid procedure is presented, together with its validation using 2D and 3D cases. First, the generation of coarser non-nested meshes using isotropic and anisotropic coarsening is discussed. Then, the spatial discretization and the implicit time integration of the modeling equations are described, which lead to solving a linearized system at each flow solver iteration. The multigrid acceleration of the Newton method is then presented, and several multigrid cycles are introduced.

In Section II, the FMG algorithm and the classical mesh adaptation are introduced, and their coupling is detailed. It consists in replacing each single-grid computation in the classical mesh adaptation loop by a multigrid procedure using the previously adapted meshes as coarser meshes. Finally, Section II.4 presents adaptative multigrid simulations, including a 3D industrial configuration.

I. Validation of the Implicit Multigrid Procedure

This section describes how we have implemented and validated the uniform multigrid procedure in Wolf, our in-house flow solver. We start by describing the generation of non-nested coarser meshes, which consists in scaling the metric field and using the usual adaptive remesher. The need for anisotropic coarsening is discussed, as well as how it can be done naturally using mesh adaptation. Then, the implementation of the different multigrid schemes used is detailed. Starting from the modeling equations, we introduce the linearized system that is solved after each time iteration and how the Newton method is accelerated and robustified using multigrid. Finally, several numerical examples are provided for validation purposes.

I.1. Uniform Coarse Meshes Generation

This section describes the generation of a hierarchical set of coarser meshes, in the uniform case (i.e. non adapted meshes). Starting from an initial finest mesh H h (whose representative edge size is h), we want to generate coarser meshes H 2h , H 4h , H 8h etc. suitable for a multigrid computation. Various algorithms exist to generate such meshes, including agglomeration techniques [START_REF] Koobus | Unstructured volume-agglomeration MG : solution of the Poisson equation[END_REF][START_REF] Francescatto | A Semi-Coarsening Strategy for Unstructured MG with Agglomeration[END_REF][START_REF] Guillard | An Aggregation Multigrid Solver for Convection-diffusion Problems on Unstructured Meshes[END_REF] and non-nested methods [START_REF] Guillard | Node-Nested Multi-Grid With Delaunay Coarsening[END_REF][START_REF] Mavriplis | Multigrid techniques for unstructured meshes[END_REF][START_REF] Mavriplis | Multigrid Strategies for Viscous Flow Solvers on Anisotropic Unstructured Meshes[END_REF] (see Figure 2). We chose to generate unstructured meshes containing simplicial elements (triangles in 2D, surface triangles and tetrahedra in 3D), which are not necessarily nested (see Section I.1.4), i.e. the set of vertices of a coarse mesh H 2h is not necessarily included in the set of vertices of H h . The use of nested meshes in a multigrid computation is investigated in Section II.4. The coarse mesh generation process is the following: (i) the geometric metric field of H h is computed, (ii) a scaling factor is applied to the metric and (iii) an adaptive remesher is used with respect to a set of constraints such as the conservation of the curvature of a surface.

I.1.1. Isotropic Scaling of the Geometric Metric of H h

At each vertex of a mesh, a metric tensor field prescribes sizes and directions of elements for the mesh generation procedure. More formally, a metric tensor M in R n is an n × n symmetric definite positive matrix, according to which the scalar product of two vectors u and v in R n is defined as:

u, v M = u, Mv = t uMv ∈ R.
So, the associated norm of a vector in R n is defined as:

u M = u, u M = √ t uMu,
which measures the length of the vector u in the metric M.

Adapting a mesh H according to a metric tensor field M consists in generating a unit mesh H new with respect to M. A mesh H new is unit with respect to M, if every one of its edge ab has a length ℓ M (ab) equal to 1 in M, and if every one of its tetrahedra K has a volume |K| M in M equal to √ 2/12. The length and the volume in M are given by:

ℓ M (ab) = 1 0 t ab M(a + tab) ab dt |K| M = K √ det M dK
The first step of mesh coarsening is to compute the metric M geo (H h ), according to which H h is unit. In order words, M geo (H h ) reflects the sizes and anisotropy of the elements of H h . M geo (H h ) is then multiplied by a scaling factor: 1/4 for H 2h , 1/16 for H 4h etc. Using a local mesh operator, the set of coarse meshes is then obtained with respect to these scaled metrics.

I.1.2. Adaptive Remesher

Given H h and a scaled metric M nh , vertex insertions and edge collapses are performed recursively until all the edges have a length in the interval [ √ 2, 1/ √ 2] with respect to M nh . Then, an optimization procedure is applied, based on vertex smoothing, edges and faces swapping. During this optimization step, we try to improve the overall quality of the mesh (as opposed to the unit mesh step where the convergence to a prescribed length distribution is the main concern). The quality in a metric M is given by:

Q M (K) = 36 3 1 3 |K| 2 3 M 6 i=1 ℓ 2 M (e i ) ∈ [0, 1],
and a "perfect" element has a quality equal to 1.

All the aforementioned mesh modification operators are based on a unique cavity-based operator revisited in an anisotropic context. [START_REF] Loseille | Serial and Parallel Mesh Modification Through A Unique Cavity-Based Primitive[END_REF] This operator also performs surface-based operations. In such a case, a local surface approximation is computed based on a user-provided background surface mesh and metric. This allows us to limit, according to an input parameter, the approximation error of the surface.

I.1.3. Anisotropic Coarsening

A breakdown in efficiency of the multigrid methods can be observed when dealing with high anisotropy. Turbulent viscous flow simulations, for instance, require mesh elements whose maximum aspect ratio can be in the order of 10 -6 in the boundary layer regions. Anisotropic elements are also necessary to accurately capture shock waves during transonic simulations. In order to prevent this breakdown in efficiency, existing isotropic coarsening techniques [START_REF] Guillard | Node-Nested Multi-Grid With Delaunay Coarsening[END_REF][START_REF] Mavriplis | Adaptive meshing techniques for viscous flow calculations on mixed element unstructured meshes[END_REF] were extended to the anisotropic case. [START_REF] Brandt | Multigrid Techniques[END_REF][START_REF] Francescatto | A Semi-Coarsening Strategy for Unstructured MG with Agglomeration[END_REF][START_REF] Mavriplis | Multigrid Strategies for Viscous Flow Solvers on Anisotropic Unstructured Meshes[END_REF][START_REF] Mavriplis | Adaptive meshing techniques for viscous flow calculations on mixed element unstructured meshes[END_REF][START_REF] Mesri | Automatic coarsening of three dimensional anisotropic unstructured meshes for multigrid applications[END_REF] The general idea is to coarsen the mesh only in the directions that are perpendicular to a direction of anisotropy. To do so, only the scaling of the geometrical metric M geo (H h ) differs from the isotropic case.

We chose to compute the anisotropic scaled metric tensor field of M geo (H h ) as follows [START_REF] Mesri | Automatic coarsening of three dimensional anisotropic unstructured meshes for multigrid applications[END_REF] (in 3D). Let M k be the metric tensor of M geo (H h ) associated to a vertex v k , and λ 1 , λ 2 , λ 2 its eigenvalues, and h 1 , h 2 , h 3 the corresponding sizes (h i = (λ i ) -1/2 ). We start by ordering these sizes.

h 1 ≤ h 2 ≤ h 3 .
Then, the new coarsened sizes are computed:

h i,coarse = max(h i , min(c × h i , h i-1 )) ,
where c is the coarsening factor. This defines the new coarsened eigenvalues:

λ i,coarse = h -1/2
i,coarse . This algorithm will produce the same coarse metric in the isotropic regions of the mesh, as in the anisotropic regions, only the directions where the mesh size is minimal are scaled.

The mesh adaptation process automatically provides a kind of anisotropic coarsening, thanks to the considered L p error estimate. This is one of the advantages of coupling it with an FMG algorithm. We recall [START_REF] Dervieux | High-order adaptive method applied to high-speed flows[END_REF] that the convergence order on a discontinuity is 1/p in L p norm. In other words, dividing the mesh size around a discontinuity by a factor s will produce an error on the discontinuity divided by s 1 p . If, however, the mesh size is divided by a factor s = 2 everywhere but around the discontinuity, and divided by a factor s = 2 2p in the direction of the discontinuity, then the error will be divided by a factor 4 and for a cost of N → 2N + 2 2p .

I.1.4. Generation of Nested Meshes

Generating embedded meshes is a simple way to build a hierarchical nested set of coarse meshes for multigrid, such as depicted in Figure 2. It consists in first generating an initial coarse mesh, which is then refined by element subdivision. This method has a major inconvenience, since the quality of the meshes generated decreases as they are refined. Indeed, as element subdivisions are iteratively performed, patterns corresponding to the coarsest mesh elements appear. These patterns may influence the computation, as they can act as artificial internal boundaries, 7 as illustrated in Figure 3.

I.2. Spatial Discretization and Implicit Time Integration in Wolf

This section introduces the linearized system that is to be solved after each time iteration. More details are given in. [START_REF] Menier | CFD validation and adaptivity for viscous flow simulations[END_REF] The modeling equations we consider are the Euler equations, which read:

             ∂ρ ∂t + ∇ • (ρu) = 0 , ∂(ρu) ∂t + ∇ • (ρu ⊗ u) + ∇p = 0 , ∂(ρe) ∂t + ∇ • ((ρe + p)u) = 0 ,
where ρ is the density, p the pressure, u = (u, v, w) the velocity vector, e the total energy per mass unit, and E = ρe the total energy per volume unit. This system can be written in vectorial form:

W t + F 1 (W ) x + F 2 (W ) y + F 3 (W ) z = 0 ( 1
)
where W is the nondimensionalized conservative variables vector:

W = (ρ, ρu, ρv, ρw, ρE) T ,
and

F(W ) = (F 1 (W ), F 2 (W ), F 3 (W )
) are the convective (Euler) flux functions:

F 1 (W ) = (ρu, ρu 2 + p, ρuv, ρuw, u(ρE + p)) T F 2 (W ) = (ρv, ρuv, ρv 2 + p, ρvw, v(ρE + p)) T F 3 (W ) = (ρw, ρuw, ρvw, ρw 2 + p, w(ρE + p)) T .
The spatial discretization of the fluid equations is based on a vertex-centered finite element/finite volume formulation on unstructured meshes. It combines HLLC upwind schemes for computing the convective fluxes and second order space accuracy is achieved through a piecewise linear interpolation based on the Monotonic Upwind Scheme for Conservation Law (MUSCL) procedure which uses a particular edge-based formulation with upwind elements. A specific slope limiter is employed to damp or eliminate spurious oscillations that may occur in the vicinity of discontinuities. [START_REF] Alauzet | High Order Sonic Boom Modeling by Adaptive Methods[END_REF] The implicit time discretization of (1) reads:

|C i | δt n i W n+1 i -W n i = R i (W n+1 ) 5 of 20
which is linearized as:

|C i | δt n i I d - ∂R i ∂W (W n ) W n+1 i -W n i = R i (W n )
where ∂Ri ∂W (W n ) contributes the i th line of the matrix and using a set of linearizations. [START_REF] Menier | CFD validation and adaptivity for viscous flow simulations[END_REF] For the sake of clarity, we rewrite the linearized system in a more compact form:

A n δW n = R n (2)
where

A n = |C| δt n I - ∂R n ∂W and δW n = W n+1 -W n .
During an implicit simulation, this linearized system (2) is being solved at each flow solver iteration. Thus, accelerating the Newton method can dramatically decrease the total CPU time of the whole simulation. Moreover, an unsufficient convergence of the Newton method can cause a bad convergence of the whole simulation. We now describe how the linear system is solved using the symmetric Gauss-Seidel (SGS) relaxation, and then we present the implementation of multigrid methods which improve both the convergence speed and the robustness.

I.3. Solving the Linear System Using the SGS Relaxation

The linearized system (2) is solved after each time advancing using the approach based on Lower-Upper Symmetric Gauss-Seidel (LU-SGS) implicit solver initially introduced by Jameson [START_REF] Jameson | Lower-Upper implicit schemes with multiple grids for the Euler equations[END_REF] and fully developed by Sharov et al. and Luo et al. [START_REF] Luo | A fast, matrix-free implicit method for compressible flows on unstructured grids[END_REF][START_REF] Luo | An accurate, fast, matrix-free implicit method for computing unsteady flows on unstructured grids[END_REF][START_REF] Sharov | Reordering of hybrid unstructured grids for Lower-Upper Symmetric Gauss-Seidel computations[END_REF][START_REF] Sharov | Implementation of unstructured grid GMRES+LU-SGS method on shared-memory, cache-based parallel computers[END_REF] The Newton method we use is the symmetric Gauss-Seidel (SGS) relaxation and it is iterated until the residual of the system is reduced by a chosen order of magnitude (usually 0.01) with respect to the initial residual.

We now briefly describe the SGS relaxation. We start by decomposing A into a strictly lower triangular component L, a diagonal component D, and a strictly upper triangular component U:

A = L + D + U .
The linear system (2) can be rewritten :

(D + L)D -1 (D + U) δW n = R n + (LD -1 U) δW n ,
and approximated by the following system:

(D + L)D -1 (D + U) δW n = R n .
Using an SGS relaxation to converge the Newton method consists in performing k max SGS sub-iterations using forward and backward sweeps:

(D + L) δW k+1/2 = R -U W k (D + U) δW k+1 = R -L W k+1/2 .
which can be rewritten point wise:

δW k+1/2 i = D -1 ii R i - j∈L(i) L ij δW k+1/2 j - j∈U (i) U ij δW k j δW k+1 i = D -1 ii R i - j∈U (i) U ij δW k+1 j - j∈L(i) L ij δW k+1/2 j
. k max SGS sub-iterations are performed, unless fewer iterations are needed to decrease the residual of the system by the desired order of magnitude. This section has described how the linear system (2) is solved using an SGS relaxation in the case of a single-grid computation. We now present how we implemented multigrid methods, which use coarser meshes to accelerate and robustify the convergence of the Newton method.

I.4. Multigrid Methods for Accelerating the Convergence of the Newton Method

This section describes the multigrid procedures for converging the Newton method. Multigrid methods require a sequence of N meshes :

H h , H 2h , H 4h , . . . , H 2N h ,
where H h is the initial (and finest) mesh, and (H 2ih ) (i=1,N ) are the coarsened versions of H h generated as previously explained in Section I.1. For the sake of clarity, we rewrite the linear system (3) :

A h δu h = F h . (3)
where A h is the matrix of the linearized system built on H h and F h is the right hand side (RHS). The residual of ( 3) is given by

r h = A h δu h -F h .
In the case of a single-grid computation, we have seen that k max (or less) SGS sub-iterations are performed on H h in order to reduce r h by a desired order of magnitude. In the case of a multigrid computation, n max multigrid cycles are performed. One multigrid cycle consists in (i) performing one SGS sub-iteration on H h , (ii) computing a correction using the coarser meshes, and (iii) adding the computed correction to the solution on H h . The way this correction is computed depends on the number of coarse meshes involved and on the type of the multigrid cycle used.

Although one multigrid cycle is more costly in terms of CPU (than one single-grid sub-iteration), the number of cycles required to reach the targeted residual is expected to be smaller, thanks to the corrections. Note that the smaller the number of vertices of the coarsest mesh is, the quicker the correction is computed. Moreover, coarser meshes have a strong smoothing property, which increases the robustness, i.e. using a multigrid procedure makes it possible to reduce the residual by some orders of magnitude that could not be reached using one single mesh. To summarize, we use a multigrid cycle to compute a correction at each sub-iteration of the Newton method, in order to (i) increase the convergence speed, while (ii) improving the robustness.

We now describe the three different types of multigrid cycles we use to compute the corrections: the V-cycle , the W-cycle and the F-cycle. We start by explaining the case of the two-grid V-cycle, which only requires one coarser mesh. Then, the three types of cycles are introduced in the general case of N meshes.

I.4.1. Basic Two-grid V-cycle

One V-CYCLE ... The two-grid V-cycle requires a mesh H h and a coarser mesh H 2h . We suppose that N time iterations were performed by the flow solver. Let A h δu h = F h be the linearized system obtained after the N -th time iteration. In order to accelerate the convergence of the Newton method for solving this linear system, a given number of multigrid cycles can be performed. The bigrid V-cycle (see Figure 4) consists in computing a correction by performing several SGS iterations on H 2h .

A h , the matrix of the linearized system, was built on H h as explained in Section I.3. A 2h was built in a similar way on H 2h after S h , the solution obtained on H h after N time iterations, was linearly interpolated to H 2h . Starting from an initial solution δu 0 h , a pre-smoothing is performed on H h , i.e. one SGS iteration. Note that when the multigrid cycle is the first of the current time iteration, δu 0 h is set to 0, and otherwise δu 0 h is the solution of the previous cycle. Le δu 1 h be the solution obtained after the pre-smoothing, the residual is computed:

r h = A h δu 1 h -F h and is restricted to H 2h .
The restriction operator R h→2h first consists in locating each vertex P h of H h in H 2h , i.e. identifying the element K 2h = (P 2h (i)) i=0,3 of H 2h containing P h . Then, the restricted residual is summed to the vertices of K 2h :

R h→2h (r h )(P 2h (i))+ = β i × r h (P h ) for i = 0, 3 ,
where β i is the barycentric coordinate of P h in K 2h associated to P 2h (i).

The correction c 2h is then computed on H 2h by using R h→2h (r h ) as the source term. The initial correction is set to 0 : c 0 2h = 0 and a given number of SGS iterations is performed:

A 2h c 0 2h = R h→2h (r h ) (4)
Then, the resulting correction c 1 2h is linearly interpolated to H h and added to the solution:

δu 2 h = δu 1 h + I 2h→h (c 1 2h
) . So, at each time iteration of the flow solver, the corrections added after each multigrid cycle are expected to improve the convergence of the Newton method, and thus to improve the convergence of the whole simulation. The number of multigrid cycles required to reach the targeted residual is expected to be smaller than the required number of SGS iterations on H h in the single-grid case. The fastest convergence in terms of the number of iterations can be reached using an ideal bigrid V-cycle.

Ideal Bigrid V-cycle. The ideal bigrid V-cycle consists in performing as many SGS iterations on H 2h as are necessary to reduce the residual of the linear system (4) as much as possible, and thus obtain a fully converged correction c 1 2h . This ideal bigrid cycle is obviously too costly in terms of CPU to be used during a real-life simulation, but since the linear system is converged to its maximum on H 2h , it provides the best correction that can be obtained using a multigrid cycle. The number of iterations needed by an ideal bigrid cycle to converge the linear system on the finest mesh H h can thus be targeted when using another multigrid cycle (using more mesh levels). In other words, a "good" multigrid cycle aims at requiring as few iterations as the ideal bigrid cycle to decrease the residual on H h , while being less costly in terms of CPU thanks to the use of more coarser mesh levels.

I.4.2. N-grid V-Cycle, W-Cycle and F-Cycle

The N-grid V-Cycle is simply the extension of the bigrid V-cycle to N grids. Only one SGS iteration is performed on H 2h and the residual r 1 2h = A 2h c 1 2h -r 2h is computed and restricted to H 4h . R 2h→4h is then used as the source term to compute a correction C 1 4h on H 4h . This is how a correction is computed on each coarse mesh. Once on the coarsest mesh H 2N h , not one but several SGS iterations are used to compute C 1 2N h , which is not costly in terms of CPU due to the low number of vertices. Then, the correction of the coarsest mesh is interpolated and added to the correction of the second coarsest mesh and so on. In the end, the final correction containing all the contributions of the coarser meshes is interpolated on the finest mesh and added to δu 1 h . A post-smoothing (i.e. one SGS iteration) can be performed on each level i after I 2(i+1)h→2ih (c 2(i+1)h ) has been added to c 2ih .

Other types of multigrid cycles may be used, such as the V-, the W-and the F-cycle. The structures of these three cycles is depicted in Figure 6 for the case of four grids, and a 5-grid W-cycle is depicted in Figure 7.

I.5. Uniform Multigrid Results

Two different comparisons were made in order to measure the impact of multigrid methods. We compared (i) the convergence of the linear system at a given time iteration, and (ii) the convergence of the residual of the whole simulation from the first to the last iteration in time. Three cycles were used (V-cycle, W-cycle and F-cycle), as well as several numbers of coarse levels.

I.5.1. Resolution of the Linear System

The convergence of the linear system at a given time iteration is compared. A solution was previously "almost" converged on the finest mesh by performing N time iterations using any method (single-grid or multigrid) and using an adequate CFL law. The evolution of the residual of the resolution of the linear system obained after the N -th iteration is then compared for (i) the single-grid method, (ii) the ideal bigrid, and (iii) V-, W-and F-cycles using various numbers of mesh levels. A high CFL is prescribed in order to evaluate the robustness of each method.

The first example considers a 2D transonic NACA 0012 airfoil. The Mach number is M = 0.8 and the angle of attack α = 1.25. A solution (see Figure 5) is precomputed on the finest mesh using 120 time iterations at CFL max = 10. Then, the resulting solution is used as a restart solution and one iteration at CFL = 1000 is performed using the different methods. Views of the four meshes used during multigrid simulations are depicted in Figure 8. Figure 9 presents the convergence rates obtained using one single-grid, an ideal bigrid, and 3 V-cycles (3, 4 and 5 meshes). All the multigrid methods manage to decrease the initial residual by twelve orders of magnitude, while in the same CPU time interval, the single-grid computation fails to decrease it by one order due to the high CFL. As expected, the ideal bigrid shows the fastest convergence in terms of the number of iterations but is also the slowest method in terms of CPU. Figure 10 presents a comparison of the three different 4 grid cycles used (V, W and F). Although both the W-cycle and the F-cycle are really close to the ideal bigrid in terms of the number of iterations, they are slower than the V-cycle. To summarize, the fastest convergence for the transonic NACA is the 4-grid V-cycle in terms of CPU, and the 4-grid F-cycle in terms of the number of iterations.

I.5.2. Impact on the Whole Simulation

The evolution of the residual after each time iteration is compared for the single-grid method and several multigrid cycles. Starting from a uniform solution, the number of time iterations needed to reach a targeted Transonic NACA: Comparison of the V-cycles for the convergence of the Newton method after 120 time iterations. residual is compared for three examples : a 2D transonic NACA 0012, a 3D subsonic NACA 0012 and a 3D transonic wing body tails (WBT) configuration.

2D transonic NACA 0012. Section I.5.1 presented how multigrid methods help to converge the Newton method at a given time iteration for the 2D transonic NACA 0012. This time, the convergence of the residual in terms of the number of flow solver iterations is considered for the same test case. Figure 11 presents a comparison between the single-grid and the multigrid methods. Multigrid methods improve the convergence rate in terms of both the number of iterations and wall clock time. As regards the number of iterations, the best convergence rate is obtained in the single-grid case by performing 40 SGS sub-iterations, and in the multigrid case using a 3-grid V-cycle. The fastest methods in terms of CPU are the 3-grid and 4-grid V-cycles. 3D subsonic NACA 0012. This example considers a subsonic flow around an extruded NACA 0012 geometry. The meshes used for the multigrid computation are depicted in Figure 12 and the final solution in Figure 14. This example is interesting because the convergence of the residual of the whole simulation greatly depends on the Newton method. As shown in Figure 13, no fewer than 25 SGS sub-iterations are required in the single-grid case to reduce the residual of the whole simulation by the desired order of magnitude (10 -9 is the target). Using a multigrid method, however, Figure I.5.2 shows that only one V-cycle is enough, and that performing two cycles is enough to obtain an optimal residual convergence in terms of the number of iterations, i.e. performing more than two cycles does not help to increase the convergence rate. Figure I.5.2 also presents a comparison between the most efficient single-grid method (i.e. 25 SGS-sub-iterations), and the multigrid. As concerns the number of iterations, the residual convergence of the optimal single-grid method (25 SGS sub-iterations) and the optimal multigrid method (2 V-cycles) are identical. The wall clock time, however, drops from 6m50s in the single grid case to 1m54s using one V-cycle. 3D transonic WBT configuration. A transonic flow is computed over a 3D wing body tails (WBT) configuration geometry depicted in Figure 17, which includes a wing, a body and two tails (horizontal and vertical). The prescribed Mach number is M = 0.8 and the angle of attack α = 1 degree. The convergence of the residual is compared for a single-grid computation and a 3-grid V-cycle. The three meshes used for the computation are presented in Figure 16.

II. Full Multigrid (FMG) Coupled with Adaptivity

II.1. Full Multigrid Validation

This section presents the full multigrid (FMG) algorithm and its validation. Figure 18. Solution computed on the finest WBT mesh using a 3-grid V-cycle.

II.1.1. Description of the FMG Method

The full multigrid algorithm is an iterative process (see Figure 23). At each stage, a solution is computed on a mesh, whose complexity of the mesh is increased from one stage to the next. The coarsest mesh is used for the first stage and the finest mesh is used for the final stage. At each stage i (i >= 2), the coarser meshes from the previous stages are used to run a multigrid simulation. At the end of each stage, the solution is linearly interpolated to the next (finer) mesh and then used as a restart solution by the flow solver. For the sake of clarity, we describe the FMG algorithm for the case of a sequence of four meshes:

H h , H 2h , H 4h , H 8h ,
where H h is the finest mesh and H 8h the coarsest. A full multigrid algorithm using these four meshes consists in computing a solution on each mesh, as follows:

1. On H 8h : starting from a uniform solution S 0 8h , a solution S 8h is computed using a single-grid method. S 8h is then interpolated to H 4h .

2. On H 4h : the interpolated of S 8h is used as a restart solution by the flow solver : S 0 4h = I 8h→4h (S 8h ). A two-grid multigrid simulation is then performed on H 4h using H 8h as the coarse mesh. S 4h is then interpolated to H 2h .

On H

2h : S 0 2h = I 4h→2h (S 4h
) is used as a restart solution and a 3-grid multigrid simulation is performed on H 2h using H 4h and H 8h as coarser meshes.

On H

h : S 0 h = I 2h→h (S 2h
) and a 4-grid multigrid simulation is performed.

This algorithm is naturally extended to N meshes. The FMG theory [START_REF] Carré | On the Application of FMG to Variational Approximation of Flow Problems[END_REF] states that if the solution is fully converged at stage 1, then it is sufficient to converge the solution by one order of magnitude at stages 2, 3, etc. in order to achieve the global convergence on the finest mesh. In other words, fully converging the solution at every stage would not improve the residual on the finest mesh (and would be more CPU comsuming). In our approach, we chose to converge the solution by two orders at stages 2, 3, etc. instead of one order, to be sure to achieve the global convergence on adapted meshes.

II.1.2. FMG Validation

The FMG method is validated using two examples introduced in Section I.5.2: a 2D transonic NACA and a 3D wing body tails (WBT) configuration. We compared the convergence rate of the residual in terms of wall clock time, for the FMG algorithm and a 4 grid V-cycle simulation on the finest mesh, which can be seen as the last stage of the FMG algorithm, but starting from a uniform solution. Results are presented in Figure 20 for the 2D NACA and in Figure 24 for the 3D transonic WBT configuration. In both cases, a faster convergence is obtained using the FMG method. 3D transonic WBT configuration: Comparison of convergence rates with respect to wall clock time. A faster convergence is obtained using an FMG algorithm.

II.2. Steady Mesh Adaptation

Mesh adaptation provides a way to control the accuracy of the numerical solution by modifying the domain discretization according to size and directional constraints. When dealing with real life flow problems, Hessian-based unstructured mesh adaptation has already proved its efficiency to improve the tradeoff between the accuracy of the solution and the number of degrees of freedom (i.e. the complexity of the problem). [START_REF] Castro-Díaz | Anisotropic Unstructured Mesh Adaptation for Flow Simulations[END_REF][START_REF] Frey | Anisotropic mesh adaptation for CFD computations[END_REF][START_REF] Gruau | 3D tetrahedral, unstructured and anisotropic mesh generation with adaptation to natural and multidomain metric[END_REF][START_REF] Li | 3D anisotropic mesh adaptation by mesh modification[END_REF] In addition, as a large number of physical phenomena are anisotropic by nature, anisotropic mesh adaptation improves this ratio even more.

The classical steady mesh adaptation scheme is a fixed-point algorithm. Starting from an initial mesh and solution (H 0 , S 0 0 ), both of them are converged to an optimal state based on metric tensor fields and on the concept of unit-mesh.

The classical steady mesh adaptation loop. The main steps of the algorithm are depicted in Figure 22. Given (H i , S i ), a metric tensor M i is computed at each vertex of the mesh H i , according to a given error estimation of the solution (Hessian-based, for instance). M i contains information on optimal sizes and

(H i , S i ) (H i , M i ) (H i , S 0 i ) (H 0 , S 0 0 ) (H i+1 , S i , H i ) S i M i H i+1 S 0 i+1 Compute Solution Compute Metric

Generate Mesh

Interpolate Solution directions. This information given by the metric tensor field M i is then used by the adaptive anisotropic remesher (described in Section I.1) to generate a new mesh H i+1 . Then, S i is interpolated on H i+1 : we obtain S 0 i+1 which is then used as a restart solution by the solver for the next stage of the mesh adaptation loop. This loop is repeated until convergence, i.e. until there is not too much variation in both the mesh and the solution from one stage to the next.

II.3. Coupling FMG and Mesh Adaptation

There are many similarities between the FMG algorithm and the classical mesh adaptation loop. In both cases, (i) we start from an initial coarse mesh, (ii) the complexity of the current mesh is increased at each stage, and (iii) the solution on the current mesh is linearly interpolated to the next mesh and then used as a restart solution by the flow solver. We start by describing the coupling of the two algorithms, then we explain how both of them can benefit from this coupling.

Description of the coupling. Coupling a FMG algorithm with adaptivity consists in modifying the solution computation step in the classical mesh adaptation loop (see Figure 22). Instead of a single-grid computation, a i-grid multigrid computation is performed at stage i, using the meshes previously adapted as coarser meshes.

Improving the FMG algorithm using a coupling. There are two reasons why the FMG algorithm can benefit from this coupling. First, it is based on uniform meshes which are generated previous to the computation. These non-adapted meshes do not take into account the characteristics of the physical phenomena, which are often anisotropic and located in small areas of the computational domain. Thus, FMG uses meshes which are not optimal in terms of both sizes and directions. As explained in Section II.2, mesh adaptation has proved its efficiency in reducing computational time while improving the accuracy of the solution. The second reason is that the anisotropic meshes generated using mesh adaptation are particularly well-suited to a multigrid computation.

Improving the mesh adaptation loop using a coupling. An example of how mesh adaptation can benefit from multigrid is presented in Figure 24. Two simulations of a transonic flow over a Falcon business jet geometry are considered. First, a classical FMG algorithm is performed using four mesh levels H h , H 2h , H 4h and H 8h . The resulting residual convergence is compared to a second simulation, which consists in replacing every multigrid simulation of the FMG algorithm by a single-grid computation. In other words, this second computation is very similar to the classical mesh adaptation loop: mesh complexity is increased at each stage but no multigrid strategy is used. The residual convergences are identical for the two first stages, but a major difference appears during the third: the single-grid computation fails to converge on H 2h , while the multigrid procedures used on H 2h and then on H h ensure clean convergences on both meshes. So, in addition to a faster convergence in terms of wall clock time, this comparison shows that multigrid procedures can improve the robustness of the mesh adaptation strategy.

II.4. Results

Two cases were considered to validate the adaptive FMG algorithm: a 3D subsonic NACA 0012 airfoil and a 3D transonic WBT configuration, which were both introduced in Section I.5.2. For each case, two simulations were performed:

• A classical mesh adaptation process: the prescribed mesh complexity was increased at each stage, and monogrid simulations were run. As is it usually done, the residual of the solution was fully converged at each stage.

• An adaptive FMG algorithm: the same mesh complexities were prescribed. The residual of the solution was fully converged for the lowest mesh complexity (i.e. stage 1), and then reduced by two orders of magnitude for the other complexities.

In both cases, the adaptive multigrid algorithm showed a significant reduction of the total wall clock time of the simulation. It was ensured that both methods converged to same final solution. To do so, one more stage of the classical mesh adaptation loop was performed using a higher mesh complexity, and the resulting couple mesh/solution was then used as a reference solution to compute spatial errors. Both algorithms showed the same mesh convergence.

II.4.1. 3D Subsonic NACA 0012

The first mesh adaptation considers the 3D subsonic NACA 0012 case that was introduced in Section I.5.2. Six stages of the classical mesh adaptation loop and of the adaptive FMG alorithm were performed, using mesh complexities leading approximately to the following numbers of vertices: 8k, 16k, 32k, 64k, 128k, 256k .

For each mesh complexity, three sub-iterations in the adaptation loop were performed. The residual of the solution was fully converged to 10 -9 at each stage of the classical adaptation. A slope limiter was used in order to avoid spurious oscillations. [START_REF] Alauzet | High Order Sonic Boom Modeling by Adaptive Methods[END_REF] A freeze of this limiter is activated in case the limiter itself is oscillating. During the adaptive FMG algorithm, the residual was fully converged at stage 1 (which does not differ from the classical algorithm), and then it was reduced by two orders of magnitude at stages 2, 3, etc. The residual convergence of both simulations in terms of wall clock time is presented in Figure 25. The total wall clock time of the simulation is dramatically improved: 13min9s for the adaptive FMG method, and 2h25min for the classical adaptation algorithm. As shown in Figure 26, the mesh convergence observed for both methods are similar. The reference couple mesh/solution was computed using one more step in the classical adaptation loop at a mesh complexity leading to ∼ 512k vertices. The final adapted mesh and the solution are depicted in Figures 27 and28 3D subsonic NACA 0012: comparison of the residual convergence in terms of wall clock time. Left: whole simulation. Right: close-up view of the first stages. Note that each stage corresponds to a mesh complexity, and that three sub-iterations were performed for each one of them.

II.4.2. 3D Transonic Wing Body Tails Configuration

The second example considers the 3D transonic WBT configuration that was also introduced in Section I.5.2. Four stages of the classical mesh adaptation loop and of the adaptive FMG algorithm were performed. The prescribed mesh complexities (corresponding to each stage) lead approximately to the following numbers of vertices:

140k, 210k, 340k, 620k .

At each stage (i.e. mesh complexity), five sub-iterations in the adaptation loop were performed. A comparison of the residual convergence of both methods in terms of wall clock time is presented in Figure 29. The total wall clock time of the simulation is reduced from 1d3h57m for the classical mesh adaptation loop, to 2h53m for the adaptive FMG algorithm. As shown in Figure 30, converging the residual by two orders of magnitude

III. Conclusion

The implementation of an implicit multigrid procedure has been described and its validation study has been carried out. A significant improvement of both the convergence speed and the robustness has been observed in the non-adapted case.

The coupling of an FMG algorithm with mesh adaptation has been detailed and compared to the classical mesh adaptation loop for a 3D subsonic and a 3D transonic case. Converging the solution by two orders of magnitude in the adaptive FMG algorithm significantly reduced the total wall clock time of the simulation, while achieving the same global convergence on adapted meshes.

We are working on mesh adaptation strategies for boundary layers, which we hope will allow to extend the adaptive FMG algorithm to turbulent viscous flow simulations.

Figure 1 .

 1 Figure 1. Example of mesh coarsening. Left: initial mesh. Middle: Isotropic coarsening. Right: Anisotropic coarsening.

Figure 2 .

 2 Figure 2. Comparison between nested and non-nested meshes. The coarse mesh H 2h (in black) and the finer mesh H h (in grey) are juxtaposed for visualization purposes. Left: H h and H 2h are nested. Right: Non-nested meshes.

Figure 3 .

 3 Figure 3. Close-up views of three nested meshes of a 2D scramjet. Left: the initial coarse mesh. Middle: all the edges of the mesh were refined. Right: only some edges were refined. Both nested meshes (middle and right) contain elements of bad quality due to the constraints from the initial discretization.

Figure 4 .

 4 Figure 4. Bigrid V-cycle Figure 5. 2D transonic NACA 0012: pressure.

Figure 6 .Figure 7 .

 67 Figure 6. Four-grid methods : V-cycle, W-cycle and F-cycle (• : 1 smoothing SGS iteration, • : Several SGS iterations)

Figure 8 .

 8 Figure 8. Close-up views of the four meshes used during the multigrid computations of the 2D transonic NACA. Number of vertices, from left to right: 29024, 7379, 2499 and 1305.

  Figure 9.Transonic NACA: Comparison of the V-cycles for the convergence of the Newton method after 120 time iterations.

Figure 10 .

 10 Figure 10. Transonic NACA: Comparison of the 3 cycles (4 grids) for the convergence of the Newton method after 120 time iterations.

Figure 11 .

 11 Figure 11. 2D transonic NACA: Convergence of the residual of the whole simulation. Left: number of time iterations. Right: wall clock time (sec).

Figure 12 .

 12 Figure 12. 3D subsonic NACA 0012 : Cuts in the volumes of the four meshes used during the multigrid computations.

Figure 13 .

 13 Figure 13. 3D subsonic NACA 0012: At least 25 SGS sub-iterations are required in the singlegrid case to converge the residual of the whole simulation.

Figure 14 .

 14 Figure 14. 3D subsonic NACA 0012: Solution (velocity).

Figure 15 .

 15 Figure 15. 3D subsonic NACA 0012 : Convergence of the residual of the whole simulation. Left: number of time iterations. Right: wall clock time (sec).

Figure 16 .

 16 Figure 16. 3D WBT configuration : Close-up views of the three meshes used for multigrid simulations.

Figure 17 .

 17 Figure 17. 3D wing body tails (WBT) configuration.

Figure 19 .

 19 Figure 19. 3D WBT : Residual convergence in terms of the number of solver iterations (left) and wall clock time (right). Comparison between a monogrid computation and a 3-grid V-cycle.

Figure 20 .

 20 Figure 20. 2D transonic NACA : Comparison of convergence rates with respect to wall clock time. A faster convergence is obtained using an FMG algorithm.

Figure 21 .

 21 Figure 21.3D transonic WBT configuration: Comparison of convergence rates with respect to wall clock time. A faster convergence is obtained using an FMG algorithm.

Figure 22 .Figure 23 .

 2223 Figure 22. The classical mesh adaptation loop

Figure 24 .

 24 Figure24. Transonic Falcon: residual convergence in terms of time iterations. Blue: classical FMG algorithm. Green: multigrid computations were replaced by single-grid computations in the FMG algorithm.

  .

  Figure 25.3D subsonic NACA 0012: comparison of the residual convergence in terms of wall clock time. Left: whole simulation. Right: close-up view of the first stages. Note that each stage corresponds to a mesh complexity, and that three sub-iterations were performed for each one of them.

Figure 26 .

 26 Figure 26. 3D subsonic NACA 0012: mesh convergence to the reference solution.

Figure 27 .

 27 Figure 27. 3D subsonic NACA 0012: Cut in the volume of the final adapted mesh.

Figure 28 .Figure 29 .

 2829 Figure 28. 3D subsonic NACA 0012 : velocity isovalues.

Figure 30 .

 30 Figure 30. 3D transonic WBT: Mesh convergence to the reference solution.Figure31. Transonic WBT configuration: cut in the trailing vortices region of the final adapted mesh.

  Figure 30. 3D transonic WBT: Mesh convergence to the reference solution.Figure31. Transonic WBT configuration: cut in the trailing vortices region of the final adapted mesh.

Figure 32 .

 32 Figure 32. Transonic WBT configuration: pressure on the wing and corresponding adapted surface mesh.

Figure 33 .

 33 Figure 33. Transonic WBT configuration: velocity isovalues.

of 20 American Institute of Aeronautics and Astronautics

IV. Acknowledgements

Thanks to Dave Marcum for providing the initial mesh of his 3D WBT configuration and to Richard James for the proofreading. This work was partially funded by the Airbus Group Foundation.