CNAM Emmanuelle Itam

Stephen Wornom

Bruno Koobus

Biot France Lemma

A Volume-agglomeration multirate time advancing for high Reynolds number flow simulation

come     A Volume-agglomeration multirate time advancing for high Reynolds number flow simulation

SummaryA frequent configuration in computational fluid mechanics combines an explicit time advancing scheme for accuracy purposes and a computational grid with a very small portion of much smaller elements than in the remaining mesh. Two examples of such situations are the travel of a discontinuity followed by a moving mesh, and the large eddy simulation of high Reynolds number flows around bluff bodies where together very thin boundary layers and vortices of much more important size need to be captured. For such configurations, multi-stage explicit time advancing schemes with global time stepping are very accurate, but also very CPU consuming. In order to reduce this problem, the multirate time stepping approach represents an interesting improvement. The objective of such schemes, which allow to use different time steps in the computational domain, is to avoid penalizing the computational cost of the time advancement of unsteady solutions which would become large due to the use of small global time steps imposed by the smallest elements such as those constituting the boundary layers. In the present work, a new multirate scheme based on control volume agglomeration is proposed for the solution of the compressible Navier-Stokes equations equipped with turbulence models. The method relies on a prediction step where large time steps are performed with an evaluation of the 0 fluxes on macro-cells for the smaller elements for stability purpose, and on a correction step in which small time steps are employed only for the smaller elements. The accuracy and efficiency of the proposed method are evaluated on several benchmarks flows: the problem of a moving contact discontinuity (inviscid flow), the computation with a hybrid turbulence model of flows around bluff bodies like a flow around a space probe model at Reynolds number 10 6 , a circular cylinder at Reynolds number 8.4 × 10 6 , and two tandem cylinders at Reynolds number 1.66 × 10 5 and 1.4 × 10 5 . computational fluid dynamics, multirate time advancing, explicit scheme, volume agglomeration, unstructured grid, hybrid turbulence model

INTRODUCTION

A frequent configuration in Computational Fluid Dynamics (CFD) calculations combines an explicit time advancing scheme for accuracy purpose and a computational grid with a very small portion of much smaller elements than in the remaining mesh. Two typical examples are the following:

A first example is the hybrid RANS/LES simulation of high Reynolds number flows around bluff bodies. In that case, very thin boundary layers need be addressed with extremely small cells. When applying explicit time advancing, the computation is penalized by the very small time-step to be applied (CFL number of order 1). But this is not the only interesting region of the computational domain. An important part of the meshing effort is devoted to large regions of medium cell size in which the motion of vortices need be accurately captured. For these vortices, the efficient and accurate time-step is of the order of the ratio of local mesh size by vortex velocity. We can apply an implicit scheme with such a time-step, which would produce a local CFL of order 1 for the vortices advection and a local CFL of order hundreds for the boundary layer. However, this may have several disadvantages. First this standpoint completely neglects a possible need of unsteady accuracy in the small-cell region, due to unsteady separation, for example. Second, considering the need of accuracy for vortices motion on medium cells, highly accurate explicit schemes are easily assembled. This includes the TVD third-order ones, and the standard fourth-order Runge-Kutta method (RK4). In contrast, high-order implicit schemes are complex and cpu consuming. More simple implicit schemes using backward differencing show much more dissipation than explicit schemes.

Our second example concerns an important complexity issue in unsteady mesh adaptation. Indeed, unsteady mesh adaptive calculations are penalized by the very small time-step imposed by accuracy requirements on regions involving small space-time scales. This small time step is for example an important computational penalty for mesh adaptive methods of AMR type [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF]. This is also the case for unsteady fixed-point mesh-adaptive methods as in [START_REF] Belme | Time accurate anisotropic goaloriented mesh adaptation for unsteady flows[END_REF]. In that latter method, the loss of efficiency is even more crucial when the anisotropic mesh is locally strongly stretched. In [START_REF] Belme | Time accurate anisotropic goaloriented mesh adaptation for unsteady flows[END_REF], due to this loss, the numerical convergence order for discontinuities is limited to 8/5 instead of second-order convergence. This limitation also applies to mesh adaptation by mesh motion. In the present work, our first numerical example concerns the computation of an isolated traveling discontinuity. The discontinuity needs to be followed by the mesh, preferably in a mesh-adaptive mode. Except if the adaptation works in a purely Lagrangian mode, an implicit scheme will smear the discontinuity of the solution. An explicit scheme will applied a costly very small time step on the whole computational domain.

In order to overcome these problems, the multirate time stepping approach represents an interesting alternative. A part of the computational domain is advanced in time with the small time-step imposed by accuracy and stability constraints. Another part is advanced with the larger time-step giving a good compromise between accuracy and efficiency.

The development of multirate schemes was first limited to ODEs with different time scales and their application was restricted to a small number of special industrial problems. The first multirate method was due to the pioneering work of Rice [START_REF] Rice | Split Runge-Kutta method for simultaneous equations[END_REF]. In this work, a system of first-order ODEs made of a latent component (slow variation) and an active component (fast variation) was considered. Runge-Kutta type integration methods were developed, in which different integration steps (large for the latent component and small for the active component) were used and appropriate extrapolations were made for the latent components when the active components are integrated. Following this early work, many other developments have been conducted in the field of multirate schemes and ODEs [START_REF] Rice | Split Runge-Kutta method for simultaneous equations[END_REF][START_REF] Andrus | Numerical solution of systems of ordinary differential equations separated into subsystems[END_REF][START_REF] Gear | Multirate linear multistep methods[END_REF][START_REF] Skelboe | Stability properties of backward differentiation multirate formulas[END_REF][START_REF] Sand | Stability of backward Euler multirate methods and convergence of waveform relaxation[END_REF][START_REF] Andrus | Stability of a multi-rate method for numerical integration of ODE's[END_REF][START_REF] Günther | Multirate row methods and latency of electric circuits[END_REF][START_REF] Engstler | Mur8: a multirate extension of the eightorder Dormer-Prince method[END_REF][START_REF] Engstler | Multirate extrapolation methods for differential equations with different time scales[END_REF][START_REF] Günther | Multirate partitioned Runge-Kutta methods[END_REF][START_REF] Günther | Multirate partitioned Runge-Kutta methods[END_REF][START_REF] Savcenco | A multirate time stepping strategy for stiff ordinary differential equations[END_REF]. These multirate strategies usually combine classic numerical integration methods (Backward Differentiation Formulas, Adams methods, Runge-Kutta schemes, Rosenbrock-Wanner methods), in which a large integration step is used for the slow subsystem of EDOs and a smaller integration step is used for the fast subsystem, with extrapolation/interpolation algorithms for coupling active to latent parts and vice versa. In order to assess their stability and efficiency, these multirate methods were applied to (more or less large) systems of stiff ODEs, often derived from electric circuits problems. As an example, in [START_REF] Günther | Multirate row methods and latency of electric circuits[END_REF], a multirate 4-steps Rosenbrock-Wanner method was implemented to solve systems of 250 -4000 ODEs which model electric circuits, leading to a gain in efficiency up to 2.8 compared to a classic 4-steps Runge-Kutta method. Fewer and more recent multirate works were conducted in the field of PDEs and hyperbolic conservation laws [START_REF] Constantinescu | Multirate timestepping methods for hyperbolic conservation laws[END_REF][START_REF] Sandu | Multirate explicit Adams methods for time integration of conservation laws[END_REF][START_REF] Mugg | Construction and analysis of multi-rate partitioned Runge-Kutta methods[END_REF][START_REF] Kirby | On the convergence of high resolution methods with multiple time scales for hyperbolic laws[END_REF][START_REF] Seny | An efficient parallel implementation of explicit multirate Runge-Kutta schemes for discontinuous Galerkin computations[END_REF][START_REF] Löhner | The use of domain splitting with an explicit hyperbolic solver[END_REF], and rare applications were performed in CFD [START_REF] Löhner | The use of domain splitting with an explicit hyperbolic solver[END_REF][START_REF] Seny | An efficient parallel implementation of explicit multirate Runge-Kutta schemes for discontinuous Galerkin computations[END_REF]. The first multirate work dealing with hyperbolic problems and applications in CFD was made by Löhner et al. [START_REF] Löhner | The use of domain splitting with an explicit hyperbolic solver[END_REF]. A domain-splitting method was developped in which the Euler equations are discretized by a second order explicit finite element scheme on a domain splitted in several subregions of different grid resolution. For the time-integration of the equations, larger time steps are used in subregions of coarser mesh fineness leading to an increase in efficiency. One and two-dimensional test-cases were considered : a shock tube problem, and supersonic inviscid flows around a circular cylinder and past a wedge. The results showed that the method can handle shocks and that a speedup of 2 between the multirate scheme and its single-rate counterpart can be reached. Another work of interest, more recent, is the one of Constantinescu et al. [START_REF] Constantinescu | Multirate timestepping methods for hyperbolic conservation laws[END_REF] which focuses on the development of multirate methods for the solution of one-dimensinal scalar hyperbolic equations. They propose strategies that are based on an appropriate transition between subregions of different local stability conditions where a classic Runge-Kutta scheme is used, by introducing buffer regions, where an adapted Runge-Kutta method is employed. For the case of the one-dimensional Burger equation, speedups up to 2.5 were obtained.

In the context of parallel computing, an interesting work on multirate methods is the recent one of Seny et al. [START_REF] Seny | An efficient parallel implementation of explicit multirate Runge-Kutta schemes for discontinuous Galerkin computations[END_REF]. This work focuses on the efficient parallel implementation of explicit multirate Runge-Kutta schemes in the framework of discontinuous Galerkin methods. The multirate Runge-Kutta scheme used is the approach proposed by Constantinescu [START_REF] Constantinescu | Multirate timestepping methods for hyperbolic conservation laws[END_REF]. In order to optimize the parallel efficiency of the multirate scheme, Seny et al. propose a solution based on multi-constraint mesh partitioning [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF]. The objective is to ensure that the workload, for each stage of the multirate algorithm, is almost equally shared by each computer core i.e. the same number of elements are active on each core, while minimizing inter-processor communications. The Metis software is used for the mesh decomposition, and the parallel programming is performed with the Message Passing Interface (MPI). The efficiency of the parallel multirate strategy is evaluated on three test cases: the wind driven circulation in a square basin and the propagation of a tsunami wave using a shallow water model (two-dimensional), and the acoustic propagation in a turbofan engine intake using the linearized Euler equations (three-dimensional). It is shown that the multi-constraint partitioning strategy increases the efficiency of the parallel multirate scheme compared to the classic single-constraint partitioning. Speed-up as high as 4.4 were obtained.

In the present work, we propose a new multirate scheme based on control volume agglomeration which is at the same time very simple to develop in an existing software relying on explicit time-advancing, and well suited to a large class of finite volume approximations. The agglomeration produces macro-cells by grouping together several neighboring cells of the initial mesh. The method relies on a prediction step where large time steps are used with an evaluation of the fluxes performed on the macro-cells for the region of smallest cells, and on a correction step advancing solely the region of small cells, this time with a small time step.

We demonstrate the method in a numerical framework using a vertex centered approximation, the mixed finite volume/finite element formulation.

Target applications are three-dimensional unsteady flows modeled by the compressible Navier-Stokes equations equipped with turbulence models and discretized on unstructured possibly deformable meshes. The numerical illustration involves the two above examples.

The proposed algorithm is described in Section 2. Section 3 provides some motivations of this construction. Section 4 gives several examples of applications.

Multirate time advancing by volume agglomeration 2.1 Finite-Volume Navier-Stokes

The proposed multirate time advancing scheme based on volume agglomeration will be denoted by MR and is developed for the solution of the threedimensional compressible Navier-Stokes equations. The main assumption is that the computational domain is split into computational finite volume cells such that cells intersect only by their boundaries and cover the whole computational domain. The discrete Navier-Stokes system is assembled by a flux summation Ψ i involving the convective and diffusive fluxes evaluated at all the interfaces separating cell i and its neighbors. More precisely, the finite-volume spatial discretization combined with an explicit forward-Euler time-advancing writes for the Navier-Stokes equations possibly equipped with a k -ε model:

vol i w n+1 i = vol i w n i + ∆t Ψ i , ∀ i = 1, ..., ncell,
where vol i is the volume of cell i, ∆t the time step, and w n i = (ρ n i , (ρu) n i , (ρv) n i , (ρw) n i , E n i , (ρk) n i , (ρε) n i ) are as usually the density, moments, total energy, turbulent energy and turbulent dissipation at cell i and time level t n , and ncell the total number of cells in the mesh.

In the examples given below, the accuracy of the initial scheme can be defined as a third-order spatial accuracy on smooth meshes, through the use of a MUSCL-type upwind-biased finite volume, combined with a third-order time accurate Shu multistage scheme, see [START_REF] Koobus | Some compressible numerical models for unstructured meshes[END_REF] for details.

Given an explicit -conditionally stable-time advancing, we assume that we can define a maximal stable time step (local time step) ∆t i , i = 1, ..., ncell on each node. For the Navier-Stokes model, the stable local time step is defined by the combination of a viscous stability limit and an advective one according to the following formula:

∆t i ≤ CF L × ∆l 2 i ∆l i (||u i || + c i ) + 2 γ ρ i µ i P r + µt i P rt (1)
where ∆l i is a local characteristic mesh size, u i the local velocity, c i the sound celerity, γ the ratio of specific heats, ρ i the density, µ i P r + µt i P rt the sum of local viscosity to Prandtl ratio, laminar and turbulent, and CFL a parameter depending of the time advancing scheme, of the order of unity. Using the local time step ∆t i leads to a stable but not consistent time advancing. A stable and consistent time advancing should use a global/uniform time step defined by: ∆t = min

1,ncell ∆t i .
For many advective explicit time advancing, in regions where ∆t is of the order of ∆t i , accuracy is quasi-optimal, and in other regions, the accuracy is suboptimal, due to the relatively large spatial mesh size.

Inner and outer zones

We first define the inner zone and the outer zone, the coarse grid, and the construction of the fluxes on the coarse grid, ingredients on which our MR time advancing scheme is based. For this splitting into two zones, the user is supposed to choose a (integer) time step factor K > 1. We define the outer zone as the set of cells i for which the explicit scheme is stable for a time step K∆t ∆t i ≥ K∆t, the inner zone is the set of cells i for which ∆t i < K∆t.

We shall build over the whole domain a coarse grid which should allow that:

• Advancement in time is performed with time step K∆t,

• Advancement in time preserves accuracy in the outer zone,

• Advancement in time is consistent in the inner zone.

A coarse grid is defined on the inner zone by applying cell agglomeration in such a way that on each macro-cell, the maximal local stable time step is at least K∆t. Agglomeration consists in considering each cell and aggregating to it neigboring cells which are not yet aggregated to an other one (Figure 1). Agglomeration into macro-cell is re-iterated until macro-cells with maximal time step smaller than K∆t have disappeared. We advance in time the chosen explicit scheme on the coarse grid with K∆t as time step. The nodal fluxes Ψ i are assembled on the fine cells (as usual). Fluxes are then summed on the macro-cells I (inner zone) :

Ψ I = k∈I Ψ k . ( 2
)
fine cells i macro-cell I

MR time advancing

The MR algorithm is based on a prediction step and a correction step as defined hereafter :

Step 1 (prediction step) :

The solution is advanced in time with time step K∆t, on the fine cells in the outer zone and on the macro-cells in the inner zone (which means using the macro-cell volumes and the coarse fluxes as defined in expression ( 2)) :

For α = 1, nstep outer zone : vol i w (α) i = a α vol i w (0) i + b α vol i w (α-1) i + K∆t Ψ (α-1) i (3)
inner zone :

vol I w I,(α) = a α vol I w I,(0) + b α vol I w I,(α-1) + K∆t Ψ I,(α-1) (4)
w (α) i = w I,(α) for i ∈ I (5)
EndFor α.

where w I,(α) denotes the fluid/turbulent variables at macro-cell I and stage α, (a α , b α ) are the multi-stage parameters, a 1 = 0, a 2 = 3/4, a 3 = 1/3 and b 1 = 1, b 2 = 1/4, b 3 = 2/3 for the three-stage Shu time advancing [START_REF] Shu | Efficient implementation of essentially nonoscillatory shock-capturing schemes[END_REF], and vol I is the volume of macro-cell I.

From a practical point of view, we do not introduce in our software a new variable w I to that already existing w i . In other words, the previous sequence ( 4) is actually inner zone :

vol I w (α) i = a α vol I w (0) i + b α vol I w (α-1) i + K∆t Ψ I,(α-1) (6)
where I is the macro-cell containing cell i.

On the other hand, the coarse fluxes Ψ I are not stored in a specific variable, which means that no extra storage is necessary than the one required to store the fluxes Ψ i . Indeed, after the computation of the fluxes Ψ (α-1) i (using the values of w (α-1) i

) at each stage α of the multi-stage time advancing scheme, the coarse flux Ψ I,(α-1) is evaluated according to expression (2) for each macro-cell I and then stored in the memory space allocated to Ψ (α-1) i for each cell i in the inner zone belonging to the macro-cell I.

Step 2 (correction step) :

• The unknowns in the outer zone are frozen at level t n + K∆t.

• The unknowns in the outer zone close to the inner zone, which are necessary for advancing in time the inner zone (which means those which are useful for the computation of the fluxes Ψ i in the inner zone), are interpolated in time.

• Using these interpolated values for the computation of the fluxes Ψ i in the inner zone (at each stage of the time-advancing sheme), the solution in the inner zone is advanced in time with the chosen explicit scheme and time step ∆t.

This time advancing writes:

For kt = 1, K For α = 1, nstep inner zone : vol i w (α) i = a α vol i w (0) i + b α vol i w (α-1) i + K∆t Ψ (α-1) i (8)
(outer zone : nothing is done)

EndFor α. EndFor kt.

The arithmetic complexity, proportional to the number of points in the inner zone, is therefore mastered.

Elements of analysis 3.1 Stability

The central question concerning the coarse grid is the stability resulting from its use in the computation. Considering (1), we expect that the viscous stability limit will improve by a factor four (1D) for a twice larger cell. The viscous stability limit can therefore be considered as more easily addressed by our coarsening. For the advective stability limit, we can be a little more precise. The coarse mesh is an unstructured partition of the domain in which cells are polyhedra. Analyses of time advancing schemes on unstructured meshes are available in L 2 norm for unstructured meshes, see [START_REF] Angrand | Some explicit triangular finite element schemes for the Euler equations[END_REF] [12] [START_REF] Giles | Energy stability analysis of multi-step methods on unstructured meshes[END_REF].

Here we solely propose a L ∞ analysis of the first order advection scheme. The gain in L ∞ stability can be analysed for a first-order upwind advection scheme. We get the following (obvious) lemma:

Lemma : The upwind advection scheme is positive on the mesh made of macro-cells as soon as for all macro-cell I :

∆t ||V I || < J∈N (I) ∂cell(I)∩∂cell(J) dΣ -1 cell(I)
dx where N (I) holds for the neighbouring macro-cells of I.

The application of an adequate neighboring-cell agglomeration, producing large macro-cells of good aspect ratio will produce a K-times larger stability limit.

Accuracy

In contrast to more sophisticated MR algorithms, the proposed method has not a rigorous control of the accuracy. Let us however remark that the generic situation involves variable-size meshes, which limits the unsteady accuracy on small scales propagation, already before applying the MR method.

However the two following remarks tend to show that the scheme accuracy is conserved:

-the predictor step involves simply a sum of the fluxes and the formal accuracy order is kept, with a coarser mesh size.

-still during the predictor step, if we assume that the mesh is reasonably smooth, then the CFL applied in the inner part near the matching zone will be close to the explicit CFL (applied on the outer part near the matching zone) and therefore accuracy is high.

Under these conditions, the effect of the corrector step will be just improving the result.

In practice, most of our experiments will involve a comparison between the three-stage Shu explicit time advancing and the MR algorithm using the same explicit scheme.

Efficiency

Two-level MR: The proposed method depends on only one parameter, the ratio K between the large and small time step. Considering a mesh with N vertices, a short loop on the mesh will produce the function K → N inner (K) ≤ N which gives the number of cells in the inner region for K.

If CP U ExpN ode (∆t) denotes the CPU per node and per time step ∆t of the underlying explicit scheme, a model for the MR cpu per ∆t would be

CP U M R(K) (∆t) = N K + N inner (K) × CP U ExpN ode (∆t) (10)
to be compared with the explicit case:

CP U Expli (∆t) = N × CP U ExpN ode (∆t).
We shall call the expected gain the ratio:

Gain = CP U Expli (∆t) CP U M R(K) (∆t) = 1 1 K + N inner (K) N .
The above formula emphasizes the crucial influence of a very small proportion of inner cells.

Remark 1: In most other multirate methods, the phase with a larger time-step does not concern the inner region and then their gain would be modelled by:

Gain = 1 1 KN (N -N inner (K)) + N inner (K) N .
Both gains are bounded by N/N small (K) and show that this ratio has to be sufficiently large.

Three-level MR: Multirate strategies are supposed to extend to more than two different time step lengths while keeping a reasonable algebraic complexity. Let us examine the case of three lengths, namely ∆t, K∆t, K 2 ∆t. It is then necessary to generate two nested levels of agglomeration in such a way that Grid1 is stable for CF L = 1, Grid2 is stable for CF L = K, Grid3 is stable for CF L = K 2 . While a two-rate calculation would involve a prediction-correction based on Grid3 and Grid2 ∪ Grid1 and a factor K: -prediction on Grid3 ∪ Grid2 ∪ Grid1, -correction on Grid2 ∪ Grid1, in a three-rate calculation, the correction step is replaced by two corrections and the whole process writes: -prediction on Grid3 ∪ Grid2 ∪ Grid1, -correction on medium part of Grid2 ∪ Grid1, -correction on inner part of Grid1, but this replacement is just the substitution (on a part of the mesh) of a single-rate advancing by a two-rate one and therefore can carry a higher efficiency (the smallest time step is restricted to a smaller inner zone Grid1). In contrast to other MR methods, we have a (second) duplication of flux assembly on the inner zone:

CP U M R(K) (∆t) = N 1 + N 2 + N 3 K 2 + N 1 + N 2 K + N 1 (K) ×CP U ExpN ode (∆t) (11
) to be compared with [START_REF] Gear | Multirate linear multistep methods[END_REF], noting that N = N 1 + N 2 + N 3 . However, this increment by a factor 1 + K -1 + K -2 remains limited (111% for K = 10).

Impact of our MR complexity on mesh adaption

An important impact of MR methods is the increment of accuracy order in unsteady mesh adaption methods. We now check that the proposed MR indeed improves mesh adaption accuracy order. Let us consider the spacetime mesh used by a time advancing method. A usual time advancing uses the Cartesian product {t 0 , t 1 , ...t N } × spatial mesh as space-time mesh. The space-time mesh is a measure of the computational cost since the discrete derivatives are evaluated on each node (t k , x k ) of the N time × N space nodes of the space-time mesh.

In [START_REF] Coudière | Accuracy barriers in mesh adaptation[END_REF] is proposed an analysis which determines the maximal convergence order (in terms of number of space-time nodes) which can be attained on a given family of mesh. This analysis is useful for evaluating mesh-adaptive methods. For example, in [START_REF] Belme | Time accurate anisotropic goaloriented mesh adaptation for unsteady flows[END_REF], a 3D mesh adaptive method for computing a traveling discontinuity has a convergence order α not better than α max = 8/5, according to

error = O(N -α/4 st ), (12)
4 being the space-time dimension. The purpose of this section is to show that replacing the usual time-advancing by the MR algorithm will indeed improve the maximal convergence order of a mesh adaption method defined as in [START_REF] Belme | Time accurate anisotropic goaloriented mesh adaptation for unsteady flows[END_REF].

Let us concentrate on the calculation of a planar horizontal discontinuity vertically travelling in a 3D spatial cubic domain and its approximation by a second-order accurate P 1 finite-element method. The approximation error is of first order (space and time) near the discontinuity. Then starting from an uniform space-time mesh, dividing mesh size and time step by a factor 4 will divide by 4 the error, at the cost of N 2 /N 1 = 256 times more space-time nodes. We verify that, according to [START_REF] Giles | Stability analysis of a Galerkin/Runge-Kutta Navier-Stokes discretisation on unstructured tetrahedral grids[END_REF], this is an order of 1.

In order not to go into many details given in [START_REF] Belme | Time accurate anisotropic goaloriented mesh adaptation for unsteady flows[END_REF] , we ask the reader to believe us if we say that an anisotropic mesh adaptor is able to improve the spatial mesh of each time level in such a way that starting from a constant time step and adapted space meshes at each time level, and dividing the time step by 4 and using adapted meshes with only 8 times more nodes than before will produce a 4 times smaller error. In short, this is due to the fact that mesh size normal to the discontinuity can be divided by 4, ensuring a 4 times smaller error with only a small fraction of N space , typically of the order N 2/3 space , the number of points lying on the planar discontinuity. The resulting performance an order of 8 5 is poor with respect to second order due to the uniform division by 4 of the time step. With a MR time advancing, the fourtimes smaller time step can be restricted to the N space for a usual MR algorithm). Both formulas, for N space large, give the second order (α = 2) convergence.

Parallelism

Implementation of MR

The proposed method is experimented with a parallel MPI software relying on mesh partitioning. In a preprocessing phase, the cell-agglomeration is applied at run time inside each partition, which saves communications. Further, it is done once for the whole computation, while fluctuations of the inner zone will be taken into account by using solely the agglomeration of the cells in inner zone at each time level. Since our purpose is to remain with a rather simplified modification of the initial software, we did not modify the communication library in order to restrict the communications to the inner zone when it is possible (i.e. in the correction step). Due to this, the complexity of the correction step is not strictly of order of the number of nodes in the inner zone, while its arithmetic complexity satisfies this condition.

To synthetize, the MR algorithm involves at each time step: -an updating of the inner zone (with a volume agglomeration done once for the computation), -a prediction step which is similar to an explicit step (with a larger time step length), but with also a local sum of the fluxes in each macro-cell, -a correction step which is similar to explicit arithmetics restricted to the inner region, and for simplicity of coding, communications which are left identical to the explicit advancing, that is communications applied to both inner and outer zones.

An intrinsic extra cost of our algorithm with respect to previous multirate algorithms is the computations on the inner zone during the predictor step. An intrinsic cost is also the local sums on macro-cell which account for a very small part of the computational cost. The correction step complexity is close to an explicit advancing one on the inner zone except the two phases of time interpolation and communications. Time interpolation can be implemented with a better efficiency by applying it only on a layer around the inner zone. However, the cost of the time interpolation is a very small part of the total cost. Global communications are less costly than 10% of the explicit time step cost. If the inner zone is 30% of the domain, developing communication restricted to inner zone will reduce the communication from 10% to 3% of the explicit time step on the whole domain, which shows that the correction step would be decreased from 40% to 33% of an explicit time stepping CPU.

Load balancing

The usual Metis software can be applied on the basis of a balanced repartion of the mesh. However, as remarked in previous works (see for example [START_REF] Seny | An efficient parallel implementation of explicit multirate Runge-Kutta schemes for discontinuous Galerkin computations[END_REF]), if the mesh partition does not take into account the inner zone, then the work effort will not be balanced during the correction step. The bad work balance for correction step can be of low impact if this step concerns a sufficiently small part of the mesh, resulting in a small part of the global work. However, a more reasonable assumption is that the correction phase represents a non-negligible part of the effort. In this section we discuss the question of a partitioning taking into account the correction phase. We observe that in the proposed method the inner zone depends on the flow through the CFL condition. This means that dynamic load balancing may be necessary. However, in the important class of flow which we consider, the change in inner zone can be neglected and we consider only static balancing. An option resulting from the work of Karypis and co-workers [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF] and available in Metis is the application of a multi-constrained communication cost minimisation, with the two constraints that: -partition is balanced for the whole computational domain, which will be favourable to the prediction step, -partition is balanced for the inner part of the computational domain, which will be favourable to the correction step. This partitioning algorithm produces a compromise between: -the number of nodes in each subdomain of the global mesh, -the number of nodes in each subdomain of the inner part of the mesh, -the communications between the subdomains which are minimized in the multi-constrained algorithm.

In some particular cases, the user can specify an evident partition which perfectly balances the number of nodes in each subdomain of the global mesh and in each subdomain of the inner part of the mesh. In our experiments, we explicitly specify when it is the case and how it is performed.

Numerical tests

The MR algorithm is implemented into the MPI-parallel code AIRONUM shared by INRIA, Lemma company and university of Montpellier. A description of this tool, which solves with a mixed element/volume method on unstructured meshes the compressible Euler and Navier-Stokes equations possibly equipped with a turbulence model, can be found in [START_REF] Moussaed | Simulation of the flow past a circular cylinder in the supercritical regime by blending rans and variational-multiscale les models[END_REF] and [START_REF] Itam | Simulation of multiple blunt-body flows with a hybrid variational multiscale model[END_REF]. The explicit time-advancing is a three-stage Shu method. The mean CPU for an explicit time step per mesh node varies between 10 -7 and 4 × 10 -7 seconds according to the partition quality and the number of nodes per subdomain.

For some of the test cases, it will be interesting to compare the efficiency and accuracy of the proposed MR time advancing with an implicit calculation of the same flow over the same time interval. The implicit algorithm (BDF2) which we use combines a second-order backward differencing formula for the time quadrature and a GMRES linear solver using a Restrictive-Additive Schwarz preconditioner and ILU(0) in each partition, see [START_REF] Koobus | Some compressible numerical models for unstructured meshes[END_REF] for further details. In the cases computed with the implicit scheme, the CFL is fixed to 30 and the total number of GMRES iterations for one time step is around 20. For this CFL, the gain of an implicit computation with respect to an explicit one at CFL 0.5 is measured between 12. and 22. depending on the number of nodes per processor. The implicit scheme scalability decreases with partitions less than 10, 000 vertices while the explicit (and MR) scheme remains scalable for partitions of 5, 000 vertices. The BDF2 algorithm is second-order accurate in time and we shall use this property when estimating which time step reduction is necessary for reducing by a given factor the deviation with respect to explicit time-advancing.

Contact discontinuity

In this first example, we consider the case of a moving contact discontinuity. For this purpose, the compressible Euler equations are solved in a rectangular parallelepiped as computational domain where the density is initially discontinuous at its middle (see Figure 2) while velocity and pressure are uniform.

The uniform velocity is a purely horizontal one. As can be seen in Figure 2, small cells are present on either side of the discontinuity. The mesh moves during the computation in such a way that the nodes located at the discon-Figure 2: ALE calculation of a traveling contact discontinuity. Instantaneous mesh with mesh concentration in the middle of zoom and corresponding advected discontinuous fluid density. tinuity are still the same, and that the number of small cells are equally balanced on either side of the discontinuity. An Arbitrary Lagrangian-Eulerian formulation is then used to solve the Euler equations on the resulting deforming mesh. Our long term objective is to combine the MR time advancing with a mesh adaptation algorithm in such a way that the small time steps imposed by the necessary good resolution of the discontinuity remain of weak impact on the global computational time.

The 3D mesh used in this simulation contains 25000 nodes and 96000 tetrahedra. The computational domain is decomposed into 2 subdomains, the partition interface being defined in such a way that it follows the center plan of the discontinuity. When integer K, used for the definition of the inner and outer zones, is set to 5, 10 and 15, the percentage of nodes located in the inner zone is always 1.3%, which corresponds to the vertices of the small cells located on either side of the discontinuity. The CFL with respect to propagation is 0.5. The MR scheme with the aforementioned values of K is used for the computation. Each simulation was run on 2 cores of a Bullx B720 cluster. In Table correction step, consists of explicit time advancing on inner zone, 1.3% of the mesh (solely 78 vertices on each partition), but, due to parallel and vector inefficiency, one Shu step of it is 39% of one Shu explicit step on the whole mesh. As a result, an improvement in the efficiency of about 1.7, 2.0 and 2.2 is observed when K is set to 5, 10 and 15, respectively, unstead of the 4.7, 8.8 and 12.5 predicted by the theory. We now consider a less academic example, the supersonic flow around a probe model for Exomars (see for example [START_REF] Guelhan | Experimental study of the dynamic stability of the exomars capsule[END_REF]). The Reynolds number is 1 million with respect to probe diameter. Delicate features in this simulation are a separation arising on a highly curved wall and a relatively large recirculation zone at afterbody. A hybrid RANS-LES calculation [START_REF] Itam | Simulation of multiple blunt-body flows with a hybrid variational multiscale model[END_REF], which brings more information than pure RANS does, is used for this simulation. The mesh involves 4, 380, 000 cells and the smallest element thickness is 2.10 -5 . The smallest cells are only a few and are concentrated near the largest radius of the probe. The mesh is partitioned with the usual option of Metis, into 192 subdomains. A sketch of the flow computed is presented in Figure 3. Two typical values of K are shown (see Table 2). For K = 10, only 657 cells are in the inner zone (located in the high curvature region of the boundary layer where few very small cells are present. The gain in efficiency is already 2.93. This illustrates the interest of the proposed MR when the mesh involves only a very small number of small cells. For K = 40, yet only 1752 cells are in the inner zone. The gain is 3.82, much smaller than the scalar theoretical evaluation. Further, for higher K the gain does not improve. In this particular case, the limitation is not related to an increased number of inner cells, but to the inefficient parallel treatment of the inner zone, due to a poor and heterogeneous partition of it, since the 1752 cells are distributed on 16 subdomains, the other 176 containing none of these cells. In order to Figure 4: Spatial probe at Reynolds number 1 million. Zoom of the lift curves obtained with explicit, implicit and MR schemes. evaluate the impact on accuracy, we consider as initial condition the result of a preliminary calculation giving an already installed vortex shedding. Then we continue the computation for one shedding cycle with either the explicit time advancing, or the MR one, or the BDF2 implicit time advancing (Figure 4). For 20 shedding cycles, the BDF2 run would produce a 0.4% deviation, which is probably reasonable for measuring bulk coefficients and probably not enough accurate for high frequency fluctuations. The deviation of the MR run is 125 times smaller for K = 40. According to the second-order accuracy of the BDF2, the CFL should be reduced to 2.7 if the user needs the same level of accuracy as MR. Then the measured gain is no longer in favor of the implicit algorithm (see Table 2).

Spatial probe

Circular cylinder at very high Reynolds number

The third application concerns the simulation of the flow around a circular cylinder at Reynolds number 8.4 × 10 6 . As for the previous benchmark, the computational domain is made of small cells around the body in order to allow a proper representation of the very thin boundary layer that occurs at such a high Reynolds number. On the other hand, the same hybrid RANS/VMS-LES model as that of the previous benchmark is used to compute this flow, which implies again that both the fluid and turbulent variables need to be advanced by the time integration scheme, and therefore also the MR method. Figure 5 depicts the Q-criterion isosurfaces and shows the very small and complex structures that need to be captured by the numerical and the turbulence models, which renders this simulation very challenging.

The mesh used in this simulation contains 4.3 million nodes and 25 million tetrahedra. The smallest cell thickness is 2.5.10 -6 . The computational domain is decomposed into 768 subdomains. When integer K, used for the definition of the inner and outer zones, is set to 5, 10 and 20, the percentage of nodes located in the inner zone is 15%, 19% and 24%, respectively (see Table 3). For each simulation, 768 cores were used on a Bullx B720 cluster, schemes with different values of K are given in Table 3. One can observe that the efficiency of the MR approach is rather moderate. The cost of the correction step is indeed relatively high compared to the prediction step. This is certainly due to an important number of inner nodes (which implies also a moderate theoritical scalar gain) and a non uniform distribution of these nodes among the computational cores. An implicit simulation, with a CFL number set to 30, was also performed. An important gain is observed compared to the MR case, but at the cost of a degradation of the accuracy (see Table 3 and Figure 6) . This is probably related to the very small scale fluctuations which arise at this Reynolds number (Figure 5). They need be captured with a rather accurate time advancing. A 1% deviation after a shedding cycle may become 20% after 20 cycles and deteriorate the prediction of bulk fluctuations. In order to obtain the same level of error, the implicit time advancing, which is second-order accurate in time, should be run with a CFL of 2, with a gain of only 1.5. 4.4 Tandem cylinders with a separation distance of 3.7D

We turn now to a rather well-known benchmark test case, the calculation of a flow around a tandem cylinders at Reynolds number 1.66 × 10 5 . This was a test case of an AIAA workshop, see [START_REF] Lockard | Summary of the tandem cylinder solutions from the benchmark problems for airframe noise computations-I[END_REF]. It is a challenging computation since several complex flow features need to be captured around and between the two bodies (stagnation zones, boundary layers, shear layers, separations, laminar-turbulent transition, recirculations, vortex sheddings, wakes). In order to illustrate the quality of resolution with 16 milion vertices, the Q-criterion isosurfaces are shown in Figure 7. It shows the complex flow features and the very small structures that need to be captured by the numerical model and the turbulence model, which renders this simulation particularly challenging. Further information concerning the comparison between computation and experiments are available in [START_REF] Itam | Simulation of multiple blunt-body flows with a hybrid variational multiscale model[END_REF]. When integer K, used for the definition of the inner and outer zones, is set to 20, the percentage of nodes located in the inner zone is 35% (see Table 5), which is a little too large for MR efficiency. The theoretical scalar gain is already rather small in this case, and the MR option turns out to show only a gain of 2 when 768 processors are used. K nproc.

Tandem cylinders with a large separation distance

We turn now to the flow around a tandem cylinder with a larger separation distance. As mentioned in the previous case, several complex flow features need to be captured around and between the two bodies. Small cells are necessary for a proper prediction of the very thin boundary layers, which implies very small global time steps so that classical explicit calculations become very costly. The case chosen here presents a supplementary difficulty since the distance between the two cylinders is 12 times a cylinder diameter. This kind of flow is of interest typically to measure the interaction between two risers in the sea. It is however a rather unusual test case and we have only measures of the horizontal mean velocity along the axis between the cylinders. The Reynolds number is 1.4 10 5 . We use a hybrid turbulence model based on RANS and VMS-LES approaches [START_REF] Itam | Simulation of multiple blunt-body flows with a hybrid variational multiscale model[END_REF].

As for numerical simulations, we did not find any similar computation in the litterature. In order to capture accurately the interaction between the cylinders, we use a superconvergent approximation of advective terms, able to provide up to fifth order accuracy on a Cartesian mesh. The local mesh between the cylinders is built as Cartesian in most of the interval. The mesh extends in the streamwise direction from x = -15. to x = 27., that is an extension of 42 diameters. The total number of vertices is 12.4 millions.

In order to illustrate the quality of resolution, vorticity contours are depicted in Figure 8. It shows the complex flow features and the very small structures that need to be captured by the numerical model and the turbulence model. The computational domain is decomposed into 192 and 768 subdomains, and as many cores on a Bullx cluster were used to perform these computations. When integer K, used for the definition of the inner and outer zones, is set 20, the percentage of nodes located in the inner zone is 35% (see Table 5). The CPU times for the explicit and MR schemes are shown in Table 5.

Conclusion

A new simplified multirate strategy for unstructured finite volume CFD is proposed in this work. The motivation of this research is two folds. First, we are interested by increasing the efficiency of accurate unsteady simulations. Indeed, the very high Reynolds number hybrid simulations can be computed with implicit time advancing for maintaining a reasonable cpu. But in many cases this is done with an important degradation of the accuracy with respect to smaller time steps on the same mesh. Second, with the arising of novel anisotropic mesh adaptation methods, the complexity of unsteady accurate computations with large and small mesh sizes needs to be mastered with new methods.

The proposed method is based on control volume agglomeration, and relies on: -a prediction step where large time steps are used and where the fluxes for the smaller elements are evaluated on macro-cells for stability purpose, -a correction step in which only the smaller elements of the so-called inner zone are advanced in time with a small time step. An important interest of the method is that the modification effort in an existing explicit unstructured code is very low.

K nproc. Theoretical Measured

Preliminary interesting results are given. They show that the proposed MR strategy can be applied to complex unsteady CFD problems such as the prediction of three-dimensional flows around bluff bodies with an hybrid RANS/LES turbulence model.

Simulations with Reynolds numbers ranging from 1.66 × 10 5 to 8.4 millions, representative of problems that can be encountered in industrial applications, are performed. For the considered flow calculations, the fully explicit option is still usable but of high computational cost.

All the numerical experiments are parallel computed with MPI. This allows to identify the main difficulty in obtaining high computational gain, which is related with the parallel efficiency of the computations restricted to the inner zone.

The case of the simplified mesh adaptive calculation of a moving shock is first studied, as a preliminary test for mesh adaptation. The particular case of a mesh of large cells involving, due to delicate mesh generation, a very small number of very small cells is then considered. We observed that the proposed MR strategy offers a superior efficiency in that case.

Thanks to the use of an explicit Runge-Kutta time advancing, the time accuracy of the MR scheme remains high and the dissipation remains low, as compared with an implicit computation. Only very small time-scales are lost with respect to a pure explicit computation.

In some numerical experiments, CPU gains are obtained with a usual mesh partition. However, the workload should be equally shared for both Step1 (prediction step) and Step2 (correction step) of the MR algorithm. This has been done either with a multi-constrained version of Metis, or, in some case, by hand in an exactly optimal way. In both case, however, communication is increased.

In all the cases, the proposed MR is very stable in practice and the loss of accuracy with respect to an explicit scheme is very low, in contrast to implicit BDF2-based calculations, although we applied the implicit scheme with a CFL of 30, not much larger than with the MR calculation (CFL= 10 for K = 20 in our simulations). Implicit accuracy is limited not only by the intrinsic scheme accuracy but also by the conditions required to achieve greater efficiency which involve a sufficiently large time-step and a short, parameter dependant, convergence of the linear solver performed in the time advancing step. In contrast, explicit and MR computations are parameter safe, and the accuracy of the MR method is optimal in regions complementary to the inner zone, that is, in the case of vortex shedding flow simulations, in regions where it can be necessary to propagate accurately vortices, for example from the first cylinder to the second one in the case of the tandem cylinders.

However, the main output of the experiments is already announced by the complexity analysis: a compulsory condition for a good speedup is that for a K ≥ 10, the inner zone is less than 10% of the complete domain. In the selected case, the speed up is about 2.

Due to its simplicity, the proposed method can be easily extended to several MR layers corresponding to different time step stability regions in order to separate, for example in the case of three layers, very small scales from intermediate ones, and intermediate scales from larger ones. In a sequel of this work, we plan to combine a higher number of layers with an unsteady mesh-adaptive algorithm.

Figure 1 :

 1 Figure 1: Sketch (in 2D) of the agglomeration of 4 cells into a macro-cell. Cells are dual cells of triangles, bounded by sections of triangle medians.

2 / 3 space

 23 smaller spatial cells and a two-times smaller time step is applied to the others. The amplification of CPU time is then 16 + 32 * N -1/3 space for the proposed algorithm and (16 * (1 -N

Figure 3 :

 3 Figure 3: Spatial probe at Reynolds number 1 million. Q criterion.

Figure 5 :

 5 Figure 5: Circular cylinder at Reynolds number 8.4 × 10 6 . Instantaneous Q-criterion isosurfaces (coloured with velocity modulus).

Figure 6 :

 6 Figure 6: Circular cylinder at Reynolds number 8.4 × 10 6 . Zoom of the lift curves obtained with explicit, implicit and MR schemes.

Figure 7 :

 7 Figure 7: Tandem cylinders with distance 3.7D: instantaneous Q-criterion isosurfaces (coloured with velocity modulus).

Figure 8 :

 8 Figure 8: DDES calculation around a tandem cylinders with distance 12D. Velocity contours.

Table 1 :

 1 1, CPU times (prediction step / correction step) are given for the MR approach and different time step factors K. The ALE propagation of a contact discontinuity: Time step factor K, CPU of the explicit scheme per explicit time-step ∆t and per node, percentage of nodes in the inner region, theoretical gain in scalar mode, CPU of the prediction step per time-step K∆t, CPU of the correction step per time-step K∆t, and measured parallel gain.

	K N small (K)/N	Expected	CPU	CPU	Measured
		(%)	gain	pred. step correc. step	gain
			(theoretical) (s/K∆t)	(s/K∆t)	(parallel)
	5	1.3	4.7	0.124	0.244	1.7
	10	1.3	8.8	0.124	0.482	2.0
	15	1.3	12.5	0.124	0.729	2.2

Table 2 :

 2 Spatial probe: Time step factor K, CPU of the explicit scheme per explicit time-step ∆t and per node, percentage of nodes in the inner region, theoretical gain in scalar mode, CPU of the prediction step per timestep K∆t, CPU of the correction step per time-step K∆t, measured parallel gain, and relative lift deviation with explicit scheme after a shedding period.

Table 3 :

 3 Circular cylinder at Reynolds number 8.4 × 10 6 Time step factor K, CPU of the explicit scheme per explicit time-step ∆t and per node, percentage of nodes in the inner region, theoretical gain in scalar mode, CPU of the prediction step per time-step K∆t, CPU of the correction step per time-step K∆t, measured parallel gain, and relative error for the explicit, MR and implicit BDF2 time advancing. UP holds for usual partition, MCP for Metis multi-constrained partition, R for analytic radial optimal partition.

Table 4 :

 4 Tandem cylinder with distance 3.7 D. Time step factor K, number of processors, percentage of nodes in the inner region, theoretical gain in scalar mode, and measured parallel gain. UP holds for usual partition and MCP for multi-constrained partition.

			N small N	Expected	Measured
			(%)	gain	gain
				(theoretical) (UP/MCP)
	20	192	35	2.50	1.12/1.77
	20	768	35	2.50	1.02/2.0

Table 5 :

 5 Tandem cylinder with distance 12D. Time step factor K, number of processors, CPU of the explicit scheme per explicit time-step ∆t and per node, percentage of nodes in the inner region, theoretical gain in scalar mode, CPU of the prediction step per time-step K∆t, CPU of the correction step per time-step K∆t, and measured parallel gain. UP holds for usual partition, R for optimal radial partition.

					Measured
			gain	gain/UP	gain/R
	40	288	3.7	1.09	2.11
	60	288	3.3	1.11	2.07
	40	384	3.7	1.10	2.12
	60	384	3.3	1.12	2.07
	40	768	3.7	1.35	2.42
	60	768	3.3	1.35	2.36

Acknowledgments

This work has been supported by French National Research Agency (ANR) through "Modèle numérique" program (projet MAIDESC n 0 ANR-13-MONU-0010). This work was granted access to the HPC resources of CINES under the allocations 2017-A0022A05067 and 2017-A0022A06386 made by GENCI (Grand Equipement National de Calcul Intensif).


Amira Filali
email: amira.filali@univ-nantes.fr

Vincent Ricordel
email: vincent.ricordel@univ-nantes.fr

Nicolas Normand
email: nicolas.normand@univ-nantes.fr

Contexte et état de l'art : Tree-Structured Point-Lattice Vector Quantization for 3-D Point Cloud Geometry Compression

This paper deals with the current trends of technologies which aim at representing efficiently 3D point cloud, they have become a popular challenge for processing, storing and conveying the data independently of how it was captured. Compressing attributes of 3D point clouds such as geometry, colours or normal directions remain challenging problem, since these signals are unstructured. We propose an adaptive Tree-Structured Point-Lattice Vector Quantization (TSPLVQ), resulting in hierarchically structured 3D content, to improve compression performance for static point cloud.The novelty of the proposed approach lies in adaptively selecting the optimal quantization scheme for the point cloud geometric data, such that its intrinsic correlations can be better exploited.Two quantization modes are dedicated to project recursively the 3D point clouds into a series of embedded truncated cubic lattices. At each step of the process, the optimal quantization mode is chosen according to a rate-distortion criterion in order to achieve the best trade-off between coding rate and geometry distortion, such that the compression flexibility and performance can be greatly improved. To achieve rate-distortion optimization, the core of the proposed method relies on an accurate estimation of the distortion inside each cubic voronoï, between fitting plane and cloud points, taking into account the visual rendering of the 3D object. We experimentally evaluated the interest of the proposed method for geometry compression and visualization tasks. 3D point cloud, Lattice Vector Quantization, Rate-Distortion Optimization, Compression Advances in 360-degree video and virtual/augmented reality fields, provide many fast and reliable cameras, laser scanners, etc. which are rapidly spreading as devices for capturing the 3D data in different ways providing high definition 3D content. The highly increased quality and quantity of the captured data with additional reflectance and color produces very dense 3-D point cloud (PC) data to realistically represent the real world enabling immersive forms of interaction, navigation and communication. The common processing such as storage, transmission from the source data to the diverse platforms, or visualization, etc. become complex, or just impossible. In many industrial domains, it is desirable to transmit this data via restricted bandwidth networks. PCs provide a significant barrier for mass market applications. Thus, MPEG has issued a Call for Proposal (January 2017) for PC compression for its intended applications. Multiple methods of compression have been developed. In [1], Morell et al. have proposed an approach based on plane extraction which represents the points of each scene plane as a Delaunay triangulation and a set of points/ area information. The concave hull of a plane is computed to find the borders and get a polygon representing it, an edge reduction is then performed to reduce the number of points needed. Many methods (available in the PCL [2]) are based on Octree decomposition to reduce the amount of data content. The principle is to build a large cuboid to contain all the PC data in a 3-D space. Then the algorithm splits the large cuboid into smaller 3D voxels of same size. For each voxel, the algorithm calculates the mean value of the 3-D coordinates of all the points. The mean value is used to represent all the points that lie in that voxel. In this way, the point cloud can be compressed while keeping the geometric details. In the context of classical images and video encoding applications, Ricordel et al. [4] have proposed a new vector quantizer (TSLVQ) based on truncated lattice embedding. They have investigated its complete design with: the lattice truncation, the multi-stage procedure of quantization, the unbalanced tree-structured codebook design according to a distortion v.s. rate trade-off. The purpose of our method is to reduce the amount of data of a

Travail proposé (Décrire l'objectif du travail et indiquer clairement le problème/défi technique étudié) :

3D points set but trying to preserve as much information as possible. To the best of our knowledge, we are the first, in the context of 3-D point clouds compression, to propose applying adaptive TSPLVQ with two quantization schemes (2x2x2 or 3x3x3 Voronoï splitting) depending on a rate-distortion criterion which takes into consideration the distortion between points and reconstructed planes associated to points in each cubic Voronoï taking into account the visual rendering of the 3D object and the tree coding cost.

Background: Tree-Structured LVQ (TSLVQ) is a gathering of many quantization approaches where the quantization is processed through a decision tree. Its benefits are reduced computation with the use of simpler sub-codebooks, and a structure adapted to progressive representation. Here in 2D, the scaling factor equals b=3 such as a cubic Voronoï is splitted in (3x3) [4].

Proposed Work: The point cloud is highly complex and unstructured, such that an organization strategy of the point cloud can directly affect the coding performance. The proposed approach, based on the embedding of truncated lattices, permits hierarchical description of the 3-D point cloud through an adaptive tree-structure codebook simplifying its representation and compression. One single splitting mode, such as octree structure (where a cubic voronoï is splitted in 2x2x2), is insufficient to account for the characteristics of all kinds of point cloud geometries. Instead of a single splitting mode, in this work two slitting modes (2x2x2 versus 3x3x3) are employed to better map the point cloud to 3-D grids. We also have to scale and truncate the lattices suitably providing the best rate-distortion trade-off. The basic quantization procedure used for our method TSPLVQ based on TSLVQ is illustrated in Fig. 2: (0): Initialization: The points set is first normalized to fit inside the root voronoï cell by centering the points and then by scaling data with the factor:

, where is the maximal euclidean distance between a point and the center.

(1) and (2): An input vector x (its components are 3D point coordinates) inside a given cubic voronoï is scaled by the b factor and then shifted, the shifting U and factor b set up the splitting of the cube either in 2x2x2 (b=2, U=(1/2,1/2,1/2)) or in 3x3x3 (b=3, U=(0,0,0)).

(3): The fast quantization algorithm is then used to produce by rounding the corresponding reproduction vector Y.

(4): The output vector is centered to permit the next quantization level. So to initialise the construction of the tree-strutured codebook that will represent the 3D point cloud geometry, all the cloud points are projected into a first cubic Voronoï (namely the tree root). Next, the basic quantization structure is repeated iteratively from each point of the cloud. At each loop of the greedy process, the cube to split and the adapted splitting method, have to be chosen according to a ratedistortion criterion. Exactly a Lagrangian optimization is then employed at this level, the multiplier associated to each Voronoï (namely tree node) is computed, such as when a node is splitted, the increase in rate is calculated in terms of tree coding cost (here: the entropic cost of the nodes is considered), and the decrease in distortion in terms of geometric distortion. To assess objectively this distortion, we put in competition 2 metrics:

-Point-to-point metric, as the square distance between the lattice point and all assigned cloud points X TSLVQ [4] aims at using a hierarchical set of embedded lattices which is obtained such as it is possible to embed a lower scale truncated lattice into a cell of the next higher scale truncated lattice. So a scaling factor b between successive truncated lattices of the hierarchy has to be set (see Fig. 1). The use of the cubic lattice permits optimal embedding and fast quantization. A ratedistortion optimization is used in order to choose which node to split in the iterative tree design.

F = 1 E ma x E ma x Y i CORESA 2018
Poitiers, 12-14 Novembre 2018

Résultats

For the experimental results, we compared different voronoï splitting strategies: for the first case only 2x2x2 splitting mode (octree) is used, for the second case only 3x3x3 mode, and for the third case the 2 splitting modes (2x2x2 and 3x3x3) have been put in competition. The performance comparisons between the 3 different schemes are shown in Fig. 3. Note that the distortion is assessed according to the 2 metrics that are considered: point-to-point or point-to-plane. For a fair comparison, we set up a parameter to our method, such as the rate in function of leaves number. It can be observed that the splitting strategies leads to different rate-distortion performances, for example, the 2x2x2 splitting scheme is more progressive and requiring more loops for the same number of obtained leaves compared to the two other schemes. The slopes of the curve is a key-point because when it becomes flat, it indicates that further splittings are not necessary. We can note that the slopes behaviors depend on the splitting method and of the used metric, and again the 2x2x2 method seems more interesting than the 3x3x3 method. When we put in competition the 2 splitting methods (called adaptive method), the metric choice is important: with P2Ppoint the 2 splitting methods are combined and the slopes are not flat, with P2Plane no combination are made. At this level of the work, it also seems that the P2Plane metric does not show advantage in term of objective quality evaluation. But it is also important to consider the rendering of the encoded point cloud, Fig. 4 presents the rendering coded 3D point cloud in order to analyze and compare their visual quality. Our adaptive method preserves more visual details in the complex areas (for instance the head of the dragon) compared to the two other splitting schemes, whereas it produces less points to represent smooth surfaces (for instance the dragon body surface). Hence our proposed adaptive method using point-to-plane distortion seems more linked with the subjective quality of the rendered decoded object. Page / 3 4

inside the corresponding Voronoï cell :

-Point-to-plane metric [3,5], as the square distance between the lattice point to normals to the plane associated to the cloud points X inside the corresponding cell :

For a given loop of the tree growing process, the rate-distortion optimization is performed locally. It aims at choosing among all the possible splitted nodes, the one which offers the maximal decrease in distortion, and the minimal increase in rate. The codebook process is stopped for instance when the bit budget is reached, and the final tree-structure is then unbalanced. Each occupied leaf node corresponds to a point output. The location of the output representant is set to the average value of the PC points contained within the leaf node. C

i D P2Point (X, Y i ) = ∑ X∈C i d (X, Y i ) Y i N X i C i D P2Plane (X, Y i ) = ∑ X∈C i | | V d (X, Y i ) ⋅ N X i | |

Conclusion et perspectives :

Références (5 max.):

Our main research investigated into a compression method of the geometric data of PCs. We have proposed a lossy compression method based on TSPLVQ. Our method innovates because it introduces three splitting strategies (2x2x2 versus 3x3x3) for the cubic voronoï, and by the use of two metrics (P2Point vs. P2Plane).

Experimental results show the different behaviours of the rate-distortion curves and different qualities of the rendered contents, their analysis explains the TPSLVQ coding improvement in this context. As future work, we would like to include more accurate geometric primitives and other attributes (as colour) to better characterise the content of a voronoï cell, in order to better render the decoded 3D objects. We also plan to improve the rate computation for the coding of the tree-structure.

Fig. 1 .

 1 Fig.1. Hierarchy principle corresponding to the cubic lattice.

Fig. 3 .

 3 Fig.3. TSPLVQ method using P2Point vs. P2Plane distortion metricOf « Dragon.ply ».

Fig. 2 .

 2 Fig.2. Initialization (left), basic procedure of TSPLVQ where X is the input vector, Y the reproduction vector, b the scaling factor (right).

Fig. 4 .

 4 Fig.4. Visual rendering result of Dragon point cloud: (a), (b) only 2x2x2 splitting using P2Point and P2Plane distortion respectively; (c), (d) only 3x3x3 splitting using P2Point and P2Plane distortion respectively; (e), (f) 2x2x2 vs.3x3x3 splitting using P2Point and P2Plane distortion respectively.