G Jannoun

E Hachem

J Veysset

J-F Zaragoci

T Coupez

FULLY SPACE-TIME METRIC BASED ANISOTROPIC MESH ADAPTATION FOR UNSTEADY PROBLEMS

Keywords: Anisotropic Mesh Adaptation, Adaptive Time-Stepping, Paradoxical meshing technique, Time-dependent problems

. The latter is developed based on the length distribution tensor approach and the associated a posteriori edge based error analysis. Then we extend the mesh adaptation technique to contain adaptive time advancing. A newly developed time error estimator is constructed and intends to homogenize the global error over space and time. The main purpose of this work is the development of a novel meshing algorithm, the paradoxical meshing, that provides optimal space and time meshes suitable for several simulation time subintervals. The advantage of the proposed method relies in its conceptual and computational simplicity as it only requires from the user a number of nodes and a frequency of adaptation according to which the mesh and the time-steps are automatically adapted. Numerical solutions on time-dependent problems demonstrate the accuracy and efficiency of the proposed space-time error estimator.

INTRODUCTION

Despite the increasing computer performances and the progress of computational fluid dynamics in modelling and simulating time dependent PDEs, numerical restrictions are still present and caused by the complexity of the numerical simulations.

Anisotropic mesh adaptation has proved to be a powerful strategy to improve the accuracy and efficiency of finite element methods. It enables the capture of multi-scaled physical or mechanical phenomena. The method, as developed in [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF], allows the creation of highly stretched and highly directional elements leading to very good capture of the gradients of the solution and the internal and boundary layers. Moreover, it provides a good level of accuracy within a reasonable degree of freedom. Another extension was proposed in [START_REF] Coupez | Adaptive Time-step with Anisotropic Meshing for Incompressible Flows[END_REF], and accounts for time-step adaptation. Based on the derived error estimator in space and the solutions at the previous times, the proposed algorithm automatically computes an appropriate time-step for the following computations.

The above described mesh adaptation is optimal for steady problems. It is theoretically valid for transient CFD applications only when applied, together with the time adaptation method, at every solver iteration. Nevertheless this would excessively increase the computational cost and lead to an accumulation of interpolation errors. Fixing a low frequency of adaptation can help solving these issues but this might lead to a time lag between the mesh and the solution. To overcome these problems, a new fully adaptive method is proposed in this paper: the paradoxical meshing. It intends to predict the solution's evolution over a period of time and to automatically generate corresponding mesh and set of time-step sizes.

This paper is structured as follows: we start section 2 with a brief description of the classical anisotropic mesh adaptation. Section 3 is dedicated to the time adaptive technique. The good performance of these methods is evaluated on a 3D example with complex geometry. The extension of these two methods into the paradoxical meshing algorithm is described in Section 4. Finally, in section 5, we test the efficiency and accuracy of the space-time adaptive algorithm on time-dependent problems.

Construction of an anisotropic mesh for stationary problems

In [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF], we have developed an a posteriori edge based spatial error estimator relying on the length distribution tensor approach. Working on a nodal based metric, an anisotropic mesh adaptation procedure is obtained under the constraint of a fixed number of nodes.

Edge based error estimation

We consider u ∈ C 2 (Ω) = V and V h a simple P 1 finite element approximation space:

V h = w h ∈ C 0 (Ω), w h | K ∈ P 1 (K), K ∈ K
where Ω = K∈K K and K is a simplex (segment, triangle, tetrahedron, ... ).

We define

X = X i ∈ R d , i = 1, • • • , N
as the set of nodes of the mesh and we denote by U i the nodal value of u at X i and we let Π h be the Lagrange interpolation operator from V to V h such that:

Π h u(X i ) = u(X i ) = U i , ∀i = 1, • • • , N
As shown in figure 1, we define the set of nodes connected to node i by

Γ(i) = j , ∃ i K ∈ K , X i , X j are nodes of K
By introducing the following notation: X ij = X j -X i and using the analysis carried in [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF], we can set the following results:

∇u h • X ij = U ij , (1)
Figure 1: length X ij of the edge joining nodes i and j (left). Varying the edge in its own direction (right).

|| ∇u h • X ij U ij -∇u(X i ) • X ij || ≤ max Y ∈[X i ,X j ] |H(u)(Y )X ij • X ij | , (2)
where H(u) = ∇ (2) u is the associated Hessian of u. A definition of the projected second derivative of u is obtained using (1) and the interpolation operator on ∇u :

∇g h X ij • X ij = g ij • X ij (3)
where

∇g h = Π h ∇u, g i = ∇u(X i ) and g ij = g j -g i .
Using a mean value argument, we set that:

∃y ∈ [x i , x j ]|g ij • X ij = H(u)(Y )X ij • X ij .
This projection is considered as an expression of the error along the edge:

e ij = g ij • X ij (4)
However a gradient recovery procedure is needed as the gradient of u is not known and is not necessarily continuous at the nodes of the mesh.

Gradient Recovery

Based on an optimization analysis, the author in [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF][START_REF] Coupez | Adaptive Time-step with Anisotropic Meshing for Incompressible Flows[END_REF] proposes a recovery gradient operator defined by:

G i = (X i ) -1 j∈Γ(i) U ij X ij (5)
where

X i = d |Γ(i)| j∈Γ(i)
X ij ⊗ X ij is the length distribution tensor at node X i . Note that this construction preserves the second order:

G i -g i • X ij ∼ H(u)X ij • X ij
where G i is the recovery gradient at node i (given by ( 5)) and g i being the exact value of the gradient at node i. The error is evaluated by substituting g by G in (4):

e ij = G ij • X ij

Metric construction from the edge distribution tensor

Taking into account this error analysis, we construct the metric for the unit mesh as follows:

M i =   d |Γ(i)| j∈Γ(i) X ij ⊗ X ij   -1
For a complete justification of this result, the reader is referred to [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF][START_REF] Coupez | Adaptive Time-step with Anisotropic Meshing for Incompressible Flows[END_REF].

Error behavior due to varying the edge length

We examine now how the error behaves when we change the length of the edges by stretching coefficients

S = s ij ∈ R + , i = 1, • • • , N , j = 1, • • • , N , Γ(i) ∩ Γ(j) = φ
Stretching factors s ∈ R are employed to link the error variations to the changes in edge lengths:

X ij = sX ij || e ij || = s 2 ||e ij || = s 2 ||G ij • X ij || (6)
where e ij and X ij are the target error at edge ij and its associated edge length.

The metric associated with S can be redefined as:

M i = |Γ(i)| d X i -1 with X i = d |Γ(i)| j∈Γ(i) s 2 ij X ij ⊗ X ij (7)
being is the length distribution tensor. Let n ij be the number of created nodes in relation with the stretching factor s ij and along the edge ij. When scaling the edges by a factor s ij , the error changes quadratically so that the number of created nodes along the edge ij is given by:

n ij = e ij e ij -1 2 = s -1 ij
Recall that e ij denotes the induced error for edge X ij . As per node i, the number of created nodes along the different edges' directions is given by the following tensor:

N i = X i -1   d |Γ(i)| j∈Γ(i) n ij 2 X ij ⊗ X ij  
So that the total number of created nodes per node i is:

n i = det   (X i ) -1   d |Γ(i)| j∈Γ(i) n ij 2 X ij ⊗ X ij    
Assuming a uniform totally balanced error along the edge e ij = e = cst, we get a direct relation between N and e as follows:

n ij (e) = s -1 ij (e) = e e ij -1 2
Hence for a node i we have

n i (e) = det   (X i ) -1   d |Γ(i)| j∈Γ(i) n ij (e) 2 X ij ⊗ X ij    
Replacing n ij (e) by its expression,

n i (e) = e -d 2 det   (X i ) -1   d |Γ(i)| j∈Γ(i) e ij X ij ⊗ X ij    
and this is equivalent to:

n i (e) = e -d 2 n i (1)
so that the total number of nodes in the adapted mesh is:

N = e -d 2 i n i (1).
Hence, the global induced error for N nodes can be determined by:

e(N ) =   N i n i (1)   -2 d
Therefore the corresponding stretching factors under the constraint of a fixed number of nodes N are given by:

s ij = e e(N ) -1 2 =   i n i (1) N   2 d e -1/2 ij
Note that the mesh does not change during time advancing but at a certain time level t n . Hence, an optimal mesh at a time level t n need not be an optimal one at t n+1 which is the case when propagating a discontinuity. This raises the question about the frequency of remeshing.

As time-dependent problems exhibit arbitrary progression with time, the duration of applicability of a mesh cannot be known apriori. When the time-step size is greater than the length of the mesh's time interval, the solution may propagate into a non preadapted region of the domain resulting in a mesh/solution lag. Adapting the mesh at every solver iteration guarantees that the spatial error remains bounded. Nevertheless, this approach increases significantly the computational cost and leads to the accumulation of interpolation errors polluting the solution.

The focus now can be waived to the choice of the time-step size computed at each solver iteration. The main objective is to produce a time-step that preserves the accuracy of the mesh adapted solution and avoids unnecessarily small time-step sizes. This method was first introduced and validated in [START_REF] Coupez | Adaptive Time-step with Anisotropic Meshing for Incompressible Flows[END_REF]. In this work, we revisit the developed method with the intention of extending it to a fully space-time adaptivity.

We apply the above described analysis in 1D where the only variable is time. Denote by T = {t 0 , • • • , t n-1 , t n , t n+1 , • • • } and let t nk = |t n -t k | , n, k ∈ T be the temporal nodes and ∆t n = t nn+1 the time increments as shown in figure 2. Assume that the solution is already computed on the whole domain up to time t n . The aim is to choose an appropriate time-step ∆t n .

Without loss of generality, the analysis will be carried on an arbitrary spatial node i. Note that at a spatial node i, we only have one time edge to be determined (t n t n+1 ). Define {τ nn+1 } to be the temporal edge scaling (stretching) factor such that:

ẽn+1,n = τ 2 n+1,n e n+1,n tn+1,n = τ n+1,n |t n+1,n | (8)
where e n+1,n is an approximation of the interpolation error from t n to t n+1 , ẽ and t are the target error at the temporal edge t n t n+1 and its associated edge length. Let u i n-1 , u i n and u i n+1 be the solutions at node i and times n -1, n, and n + 1, respectively. Using a forward difference approximation, we have that u n+1 -u n = un ∆t n and u n-1 -u n = -un-1 ∆t n-1 . Then applying the recovery gradient in 1D, we get:

un = u i n,n+1 ∆t n + u i n,n-1 ∆t n-1 ∆t 2 n + ∆t 2 n-1
and the quadratic interpolation error:

e i n,n-1 = ui n,n-1 ∆t n-1 and e i n,n-1 = τ i n-1 2 ui n,n-1 ∆t n-1 where ui n,n-1 = ui n -ui n-1
. Now using the equidistribution error argument, we write

ẽi n,n-1 = e n (N )
with e n (N ) being the maximal error in space for a total number of nodes N . Hence the stretching factor of the time-step size is given by:

τ i n-1 = e(N, t n ) e i n,n-1 1 2
and the optimal time-step is determined by:

∆t n = min i τ i n-1 ∆t n (9)
Looking closely at this formula, we notice that it requires the solution at time t n+1 which is not yet computed. Therefore instead of computing the optimal time-step ∆t n we calculate ∆t n-1 and we let ∆t n = ∆t n-1 .

Application to 3D heat transfer and turbulent flow inside an industrial furnace

In this section, we will apply the classical space and time-adaptive methods to simulate the heat transfer and fluid flows inside an industrial furnace with complex geometry. The objective of this test case is to show the applicability and the potential of the developed algorithm in simulating long time heating inside large scale furnaces.

The furnace is modelled as a hexagonal section duct of 2.7 × 8.1 × 5.3m 3 forming one heat transfer zone. All computations have been conducted by starting with a gas at rest with a constant temperature of 1463 • C. Adiabatic temperature is considered at all other boundaries for sake of simplicity. The heated air is pumped into the furnace at a velocity 14.3m/s by a circular burner with 6m diameter and located at the left vertical wall. The air is vented out of the furnace through two outlets positioned at the bottom of the left vertical wall. The 3D computations aim at simulating an hour of heating and have been conducted in parallel on 16, 2.4Ghz Opteron cores.

Figure 3 (top) shows the isothermal distribution at different time-steps. When the hot fluid spreads along the volume of the furnace, it induces a turbulent motion within the geometry. This forced convection is caused by the interaction of the moving stream and the G. JANNOUN, E. HACHEM, J. VEYSSET, J-F. ZARAGOCI, AND T. COUPEZ 3) clearly reflects the expected flow pattern. A number of small vortices inside different buffer zones can be observed. The latter are due to the turbulence dissipation and mixing of the hot and cold air. The numerical results obtained with our space and time adaptive algorithm reflect well the efficiency and potential of the methods. Figure 3 (bottom) highlights how well the mesh is adapted to sharply capture, with highly stretched elements, the gradients of the solution, the boundary layers and the emerging vortices. The algorithm builds up the mesh in a way to maximize the accuracy of the numerical solution with a fixed number of nodes ( 100, 000). Note that the mesh is adapted according to the velocity components and its norm as in this test case it is the motor of the induced airflow and spread of the temperature. The results describing only one hour of the heating process required 100 days of computations with a fixed time-step equal to 0.005s. Significant CPU time and computational cost were saved by applying our time adaptive procedure as it required only 2 days of calculations (see table 1). Figure 4 shows the evolution of the time-step sizes allowing at the same time a certain level of accuracy and an acceleration of the computations.

G. JANNOUN, E. HACHEM, J. VEYSSET, J-F. ZARAGOCI, AND T. COUPEZ 7(left), in order to reduce the temporal error, preventing the progress of the solution in time so that it remains in phase with the adapted mesh. Recall that the method aims at giving a better efficiency than the classical algorithm. This is exactly what we notice when comparing figures 7(left) and 7(right); the paradoxical meshing technique produces time-step sizes that are almost 200 times larger than those generated by the classical algorithm. Hence, the resulting computations will be about 200 times faster reflecting the high efficiency of the novel method. Note that the inner loop of the algorithm is repeated only two times to get the optimal meshes and time-step size for which the solution remains bounded.

CONCLUSIONS

In this paper, we have presented a classical anisotropic mesh adaptation that showed very good performance when applied together with the new time adaptive technique for resolving time dependent problems. An extension of these algorithms lead to a novel and very powerful method for full adaptation known as the paradoxical meshing. This method demonstrated its efficiency in generating meshes and time-step sizes that guarantee the convergence of the solution all over computations for a limited number of nodes and a fixed frequency of adaptation.

Figure 2 :

 2 Figure 2: Temporal discretization at the spatial node i.

Figure 3 :

 3 Figure 3: Isotherms (Top) and corresponding adapted meshes (bottom) at three different time-steps.

Table 1 :

 1 CPU time for computing the solution with non-adaptive and a space/time adaptive methods stationary fluid inside the furnace. The numerically obtained temperature distribution (fig.

Acknowledgement:

The authors gratefully acknowledge the support from the "ANR: Agence nationale de la recherche".

paradoxical meshing: Full adaptivity algorithm

When dealing with steady state problems, the classical mesh adaptation technique, presented in section 2, performs pretty well for converging the mesh-solution couple. Nevertheless, this method is no more optimal when applied to unsteady problems as the physical solution evolves in time. Together with the time-adaptive method introduced in section 3, the classical mesh adaptation can be efficiently adapted to time dependent problems.

In this paper we are interested in developing a space and time fully adaptive algorithm. The latter aims at anticipating the solution progress over a period of time and generating the optimal mesh that is adequately adapted, for a fixed number of nodes, to the evolving solution along that time interval. The analysis is carried out on a (3D+1D) mesh, i.e. the computations are performed synchronously on a 3D spatial mesh and a 1D temporal mesh. We aim at generating a mesh that holds for several solver iterations together with the corresponding optimal set of time-step sizes. Note that the user can assign a frequency of adaptation and the algorithm will accordingly adapt the meshes.

The principle consists of dividing the simulation time [0, T ] into n SI subintervals:

that will in turn be divided into n freq -1 subintervals where n freq is the frequency of adaptation assigned by the user. We call the adaptation method a paradoxical meshing as the resulting mesh is being adapted to n freq time-steps while adapting every n freq steps. The mesh and the set of time-step sizes are computed through an iterative process along which we try to converge both meshes (the spatial and the temporal one) to the optimal configurations that give the most accurate solution for the corresponding interval of time. At every iteration, we consider each of the n SI intervals at a time and divide it into subintervals

The solution is predicted at each of the {t k l } temporal nodes using the numerical scheme. From these solutions we construct a vector field

and we estimate the edge based spatial errors:

We compute a global error e(N, T k ), as in the classical approach, for equidistributing the error on the edges of the discrete domain. An optimal metric is deduced controlling the spatial error over [T k , T k+1 ]. We also compute the temporal error at the nodes t k l as in section 3. Once this is done we optimize the time-step sizes ∆t k,l by equidistributing the error in space and time:

An optimal size δ t k of the time interval [T k , T k+1 ] is recomputed as follows:

These δt k are given to the 1D mesher that will generate a new discretization of the interval [0, T ] as well as its corresponding temporal nodes T k . Notice that the number of subintervals n SI will be automatically changed due to this remeshing. The above described algorithm is repeated iteratively until convergence of the metrics and the set of time-step sizes. Therefore, for each interval [T k , T k+1 ] a metric is computed accounting for the solution's transient evolution. At convergence, a mesh is generated from this metric field. Computations are then resumed on the predicted optimal set of meshes with the corresponding set of time-step sizes.

The novel method that we presented herein is perceived not only as a fully adaptive technique but also as a space and time accurate way of solving time-dependent problems within reasonable computational costs.

Numerical example: Rotating Circle

In this section, we assess the performance of the newly developed fully-adaptive algorithm on a numerical test case and compare the result with the classical space and time adaptive techniques. We consider a 2D circle of radius 0.3 located at (0.5, 0.5) in the computational domain [-1.5, 1.5] × [-1.5, 1.5]. The simulation consists of rotating the circle in the counterclockwise direction at the rate θ = 1rad/s. The objective is to test the capability of the anisotropic paradoxical meshing technique to accurately capture the dynamically evolving interface. Figure 6 shows the adapted mesh, made up of 10, 000 nodes, for the time interval [t n , t n+1 ] made up of 10 sub-intervals of time. We can clearly see how refined the mesh is at the location of high gradients of the solution and how accurate is the capture of the interface as it rotates from t n (left) to time t n+1 (right). The elements all along the interface are isotropic yielding a well respected curvature. The time intervals' lengths

] are automatically generated to guarantee the validity of the mesh for 10 consecutive time-steps. Figure 7 presents the time-step sizes for the first few iterations of the algorithm. The periodic variation of the time-steps is in good agreement with the nature of the problem, as the circle rotates at a constant rate and maintains the same behavior all over computations.

Using the classical mesh adaptation with the time adaptive technique and adapting the mesh every 10 iterations, the generated time-step sizes will be too small, as seen in figure


J P Moitinho De Almeida

P Díez

C Tiago

Eds N Parés

H Digonnet

L Silva

T Coupez

VI International Conference on Adaptive Modeling and Simulation ADMOS 2013

Keywords: Massively parallel computation, Anisotropic mesh adaptation, Multigrid

In this paper, we present developments done to obtain efficient parallel computations on supercomputers up to 8192 cores. While most massively parallel computation are shown using regular grid it is less common to see massively parallel computation using anisotropic adapted unstructured meshes. We will present here two mains components done to reach very large scale calculation up to 10 billions unknowns using a muligrid method over unstructured mesh running on 8192 cores. We firstly focus on the strategy used to generate computational meshes and in particular anisotropic ones adapted to capture a quite complicated test function. Then we will briefly describe a parallel multigrid method. Performance test over a large range of cores from 512 to 8192 cores is then presented using the French national supercomputers Jade and Curie. The last section will present a calculation done on smallest number of cores on our own cluster, but using more realistic data obtain directly from experimentation. The goal is to be able to realize such kind of simulation on really complex micro structure obtain by tomography at a larger scale.

INTRODUCTION

From the last years computers power increase mainly by cores multiplication (rather than clock rate rise) inside each CPU (classically 8 or 16 cores in 2012) and also for the supercomputers of the top500 that contain several thousand to more than one million cores. this context impose to develop fully parallel softwares to at least be able to take advantage of new computer hardware. If we look at numerical simulation, computation times are always to important and has soon as they decrease the user want to have more precise computation by introducing more physical properties or by looking in computing lower scale calculation to improve the material behavior. One of the mains way to reduce computation time is to use parallel computers and in the ideal case divide the CPU time by the number of cores used. For that it is then necessary to develop algorithms that could run on massively parallel computer containing hundred, thousand or even more cores like ones present in the top500 list of supercomputers. We present here some algorithms done to allow our application CimLib to run over massively parallel supercomputer using up to 8192 cores. The first section will focus on parallel mesh adaptation with anisotropic elements taht lead to generate smaller meshes for a given precision. A second part will briefly present parallel multigrid method implemented to reduce the complexity of the algorithm used to solve linear systems and fully utilize the power and memory given on a computer. We present parallel performance on solving incompressible Stokes equation using from 8 to 8192 cores and enable resolution of a linear system containing more that 10 billions unknowns by using 8192 cores. The last section will present more reasonable simulation but on a more realistically and complex case. It consists in computing a flow through a micro structure given directly from a real one using tomography.

parallel anisotropic mesh adaptation

parallel mesh adaptation strategy

We briefly present here the methodology used to parallelize our unstructured and non hierarchic mesher, MTC [START_REF] Coupez | Génération de maillage et adaptation de maillage par optimisation locale[END_REF]. At the beginning of this work, the sequential mesher already exist and is likely to evolve (we recently add anisotropic mesh size [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF]), so we want to keep it as it was, by not including parallel instructions inside. So, we have not parallelized directly it, but only give the oprtunity to use it inside a parallel context. The strategy is then to combine local remeshing (inside each processor domain, freezing the interfaces between the partitions) [START_REF] Coupez | Parallel meshing and remeshing[END_REF] and parallel mesh repartitioning [START_REF] Basermann | Dynamic load-balancing of finite element applications with the DRAMA library[END_REF] to move unremeshed interfaces inside domains for the next step, in order to be able to apply the remeshing procedure on theses zones during the next phase. Figure 1 illustrate this strategy on a simple geometry using seven processors.

optimization and performance

This strategy of parallelization leads to several iterations (depending of the space dimension) between the mesher and the repartitioner, but the work to be done at each iteration decreases quickly. For example, in a 2d space : the first remeshing step is proportional to a surface, the second to a line and the last one to a point. For the repartitionner, as we only need to move bad quality zones inside the domain in order to remesh it, the proportion of elements and nodes to be migrate across the processors decreases in the same way. For that reason the cost of using a multi steps strategy must be really close to a single remeshing step: the global work done by the mesher is approximately the same, and only some few mesh migrations are done during repartitioning steps. For that purpose an optimization of the update procedures in case of small changes has been done (the time spent to remesh or migrate 10 elements is then close to zero). This optimization is based on a permute-copy-past algorithm that leads to reduce the complexity from N (the data size) to m (the moving size) with m ≪ N. Permutation removes any copy of all the data by making instead some few permutations. After that, work is performed only in the cutting zone and pasted back after being done. This optimization was essential to maintain the iterative strategy costless, and so provide a good parallel efficiency to the mesher.

anisotropic error estimator

Anisotropic meshes are used to keep the same accuracy of isotropic meshes but with a smaller number of nodes and elements by allowing different mesh size depending of the direction. The use of such anisotropic meshes could reduce drastically the size of the mesh and in particular in 3d where a uniform refinement of a factor 2 increase the number of nodes by a factor 8 that has to be compared to only a factor 2 if the refinement is only needed in one direction. For anisotropic mesh adaptation we need use a mesh size definied by a scalar for isotropic mesh and a symmetric positive defined tensor (named metric) for anisotropic one. In order to enable mesh adaptation we need to provide to the mesher a continuous field of mesh size one simple solution is to give a P1 interpolate field where a local mesh size is defined at each mesh nodes. The way to compute this field could be a average of the elements metrics containing this node with the difficulty of choosing a good average value for tensors. A more direct way to defined this metric consist in using the distribution tensor describe in [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF] and given by :

X i = 1 |Γ(i)| j∈Γ(i) X ij ⊗ X ij
where X ij is the vector ---→ X i X j between two neighbor nodes i and j and Γ(i) the neighbors nodes of i.

This tensor give the anisotropic size of the envelope containing of the elements belonging to node i, and the natural metric for the mesh is then simply defined by M i = X -1 i to insure that using this metric this envelope has an unit length.

The anisotropic error estimator is well presented in [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF] and we will here only present the mains ideas. To be able to build a anisotropic error estimator we first need to be able to define an anisotropic error and for that we will use the same technique of the distribution tensor. For a given Level Set we will compute the error over nodes edges as :

e ij = G ij X ij
where G ij is the gradient variation of the level set function across the edge ij. With this anisotropic error and given a number of nodes mesh we could define (for demonstration refer to paper [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF]) a metric tensor to generate an adapted mesh that will uniformly distribute the error across its edges. The metric is given by :

M i = 1 |Γ(i)|   j∈Γ(i) s 2 ij X ij ⊗ X ij   -1 with p ∈ [1, d] λ =           i j∈Γ(i) e p p+2 ij   A         p+2 2
where d is the space dimension, A is the global number of edges in the mesh and p is the power of the power law estimation for the number of edges generated in the mesh when we want to divide one edge. For example: for a one direction anisotropic mesh divide an edges in that direction by a factor 2 will only multiply by 2 the global number of edges in the mesh and then p = 1; in case of an isotropic mesh divide an edge by two will generate eight times more edges so p = 3.

application test case

We present in this section an application test case adapted to massively parallel computation. It consist in computing a adapted mesh to allow a good representation of a complicated function defined by :

f (x) = a • a(x -x 0 ) + a • a(x -x 1 )
where

a(x) = tanh Esin 4 * N + 1 2 Π x and x 0 = (0, .., 0) x = (1, ..., 1)
E and N are two parameters that can be adjusted. They respectively influence the thickness and the complexity of the details present in the function: if E increase the size of the detail decrease and if N increase then each detail contain become more and more complex. The figure 2 represent the test function on an uniform square mesh containing around 50000 nodes for E,N parameters equal to 1, 2 and 4. We notice that for E=4 and N=4 the uniform mesh start to be to small to obtain a good representation of the function f.

As soon as these parameters increase the function become really difficult to capture and we need to activate both the anisotropic mesh adaptation and parallel computation to generate meshes fine enough to represent the test function. The bottom of figure 3represent the 23 millions nodes anisotropic mesh generate to represent the test function with E=16 and n=6 on the unit square. It clearly illustrate the benefit of being able to execute anisotropic mesh adaptation in parallel, here using 256 cores.

The figure 4 represent the same test case in a 3d context using smaller parameters N = 1 and E = 4. E and N parameters. The 3d adapted mesh contain 2 millions nodes and have been generate using 50 mesh adaptation iterations using 64 cores of Jade supercomputer in 15.5 hours .

Multigrid solver

Using powerful computers containing a large numbers of cores could lead using very big meshes with several millions nodes. Solving physical equation using for example the finite element method will ends in solving very large linear systems. Using iterative methods like Krylov ones will parallel preconditioner give good parallel performance but the algorithm complexity (around O(n 3 /2 ) in 3d) become a bottleneck as soon as the number of unknowns start to be important (taking a problem 2 times bigger and using 2 times more cores the time spent to solve will be greater than 2). To over come this difficulty we have implemented a parallel multigrid solver using the framework give by PETSc [START_REF] Balay | PETSc Users Manual[END_REF]. To do this we need to provide for each level the system to be solve and also the interpolation/restriction matrix operator. Of course all these matrices must be build and store in a distributed way to allow using a large number of cores. More details are presented in [START_REF] Digonnet | Multigrid using Adaptive Unstructured Meshes for Massively Parallel Computation[END_REF] but a particular attention have been given to reduce the communication during the construction of the interpolation operator using a parallel octree localization algorithm and pixel mask filter distributed over processors to limit false positive external nodes detection. The figure 5 present parallel performance of the multigrid resolution for Stokes equation over a 2d square mesh of 800 millions nodes and using 512 to 4096 cores. Time spent to solve the system decrease from 96.7 seconds with 512 cores to 17.7 seconds with 4096 cores, in the same time the number of multigrid iteration only increase from 11 to 13 and leads to a speed-up of 5.46.

Even if a good parallel performance of the multigrid solver is important, the main interest is the scalability of this approach that make us able to solve very very large systems. An ideal scalability will be give by a algorithm that will use exactly the same computing resources for a given per core problem size independently of the number of cores uses. A scalability test have been done from 8 to 8192 cores on the multigrid solver. The results are very interesting : the memory used per core during the process stay nearly constant (from 1.94 to 2.14 GB) as well as the time spent to assemble all levels systems (from 9.7 to 11 seconds) and the multigrid iterations (from 13 to 11); only the resolution time increase a beat more from 90 to 148 seconds. This augmentation of 50% of the time spent to solve the problem between 8 and 8192 cores is still reasonable and may also be reduce using more multigrid level over 8192 cores as the coarse level system size may become to big (around 600 000 unknowns) to be solve in a efficient way. To conclude this multigrid solver implementation have been able to solve a global system (based on the Stokes equation) containing more than 10 billions unknowns using 8 levels in 158 seconds using 8192 cores and 17.5 PB of memory.

From real to virtual

In this section we present some work done to enable making simulation using real micro structure. Today tomography could produce some really nice (well defined but also heavy) image of real material as shown on the left of figure 6 [START_REF] Orgéas | In-plane conduction of polymer composite plates reinforced with architectured networks of copper fibres[END_REF]. This image is constituted of a large number of gray level voxel (equivalent to pixel in 2d), the back ones represent the polymer matrix and white one solid fibers. From this image we extract the an isosurface mesh (shown on the right of figure 6). Depending of the image definition the number of voxel and thus the number of faces in the isosurface mesh could become really important and may need to be executed in parallel using some small enough part of the whole domain using some similar technique as octree. Once we have the isosurface mesh, it could be immerse into the computational mesh by computing for all the mesh nodes the distance to the isosurface that give a level set representation of the micro structure [START_REF] Laure | Permeability computation on a REV with an immersed finite element method[END_REF]. To improve the 6: From real to virtual: a large 3d voxel image of the microstructure obtain by tomography and treated with a segmentation algorithm, its surface mesh given by the isosurface 127, and finally the computational domain with an anisotropic adapted mesh.

immersion of the microstrure we could lauch some anisotropic mesh adaptation procedures that will provide the computational adapted mesh as shown in bottom of figure 6.

If the previous figures 6were done to illustrate the mesh adaptation to the real microstructure using anisotropic mesh, using bigger image for the microstructure generate very big mesh and justify using massively parallel computer for both anisotropic mesh adaptation but also for computing the solution. Figure 7 represent the computed flow using Stokes equation across the immerged microstucture. Computation was done using 96 cores on a 10 millions nodes adapted mesh to a microsuture image containing 900x900x220 voxels. This is this type of computation than we plan to do in a close future but over much more larger image with around 4000x4000x4000 voxels given by new tomographic acquisition. For being able to do that all the computational chain need to be run in parallel from the image generation to the visualization of the computed results.

CONCLUSIONS

In this paper, we have presented the parallelization strategy for the mesher that consist in iterate between a remeshing step with frozen interfaces and repartitionning step used to migrate interfaces. In that way the mesher engine stay sequential and do not contain any parallel instruction so it could still be developed (introduce anisotropic mesh) without any specific knowledge. The implementation of a permutation cut and past optimization gives excellent parallel performance as well as a better quality control for meshes generated. Some examples of anisotropic distributed meshes adapted to well represent a complicated test function are give and clearly show the interest of using anisotropic mesh rather than isotropic ones : a 23 millions nodes anisotropic mesh adapted using 128 cores can capture the test function with parameters N=6 and E=16 leading to a very complex function, a equivalent representation will need a picture containing around one thousand billions pixels.

A massively parallel multigrid method working on unstructured meshes have been also presented. The construction of interpolation/restriction operators between very large distributed meshes and using a large number of cores need a powerful localization algorithm used to determine which element contain a node. For that we use a parallel octree localization method improve by a domain filter define by a pixel mask to reduce false detection of non local nodes. By using this filter we have a important reduction of communication even when domains are unrelated???. Very good parallel performance have been presented for the multigrid approach use to solve Stokes equation on very fine meshes over a large range of cores from 8 up to 8192 cores. The multigrid method developed allow to better use all the resources of massively parallel computers (CPU and memory). Over 8192 cores a global system containing more than 10 billions nodes have been solved with an 8 levels multigrid solver in 148 seconds, using a total of 17,5Po of memory. It is important to notice that for such a resolution we have been oblige to remove the 32 bits integer limit in our code.

Finally we presented some more realistic simulation done on real micro structure obtained by tomography. Tomography of a real material gives a 3d image (voxels) segmented to represent the micro structure. This image is then used to build a Level-Set function. A anisotropic adapted computational mesh is build in parallel to fit the micro structure. Once this adaptation is done we could compute the flow through the micro structure. We have presented here one simulation done over a 10 millions nodes 3d adapted mesh using a micro structure image with 900x900x200 voxels. Visualization of the results has been done using a parallel version of the Paraview [START_REF] Henderson | The ParaView Guide: A Parallel Visualization Application[END_REF]software, coupled with distributed results files and using python scripts [START_REF] Digonnet | Making Massively Parallel Computations Available for End Users[END_REF].

Future work we be done to improve and validate our algorithms and be able to scale over tens to one hundred of thousands cores as the 32 bits integer limits have already been removed. We will also look to simulate non linear materials on a larger VER given by new tomographic acquisition leading to 4000x4000x4000 voxels images.

Figure 1 :

 1 Figure 1: Illustration of the parallelization strategy used for the mesher in a 2d case with a refinement by 4, over 7 processors. Each image represent one step of the iterative methodology consisting in successively call parallel remeshing with frozen interfaces and parallel repartitioning, until convergence to good quality mesh. The last repartitionning phase is done to load balance the work by considering finite element resolution to come.

Figure 2 :

 2 Figure 2: Influence of the two parameters N and E on the complexity of the test function. This test function is plot on an uniform 2d mesh containing 50 000 nodes. From the top to the bottom the parameter N = 1, 2, 4 and from the left to the right the parameter E = 1, 2, 4. On the left bottom the uniform mesh used to plot the test function.

Figure 3 :

 3 Figure 3: Anisotropic mesh adaptation around the test function f (x) = g • g( x -0 + g • g( x -1 ) with g(x) = tanh(Esin( 4 * N +1 2 πx)) with 2 sets of parameters : on the top N=3, E=16 and at the bottom N = 6, E = 16. For these two sets we present : the function on the unit square, one cross detail and a zoom on the mesh for a deep detail inside the cross detail. The 2d adapted meshes containing respectively 1 million and 23 millions nodes. They where build in respectively 576 seconds over 128 cores and 3 hours on 256 cores of Jade supercomputer. The smallest mesh size is around receptively 10 -5 and 10 -6 .

Figure 4 :

 4 Figure 4: Anisotropic mesh adaptation in 3d around the test function with N = 1 and E = 8 for a final mesh containing 2 millions nodes.

Figure 5 :

 5 Figure 5: Parallel efficiency for the multigrid resolution using 512 to 4096 cores. The number of multigrid iteration needed to converge stay almost constant (efficiency close to one) has well as the parallel efficiency per iteration. At the global point of view the parallel efficiency is a bit worse due to the combination of per iteration parallel performance and increasing number of iteration but still very good (close to one) has the problem size per core decrease as number of cores used increase.

Figure

  Figure6: From real to virtual: a large 3d voxel image of the microstructure obtain by tomography and treated with a segmentation algorithm, its surface mesh given by the isosurface 127, and finally the computational domain with an anisotropic adapted mesh.

Figure 7 :

 7 Figure 7: Large numerical simulation of the 3d flow across a microstructure done using a big 900x900x220 tomographic image. Flow is compute over a 10 millions nodes adapted mesh using 96 cores. The isosurface represent fibers microstructure, vector the velocity field.

Acknowledgments

This work was granted access to the HPC resources of [CCRT/TGCC/CINES/IDRIS] under the allocation 2012-[x2012066109] made by GENCI (Grand Equipement National de Calcul Intensif)