© F. De Vuyst et al., published by IFP Energies nouvelles, 2016

Motivation and Introduction

For complex compressible flows involving multiphysics phenomena like e.g. high-speed elastoplasticity, multimaterial interaction, plasma, gas-particles, etc., a Lagrangian description of the flow is generally preferred. To achieve robustness, some spatial remapping on a regular mesh may be added. A particular case is the family of the so-called Lagrange+remap schemes [1–3], 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 centers and velocity at cell nodes (Fig. 1). In Poncet et al. [4], we have achieved a multicore node-based performance analysis of a reference Lagrange-remap hydrodynamics solver used in industry. By analyzing each kernel (i.e. elementary computing function) 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” Lagrange-remap schemes, with possibly modifying some aspects of the solver in order to improve node-based performance of the hydrocode. 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 1, we first formulate the requirements for the design of better Lagrange-remap schemes. In Section 2 we give a description of the Lagrange step and formulate it under a finite volume form. In Section 3 we focus on the remap step which is reformulated as a finite volume scheme with pure convective fluxes. This interpretation is applied in Section 4 to build the so-called Lagrange-flux schemes. We also discuss the important issue of achieving second order accuracy (in both space and time). Performance comparisons between the reference Lagrange-remap and the proposed Lagrange-flux scheme are presented in Section 5. In Section 6 we comment the possible extension to multimaterial flow with the use of low-diffusive and accurate interface-capturing methods. We will close the paper by some concluding remarks, work in progress and perspectives.

1 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 their features 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 split 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 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.

2 Lagrange Step

As examples, let us consider the compressible Euler equations for two-dimensional problems. Denoting i ∈ {1,2}, p and E the density, velocity, pressure and specific total energy respectively, the mass, momentum and energy conservation equations read(1)where , , π ₂ = (p, 0) ^T , π ₃ = (0, p) ^T and π ₄ = p u . For the sake of simplicity, we will use a perfect gas equation of state . The speed of sound c is given by .

For any volume V _t that moves with the fluid, from the Reynolds transport theorem we havewhere ν is the normal unit vector exterior to . This leads to a natural explicit finite volume scheme in the form(2)

In (2), the superscript “L” indicates the Lagrange evolution of the quantity. Any Eulerian cell K is deformed into the Lagrangian volume at time , and into the Lagrangian volume K^n+1,L at time tⁿ⁺¹. The pressure flux terms through the edges A are evaluated at time in order to get second-order accuracy in time. Of course, that means that we need a predictor step for the velocity field at time (not written here for simplicity).

Notations

To simplify, in all what follows we will use the notation .

3 Rethinking the Remapping Step

The remapping step considers the remapping fields defined at cell centers K ^n+1,L on the initial (reference) Eulerian mesh with cells K. Let us denote a linear operator that reconstructs piecewise polynomial functions from a discrete field U ^n+1,L defined at cell centers of the Lagrangian mesh at time t ⁿ⁺¹. The remapping process consists in projecting the reconstructed field on piecewise-constant functions on the Eulerian mesh, according to the integral formula(3)

Practically, there are many ways to consider the projection operation (3). One can assemble elementary projection contributions by computing the volume intersections between the reference mesh and the deformed mesh. But this procedure requires the computation of all the geometrical elements. Moreover, the projection needs local tests of projection with conditional branching (think about the very different cases of a compression and expansion). Thus the procedure is not SIMD and potentially leads to a loss of performance. The incremental remapping 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:

1. Backward convection:

   (4)

2. Forward convection:

   (5)

Each convection problem is well-posed on the time interval [0, Δtⁿ]. Let us now focus into these two steps and the way to solve them.

3.1 Backward Convection in Lagrangian Description

After the Lagrange step, if we solve the backward convection problem (4) over a time interval Δt ⁿ using a Lagrangian description, we have(6)

Actually, from the cell K^n+1,L we go back to the original cell K with conservation of the conservative quantities. For (conservation of mass), we haveshowing the variation of density by volume variation. For , it is easy to see that both velocity and specific total energy are kept unchanged in this step:

Thus, this step is clearly computationally inexpensive.

3.2 Forward Convection in Eulerian Description

From the discrete field defined on the Eulerian cells K, we then solve the forward convection problem over a time step Δtⁿ under an Eulerian description. A standard finite volume discretization of the problem will lead to the classical time advance scheme(7)for some interface values defined from the local neighbor values . We finally get the expected Eulerian values at time tⁿ⁺¹.

Notice that from (6) and (7) we have also(8)thus completely defining the remap step under the finite volume scheme form (8). One can notice that neither mesh intersections nor geometric considerations are required to achieve this 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).

3.3 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 :(9)that can also be written(10)

We recognize in (10) pressure-related fluxes and convective fluxes that define the whole numerical flux.

3.4 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 Lagrangian velocity vector field 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 . 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 Lagrange-remap 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.

4 Derivation of a Second-Order Accurate Lagrange-Flux Scheme

From the intermediate conclusions of the discussion 3.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 at time .

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 fulfills the usual stability CFL condition. We have

By making t tend to zero, (t > 0), we have A^L(t/2) → A, , , , then we get a semi-discretization in space of the conservation laws. That can be seen as a particular method of lines [13]:(11)

We get a classical finite volume methodwith a numerical flux Φ _A whose components are(12)

In (11), pressure fluxes and interface normal velocities can be computed from an approximate Riemann solver in Lagrangian coordinates (for example the Lagrangian HLL solver, [14]). Then, the interface states should be computed from an upwind process according to the sign of the normal velocity . 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 [16]:(13)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 (Runge-Kutta 2 (RK2), etc.). For the second-order Heun scheme for example, we have the following algorithm:

1. Compute the time step Δt ⁿ subject to some CFL condition;

2. Predictor step. MUSCL reconstruction + slope limitation on primitive variables ρ, u and p. From the discrete states , compute a discrete gradient for each cell K and interpolated values at cell interfaces;

3. Use a Lagrangian approximate Riemann solver to compute pressure fluxes and interface velocities from the MUSCL reconstructed values at the interface;

4. Compute the upwind edge values according to the sign of ;

5. Compute the numerical flux as defined in (12);

6. Compute the first order predicted states :

   

7. Corrector step. MUSCL reconstruction + slope limitation: from the discrete values , compute a discrete gradient for each cell K and interpolated values at cell interfaces;

8. Use a Lagrangian approximate Riemann solver to compute pressure fluxes and interface velocities from the MUSCL reconstructed values at the interface;

9. Compute the upwind edge values according to the sign of ;

10. Compute the numerical flux as defined in (12);

11. Compute the second-order accurate states at time t ⁿ⁺¹:

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 [4]. Here the predictor and corrector kernel functions have similar programming codes and there is no intermediate variables to save in memory.

4.1 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(14)and the normal contact velocity u ^⋆ by(15)leading to simple formulas easily implementable in the second-order Heun time integration scheme.

4.2 Numerical Experiments

As example, we test the Lagrange-flux scheme presented in Section 4 on few one-dimensional shock tube problems. We use a RK2 time integrator and a MUSCL reconstruction with the Sweby slope limiter given in (13).

4.2.1 Sod’s Shock Tube [17]

The initial data defined on space interval [0, 1] is made of two constant states and with initial discontinuity at x = 0.5. We successively test the method on two uniform mesh grids made of 100 and 400 cells, respectively. The final computational time is T = 0.23 and we use a CFL number equal to 0.25 and a Sweby limiter with coefficient β = 1.5. On Figure 2, one can observe a nice behavior of the Euler solver, with sharp discontinuities and a low numerical diffusion into the rarefaction fan even for the coarse grid.

4.2.2 Two-Rarefaction Shock Tube

The second reference example is a case of two moving-away rarefaction fans under near-vacuum conditions [14]. It is known that the Roe scheme breaks down for this case. The related Riemann problem is made of the left state and right state . The final time of T = 0.16. We again test the method on a coarse mesh (200 points) and a fine mesh (2000 points). Numerical results are given in Figure 3. The numerical scheme appears to be robust especially in near-vacuum zones where both density and pressure are close to zero.

4.2.3 Case With Sonic Rarefaction and Supersonic Contact

The following shock tube case with initial data and generates a sonic 1-rarefaction, a supersonic 2-contact discontinuity and a 3-shock wave. The final time is T = 0.16 and we use 400 mesh points, CFL = 0.25. Numerical results show a good capture of the rarefaction wave, without any non-entropic expansion-shock (Fig. 4).

4.2.4 Case of Shock-Shock Hypersonic Shock Tube

This last shock tube problem is a violent flow case made of two hitting fluids with and . Both 1-wave and 3-wave are shock waves, and the right state has a Mach number of order 40. Final time is T = 0.16, we use 400 grid points and the limiter coefficient β is here 1 (equivalent to the minmod limiter). One can observe a nice behavior of the solver: there is no pressure or velocity oscillations at the contact discontinuity, and the numerical scheme preserves the positivity of density, pressure and internal energy (Fig. 5).

5 Performance Results

In this section, we compare the new Lagrange-flux scheme to the reference (staggered) Lagrange-remap 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 [4], 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 fluxes are computed and variables are updated for the prediction. The second kernel, named 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 Table 1 a performance comparison study. First, we notice that the baseline performance – i.e. the single core absolute performance without vectorization – is quite similar for the two schemes, as can be seen in the first column. However, we the Lagrange-flux scheme has a better scalability, due to both vectorization and multithreading: our Lagrange-flux implementation achieves a speedup of 31.1X with 16 cores and AVX vectorization (a 2X speedup from AVX vectorization, which is ideal for this solver with many divisions instructions, and an almost 16X perfect speedup from the multicore usage) whereas the reference Lagrange-remap algorithm reaches a speed-up of only 14.8X. This difference is mainly due to the memory-bound kernels composing the reference Lagrange-remap scheme. Indeed, speedups due to AVX vectorization and multithreading are not ideal for kernels with relatively low intensity since memory bandwidth is shared between cores.

6 Dealing With Multimaterial Flows

Although this is not the aim and the scope of the present paper, we would like to give an outline of the possible extension of Lagrange-flux schemes to compressible multimaterial/multifluid flows, i.e. flows that are composed of different immiscible fluids and separated by free boundaries.

For pure Lagrange+remap schemes, usually VOF-based Interface Reconstruction (IR) algorithms are used (Young’s PLIC, etc.). After the Lagrangian evolution, for the cells that host more than one fluid, fluid interfaces are reconstructed. During the remapping step, one has to evaluate the mass fluxes per material. From the computational point of view and computing performance, this process generally slows down the whole performance because of many array indirections in memory and specific treatment into mixed cells along with the material interfaces.

If the geometry of the Lagrangian cells is not completely known (as in the case of Lagrange-flux schemes), 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 next two 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.

6.1 Interface Capturing for Pure Advection Problems

Let us first present numerical results on a pure scalar linear advection problem. The forward-backward advection case proposed by Rider and Kothe [22] is a hat-shaped function which is advected and stretched into a rotating velocity 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² 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.

6.2 Three-Material Hydrodynamics Problem

We then consider our interface capturing method for multifluid hydrodynamics. Because material mass fractions are advected, that isone can use the advection solver of these variables but we prefer dealing with the conservative form of the equationsin order to enforce mass conservation [23]. It is known that Eulerian interface-capturing 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 as described in Figure 7. The domain is split 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.5 and γ₂ = 1.4. Due to the density differences, two shocks in sub-domains Ω₂ and Ω₃ propagate with different speeds. This creates 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.8 M 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 Lagrange-remap 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.

References

All Tables

All Figures