Tetrahedral Remeshing in the Context of Large-Scale Numerical Simulation and High Performance Computing

. The purpose of this article is to discuss several modern aspects of remeshing, which is the task of modifying an ill-shaped tetrahedral mesh with bad size elements so that it features an appropriate density of high-quality elements. After a brief sketch of classical stakes about meshes and local mesh operations, we notably expose (i) how the local size of the elements of a mesh can be adapted to a user-deﬁned prescription (guided, e


Introduction 2. A few aspects of tetrahedral remeshing
2.1.General facts about meshes and notation 2.2.Basic stakes about remeshing 2.3.Goal-oriented mesh adaptation 2.4."Lagrangian" mesh deformation 2.5.Implicit domain remeshing 2.6.Parallel remeshing 3. h-adaptive RANS and hybrid RANS/LES simulations of a nozzle with the Discontinuous Galerkin method using unstructured meshes and isotropic mesh adaptation 3.1.Context and motivation 3.2.Mesh adaptation algorithm 3.3.Numerical example: simulation of the subsonic turbulent nozzle jet flow 4.An application of implicit domain meshing in shape optimization 15 4.1.A few generalities about shape and topology optimization 16 4.2.Body-fitted shape and topology optimization 16 4.3.Numerical example: optimization of an elastic mechanical system 17 5.Conforming remeshing using level set discretization for geophysical inverse problem 18 5.1.Context and motivation 18 5.2.Numerical example: sensitivity of a two-phase Darcy flow with respect to the existence and the geometry of fractures 21 6.Dynamic and parallel mesh adaptation for premixed flame fronts and liquid/gas interfaces 23 6.1.Dynamic mesh adaptation for material interfaces 23 6.2.Numerical example: application to the lean-premixed PRECCINSTA burner 25 6.3.Numerical example: application to the droplet impact on a cone 27 References 28

Introduction
Since the early days of scientific computing, simplicial meshes (that is, meshes composed of triangles in 2d, or tetrahedra in 3d) have been raising a constantly renewed interest as a prominent means to represent and process complex domains.For instance, they have been used extensively to account for a shape in the perspective of its visualization; more recently, STL mesh files have become one of the most widespread formats under which shapes are supplied to 3d printing machines.Furthermore, and closer to the scope of the present article, the most popular numerical simulation frameworks for physical phenomena (namely, the finite element method, or the finite volume method) crucially rely on a meshed discretization of the computational domain.We refer for instance to [21] for a general presentation of several stakes and applications of meshing.
In line with this omnipresence of meshes in the numerical treatment of shapes, the issues of mesh generation and mesh processing have received a tremendous amount of attention in the literature, resulting in various efficient algorithms and software packages, such as Meshlab [32] Gmsh [59], CGAL [97], Tetgen [93], to name a few free and open-source instances.
Despite these numerous investigations and achievements, meshing is still a thriving field for academic and industrial research: some old and major challenges remain -in particular, it is still unfortunately difficult to create a valid surface triangulation of the boundary of a complex 3d domain, then to fill the volume with a valid tetrahedral mesh -while modern applications suggest new and promising directions for investigations.
Many such applications fit in with the context of remeshing, where one aims to modify an existing, valid, albeit ill-shaped and badly sampled mesh, so that it better comply with given requirements, such as a good element quality, an adapted element density, etc.The purpose of this article is to discuss and illustrate several such aspects of tetrahedral remeshing.These, as well as the proposed means to address them, are suggested by the knowledge of the authors, without ambition of exhaustivity.
(1) The primary target of remeshing is to adapt the local size of a mesh (i.e. its element density) to the geometric features of the represented domain, or to an error estimate associated to a numerical simulation of interest.(2) When a physical simulation implies an evolution of the underlying domain (as it is often the case for instance in computational fluid dynamics, or in shape optimization), suitable combinations between Lagrangian strategies and remeshing techniques make it possible to efficiently account for this motion.(3) The recent development of Eulerian interface-capturing methods -such as the celebrated level set method -has made it quite popular to represent a shape Ω implicitly, i.e. as the negative subdomain of a scalar function φ : D → R, defined on (a fixed mesh of) a given computational domain D. In this context, one may need an exact mesh of Ω, which calls for a methodology for constructing a mesh of such an implicit domain.(4) When the size of the computational mesh is so large that it can barely be stored in memory, it is crucial that the remeshing process be carried out in a parallel fashion.
Understandably enough, these issues find particularly relevant applications in the recent context where the sustained increase in computational power and the advent of high-performance computing herald high resolution, very accurate numerical simulations of complex physical phenomena.The presentation of this article is guided by this perspective of large-scale simulations.In particular, it mainly takes place in the three-dimensional context, which is on every aspect more involved than its two-dimensional counterpart, to which we shall occasionally refer for pedagogical purposes, or as a preliminary step towards a future three-dimensional implementation.We discuss the use of the presented remeshing features in the context of applications involving high performance computing; all these features are implemented in the open-source software mmg [1], which is used consistently in this exposition.A thorough technical description of this environment can be found in the article [37], see also [36,43].This article is organized as follows.In the next Section 2, we briefly review some basic material about meshes and remeshing, before describing a little more precisely some modern applications in different contexts of use.We then present several applications of these methods, and variants of them dedicated to the context of high-performance computing: in Section 3, an error estimate based adaptive remeshing strategy is implemented in the context of unsteady fluid mechanics simulations.The next Section 4 arises in the context of shape and topology optimization, where an implicit domain meshing methodology is carried out to track the motion of the optimized domain, while allowing for an explicit, meshed description of the latter throughout the optimization iterations.The application of this idea of meshing implicitly-defined interfaces is then applied in the field of geoscience in Section 5, to the construction of an explicit and accurate meshed representation of geophysical interfaces.Finally, in Section 6, we exemplify how a parallel implementation of mesh adaptation techniques makes it possible to track efficiently the motion of complex fluid interfaces, such as a turbulent flame front or a liquid droplet, with a very high resolution.

A few aspects of tetrahedral remeshing
In this section, we present the remeshing features at stake in this article.After a short summary of basic concepts about meshes in Section 2.1 and a general presentation of remeshing in Section 2.2, we describe in Section 2.3 the issue of mesh adaptation.In the next Section 2.4, the problem of mesh displacement is discussed; in Section 2.5 we broach the idea of isosurface discretization before finally dealing with parallel remeshing in Section 2.6.

General facts about meshes and notation
Let us start by introducing briefly the needed background material about tetrahedral meshes in this article; for these issues, we refer to e.g.[20,27,55,68,98] for further details.
Formally, let Ω ⊂ R 3 be a bounded domain; a mesh of Ω is a collection T = {T i } i=1,...,N T of open tetrahedra, accounting for a covering of Ω, in the sense that In addition, it is assumed that • T is valid : the open simplices T i are mutually disjoint: • T is conforming: each intersection T i ∩ T j , i = j, is reduced to either a vertex, an edge, or a face of the mesh.
A mesh T additionally bears information about a surface triangulation S, that is, a collection {S j } j=1,...,N S of triangles S j ⊂ R 3 accounting for one or several pieces of surface.In the most simple instances, S is made of the external faces of the tetrahedra T ∈ T , as a (approximate) representation of the boundary of Ω.However, in some applications, Ω comprises several subdomains, whose tetrahedra are identified with different labels.In such cases, S also contains the triangles of the corresponding inner boundaries, see Fig. 5 (right) below for an example.
Usually, the creation of a mesh T of Ω starts from the datum of a surface triangulation S of the boundary ∂Ω (which is often supplied by a CAD software).Thence, several strategies are available to fill the volume of Ω with tetrahedra conforming to this triangulation, the perhaps most efficient ones being based on the constrained Delaunay algorithm.Let us mention that, in spite of having been extensively addressed for decades, mesh generation is still a delicate issue when very intricate shapes are considered.
Beyond the aforementioned mild requirements, it is often desirable to evaluate more closely "how well" a mesh T lends itself to accurate numerical simulations.This feature can be appraised in at least two independent ways, which are illustrated on Fig. 1.
• Most numerical methods experience trouble when the mesh T contains very flat, nearly degenerate elements.For instance, such configurations are well-known to increase dramatically the condition number of finite element systems based on this mesh, thus slowing down iterative matrix solvers, see e.g. the books [31,49] about this classical issue.Several quantities are used in the literature to discriminate "well-shaped elements" T i (i.e.those that are close from being regular tetrahedra), from "ill-shaped" ones (i.e. that are nearly degenerate); for instance, the following measure of the quality of a tetrahedron T is quite popular in the literature: , where e i , i = 1, . . ., 6 are the edges of T , Vol(T ) is its volume, and α is a normalization factor.This quantity Q(T ) equals 1 when T is regular, and it is close to 0 when T is nearly degenerate.
• Another crucial feature in the evaluation of the quality of T is related to an issue which we have overlooked so far.Often, the domain Ω of interest is smooth, and has a curved boundary ∂Ω, while the faces of the surface triangulation S, intended as a discrete approximation of the latter, are planar.Hence, it should be required from S that it be a close approximation of the smooth surface ∂Ω up to a certain tolerance, for instance in terms of the Hausdorff distance between S and ∂Ω.
Summarizing, we shall expect from a mesh T of Ω that each tetrahedron T i have quality Q(T i ) close to 1, and that the surface triangulation S be a "close" approximation of the continuous surface ∂Ω.

Basic stakes about remeshing
As suggested by the name, remeshing assumes the datum of a tetrahedral mesh T of Ω which is valid and conforming, but may still be unsuitable for computation: as discussed in Section 2.1, T may suffer from poor element quality, or its element density may be inadequate, in the sense that curved regions of the boundary ∂Ω (or other, internal surfaces) are represented by too few elements.Remeshing aims to modify T into a high-quality mesh T of Ω, whose element density is well-tailored to its geometric features.
This objective is usually achieved thanks to a series of local operations, which are applied iteratively; these are briefly described below, and illustrated on Fig. 2 in the (simpler) two-dimensional context.
• Edge split: When an edge pq is "too long", a new vertex m is inserted in the mesh T ; pq is replaced by the two edges pm and mq and the connections of T are updated accordingly.• Edge collapse: When an edge pq is "too short", its endpoints p and q are merged and the connections of the mesh T are updated accordingly.• Edge swap: An edge pq is removed from the mesh T and the "shell" of pq, consisting of all the tetrahedra sharing this edge, is adequately reconnected.• Vertex relocation: A vertex p is moved slightly, while all its connections remain unaltered.
Each of these operators exists under two different versions, depending on whether it is applied to a surface configuration, involving tetrahedra bearing surface triangles S j ∈ S, or to an internal one.For instance, when an internal edge pq of T is split, the inserted point m is usually chosen as the midpoint of p and q; on the contrary, when pq ∈ S is a boundary edge, m is rather placed on the continuous surface ∂Ω.In practice, since this "ideal" surface is unknown, the position of m is inferred from those of p and q, and other associated geometric quantities (such as the normal vectors at p and q).
In the above description, the criterion whereby an edge pq is deemed to be "too long" (or "too short") may depend on the situation.For instance, pq may be "too long" with respect to a user's prescription for the size of the elements of the mesh (see the next Section 2.3 about this practice) or, when pq ∈ S is a surface edge, it may be considered to be "too long" with respect to the curvature of ∂Ω, as its size entails a too coarse geometric approximation of ∂Ω.
Last but not least, let us emphasize that the use of the above remeshing operators should be carefully monitored: several checks are in order so as to prevent the emergence of invalid configurations.

Goal-oriented mesh adaptation
As we have seen, the primary objective of remeshing T is to produce a mesh T with fine element quality, which is a close approximation of the underlying domain Ω.In addition to these requirements, one usually demands that T comply with a user-defined local size prescription, which may be either isotropic or anisotropic: • In the former case, the size prescription is encoded as a size map h : Ω → R, which is often supplied at the vertices of T and interpolated from these data when necessary: h(x) is the desired size for the edges near the point x ∈ Ω. • In the latter case, the local size is imposed under the form of a metric tensor M : Ω → R 3×3 : the eigenvalues of the symmetric, positive definite matrix M (x) encode the desired size near x ∈ ∂Ω, in the directions of the corresponding eigenvectors; see [99] about this Riemannian framework and [70] for an interesting continuous paradigm based on this idea.This extra ingredient is incorporated into the general remeshing framework of Section 2.2 via the calculation of the length of an edge pq when deciding whether it is "too long" or "too short".For instance, the length r o i b 6 5 / q 0 q Q Q 0 S Y j A e U j y l m S j n r s 6 4 0 i a p d 9 t Z W + T f F l K j c s 4 y b 4 l 3 e k g Y 8 N 8 7 5 o F 0 2 r B O j 0 q y W 6 g f Z q P P Y w z 6 O a J 6 n q O M S D b S U 9 y O e 8 K x d a J 6 W a O k n V c t l m l 3 8 W N r D B 6 U N j 4 M = < / l a t e x i t > m < l a t e x i t s h a 1 _ b a s e 6 4 = " I o m 4 R N T n i S s T 7 N N P D P O m v 3 N e 2 I 8 9 d 3 G 9 P d z r 5 B Y h V t i / 9 J N M / + r 0 7 U o D H B m a g i o p s Q w u j q W u 2 S m K / r m 9 p e q F D k k x G n c p 7 g k z I x y 2 m f b a F J T u + 6 t Z + J v J l O z e s / y 3 A z v + p Y 0 Y P f n O G d B 6 7 D q n l S P G s e V 2 l 4 + 6 i J 2 s I s D m u c p a r h E H U 3 j / Y g n P F s X l r B S K / t M t Q q 5 Z h v f l v X w A V G y j 1 8 = < / l a t e x i t > p < l a t e x i t s h a 1 _ b a s e 6 4 = " n u U S J t T l m M M j x i b + g 7 p F 3 X s A n t l a f Q 6 o B O i e j l p L S x S 5 q U 8 j h h d Z q t 4 7 l 2 V u x v 3 m P t q e 4 2 o r 9 v v G J i J a 6 J / U s 3 y f y R N T n i S s T 7 N N P D P O m v 3 N e 2 I 8 9 d 3 G 9 P d z r 5 B Y h V t i / 9 J N M / + r 0 7 U o D H B m a g i o p s Q w u j q W u 2 S m K / r m 9 p e q F D k k x G n c p 7 g k z I x y 2 m f b a F J T u + 6 t Z + J v J l O z e s / y 3 A z v + p Y 0 Y P f n O G d B 6 7 D q n l S P G s e V 2 l 4 + 6 i J 2 s I s D m u c p a r h E H U 3 j / Y g n P F s X l r B S K / t M t Q q 5 Z h v f l v X w A V G y j 1 8 = < / l a t e x i t > p < l a t e x i t s h a 1 _ b a s e 6 4 = " n u U S J t T l m M M j x i b + g 7 p F 3 X s A n t l a f Q 6 o B O i e j l p L S x S 5 q U 8 j h h d Z q t 4 7 l 2 V u x v 3 m P t q e 4 2 o r 9 v v G J i J a 6 J / U s 3 y f y j O j L m I P + z i k e Z 6 i j k s 0 0 N L e j 3 j C s 3 V h h Z a w 8 s 9 U q 2 A 0 u / i 2 r I c P X a K P c g = = < / l a t e x i t > s < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 X e v 0 j O j L m I P + z i k e Z 6 i j k s 0 0 N L e j 3 j C s 3 V h h Z a w 8 s 9 U q 2 A 0 u / i 2 r I c P X a K P c g = = < / l a t e x i t > s < l a t e x i t s h a 1 _ b a s e 6 4 = " X o d 6 H j e N y 5 a x 8 c n N a q h 5 k V 5 3 H H v Z x R P d 5 j i q u U E P d d v m M F 7 w 6 1 4 5 y x s 7 k W + r k M s 8 u f g z n 6 Q s S / J I Q < / l a t e x i t > • < l a t e x i t s h a 1 _ b a s e 6 4 = " w Z U 2 C 5 s U U W t D 8 + 8 l / A S Z 9 n P l L e s = " > A A A C y n i c j V H L S s N A F D 3 G V 6 2 v q k s 3 w S K 4 C k l b t e 4 K b l y 4 q G A f 0 B Z J 0 m k N n S Z h M h F K c e c P u N U P E / 9 A / 8 I 7 Y y p K E Z 2 Q 5 M 6 5 5 5 y Z e 6 8 X 8 y C R t v 2 6 Y C w u L a + s 5 t b y 6 x u b W 9 u F n d 1 m E q X C Z w 0 / 4 p F o e 2 7 C e B C y h g w k Z + 1 Y M H f s c d b y R u c q 3 7 p j I g m i 8 F e e i V b 7 d A q n V 5 D S x C F p I u I J i t V p p s 6 n 2 l m h v 3 l P t a e 6 2 4 T + X u Y 1 J l T i l t C / d D P m f 3 W q F o k B q r q G g G q K N a K q 8 z O X V H d F 3 d z 8 V p U k h 5 g w F f c p L y j 2 t X L W Z 1 N r E l 2 7 6 q 2 r 8 2 + a q V C 1 9 z N u i n d 1 S x r w 3 D j n g 2 b J c k 6 s 8 l W l W K t k o 8 5 h H w c 4 o n m e o o Y L 1 N H Q V T 7 i C c / G p S G M i T H 9 p B o L m W Y P P 5 b x 8 A F v 6 5 J G < / l a t e x i t > • < l a t e x i t s h a 1 _ b a s e 6 4 = " p 3 K u B I 9 j W D 1 q S q 7 e N y 7 Figure 2. (a) Split of the "long" edge pq: a new point m is inserted in T and the triangles sharing this edge are subdivided accordingly; (b) collapse of point q onto p: the two red triangles disappear, and q is replaced by p in all the other triangles connected to q (in blue); (c) swap of the edge pq (in red): the connection pq is traded for the alternative configuration, featuring the edge rs (in blue); (d) relocation of the vertex p while maintaining all the connections in the mesh.
Whether it is isotropic or anisotropic, the size prescription may stem from various considerations.Here are two examples: • Geometric error estimate: The local size of the edges of T should guarantee that the geometric approximation of the smooth boundary ∂Ω of Ω by the surface triangulation S is accurate enoughfor instance, that the Hausdorff distance between ∂Ω and S is small; see Fig. 3 for an example.• A priori or a posteriori error estimate: In the perspective of reducing the computational effort of physical simulations, the local size should be adapted to a surrogate quantity for the error ||u − u h || between the solution u to a partial differential equation posed on Ω and its finite element approximation u h .This surrogate quantity may depend on u (a priori estimate) or u h (a posteriori estimate); see for instance [49] for an introduction to this practice.Eventually, let us observe that while isotropic remeshing is by now a very popular idea for tackling realistic and industrial applications, the use of an anisotropic size prescription seems more limited, since even a slight misalignment of the resulting, very stretched elements may jeopardize with the accuracy of the numerical computation, see e.g.[84,92] about this issue.

"Lagrangian" mesh deformation
Many time-dependent physical phenomena imply an evolution of the simulation domain Ω; for instance, Ω may represent a volume filled by a fluid whose velocity is predicted by the resolution of the Stokes, or the Navier-Stokes equations.
In this spirit, various applications demand to realize the motion of a domain Ω n , equipped with a mesh T n , according to a velocity field V n : Ω n → R 3 , and to obtain a mesh T n+1 of the next domain Ω n+1 := (Id + V n )(Ω n ).Here and throughout this section, the superscript n refers to discrete time.As we have (left) initial mesh with size map, calculated on the basis of the geometric approximation; (middle) adapted mesh with respect to this size prescription; (right) adapted mesh with respect to a finer geometric approximation requirement.mentioned, in practice, V n results from a physical simulation, conducted on the mesh T n ; for the purpose of this section, we suppose that it is given, as a (time independent) vector field defined at the vertices of T n .
Realizing this motion is a highly challenging task, which has been extensively considered in the literature, see for instance [4,11,35,44,72].The perhaps most intuitive strategy consists in intertwining steps where each vertex p ∈ T n is pushed in the direction of the vector V n (p) as long as the resulting mesh is validwhich quickly deteriorates the quality of T n -, with remeshing steps for improving the quality of the resulting mesh, thereby allowing to reiterate the process.One possible procedure based on this principle outlines as follows, see also Fig. 4: (1) Find the largest number 0 < t * ≤ 1 such that moving each vertex p ∈ T n to p + t * V n (p) results in a valid mesh; (2) Apply all or just one subset of the remeshing operations of Section 2.2 in order to improve the quality and the density of the resulting mesh; (3) If t * < 1, go back to (1) by replacing the vector field A number of additional heuristics may help this process, postponing the emergence of overlapping elements, i.e. allowing t * to be as close to 1 as possible during the first step.One popular practice relies on the fact that only the values of the velocity field V n (p) at those vertices p in the surface part S n of T n have an influence on the geometry of Ω n+1 , so that the motion of the internal vertices of T n can actually be chosen arbitrarily.Among such possibilities, extending the values of V n on S n to the inner nodes of T n by solving a linearized elasticity system with Dirichlet data V n on S n tends to produce a motion with little compression, thus mitigating the trend of elements to degenerate."Reasonably large" displacements of Ω can be realized owing to such Lagrangian strategies.Yet, a number of situations are difficult to handle, especially when the considered motion implies topological changes (e.g. the merger of two holes).A more robust strategy, combining remeshing techniques with an implicit representation of the evolving domain, is presented in the next Sections 2.5 and 4.

Implicit domain remeshing
In several applications, the domain Ω of interest is not readily defined by a meshed representation.Rather, a large computational box D is introduced, which is equipped with a mesh T , and Ω is defined as the negative subdomain of a scalar "level set" function φ : D → R (in practice, defined at the vertices of T ), that is: This type of representation is ubiquitous in e.g.geometric modeling or shape reconstruction.It has raised a tremendous enthusiasm within the numerical simulation community as the pivotal ingredient of the celebrated level set method [82], devoted to the description of the motion of a domain Ω(t) ⊂ D according to a velocity field V (t, x).Indeed, introducing a level set function φ(t, •) for Ω(t) (i.e.(2.1) holds at any time Figure 4. Evolution of a domain Ω (an elastic cantilever beam) via a velocity field V n (supplied by the resolution of a shape optimization problem); the vertices of the mesh T n are displaced according to V n .(Left) initial shape and associated deformation field; (middle) at iteration 150, the mesh becomes very stretched; the resolution of the linear elasticity equation is very inaccurate, and the deformation cannot be applied lest that the mesh becomes invalid; (right) after quality-oriented remeshing, the deformation process is able to go on.t > 0), this motion translates into the following advection-like equation: Hence, a geometric evolution problem is reformulated as a partial differential equation posed on a fixed domain D, equipped with the fixed mesh T ; see for instance [81,91] for comprehensive introductions to the level set method.Depending on the application, it may be of utmost importance to recover a mesh of Ω from the data of D and φ; this is the case when for instance Ω is intended as the domain of physical phenomenon, whose accurate numerical simulation is desired.
It is actually possible to modify the mesh T of D into a new (valid, conforming) "body-fitted" mesh T of D where both Ω and its complement D \ Ω are explicitly discretized.This may be achieved within two steps, as illustrated on Fig. 5: (1) Discretize the isosurface {x ∈ D, φ(x) = 0} into the mesh T .This stage is simple and relies the marching tetrahedra algorithm [45], as a variant of the well-known marching cubes algorithm [69]: in a nutshell, each tetrahedron T ∈ T crossed by this isosurface is subdivided according to a pattern.This ends up with a valid, conforming mesh T temp of D, where Ω is explicitly discretized, but which generally has very bad quality.(2) Remesh T temp into a fine quality mesh T of D, where Ω is explicitly discretized.

Parallel remeshing
As the range of applications of physical simulations grows wider, and more and more realistic applications are targeted, it is natural that very large meshes need to be handled and processed.This raises the issue of modifying such a large mesh T in a parallel fashion.
In this direction, two quite different paradigms can be thought of.On the one hand, according to shared memory strategies, the whole mesh data are stored on one single node; they are shared by the different processes or threads, and the remeshing operations are carried out in parallel.On the other hand, distributed memory strategies advocate to divide T into several regions with shared interfaces, which are independently remeshed on the different processors.Each of these two strategies faces new, specific issues with respect to sequential remeshing algorithms, such as the prevention of "parallel data races" (i.e.multiple processes trying to adapt the same mesh entities) in the implementation of shared memory algorithms, and the preservation of mesh conformity at the interface between regions treated on different processes in distributed memory strategies.
In practice, physical simulation solvers are often parallelized over distributed memory architectures.Hence, their combination with a sequential or shared memory implementation of remeshing would imply to gather the whole mesh on one process, call the remesher and last, redistribute the mesh onto each parallel process before calling the solver.The intense exchanges of memory entailed in doing so are a well-known performance bottleneck for remeshing-based numerical simulation strategies, see [84] for a discussion.Moreover, the considered mesh is sometimes so large that it cannot be stored in the memory of one single node, thus ruling out the very possibility to use a sequential, or a shared memory implementation.These major drawbacks plead in favor of distributed memory parallel strategies for remeshing.
Conducting the remeshing of T over a distributed memory architecture first requires to partition T into disjoint submeshes T k , k = 1, . . ., K: no tetrahedra T ∈ T belongs to two different T k .The surface triangles at the interface between two of these submeshes constitute a surface mesh O, also referred to as parallel interface.The T k are sent and treated on different processes, or ranks; this raises the following issues: • The imbalance between the amounts of work conducted on each rank can hardly be predicted.Most often, the modifications involved in the remeshing process are non uniformly distributed over T (this is the case when, for instance, goal-oriented mesh adaptation is performed, see Section 2.3) and it is difficult to propose a priori an even repartition of the amount of operations across the T k .For instance, balancing the number of elements between ranks does not guarantee a balance of the work loads.• Enforcing the conformity of the reunion of the submeshes T k requires specific parallel data structures (e.g. a table maintaining the correspondence between the shared surface triangles S ∈ O and the supporting tetrahedra pertaining to different ranks), whose management is intricate.Notably, their update requires some parallel communication and synchronization between ranks when one tetrahedron T bearing an interface entity S ∈ O is modified, and a careful reconstruction of the parallel interface O is in order when the partition of T is modified (as is the case in the remeshingrepartitioning parallel strategies broached below).These concerns open the way to two strategies for conducting remeshing in parallel over distributed memory architectures, which essentially differ by their treatment of the entities S ∈ O at the interface between two of the submeshes T k .
• On the one hand, each remeshing operator could be applied in parallel to the entities in O, see e.g.[23] [39] [80] [30].Doing so requires a tight communication between the various processes sharing the considered configuration, in order to check the validity of the realized operation and to update consistently the mesh connectivity; • On the other hand, an iterative remeshing-repartitioning (or moving-interface) strategy could be used, which intertwines -A parallel remeshing of each domain T k on the associated process, while keeping the entities in O unmodified; -A new subdivision of the mesh T , creating new submeshes T k and interface triangulation O; see [52] [25] [42] [14] about this approach, and Fig. 6 for a 2d illustration.
It is expected that the repetition of this procedure allows each region of the mesh to be modified.In this direction, the most crucial issues lie in the repartitioning operations.Indeed, the remeshing of each region T k simply makes use of a sequential algorithm such as those presented in Section 2.2, but a fluid migration of the domain interfaces from one rank to the other is mandatory to ensure that no element in the mesh stay at the interface between submeshes, lest that it would stay unmodified.Moreover, this interface migration procedure must ensure a fair balance of the work load between ranks in terms of CPU cost, while keeping the amount of migrating data low, insofar as possible.

h-adaptive RANS and hybrid RANS/LES simulations of a nozzle with the Discontinuous
Galerkin method using unstructured meshes and isotropic mesh adaptation In this section, we present a mesh adaptation strategy devoted to the solution of the compressible steady Reynolds-Averaged Navier-Stokes (RANS) and the unsteady Zonal Detached Eddy Simulation (ZDES) equations on hybrid prismatic/tetrahedral grids using Discontinuous Galerkin (DG) methods, in the context of realistic, industrial applications.The developed mesh adaptation algorithm is applied to RANS and hybrid RANS/LES simulations on the PPRIME nozzle configuration for a DG formulation featuring elements with polynomial degree p = 1.

Context and motivation
The work presented in this section is motivated by the search for a simple mesh adaptation algorithm to improve the computational efficiency of DG simulations for complex flow configurations relevant in an industrial context, where the necessity of keeping the simulation cost "reasonable" is of utmost importance.Steady RANS simulations are well established, and extensively used for industrial purposes.Nevertheless, they may fail to capture the turbulence and noise generation mechanisms which are often required in the design process.On the other hand, scale-resolving simulations (relying for instance on the DNS, LES, or hybrid RANS/LES methods [89]), are capable of capturing the unsteady features in transitional flows, gas turbine combustors and nozzles.In this work, we assess both the steady RANS [95] and the ZDES hybrid RANS/LES models [40], in order to reduce the computational cost of the DNS and wall-resolved LES turbulence modeling approaches, which require a very large (often intractable in practice) number of degrees of freedom in space and time for capturing the smaller structures developed in the boundary layer.
DG methods are particularly suited for turbulent flow simulations thanks to their good dispersion and dissipation properties.In addition, these methods provide a high order of accuracy on arbitrary unstructured meshes, are suitable for parallel computing thanks to their compact stencil and provide a natural framework for hp-adaptation -not only the size h of the elements of the mesh can be adapted but also their polynomial degree p [46,66,100].
The new generation compressible flow solver CODA is employed to compute the flow solution and the error estimators on an unstructured mesh.The CODA solver, designed for an efficient use on current and future parallel HPC systems, is developed in partnership by Airbus, ONERA and DLR [67] and targets academic and industrial aerodynamic problems.The object-oriented CODA framework permits the integration of two spatial discretizations: a second-order Finite Volume (FV) and a DG discretization with variable order, applied to the Euler, Navier-Stokes, RANS and hybrid RANS/LES equations.

Presentation of the mesh adaptation strategy
The general strategy is to adapt the computational mesh with respect to a relatively affordable steady RANS simulation, which already captures some important features of the targeted unsteady simulation.This mesh is then used as the computational support of an accurate, albeit more expensive hybrid RANS/LES simulation.The mesh adaptation process implemented in this work outlines as follows.A first steady RANS simulation is carried out on a very coarse initial mesh.An a posteriori error estimator aimed at controlling the accuracy of the solution in the elements of the mesh, is then computed; the latter is described in the next Section 3.2.2.It is used to define a new mesh size prescription guiding the refinement of the elements where the error is high, as we discuss in Section 3.2.4.Remeshing is then conducted on the basis of this datum.The previous flow solution is projected onto the obtained mesh, and the process is iterated: a new steady RANS simulation is performed, etc.
Five such RANS adaptation steps are carried out, and a hybrid RANS/LES simulation is eventually realized on the finest mesh adapted to the RANS equation.
Note that the present strategy, whereby the accurate and expensive hybrid RANS/LES analysis is performed only on the finest mesh, which is adapted with respect to the relatively cheap RANS simulation, is preferred over adapting the mesh from the beginning with respect to the hybrid RANS/LES simulation.Indeed, an initial very coarse mesh would prevent the turbulent structures of the flow from developing and yield numerical instabilities as well as a dramatic increase in computational time for the whole adaptation process.
Note that this idea of employing adapted meshes obtained from steady simulations at a low computational cost, as the starting point of an unsteady turbulent adaptation procedure is expected to be a robust and reasonable cheap means to achieve a highly accurate hybrid RANS/LES simulation.The evaluation of this methodology will be addressed in future research.

Description of the error estimator for the steady RANS simulation
Several indicators based on the discretization error have been developed in the DG simulation literature.These types of error indicators are convenient thanks to their efficiency, locality, simplicity and low computational cost [71,85,58,87].The error estimator employed in this work is made of two contributions, similarly to what has been proposed in [34], [13] and [12]: the first one is based on the measure of the energy associated with the highest-order polynomial modes, the Small Scale Energy Density (SSED) [78], while the second one relies on the inter-element jumps (JUMP) [17] of the momentum.The resulting error indicator is local, inexpensive and flexible, in the sense that an error indicator based on the highest order modes of the solution is more reliable for a high polynomial degree p, while a jump-based error estimator is accurate for every value of p.The two aforementioned contributions are normalized by their respective maximum and minimum values over the tetrahedra T ∈ T , so that the considered error estimator finally reads: The quantity T is naturally defined at the level of the tetrahedra T ∈ T .In practice, most remeshing strategies rely on error estimators attached to the vertices of the computational mesh.A consistent value x may be attached to a vertex x by using a volume-weighted average of the values T at the elements sharing x as vertex.This indicator is the key ingredient in the definition of the imposed size h * x , as we describe next in Section 3.2.4.

Description of the mesh adaptation loop
As we have mentioned, five mesh adaptation steps are performed.At the end of the RANS simulation of step n, the element-based and vertex-based error estimators n T and n x are computed, and the next imposed size h n+1 x is defined (see Section 3.2.4).The computational mesh T n available at step n is eventually remeshed, resulting in the new mesh T n+1 .
The meshes T n are made of 2 regions: while most of the computational domain is filled with tetrahedra, another region, associated to the boundary layer of the physical phenomenon under scrutiny is composed of structured or pseudo-structured (i.e.prismatic) elements.At each step n, only the tetrahedral part of T n is remeshed.Moreover the user may decide not to remesh some of the tetrahedra in the computational domain.Summarizing, the computational mesh T n is split into two parts: where T n free is the tetrahedral zone subject to remeshing at step n, and T n fixed contains the constituent prisms of the boundary layer mesh, as well as the fixed tetrahedra.

Size prescription for the steady RANS simulation
The size prescription is based on the idea, already present in [17] and [88], that the error T attached to each element T ∈ T converges to zero with the rate p + 1, when no geometrical or physical discontinuities are observed.
At each adaptation step n, the imposed size h n+1 x at the vertex x ∈ T is updated according to the rule: , where n+1, * is a maximum value for the error, which is imposed globally to all the elements of the mesh, with the aim to provide a fixed increase in the number of elements of the mesh at each adaptation step [12,17].Adopting this strategy, the regions showing a larger value of the error estimator than the imposed target n+1, * are refined by a factor depending on the ratio n+1, * / n x between this target and the actual value of the error, while regions characterized by an error lower than n+1, * are left unchanged.

Numerical example: simulation of the subsonic turbulent nozzle jet flow
We appraise the above adaptive remeshing strategy in the physical context described in [22], for which experimental and numerical data are available in the literature, see also [57] and [79] for further investigations aimed at predicting turbulence generation and jet noise.
The main goal of this study is to simulate the turbulent isothermal subsonic jet issued from a nozzle with exit diameter d N = 0.05m.The nozzle has an axial symmetry with respect to the first coordinate axis e 1where (e 1 , e 2 , e 3 ) stands for the canonical basis of R 3 and we denote by (x 1 , x 2 , x 3 ) the coordinates of a point x ∈ R 3 in this frame.The operating conditions are defined in terms of the total pressure ratio P t /P ∞ = 1.7 and total temperature ratio T t /T ∞ = 1.15,where the t and ∞ subscripts refer to the stagnation and freestream states, respectively.The jet is assumed to be isothermal (T jet /T ∞ = 1.0), the jet Mach number is M jet = U jet /c jet = 0.9, and the Reynolds number equals Re , where U jet is the mean jet exit streamwise velocity, d N = 0.05m is the exit diameter of the nozzle, c jet is the speed of sound, ρ jet is the density and µ jet is the dynamic viscosity.
As we have mentioned, five mesh adaptation steps based on steady RANS simulations are performed in the present study, and a final unsteady hybrid RANS/LES simulation is performed on the finest RANS adapted mesh.
The initial mesh is generated with the pre-processing software ANSA [3]: the geometry of the body and the far field boundaries are created and meshed with surface triangles; then, the prismatic boundary layer surrounding the surface of the body (the internal and external walls of the nozzle) is obtained thanks to a normal extrapolation of these surface triangles.The remaining volume of the computational domain is eventually filled with tetrahedra.
Not only the boundaries of the computational domain and the prismatic elements of the boundary layer mesh, but also the tetrahedral elements which are internal to the nozzle, are fixed throughout the computation, and are not subject to remeshing.Indeed, tetrahedra free to change size inside the nozzle but constrained by the fixed size of the surface could severely deteriorate the quality of the mesh.
The initial mesh, containing around 1.5 million degrees of freedom is represented on Fig. 7, as well as the final mesh, resulting from the five-step adaptation process, containing around 10 millions of degrees of freedom.Remarkably, the mesh adaptation algorithm is capable of detecting the flow regions of interest, leading to a concentration of the elements around the potential core, and in the shear layer of the jet.The streamwise velocity profiles for four stations of the jet issued from the nozzle and the mean streamwise velocity profiles along the centerline are compared to the experimental and LES results obtained in [22].In order to assess our DG h-adaptive results obtained with a RANS model, we use for reference simulation a RANS second-order FV computation on a hexahedral structured mesh counting 48 millions of elements.
The meshes resulting from the third (4.5 millions of degrees of freedom), the fourth (6.9 millions of degrees of freedom) and the fifth adaptation steps (10.1 millions of degrees of freedom) are very similar, indicating that mesh convergence for the RANS case has been attained.
The results reported on Fig. 8 reveal that the RANS computations, conducted with either the adaptive DG and structured FV methods, show a reasonable agreement with the LES results inside the potential core of the jet and in its vicinity, but that they tend to underestimate the centerline velocity of the jet for x 1 /d N > 6, leading to an earlier dissipation of the streamwise axial velocity.This underestimation is a typical behavior of RANS models applied to jet flows, see [2].Nevertheless, the results of our adaptive simulation show a similar overall behavior as those of the FV reference computation, achieving closer results to the LES and experimental computations, with around 20% the number of degrees of freedom of the structured FV simulation (10.1 vs 48 millions).Moreover, the use of an adaptive process circumvents the difficulties posed by a classical, manual structured meshing process, especially when complex geometries are considered.
In Fig. 9, we report radial velocity profile cuts at four locations in the mesh.At x 1 /d N = 1 the adaptive simulations match almost perfectly the LES and experimental results, and show a large increase in accuracy with respect to the fine FV structured simulation.At x 1 /d N = 5, our adaptive simulations are still very close to the LES and experimental results, but show a slight over-prediction of the spreading rate of the jet, while at x 1 /d N = 10 and x 1 /d N = 15 the radial velocity profiles of RANS simulations show significant discrepancies with respect to the experimental and LES reference results, but still perform slightly better than FV computations.
Eventually, we investigate the improvement of the jet flow prediction entailed by the use of a hybrid RANS/LES approach on RANS-adapted hybrid prismatic/tetrahedral meshes.The formulation of ZDES used in this work forces the whole interior part of the nozzle to act in RANS mode, while the DES equations are solved in the rest of the domain, see [96].
The benefits that the unsteady simulation can bring to the finest mesh solution are shown with red dotted lines in Figs. 8 and 9.The velocity field of the unsteady simulation is averaged over a period of 150 times the characteristic time of the flow (d N /U jet ).Despite being a relatively coarse mesh with respect to the standards of unsteady turbulent computations, the velocity profiles match more closely the experimental and LES references, thanks to the lower degree of modeling of turbulence that hybrid RANS/LES and classic LES models introduce with respect to RANS.In Fig. 10 we represent the steady Mach RANS solution on the initial mesh and on the finest mesh, and an instantaneous snapshot of the Mach field on the finest mesh 4.1.A few generalities about shape and topology optimization Shape and topology optimization is about finding the "best" shape with respect to mechanical, geometrical and manufacturing specifications.Over the last decades, this discipline has been raising an increasing enthusiasm within the academic and industrial communities, as the dramatic increase in the cost of raw materials calls for the optimization of the shape of mechanical parts since the early stages of design; we refer to e.g.[6,7,8,16] for further motivations and comprehensive introductions to this field.
A typical shape optimization problem is of the form: where • Ω is the optimized shape, which is sought within the fixed computational domain D; • J(Ω) is a given performance criterion; • C(Ω) stands for a (collection of) constraint functional.Usually, the functionals J(Ω) and C(Ω) depend on the shape Ω via its mechanical behavior, that is, mathematically, via a state u Ω , solution to a "physical" partial differential equation posed on Ω (for instance, the heat equation, the elasticity system, etc.).
The treatment of shape optimization problems of the form (4.1) by constrained optimization algorithms requires a notion of derivative with respect to the domain.Here, we rely on Hadamard's boundary variation method [9,62,76,94]; when combined with the so-called adjoint technique from optimal control theory, it allows to calculate a "shape gradient", that is, a descent direction, for a function F (Ω) of the domain.The latter arises as a vector field θ : R 3 → R 3 such that, for small t > 0: intuitively, a "small" displacement of the boundary ∂Ω in the direction pointed by the vector field θ results in a new domain (Id + tθ)(Ω) with a slightly lower value of F (Ω).

Body-fitted shape and topology optimization
Usually, large modifications of the shape are expected in the course of the resolution of a shape and topology optimization problem of the form (4.1), making the use of Lagrangian strategies such as those described in Section 2.4 unrealistic in this context.
Multiple alternatives have been thought of in order to circumvent this problem, and notably the resort to the level set method outlined in Section 2.5, see [8] for the seminal contribution in shape optimization.
Recently, this framework has been successfully combined with remeshing techniques, resulting in a bodyfitted shape optimization method, which features an exact mesh of the shape Ω (and of the complement D \ Ω), see [5,37,50,51].In a nutshell, two complementary representations of a shape Ω ⊂ D are available at each stage of the process: • Meshed representation: On the one hand, Ω (and thus D \ Ω) is explicitly discretized in the mesh T of D: T is a valid, conforming mesh of the total computational domain D, which can be divided into two submeshes T Ω and T D\Ω of Ω and D \ Ω, respectively.• Level set representation: On the other hand, Ω is described as the negative subdomain of a level set function φ : D → R defined at the vertices of T , see (2.1).Dedicated numerical algorithms are then used to switch from one representation to the other.When a meshed representation of Ω is available, a corresponding level set function φ : D → R is calculated as the signed distance function to Ω at the vertices of the mesh T , for instance by using the fast marching method, see [90].Conversely, when Ω is described by a level set function φ : D → R on a mesh T of D, a new mesh T of D where Ω is explicitly discretized can be obtained thanks to the methodology described in Section 2.5.
In the course of the numerical resolution of a shape optimization problem of the form (4.1), each representation of the shape Ω is used depending on the performed operation, see Fig. 11: • The finite element analyses needed to calculate the state u Ω , and thereby a descent direction θ for the problem (4.1) can be conducted from the meshed representation of Ω. • The motion of Ω to the next iterate (Id + tθ)(Ω) can be realized with the level set representation φ of Ω, by solving the evolution equation (2.2) on the mesh T .This approach has the following benefits: • Accurate mechanical calculations can be realized on the exact mesh T Ω of the shape Ω, via the finite element method for instance.This mesh can be readily exported, and finite element analyses can be carried out by using an external software application in a black-box fashion; hence, this strategy totally decouples the update of Ω from the mechanical analyses needed to evaluate its performance and the sensitivities of the optimization functionals.• Since the structural interface ∂Ω is explicitly discretized at each step of the iterative procedure, this body-fitted approach simplifies the evaluation of geometric and mechanical quantities of interest near ∂Ω (such as the perimeter, or the curvature of ∂Ω, or the normal stresses applied on ∂Ω).• Since the level set method is used to realize the motion of the shape, dramatic evolutions of shapes can be accounted for, including topological changes.

Numerical example: optimization of an elastic mechanical system
We illustrate the proposed approach with the optimization of the design of an elastic two-component mechanical system.This numerical example was treated by using PISCO, a Research and Development software platform devoted to topology optimization that is under active development at IRT SystemX and Safran Tech.The industrial-grade solver Code Aster1 is used to perform the finite element analyses needed to evaluate the physical criteria and the corresponding sensitivities.
The considered setting is depicted on Fig. 12.The computational domain D is composed of two regions: an L-shaped support corresponding to a non optimizable region, which supports the boundary conditions of the physical analysis, and a structural block, which is the optimized domain, properly speaking.These are filled with two different linearly elastic materials -the support being characterized by a higher Young's modulus than the block -and they are linked by rigid mechanical connections, represented on Fig. 12, (b).Both regions are discretized using linear tetrahedral elements while the rigid connections are accounted for by four discrete 1-order elements.Dirichlet and linear sliding contact boundary conditions are prescribed on two regions of the boundary ∂D located respectively on the support and the block boundaries.Eventually, two load cases are considered, which are applied on another region of ∂D see Fig. 12 (c,d).
The goal of the present experiment is to minimize the volume of the structural block (so as to design a lighter structure) under four constraints: one on the compliance of the total structure Ω and one on the Von Mises stress field for each of the two considered load cases.The Von Mises constraints are formulated using a p-norm aggregation technique which is quite classical in stress-constrained topology optimization, see e.g.[48].
We rely on the body-fitted shape and topology optimization strategy introduced in Section 4.2 to solve this problem, starting with the full structural block of Fig. 12 (a) as initial shape.The optimized design is depicted in Fig. 13, and the associated displacement and Von Mises stress fields for both considered load cases are reported in Fig. 14.At convergence, all the constraints are fulfilled and the weight of the assembly has been reduced by 52%.
Note that a full remeshing of the structural system is carried out at each step of the optimization process.As we have already emphasized, the employed body-fitted approach simplifies greatly the definition and the evaluation of the physical criteria and sensitivities, since any finite element solver can be used in a non intrusive fashion.On the other hand, the remeshing stages involved in this strategy have admittedly a huge impact over the total computational cost of the process.Thus, topology optimization would greatly benefit from the development of dedicated parallel remeshing routines to target large-scale structural optimization challenges.

Context and motivation
This section illustrates the contributions of local remeshing algorithms in geophysical modeling, and more particularly in the management and the reduction of uncertainty over the identification of geological interfaces.The underground properties of the earth govern many natural or anthropic phenomena occurring over multiple time scales, such as mineral element migration and concentration processes, earthquake initiation and propagation, landslides, groundwater flow, hydrocarbon recovery, CO 2 sequestration and geothermal heat recovery, etc.In particular, the representation of geological interfaces is a crucial aspect of the numerical simulation of these processes, as these surfaces delimit regions containing different geological materials (rock units) characterized by specific ranges of mineralogical and physical properties.For example, in sedimentary basins, these rock units arise as layers separated by horizons and can be described at multiple nested scales, as exemplified on Fig. 15.The geometry of these layers may contain challenging geometric features such as sharp creases or low-angle contacts between surfaces, for instance when geological layers vanish laterally because some previously deposited material has been eroded (see Fig. 15 (c)).Under the effect of tectonic forces, cracks naturally present in the rock may grow as fractures, and eventually become faults under the accumulation of tangential displacement, which results in possibly complex juxtapositions of materials on either side of the fault, see e.g.Fig. 15 (d).
From the numerical point of view, the simulation of subsurface phenomena such as those mentioned above relies on a model describing the features of the underground medium (and notably the geophysical interfaces) at the scale of interest, which is then meshed and populated with physical rock properties before the physical equations of interest are solved.
The aforementioned complex underground properties of the earth are therefore pivotal in these computations.Unfortunately, they are generally inaccessible to direct observation.Rather, they must be inferred from indirect data, such as surface measurements, borehole sample analyses or in situ geophysical measurements and geophysical images.This is a highly challenging task, since such data are sparse, ambiguous and insufficient to precisely characterize the domain.Hence, geological knowledge is a key ingredient of the interpretation process and significant uncertainty exists about the existence, location, connectivity and geometry of interfaces between different rock units [101].
Numerical approaches based on the level set method (see Section 2.5) have raised a significant interest within the geoscience modeling community, as a convenient tool to cope with these uncertainties.One category of methods, referred to as "implicit structural modeling", consists in creating level set functions from typical interpretation points.In tetrahedral formulations [54,24], this paves the way to the generation of meshes where the discontinuities under scrutiny are explicitly discretized, as described in Section 2.5, which allows to conveniently manage faults and other geological interfaces.
Other investigations have used perturbations of level set functions for generating alternative geological models representing uncertainties [33,103].In the perspective of reducing the related uncertainties, several deterministic or stochastic inverse methods have been formulated to take into account indirect geophysical data [83,63,56], but most of them feature one single level set function during the inversion process.Going further calls for inverting several and possibly finite open interfaces while preserving their geological consistency [101,60].Moreover, even though some Monte Carlo methods are available to address the uncertainty in the geometric configuration [18,28], their dissemination has been hampered by the difficulty to robustly automate mesh generation for arbitrary configurations.
Using the level set method in combination with body-fitted tetrahedral meshes in the spirit of Section 2.5 is a very attractive strategy to handle these requirements: it indeed combines geometric adaptivity, accuracy and versatility to several PDE discretization schemes such as the finite element method or the control volume finite element method.5.2.Numerical example: sensitivity of a two-phase Darcy flow with respect to the existence and the geometry of fractures In this section, we consider the porous fractured network configuration depicted on Fig. 16: we assume 6 fracture intersections along two distinct observation lines have been located.As often in geoscience practice, the performed observations do not contain the same amount of information, neither are they equally reliable: some of these data provide the fracture position and orientation while others provide at best an approximate insight into its position.From these incomplete observations, spatially exhaustive scenarii are generated, representing possible fracture configurations compatible with the observations.In the present 2D case, this task is easily achieved by hand, but many sampling methods exist to generate this type of model, based on statistical and/or physical reasoning, see for instance [19,38,61,29,102].
In each situation, the initial porous region is assumed to be filled with oil.We then simulate the injection of water in the underground, from the lower left corner of the computational domain D; both fluids are recovered in the upper right corner of D. Assuming the fluids to be incompressible and immiscible, the result of this experiment is predicted by the resolution of the equation for the pressure P : which is complemented with the saturation equation for each phase Here, ϕ stands for the rock porosity; k is the rock permeability, S α , µ α , k rα and f α are respectively the saturation, viscosity, relative permeability and fractional velocity of the phase α, and q α is the volumetric source term.The numerical simulation solver is based on a control-volume finite element approach, where fractures are represented as lower-dimensional mesh elements, see [65,73].In this experiment, we assume for simplicity that the relative water permeability is given by the square of the water saturation.The intact rock has a porosity of 0.2.Fractures are assumed to be open (porosity of 1).The matrix and fracture permeabilities are taken equal to 10 −15 m 2 and 10 −9 m 2 , respectively.The viscosity for oil and water are 0.45 and 1 mPa.s,respectively.Each fracture is parametrized by its center, size and azimuth.Starting from an initial mesh of the computational domain D, the various fractures at stake are inserted one after another by relying on the implicit surface meshing technology described in Section 2.5: briefly, the signed distance field φ to the fracture is computed by a least-square optimization method [54].A subregion of the mesh of D is then defined from the fracture center and size to bound the editing.Finally, the methodology of Section 2.5 is used to insert the zero level set of φ into the portion of the mesh corresponding to this region.
This strategy presents numerous assets, when it comes to exploring multiple possible fracture geometries and to appraising their impact on the flow.In particular, this iterative approach allows to locally perturb an existing scenario (via a perturbation on the corresponding level set function) in order to explore uncertainties over its geometry (such as the fracture shape or position), while ensuring that the mesh features adequate element size and quality.
On the contrary, the opposite strategy consisting in meshing the fracture network at once is highly challenging, since its very thin features inevitably lead to meshes containing a large number of elements, or ill-shaped mesh elements, which either slows down simulations or affect their stability.Admittedly, several simplification approaches have been proposed in this context, with the aim to approximate the target fracture geometry while preserving the mesh quality, while retaining an affordable amount of elements, see [77,53,10,64].
This simplification problem is significantly easier in our framework where each fracture is sequentially inserted as the 0 isosurface of an adequate level set function.In order to prevent the creation of "ill-shaped" elements in the remeshed model when the new fracture nodes are located on another interface, function values that are "too close to 0" on these nodes are set to 0.
In order to illustrate the proposed workflow, we consider the example depicted on Fig. 16: the initial scenario depicted on Fig. 16 (a) includes 6 and 4 fractures belonging to two fracture sets of azimuth ca.45 • and 105 • , respectively.In the configuration of Fig. 16 (b), only one fracture has been extended with respect to that in (a), leaving the rest of the model unchanged.This minor perturbation has a limited effect on the production curves, and induces only a moderate change in the saturation field.The scenario of Fig. 16 (c) involves a more drastic change with respect to (a): one fracture is replaced by two fractures.This type of operation could occur for instance in a reversible jump Monte Carlo Markov Chain transition [18].As this change affects the fracture connectivity, the impact on the production curves and on the saturation field is significant: the water breakthrough at the production well occurs later, meaning that this fracture configuration yields a better sweeping efficiency of the injected water.The scenario of Fig. 16 (d) represents what could occur after multiple transitions from the initial model (a): although the orientations and sizes are comparable, most fractures have been moved, removed or replaced.Interestingly, the production curves in the scenario (d) are similar to those of (a) and (b).This illustrates the ill-posedness of this particular type of problems.This motivates the use of stochastic inverse methods rather than deterministic optimization to address subsurface uncertainty.
6. Dynamic and parallel mesh adaptation for premixed flame fronts and liquid/gas interfaces 6.1.Dynamic mesh adaptation for material interfaces

Motivation
The numerical simulation of dynamic material surfaces in flows such as turbulent flame fronts and liquid/gas interfaces has long been a thorny issue.One reason is that the reaction or interface forces occur in very thin and moving layers around the surface, so that an accurate capture of the latter requires a very fine mesh resolution in its vicinity.On the other hand, such interfaces typically undergo dramatic displacements, spanning the whole computational domain, which makes it impossible to adapt the mesh once for the whole simulation.Dynamic mesh adaptation is a more realistic alternative to this unfeasible, brute-force approach.
It consists in modifying the mesh regularly so that it features a fine mesh size in the vicinity of the interface and a coarser one away from it.Such a strategy is usually hindered by two factors: (1) The huge increase in CPU cost entailed by the mesh adaptation process, which cancels out the potential gain allowed by applying the numerical solver with a reduced number of elements.(2) Even though the cost of mesh adaptation is affordable, another challenge comes from the compromise between the thickness of the refined layer around the interface (and thereby, the size of the mesh) and the frequency at which remeshing is conducted.This section presents solutions to reduce the overhead of dynamic mesh adaptation by relying on a parallel implementation, and to find a compromise between the thickness of the refined layer and the remeshing frequency.From the algorithmic viewpoint, we use the YALES2 flow solver [74] for the numerical simulation of the physical problems at stake.This solver also implements a version of the iterative remeshing-repartitioning algorithm relying on mmg to perform parallel mesh adaptation.As was presented in Section 2.6, the efficiency of parallel remeshing strategies carried out on distributed architectures (in terms of the quality of the resulting mesh, notably) crucially depends on the ability of the procedure to modify the whole mesh T .More precisely, in each submesh T k (processed on rank k), not only the regions lying "far" from the interface O between ranks, but also those lying close to O have to be remeshed.While the modifications of the former regions can be done in parallel by independent calls to the mesh adaptation library, the treatment of the interface between ranks is much more delicate.As we have outlined in Section 2.6, several strategies are available to tackle this issue, and that retained in the present work is the so-called moving interface method.In a first step, the mesh is refined inside each rank, with respect to a metric whose definition is detailed in Section 6.1.3,while leaving the regions near O untouched.During this step, the target metric is interpolated from the old to the new mesh.This adaptation step is therefore local and does not require any communication between the ranks.In a second step, low-quality elements are sent from one rank to another.The scheduling of this movement is driven by a constrained load balancing algorithm.This operation has to ensure that the grid remains properly balanced from one rank to another, i.e. that all ranks have approximately the same number of tetrahedra, and it has to enforce the movement of the interface.To achieve this, one solution is to impose strong weights at the connections (or arcs) between the elements at the interface.Thus, when minimizing the edge cut of the connectivity graph of the mesh T , the load balancing algorithm will be strongly enticed to place the new interface away from the previous interface so as to avoid cutting heavy-weight connections.Once the elements have been moved from one rank to another, the process can be repeated until convergence.This convergence is achieved when the difference between the target and resulting metrics is below a threshold value and when the element skewness does not exceed a limit value.The performances of this method depend on the efficiency of each of the aforementioned operations, namely:

A two-level parallel mesh adaptation method
(1) The sequential mesh adaptation library; (2) The interpolation of the data stored on the grid.
(3) The load balancing algorithm; (4) The cell and data transfer from one rank to another; (5) The parallel connectivity reconstruction.
All these steps may limit the performance of the overall process when a large number of cores is involved.In order to alleviate all the above potential limitations (except those concerning the sequential mesh adaptation algorithm, which we assume to be sufficiently efficient in the following), we rely on a two-level domain decomposition method.The latter is based on a slight refinement of the general parallel remeshing principles exposed in Section 2.6, which is illustrated in Fig. 17: in a nutshell, most of the above tasks are actually conducted at a coarser scale than that of the tetrahedra T ∈ T , that of groups (hereafter referred to as cell groups) composed of several thousands of tetrahedra.For instance, the load balancing algorithm is applied to the connectivity graph of the cell groups instead of that of the cells, which allows for a significant gain in CPU time as the latter graph is much lighter and easier to partition.Likewise, the cell migration operation is performed cell group by cell group and it is combined to automatic packing/unpacking with non-blocking send/receive calls to speed-up the transfer.Finally, the interpolation of numerical quantities from the initial mesh to the modified one is done at each sub-step, benefiting from the 2-level domain decomposition to locate the nearest source tetrahedron from a given destination element.The only drawback of this two-scale strategy is that it requires additional merging and splitting algorithms in order to provide a large block to the sequential remeshing software instead of small cell groups.However, this drawback does not alter the performances of the complete workflow, and the overall gain far compensates this overhead.

Choice of the metric
The dynamic mesh adaptation of a moving material interface such as a flame or a liquid/gas interface requires frequent remeshing steps.When using low-dissipation and low-dispersion numerical schemes such as those used in Large-Eddy Simulation and Direct Numerical Simulation of turbulent flows, the remeshing steps have to introduce as small numerical perturbations as possible.The numerical perturbations entailed by the adaptation of a mesh T into a new mesh T better suited to the further evolution of the interface mainly come from the interpolation errors of numerical quantities defined on T onto T , whose values near the interface are of critical importance.These perturbations can be made negligible if the remeshing occurs away from the interface, i.e. if the mesh metric is kept constant near the interface under scrutiny and if the tetrahedra nearby are frozen during this operation.Keeping the size prescription constant near the interface and controlling its variation away from it leads to define the metric as in Fig. 18: it is set to ∆x min in a narrow band with thickness 2N p ∆x min around the interface, and its gradient is bounded outside this region, so that the metric attains that of the initial coarse mesh "far" from the interface.With this definition and as illustrated on the left-hand side of Fig. 18, when the interface moves of a distance less than N p ∆x min , the size prescription changes where either the old or the new metric shows linear variations, while keeping a constant value at the interface location; this ensures an exact interpolation of the data within this protected zone.The remeshing effort is therefore concentrated in the region where the metric changes most, i.e. in the linear metric gradient zone.6.2.Numerical example: application to the lean-premixed PRECCINSTA burner Figure 19.Dynamic mesh adaptation in the PRECCINSTA burner test case of Section 6.2: the instantaneous CO mass fraction and the heat release fields are displayed, together with a slide of the 3D tetrahedral mesh.
The computational benefits of the proposed dynamic mesh adaptation strategy are first illustrated with the example of a semi-industrial lean-premixed burner, operating with a methane-air premixing under atmospheric conditions.This burner, called PRECCINSTA, has been widely used for model development and validations [15,75].In this burner, the turbulent flame is anchored to the injector thanks to a swirl motion, which creates large inner and outer recirculation zones.A Large-Eddy Simulation with the same models  and numerics as in [15] is performed with the help of the dynamic mesh adaptation strategy presented in Section 6.1.The results are depicted in Fig. 19; in there, a mesh resolution of 300 microns in the flame front and of 600 microns in the primary combustion zone are used.These coarse resolutions are chosen to better highlight the metric definition.As the flame is highly unsteady due to flame/turbulence interactions, the flame front displacement is tremendous and it eventually features topological changes, as isolated burning pockets are created.The remeshing frequency is therefore piloted by a Courant-Friedrichs-Lewy condition based on the displacement velocity, which ensures that the flame front remains in the refined zone.The flame topology is rendered as a progress variable iso-surface colored by the temperature in Fig. 20.In there, the left-and right-hand side subfigures correspond to two treatments of this simulation using static meshes (i.e. the computational mesh is adapted at the beginning of the process, and is then left untouched) containing 110 millions and 878 millions tetrahedra, respectively.In the center, the same situation is treated with our dynamic mesh strategy; the computational mesh is obtained with the background mesh of the 110 million cell mesh while the resolution at the flame front is the same as that of the 878 million cell mesh.The flame topology of the front obtained thanks to the dynamic strategy features small-scale wrinkles similar to the refined static case (878M) while the total CPU cost is around 4 times smaller.This observation is confirmed by the plot of the flame curvature distribution in Fig. 21.The distributions of the dynamic and of the refined static cases are the same while the simulation based on the coarser static mesh exhibits smaller values of the curvature.
6.3.Numerical example: application to the droplet impact on a cone  In this section, the benefits of dynamic mesh adaptation are illustrated with the simulation of the impact of a water drop on a superhydrophobic cone characterized by its semi-angle 50°-a situation which has been studied experimentally in [47,26].It is an interesting case as the drop undergoes various topological changes after its impact on the cone: it subsequently takes the shape of an attached ring before bouncing and breaking into a number of smaller droplets.The modeling of the liquid/gas interface and that of the contact angle of 163°are detailed in [86].The liquid/gas interface is captured with a conservative level set algorithm described in [41] and the contact angle is imposed through a curvature modification at the contact line.The metric definition follows the strategy exposed in Section 6.1.3,with a thickness of the refined zone tailored so that 30 cells are in the initial drop radius.The results and the mesh are presented in Figs.22 and 23.The experimental results are very well reproduced in terms of topology and of break-up dynamics.Such simulation would be intractable without dynamic mesh adaptation.

Figure 1 .
Figure 1.(a) Fine geometric approximation of a domain Ω 1 , with a mesh containing bad quality elements; (b) good quality mesh of Ω 1 ; (c) good quality mesh of a domain Ω 2 , accounting for a poor geometric approximation of Ω 2 ; (d) mesh of Ω 2 with a good geometric approximation property.
w s w 4 p + d s G 0 9 i e t d n 6 5 r 6 p 1 F q V s 9 Z p k 3 x p X c 5 e 8 H z Q e u 4 X D k v n z R O S 7 W D 7 K r z 2 M M + j u g + L 1 D D D e p o m u w X v O L N u r Z 8 K 7 H S H 6 m V y z y 7 + D W s 5 2 + I Q Y 9 3 < / l a t e x i t > p < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 Y D 6 o f 0 g y 5 w b H r N r P Y n t 3 Z y t b + s f V m l Y M 2 e Z N s W n 2 e X s B c 8 H j e N y 5 a x 8 c n N a q h 5 k V 5 3 H H v Z x R P d 5 j i q u U E P d d v m M F 7 w 6 1 4 5 y x s 7 k W + r k M s 8 u f g z n 6 Q s S / J I Q < / l a t e x i t > • < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 Y D 6 o f 0 g y 5 w b H

Figure 3 .
Figure 3. (Upper row) Adaptation of a mesh with respect to a geometric size prescription;(left) initial mesh with size map, calculated on the basis of the geometric approximation; (middle) adapted mesh with respect to this size prescription; (right) adapted mesh with respect to a finer geometric approximation requirement.

Figure 5 .
Figure 5. (Left) One "level set" function φ is defined at the vertices of the mesh T of the square-shaped computational domain D; its 0 level set (which is not explicitly discretized in T ) is depicted in red; (middle) explicit discretization of the 0 level set of φ into the mesh T ; the intermediate, low-quality mesh T temp is obtained; (right) high-quality mesh T obtained after remeshing T temp .

Figure 6 .
Figure 6.Illustration of an iterative remeshing-repartitioning parallel mesh adaptation strategy performed over 3 processors: at each iteration, the mesh is divided into 3 submeshes T k , k = 1, 2, 3, which are distributed over as many processors; each T k is then remeshed, independently of the other two, while the entities at the rank interfaces (cyan edges) are preserved.

Figure 7 .
Figure 7. Illustration of the PPRIME nozzle configuration studied in Section 3.3; (left) zoom of the initial mesh; (right) computational mesh after 5 steady RANS adaptation steps in the jet zone.

Figure 8 .
Figure 8. Velocity profiles (left) on the centerline and (right) on the lipline for the PPRIME nozzle example of Section 3.3.

Figure 11 .
Figure 11.Level set advection and body-fitted remeshing for interface tracking.(Left) initial level set defined at the vertices of a large conformal mesh and associated velocity field; (middle) level set update using the Hamilton-Jacobi advection equation; (right) domain remeshed to fit the zero isovalue of the advected level set function.

Figure 12 .
Figure 12.(a) The two-component mechanical system test case considered in Section 4.3; the non optimizable L-shaped support region is represented in orange and the optimized block is in blue; (b) rigid connections linking both regions; (c) clamped region of the device; (d) the surface supporting contact boundary conditions is represented in green, and that submitted to loads is in light orange.

Figure 13 .Figure 14 .
Figure 13.(Left) optimized shape of the structural block in the example of Section 4.3; (right) cut in the mesh of the computational domain D, divided into material and void parts, in red and blue, respectively.

Figure 15 .
Figure 15.Example of a challenging geometrical configuration in a 3D geological model, as discussed in Section 5. (a) Global view of a faulted stratigraphic model; dark lines account for depth contours and colored lines in the upper horizon represent stratigraphic unconformities; (b) cut-off view; (c) the low angle contacts lines at the stratigraphic unconformities caused by erosion are displayed in red; (d) view of how horizons (in transparency) are offset by faults, forming small geometric features and low-angle intersections on cut-off lines on either sides of each fault.

Figure 16 .
Figure 16.Several fracture network scenarii in the numerical example of Section 5.2, together with the associated multiphase porous medium behaviors.The spatial fracture observations and the associated uncertainties are displayed as grey rectangles and ellipses along thin dotted lines.The fractures are bold dark lines, which are explicitly discretized in the computational mesh.The color scale shows the water saturation.Each graph shows the water fraction and produced volume at the upper right corner, after injection of one pore volume (PVI) in the lower left corner.The dashed ellipses highlight the changes between (a) and each alternative scenario.

Figure 17 .
Figure 17.Schematic of the two-level moving interface parallel adaptation strategy developed in Section 6.The thin black lines represent the interfaces between cell groups and the different colors correspond to as many ranks.

Figure 18 .
Figure 18.Illustration of the construction of the size map involved in the remeshing strategy of Section 6.1 (left) and its displacement (right).

Figure 20 .
Figure 20.Comparison of the numerical simulations using static and dynamic meshes at several resolutions in the PRECCINSTA burner example of Section 6.2.

Figure 21 .
Figure 21.Comparison of flame front curvature distribution for static and dynamic meshes in the PRECCINSTA burner.

Figure 22 .
Figure 22.Impact of a water droplet on a superhydrophobic cone as considered in Section 6.3.

Figure 23 .
Figure 23.Dynamic mesh adaptation in the example of the impact of a water droplet on a superhydrophobic cone in Section 6.3.