A high performance and portable all-Mach regime flow solver code with well-balanced gravity. Application to compressible convection

TL;DR

This paper introduces ARK, a high-performance, portable flow solver for compressible convection that accurately handles low and high Mach regimes, with applications in astrophysics such as stellar atmospheres and accretion shocks.

Contribution

The paper presents a new numerical solver, ARK, which combines a well-balanced gravity scheme and low Mach correction for accurate, high-performance simulations of stratified compressible flows.

Findings

Achieves machine precision hydrostatic balance.

Improves accuracy in low Mach regimes over standard methods.

Enables detailed 3D turbulent convection simulations.

Abstract

Convection is an important physical process in astrophysics well-studied using numerical simulations under the Boussinesq and/or anelastic approximations. However these approaches reach their limits when compressible effects are important in the high Mach flow regime, e.g. in stellar atmospheres or in the presence of accretion shocks. In order to tackle these issues, we propose a new high performance and portable code, called "ARK" with a numerical solver well-suited for the stratified compressible Navier-Stokes equations. We take a finite volume approach with machine precision conservation of mass, transverse momentum and total energy. Based on previous works in applied mathematics we propose the use of a low Mach correction to achieve a good precision in both low and high Mach regimes. The gravity source term is discretized using a well-balanced scheme in order to reach machine precision hydrostatic balance. This new solver is implemented using the Kokkos library in order to achieve high performance computing and portability across different architectures (e.g. multi-core, many-core, and GP-GPU). We show that the low-Mach correction allows to reach the low-Mach regime with a much better accuracy than a standard Godunov-type approach. The combined well-balanced property and the low-Mach correction allowed us to trigger Rayleigh-B\'enard convective modes close to the critical Rayleigh number. Furthermore we present 3D turbulent Rayleigh-B\'enard convection with low diffusion using the low-Mach correction leading to a higher kinetic energy power spectrum. These results are very promising for future studies of high Mach and highly stratified convective problems in astrophysics.

Figures

Tables

Table 1. Table 1: Isothermal atmospheres at rest. Table shows for different spatial resolutions the maximum velocity in the domain. We see that the velocity is maintained around zero up to the machine precision , thus illustrating the well-balanced property.

Number of cells	velocity ($\theta =1$)	velocity ($\theta =Ma$)
128	${2.910}^{-15}$	${1.410}^{-13}$
256	${8.110}^{-15}$	${5.710}^{-13}$
512	${1.510}^{-14}$	${1.110}^{-12}$
1024	${2.210}^{-14}$	${2.210}^{-12}$
2048	${4.710}^{-14}$	${1.610}^{-12}$
4096	${1.110}^{-13}$	${4.010}^{-12}$

Equations

$${\partial }_t\rho$$

$${\partial }_t\left(\rho \mathbf{u}\right)$$

$${\partial }_t\left(\rho E\right)$$

$${\mathbf{\tau }}_{\text{visc}}=\mu \left(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T\right)+\eta \left(\nabla \cdot \mathbf{u}\right)\mathbf{I},$$

$$p$$

$${\mathbf{q}}_{\text{heat}}$$

$$\eta$$

$${\partial }_t\left(\rho \Phi \right)+\nabla \cdot \left(\rho \Phi \mathbf{u}\right)=\rho \mathbf{u}\cdot \nabla \Phi .$$

$${\partial }_t\left(\rho \mathcal{E}\right)+\nabla \cdot \left(\rho \mathcal{E}\mathbf{u}-{\mathbf{\sigma }}_{\text{stress}}\mathbf{u}-{\mathbf{q}}_{\text{heat}}\right)=0$$

$$\nabla p=-\rho \nabla \Phi ,\mathbf{u}=0,$$

$${\partial }_t\rho$$

$${\partial }_t\left(\rho \mathbf{u}\right)$$

$${\partial }_t\left(\rho \mathcal{E}\right)$$

$${\partial }_t\left(\rho \Phi \right)$$

$$\rho \mathcal{E}=\rho e+\frac{1}{2}\rho {\mathbf{u}}^2+\rho \Phi .$$

$${\partial }_t\rho$$

$${\partial }_t\left(\rho \mathbf{u}\right)$$

$${\partial }_t\left(\rho \mathcal{E}\right)$$

$${\partial }_t\left(\rho \Psi \right)$$

$$\rho \mathcal{E}=\rho e+\frac{1}{2}\rho {\mathbf{u}}^2+\rho \Psi .$$

$${\partial }_t\rho$$

$${\partial }_t\left(\rho \mathbf{u}\right)$$

$${\partial }_t\left(\rho \mathcal{E}\right)$$

$${\partial }_t\left(\rho \Psi \right)$$

$${\partial }_t\rho$$

$${\partial }_t\left(\rho \mathbf{u}\right)$$

$${\partial }_t\left(\rho \mathcal{E}\right)$$

$${\partial }_t\left(\rho \Psi \right)$$

$${\partial }_t\tau$$

$${\partial }_{t}u$$

$${\partial }_{t}v$$

$${\partial }_t\mathcal{E}$$

$${\partial }_t\Pi$$

$${\partial }_t\Psi$$

$$\mathcal{E}=e+\frac{1}{2}(u^2+v^2)+\Psi .$$

$${\tilde{\tau }}_i$$

$${\tilde{u}}_i$$

$${\tilde{v}}_i$$

$${\tilde{\mathcal{E}}}_i$$

$${\tilde{\Pi }}_i$$

$${\tilde{\Psi }}_i$$

$$u_{i+1/2}^{\ast }$$

$${\Pi }_{i+1/2}^{\ast }$$

$$a_{i+1/2}^n$$

Full text

In order to tackle these issues, we propose a new high performance and portable code, called “ARK” with a numerical solver well-suited for the stratified compressible Navier-Stokes equations. We take a finite volume approach with machine precision conservation of mass, transverse momentum and total energy. Based on previous works in applied mathematics we propose the use of a low Mach correction to achieve a good precision in both low and high Mach regimes. The gravity source term is discretized using a well-balanced scheme in order to reach machine precision hydrostatic balance. This new solver is implemented using the Kokkos library in order to achieve high performance computing and portability across different architectures (e.g. multi-core, many-core, and GP-GPU).

We show that the low-Mach correction allows to reach the low-Mach regime with a much better accuracy than a standard Godunov-type approach. The combined well-balanced property and the low-Mach correction allowed us to trigger Rayleigh-Bénard convective modes close to the critical Rayleigh number. Furthermore we present 3D turbulent Rayleigh-Bénard convection with low diffusion using the low-Mach correction leading to a higher kinetic energy power spectrum. These results are very promising for future studies of high Mach and highly stratified convective problems in astrophysics.

\newdateformat

myformat\monthname[\THEMONTH] \THEYEAR

Introduction

The study of convection is an active topic of research in the astrophysics community because of its major role in different mecanisms such as heat transport in solar and stellar interiors (Spruit et al., 1990), mixing of elements (Pinsonneault, 1997) and dynamo (Charbonneau, 2014). As these mecanisms play a role in the estimation of the lifetime of these objects it is of great importance for stellar evolution theory.

Different approximations have been developed to ease the study of convection. The Boussinesq and the anelastic approximations simplify the Navier-Stokes system by getting rid of acoustic waves and keeping buoyancy effects. In practice these approximations are derived by looking at the equations satisfied by small perturbations near a reference state (Spiegel & Veronis, 1960). The Boussinesq approximation is quite restrictive as it is valid for a small layer of the reference state, such that the flow can be considered incompressible. On the other hand the anelastic approach allows to have a larger scale height by keeping the density stratification of the reference state (Gilman & Glatzmaier, 1981). Another way to understand these approximations is to consider the flow regime in terms of the Mach number $\mathrm{Ma}$. As it is shown in Mentrelli (2018), these approximations can be recovered by considering low-Mach asymptotic limits of the Navier-Stokes system. The Froude number, defined as the non-dimensional ratio of kinetic energy to gravitational energy, characterizes the influence of gravity in the flow. By taking into account different Froude regimes, they recover the incompressible, the Boussinesq and the anelastic models. From a numerical point of view the removal of the acoustics waves in these models is quite attractive because it allows to have larger time steps. The anelastic model has been successfully implemented in different codes like Rayleigh (Featherstone & Hindman, 2016) or Magic (Gastine & Wicht, 2012) and it is widely used in the community (see Glatzmaier, 2017). We can also mention the MAESTRO code (see Nonaka et al., 2010) which uses an extended version of the anelastic model. The velocity constraint takes into account the time variation of pressure. However these approaches present some drawbacks. The addition of new physics and source terms to the model is difficult, one has to derive another asymptotic model to take the new physics into account in the anelastic regime (see Mentrelli, 2018). Furthermore one has to be careful that the simulation stays in the regime of validity of the model (especially in the Boussinesq regime). Finally a numerical difficulty is the parallelization of those codes. They usually use pseudo-spectral methods for which it is more difficult to achieve a good scalability (e.g. need to use pencil-type domain decomposition Featherstone & Hindman, 2016).

We chose to take a more flexible approach by solving the full compressible Navier-Stokes system, as in the MUSIC code (Viallet et al., 2011; Goffrey et al., 2017) but with a collocated finite volume solver instead of using a staggered grid. Different discretization techniques of the Euler system are used in the astrophysics community. We can classify them in various ways. One way is to separate SPH techniques from grid-based techniques. Furthermore grid-based approaches can be divided in different families, finite difference, finite element and finite volume. The finite volume method is of particular interest because of its natural property of being conservative and to capture shocks and discontinuities. Designing a finite volume scheme essentially resides in the definition of a numerical flux, numerical counterpart of the physical flux. A widely used family of fluxes is the Godunov (see Godunov, 1959) flux which is the flux of the — usually approximate — Riemann problem between two neighbour cells.

However we have to face multiple numerical difficulties with this approach. Compressible solvers and mainly Godunov-type solvers are known to have an excessive amount of numerical diffusion in the low-Mach regime which make them unusable in this regime (see Guillard & Viozat, 1999; Dellacherie, 2010; Miczek et al., 2015; Chalons et al., 2016; Barsukow et al., 2017). In this regime, in which flows are smoother, considering Riemann problems at interfaces is not adapted. Indeed in the work of Miczek et al. (2015) they show that part of the kinetic energy is dissipated into internal energy whereas it should be conserved. To tackle this issue they propose a preconditionned Roe scheme to remove the numerical diffusion. Secondly, hydrodynamics and gravity are usually discretized independently from each other. In the case of highly stratified medium, the numerical scheme does not maintain the hydrostatic equilibrium and produces spurious flows that pollutes the simulation. Different approaches have been investigated to solve this issue both for the Euler and the shallow water equations. In Leroux & Cargo (1994), they rewrite the Euler system as a fully conservative system by defining an hydrostatic pressure satisfying a conservation law. In Chandrashekar & Klingenberg (2015) they use a variable reconstruction by taking advantage of the equilibrium profile. In Chalons et al. (2010); Vides et al. (2014); Chalons et al. (2017) they incorporate the source term in the Riemann problem itself allowing to compensate pressure gradients at the interface. As in Leroux & Cargo (1994), authors from Chertock et al. (2018) also propose to discretize the Euler system with gravity as a fully conservative system but using global fluxes and a reconstruction on equilibrium variables. Finally the last numerical difficulty is the time step in the low-Mach regime. Because of the stability condition involving the fast acoustic waves, the time step becomes very small compared to the material transport timescale. It can either be resolved using a full implicit approach as in the MUSIC code Viallet et al. (2011); Goffrey et al. (2017), or by using an implicit-explicit (IM-EX) approach in which only the system with fast acoustic waves is solved implicitly (Chalons et al., 2016, 2017).

Following the original work of Chalons et al. (2016) and Chalons et al. (2017) we use an acoustic-transport splitting. In Chalons et al. (2016) they derive a finite volume scheme of the Euler system on unstructured mesh. This scheme uses an acoustic splitting to separate acoustic waves from material ones. In the low Mach regime, this translates to a splitting between fast waves and slow waves. In the low Mach regime, the fast waves can be treated with an implicit solver to get rid of the restrictive stability condition. Then in the work of Chalons et al. (2017), the scheme has been adapted to shallow water equations with a source term which is the topography. This source term is added in the equivalent acoustic subsystem to obtain a well-balanced scheme. In this paper we adapt their approach for the Euler system by taking care of the discretization of the energy equation.

The paper is organized as follows. In Section 1 we briefly recall the compressible model we use to study convection, i.e. the Navier-Stokes equations with gravity. In Section 2 we present the derivation of the well-balanced and all-regime numerical scheme using a splitting approach between an acoustic step and a transport step both solved explicitly in this work. In Section 3 we present some implementation features about the “ARK” 111https://gitlab.erc-atmo.eu/erc-atmo/ark, version v1.0.0 code in particular the Kokkos library used for the shared memory parallelization. We also give some performance results. Finally in Section 4 we present different numerical test cases illustrating the importance of the low-Mach correction and the well-balanced discretization of gravity.

1 Navier-Stokes equations

We want to solve Navier-Stokes equations expressing conservation of mass, balance of momentum and balance of energy, respectively written as follows

[TABLE]

where $\rho$ is the density, $\mathbf{u}$ the material velocity, $p$ the pressure, $\mathbf{g}$ the external gravitational field, $\rho E=\rho e+\frac{1}{2}\rho\mathbf{u}^{2}$ the density of total energy with $e$ the specific internal energy, $\mathbf{q}_{\text{heat}}$ the heat flux and $\bm{\tau}_{\text{visc}}$ the viscous tensor satisfying

[TABLE]

where $\mu$ is the dynamic viscosity and $\eta$ the bulk viscosity. We use $\cdot$ as a scalar product and thus $\bm{\nabla\cdot}{}$ represents the divergence operator. In order to close Navier-Stokes system (1) we add constitutive equations namely a pressure law $p^{\text{EOS}}$ (3a), the Fourier’s law (3b) and the Stokes hypothesis (3c)

[TABLE]

We recall that the gravitational field is derived from a gravitational potential $\Phi$ for which $\mathbf{g}=-\bm{\nabla}{\Phi}$. Dealing with a constant in time external gravity field, $\partial_{t}\Phi=0$ and using the conservation of mass we get (4)

[TABLE]

Let us emphasize that in this equation, the gravitational energy $\rho(\mathbf{x},t)\Phi(\mathbf{x})$ is time dependent only through the density $\rho(\mathbf{x},t)$. Hence the energy equation (4) can be rewritten in the following conservative form

[TABLE]

where we define $\rho\mathcal{E}=\rho e+\frac{1}{2}\rho\mathbf{u}^{2}+\rho\Phi$. Equation (5) expresses the local conversion between three different energy reservoirs, as depicted in figure 1: internal, kinetic and gravitational.

There can be a direct transfer between gravitational energy and kinetic energy through the work of gravitational forces, from kinetic energy to internal energy because of the second law of thermodynamics but no direct transfer between gravitational energy and internal energy, see also Section 5 of Springel (2010) and Section 2.2 of Marcello & Tohline (2012) for a discussion on energy conservation for both external and self-gravity cases.

Because of this conservation of energy including gravitational energy we will use the formulation (5) of the energy equation and we will use the gravitational potential instead of the usual gravitational field $\mathbf{g}$. To our knowledge this approach is quite rare, see Graham (1975) or Chertock et al. (2018) where they use global fluxes to have a well-balanced and conservative scheme.

An important steady state solution of this system for stratified objects is the hydrostatic balance. The flow is static and the gravitational force is balanced by the pressure forces, i.e. following equation (6)

[TABLE]

As we mentioned in the introduction, convective flows can be considered as a perturbation flow of the hydrostatic equilibrium. Thus this steady state is particularly important in order to study convection problems in stratified flows.

2 Numerical scheme

2.1 Euler system — Hyperbolic system

Before going into the derivation of the scheme we introduce the notations. We define by $\Delta x$ (resp. $\Delta y$ and $\Delta z$) the step along the x-direction (resp. the y and z-direction). We note by $\Delta t$ the time interval between current time $t^{n}$ and $t^{n+1}$. We use the notation $q^{n}_{i}$ (resp. $q^{n}_{i,j,k}$) to represent the averaged quantity associated to the field $q$ at time $t^{n}$ and in the cell $i$ (resp. $i,j,k$) in the one-dimensional case (resp. the three-dimensional case). We use the notation $q^{n}_{i+1/2}$ (resp. $q^{n}_{i+1/2,j,k}$) to represent the quantity associated to the field $q$ at time $t^{n}$ and at the interface between cells $i$ and $i+1$ (resp. $i,j,k$ and $i+1,j,k$) in the one-dimensional case (resp. the three-dimensional case). Finally we define the notation ${[q]}_{i}=q_{i+1/2}-q_{i-1/2}$ in the one-dimensional case.

2.1.1 Acoustic-Transport splitting approach

Following Chalons et al. (2017) we use a splitting strategy that separates acoustic terms and transport terms and we choose to add the gravitational source terms to the acoustic part. This way, pressure gradient can be balanced by the gravity source term.

However we have another equation compared to the shallow water system that is the energy equation. As in Chalons et al. (2017), we want an isentropic acoustic step for smooth solutions. Thereby we choose to solve the equation on the gravitational energy,

[TABLE]

[TABLE]

However, this leads to a non constant gravitational potential in the acoustic step whose time variations are compensated in the transport step in order to have a constant potential in the full step. The potential is constant in the full step at the continuous level, but discretization errors with the splitting can lead to a non-constant discretized potential. Thus we choose to introduce an approximation of the gravitational called $\Psi\approx\Phi$ and a relaxation parameter $\lambda$

[TABLE]

[TABLE]

We consider the relaxation system (8) to be an approximation of the original system (7) that we formally recover in the limit $\lambda\to 0$. System (8) is solved by first solving the system in the limit $\lambda\to\infty$ and then in the limit $\lambda\to 0$ in which $\Psi$ is projected onto $\Phi$, the initial condition. This way, the evolution of the gravitational potential $\Psi$, consistent with zero, is forced to be constant. The relaxation technic used here for the gravitational potential is similar to what is done for pressure relaxation in many approximate Riemann solvers and we emphasize that $\Psi$ is just an intermediate used to design the scheme and can be removed when writing the final scheme (see 2.1.5) .

We now turn to the discretization of the system (8) in the limit $\lambda\to\infty$. Transport phenomena of the form $\mathbf{u}\cdot\bm{\nabla}{}$ are separated from the other terms to give two subsystems, first the acoustic subsystem (9)

[TABLE]

then the transport subsystem (10)

[TABLE]

We now briefly study the eigenstructure of systems (9)-(10). Let $\mathbf{n}$ be any unit normal vector, the acoustic system (9) involves seven eigenvalues: $-c,0,c$. The fields associated with 0 (resp. $\pm c$) are linearly degenerate (resp. genuinely nonlinear), see Appendix A for more details. The eigenvalues for transport system (10) are given by $\mathbf{u}\cdot\mathbf{n}$. Both systems (9)-(10) are hyperbolic. We emphasize here that the choice of using a relaxation procedure for the gravitational potential by introducing the equation on the gravitational potential energy $\rho\Psi$ has been made to obtain this simple wave pattern for the splitted Euler system with gravity. (i.e. the same pattern as without gravity). Other choices for the relaxation procedure (e.g. $\partial_{t}\Psi=0$ in both steps) would either lead to the introduction of $\mathbf{u}\cdot\mathbf{n}$ in the eigenvalues of the acoustic subsystem or would significantly complexify the relaxation procedure for the pressure..

To summarize our numerical procedure, we propose to define a flux interface by approximating system (7) with a three-step procedure that involves solving the acoustic system (9) (acoustic step), the transport system (10) (acoustic step) and finally project $\Psi$ onto $\Phi$ (relaxation step). We detail each step in the next sections using the one-dimensional equations.

2.1.2 Acoustic step

We follow the idea of Chalons et al. (2016) to discretize the acoustic subsystem. They introduce a pressure relaxation $\Pi\approx p$, an acoustic impedance $a\approx\rho c$ and a relaxation parameter $\nu$ to get a fully linearly degenerated system. It is then written using Lagrangian variables $(\tau,u,v,\mathcal{E},\Psi)$ where $u$ represents the normal velocity component at an interface and $v$ a transverse component. We also use a mass variable $\mathrm{d}m=\rho(t^{n},x)\mathrm{d}x$ where time is frozen at instant $t^{n}$

[TABLE]

where

[TABLE]

The discretization of this system is realized with an approximate Riemann solver that accounts for the source term by means of integral consistency and composed by three waves $-a,0,a$, see Gallice (2002); Chalons et al. (2013, 2017). After the relaxation, in which $\nu\to 0$, it gives

[TABLE]

where

[TABLE]

and $a^{n}_{i+1/2}\geq\max{\left(\rho^{n}_{i}c^{n}_{i},\rho^{n}_{i+1}c^{n}_{i+1}\right)}$ which is the so-called sub-characteristic condition (see Chalons et al., 2013).

The update of the conservative variables is then

[TABLE]

where $\widetilde{L}_{i}=1+\frac{\Delta t}{\Delta x}{\left[u^{*}\right]}_{i}$.

2.1.3 Transport step

The transport subsystem can be written in the following form, for $b\in\left\{\rho,\rho u,\rho v,\rho\mathcal{E},\rho\Psi\right\}$

[TABLE]

that is discretized as follows

[TABLE]

The interface term $\widetilde{b}_{i+1/2}$ is defined by the upwind choice with respect to the velocity $u^{*}_{i+1/2}$

[TABLE]

2.1.4 Relaxation step

At this stage, the relaxed gravitational potential $\Psi$ still evolves in time. So we perform the relaxation $\lambda\to 0$ that boils down to set $\Psi^{n+1}_{i}=\Phi_{i}$.

2.1.5 Overall algorithm

Gathering the previous steps and intermediate variables, the overall scheme reads

[TABLE]

It may also be expressed as a first-order classic finite-volume scheme involving flux terms for the conservative part for energy $\rho E=\rho e+\frac{1}{2}\rho u^{2}$ and source terms for gravity

[TABLE]

where

[TABLE]

We emphasize that both formulations are equivalent and conservative for the energy $\rho\mathcal{E}$. A non-conservative energy approach is also detailed in Appendix B.

We can notice that in the case of a constant gravitational potential, we recover the original scheme derived in Chalons et al. (2016).

2.1.6 On the low-Mach correction

As for the scheme of Chalons et al. (2016) and as explained in Dellacherie (2010), the numerical scheme defined by (11) poorly performs in the low Mach regime due to truncature error of magnitude $\frac{\Delta x}{\mathrm{Ma}}$ that comes from the term $\Pi^{*}_{i+1/2}$. To tackle this issue, following Chalons et al. (2016) we modify the upwinding part of $\Pi^{*}_{i+1/2}$ thanks to an extra parameter $\theta_{i+1/2}$ by setting

[TABLE]

[TABLE]

Using a truncation analysis in dimensionless form it can be shown that this correction acts like a rescaling of the numerical diffusion induced by the pressure discretization (see Chalons et al., 2016).

As we can see, the low-Mach correction does not directly come from the derivation of the numerical scheme 11. Some ongoing works are trying to derive directly all-Mach schemes using more sophisticated relaxation schemes (see Bouchut et al., 2017).

2.1.7 On the well-balanced property

A numerical scheme is said to be well-balanced for equilibrium states satisfying equation (6), if it exists a discrete counterpart of equation (6) in which solutions are preserved by the numerical scheme.

The discrete counterpart of equation (6) for scheme (11) is given by

[TABLE]

Let us now verify that we have obtained a well-balanced scheme. If at time $t^{n}$, for some density profile the initial state reads as in (16) then fluxes from the acoustic step reduce to

[TABLE]

Then we have for the acoustic step

[TABLE]

Finally, because $u^{*}_{i+1/2}$ vanishes, transport step is trivial and the initial state remains unchanged. Once we have made the appropriate choice for the discretization of the gravitational source term in the acoustic step, the well-balanced property is automatically verified without the need to introduce an other algorithmic correction.

2.2 Dissipative fluxes — Parabolic

system

We now turn to the discretization of dissipative fluxes (2)-(3b). They are discretized using first order discrete fluxes

[TABLE]

where $\mathbf{f}^{dissipative}$ is either the heat flux $\mathbf{q}_{\text{heat}}$ or the viscous flux $\bm{\tau}_{\text{visc}}$. In the case of the heat flux we have

[TABLE]

With the addition of the viscous terms and the heat flux, this all-regime well-balanced scheme is now well-suited for the study of convection problems in highly stratified flows in both low Mach and high Mach regimes. Before showing validating numerical tests, we present some specificities about the numerical implementation and parallelization used in this work.

3 Implementation and parallelization

In this section we describe the implementation of the scheme using Kokkos library. We begin by giving a brief overview of the Kokkos library.

3.1 Exascale computing

To reach the exascale, the distributed memory model is not sufficient to take advantage of all the computing power of new architectures. There are mainly two reasons for this. First, nodes of supercomputers tend to grow more and more and hence are more suited to a shared memory model (Sunderland et al., 2016). Secondly, nodes tend be more and more heterogeneous by using multi-core, many-core and/or accelerators like GP-GPUs. So it means that even if shared memory is exposed, it needs to be handled differently from one architecture to another. For example we can think of OpenMP or C++11 threads for multi-core and many-core processors, and CUDA or OpenACC for GP-GPUs.

Moreover this architecture heterogeneity raises a performance portability issue. Currently, many HPC codes are optimized for some specific architectures to get the maximum computing power. However this optimization process couples the numerical scheme to its implementation details like the memory management, the loop ordering, cache blocking and so on. Hence running a code on a different architecture results in bad performance.

We propose to use the recent C++ library Kokkos (see Carter Edwards et al., 2014) that implements a new shared memory model. Using abstract concepts such as execution spaces (where a function is executed), data spaces (where data resides) and execution policies (how the function is executed) the library is able to efficiently take advantage of multi-core many-core processors and GP-GPUs. This way the portability relies on the library and no more on the numerical code.

3.2 Implementation

Following the work of Kestener (2017), the code is then organized with computation kernels:

• •

Acoustic and transport kernels,

• •

Viscous and heat diffusion operator kernels,

• •

Conservative variables to primitive variables kernel,

• •

Time step kernel.

Each kernel is a C++ functor. They are given to Kokkos through the function Kokkos::parallel_for. Internally, depending on the device chosen at compile-time, it hides a parallelized one-dimensional loop where the current index is given as an argument to the functor. This index is then interpreted as a cell index in the domain.

Kokkos only deals with shared memory systems. We use the Message Passing Interface (MPI) programming model with a regular domain decomposition to take advantage of distributed memory machines across multiples computing nodes. Kokkos is then used as a shared memory programming model inside each node. These domains are endowed with ghost zones which are used to both implement physical boundary conditions and to contain values from neighbour domains. Communications are handled through the ghost cell pattern (see Kjolstad & Snir, 2010). Thus for a given direction X, Y (or Z) and a given side, left or right, one MPI process sends data from its domain to its neighbour’s ghost zone and receives data into its own ghost zone.

3.3 Performance results

Thanks to Kokkos, we were able to use the same code on different architectures like Intel Skylake, Intel Knights Landing (KNL) and NVIDIA GP-GPUs (K80, P100, V100). We measured performance on the Intel Skylake and the Intel KNL partition of the Joliot Curie machine at TGCC. Figure 2 shows the results. We see that the Kokkos library is able to provide good performance on the different tested architectures. Nevertheless, even if the peak performance of the Intel KNL architecture is higher than the Intel Skylake one we have better performance on the Intel Skylake architecture. We also notice the important speed-up (around five) between the Intel Skylake architecture and the NVIDIA V100 GP-GPU.

Figure 3 shows a weak scaling test performed with a hybrid configuration OpenMP/MPI. We went up to 512 MPI processes, one MPI process per Intel Skylake socket to avoid NUMA effects. It results in a total of 12288 cores at 512 MPI processes. Each MPI process is getting a piece of the whole domain of $128^{3}$, so a domain of $44^{3}$ per core. We can see that we obtain a plateau of 85% of maximum performance from 128 MPI processes.

The performances obtained with the use of the Kokkos library are encouraging for the study of convection problems with the ARK code on massively parallel present and future architectures. In the next section, we use several numerical tests to show that the numerical scheme used in the ARK code is indeed very well suited for the study of convection.

4 Numerical results

In this section we specialize the equation of state 3a. We will use an ideal gas satisfying

[TABLE]

where $\gamma$ is the adiabatic index of the gas. The speed of sound satifies the following simple relation

[TABLE]

We emphasize that it is possible to use a different equation of state with the all-regime well-balanced numerical scheme. Moreover we consider two versions of the all-regime scheme depending on the low-Mach correction. We will refer to the disabled low-Mach correction scheme when $\theta=1$ and to the enabled one when $\theta$ follows equation 13.

We will test different properties of the scheme with different test cases: wave speeds with the Sod test (no gravity), low-Mach accuracy with the Gresho vortex test (no gravity), hydrostatic balance with the test of an atmosphere at rest and out of equilibrium behavior with the Rayleigh-Taylor test. We then use the ARK code for the study of Rayleigh-Bénard convection.

4.1 Shock tube test

The Sod shock tube (Sod, 1978) is a classical test for compressible solvers. It tests the ability of the solver to have correct wave speeds and its numerical diffusion near discontinuities.

The computational domain is the interval [0,1], the initial condition is defined by

[TABLE]

Results are shown in figure 4 for simulations with nx = 100.

First we can observe that the solver is as good as a first order Godunov-type scheme with a HLLC approximate Riemann solver around the contact discontinuity and the shock. However the rarefaction wave is a bit more diffused. We also notice that the low-Mach correction does not influence the behavior of the scheme for this test case. However we want to stress out some instability near discontinuities, as shown in Chalons et al. (2016). This can also be seen in a double shock waves test case.

4.2 Gresho vortex test case

The Gresho vortex (Gresho & Chan (1990); Miczek et al. (2015)) is a test case that has already been used to test numerical schemes in the low Mach regime. It is a two dimensional stationary test case that can be parameterized by the maximum value of the Mach number. It is thus well-suited to study the behavior of the scheme in the low Mach regime. We recall that the test case is defined using polar coordinates $(r,\theta)$ defined with respect to the center of the vortex as follows

[TABLE]

where $p_{0}$ satisfies $p_{0}=\frac{1}{\gamma\mathrm{Ma}^{2}}$. In this case $\mathrm{Ma}$ is a parameter and $\gamma$ is the adiabatic index of the ideal gas. The velocity is normalized so a particle placed at the peak of velocity ($u=1.0$ at location $r=0.2$) make a full rotation in $\Delta t=\frac{2}{5}\pi\approx 1.26$.

We ran a serie of simulations with different solvers where we explored parameter space nx and Ma from 32 to 2048 and from 1.0 to $1.0\times 10^{-5}$ respectively. Final time is set to $t_{f}=1.0\times 10^{-3}$, which has been chosen sufficiently small such that the error doesn’t saturate.

Figure 5 shows snapshots of the the magnitude of the velocity field at the final time and at resolution $512^{2}$. We see that when the Mach number decreases the velocity field becomes more and more degraded when the low-Mach correction is disabled. At $\mathrm{Ma}=10^{-5}$, the vortex has completely disappeared. Figures 6 and 7 show more quantitative results where we show absolute $L^{1}$ error on velocity in function of the Mach number Ma and the spatial resolution dx respectively. Figure 6 shows that $L^{1}$ error on velocity depends on the Mach number. More precisely we measure a slope of -1 on schemes or order 1 and a slope of -0.5 on scheme of order 2. On the other hand the low Mach correction of the all-regime scheme gives a uniform error with respect to the Mach number.

Figure 7 shows convergence curves at $\mathrm{Ma}=1.0\times 10^{-3}$. We see that both Godunov-type and all-regime without the low Mach correction converge at order 1 as expected. Nevertheless Godunov-type with Muscl-Hancock reconstruction converges only at order 1.5. It may be due to the lack of regularity of the velocity field as it can be observed in the case of a contact discontinuity (see Springel, 2010). All-Regime scheme shows two different behaviors, at first it converges at order 1.5 then around $\mathrm{nx}=1024$ the slope changes and it converges at order 1.2. We assume that at higher resolution we would recover order

1. We see that at low Mach number the precision, independently of the order, is better than the one of a Godunov-type scheme.

4.3 Well-balanced test case

The well-balanced test case is a simple isothermal column of atmosphere at equilibrium. This column of atmosphere is in a stable equilibrium state. The test allows us to measure the ability of the scheme to preserve this equilibrium. After normalization, it is given by

[TABLE]

which is the solution of the following system

[TABLE]

We take advantage of the formula (16) and we initialize the test case with the following formula

[TABLE]

The computational domain used is the interval $[0,3]$. Results are displayed in table 1 at time $t=10$, more than three times the sound crossing time in the box. We see that we stay near machine precision at the end of the simulation. We see a shift of two orders of magnitude in the error when using the low-Mach correction. The reason of this shift is not entirely clear and is difficult to interpret as it involves truncature errors. Looking at the spatial pattern of the error in the simulation, it does seem to come from the boundary conditions (extrapolation of the hydrostatic balance for pressure and density and reflexive conditions for the velocity) with the use of the low-Mach correction. A more appropriate boundary condition might remove this shift in the error (which is in any case sufficiently small and stable to allow the use of controlled seeded perturbations).

4.4 Rayleigh-Taylor instability test

case

The Rayleigh-Taylor test case is a two dimensional test case where two fluids of different densities are superposed and are at equilibrium. The denser one is on top. A small perturbation is introduced to break equilibrium.

The full setup is as follow, for a domain $\left[-0.25,0.25\right]\times\left[-0.75,0.75\right]$:

[TABLE]

Where $C=0.01$ is the magnitude of the velocity perturbation, $L_{x}=0.5$ and $L_{y}=1.5$ are the size of the domain in each direction. We do not need to use the well-balanced formula (16), the equilibrium is preserved in the case $A=0$.

Figure 8 shows two simulations of the Rayleigh-Taylor test case, one with the low Mach correction and the other without it ($\theta^{n}_{i+1/2}=1$). Both simulations are at the same time $t=12.4$ and the same resolution $200\times 600$. The yellow part is at density $2$ and the purple is at density $1$. We see that we recover the classical linear growing mode. Moreover the simulation with the low Mach correction is able to capture secondary instabilities in the non linear regime. They are closer to the second order Godunov-type simulation than the order one. However the low Mach correction does not help on the interface diffusion between the two mediums. It also shows a peak that is not present without the correction at the same resolution. This spurious behavior is therefore caused by the low Mach correction that removes some numerical diffusion in the scheme. By looking at higher resolutions, we identify that this peak is a grid-seeded secondary RT unstable mode that appears at the top of the large scale seeded mode. This type of secondary modes are not unexpected and can be seen for example in Fig. 9 of Almgren et al. (2010). This peak disappears with the addition of some physical viscosity in the simulation.

4.5 Rayleigh-Bénard instability test

case

This last test case is about compressible convection simulations both in 2D and 3D. In this test case there are different important parameters. First, from stability analysis we know that the Rayleigh number $\mathrm{Ra}$ is an important non-dimensional number. Beyond a threshold, called the critical Rayleigh $\mathrm{Ra}_{c}$, the convection process starts and efficiently transports the heat (see Figures 1 and 3 in Hurlburt et al., 1984). Below this threshold, diffusion processes are sufficient to transport heat and no material displacement is necessary. Then another important parameter is the density stratification $\chi$ which the ratio between the density at the bottom of the domain and the density at the top. In the highly stratified case, study of convection becomes more difficult as there is not a unique Rayleigh number but more a whole range of values extending on the scale height. Notice that when $\chi\to 1$ we recover the Boussinesq-like situation. Finally the last parameter is the polytropic index $m$ which is a measure of how close is the initial temperature gradient from the adiabatic gradient. One can show that the Schwarzschild criterion writes

[TABLE]

The initial setup is inspired from Hurlburt et al. (1984); Toomre et al. (1990). Following their notation, the initial state is given by a polytropic profile of polytropic index $m$

[TABLE]

where $z$ is the vertical variable. It is initialized using to the recursive formula (16). So we begin with a hydrostatic equilibrium that we destabilize whether with a velocity mode perturbation or with a temperature random perturbation.

4.5.1 2D case

We begin with 2D simulations in a weak stratification setup where $\chi=1.1$ and $m=1.3$ in order to be close to the adiabatic gradient. The initial perturbation is close to the fundamental velocity mode. The spatial resolution is set to $128^{2}$, and we impose the temperature flux on the bottom boundary. We then obtain stationary symmetric convective rolls. We study the effect of the low-Mach correction on the onset of the Rayleigh-Bénard instability by varying the initial Rayleigh number. Figure 9 shows the evolution of the mean absolute velocity. The linear phase, in logarithmic scale, corresponds to the exponential growth of modes. We can see that without the low-Mach correction we have an effective critical Rayleigh number between 10 and 15. Whereas with the low-Mach correction we recover an effective critical Rayleigh number close to the theoretical critical one.

If we now turn to a stronger stratification, the convective rolls pattern change. We increase the density ratio to $\chi=21$. Figure 10 shows a snapshot of the local Mach number field with the velocity field, low-Mach correction enabled. As observed in Hurlburt et al. (1984) we see a downward shift of the center of mass of convective rolls compared to the weak stratification case. By conservation of mass, the upper part of the convective roll has to be larger. The strong stratification case also exhibits higher Mach flows, around $\mathrm{Ma}\approx 0.5$ at the top of the box due to the low density. The all-regime well balanced scheme is indeed able to capture properly convection in highly stratified and high Mach flows.

4.5.2 3D case

We now turn to 3D simulations in a weak stratification situation. In this setup we want to look at the effect of the low-Mach correction on the kinetic energy spectrum in a more turbulent situation. So we change the polytropic index to $m=0.1$ and increase the initial Rayleigh number to $\mathrm{Ra}\approx 650000$. We also change the boundary condition to a fixed temperature for both boundaries in order to continuously force a large Rayleigh number in the simulation. By using different upscaling, from $128^{3}$ to $512^{3}$ we reach a stationary state 222The simulation outputs are available at http://opendata.erc-atmo.eu. Figure 12 shows a snapshot of the velocity in the box. We see large and structured vertical flows whereas in horizontal plans the flow is more turbulent. In order to study the different scales and the energy in this turbulent state we compute power spectrum of the kinetic energy of the horizontal middle plane. Figure 11 shows the results, the orange curve corresponds to the simulation performed with the low-Mach correction and the blue one without it. We see a net difference in the overall kinetic energy due to a lower dissipation into the internal energy. We notice that we recover higher kinetic energies at all scales showing that the low Mach correction is important to properly capture the power spectrum of turbulent convection.

Conclusion

We have presented a new numerical code that is able to perform simulations of convection without any approximation of Boussinesq nor anelastic type. To do so we have adapted an all-Mach number scheme into a well-balanced scheme for gravity. We have been able to show that it preserves arbitrary discrete equilibrium states up to the machine precision. Moreover the low-Mach correction in the numerical flux allows to be more precised in the low-Mach regime. This new scheme is well suited to properly study highly stratified and high Mach convective flows. The low Mach correction is important to properly capture convection modes in the laminar low Mach regime and the kinetic energy power spectrum in the turbulent regime. This code has been parallelized using a hybrid approach MPI+Kokkos in order to be well prepared for running on forthcoming exascale machines.

Further work will consist in using the implicit-explicit approach to reach very low Mach number simulations, see Chalons et al. (2016), and still keeping the well-balanced property for the gravity source term. Indeed by solving the acoustic part implicitly we avoid the restrictive CFL condition due to the fast acoustic waves. With both the explicit-explicit and implicit-explicit approach, this numerical scheme will be able to efficiently study convection problems in all regimes, low Mach and high Mach on the largest next generation massively parallel architectures.

Acknowledgement

P. Tremblin acknowledges supports by the European Research Council under Grant Agreement ATMO 757858. This work was granted access to the HPC resources of TGCC under the allocation A0040410097 attributed by GENCI (Grand Equipement National de Calcul Intensif). The authors acknowledge IDRIS (Institut du Développement et des Ressources en Informatique Scientifique) center to allow access to the Ouessant supercomputer. The authors would also like to thank Martial Mancip (CEA Saclay, Maison De La Simulation) for helping in remote visualization, Maxime Stauffert for his insight in the development of the well-balanced scheme and G. Grasseau (LLR, Polytechnique, IN2P3) to allow access to a NVIDIA V100 GPU.

Appendix A Eigenstructure of the acoustic system

For the sake of simplicity, the eigenstructure analysis of the acoustic system (9) is made in the one-dimensional case. We use the following change of variables, valid for smooth flows

[TABLE]

where $s$ is the specific entropy. By using equation of mass, one obtains

[TABLE]

By using the second law of Thermodynamics and the equation on the specific internal energy, one can show that $\partial_{t}s=0$ (see Godlewski & Raviart, 1996). Thus the acoustic system (9) writes equivalently

[TABLE]

The Jacobian matrix associated to the quasi-linear system A1 involves five eigenvalues: $-c<0<c$ where 0 is degenerated three times and $c$ satifies $c^{2}=\partial_{\rho}p^{\text{EOS}}\left(\rho,s\right)$. It is then hyperbolic. The four eigenvectors are given by

[TABLE]

Clearly the field associated to the stationary wave is linearly degenerated. The fields associated to $\pm c$ are genuinely non-linear under the condition that the following quantity does not vanish

[TABLE]

Appendix B Non-conservative energy

scheme

To obtain the non-conservative scheme, we do not need anymore the relaxation on the gravitational potential. This scheme is then obtained through the following splitting, for the acoustic subsystem

[TABLE]

followed by the transport subsystem

[TABLE]

then we use the same techniques for the acoustic system as in 2.1.2, in other words the use of the mass variable and the Lagrangian variables. The acoustic system in these variables writes

[TABLE]

[TABLE]

Using a pressure relaxation, an approximate Riemann solver with source term, see Gallice (2002), and the same upwind scheme for the transport system as in 2.1.3 we obtain the following non-conservative counterpart scheme

[TABLE]

where

[TABLE]

This scheme is not conservative for the whole energy but is closer to the scheme proposed for the shallow water equations in Chalons et al. (2017), for which the authors have obtained a discrete entropy inequality. It seems therefore possible to obtain a similar inequality for this non-conservative scheme, but this demonstration is beyond the scope of this paper.