Lagrange-flux schemes: reformulating second-order accurate Lagrange-remap schemes for better node-based HPC performance

In a previous paper, we have achieved the performance analysis of staggered Lagrange-remap schemes, a class of solvers widely used for Hydrodynamics applications. This paper is devoted to the rethinking and redesign of the Lagrange-remap process for achieving better performance using today's computing architectures. As an unintended outcome, the analysis has lead us to the discovery of a new family of solvers -- the so-called Lagrange-flux schemes -- that appear to be promising for the CFD community.

: "Legacy" staggered Lagrange-remap scheme: thermodynamical variables are located at cell centers (circles) whereas velocity variables are located at cell nodes (squares). [2,3,4], also referred to as remapped Lagrange solvers that apply a remap step on a reference (Eulerian) mesh after each Lagrangian time advance. Legacy codes implementing remapped Lagrange solvers usually define thermodynamical variables at cell centres and velocity at cell nodes (see figure 1). In Poncet et al. [1], we have achieved a multicore node-based performance analysis of a reference Lagrange-remap hydrodynamics solver used in industry. By analyzing each kernel of the whole algorithm, using roofline-type models [5] on one side and refined Execution Cache Memory (ECM) models [6], [7] on the other side, we have been able not only to quantitatively predict the performance of the whole algorithm -with relative errors in the single digit range -but also to identify a set of features that limit the whole performance.
This can be roughly summarized into three points: 1. For typical mesh sizes of real applications, spatially staggered variables involve a rather big amount of communication to/from CPU caches and memory with low arithmetic intensity, thus lowering the whole performance; 2. Usual alternating direction (AD) strategies (see the appendix in [8]) or AD remapping procedures also generate too much communication with a loss of CPU occupancy. 3. For multimaterial flows using VOF-based interface reconstruction methods, there is a strong loss of performance due to some array indirections and noncoalescent data in memory. Vectorization of such algorithms is also not trivial.
From these observations and as a result of the analysis, we decided to "rethink" Lagrangeremap schemes, with possibly modifying some aspects of the solver in order to improve nodebased performance of the hydrocode solver. We have looked for alternative formulations that reduce communication and improve both arithmetic intensity and SIMD (Single Instruction -Multiple Data) property of the algorithm. In this paper, we describe the process of redesign of Lagrange-remap schemes leading to higher performance solvers. Actually, this redesign methodology also gave us ideas of innovative eulerian solvers. The emerging methods, named Lagrange-flux schemes, appear to be promising in the extended computational fluid dynamics community.
The paper is organized as follows. In section 2, we first formulate the requirements for the design of better Lagrange-remap schemes. In section 3 we give a description of the Lagrange step and formulate it under a finite volume form. In section 4 we focus on the remap step which is reformulated as a finite volume scheme with pure convective fluxes. This interpretation is applied in section 5 to build the so-called Lagrange-flux schemes. We also discuss the important issue of achieving second order accuracy (in both space and time). In section 7 we comment the possible extension to multimaterial flow with the use of low-diffusive and accurate interfacecapturing methods. We will close the paper by some concluding remarks, work in progress and perspectives.

Requirements
Starting from "legacy" Lagrange-remap solvers and related observed performance measurements, we want to improve the performance of these solvers by modifying some of the features of them but under some constraints and requirements: 1. A Lagrangian solver (or description) must be used (allowing for multiphysics coupling).
2. To reduce communication, we prefer using collocated cell-centered variables rather than a staggered scheme.
3. To reduce communication, we prefer using a direct multidimensional remap solver rather than splitted alternating direction projections. 4. The method can be simply extended to second-order accuracy (in space and time).
5. The solver must be able to be naturally extended to multimaterial flows.
Before going further, let us comment the above requirements. The second requirement should imply the use of a cell-centered Lagrange solver. Fairly recently, Després and Mazeran in [9] and Maire and et al. [10] (with high-order extension in [11]) have proposed pure cell-centered Lagrangian solvers based on the reconstruction of nodal velocities. In our study, we will examine if it is possible to use approximate and simpler Lagrangian solvers in the Lagrange+remap context, in particular for the sake of performance. The fourth assertion requires a full multidimensional remapping step, probably taking into account geometric elements (deformation of cells and edges) if we want to ensure high-order accuracy remapping. To summarize, our requirements are somewhat contradictory, and we have to find a good compromise between some simplifications-approximations and a loss of accuracy (or properties) of the numerical solver.

Lagrange step
As example, let us consider the compressible Euler equations for two-dimensional problems.
For any volume V t that moves with the fluid, from the Reynolds transport theorem we have where ν is the normal unit vector exterior to V t . This leads to a natural explicit finite volume scheme in the form In (2), the superscript "L" indicates the Lagrange evolution of the quantity. Any Eulerian cell K is deformed into the Lagrangian volume K n+ 1 2 ,L at time t n+ 1 2 , and into the Lagrangian volume K n+1,L at time t n+1 . The pressure flux terms through the edges A are evaluated at time t n+ 1 2 in order to get second-order accuracy in time. Of course, that means that we need a predictor step for the velocity field u n+ 1 2 ,L at time t n+ 1 2 (not written here for simplicity).
Notations. To simplify, in all what follows we will use the notation v n+ 1 2 = u n+ 1 2 ,L .
The remapping step considers the remapping the fields U defined at cell centers K n+1,L on the initial (reference) Eulerian mesh with cells K. Let us denote R n+1,L a linear operator that reconstructs piecewise polynomial functions from a discrete field U n+1,L defined at cell centers of the Lagrangian mesh M N +1,L at time t n+1 . The remapping process consists in projecting the reconstructed field on piecewise-constant function on the Eulerian mesh, according to the intergral formula Practically, they are many ways to consider the projection operation ( ping can also be interpreted as a transport/advection algorithm, as emphasized by Dukowicz and Baumgardner [12] that appears to be better suited for SIMD treatments.
Let us now write a different original formulation of the remapping step. In this step, there is no time evolution of any quantity, and in some sense we have ∂ t U = 0, that we write We decide to split up this equation into two substeps: i) Backward convection: ii) Forward convection: Each convection problem is well-posed on the time interval [0, ∆t n ]. Let us now focus into these two steps and the way to solve them.

Backward convection in Lagrangian description
After the Lagrange step, if we solve the backward convection problem (4) over a time interval ∆t n using a Lagrangian description, we have conservative quantities. For = 1 (conservation of mass), we have showing the variation of density by volume variation. For = 2, 3, 4, it is easy to see that both velocity and specific total energy are kept unchanged is this step: Thus, this step is clearly computationally inexpensive.

Forward convection in Eulerian description
From the discrete field (U n, K ) K defined on the Eulerian cells K, we then solve the forward convection problem over a time step ∆t n under an Eulerian description. A standard Finite Volume discretization of the problem will lead to the classical time advance scheme for some interface values U n+ 1 2 , A defined from the local neighbor values U n, K . We finally get the expected Eulerian values U n+1 K at time t n+1 .
Notice that from (6) and (7) we have also thus completely defining the remap step under the finite volume scheme form (8). We find that we no more need any mesh intersection or geometric consideration to achieve the remapping process. The finite volume form (8) is now suitable for a straightforward vectorized SIMD treatment. From (8) it is easy to achieve second-order accuracy for the remapping step by usual finite volume tools (MUSCL reconstruction + second-order accurate time advance scheme for example).

Full Lagrange+remap time advance
Let us note that the Lagrange+remap scheme is actually a conservative finite volume scheme: putting (2) into (8) gives for all : that can also be written We recognize in (10) pressure-related fluxes and convective fluxes that define the whole numerical flux.

Comments
The finite volume formulation (10) is attractive and seems rather simple at first sight. But we should not forget that we have to compute a velocity Lagrange vector field v n+ 1 where the variables should be located at cell nodes to return a well-posed deformation. Moreover, expression (10) involves geometric elements like the length of the deformed edges A n+ 1 2 ,L .
Among the rigorous collocated Lagrangian solvers, let us mention the GLACE scheme by Després-Mazeran [9] and the cell-centered EUCCLHYD solver by Maire et al [10]. Both are rather computationally expensive and their second-order accurate extension is not easy to achieve.
Although it is possible to couple these Lagrangian solvers with the flux-balanced remapping formulation, it is also of interest to think about ways to simplify or approximate the Lagrange step without losing second-order accuracy. One of the difficulty in the analysis of Lagrangeremap schemes is that, in some sense, space and time are coupled by the deformation process.
Below, we derive a formulation that leads to a clear separation between space and time, in order to simply control the order of accuracy. The idea is to make the time step tend to zero in the Lagrange-remap scheme (method of lines [13]), then exhibit the instantaneous spatial numerical fluxes through the Eulerian cell edges that will serve for the construction of an explicit finite volume scheme. Because the method needs an approximate Riemann solver in Lagrangian form, we will call it a Lagrange-flux scheme.

Derivation of a second-order accurate Lagrange-flux scheme
From the intermediate conclusions of the discussion 4.4 above, we would like to be free from any rather expensive collocated Lagrangian solver. However, such a Lagrangian solver seems necessary to correctly and accurately define the deformation velocity field v n+ 1 2 at time t n+ 1 2 .
In what follows, we are trying to deal with time accuracy in a different manner. Let us come back to the Lagrange+remap formula (10). Let us consider a "small" time step t that fulfils the usual stability CFL condition. We have By making t tend to zero, (t > 0), we have A L (t/2) → A, (π ) L (t/2) → π , v(t/2) → u, (U ) → U , then we get a semi-discretization in space of the conservation laws. That can be seen as a particular method of lines ( [13]): We get a classical finite volume method with a numerical flux Φ A whose components are In (11), pressure fluxes (p ) A and interface normal velocities (u A · ν A ) can be computed from an approximate Riemann solver in Lagrangian coordinates (for example the Lagrangian HLL solver, see [14]). Then, the interface states (U ) A should be computed from an upwind process according to the sign of the normal velocity (u A · ν A ). This is interesting because the resulting flux has similarities with the so-called advection upstream splitting method (AUSM) flux family proposed by Liou [15], but the construction here is different and, in some sense, justifies the AUSM splitting.
To get higher-order accuracy in space, one can use a standard MUSCL reconstruction + slope limiting process involving classical slope limiters like for example Sweby's limiter function [17]: with β ∈ [1, 2] for achieving second order accuracy. At this stage, because there is no time discretization, everything is defined on the Eulerian mesh and fluxes are located at the edges of the Eulerian cells. This is one originality of this scheme compared to the legacy staggered Lagrange-remap scheme that has to use variables defined on the Lagrangian cells.
To get high-order accuracy in time, one can then apply a standard high-order time advance scheme (RK2, etc.). For the second-order Heun scheme for example, we have the following algorithm: 1. Compute the time step ∆t n subject to some CFL condition; 2. Predictor step. MUSCL reconstruction + slope limitation on primitive variables ρ, u and p. From the discrete states U n K , compute a discrete gradient for each cell K.
3. Use a Lagrangian approximate Riemann solver to compute pressure fluxes π n A and interface velocities u n A 4. Compute the upwind edge values (U ) n A according to the sign of (u n A · ν A ); 5. Compute the numerical flux Φ n A as defined in (12) One can appreciate the simplicity of the numerical solver compared to the legacy staggered Lagrange-remap algorithm. The complexity of the latter mainly due to various kernel (function) calls and too much communications is detailed in [1]. Here the predictor and corrector kernel functions have similar programming codes and there is no intermediate variables to save in memory.

Lagrangian HLL approximate solver
A HLL approximate Riemann solver [14] in Lagrangian coordinates can be used to easily compute interface pressure and velocity. For a local Riemann problem made of a left state U L and a right state U R , the contact pressure p is given by the formula and the normal contact velocity u by leading to simple formulas easily implementable in the second-order Heun time integration scheme.

Numerical experiments
As example, we test the Lagrange-flux scheme presented in section 5 on few one-dimensional shock tube problems. We use a Runge-Kutta 2 (RK2) time integrator and a MUSCL reconstruction with the Sweby slope limiter given in (13).
Sod's shock tube [16] The initial data defined on space interval  Two-rarefaction shock tube The second reference example is a case of two moving-away rarefaction fans under near-vacuum conditions (see Toro [14]). It is known that the Roe scheme breaks down for this case. The related Riemann problem is made of the left state   preserves the positivity of density, pressure and internal energy (see figure 5).

Performance results
In this section, we compare the new Lagrange-flux scheme to the reference (staggered) Lagrangeremap scheme in terms of millions of cell updates per second (denoted hereafter as MCUPs).
Tests are performed on a standard 2 × 8 cores Intel Sandy Bridge server E5-2670. Each core has a frequency of 2.6 GHz, and supports Intel's AVX (Advanced Vector Extension) vector instructions. For multicore support, we use the multithreading programming interface OpenMP.
In the reference staggered Lagrange-remap solver (see [1]), thermodynamic variables are defined at grid cell centers while velocity variables are defined at mesh nodes. Due to this staggered discretization and the alternating direction (AD) remapping procedures, this solver is decomposed into nine kernels. This decomposition mechanically decreases the mean arithmetic intensity (AI) of the solver.
On the other hand, the Lagrange-flux algorithm consists in only two kernels with a relative high arithmetic intensity which leads to two compute-bound (CB) kernels. In the first kernel, named PredictionLagrangeFlux(), an appropriate Riemann solver is called, face CorrectionLagrangeFlux(), is close in terms of algorithmic steps, since it also uses a Riemann solver, computes fluxes and updates the variables for the correction part of the solver.
In order to assess the scalability and absolute performance of both schemes, we present in

Dealing with multimaterial flows
Although this is not the aim and the scope of the present paper, we would like to give an out- If the geometry of the Lagrangian cells is not completely known (as in the case of Lagrangeflux schemes), anyway we have to proceed differently. A possibility is to use interface capturing (IC) schemes, e.g. conservative Eulerian schemes that evaluate the convected mass fluxes through Eulerian cell edges. This can be achieved by the use of antidiffusive/low-diffusive advection solvers in the spirit of Després-Lagoutière's limited-downwind scheme [18] of VoFire [19].
In a recent work [20], we have analyzed the origin of known artifacts and numerical interface instabilities for this type of solvers and concluded that the reconstruction of fully multidimensional gradients with multidimensional gradient limiters was necessary. Thus, we decided to use low-diffusive advection schemes with a Multidimensional Limiting Process (MLP) in the spirit of [21]. The resulting method is quite accurate, shape-preserving and free from any artifact. We show some numerical results in the two next subsections. Let us emphasize that the interface capturing strategy perfectly fits with the Lagrange-flux flow description, and the resulting schemes are really suitable for vectorization (SIMD feature) with data coalescence into memory.

Interface capturing for pure advection problems
Let us first present numerical results on a pure scalar linear advection problem. The forwardbackward advection case proposed by Rider and Kothe [22] is a hat-shaped function which is advected and stretched into a rotating vector field, leading to a filament structure. Then by applying the opposite velocity field, one have to retrieve the initial disk shape. In figure 6, we show the numerical solutions obtained on a grid 500 2 for both the passive scalar field of variable z ∈ [0, 1] and the quantity z(1 − z) that indicates the numerical spread rate of the diffuse interface. One can conclude the good behavior of the method, providing both stability, accuracy and shape preservation.

Three-material hydrodynamics problem
We then consider our interface capturing method for multifluid hydrodynamics. Because material mass fractions are advected, that is ∂ t y k + u · ∇y k = 0, k = 1, ..., N, one can use the advection solver of these variables but we prefer dealing with the conservative form of the equations ∂ t (ρy k ) + ∇ · (ρy k u) = 0 in order to enforce mass conservation (see also [23]). It is known that Eulerian interfacecapturing schemes generally produce spurious pressure oscillations at material interfaces ( [24,25]). Some authors propose locally non conservative approaches [26,27] to prevent from any pressure oscillations. Here we have a full conservative Eulerian strategy involving a specific limiting process which is free from any pressure oscillation at interfaces, providing strong robustness. This will be explained in a next paper.
The multimaterial Lagrange-flux scheme is tested on the reference "triple point" test case, found e.g. in Loubère et al. [28]. This problem is a three-state two-material 2D Riemann problem in a rectangular vessel. The simulation domain is Ω = (0, 7) × (0, 3) as described in figure 7. The domain is splitted up into three regions Ω i , i = 1, 2, 3 filled with two perfect gases leading to a two-material problem. Perfect gas equations of state are used with γ 1 = γ 3 = 1.5 and γ 2 = 1.4. Due to the density differences, two shocks in sub-domains Ω 2 and Ω 3 propagate with different speeds. This create a shear along the initial contact discontinuity and the formation of a vorticity. Capturing the vorticity is of course the difficult part to compute.
We use a rather fine mesh made of 2048 × 878 points (about 1.8M cells). On figure 8, we plot the density, pressure, temperature fields respectively and indicate the location of the three material zones. One can observe a nice capture of both shocks and contact discontinuities. The vortex is also captured accurately.

Concluding remarks and perspectives
This paper is primarily focused on the redesign of Lagrange-remap hydrodynamics solvers in order to achieve better HPC node-based performance. We have reformulated the remapping step under a finite volume flux balance, allowing for a full SIMD algorithm. As an unintended outcome, the analysis has lead us to the discovery of a new promising family of Eulerian solvers -the so-called Lagrange-flux solvers -that show simplicity of implementation, accuracy, and flexibility with a high-performance capability compared to the legacy staggered Lagrangeremap scheme. Interface capturing methods can be easily plugged for solving multimaterial flow problems. Ongoing work is focused of the effective performance modeling, analysis and measurement of Lagrange-flux schemes with comparison of reference "legacy" Lagrange-remap solvers including multimaterial interface capturing on different multicore processor architectures. Because of the multicore+vectorization scalability of Lagrange-flux schemes, one can also expect high-performance on manycore co-processors like graphics processing units (GPU) or Intel MIC. This will be the aim of next developments.