N Barral

G Olivier

F Alauzet

Time-Accurate Anisotropic Mesh Adaptation for Three-Dimensional Time-Dependent Problems with Body-Fitted Moving Geometries

Keywords: Metric-based mesh adaptation, anisotropy, unsteady, moving geometry, moving mesh, connectivity change

Anisotropic metric-based mesh adaptation has proved its efficiency to reduce the CPU time of steady and unsteady simulations while improving their accuracy. However, its extension to time-dependent problems with body-fitted moving geometries is far from straightforward. This paper establishes a well-founded framework for multiscale mesh adaptation of unsteady problems with moving boundaries. This framework is based on a novel space-time analysis of the interpolation error, within the continuous mesh theory. An optimal metric field, called ALE metric field, is derived, which takes into account the movement of the mesh during the adaptation. Based on this analysis, the global fixed-point adaptation algorithm for time-dependent simulations is extended to moving boundary problems, within the range of body-fitted moving meshes and ALE simulations. Finally, three dimensional adaptive simulations with moving boundaries are presented to validate the proposed approach.

Introduction

Simulating complex moving geometries evolving in unsteady flows in three dimensions, which is more and more required by industry, still remains a challenge because it is very time consuming. To reduce the CPU time of numerical simulations while preserving their accuracy, anisotropic metric-based mesh adaptation has already proved its efficiency for steady problems [START_REF] Alauzet | High order sonic boom modeling by adaptive methods[END_REF][START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF][START_REF] Loseille | Achievement of global second-order mesh convergence for discontinuous flows with adapted unstructured meshes[END_REF][START_REF] Pain | Tetrahedral mesh optimisation and adaptivity for steady-state and transient finite element calculations[END_REF], and appears as a promising way to reduce the complexity of such simulations. However, its extension to the unsteady case with body-fitted moving geometries is not straightforward. Indeed, time-dependent simulations combine the difficulties arising from unsteadiness and geometrical complexity: global time-step driven by the mesh smallest height, evolution of physical phenomena in the whole domain, solution interpolation spoiling, and also three-dimensional meshing and remeshing issues with an imposed discretized surface. The introduction of body-fitted moving geometries in this process even raises new difficulties, due to the handling of the mesh movement and the deterioration of its quality, the specific numerical schemes imposed by moving geometries as well as fluid/structure coupling and contact handling.

There exist two main branches of mesh adaptation for unsteady problems: r-adaptation and adaptive remeshing. In the case of r-adaptation also called adaptive moving mesh methods, the mesh is continuously moved to an adapted configuration. The movement of the mesh vertices is governed by an extra equation that is strongly coupled to the underlying physics equations and that is solved at the same time. Since the mesh movement is explicitly taken into account in the equations, there are no spoiling interpolation steps. Several approaches can be found in the literature, that differ by the adaptation criteria, the way the adaptation governs the mesh movement and the way the mesh movement is coupled with the physics equations. Recent works on the topic include Moving Mesh PDEs [START_REF] Arpaia | Mesh adaptation by continuous deformation. basics: accuracy, efficiency, well balancedness[END_REF][START_REF] Huang | Adaptive mesh movement -the MMPDE approach and its applications[END_REF], and Lagrangian methods with ALE rezoning [START_REF] Loubère | ReALE: A reconnection-based arbitrary-Lagrangian-Eulerian method[END_REF]. A Monge-Ampère equation has also been used to drive the adaptation and has shown interesting results in 3D in [START_REF] Browne | Fast three dimensional r-adaptive mesh redistribution[END_REF][START_REF] Chacón | Robust, multidimensional mesh-motion based on Monge-Kantorovich equidistribution[END_REF]. However, the solution of the extra moving mesh equation is often costly in terms of CPU in 3D, and the optimality of the adapted mesh is not achieved temporally as the same number of degrees of freedom is used through the whole simulation. Moreover, it is unsure how these methods can be coupled to moving boundary problems.

The other branch consists in performing frequent adaptive remeshings. Three different approaches can be distinguished in the literature. First [START_REF] Löhner | Three-dimensional fluid-structure interaction using a finite element solver and adaptive remeshing[END_REF][START_REF] Löhner | Adaptive h-refinement on 3D unstructured grids for transient problems[END_REF][START_REF] Rausch | Spatial adaptation procedures on tetrahedral meshes for unsteady aerodynamic flow calculations[END_REF][START_REF] Speares | A 3D unstructured mesh adaptation algorithm for time-dependent shock-dominated problems[END_REF], an isotropic mesh is adapted frequently in order to maintain the solution within refined regions and introduce a safety area around critical regions. Another approach is to use an unsteady mesh adaptation algorithm [START_REF] Bourgault | Anisotropic error estimates and space adaptivity for a semidiscrete finite element approximation of the transient transport equation[END_REF][START_REF] Carpio | An anisotropic, fully adaptive algorithm for the solution of convection-dominated equations with semi-Lagrangian schemes[END_REF][START_REF] De Sampaio | Petrov-Galerkin solutions of the incompressible Navier-Stokes equations in primitive variables with adaptive remeshing[END_REF][START_REF] Porta | A space-time adaptation scheme for unsteady shallow water problems[END_REF][START_REF] Wu | Error estimation and adaptivity in Navier-Stokes incompressible flows[END_REF] based on local or global remeshing techniques and the estimation of the error every n flow solver iterations. If the error is greater than a prescribed threshold, the mesh is re-adapted. In [START_REF] Hay | Computation of free-surface flows with local mesh adaptation[END_REF], the authors combine both approaches mentioned above. Local adaptive remeshing enabling the construction of anisotropic meshes has been also considered. In this case [START_REF] Pain | Tetrahedral mesh optimisation and adaptivity for steady-state and transient finite element calculations[END_REF][START_REF] Remacle | Anisotropic adaptive simulation of transient flows using discontinuous Galerkin methods[END_REF], the mesh is frequently adapted in order to guarantee that the solution always evolves in refined regions. All these approaches involve a large number of mesh adaptations leading to a large increase in CPU time due to the generation of many meshes 1 , and the numerous IOs and database initializations 2 . This large number of mesh adaptations also introduces unquantified errors due to the transfer of the solution from the old mesh to the new one. In the context of conservative equations, this can result in a time-shift of the solution if a non-conservative solution transfer is considered [START_REF] Alauzet | A parallel matrix-free conservative solution interpolation on unstructured tetrahedral meshes[END_REF]. Moreover, these methods can result in a time shift between the mesh and the solution, since the adapted mesh is generated at one time-step, from the solution at that time-step, and is kept for the next few time-steps. Finally, none of them consider the inherent non-linear nature of the mesh adaptation problem: the convergence of the mesh adaptation process is never addressed and therefore obtaining the optimal mesh cannot be expected.

A first answer to these issues has already been proposed [START_REF] Alauzet | 3D transient fixed point mesh adaptation for time-dependent problems: Application to CFD simulations[END_REF][START_REF] Ducrot | Anisotropic level set adaptation for accurate interface capturing[END_REF][START_REF] Guégan | An L ∞ -L p mesh adaptive method for computing unsteady bi-fluid flows[END_REF][START_REF] Habashi | Improved transient-fixed-point mesh adaptation using orthogonality-preserving metric intersection[END_REF]. It is based on a control of the space-time interpolation error in L ∞ norm, the subdivision of the time interval into large sub-intervals on which the mesh is kept constant and on a local fixed-point algorithm to address the prediction of the solution and the convergence of the non-linear mesh adaptation problem. Since then, multiscale anisotropic mesh adaptation [START_REF] Alauzet | High order sonic boom modeling by adaptive methods[END_REF][START_REF] Loseille | Achievement of global second-order mesh convergence for discontinuous flows with adapted unstructured meshes[END_REF] (i.e., the control of the interpolation error in L p norm) and goal-oriented anisotropic mesh adaptation [START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF] (i.e., the control of a goal-oriented error estimate in L 1 norm) have proved to be a lot more efficient for steady CFD computations. Therefore, it seems relevant to extend these approaches to time-dependent simulations. Time-accurate multiscale and goal-oriented anisotropic mesh adaptation have been derived in [START_REF] Alauzet | Multi-scale anisotropic mesh adaptation for time-dependent problems[END_REF] and [START_REF] Belme | Time accurate anisotropic goal-oriented mesh adaptation for unsteady flows[END_REF], respectively, but for fixed meshes. They are based on a space-time interpolation error analysis and a global fixed-point mesh adaptation algorithm. The present paper extends the time-accurate multiscale anisotropic mesh adaptation of [START_REF] Alauzet | Multi-scale anisotropic mesh adaptation for time-dependent problems[END_REF] to the case of moving geometries with dynamic meshes.

Three classes of methods exist to deal with moving geometries: immersed/embedded methods, chimera/overset methods, and body-fitted methods. In this work, we select the last class. If very frequent remeshings are performed (each n flow solver time-steps), the movement of the geometries and thus of the meshes can be dealt with by the remeshings [START_REF] Compere | A mesh adaptation framework for dealing with large deforming meshes[END_REF][START_REF] Dobrzynski | Anisotropic Delaunay mesh adaptation for unsteady simulations[END_REF]. But, this leads to the issues stated above. In our case, the adaptive remeshings are much less frequent as only one adapted mesh is generated for each sub-interval, so the handling of the moving geometries is crucial, both from the purely moving mesh point of view and from the solver point of view. Indeed, algorithmically (to keep a simple mesh adaptation scheme) and theoretically (to derive the error analysis), it is fundamental to only remesh when it is required by the adaptation, i.e., for each sub-interval, and to never remesh because the mesh becomes too distorted due to the mesh deformation. Many works exist on mesh deformation strategies and Arbitrary-Lagrangian-Eulerian (ALE) numerical schemes [START_REF] Baum | A new ALE adaptive unstructured methodology for the simulation of moving bodies[END_REF][START_REF] Hassan | A method for time accurate turbulent compressible fluid flow simulation with moving boundary components employing local remeshing[END_REF][START_REF] Löhner | Adaptive embedded unstructured grid methods[END_REF][START_REF] Murman | Simulation of 6-DOF motion with cartesian method[END_REF][START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF], but generally they are not able to handle large geometry displacement without any remeshing throughout the simulation. Here, we consider the connectivity-change moving mesh algorithm proposed in [START_REF] Alauzet | A changing-topology moving mesh technique for large displacement[END_REF][START_REF] Barral | Two mesh deformation methods coupled with a changing-connectivity moving mesh method for CFD applications[END_REF] which proves to be very efficient to handle any 3D geometry displacements thanks to mesh optimizations based on connectivity changes that manage efficiently any shearing inside the mesh deformation.

Regarding adaptive strategies for moving mesh simulations, only a few attempts can be found in the literature among which [START_REF] Compere | A mesh adaptation framework for dealing with large deforming meshes[END_REF][START_REF] Hassan | A method for time accurate turbulent compressible fluid flow simulation with moving boundary components employing local remeshing[END_REF][START_REF] Löhner | Adaptive remeshing for transient problems[END_REF][START_REF] Löhner | Numerical simulation of a shock interaction with complex geometry three-dimensional structures using a new adaptive Hrefinement scheme on unstructured grids[END_REF][START_REF] Löhner | Three-dimensional store separation using a finite element solver and adaptive remeshing[END_REF][START_REF] Saksono | An adaptive remeshing strategy for flows with moving boundaries and fluid-structure interaction[END_REF]. As impressive as they can be, these results nevertheless still suffer from some of the weaknesses described above, i.e., mainly very frequent remeshing and spoiling interpolation stages. Moreover, there is no well-established framework for the space-time error estimate and for the consideration of the dynamic of the mesh. The contribution of this paper is to propose an extension of the time-accurate multiscale anisotropic mesh adaptation to moving-geometry problems, and to demonstrate numerically that three-dimensional moving mesh adaptation simulations can actually be run with this algorithm.

To this aim, we first (Section 2) recall the extension of the space-time error analysis of [START_REF] Belme | Time accurate anisotropic goal-oriented mesh adaptation for unsteady flows[END_REF] to time-accurate multiscale anisotropic mesh adaptation, which provides a control of the space-time interpolation error in L p norm and is thoroughly detailed in [START_REF] Alauzet | Multi-scale anisotropic mesh adaptation for time-dependent problems[END_REF]. Section 3 describes the global fixed-point mesh adaptation algorithm and its practical implementation. The extension of the space-time analysis and of the mesh adaptation algorithm to moving mesh simulations is then addressed in Sections 4 and 5. The contributions of this work are the following:

• propose an ALE metric field formulation which takes into account the displacement of the mesh. This ALE metric field formulation is validated on several analytic examples by analyzing the quality of the resulting meshes and by performing a convergence study of the error.

• extend the space-time interpolation error analysis to dynamic meshes using the ALE metric field formulation. The idea is to match instantaneous interpolation error using the ALE metric field in order to map the error onto the initial continuous mesh configuration, this removes the dependence of the continuous mesh in time.

• modify the connectivity-change moving mesh algorithm [START_REF] Alauzet | A changing-topology moving mesh technique for large displacement[END_REF] to take into account anisotropic adapted dynamic meshes. In particular, a dynamic metric field is considered to perform the mesh optimizations.

• demonstrate that body-fitted ALE numerical simulations involving large displacement of complex geometries and coupled with anisotropic mesh adaptation is now achievable in a very robust and efficient way. To this end, many numerical examples of 3D unsteady adaptive moving mesh simulations are presented and analyzed (Section 6).

Time-accurate multiscale anisotropic mesh adaptation for unsteady problems

This section first recalls the context of metric-based mesh adaptation, then a space-time interpolation error analysis is presented. Presenting the context and the time-accurate anisotropic mesh adaptation framework is essential to understand the analysis with dynamic meshes. Note that the results are only presented in dimension three.

Metric-based generation of anisotropic adapted meshes

Metric-based generation of anisotropic adapted meshes uses the notion of Riemannian metric space [START_REF] Frey | Mesh generation. Application to finite elements[END_REF][START_REF] Loseille | Continuous mesh framework. Part I: well-posed continuous interpolation error[END_REF][START_REF] Loseille | Continuous mesh framework. Part II: validations and applications[END_REF]. For a computational domain Ω ⊂ R 3 , a Riemannian metric space (M(x)) x∈Ω is a spatial metric field that defines at any point of Ω a metric tensor M(x), e.g. a 3 × 3 symmetric positive definite matrix. It is then possible for a mesh generator to work, i.e., to evaluate all geometric quantities, in this Riemannian metric space instead of working in the canonical Euclidean space. In a Riemannian metric space, the dot product is defined locally by metric tensor M: u, v M = u, Mv for (u, v) ∈ R 3 × R 3 . Thus, the computation of geometric quantities requires integral formulae to take into account the variation of the metric field. In that case, the length of edge ab is computed using the straight line parameterization γ(t) = a + t ab, where t ∈ [0, 1]:

M (ab) = 1 0 γ (t) M dt = 1 0 ab T M(a + t ab) ab dt ,
and the volume of element K is:

|K| M = K det M(x) dx .
The main idea of metric-based mesh adaptation, initially introduced in [START_REF] George | Creation of internal points in Voronoi's type method. Control and adaptation[END_REF], is to generate a unit mesh in the prescribed Riemannian metric space, e.g. a mesh of Ω ⊂ R 3 such that each edge has a unit length and each tetrahedron is regular (or equilateral) with respect to (M(x)) x∈Ω :

∀e, M (e) = 1 and ∀K, |K| M = √ 2 12
.

The resulting mesh in the canonical Euclidean space will be anisotropic and adapted. As it is not possible to tesselate R 3 with the regular tetrahedron, we seek for a mesh such that all its edges have a length close to unity and such that all the elements are almost regular in the considered Riemannian metric space. Following the duality introduced in [START_REF] Loseille | Continuous mesh framework. Part I: well-posed continuous interpolation error[END_REF][START_REF] Loseille | Continuous mesh framework. Part II: validations and applications[END_REF], a metric tensor (when it is defined point-wise) is a continuous element, and a Riemannian metric space (when it is defined all-over the domain) -also called metric field -is a continuous mesh. One or the other of these names will be used thereafter.

Space-time L p interpolation error analysis

The space-time error analysis is based on the continuous mesh framework presented in [START_REF] Loseille | Continuous mesh framework. Part I: well-posed continuous interpolation error[END_REF][START_REF] Loseille | Continuous mesh framework. Part II: validations and applications[END_REF]. However, the multiscale mesh adaptation described in these papers only controls spatial errors. In the context of time-dependent problems, temporal errors must be controlled as well. In [START_REF] Belme | Time accurate anisotropic goal-oriented mesh adaptation for unsteady flows[END_REF], a space-time error analysis is performed for a goal-oriented error estimate. Here, we recall the extension of this analysis to the control of the space-time interpolation error in L p norm. In what follows, we do not account for time discretization errors but we focus on a space-time analysis of the spatial errors in unsteady simulations. In other words, we seek for the optimal space-time mesh controlling the space-time spatial discretization error.

The following assumption is made: as an explicit time scheme is used for time advancing, the error in time is controlled by the error in space under the CFL condition. The proof of this assumption for the 1D scalar advection equation is provided in [START_REF] Alauzet | 3D transient fixed point mesh adaptation for time-dependent problems: Application to CFD simulations[END_REF] and the extension of this proof to the multi-dimensional linear advection problem is given in [START_REF] Olivier | Anisotropic metric-based mesh adaptation for unsteady CFD simulations involving moving geometries[END_REF]. However, there is no proof of this fact for multi-dimensional non-linear problems. This assumption is illustrated on Sod's shock tube problem in [START_REF] Olivier | Anisotropic metric-based mesh adaptation for unsteady CFD simulations involving moving geometries[END_REF] where, for explicit and implicit schemes with a CFL number less than one, we observe that no error in time appears while for implicit schemes with a CFL number larger than one large error in time appears. Therefore, if implicit schemes with arbitrary large time-steps are considered, an adaptive time-stepping method is required to control the error in time such as the analysis proposed in [START_REF] Coupez | Adaptive time-step with anisotropic meshing for incompressible flows[END_REF] for incompressible unsteady simulations and [START_REF] Olivier | Anisotropic metric-based mesh adaptation for unsteady CFD simulations involving moving geometries[END_REF] for the multi-dimensional linear advection problem. As far as the above assumption is true, the spatial interpolation error is a good measure of the total space-time error of the discretized unsteady system.

The space-time L p interpolation error analysis is thoroughly detailed in [START_REF] Alauzet | Multi-scale anisotropic mesh adaptation for time-dependent problems[END_REF]. However, it is absolutely essential to understand well the reasoning to understand its extension to the moving mesh case, which is why some of it is recalled in the present paper. Our goal is to solve an unsteady PDE which is set in the computational space-time domain Q = Ω × [0, T ] where Ω ⊂ R 3 is the spatial domain and T is the (positive) maximal time. An essential ingredient of our discretization and of our analysis is the element-wise linear interpolation operator. We define our working functional space which is the set of measurable functions that are continuous with square integrable gradients:

V = H 1 (Ω) ∩ C( Ω) , and V = H 1 {[0, T ]; V} .
We assume that Ω is covered by a finite-element partition made of simplicial elements denoted K. The mesh, denoted by H, is the set of the elements. Let us introduce the following approximation spaces:

V h = ϕ h ∈ V ϕ h|K is affine ∀K ∈ H , and V h = H 1 {[0, T ]; V h } ⊂ V.
Let Π h be the usual P 1 projector:

Π h : V → V h such that Π h ϕ(x) = ϕ(x), ∀ x vertex of H.
We extend it to time-dependent functions:

Π h : V → V h such that (Π h ϕ) (t) = Π h (ϕ(t)) , ∀ t ∈ [0, T ].
The problem of mesh adaptation consists in finding the space-time mesh H of Q that minimizes the space-time linear interpolation error u -Π h u in L p norm, for a given sensor function u and for a given number of space-time mesh vertices N st . The problem is thus stated in an a priori way:

Find H opt having N st space-time vertices such that E L p (H opt ) = min H u -Π h u L p (Ω h ×[0,T ]) . (1)
As it, this problem is a global combinatorial problem which turns out to be intractable in practice. This ill-posed problem can be reformulated in the continuous mesh framework [START_REF] Loseille | Continuous mesh framework. Part I: well-posed continuous interpolation error[END_REF][START_REF] Loseille | Continuous mesh framework. Part II: validations and applications[END_REF]. In this framework, we propose the following continuous model of a mesh. A continuous mesh M of a domain Ω is identified to a Riemannian metric space M = (M(x)) x∈Ω that, at each point x of Ω, prescribes a density, a set of anisotropy directions and the stretching along these directions, this information being locally contained in metric tensor M(x). The spatial size of the continuous mesh is given by its spatial complexity:

C(M) = Ω det M(x)
dx which is the continuous counterpart of the number of vertices. Continuous mesh M defines a class of equivalence of discrete meshes, which are all unit meshes with respect to M. We also define a continuous model of the linear interpolation operator Π h denoted π M . It is then possible to recast (1) into a well-posed continuous global optimization problem of finding the optimal space-time continuous mesh minimizing the space-time continuous interpolation error in L p norm, for a given sensor u and for a given space-time complexity N st :

Find M L p = (M L p (x, t)) (x,t)∈Q such that E L p (M L p ) = min M u -π M u L p (Ω×[0,T ]) , (2)
under the space-time complexity constraint:

C st (M) = T 0 τ(t) -1 Ω det(M(x, t)) dx dt = N st . (3)
The continuous mesh space-time complexity N st enables the user to control the level of accuracy of the mesh, and thus, to implicitly control the number of space-time vertices of the resulting discrete mesh. In its definition, τ(t) is the time-step used at time t of interval [0, T ] and det(M(x, t)) represents the continuous mesh local density. The continuous space-time error model can be written as follows (see [START_REF] Loseille | Continuous mesh framework. Part I: well-posed continuous interpolation error[END_REF]):

E L p (M) = T 0 Ω trace M -1 2 (x, t)|H u (x, t)|M -1 2 (x, t) p dx dt 1 p . ( 4
)
where H u is the Hessian of sensor u. To find the optimal space-time continuous mesh (see [START_REF] Alauzet | Multi-scale anisotropic mesh adaptation for time-dependent problems[END_REF][START_REF] Belme | Time accurate anisotropic goal-oriented mesh adaptation for unsteady flows[END_REF]), Problem (2-3) is solved in two steps: first a spatial minimization is done for a fixed t, then a temporal minimization is performed.

Spatial minimization for a fixed t. At time t, we seek for the optimal continuous mesh M L p (t) which minimizes the instantaneous interpolation error E L p in which the integral in time of Relations (3) and (4) was removed. Similarly to [START_REF] Loseille | Continuous mesh framework. Part I: well-posed continuous interpolation error[END_REF], solving the optimality conditions provides the optimal instantaneous continuous mesh in L p norm M L p (t) = (M L p (x, t)) x∈Ω at time t defined by:

M L p (x, t) = N(t) 2 3 Ω (det |H u (x, t)|) p 2p+3 dx -2 3 (det |H u (x, t)|) -1 2p+3 |H u (x, t)| . (5)
We set:

K(t) = Ω (det |H u (x, t)|) p 2p+3 dx.
The corresponding optimal instantaneous interpolation error in L p norm at the power p at time t is:

E L p (M L p (t)) = 3 p N(t) -2p 3 Ω (det |H u (x, t)|) p 2p+3 dx 2p+3 3 = 3 p N(t) -2p 3 K(t) 2p+3 3 . ( 6
)
Temporal minimization. To complete the resolution of optimization Problem (2-3), we perform a temporal minimization. We need to find the optimal time law t → N(t) for the instantaneous mesh size. We consider the case where the time-step τ is specified by the user as a function of time t → τ(t). A similar analysis is done in [START_REF] Alauzet | Multi-scale anisotropic mesh adaptation for time-dependent problems[END_REF] to deal with the case of an explicit time advancing solver subject to Courant time-step condition. After the previous spatial minimization, optimization Problem (2-3) becomes: find space-time continuous mesh M L p such that:

(E L p (M L p )) p = min M T 0 E L p (M L p (t)) dt = min M 3 p T 0 N(t) -2p 3 K(t) 2p+3 3 dt under constraint T 0 τ(t) -1 N(t) dt = N st .
After computations detailed in [START_REF] Alauzet | Multi-scale anisotropic mesh adaptation for time-dependent problems[END_REF], the expression of the optimal space-time continuous mesh M L p = (M L p (x, t)) (x,t)∈Q for a prescribed time-step is:

M L p (x, t) = N 2 3 st T 0 τ(t) -2p 2p+3 K(t)dt -2 3 τ(t) 2 2p+3 (det |H u (x, t)|) -1 2p+3 |H u (x, t)| . (7)
The following optimal error is obtained:

E L p (M L p ) = 3 N -2 3 st T 0 τ(t) -2p 2p+3 K(t)dt 2p+3 3p . (8)

Space-time L p error analysis with time sub-intervals

The previous analysis provides the optimal adapted meshes for each time level, it thus requires the mesh to be adapted at each flow solver time-step which is inconceivable in practical applications. In [START_REF] Alauzet | 3D transient fixed point mesh adaptation for time-dependent problems: Application to CFD simulations[END_REF][START_REF] Alauzet | Multi-scale anisotropic mesh adaptation for time-dependent problems[END_REF], the use of a coarse adapted discretization of the time axis is proposed. The basic idea consists in splitting the simulation time frame [0, T ] into n adap adaptation sub-intervals:

[0, T ] = [0 = t 1 , t 2 ] ∪ . . . ∪ [t i , t i+1 ] ∪ . . . ∪ [t n adap , t n adap +1 = T ] ,
and to keep the same adapted spatial mesh for each time sub-interval. On each sub-interval, the mesh is adapted to control the solution accuracy from t i to t i+1 . Consequently, the time-dependent simulation is performed with n adap different adapted meshes. This drastically reduces the number of meshes generated during the simulation, hence the number of solution transfers. Moreover, the flow solver performs many iterations (hundreds to thousands) on the same fixed spatial mesh. This provides a first answer to the adaptation of the whole space-time mesh, the spatial mesh being kept constant for each sub-interval when the global space-time mesh is visualized.

In the following, we extend the previous error analysis to the fixed-point unsteady mesh adaptation algorithm context where the simulation time interval [0,T] is split into n adap sub-intervals [t i , t i+1 ] for i = 1, .., n adap . Each spatial continuous mesh M i is kept constant during each sub-interval [t i , t i+1 ], thus it has no more dependency in time. As previously, a spatial minimization problem is solved first, followed by a temporal minimization resolution (see [START_REF] Alauzet | Multi-scale anisotropic mesh adaptation for time-dependent problems[END_REF][START_REF] Belme | Time accurate anisotropic goal-oriented mesh adaptation for unsteady flows[END_REF] for details of the demonstration).

Spatial minimization for a given sub-interval. Given the continuous mesh spatial complexity N i for the single adapted mesh used during sub-interval i, we seek for the optimal continuous mesh M i L p = (M i L p (x)) x∈Ω solution of the following problem:

E i L p (M i L p ) = min M i t i+1 t i Ω trace (M i ) -1 2 (x) |H u (x, t)| (M i ) -1 2 (x) p dx dt = min M i Ω trace (M i ) -1 2 (x) H i u (x) (M i ) -1 2 (x) p dx such that C(M i ) = N i , (9)
where matrix H i u on the sub-interval is defined by:

H i u (x) = t i+1 t i |H u (x, t)| dt . ( 10
)
This comes to moving the integral over time into the trace in Expression [START_REF]Time-accurate anisotropic mesh adaptation for three-dimensional moving mesh problems[END_REF]. This can be done because M i does not depend anymore on t, i.e., the mesh is static. Note that this cannot extend directly when meshes become dynamic. The spatial minimization gives the solution of Problem [START_REF]Time-accurate anisotropic mesh adaptation for three-dimensional moving mesh problems[END_REF] where the expressions of the optimal continuous mesh and of the optimal error are identical to Relations ( 5) and ( 6) where H u is substituted with H i u .

Temporal minimization. Again, we consider the case where the time-step τ is specified by a function of time. After the spatial minimization, the temporal minimization problem becomes:

(E L p (M L p )) p = min M n adap i=1 E i L p (M i L p ) such that n adap i=1 N i        t i+1 t i τ(t) -1 dt        = N st .
After the temporal minimization, the expression of the optimal continuous mesh M L p = {M i L p } i=1,..,n adap = (M i L p (x)) x∈Ω i=1,..,n adap and error are:

M i L p (x) = N 2 3 st n adap j=1 K j t j+1 t j τ(t) -1 dt 2p 2p+3 -2 3 t i+1 t i τ(t) -1 dt -2 2p+3 (det H i u (x)) -1 2p+3 H i u (x) (11)
E L p (M L p ) = 3 N -2 3 st n adap i=1 K i t i+1 t i τ(t) -1 dt 2p 2p+3 2p+3 3p , (12)
where

K i = Ω det H i u (x) p 2p+3 dx.

Global fixed-point mesh adaptation algorithm

In [START_REF] Alauzet | 3D transient fixed point mesh adaptation for time-dependent problems: Application to CFD simulations[END_REF], a local transient fixed-point mesh adaptation algorithm for time-dependent simulations has been proposed. However, this algorithm cannot be used here as the computation of the optimal continuous mesh given by Relation (11) involves a global normalization term which requires the knowledge of quantities over the whole simulation time frame. Thus, the complete simulation must be performed before evaluating any space-time continuous mesh. To this end, we consider a global fixed-point mesh adaptation algorithm covering the whole time frame [0, T ]. This iterative algorithm is used to converge the non-linear mesh adaptation problem, and also to predict the solution evolution.

The global fixed-point unsteady mesh adaptation algorithm is schematized in Algorithm 1 where H, S and M denote respectively meshes, solutions and metric fields (i.e., continuous meshes). H is the Hessian-metric given by Relation [START_REF] Barral | Large displacement body-fitted FSI simulations using a mesh-connectivity-change moving mesh strategy[END_REF]. Several steps of the previous algorithm are not straightforward, and need to be further detailed.

Algorithm 1 Mesh Adaptation Loop for Unsteady Flows

Initial mesh and solution(H 0 , S 0 0 ) and set targeted space-time complexity N st // Fixed-point loop to converge the global space-time mesh adaptation problem For j = 1, n pt f x 1. // Adaptive loop to advance the solution in time on time frame [0, T ] For i = 1, n adap (a) S i, j 0 = Interpolate conservatively next sub-interval initial sol. from (H i-1, j , S i-1, j , H i, j ); (b) S i, j = Compute solution on sub-interval from pair (S i, j 0 , H i, j ); (c) |H| i, j = Compute sub-interval Hessian-metric from sol. sample (H i, j , {S i, j (k)} k=1,nk ); EndFor 2. C j = Compute space-time complexity from all Hessian-metrics ({|H| i, j } i=1,n adap );

3. {M i, j L p } i=1,n adap = Compute all sub-interval unsteady metric fields (C j , {|H max | i, j } i=1,n adap ); 4. {H i, j+1 } i=1,n adap = Generate all sub-interval adapted meshes ({H i, j , M i, j
L p } i=1,n adap ); EndFor

Computation of the optimal continuous mesh

Computation of the Hessian-metric. The optimal L p metric involves an averaged Hessian-metric H i u on sub-interval i given by Relation [START_REF] Barral | Large displacement body-fitted FSI simulations using a mesh-connectivity-change moving mesh strategy[END_REF], but it still remains to know how to compute it practically, i.e., how it is discretized. The strategy adopted in [START_REF] Alauzet | 3D transient fixed point mesh adaptation for time-dependent problems: Application to CFD simulations[END_REF] is to sample the solution on the time sub-interval. More precisely, n k solutions equally distributed on the sub-interval time frame are saved, including the initial solution at t i and the final solution at t i+1 . Positive Hessian |H u (x, t k )| is evaluated for each sample. The Hessian matrices are computed using a double least-square procedure: first the gradients are computed using a linear least-square reconstruction procedure [START_REF] Mavriplis | Revisiting the least-square procedure for gradient reconstruction on unstructured meshes[END_REF], and a similar procedure is applied to recover the Hessians. The error analysis leads to write H i u as the integral over time of the Hessian matrices, and thus the following discretization is used:

H i u (x) = 1 2 ∆t i n k -1 |H u (x, t i )| + ∆t i n k -1 n k -1 k=2 |H u (x, t k )| + 1 2 ∆t i n k -1 |H u (x, t i+1 )| = ∆t i |H i avg (x)| ,
where ∆t i = t i+1t i is the sub-interval time length and

t k = t i + k-1 n k -1 ∆t i .
Choice of the optimal continuous mesh. The optimal adapted mesh for each sub-interval is obtained according to analysis of Section 2.3. For the numerical results presented below, we select the optimal mesh given by Relation [START_REF] Barral | Two mesh deformation methods coupled with a changing-connectivity moving mesh method for CFD applications[END_REF] and the following particular choices have been made:

• the space-time error is controlled in L 2 norm (p = 2)
• all sub-intervals have the same time length ∆t = T/n adap

• function τ : t → τ(t) is constant and equal to ∆t, thus

t i+1 t i τ(t) -1 dt = 1.
In that case, we can just consider the time-step τ(t) constant and equal to ∆t (thanks to the global normalization term) so the integral t i+1 t i τ(t) -1 dt reduces to 1 (see [START_REF] Alauzet | Multi-scale anisotropic mesh adaptation for time-dependent problems[END_REF] for details). With these choices, the optimal continuous mesh M

L 2 = {M i L 2 } i=1,..,n adap = {(M i L 2 (x)) x∈Ω } i=1,.
.,n adap simplifies to:

M i L 2 (x) = N 2 3 st         n adap j=1 Ω (det |H j u (x)|) 2 7 dx         -2 3 det |H i u (x)| -1 7 |H i u (x)| . (13)
In that case, as we assume that theoretically one time-step is done by sub-interval, N avg = N st /n adap represents the average spatial complexity by sub-interval. We do not prescribe the temporal complexity, i.e., we do not control the number of time-steps done at each sub-interval.

In practice, the user prescribes the number of sub-intervals n adap and the sub-interval average spatial complexity N avg leading to a total space-time complexity of

N st = n adap × N avg . ( 14
)
This prescription is the theoretical complexity. In practice, the total number of space-time vertices N st of the simulation discrete meshes is directly proportional to the prescribed total space-time complexity N st = c N st , where coefficient c depends on the geometry of the problem, the mesh gradation, and the local remesher. The total number of space-time vertices of the simulation run on n adap adapted meshes is computed as follow:

N st = n adap i=1 n i iter × N i ,
where the time discretization corresponds to the number of time-steps of the flow solver at each sub-interval, e.g. n i iter is the number of time-steps for the i th sub-interval, and N i is the number of spacial vertices for the i th sub-interval.

Comments on parameters N st , N avg and n adap . A very detailed analysis and convergence study of these parameters has been done in [START_REF] Alauzet | Multi-scale anisotropic mesh adaptation for time-dependent problems[END_REF]. The total space-time complexity N st can be increased by fixing n adap and increasing N avg or by fixing N avg and increasing n adap . The finer adapted meshes obtained with these two approaches have the same theoretical complexity but increasing n adap leads to a mesh having less space-time vertices than increasing N avg . In other words, a better adapted space-time mesh is obtained by increasing n adap . Therefore, when convergence studies are performed, the number of sub-intervals should be increased to increase the space-time complexity in order to obtain a better rate of convergence.

Matrix-free P 1 -exact conservative solution transfer

Between each sub-interval, the solution needs to be transferred from the previous mesh to the next one to pursue the computation. This stage becomes critical in the context of unsteady problems and even more if a large number of transfers is performed, as the error introduced by this stage can spoil the overall accuracy of the solution. In the context of the resolution by a second order numerical scheme of a PDE system of conservation laws, such as the compressible Euler system, it is crucial for the interpolation method to satisfy the following properties in order to obtain a consistent mesh adaptation scheme: (i) mass conservation, (ii) P 1 -exactness preserving the second order of the adaptive strategy and (iii) verify the maximum principle. This matrix-free P 1 -exact conservative solution transfer method, and its properties, is thoroughly described in [START_REF] Alauzet | A parallel matrix-free conservative solution interpolation on unstructured tetrahedral meshes[END_REF].

The remeshing step

The remeshing step is also crucial in the adaptation process, as a poorly adapted mesh or a mesh with bad quality elements will impact the accuracy of the solution and spoil the efforts on all other parts of the process. In this paper, the remesher Amg [START_REF] Loseille | Cavity-based operators for mesh adaptation[END_REF] was used. Amg belongs to the class of 3D anisotropic local remeshers that aim at generating a unit mesh with respect to a prescribed metric field. Its main particularity is to adapt the volume and the surface mesh in a coupled way so that a valid 3D mesh is always guaranteed in output. It uses a unique cavity-based operator (that generalizes standard operators: point insertion, edge removal, edges swapping, point smoothing). This operator is governed by dedicated metric-based quality functions. This allows us to reach a good balance between high level of anisotropy, necessary to capture the physical features of the solution, and mesh quality, necessary to ensure the stability and accuracy of the numerical scheme. Amg has several enhancements designed for CFD computations, including explicit control and optimization of tetrahedra height to ensure a maximal time-step for unsteady simulations, and surface mesh re-projection based either on CAD or on background discrete meshes.

The ALE flow solver

We consider the compressible Euler equations for a Newtonian fluid in their Arbitrary-Lagrangian-Eulerian (ALE) formulation. The ALE formulation allows the equations to take arbitrary motion of the mesh into account. Assuming that the fluid is a perfect gas and that there is no thermal diffusion, the ALE formulation of the Euler equations is written for any arbitrary closed volume C(t) of boundary ∂C(t) moved with mesh velocity w:

d dt C(t) Wdx + ∂C(t) (F (W) • n -W(w • n)) ds = C(t) F ext dx (15)
where

             W = (ρ, ρu, ρe)
T is the conservative variables vector F (W) = ρu, ρu x u + pe x , ρu y u + pe y , ρu z u + pe z , ρuh is the flux tensor

F ext = (0, ρ f ext , ρu • f ext )
T is the contribution of the external forces and we have noted ρ the density of the fluid, p the pressure, u = (u x , u y , u z ) its Eulerian velocity, q = u , ε the internal energy per unit mass, e = 1/2 q 2 + ε the total energy per unit mass, h = e + p/ρ the enthalpy per unit mass of the flow, f ext the resultant of the volumic external forces applied on the particle and n the outward normal to interface ∂C(t) of C(t). The fixed mesh case is obtained by simply taking w = 0.

The following examples were run using our in-house flow solver, Wolf. The spatial discretization relies on an edge-based vertex-centered Finite Volume numerical scheme using an HLLC Riemann approximate solver to compute the numerical fluxes. A second-order scheme is derived according to a MUSCL type method with a low dissipation gradient reconstruction combined with a generalization of the Superbee limiter with three entries. The main difference when translating these schemes from the standard formulation to the ALE formulation is the addition of the mesh velocities w in the wave speeds of the Riemann problem and the re-evaluation of the finite volume cells at each time-step [START_REF] Barral | Large displacement body-fitted FSI simulations using a mesh-connectivity-change moving mesh strategy[END_REF]. In this work, the temporal discretization is explicit and based on strong stability preserving Runge-Kutta (SSPRK) schemes [START_REF] Spiteri | A new class of optimal high-order strong-stability-preserving time discretization methods[END_REF]. The 5-stage order-2 SSPRK is considered with a CFL number equal to 3. The discrete geometric conservation law (DGCL) is strictly enforced at the discrete level by determining when and how geometrical parameters that appear in the fluxes should be computed [START_REF] Olivier | A new changing-topology ALE scheme for moving mesh unsteady simulations[END_REF].

Parallelization of the mesh adaptation loop

All the steps of the adaptation loop have been parallelized. Two different approaches are used for the two big parts of the loop. The solution computation and the solution transfer procedure are parallelized using a p-thread paradigm at the element loop level [START_REF] Alauzet | On the use of space filling curves for parallel anisotropic mesh adaptation[END_REF]. As regards the computation of the metrics, the metric gradation and the generation of the adapted meshes, a pipeline approach was used. These operations have to be done at the end of the loop, as many times as there are sub-intervals, but the operations for one sub-interval are totally independent from the operations for another. Consequently, they can all be run in parallel: if N processors are available, N metrics can be computed and their associated meshes generated simultaneously in serial on one processor.

Example of a 3D blast on the London Tower Bridge

Let us illustrate the efficiency of the time-accurate multiscale anisotropic mesh adaptation method on an example involving a very complex geometry. We consider a blast on the London Tower Bridge, whose geometry was the object of the meshing contest of the 23rd International Meshing Roundtable. A CAD description of the bridge was provided, from which an initial surface (Figure 1, left) and volume mesh were generated. This very detailed geometry is a challenge for mesh adaptation, because it results in complex pattern of the solution that the error estimate has to be able to capture, and because the adaptive remesher has to be able to preserve them and preserve the anisotropy in their vicinity.

The gas is initially at rest in a box of dimension [-130, 130]×[-72, 77]×[0, 100]. The gas is in a constant state W 0 = (1, 0, 2.5) everywhere, except for a sphere of radius 1 centered in (3.5, 20, 10) containing a high density and high energy gas defined by W blast = (10, 0, 250). The gas is then left to evolve freely until dimensionless time T = 25, and wall conditions are imposed on the boundaries. The spherical shock expands in the volume, hits the walls of the towers, and insinuates itself in any hole of the geometry, see Figure 3 (right). We observe numerous shock waves reflected by the walls that then interact with each other, and even recirculations around the pillars of the bridge. Many reflected shock waves impact the low-density bubble region located at the center of the explosion. A mushroom-shape instability develops from all these shock-bubble interactions as can be seen in Figure 2. We observe that all these phenomena -the shocks and as many of their reflections as possible as well as the smaller scale instabilities of the mushroom -are automatically captured and highly resolved thanks to the time-accurate multiscale anisotropic mesh adaptation method.

To obtain that highly resolved result, the simulation interval was split into n adap = 80 sub-intervals and a theoretical average spatial complexity of N avg = 800, 000 per sub-interval has been set. This represents a total space-time complexity N st equal to 64 million. The density field of the solution was chosen as sensor for the adaptation. Five fixed-point iterations were used to converge the non-linear mesh adaptation problem. Practically, at the last fixed-point iteration, the discrete meshes have an average number of spatial vertices of N avg = 2, 970, 984 vertices and a total of n iter = 16, 353 time-steps have been performed. The total number of space-time vertices used to compute that solution is N st = 47.7 billion. To give an idea of the accuracy of the mesh, at t = 3.125, the mesh accuracy is around 1.4 e -3 , and 1.1 e -3 at t = 25. The CPU time of the last fixed-point iteration (error estimate, adaptive mesh generation, solution interpolation, and flow solver) is 21 hours on 20 cores using a 2 Xeon E5-2670 v2 chip processor (10 cores at 2.5 GHz for each chip), both chips being connected by 2 QPI links with a speed of 16 GB/s. Snapshots of the adapted meshes and solutions are shown in Figure 3. The mesh adaptation for the whole sub-interval is clearly illustrated. Indeed, the mesh refinement along band-shaped regions, which is typical of the fixed-point algorithm, is clearly visible. These band-shaped areas correspond to the evolution zone of the physical phenomena during an adaptation sub-interval.

Extension of the space-time L p error analysis to moving mesh simulations

In the previous sections, we have only considered fixed meshes in the adaptation process. Due to its efficiency in those cases, it seems natural to try to extend the global fixed-point mesh adaptation algorithm (Algorithm 1) to moving mesh simulations. In the context of body-fitted moving mesh simulations, the mesh is deformed to follow the displacement of the moving geometry inside the domain. The mesh deformation consists in prescribing a displacement to each mesh vertex by solving a PDE or by using an interpolation method where the boundary conditions are the geometry displacement. Thus, the general ideas driving the global fixed-point mesh adaptation algorithm are still valid: iterating the process to converge the couple mesh/solution and splitting the simulation time interval into sub-intervals. By design, we are only remeshing between two adaptation sub-intervals, which leads to a unique dynamic mesh on each sub-interval that is deformed because of the moving geometry. But, the previous adaptation algorithm cannot directly be extended to moving mesh simulations because it does not take into account the movement of the mesh. The error analysis of Section 2.3 from which the optimal continuous mesh is deduced, assumes that the mesh is constant in time inside a sub-interval, therefore it does not take into account the local deformation of the mesh which necessarily impacts local errors. Indeed, the mesh deformation is driven by the boundary displacement, so the mesh distortion does not necessarily follow the physical phenomena. This may in particular interfere with the adaptation of the mesh.

The previous error analysis has to be modified to take into account the movement of the mesh. As regards the discrete representation of the metric field for moving meshes, two options are available. The first approach is a Eulerian approach: a metric is associated with a fixed position in space, thus the metric field is a "background" field evolving in space and time but independently from the moving mesh. To know the value of the metric at a vertex, an interpolation has to be performed. But, the representation of this metric on a background mesh is problematic due to the displacement of the boundaries and to the definition of an adequate background mesh to support the metric field. That is why the approach considered in this work is a Lagrangian approach: a metric is attached to a moving vertex and moves with it, i.e., we have M(x(t), t) which we will write M(x(t)) in what follows.

In the context of dynamic meshes, the optimization problem is still given by Relations ( 2) and (3) where the space-time error model is given by Relation (4). As in Section 2.3, to perform the space-time L p error analysis with time sub-intervals, we first do a spatial minimization for a given sub-interval. We consider the i th sub-interval [t i , t i+1 ]. Given the continuous mesh spatial complexity N i , we seek for the optimal Arbitrary-Lagrangian-Eulerian (ALE) continuous mesh

M i,ALE L p = M i,ALE L p (x(t)) x∈Ω(t)
for the whole sub-interval which is solution of the following problem:

E i,ALE L p (M i,ALE L p ) = min M i t i+1 t i Ω(t) trace (M i ) -1 2 (x(t)) |H u (x(t), t)| (M i ) -1 2 (x(t)) p dx(t) dt such that C(M i ) = N i .
Note that the continuous mesh spatial complexity is constant on this sub-interval. We name it optimal ALE continuous mesh to stress that the continuous mesh is dynamic (i.e., it evolves in time), it represents a moving adapted mesh and it is coupled with an ALE flow solver. Now M i has a dependence in time, thus we cannot move the integral over time into the trace to make the mean Hessian metric appear like in the fixed mesh analysis. The spatial minimization for a given sub-interval cannot be done directly.

To solve this issue, we first exhibit the optimal instantaneous ALE continuous mesh minimizing the interpolation error in L p norm at a given time t, also shortly called ALE metric field, that takes into account the mesh deformation (Section 4.1) and then we validate it on several analytical examples (Section 4.2). This ALE metric field makes it possible to map a continuous moving mesh of Ω(t), t ∈ [t i , t i+1 ], onto a continuous mesh of Ω(t i ). Then, the error model can be written only on Ω(t i ) and the spacial minimization can be performed (Section 4.3). In other words, to preserve the adaptation of the mesh despite the mesh deformation, we are looking for a mesh at the beginning of each sub-interval that will be adapted to the solution at each time-step of the sub-interval once moved with the prescribed mesh displacement.

Optimal instantaneous ALE continuous mesh minimizing the interpolation error in L p norm: the ALE metric

Before dealing with the case of the mesh adaptation for a sub-interval, let us solve the previous mentioned problem for one time step. The computational domain is time-dependent Ω(t) ⊂ R 3 × [0, T ]. The problem for one time-step is the following. t n and t n+1 are two times, Ω n and Ω n+1 are the spatial domains at t n and t n+1 respectively, and generally we have Ω n Ω n+1 . We denote by d the given mesh displacement field. Then, we want to find the optimal continuous mesh M n,ALE L p = (M n,ALE L p (x)) x∈Ω n defined on Ω n from which we will generate a mesh at time t n that, once moved with displacement d, will be adapted to a sensor u n+1 at time t n+1 on domain Ω n+1 . We consider that u n+1 is a scalar sensor function, the extension to vector sensor functions being straightforward. The resolution of this problem leads to the optimal instantaneous ALE continuous mesh defined on Ω n minimizing the interpolation error of sensor u n+1 in L p norm at time t n+1 on Ω n+1 after being deformed by displacement d. This ALE metric field involves the gradient of the mesh transformation between t n and t n+1 to take into account the mesh deformation. A first draft of this analysis was given in [START_REF] Alauzet | Extension of metric-based anisotropic mesh adaptation to time-dependent problems involving moving geometries[END_REF].

For the coming mathematical analysis involving spatial derivatives, it is very important to state on which domain/mesh these derivatives are computed. To this end, the following notations will be used:

• ∇ n denotes the gradient operator3 performed on domain Ω n , i.e., computed on mesh H n • H n+1 denotes the Hessian operator performed on domain Ω n+1 , i.e., computed on mesh H n+1 . As operator H n+1 is always applied to sensor u n+1 at time t n+1 , simplified notation H n+1 u will stand for

H n+1 [u n+1 ] • M n+1 L p [u n+1
] denotes the point-wise optimal L p metric for sensor u n+1 computed on domain Ω n+1 , i.e., on mesh H n+1 . Again, simplified notation M n+1 L p will stand for

M n+1 L p [u n+1 ] • M n+1 L p = (M n+1 L p (x)) x∈Ω n+1
is the associated continuous mesh with complexity C(M n+1 L p ) = N n+1 . Now, let us introduce φ the mapping between domains4 Ω n and Ω n+1 :

φ : Ω n -→ Ω n+1 x n -→ x n+1 = φ (x n ) , (16)
and d the corresponding mesh displacement field, such that:

x n+1 = φ(x n ) = x n + d(x n ) . ( 17
)
Since φ is a diffeomorphism, we have, for any infinitesimal vector δx n ∈ Ω n :

δx n+1 = ∇ n φ(x n ) T δx n with ∇ n φ(x n ) = I + ∇ n d(x n ) , ∀x n ∈ Ω n . (18)
Finally, a . operator is defined which transports a quantity from Ω n+1 to Ω n . We note H n+1 u the Hessian of u n+1 computed on Ω n+1 and transported on domain Ω n . This mathematically writes:

H n+1 u : Ω n -→ R x n -→ H n+1 u (φ(x n )) = H n+1 u (x n ) .
According to error analysis of Section 2.2, the optimal instantaneous continuous mesh M n+1 L p = (M n+1 L p (x)) x∈Ω n+1 at time t n+1 is given by Relation [START_REF] Alauzet | High order sonic boom modeling by adaptive methods[END_REF]. Thus, an optimal mesh of Ω n+1 adapted to u n+1 can be built by generating a unit mesh H n+1 with respect to M n+1 L p . This means that the length of any arbitrary edge e n+1 of H n+1 is one in metric M n+1 L p (see Section 2.1):

1 = (e n+1 (x n+1 )) T M n+1 L p (x n+1 ) e n+1 (x n+1 ) , ∀e n+1 (x n+1 ) ∈ H n+1 . ( 19
)
Let e n be the edge of H n having e n+1 as image by φ in H n+1 . As we are only interested in controlling the prevailing term of the interpolation error, we can use the following relation:

e n+1 (x n+1 ) = e n+1 ( φ(x n ) ) = ∇ n φ(x n ) T e n (x n ) , (20)
which is true at first order in the demonstration. The idea is to unravel how Condition [START_REF] Coupez | Adaptive time-step with anisotropic meshing for incompressible flows[END_REF] writes when transposed onto mesh H n . Using Relations ( 5) and [START_REF] De Sampaio | Petrov-Galerkin solutions of the incompressible Navier-Stokes equations in primitive variables with adaptive remeshing[END_REF], unit length Equation ( 19) can be re-written:

1 = e n+1 ( φ(x n ) ) T M n+1 L p ( φ(x n ) ) e n+1 ( φ(x n ) ) = ∇ n φ(x n ) T e n (x n ) T         N n+1 K n+1 2 3 det |H n+1 u (φ(x n ))| -1 2p+3 H n+1 u (φ(x n ))         ∇ n φ(x n ) T e n (x n ) = e n (x n ) T         N n+1 K n+1 2 3 det | H n+1 u (x n )| -1 2p+3 ∇ n φ(x n ) H n+1 u (x n ) ∇ n φ(x n ) T         e n (x n ) = e n (x n ) T M n,ALE L p (x n ) e n (x n ) .
Consequently, if we create a unit mesh of Ω n with respect to ALE continuous mesh M n,ALE

L p = M n,ALE L p (x)
x∈Ω n , rewriting the above calculus upside down, we get the following implication:

e n (x n ) T M n,ALE L p (x n ) e n (x n ) = 1 =⇒ e n+1 (x n+1 ) T M n+1 L p (x n+1 ) e n+1 (x n+1 ) = 1 ,
for all edges e n of mesh H n having e n+1 as image by φ in mesh H n+1 . Therefore, if a unit mesh is generated with respect to M n,ALE L p at t n , then the deformed mesh moved with displacement d until time t n+1 is unit for the optimal metric associated with sensor u n+1 , meaning that it is optimal to control the interpolation error in L p norm of the sensor u at t n+1 .

The optimal instantaneous ALE continuous mesh reads:

M n,ALE L p (x n ) = N n+1 K n+1 2 3 det | H n+1 u (x n )| -1 2p+3 ∇ n φ(x n ) H n+1 u (x n ) ∇ n φ(x n ) T , with K n+1 = Ω n+1 det |H n+1 u (x n+1 )| p 2p+3 dx n+1 .
The number of vertices of the mesh remains the same when it is moved with displacement d. Indeed, for our moving mesh simulations, connectivity-changes can occur but no vertex addition or deletion is allowed (see Section 5.1). Consequently, the continuous mesh complexity remains constant in time:

C M n+1 L p = C M n,ALE L p ⇐⇒ Ω n+1 det M n+1 L p (x n+1 )dx n+1 = Ω n det M n,ALE L p (x n )dx n ,
thus, according to the definition of K and using det(αA) = α 3 det A (we are in R 3 ), we have:

K n+1 = K n,ALE = Ω n det det | H n+1 u (x n )| -1 2p+3 ∇ n φ(x n ) H n+1 u (x n ) ∇ n φ(x n ) T 1 2 dx n = Ω n det ∇ n φ(x n ) det | H n+1 u (x n )| p 2p+3 dx n Moreover, M n,ALE L p
can be re-written to recover the same form as Relation (5) which will be useful for the error analysis involving dynamic meshes. To this end, we set:

|H n,ALE u (x n )| = det ∇ n φ(x n ) 1 p ∇ n φ(x n ) | H n+1 u (x n )| ∇ n φ(x n ) T , (21)
from which we deduce:

det |H n,ALE u (x n )| = det ∇ n φ(x n ) 3 p +2 det | H n+1 u (x n )| and K n,ALE = Ω n det |H n,ALE u (x n )| p 2p+3 dx n .
After some algebra, we finally get:

M n,ALE L p (x n ) = N n+1 K n,ALE 2 3 det |H n,ALE u (x n )| -1 2p+3 |H n,ALE u (x n )| .
According to Section 2.2, this expression of the ALE metric field represents the optimal instantaneous continuous mesh minimizing the interpolation error in L p norm for the sensor having as second derivatives H n,ALE u on Ω n . Therefore, we can derive the expression of the interpolation error associated with continuous mesh M n,ALE L p . The instantaneous interpolation error in L p norm at the power p is given by Relation ( 6):

E L p (M n,ALE L p ) = Ω n trace (M n,ALE L p ) -1 2 (x n ) |H n,ALE u (x n )| (M n,ALE L p ) -1 2 (x n ) p dx n = 3 p (N n+1 ) 2p 3 (K n,ALE ) 2p+3 3 = 3 p (N n+1 ) 2p 3 (K n+1 ) 2p+3 3 = Ω n+1 trace (M n+1 L p ) -1 2 (x n+1 ) |H n+1 u (x n+1 )| (M n+1 L p ) -1 2 (x n+1 ) p dx n+1 = E L p (M n+1 L p ) .

Now, we can state the main result:

Proposition 1 (Optimal instantaneous ALE continuous mesh). Let Ω n and Ω n+1 be two spatial domains at t n and t n+1 , and φ the mapping between these domains. The optimal instantaneous ALE continuous mesh M n,ALE

L p = M n,ALE L p (x)
x∈Ω n defined on Ω n reads:

M n,ALE L p (x) = N n+1 2 3 Ω n det |H n,ALE u (x)| p 2p+3 dx -2 3 det |H n,ALE u (x)| -1 2p+3 |H n,ALE u (x)| , (22)
with the ALE Hessian-metric given by

|H n,ALE u (x)| = det ∇ n φ(x) 1 p ∇ n φ(x) | H n+1 u (x)| ∇ n φ(x) T .
The following properties hold: i) Let us assume continuous mesh M n,ALE L p is used to generate a unit mesh H n ALE of Ω n and let us denote by H n+1 ALE the mesh of Ω n+1 which is the image of mesh H n ALE by mapping φ. Then, mesh H n+1 ALE is optimal to control the interpolation error in L p norm of sensor u n+1 on Ω n+1 . ii) Continuous mesh M n,ALE L p has the same complexity N n+1 as continuous mesh M n+1 L p . iii) Continuous mesh M n,ALE L p achieves on Ω n the same level of interpolation error as continuous mesh M n+1 L p on Ω n+1 :

E L p (M n,ALE L p ) = E L p (M n+1 L p
) where E L p is given by Relation (6).

iv) There is no reason for adapted mesh H n ALE , which is a unit mesh for M n,ALE L p , to be optimal for the control of the interpolation error of sensor u n at t n .

Validation on two-dimensional analytic examples

We validate the result of Proposition 1 on analytic examples in two dimensions, three-dimensional examples being available in [START_REF]Time-accurate anisotropic mesh adaptation for three-dimensional moving mesh problems[END_REF][START_REF] Olivier | Anisotropic metric-based mesh adaptation for unsteady CFD simulations involving moving geometries[END_REF]. A uniform mesh H n 0 of domain Ω = [-1, 1] 2 , a displacement field d between two times t n and t n+1 , and a sensor u n+1 at time t n+1 are given. We choose to control the interpolation error in L 1 norm and theoretical spatial complexities from 10, 000 to 400, 000 are prescribed to perform a convergence study. We generate a unit mesh H n,ALE with respect to M n,ALE L 1 (given by Relation ( 22)) at time t n that, once moved into H n+1,ALE = φ(H n,ALE ) is expected to be adapted to the sensor u n+1 according to the above developments. At the same time, a reference adapted mesh H n+1 is generated that is directly adapted to sensor u n+1 using optimal continuous mesh M n+1 L 1 (given by Relation ( 5)). To validate the optimality of H n+1,ALE , we are going to compare the interpolation error in L 1 norm of sensor u n+1 on both meshes and analyse the quality of the meshes elements. Algorithm 2 describes the procedure to generate both meshes. Note that, for the ALE cases, the meshes are directly moved from their initial position at t n to their final position at t n+1 in one step. No mesh optimizations are performed after the displacement.

Algorithm 2 Procedure to generate adapted meshes with the ALE L 1 metric field and the regular L 1 metric field.

• Generation of an adapted mesh with the optimal ALE L 1 metric field

1. ∇ n φ = Compute transformation gradient H n 0 , d 2. H n+1 0 = Move mesh H n 0 , d 3. u n+1 = Compute target sensor H n+1 0 4. M n,ALE L 1 = Compute ALE metric field H n+1 0 , u n+1 , H n 0 , ∇ n φ 5. H n,ALE = Adapt mesh H n 0 , M n,ALE L 1 6. H n+1,ALE = Move mesh H n,ALE , d
• Generation of an adapted mesh with the optimal L 1 metric field

1. u n+1 = Compute target sensor H n 0 2. M n+1 L 1 = Compute metric field H n 0 , u n+1 3. H n+1 = Adapt mesh H n 0 , M n+1 L 1
The following two-dimensional analytic functions are used as sensor functions:

u n+1 1 (x, y) =            0.01 sin (50xy) if |xy| ≥ 2π 50 sin (50xy) if |xy| < 2π 50
and u n+1 2 (x, y) = 0.1 sin(50 x) + arctan 0.1 sin(5 y) -2 x .

They are shown in Figures 5 and6. First sensor exhibits features of different scales, small oscillations of amplitude 0.01 and large oscillations of amplitude 1, which makes it a suitable example for our multiscale mesh adaptation process. The second sensor exhibits a strong discontinuity-like sinusoidal feature -through the arctangent -crossing small amplitude waves.

Two displacements d i : Ω n -→ Ω n+1 are considered:

d 1 (x, y) =                                       -0.3 (x + 1) y 2 -1 exp -5x 2 , if x ≥ 0 0.3 (x -1) y 2 -1 exp -5x 2 , if x < 0          -0.3 x 2 -1 (y + 1) exp -5y 2 , if y ≥ 0 0.3 x 2 -1 (y -1) exp -5y 2 , if y < 0                              and d 2 (x, y) =         0.5 (x 2 -1) (y 2 -1) 0         .
They are illustrated in Figure 4 where the displacements have been applied to a uniform-size mesh composed of 53, 262 vertices (right) to the initial uniform-size mesh. These displacements are used to validate the ALE metric formulation.

Figure 5: Mesh adapted with the ALE metric for a spatial complexity of 50, 000 at time t n (left) and t n+1 (middle)i.e., after moving the mesh -for analytical function u n+1 to emphasize the mesh distortion induced by these mesh deformations. It is important to note that these functions lead to very large displacement: d 1 max ≈ 0.1944 and d 2 max = 0.5 which are a lot larger than the mesh size in all cases. Hence, even if the theory considers infinitesimal displacement, these test cases validate the result in a more realistic context.

Figures 5 and6 show the resulting meshes for both sensors and both displacements for a prescribed theoretical spatial complexity of 50, 000 leading to discrete adapted meshes composed of almost 60, 000 vertices. Initial adapted meshes H n,ALE at time t n are displayed on the left picture, and the resulting adapted mesh H n+1,ALE at time t n+1 after applying the mesh displacement are shown on the middle picture. Note that the optimal meshes H n+1 obtained directly from M n+1 L 1 are very similar to the meshes shown on the middle picture.

A quantitative analysis of these examples is required to make sure our ALE procedure results in effectively adapted meshes. To this end, we compare the mesh generated with a classic adaptation procedure and the mesh generated with the ALE procedure through two quantities:

• The interpolation error u n+1 -Π h u n+1 Ω committed when the sensor function is projected on the mesh remains the best way to evaluate if the mesh is really adapted to the considered sensor.

• To analyse any possible distortion of the mesh due to the mesh deformation, we can check that the shape of the elements is in accordance with the reference adapted metric. To this end, we compute the mesh elements quality with respect to metric field M n+1 L 1 used to generate the reference adapted mesh. The 2D element quality in the metric is given by:

Q M n+1 L 1 (K) = √ 3 12 3 i=1 2 M n+1 L 1 (e i ) |K| M n+1 L 1 ∈ [1, +∞]
where the e i are the three edges of the triangle, and length and area formula in the metric are given in Section 2.1. To give an idea, Q M (K) = 1 corresponds to a perfectly regular element in the metric, Q M (K) < 2 correspond to excellent quality elements, while a high value (> 100) of Q M (K) indicates a nearly degenerate element.

Figure 6: Mesh adapted with the ALE metric for a spatial complexity of 50, 000 at time t n (left) and t n+1 (middle)i.e., after moving the mesh -for analytical function u n+1

The interpolation error convergence curves are represented in Figure 7 where the interpolation error in L 1 norm is plotted with respect to the number of vertices. We observe that, for both sensors and both displacements, the interpolation error obtained with the ALE continuous mesh is almost similar or only slightly higher than the interpolation error obtained with the direct adaptation. This slight difference is due notably to approximations made when establishing the formulation of the ALE continuous mesh whereas large displacement are performed between t n and t n+1 . This confirms the optimality of the ALE continuous mesh, i.e., the property i) of Proposition 1. Note that the interpolation error is one order of magnitude smaller on adapted meshes than on uniform meshes, or that ten times more vertices are required on uniform meshes to achieve the same level of interpolation error. This clearly establishes the interest of the adaptation.

In regard to mesh elements quality, results are given in Tables 1 and2. The reference adapted meshes have an almost perfect quality with respect to the prescribed metric M n+1 L 1 , which was expected and is achieved thanks to the efficient local remesher Amg [START_REF] Loseille | Cavity-based operators for mesh adaptation[END_REF]. The quality of the meshes adapted with the ALE procedure at time t n+1 with respect to metric M n+1 L 1 is very good, which shows that resulting adapted meshes from both standard and ALE continuous mesh are very similar in terms of element-shape quality. This points out that the mesh deformation is perfectly taken into account in the ALE metric field formulation. Again, this validates property i) of Proposition 1.

Space-time error analysis for dynamic meshes

Section 4.1 provides the optimal instantaneous ALE continuous mesh which takes into account the mesh deformation. Now, we can extend the space-time error analysis with time sub-intervals done for fixed meshes -Section 2.3 -to the case of dynamic meshes. The simulation time interval is split into n adap sub-intervals. On each sub-interval, the mesh size (number of vertices) remains constant, but the mesh is deformed to follow the geometry displacement. At each time-step of the sub-interval, we want the moved mesh to be adapted to the current sensor. The key idea to perform the error analysis is to seek for the optimal dynamic continuous mesh at the beginning of the sub-interval, this continuous mesh being optimal for the whole sub-interval when deformed, instead of seeking for the expression of the optimal continuous mesh at each instant, i.e., as a function of the time. To perform the spatial minimization for a sub-interval, we consider the i th sub-interval [t i , t i+1 ]. Given the continuous mesh spatial complexity N i , we seek for the optimal ALE continuous mesh M i,ALE

L p = M i,ALE L p (x(t)) x∈Ω(t)
for the whole sub-interval which is solution of the following problem:

E i,ALE L p (M i,ALE L p ) = min M i t i+1 t i Ω(t) trace (M i ) -1 2 (x(t)) |H u (x(t), t)| (M i ) -1 2 (x(t)) p dx(t) dt , such that C(M i ) = N i . ( 23
)
The continuous mesh spatial complexity is constant on this sub-interval. To remove the time dependency of continuous mesh M i , we propose to use the optimal instantaneous ALE continuous mesh and property iii) of Proposition 1 which provides equality between interpolation errors. Indeed, continuous mesh M i (x(t))

x∈Ω(t)
can be mapped back to Ω(t i ) using M i,ALE (x(t i ))

x∈Ω(t i )
where φ maps Ω(t i ) onto Ω(t), and property iii) of Proposition 1 states that

E L p (M i,ALE ) = E L p (M i ) ,
thus expression of the interpolation error for each time t can be re-written at time t i . As we seek for the dynamic continuous mesh at time t i which is optimal to control the interpolation error for the whole sub-interval, we can recast Error Model [START_REF] Frey | Mesh generation. Application to finite elements[END_REF] into the following error model where the metric is independent of the time:

E i,ALE L p (M i ) = t i+1 t i Ω(t i ) trace (M i,ALE ) -1 2 (x(t i )) |H i,ALE u (x(t i ))| (M i,ALE ) -1 2 (x(t i )) p dx(t i ) dt ,
where the x(t i ) are in domain Ω(t i ). The time dependency in H i,ALE u is hidden in mapping φ and operator . . The previous expression can now be written on Ω(t i ):

E i,ALE L p (M i ) = Ω(t i ) trace (M i,ALE ) -1 2 (x(t i )) H i,ALE u (x(t i )) (M i,ALE ) -1 2 (x(t i )) p dx(t i ) ,
where the mean ALE Hessian metric is:

H i,ALE u (x(t i )) = t i+1 t i |H i,ALE u (x(t i ))| dt.
The expression of the error has exactly the same form as in Problem ( 9) in Section 2.3, thus the spatial minimization gives the same optimal metric where H i u is replaced by H i,ALE u . Then, the temporal minimization leads to the following optimal space-time ALE continuous mesh

M ALE L p = {M i,ALE L p } i=1,..,n adap : M i,ALE L p (x(t i )) = N 2 3 st n adap j=1 K j,ALE t j+1 t j τ(t) -1 dt 2p 2p+3 -2 3 t i+1 t i τ(t) -1 dt -2 2p+3 (det H i,ALE u (x(t i ))) -1 2p+3 H i,ALE u (x(t i )) ,
where the x(t i ) are in domain Ω(t i ). The ALE continuous mesh dependence in time is hidden in H i,ALE u by means of mapping φ. This way, the mesh generated at time t i is adapted to the solution at any time t > t i within the sub-interval [t i , t i+1 ] once moved with the mesh deformation displacement. Once again, we stress that preserving the number of degrees of freedom when moving the mesh is essential in this analysis.

Global fixed-point mesh adaptation algorithm for moving meshes

The optimal space-time ALE continuous mesh is designed to fit in the global fixed-point unsteady mesh adaptation algorithm described in Algorithm 1. Nevertheless, a few things need to be modified to extend this algorithm to moving mesh ALE simulations. In this section, we first recall the connectivity-change moving mesh strategy, how connectivity-changes are handled in the ALE numerical scheme, and how the ALE metric field is updated in time. Then, we specify the modifications to Algorithm 1.

Connectivity-change moving mesh strategy

The moving mesh strategy used in this paper is detailed in [START_REF] Alauzet | A changing-topology moving mesh technique for large displacement[END_REF]. The mesh adaptation framework makes us consider body-fitted simulations. Our strategy is designed to deal with large displacement of the boundaries, when the volume mesh quickly gets too distorted. The aim of this strategy is to preserve a good mesh quality all along the movement without costly remeshings, thanks to an improved mesh deformation step and an improved mesh optimization step.

The aim of the mesh deformation step is to compute a displacement for all inner vertices from the displacement of the boundaries, so that the volume mesh follows the moving geometries and thus remains valid. This is achieved by solving a linearelasticity-like PDE on the whole domain, the inside of the meshed domain being taken for a nearly incompressible elastic material and the displacement of the boundaries being the boundary conditions of the elasticity problem. To improve the CPU efficiency of this step, the number of such solutions is significantly reduced: the mesh deformation problem is solved for large time-steps, and the trajectories are constant over these large time-steps. To improve the precision of this step, higher-order trajectories are computed: two elasticity problems are solved, thus providing a speed and an acceleration to the inner vertices. In the present strategy, boundary vertices do not slide on the surfaces (mainly because it requires an accurate geometric re-projection on a CAD model), except in the specific case of simple planar surfaces such as a symmetry plane.

To preserve a good mesh quality, local mesh optimizations are performed between two mesh deformation steps, using only vertex smoothing and generalized edge/face swapping. No vertices are added or removed. Vertex smoothing consists in moving vertices close to the center of gravity of their vertex ball, and helps recovering nicely shaped elements. The swap operator changes the connectivity of the mesh, and is especially powerful in handling shear and large deformation movement. These mesh optimizations are simulated on the current mesh considering the current metric field to evaluate quality criteria, and we also analyse their impact on the mesh at the end of the mesh deformation (i.e., in the future) considering the metric field at the end of the mesh deformation. A node repositioning or an edge/face swapping is accepted and performed if and only if it satisfies the quality criteria on the current mesh and the quality criteria on the mesh at the end of the mesh deformation, see [START_REF] Alauzet | A changing-topology moving mesh technique for large displacement[END_REF] for details. For instance, a face swap operation that improves the quality of the current mesh but invalidate an element at the end of the mesh deformation will be rejected. The definition of the appropriate metric field -as the metric field evolves in time -is discussed in Section 5.3. This way, the mesh optimizations cannot affect the accuracy and the directional features of the mesh, otherwise the quality check will reject such optimizations.

These optimization steps are not performed at every solver time-step, but only when vertices have crossed a certain predetermined number of elements, considering a geometric CFL-like condition (generally every few tens or hundreds of flow solver time-steps). Moreover, quality optimizations are simulated only on elements having a quality greater than a given threshold which is generally set to 3, and they are applied only if they fulfill the considered quality criteria. Therefore, most of the time, mesh optimizations impact less that one or just a few percents of the mesh elements.

Remark 1. The connectivity-change moving mesh algorithm has been designed such that there is uniqueness of the result. If a simulation is re-run with the same data and on the same computer, then the result will be the same (the same smoothing and swapping operations will be performed). This uniqueness is also preserved in parallel.

Dealing with connectivity-change in the ALE numerical scheme

In [START_REF] Barral | Large displacement body-fitted FSI simulations using a mesh-connectivity-change moving mesh strategy[END_REF], this moving mesh algorithm has been successfully coupled to the flow solver described in Section 3.4. The numerical flow solver uses an Arbitrary Lagrangian Eulerian (ALE) numerical schemes to deal with flow simulations where the mesh is moving independently from the physical phenomena. In the ALE framework, the number of degrees of freedom and the connectivity of the mesh are supposed to be constant [START_REF] Baum | A new ALE adaptive unstructured methodology for the simulation of moving bodies[END_REF][START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF]. However, the efficiency of the moving mesh algorithm relies on the connectivity-change optimization. To devise a connectivity-change moving mesh algorithm fitting within this framework, the following choices have been made. First, the number of degrees of freedom of the simulation is kept constant. In two dimensions, an ALE formulation of the swap operator was proposed in [START_REF] Olivier | A new changing-topology ALE scheme for moving mesh unsteady simulations[END_REF], but its extension to the three dimensions case is rather tedious. In [START_REF] Barral | Large displacement body-fitted FSI simulations using a mesh-connectivity-change moving mesh strategy[END_REF], we have shown that the impact of swaps on the solution accuracy is very small when a linear interpolation method is considered. Nevertheless, a conservative P 1 -interpolation (see Section 3.2) is applied on each swap configuration in order to conserve the mass of the problem variables. Note that the CPU time over-cost of this conservative interpolation on each swap configuration is negligible because the number of swaps is small (it impacts only a small number of elements) at each mesh optimization and that mesh optimizations are only performed each tens or hundreds flow solver time-steps.

Update of the ALE metric field for optimizations

The mesh optimization procedure (smoothing and swapping) of the connectivity-change moving mesh strategy requires a proper metric field to evaluate geometric quantities and elements qualities of anisotropic adapted meshes. But, at a given fixedpoint iteration j and sub-interval i, the input is the optimal ALE continuous mesh M i,ALE L p defined on Ω(t i ). Therefore, when the mesh is deformed during the sub-interval due to the geometry displacement, the input metric is no more compatible with the moved mesh at a given time t. The use of an incorrect metric field will drive the smoothing to move vertices to a wrong location and the swaps to break the anisotropy, thus spoiling the adaptation of the mesh. As the mesh is evolving in time, the input metric field (i.e., the input continuous mesh) must also evolve in time. Consequently, we have to apply the deformation of the adapted mesh to the continuous mesh in order to maintain the consistency between the discrete and the continuous meshes.

But, in the theory developed in the previous sections, the mesh deformation correction is applied to the Hessian-metric |H i,ALE u | and not directly to the continuous mesh. Therefore, from the optimal ALE continuous mesh M i,ALE L p (x(t i ))

x∈Ω(t i )
at the beginning of sub-interval [t i , t i+1 ], we cannot directly find the deformed optimal ALE continuous mesh M i,ALE L p (x(t))

x∈Ω(t) for t ∈ [t i , t i+1 ].
Two possibilities arise from these remarks. First, several ALE metric field samples

M i,ALE L p (x(t k )) x∈Ω(t k ) k
can be precalculated for each sub-interval at the previous fixed-point iteration, and when an ALE metric field is required at time t, it is linearly interpolated in time on the current mesh from two of these pre-calculated ALE metric fields. This method has two drawbacks: it requires computing and saving all these ALE metric fields and it requires a metric interpolation stage at each mesh optimization which leads to consequent memory, I/Os and CPU time overhead. The second possibility is to modify the input ALE metric field according to the current mesh deformation to get a consistant metric field at time t. To deform the continuous mesh defined at the beginning of the sub-interval, we adopt the same reasoning as the one that leads to the optimal instantaneous ALE continuous mesh in Section 4.1, but on a reverse time frame.

Let us consider sub-interval [t i , t i+1 ], t k a time in this sub-interval and d the displacement field of the mesh between t i and t k . Here, the problem is reversed, we seek for the continuous mesh at time t k that will result in the continuous mesh at time t i once moved with displacement -d. At time t i , the continuous mesh is M i,ALE L p (the adapted mesh at time t i is generated directly using this metric field). Displacement d is the displacement of vertices between t i and t k :

x k = x i + d(x i ) and we write d (x k ) = x i -x k = -d(x i ). We define φ : Ω k -→ Ω i x k -→ x i = φ x k = x k + d (x k ) . (24)
We start from an adapted mesh which is unit for metric field M i,ALE L p at time t i and we search for the expression of the metric field at time t k for which the deformed adapted mesh at time t k is unit. Following Section 4.1, we consider an edge e k at time t k having edge e i as image by φ at time t i . They verify:

1 = e i (x i ) T M i,ALE L p x i e i (x i ) = e i φ (x k ) T M i,ALE L p x i e i φ (x k ) = ∇ k φ (x k ) T e k (x k ) T M i,ALE L p (x i ) ∇ k φ (x k ) T e k (x k ) = e k (x k ) T ∇ k φ (x k ) M i,ALE L p (x i ) ∇ k φ (x k ) T e k (x k ) .
In other words, the metric field M i,optim L p we are looking for at time t k is:

M i,optim L p (x(t k )) = ∇ k φ x(t k ) M i,ALE L p x(t i ) ∇ k φ x(t k ) T . ( 25
)
This metric field is easy to compute on the fly, because the central term is the input metric field M i,ALE L p

x(t i ) which is known, and the gradients of φ are easy to compute, since the displacement considered in φ is the opposite of the current displacement at each vertex. Note that the gradients of φ should be computed on Ω(t k ) that is to say on the current mesh.

Modifications of the global fixed-point mesh adaptation algorithm

The overall global fixed-point mesh adaptation algorithm remains unchanged: the simulation time frame is still divided into sub-intervals, and a global fixed-point strategy is used to converge the meshes and the solutions. On each sub-interval, the mesh number of vertices remains the same. However, and it is the first difference, within a sub-interval, the mesh is moved using the connectivity-change moving mesh algorithm presented in Section 5.1 and fully detailed in [START_REF] Alauzet | A changing-topology moving mesh technique for large displacement[END_REF]. One or several mesh deformation steps are performed during a sub-interval. The connectivity-change moving mesh algorithm uses mesh optimizations to make the mesh deformation efficient and robust. These mesh optimizations, smoothing and swapping operators, use the dynamic (corrected) metric field presented in Section 5.3 to be compliant with the current dynamic adapted mesh.

The second difference concerns the computation of the ALE Hessian-metric for the next adaptation step. ALE Hessianmetrics |H i,ALE u | need to be sampled in time. To compute each sample, we use the displacement given by the mesh deformation step without taking into account the smoothing optimization. In our numerical experiments, we did not see any significant difference on the solution accuracy and the quality of the adapted meshes by considering or not the smoothing in the ALE Hessian-metrics computations. This is due mainly to the fact that the mesh smoothing has only a slight influence in 3D, the mesh displacement is predominantly governed by the mesh deformation [START_REF] Alauzet | A changing-topology moving mesh technique for large displacement[END_REF]. The mesh smoothing is only a slight perturbation around the mesh deformation displacement to get better shaped elements. Nevertheless, we prefer not to consider the smoothing displacement when computing ALE Hessian-metrics because:

• the geometry displacement is clearly defined, thus the mesh deformation is converging toward a unique mesh deformation and therefore we converge toward a unique optimal ALE continuous mesh

• the smoothing correction for one node may change from one fixed-point iteration to the other, thus for each node it does not converge to a fixed correction.

The gradient of the mesh deformation displacement field is computed on the mesh position at the beginning of the sub-interval. The Hessian of the sensor is computed on the mesh current position and is mapped back to the mesh position at the beginning of the sub-interval. Then, as in the standard algorithm, the ALE Hessian-metrics are computed at each vertex.

Note that the connectivity-change moving mesh algorithm allows the mesh vertices to move to their final position given by the mesh deformation solution while it would be impossible to reach such positions without skewing mesh elements with a classical method. Thanks to the robustness of this algorithm, we are always able to calculate the ALE hessian-metrics.

Numerical examples

We now present several three-dimensional CFD simulations to illustrate the efficiency of the moving-mesh unsteady fixedpoint adaptation algorithm. The presented examples were run on 20 cores using a 2 Xeon E5-2670 v2 chip processor (10 cores at 2.5 GHz for each chip), both chips being connected by 2 QPI links with a speed of 16 GB/s.

A shock tube in expansion

The first example is a variation of the classic Sod shock-tube problem, with a tube being homogeneously expanded in one direction. This a priori simple test case shows that the refined regions follow the physical phenomena of interest, and that the anisotropy is well preserved despite the mesh being moved.

The tube initial dimensions are [0, 1] × [-0.15, 015] × [0, 1], and the gas is split in two regions: for x ≤ 0.82 the state is (ρ left , u left , p left ) = (1, 0, 1) whereas for x > 0.82 the state is (ρ right , u right , p right ) = (0.125, 0, 1). The gas is then left free to evolve, and the classic rarefaction wave, contact discontinuity and shock discontinuity appear and evolve in the tube. The tube is expanded to the right with the following movement imposed to the mesh vertices:

             x(t) = x 0 + 1.5 x 0 t y(t) = y 0 z(t) = z 0 .
The simulation is run until time T = 0.6 so that the size of the tube doubles in the x-direction. For this case, 5 fixed-point iterations and 20 sub-intervals were prescribed, and the density field is used as sensor for the adaptation. A space-time complexity N st = 80, 000 was prescribed leading to a sub-interval average spatial complexity N avg = 4, 000. The resulting adapted meshes at the last fixed-point iteration have an average number of spatial vertices equal to 24, 500. A total of n iter = 7, 075 time-steps have been performed in the last fixed-point iteration. The total number of space-time vertices at the last fixed-point iteration is N st = 174 million. For this simulation, the CPU time to compute the solution at the last fixed-point iteration (flow solver and interpolation steps on the last set of sub-interval adapted meshes) represents 43% of the total CPU time. In other words, the overhead of the iterative process to obtain the optimal adapted space-time mesh is slightly higher than a factor two. This is mainly due to the fact that the accuracy of the mesh is increased at each fixed-point iteration as the mesh-solution convergence is being established.

The meshes at different time-steps and the corresponding solution are shown in Figure 8. Not only do the waves evolve in the refined bands, but these bands also move as the tube is expanded. The adapted regions are transported together with the mesh in accordance with our algorithm, that associates metrics to moving vertices rather than to fixed positions in space. In Figure 9, we 21

show the mesh at the beginning and at the end of the fourteenth sub-interval, to emphasize the movement of the adapted regions throughout a sub-interval, together with the movement of the tube boundaries, still preserving the anisotropy. In Figure 10, we zoom on the solution at the beginning and the end of the fourteenth sub-interval to stress that the evolving physical phenomena remain inside the adaptation bands during a sub-interval, despite the mesh deformation and the fact that bands themselves are moving. We notice that the region where the shock wave evolves during the sub-interval is a lot larger than the refined band size because the refined band moves forward together with the shock wave. Therefore, the size of refined band has been reduced thanks to the ALE metric.

This case has the advantage of having a known analytic solution, which makes it possible to analyze quantitatively the effects of the adaptation. Three versions of this case are compared: the first one is a moving mesh ALE simulation without adaptation, the second one is an adaptive simulation with no moving mesh, the tube being fixed in expanded position, and the third one is a fully moving mesh adaptive simulation. In all three simulations, the same space-time complexity is prescribed. The results are presented in Figure 11 and point out that the moving mesh adapted solution matches the analytic solution as well as the fixed mesh adapted solution, both being much more accurate than the non-adapted ALE simulation in the shock, the contact discontinuity and the expansion regions.

To emphasize the efficiency of the proposed approach, a quantitive CPU analysis of this simulation is carried out by comparing the adaptive simulations with non-adaptive simulations using uniform meshes. Then, we quantify the CPU overhead due to moving meshes by comparing the adaptive simulations with and without moving meshes. Please note that in general, one does not run a moving mesh simulation when a fixed mesh simulation is possible.

Adaptive mesh versus uniform mesh simulations CPU efficiency. We compare the CPU time obtained with the global fixed-point mesh adaptation for moving mesh simulations (Algorithm 1 and Section 5.4) and simulations on uniform meshes on this shock tube in expansion problem. The CPU time in each case is analyzed with respect to the solution accuracy at final time T = 0.6, i.e., the spatial error in L 1 -norm of the solution at time T = 0.6 with respect to the analytical solution. For a reliable analysis, the comparison is done for many simulations to observe the convergence of the error with respect to the CPU time. This CPU time analysis was run on 12 cores using a Xeon E5-2697 v2 chip processor at 2.7 GHz.

All parameters of the ALE flow solver are identical for all simulations. The CPU time for the uniform mesh simulations corresponds to the CPU time of the ALE flow solver to run such simulations. Four uniform mesh simulations have been run with a number of vertices ranging from 450, 000 to 3.5 million. The CPU time in the adaptive case corresponds to the sum of all the CPU times for each stage of the adaptive algorithm (flow solver, interpolation, metric fields construction, and remeshing) for all fixed-point iterations. Fifteen adaptive mesh simulations have been run using 6, 12 and 18 sub-intervals and a theoretical average complexity of 1, 000, 2, 000, 4, 000, 6, 000 and 8, 000 leading to sub-interval adapted meshes with a size ranging from 6, 420 to 46, 204 vertices. The convergence in the adaptive case can be either observed by fixing the number of sub-intervals and increasing the theoretical complexity, see Figure 12 (left), or by fixing the theoretical complexity and increasing the number of sub-intervals, see Figure 12 (right). In these plots, we clearly see the superiority of the adaptive simulations. Indeed, the error is reduced by a factor 4 for the more accurate simulations, and we also observe a higher order of convergence for adaptive simulations. To reach the accuracy of the adaptive simulations, the uniform mesh simulations would need a lot more vertices.

A more detailed error analysis has been done in the fixed-boundary case [START_REF] Alauzet | Multi-scale anisotropic mesh adaptation for time-dependent problems[END_REF] showing the great improvement of the adaptive simulations versus non-adaptive simulations. This work also explains why, in terms of convergence, it is more advantageous to increase the number of sub-intervals 5 while performing a convergence study instead of increasing the average complexity as can be seen in Figure 12. creating complex patterns both in front of it and in its wake. We expect all these physical phenomena to be captured correctly by the multiscale adaptation algorithm. For this case, 80 sub-intervals and 5 fixed-point iterations were prescribed, and the density field is used as sensor for the adaptation. A space-time complexity of 12 million was prescribed leading to a sub-interval average spatial complexity of N avg = 150, 000. The resulting adapted meshes at the last fixed-point iteration have an average number of spatial vertices equal to 265, 000 with sizes ranging from 100, 000 to 450, 000 vertices. A total of n iter = 1, 854 time-steps have been performed in the last fixed-point iteration. The total number of space-time vertices N st at the last fixed-point iteration is 500 million. For this simulation, the CPU time to compute the solution at the last fixed-point iteration (flow solver and interpolation steps on the last set of sub-interval adapted meshes) represents 20% of the total CPU time.

The adapted meshes at different time-steps and the corresponding solutions are shown in Figure 13. Since we are interested in the interaction of the ball with the shock and the contact discontinuity, we only show the part x > 0.5 of the tube. The meshes shown are the ones at the end of the sub-intervals, i.e., the meshes that have been moved. Thanks to the anisotropic mesh adaptation, a highly resolved solution is obtained: the shock wave and the contact discontinuity are accurately captured all along Concerning the fluid simulation, the aircraft are moved at a speed of Mach 0.4, in an initially inert uniform fluid: at t = 0 the speed of the air is zero everywhere. Transmitting boundary conditions are used on the sides of the surrounding box, while slipping conditions are imposed on the two F117 bodies. After a short phase of initialization, the flow is established when the two F117s pass each other, and the density fields around the aircraft and in their wake interact. Acoustic waves are created in front of the F117s due to the instantaneous setting in motion of the aircraft.

The adaptation parameters are the following: the adaptation is performed on the density field, the time interval is divided into 50 sub-intervals, and 4 fixed-point iterations are performed. A space-time complexity 20 million was prescribed leading to a sub-interval average spatial complexity N avg = 400, 000. The adapted meshes at the last fixed-point iteration have between 360, 000 and 820, 000 vertices, for an average number of vertices equal to 745, 000. A total of n iter = 93, 875 time-steps have been performed in the last fixed-point iteration. The total number of space-time vertices N st at the last fixed-point iteration is equal to 70 billion. The total CPU time of this simulation is 130 hours, 58 hours being spent to compute the solution in the last fixed-point iteration (flow solver and interpolation steps on the last set of sub-interval adapted meshes) which represents 44% of the total CPU time. The CPU distribution between the different stages of the fixed-point mesh adaptation algorithm are given in Table 5.

Note that this simulation involves only 50 remeshing and solution interpolation between meshes at each fixed-point iteration. If we had considered a strategy where the mesh is frequently remeshed, for instance each 10 time-steps, the generation of 9, 387 adapted meshes and 9, 387 solution interpolation between meshes would have been required. It is clear that such a strategy would have been a lot less efficient and a lot less accurate that the method proposed in this paper.

The meshes and density solutions at different time-steps are shown in Figure 16. The meshes are well adapted to the solution. Thanks to the adaptation, the wakes of the F117s are captured far from the aircraft, and their interactions with the wakes are clearly visible. The anisotropy of the meshes is obvious, in front of the aircraft (corresponding to the acoustic waves), and in their wake. In Figure 15, a close up on the two aircraft is made, that shows the adapted surface mesh and solution. In this case again, a correct handling of the mesh optimization with a dynamic ALE metric is necessary to preserve the anisotropic mesh adaptation.

One can conclude from these examples results that the adaptive connectivity-change moving mesh algorithm has proved to be very efficient in deforming anisotropic adapted meshes containing anisotropic elements with high aspect ratio when large geometry displacements are involved.

Conclusion

This paper has addressed time-accurate anisotropic mesh adaptation in the case of dynamic meshes within the range of bodyfitted ALE simulations. The multiscale space-time interpolation error analysis for time-dependent problems has been extended to the case of dynamic meshes leading to the expression of the optimal ALE continuous mesh. The mesh deformation is taken into account in this expression: the ALE metric field is used to generate adapted meshes that will remain adapted once moved as required by the geometry movements. Then, the global fixed-point unsteady mesh adaptation algorithm has been coupled to a connectivity-change moving mesh algorithm. The ideas driving this algorithm remain the same, and only few modifications are required: use the connectivity-change moving mesh algorithm to move the vertices, compute and use the ALE metric field and use a dynamic ALE metric field that takes into account the mesh deformation for the mesh optimization.

Several numerical examples were presented in three dimensions, that validate the proposed approach. Adapted regions remain around the solution flow features, even if the mesh vertices are moved in other directions to follow the moving boundaries, and we have successfully managed to avoid the frequent remeshings usually performed for this kind of simulations. This results in a better accuracy at a lesser cost. Indeed, we have shown the superiority in CPU time efficiency of the adaptive algorithm with respect to uniform mesh simulations. These examples show both that our theoretical optimal ALE continuous mesh effectively leads to highly resolved solutions, and that it can be put into practice at a reasonable time cost for complex fluids problems: in 3D, with complex boundaries undergoing large displacements.

Several improvement tracks are considered. First, it is now within our reach to extend goal-oriented unsteady error estimates of [START_REF] Belme | Time accurate anisotropic goal-oriented mesh adaptation for unsteady flows[END_REF] to moving boundary problems. Second, the mesh deformation due to the moving geometry may lead us to over-refine some regions of the computational domain. For instance, if we consider an expanded tube that grows by a factor n on a subinterval, if an element size of h is expected at the end of the sub-interval, we may have to generate an element size of h/n at the beginning of the sub-interval. This issue is recurrent, for instance in piston engine simulations. In that case, we may want to adapt the sub-interval time length ∆t in order to control this over-refinement by a given factor. This will also provide a control on the time step. Finally, this paper has only discussed the case of an explicit temporal scheme. If an implicit scheme is considered with possibly arbitrary large time-steps to increase the flow solver efficiency, the size of the time-steps should be adapted in order to control the error introduced by the temporal scheme. To this end, the time step has to be introduced in the space-time error analysis. A first work in this direction has been proposed in [START_REF] Coupez | Adaptive time-step with anisotropic meshing for incompressible flows[END_REF].

Acknowledgments

This work was partially funded by the Airbus Group Foundation and was partly done in the MAIDESC ANR project which is supported by the french ministry of Research under contract ANR-13-MONU-0010.

Figure 1 :

 1 Figure 1: London tower bridge case: initial surface mesh and adapted surface mesh at dimensionless time t = 21.

Figure 2 :

 2 Figure 2: London tower bridge case: zoom on details of the mesh and the solution at dimensionless time t = 21 in the region where the explosion originated.

Figure 3 :

 3 Figure 3: London tower bridge case: snapshots of adapted meshes and solution at dimensionless times t = 0.63, t = 8.6, t = 17.2, t = 25.

10

 10

Figure 4 : 1 (

 41 Figure 4: Left, an initial uniform-size mesh. Middle and right, resulting meshes after applying displacements d n+1 1 (midlle) and d n+1

2

 2

Figure 7 : 1 (left) and u n+1 2 (

 712 Figure 7: Interpolation error convergence curves for sensor u n+1 1 (left) and u n+1 2 (right).

H

Figure 8 :

 8 Figure 8: Adapted meshes and density solutions for the expanding shock tube case at time 0, 0.5 and 0.66. Cuts into the volume mesh are made along the plane y = 0.

Figure 9 :

 9 Figure 9: Mesh at the beginning (top)t = 0.42 -and end (bottom)t = 0.45 -of the fourteenth sub-interval. One can see the refined bands are moved forward as the tube is being expanded.

Figure 10 :

 10 Figure 10: Zoom on the mesh and solution at the beginning (top)t = 0.42 -and at the end (bottom)t = 0.45 -of the fourteenth sub-interval. One can see the two waves move forward inside the adapted bands while the refined bands are moving forward.

Figure 11 :

 11 Figure11: Density solution extraction along a line inside the domain for a moving mesh ALE simulation without adaptation, an adaptive simulation with no moving mesh, and moving mesh adaptive simulation. The extractions are compared to the analytic solution on the whole simulation interval (first graph) and on zooms on the different features.

Figure 13 :

 13 Figure 13: Adapted meshes and density solutions for the moving ball in a shock tube test case, at times 0.07, 0.12, 0.19 and 0.44. Cuts into the volume mesh are made along the plane y = 0.

Figure 14 :

 14 Figure 14: Nosing-up F117 test case: adapted meshes (view from the top) and local Mach number isolines at different time steps.

Figure 15 :

 15 Figure 15: Close-ups on the surface mesh and solution of the two F117s.

Figure 16 :

 16 Figure 16: Adapted meshes and density solutions for the two F117s test case at different time-steps. Cuts into the volume mesh are made along the plane y = 0.

Table 1 :

 1 Adapted meshes elements quality for sensor u n+1

			n+1 re f	H n+1 ALE with d 1	H n+1 ALE with d 2
	Average quality		1.07		1.20		1.23
	Worst quality		3.28		23.7		66.0
	1.00 < Q < 2.00	127303 99.92 %	122473 96.69 %	121585 96.02 %
	2.00 < Q < 3.00	97	0.08 %	3283	2.59 %	3641	2.88 %
	3.00 < Q < 4.00	6	0.00 %	611	0.48 %	843	0.67 %
	4.00 < Q < 5.00	0	0.00 %	170	0.13 %	244	0.19 %
	5.00 < Q < 10.00	0	0.00 %	121	0.10 %	258	0.20 %
	10.00 < Q < 50.00	0	0.00 %	7	0.01 %	49	0.04 %
	50.00 < Q < 100.00	0	0.00 %	0	0.00 %	2	0.00 %
				1	and a spatial complexity of 50, 000.

Table 2 :

 2 Adapted meshes elements quality for sensor u n+1

			H n+1 re f	H n+1 ALE with d 1	H n+1 ALE with d 2
	Average quality		1.06		1.19		1.21
	Worst quality		13.8		20.5		24.4
	1.00 < Q < 2.00	126406 99.95 %	121522 96.87 %	121146 96.36 %
	2.00 < Q < 3.00	59	0.05 %	3054	2.43 %	3631	2.89 %
	3.00 < Q < 4.00	5	0.00 %	594	0.47 %	652	0.52 %
	4.00 < Q < 5.00	0	0.00 %	152	0.12 %	177	0.14 %
	5.00 < Q < 10.00	0	0.00 %	118	0.09 %	109	0.09 %
	10.00 < Q < 50.00	1	0.00 %	4	0.00 %	5	0.00 %
				2	and a spatial complexity of 50, 000.

Table 5 :

 5 Percentage of the total CPU time for the different stages of the fixed-point mesh adaptation algorithm on the last three examples.

The generation of a mesh is usually more costly than a few flow solver time-steps.Preprint submitted to Journal of Computational PhysicsDecember 4,

Unless a flow solver with dynamic data fully coupled with the adaptive remesher is used.

Here, the gradient is not the Jacobian, i.e., for an arbitrary vector field f = ( f 1 , .., f k ), its gradient matrix is ∇f = ∂ f j ∂x i i j .

Even if Ω n Ω n+1 in general, we assume there is diffeomorphism mapping φ from Ω n to Ω n+1 , .i.e., these two spatial domains can be mapped one onto the other.

(right), and displacements d 1 (top) and d

(bottom).

(right), and displacements d 1 (top) and d 2 (bottom).

Changing the number of sub-intervals optimizes (adapts) the space-time mesh in the time direction.

Convergence curves of the spatial error of the final solution with respect to the CPU time for uniform mesh simulations (black) and adaptive mesh simulations (red, green, blue, pink and cyan). The convergence in the adaptive case can be either observed by fixing the number of sub-intervals and increasing the theoretical complexity (left) or by fixing the theoretical complexity and increasing the number of sub-intervals (right).

CPU overhead of moving mesh simulations. We analyze the CPU overhead of the ALE solver with respect to the solver on fixed meshes. This overhead is due to the update of the finite volume cells, the computation of the mesh velocities, and the extra terms in the ALE fluxes. When an explicit time-stepping is considered, if the cost of the solver on fixed mesh is 1 then the cost of the solver on moving meshes is around 1.75. There is also an overhead due to mesh deformation: the elasticity solutions, the mesh optimizations, moving the mesh and checking its validity. For this case, if the cost of the solver on fixed mesh is 1 then the cost of the mesh deformation is 0.2. Moreover, in this very specific case, a total of 7, 075 time-steps have been performed for the moving mesh adaptive simulation whereas 3, 721 time-steps have been performed for the fixed mesh adaptive simulation. Indeed, as the tube is expanding and to preserve the solution accuracy, the ALE metric field prescribes a mesh with a size smaller than required in the x-direction at the beginning of the sub-intervals because this size is going to grow during the sub-intervals. Thus, a larger number of time-steps are performed for the moving mesh adaptive simulation. For the whole simulation, adding all these overheads, the cost of the ALE flow solver is 3.31 times the cost of the flow solver on fixed meshes, see Table 3.

We can then analyze the global cost of the adaptive simulations. In Table 3, we observe that there is no overhead for the interpolation or the remeshing parts. Thus, the adaptive moving mesh case is 2.85 times the cost of the adaptive fixed mesh case. The percentage of the total CPU time for the different stages of the fixed-point mesh adaptation algorithm for both cases are given in Table 4. It is clear that for these cases most of the CPU time is spent in the flow solver.

A moving ball in a shock tube

In the second example, we consider again a Sod shock tube, but instead of expanding the tube, we add a moving ball inside it, that interacts with the shock and the contact discontinuity.

The dimensions of the tube are [0, 1] × [-0.15, 0.15] × [-0.1, 0.1]. At initial time, the gas is split in two states: for x ≤ 0.5 the state is (ρ left , u left , p left ) = (1, 0, 1) whereas for x > 0.5 the state is (ρ right , u right , p right ) = (0.125, 0, 1). A ball of radius r = 0.02 is immersed in the gas, its center being in position (0.75, 0, 0). The tube is fixed, while the ball has a constant speed v ball = -0.3 e x . The simulation is run until final time T = 0.5. After a while, the ball goes through the shock and the contact discontinuity, the simulation, and more complex features around the ball and in its wake are represented in detail. It is interesting to note that the shock is perfectly meshed as soon as the ball has gone through it. A closer look at the surface of the ball also shows that it is adapted with highly anisotropic elements, and the anisotropy close to the surface is well preserved.

For this simulation, allowing mesh optimization (smoothing and edge/face swapping) is mandatory, otherwise, the shearing created by the advance of the ball in the volume mesh quickly creates invalid elements. Consequently, it is crucial to compute the quality criteria that trigger optimization with regards to the actual metric of the mesh. In this case, it is clear in figures shown that the connectivity-change optimization did not break the anisotropy in the adapted regions, although nearly 4.5 millions of swaps were performed during the last fixed-point iteration, i.e., on average 55, 000 swaps per sub-interval are performed. To give an idea of the impact of the swaps on the numerical scheme, it represents nearly 2, 500 swaps per time-step as 1, 854 time-steps have been performed at the last fixed-point iteration. These 2, 500 swaps have to be compared to the number of elements in the meshes which is 1.6 millions on average for each sub-interval, i.e., one swap for 640 elements per time-step. The CPU distribution between the different stages of the fixed-point mesh adaptation algorithm are given in Table 5.

Nosing-up F117 aircraft

The third example is a subsonic notional F117 aircraft geometry nosing up, that creates a vortical wake. An inflow of air at Mach 0.4 arrives in front of the aircraft, initially in horizontal position, that noses up, stays up for a while, then noses down. In this example, the aircraft rotates around its center of gravity. Let T = 1s be the characteristic time of the movement and θ max = 20 o the maximal angle reached, the movement is defined by its angle of rotation, of which the evolution is divided in 7 phases:

Phase (i) is an initialization phase, during which the flow around the aircraft is established. Phases (ii) and (iii) are respectively phases of accelerated and decelerated ascension. Vortices start to grow behind the aircraft, and they expand during phase (iv), where the aircraft stays in upward position. Phases (v) and (vi) are phases of accelerated and decelerated descent, the vortices start to move away and they slowly disappear in phase (vii). Free-stream conditions are imposed on the faces of the surrounding box, and slipping conditions on the aircraft. The time simulation interval [0, 5] was divided into 96 sub-intervals, and 5 fixed-point iterations are performed. The adaptation sensor is the local Mach number. A space-time complexity of 48 million was prescribed leading to a sub-interval average spatial complexity N avg = 500, 000. The adapted meshes at the last fixed-point iteration have between 100, 000 and 1, 400, 000 vertices, for an average number of vertices equal to 841, 000. A total of n iter = 38, 980 time-steps have been performed in the last fixed-point iteration. The total number of space-time vertices N st at the last fixed-point iteration is equal to 35 billion. The total CPU time of this simulation is 95 hours, 31 hours being spent to compute the solution in the last fixed-point iteration (flow solver and interpolation steps on the last set of sub-interval adapted meshes) which represents 32% of the total CPU time. The CPU distribution between the different stages of the fixed-point mesh adaptation algorithm are given in Table 5.

Views of the adapted meshes and the corresponding solutions are displayed in Fig. 14. The vortical wake is propagated far from the aircraft, and the patterns of the vortices are highly resolved. The observation of the adapted meshes shows that they are actually refined only in the vicinity of the wake, and with such a precision that one can see the vortices being created, evolving and vanishing just by looking at the meshes. A view at several meshes within the same sub-interval shows that the mesh evolves continuously following the physical phenomena. In the wake far from the aircraft, the elements are highly anisotropic, whereas they have to be isotropic in the area of vorticity due to the characteristic of the physical phenomena.

Two F117 aircraft flight paths crossing

This case is an example of moving mesh simulation, described in [START_REF] Barral | Two mesh deformation methods coupled with a changing-connectivity moving mesh method for CFD applications[END_REF]. It models two notional F117 aircraft geometries having crossing flight paths, translating and rotating. This problem is difficult in terms of mesh movement and it illustrates the efficiency of the connectivity-change moving mesh algorithm in handling large displacements of complex geometries without any remeshing. When both aircraft cross each other, the mesh encounters large shearing due to the opposite flight directions. The connectivity-change mesh deformation algorithm handles easily this complex displacement thanks to the mesh local reconnections.