APSARA: A multi-dimensional unsplit fourth-order explicit Eulerian hydrodynamics code for arbitrary curvilinear grids

A. Wongwathanarat; H. Grimm-Strele; E. Müller

doi:10.1051/0004-6361/201628205

Home

All issues

Volume 595 (November 2016)

A&A, 595 (2016) A41

Full HTML

Free Access

Issue		A&A Volume 595, November 2016


Article Number		A41
Number of page(s)		16
Section		Numerical methods and codes
DOI		https://doi.org/10.1051/0004-6361/201628205
Published online		26 October 2016

A&A 595, A41 (2016)

APSARA: A multi-dimensional unsplit fourth-order explicit Eulerian hydrodynamics code for arbitrary curvilinear grids

A. Wongwathanarat¹^,2, H. Grimm-Strele³^,2^,⋆ and E. Müller²

¹ RIKEN, Astrophysical Big Bang Laboratory, 2-1 Hirosawa, Wako, 351-0198 Saitama, Japan
e-mail: annop.wongwathanarat@riken.jp
² Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85748 Garching, Germany
³ Institute of Mathematics, University of Vienna, Oskar-Morgenstern-Platz 1, 1090 Vienna, Austria

Received: 28 January 2016
Accepted: 11 July 2016

Abstract

We present a new fourth-order, finite-volume hydrodynamics code named Apsara. The code employs a high-order, finite-volume method for mapped coordinates with extensions for nonlinear hyperbolic conservation laws. Apsara can handle arbitrary structured curvilinear meshes in three spatial dimensions. The code has successfully passed several hydrodynamic test problems, including the advection of a Gaussian density profile and a nonlinear vortex and the propagation of linear acoustic waves. For these test problems, Apsara produces fourth-order accurate results in case of smooth grid mappings. The order of accuracy is reduced to first-order when using the nonsmooth circular grid mapping. When applying the high-order method to simulations of low-Mach number flows, for example, the Gresho vortex and the Taylor-Green vortex, we discover that Apsara delivers superior results to codes based on the dimensionally split, piecewise parabolic method (PPM) widely used in astrophysics. Hence, Apsara is a suitable tool for simulating highly subsonic flows in astrophysics. In the first astrophysical application, we perform implicit large eddy simulations (ILES) of anisotropic turbulence in the context of core collapse supernova (CCSN) and obtain results similar to those previously reported.

Key words: methods: numerical / hydrodynamics / turbulence

^⋆

Present address: NUMECA International, Chaussée de la Hulpe, 189, Terhulpsesteenweg, 1170 Brussels, Belgium.

© ESO, 2016

1. Introduction

In many astrophysical simulations, such as simulations of convection inside a star and stellar explosions, spherical coordinates are often a preferable choice for integrating the respective partial differential equations. However, the spherical polar coordinates possess coordinate singularities at the coordinate origin and along the north and south poles, which prevent an easy and straightforward implementation of numerical methods. These coordinate singularities result in smaller time steps and, hence, special numerical treatments and/or boundary conditions have to be applied. This leads to a degradation of the efficiency of the employed numerical schemes and may introduce numerical artifacts.

For example, in recent state-of-the-art, three-dimensional (3D) simulations of core collapse supernovae (CCSN) performed with time-explicit, finite-volume hydrodynamic codes (e.g., Melson et al. 2015a; Lentz et al. 2015), the flow is modeled in spherical symmetry inside a specified sphere representing the inner core of the proto-neutron star, where the symmetry assumption is justified, to alleviate the restriction of the time step due to the Courant-Friedrichs-Lewy (CFL) condition. In addition, a reflecting boundary condition is usually applied at the coordinate origin. Müller (2015) applied the same spherical symmetry assumption, but he also included the effect of proto-neutron star convection by means of a mixing-length theory (e.g., Wilson & Mayle 1988). Lentz et al. (2015) avoided time steps that were too small by using nonuniform angular zones in the polar direction and an azimuthal averaging procedure for grid zones near the two poles. Similarly, Müller (2015) circumvented this problem using a mesh coarsening scheme in 30° cones around the poles whereby the short wavelength noise in the azimuthal direction is filtered out.

The problem of the severe time step restriction near the singular points at the poles of a sphere, also known as the pole-problem, has received much attention in various fields of research over the past few decades. In the particular case of finite-volume methods on a structured mesh, two popular solutions to the pole-problem are the cubed sphere grid (Ronchi et al. 1996) and the Yin-Yang grid (Kageyama & Sato 2004). The cubed sphere grid is based on a projection of the six sides of a cube onto the surface of a sphere and thus consists of six equidistant (in polar and azimuthal direction) identical grid patches. Except for the one point at the middle of each grid patch, the cubed sphere mesh is nonorthogonal. On the other hand, the Yin-Yang grid is formed by combining two identical low-latitude parts of a spherical polar grid. Hence, it is orthogonal everywhere, allowing for an easy extension of existing numerical codes that are based on an orthogonal mesh. Nevertheless, overlapping grids like the Yin-Yang grid have a drawback. Even if the numerical scheme employed on each grid patch is conservative, this does not ensure global conservation in the computational volume unless one applies a flux correction algorithm at the boundaries between the two grid patches (see, e.g., Peng et al. 2006; Wongwathanarat et al. 2010). In contrast, the boundaries of each patch of the cubed sphere mesh coincide perfectly with its neigboring patches, i.e., there is no overlap between grid patches. Thus, boundary flux corrections can be implemented in a more straightforward manner than in the case of the Yin-Yang grid.

Both the cubed sphere grid and the Yin-Yang grid have been used in various astrophysical applications. For instance, the cubed sphere grid was applied in simulations of accretion flows onto magnetized stars (e.g., Koldoba et al. 2002; Romanova et al. 2012) and accretion disks around rotating black holes (Fragile et al. 2009). The Yin-Yang grid was used, for example, in CCSN simulations (e.g., Wongwathanarat et al. 2015; Melson et al. 2015b), simulations of type Ia supernova remnants (Warren & Blondin 2013), and calculations of the solar global convection (Hotta et al. 2014), solar corona (e.g., Jiang et al. 2012; Feng et al. 2012), and coronal mass ejections (Shiota et al. 2010).

While the cubed sphere grid and the Yin-Yang grid circumvent the pole-problem, the singularity at the coordinate origin remains in both cases. One possible solution is to supplement the cubed sphere mesh or the Yin-Yang mesh with a Cartesian mesh at small radii to cover the central part of the sphere. An example of such a grid arrangement, called the Yin-Yang-Zhong, was recently developed by Hayashi & Kageyama (2016). Using this approach, the numerical code must be able to deal with different coordinate systems on different grid patches. In addition, global conservation is hampered by the overlap of the Cartesian grid patch with the spherical grid patches. On the other hand, the GenASiS code developed by Cardall et al. (2014) handles both the pole-problem and the singularity at the coordinate origin by employing the adaptive mesh refinement (AMR) approach. Instead of a block-structured AMR framework, GenASiS uses a more flexible cell-by-cell refinement in the Cartesian coordinate system. The cell-by-cell refinement can generate a centrally refined mesh, achieving higher and higher resolutions toward the origin of the sphere, which is a desirable grid property for many astrophysical applications. Nevertheless, a possible disadvantage of the cell-by-cell refinement approach is the cost of a complicated data communication on large-scale machines.

Calhoun et al. (2008) proposed grid mappings for circular or spherical domains that work on a logically rectangular mesh. These grid mappings can be used together with the mapped-grid technique in which one formulates the governing equations in an abstract (singularity-free) computational space instead of in a physical space and applies a coordinate transformation between them. A great advantage of the mapped-grid approach is that the computational domain can always be discretized by an equidistant Cartesian mesh. Therefore, any numerical scheme formulated for a Cartesian equidistant grid can be applied in a straightforward manner. The mapped-grid technique is widely used in engineering applications, where complicated geometries are common (see, e.g., LeVeque 2002). In the field of astrophysics, the mapped-grid approach is not very well known and has only been used up to now in a few numerical codes (Kifonidis & Müller 2012; Grimm-Strele et al. 2014; Miczek et al. 2015).

Colella et al. (2011) introduced a new class of high-order (better than second-order), finite-volume methods for mapped coordinates. The methods are based on computing face-averaged fluxes at each face of a control volume using high-order quadrature rules rather than the mid-point rule as in standard second-order accurate schemes. They demonstrated the capability of their approach to solve an elliptic equation and a scalar, linear hyperbolic equation up to fourth-order accuracy. An extension of the method to solve nonlinear hyperbolic equations was presented for Cartesian coordinates by McCorquodale & Colella (2011) and for mapped coordinates by Guzik et al. (2012). This extension is nontrivial because it is necessary to perform nonlinear transformations between point and (zone and face) averaged values of the conserved variables and fluxes ensuring fourth-order accuracy.

Motivated by the works of Calhoun et al. (2008) and Colella et al. (2011), we have developed a new numerical code for astrophysical applications. This code, which we named Apsara, is based on the fourth-order implementation of the high-order, finite-volume method in mapped coordinates by Colella et al. (2011). The code is extended to solve the Euler equations of gas dynamics using a high-resolution, shock-capturing (HRSC) method on arbitrary, structured curvilinear grids in three spatial dimensions as described in Guzik et al. (2012). This code robustly captures shock waves and discontinuities, while obtaining high-order accurate results in smooth flow regions. On the other hand, simulations of flows on complex grid geometry are possible thanks to the mapped-grid technique. These combined features make the code suitable for simulating a wide range of astrophysical problems.

The paper is organized as follows. In Sect. 2, we summarize the methods by Colella et al. (2011), McCorquodale & Colella (2011), and Guzik et al. (2012) for simulating 3D compressible flows in mapped coordinates. In Sect. 3, examples of mapping functions are introduced, which are employed in our test calculations. In Sect. 4, using a set of hydrodynamic tests, we demonstrate the capability of Apsara to compute fourth-order accurate solutions for smooth flows without discontinuities. In addition, we show the ability of Apsara to accurately calculate low-Mach number flows. As an astrophysical application, we also performed turbulence simulations in the context of CCSN with Apsara. We discuss the results of these simulations in Sect. 5. Finally, we summarize the code properties in Sect. 6, and mention some plans for future extensions of Apsara.

2. Numerical methods

2.1. Governing equations

The Euler equations written in Cartesian coordinates (x,y,z) in three spatial dimensions read $\frac{\partial U}{∂t} + \frac{\partial F (U)}{∂x} + \frac{\partial G (U)}{∂y} + \frac{\partial H (U)}{∂z} = 0,$ $\begin{equation} \frac{\partial\mbf{U}}{\partial t} +\frac{\partial\mbf{F(U)}}{\partial x} +\frac{\partial\mbf{G(U)}}{\partial y} +\frac{\partial\mbf{H(U)}}{\partial z} = 0, \label{eq:eulereqs} \end{equation}$ (1)where the vector of conserved variables U and the Cartesian flux functions F, G, and H are given by $\begin{matrix} U & = & (\begin{matrix} \end{matrix}), F (U) = (\begin{matrix} \end{matrix}), \\ G (U) & = \end{matrix}$ $\begin{eqnarray} \mbf{U}&=&\left(\begin{array}{c} \rho \\ \rho u \\ \rho v \\ \rho w \\ \rho E \end{array}\right),~~~~~~~~~~~ \mbf{F(U)}=\left(\begin{array}{c} \rho u \\ \rho u^2+p \\ \rho uv \\ \rho uw \\ \rho uH \end{array}\right), \nonumber\\ \mbf{G(U)}&=&\left(\begin{array}{c} \rho v \\ \rho uv \\ \rho v^2+p \\ \rho vw \\ \rho vH \end{array}\right),~~~~~~~ \mbf{H(U)}=\left(\begin{array}{c} \rho w \\ \rho uw \\ \rho vw \\ \rho w^2+p \\ \rho wH \end{array}\right). \label{eq:def-UFGH} \end{eqnarray}$ (2)The hydrodynamic state variables in Eq. (2) are the density ρ, the Cartesian components of the velocity field v = (u,v,w)^T, the pressure p, the specific total enthalpy $H = e + \frac{1}{2} (u^{2} + v^{2} + w^{2}) + \frac{p}{ρ},$ $\begin{equation} H = e+\frac{1}{2}\left(u^2+v^2+w^2\right)+\frac{p}{\rho}, \label{eq:enthalpy} \end{equation}$ (3)and the specific internal energy e. The pressure p is related to ρ and e via the equation of state (EOS) of a perfect gas $p = ρe (γ - 1),$ $\begin{equation} p = \rho e(\gamma-1), \label{eq:eos} \end{equation}$ (4)where γ is the adiabatic index.

We solve Eq. (1) using a finite-volume approach. The domain of interest is discretized into small subvolumes V_i, where the subscript i = (i,j,k) are indices referring to the ith, jth, and kth cell in each coordinate direction. We integrate Eq. (1) for each of the control volumes V_i to obtain the Euler equations in their integral form. For each component of the vector of conserved variables U we have $\frac{\partial}{∂t} \int^{V_{i}} U^{s} d V + \int^{V_{i}} \nabla \cdot F^{s} d V = 0,$ $\begin{equation} \frac{\partial}{\partial t}\int\limits_{V_{\bs{i}}}\, {U}^s\,\mrm{d}V + \int\limits_{V_{\bs{i}}}\, \bs{\nabla}\cdot\bs{\mfk{F}}^s\,\mrm{d}V = 0, \label{eq:inteuler} \end{equation}$ (5)where U^s denotes the sth component of U, and F^s = (F^s,G^s,H^s)^T is the corresponding flux vector. We then define a mapping function M, which maps coordinates ξ = (ξ,η,ζ)^T one-to-one in an abstract computational space into coordinates x = (x,y,z)^T in physical space, i.e., $M (ξ) = x, M : R^{3} \to R^{3} .$ $\begin{equation} \mbf{M}(\bs{\xi}) = \mbf{x},\; \mbf{M}\!\!:\mathbb{R}^3\rightarrow\mathbb{R}^3. \end{equation}$ (6)For simplicity we discretize the computational space with a Cartesian mesh consisting of computational cells Ω_i, which are unit cubes. The mapping function M is defined such that the image M(Ω_i) of each computational cell Ω_i corresponds to the control volume V_i in physical space. By applying the change of variables theorem (see, e.g., Appel 2007) to Eq. (5), the integrals over a control volume V_i in physical space are transformed into integrals over a computational cell Ω_i according to $\frac{\partial}{∂t} \int^{Ω_{i}} J U^{s} dΩ + \int^{Ω_{i}} J (\nabla \cdot F^{s}) dΩ = 0,$ $\begin{equation} \frac{\partial}{\partial t} \int\limits_{\Omega_{\bs{i}}}\, J{U}^s\,\mrm{d}\Omega + \int\limits_{\Omega_{\bs{i}}}\, J(\bs{\nabla}\cdot\bs{\mfk{F}}^s)\,\mrm{d}\Omega = 0, \label{eq:inteuler2} \end{equation}$ (7)where $J = | \frac{\partial (x,y,z)}{\partial (ξ,η,ζ)} |$ $\begin{equation} J=\left|\frac{\partial(x,y,z)}{\partial(\xi,\eta,\zeta)}\right| \label{eq:detj} \end{equation}$ (8)is the Jacobian determinant of the inverse mapping M^-1. The divergence of F^s can be expressed in terms of derivatives in the computational space (Colella et al. 2011), $\nabla \cdot F^{s} = \frac{1}{J} \nabla_{ξ} \cdot (N F^{s}),$ $\begin{equation} \bs{\nabla}\cdot\bs{\mfk{F}}^s=\frac{1}{J}\bs{\nabla_\xi} \cdot({\bf N}\bs{\mfk{F}}^s), \label{eq:divergence} \end{equation}$ (9)the matrix N being defined as $N = (\begin{matrix} | \frac{\partial (y,z)}{\partial (η,ζ)} | & | \frac{\partial (z,x)}{\partial (η,ζ)} | & | \frac{\partial (x,y)}{\partial (η,ζ)} | \\ | \frac{\partial (y,z)}{\partial (ζ,ξ)} | & | \frac{\partial (z,x)}{\partial (ζ,ξ)} | & | \frac{\partial (x,y)}{\partial (ζ,ξ)} | \\ | \frac{\partial (y,z)}{\partial (ξ,η)} | & | \frac{\partial (z,x)}{\partial (ξ,η)} | & | \frac{\partial (x,y)}{\partial (ξ,η)} | \end{matrix}) .$ $\begin{equation} {\bf N}=\left( {\def\arraystretch{2.5} \begin{array}{ccc} \left|\dfrac{\partial(y,z)}{\partial(\eta, \zeta)}\right| & \left|\dfrac{\partial(z,x)}{\partial(\eta, \zeta)}\right| & \left|\dfrac{\partial(x,y)}{\partial(\eta, \zeta)}\right| \\ \left|\dfrac{\partial(y,z)}{\partial(\zeta,\xi) }\right| & \left|\dfrac{\partial(z,x)}{\partial(\zeta,\xi) }\right| & \left|\dfrac{\partial(x,y)}{\partial(\zeta,\xi) }\right| \\ \left|\dfrac{\partial(y,z)}{\partial(\xi, \eta) }\right| & \left|\dfrac{\partial(z,x)}{\partial(\xi, \eta) }\right| & \left|\dfrac{\partial(x,y)}{\partial(\xi, \eta) }\right| \\ \end{array} }\right). \label{eq:n-matrix} \end{equation}$ (10)Using Eq. (9) we obtain from Eq. (7) $\frac{\partial}{∂t} \int^{Ω_{i}} J U^{s} dΩ + \int^{Ω_{i}} \nabla_{ξ} \cdot (N F^{s}) dΩ = 0,$ $\begin{equation} \frac{\partial}{\partial t} \int\limits_{\Omega_{\bs{i}}}\, J{U}^s\,\mrm{d}\Omega + \int\limits_{\Omega_{\bs{i}}}\, \bs{\nabla_\xi}\cdot({\bf N}\bs{\mfk{F}}^s)\,\mrm{d}\Omega = 0, \label{eq:transeuler} \end{equation}$ (11)which resembles Eq. (5). Applying the divergence theorem Eq. (11) becomes $\frac{\partial}{∂t} \int^{Ω_{i}} J U^{s} dΩ + \int^{\partial Ω_{i}} N F^{s} \cdot d A_{ξ} = 0.$ $\begin{equation} \frac{\partial}{\partial t} \int\limits_{\Omega_{\bs{i}}}\,J{U}^s\,\mrm{d}\Omega + \int\limits_{\partial\Omega _{\bs{i}}}\,{\bf N}\bs{\mfk{F}}^s\cdot\,\mrm{d}\mbf{A}_{\bs{\xi}}=0. \label{eq:transeuler2} \end{equation}$ (12)The second integral in Eq. (12) is performed over all faces ∂Ω_i of a cell Ω_i, where A_ξ is the surface vector of each cell face pointing in the outward normal direction. The integrals in Eq. (12) are approximated to fourth-order accuracy as $\frac{\partial}{∂t} ⟨ J U^{s} ⟩_{i} + \sum_{d = 1}^{3} (⟨ N_{d} F^{s} ⟩_{i + \frac{1}{2} e^{d}} - ⟨ N_{d} F^{s} ⟩_{i - \frac{1}{2} e^{d}}) = 0,$ $\begin{equation} \frac{\partial}{\partial t}\avg{J{U}^s}_{\bs{i}} + \sum_{d=1}^3\left( \avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}+\frac{1}{2}\bs{e}^d} - \avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}-\frac{1}{2}\bs{e}^d} \right)=0, \label{eq:avgeuler} \end{equation}$ (13)where e^d is the unit vector in the dth direction and N_d is the dth row of the matrix N. The operators ⟨ ·⟩ _i and $⟨ \cdot ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{\cdot}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ denote the fourth-order accurate approximations of volume averages and face averages, respectively. We omit factors of h³ and h² multiplying the first and second term, respectively, on the LHS of Eq. (13) owing to the cell spacing h = 1.

2.2. Temporal discretization

In accordance with our fourth-order accurate spatial scheme, we use a fourth-order accurate Runge-Kutta (RK4) scheme for the time integration. Discretizing the time derivative in Eq. (13) one obtains the following update formula for the cell averaged conserved quantities ⟨ JU^s ⟩ _i, $⟨ J U^{s} ⟩_{i}^{n + 1} = ⟨ J U^{s} ⟩_{i}^{n} + Δ t^{n} \sum_{d = 1}^{3} (⟨ N_{d} F^{s} ⟩_{i - \frac{1}{2} e^{d}}^{tot} - ⟨ N_{d} F^{s} ⟩_{i + \frac{1}{2} e^{d}}^{tot}),$ $\begin{equation} \avg{J{U}^s}^{n+1}_{\bs{i}} = \avg{J{U}^s}^n_{\bs{i}} + \Delta t^n\sum_{d=1}^3 ( \avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}-\frac{1}{2}\bs{e}^d}^\mrm{tot} -\avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}+\frac{1}{2}\bs{e}^d}^\mrm{tot} ), \label{eq:update} \end{equation}$ (14)where $⟨ J U^{s} ⟩_{i}^{n}$ $\hbox{$\avg{J{U}^s}^n_{\bs{i}}$}$ and $⟨ J U^{s} ⟩_{i}^{n + 1}$ $\hbox{$\avg{J{U}^s}^{n+1}_{\bs{i}}$}$ represent the cell averaged conserved quantities at times tⁿ and t^{n + 1} = tⁿ + Δtⁿ, respectively. The total RK4 face averaged fluxes $⟨ N_{d} F^{s} ⟩_{i \pm \frac{1}{2} e^{d}}^{tot}$ $\hbox{$\avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}^\mrm{tot}$}$ are given by the linear combination $\begin{matrix} ⟨ N_{d} F^{s} ⟩_{i \pm \frac{1}{2} e^{d}}^{tot} = \\ \frac{1}{6} (⟨ N_{d} F^{s} ⟩_{i \pm \frac{1}{2} e^{d}}^{(0)} + 2 ⟨ N_{d} F^{s} ⟩_{i \pm \frac{1}{2} e^{d}}^{(1)} + 2 ⟨ N_{d} F^{s} ⟩_{i \pm \frac{1}{2} e^{d}}^{(2)} + ⟨ N_{d} F^{s} ⟩_{i \pm \frac{1}{2} e^{d}}^{(3)}), \end{matrix}$ $\begin{eqnarray} &&\avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}^\mrm{tot} = \nonumber\\ &&\,\,\,\,\frac{1}{6}\left( \avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}^{(0)} +2\avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}^{(1)} +2\avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}^{(2)} + \avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}^{(3)} \right), \end{eqnarray}$ (15)where $\begin{matrix} ⟨ J U^{s} ⟩ \begin{matrix} (0) \\ i \end{matrix} = ⟨ J U^{s} ⟩_{i}^{n}, \\ ⟨ J U^{s} ⟩ \begin{matrix} (1) \\ i \end{matrix} = ⟨ J U^{s} ⟩ \begin{matrix} (0) \\ i \end{matrix} + \frac{Δ t^{n}}{2} \sum_{d = 1}^{3} (⟨ N_{d} F^{s} ⟩ \begin{matrix} (0) \\ i - \frac{1}{2} e^{d} \end{matrix} - ⟨ N_{d} F^{s} ⟩ \begin{matrix} (0) \\ i + \frac{1}{2} e^{d} \end{matrix}), \\ ⟨ J U^{s} ⟩ \begin{matrix} (2) \\ i \end{matrix} = ⟨ J U^{s} ⟩ \begin{matrix} (0) \\ i \end{matrix} + \frac{Δ t^{n}}{2} \sum_{d = 1}^{3} (⟨ N_{d} F^{s} ⟩ \begin{matrix} (1) \\ i - \frac{1}{2} e^{d} \end{matrix} - ⟨ N_{d} F^{s} ⟩ \begin{matrix} (1) \\ i + \frac{1}{2} e^{d} \end{matrix}), \\ ⟨ J U^{s} ⟩ \begin{matrix} (3) \\ i \end{matrix} = ⟨ J U^{s} ⟩ \begin{matrix} (0) \\ i \end{matrix} + Δ t^{n} \sum_{d = 1}^{3} (⟨ N_{d} F^{s} ⟩ \begin{matrix} (2) \\ i - \frac{1}{2} e^{d} \end{matrix} - ⟨ N_{d} F^{s} ⟩ \begin{matrix} (2) \\ i + \frac{1}{2} e^{d} \end{matrix}), \end{matrix}$ $\begin{eqnarray} &&\avg{J{U}^s}^{\left(0\right)}_{\bs{i}} = \avg{J{U}^s}^n_{\bs{i}},~~~~~~~ \\ &&\avg{J{U}^s}^{\left(1\right)}_{\bs{i}} = \avg{J{U}^s}^{\left(0\right)}_{\bs{i}} + \frac{\Delta t^n}{2}\sum_{d=1}^3 \left( \avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}-\frac{1}{2}\bs{e}^d}^{\left(0\right)} -\avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}+\frac{1}{2}\bs{e}^d}^{\left(0\right)} \right), \\ &&\avg{J{U}^s}^{\left(2\right)}_{\bs{i}} = \avg{J{U}^s}^{\left(0\right)}_{\bs{i}} + \frac{\Delta t^n}{2}\sum_{d=1}^3 \left( \avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}-\frac{1}{2}\bs{e}^d}^{\left(1\right)} -\avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}+\frac{1}{2}\bs{e}^d}^{\left(1\right)} \right), \\ &&\avg{J{U}^s}^{\left(3\right)}_{\bs{i}} = \avg{J{U}^s}^{\left(0\right)}_{\bs{i}} + \Delta t^n\sum_{d=1}^3 \left( \avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}-\frac{1}{2}\bs{e}^d}^{\left(2\right)} -\avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}+\frac{1}{2}\bs{e}^d}^{\left(2\right)} \right), \end{eqnarray}$ and $⟨ N_{d} F^{s} ⟩_{i \pm \frac{1}{2} e^{d}}^{(m)} = ⟨ N_{d} F^{s} (U^{(m)}) ⟩_{i \pm \frac{1}{2} e^{d}}, for m = 0, 1, 2, 3.$ $\begin{equation} \avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}^{(m)} = \avg{{\bf N}_d\bs{\mfk{F}}^s(\mbf{U}^{(m)})}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}, \quad \mrm{for} \quad m=0,1,2,3. \end{equation}$ (20)The size of the time step Δt is given by the condition (McCorquodale & Colella 2011) $\frac{Δ t}{h} \max_{i} (\sum_{d = 1}^{3} J^{-1} (| N v \cdot e^{d} | + c_{s} | {(N_{d})}^{T} |^{)}) ⪅ 1.3925,$ $\begin{equation} \frac{\Delta t}{h}\max\limits_{\bs{i}} \left(\sum_{d=1}^3 J^{-1}\left(\lvert{\bf N}\mbf{v}\cdot\bs{e}^d\rvert + c_{\rm s}\lvert\left({\bf N}_d\right)^T\rvert\right) \right) \lessapprox 1.3925, \label{eq:tcfl} \end{equation}$ (21)where c_s is the speed of sound.

2.3. Cell and face averages

A fourth-order approximation of the cell average of a function f on a mapped grid can be computed by performing a Taylor expansion of the function about the cell center in the computational space. The expansion yields $⟨ f ⟩_{i} = f_{i} + \frac{h^{2}}{24} \nabla_{ξ}^{2} f,$ $\begin{equation} \avg{f}_{\bs{i}}=f_{\bs{i}}+\frac{h^2}{24}\nabla_{\bs{\xi}}^2 f, \label{eq:cellavg} \end{equation}$ (22)where f_i denotes the pointwise value of f at the center of Ω_i. Because the Laplacian operator in Eq. (22) only needs to be evaluated to second-order accuracy, we use the second-order accurate central difference formula to calculate it. On the equidistant Cartesian grid in the computational space, we define the quantities $𝒞_{i}^{(2)} (a,q) = \sum_{d = 1}^{3} (a_{i + q e^{d}} - 2 a_{i} + a_{i - q e^{d}}),$ $\begin{equation} \mathcal{C}^{(2)}_{\bs{i}}(a,q) = \sum_{d=1}^3 (a_{\bs{i}+q\bs{e}^d}-2a_{\bs{i}}+a_{\bs{i}-q\bs{e}^d}), \label{eq:c2} \end{equation}$ (23)to calculate the cell averages of a function f as $⟨ f ⟩_{i} = f_{i} + \frac{1}{24} 𝒞_{i}^{(2)} (f, 1) .$ $\begin{equation} \avg{f}_{\bs{i}} = f_{\bs{i}} + \frac{1}{24}\mathcal{C}^{(2)}_{\bs{i}}(f,1). \label{eq:cellfto<f>} \end{equation}$ (24)Conversely, we can obtain a fourth-order approximation of the f_i from cell averages ⟨ f ⟩. In this case, the Laplacian is calculated using cell averages ⟨ f ⟩, which yields $f_{i} = ⟨ f ⟩_{i} - \frac{1}{24} 𝒞_{i}^{(2)} (⟨ f ⟩, 1) .$ $\begin{equation} f_{\bs{i}} = \avg{f}_{\bs{i}} - \frac{1}{24} \mathcal{C}^{(2)}_{\bs{i}}(\avg{f},1). \label{eq:cell<f>tof} \end{equation}$ (25)Similarly, expanding f about the center of a cell face in the dth direction, we obtain fourth-order accurate face values $⟨ f ⟩_{i \pm \frac{1}{2} e^{d}} = f_{i \pm \frac{1}{2} e^{d}} + \frac{h^{2}}{24} \nabla_{ξ, ⊥}^{2} f,$ $\begin{equation} \avg{f}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} = f_{\bs{i}\pm\frac{1}{2}\bs{e}^d} + \frac{h^2}{24}\nabla_{\bs{\xi},\perp}^2 f, \end{equation}$ (26)where $\nabla_{ξ, ⊥}^{2}$ $\hbox{$\nabla_{\bs{\xi},\perp}^2$}$ denotes the transverse Laplacian operator. As before, we replace the transverse Laplacian by $𝒞_{i \pm \frac{1}{2} e^{d}, ⊥}^{(2)} (a,q) = \begin{matrix} 3 \\ \sum \\ d^{'} = 1 d^{'} \neq d \end{matrix} (a_{i \pm \frac{1}{2} e^{d} + q e^{d^{'}}} - 2 a_{i \pm \frac{1}{2} e^{d}} + a_{i \pm \frac{1}{2} e^{d} - q e^{d^{'}}})$ $\begin{equation} \mathcal{C}^{(2)}_{\bs{i}\pm\frac{1}{2}\bs{e}^d,\perp}(a,q)= \sum_{\substack{d^\prime=1\\d^\prime\ne d}}^3 (a_{\bs{i}\pm\frac{1}{2}\bs{e}^d+q\bs{e}^{d^\prime}} -2a_{\bs{i}\pm\frac{1}{2}\bs{e}^d} +a_{\bs{i}\pm\frac{1}{2}\bs{e}^d-q\bs{e}^{d^\prime}}) \label{eq:c2p} \end{equation}$ (27)to obtain $⟨ f ⟩_{i \pm \frac{1}{2} e^{d}} \approx f_{i \pm \frac{1}{2} e^{d}} + \frac{1}{24} 𝒞_{i \pm \frac{1}{2} e^{d}, ⊥}^{(2)} (f, 1) .$ $\begin{equation} \avg{f}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} \approx f_{\bs{i}\pm\frac{1}{2}\bs{e}^d} + \frac{1}{24}\mathcal{C}^{(2)}_{\bs{i}\pm\frac{1}{2}\bs{e}^d,\perp}(f,1). \label{eq:facefto<f>} \end{equation}$ (28)To calculate pointwise values of f at face centers from face averages we use $f_{i \pm \frac{1}{2} e^{d}} = ⟨ f ⟩_{i \pm \frac{1}{2} e^{d}} - \frac{1}{24} 𝒞_{i \pm \frac{1}{2} e^{d}, ⊥}^{(2)} (⟨ f ⟩, 1) .$ $\begin{equation} f_{\bs{i}\pm\frac{1}{2}\bs{e}^d} = \avg{f}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} - \frac{1}{24} \mathcal{C}^{(2)}_{\bs{i}\pm\frac{1}{2}\bs{e}^d,\perp}(\avg{f},1). \label{eq:face<f>tof} \end{equation}$ (29)In addition, we also need expressions for cell averages and face averages of a product of two functions, for example, for the products JU^s and N_dF^s that appear in the update formula (Eq. (13). The Taylor expansion of the product of two functions f and g about the cell center gives $⟨ fg ⟩_{i} = ⟨ f ⟩_{i} ⟨ g ⟩_{i} + \frac{h^{2}}{12} \sum_{d = 1}^{3} \frac{∂f}{\partial ξ_{d}} \frac{∂g}{\partial ξ_{d}},$ $\begin{equation} \avg{fg}_{\bs{i}} = \avg{f}_{\bs{i}}\avg{g}_{\bs{i}} + \frac{h^2}{12} \sum_{d=1}^3 \frac{\partial f}{\partial\xi_d} \frac{\partial g}{\partial\xi_d}, \label{eq:prodavg} \end{equation}$ (30)where the first-order partial derivatives can be replaced again by second-order accurate central differences. We define $𝒞_{i}^{(1)} (a,b,q) = \sum_{d = 1}^{3} (a_{i + q e^{d}} - a_{i - q e^{d}}) (b_{i + q e^{d}} - b_{i - q e^{d}}),$ $\begin{equation} \mathcal{C}^{(1)}_{\bs{i}}(a,b,q) = \sum_{d=1}^3 (a_{\bs{i}+q\bs{e}^d} - a_{\bs{i}-q\bs{e}^d}) (b_{\bs{i}+q\bs{e}^d} - b_{\bs{i}-q\bs{e}^d}), \label{eq:c1} \end{equation}$ (31)and compute ⟨ fg ⟩ _i as $⟨ fg ⟩_{i} = ⟨ f ⟩_{i} ⟨ g ⟩_{i} + \frac{1}{48} 𝒞_{i}^{(1)} (⟨ f ⟩, ⟨ g ⟩, 1)$ $\begin{equation} \avg{fg}_{\bs{i}} = \avg{f}_{\bs{i}}\avg{g}_{\bs{i}} + \frac{1}{48} \mathcal{C}^{(1)}_{\bs{i}}(\avg{f},\avg{g},1) \label{eq:cell<fg>} \end{equation}$ (32)Similarly, the face averages are given by $⟨ fg ⟩_{i \pm \frac{1}{2} e^{d}} = ⟨ f ⟩_{i \pm \frac{1}{2} e^{d}} ⟨ g ⟩_{i \pm \frac{1}{2} e^{d}} + \frac{1}{48} 𝒞_{i \pm \frac{1}{2} e^{d}, ⊥}^{(1)} (⟨ f ⟩, ⟨ g ⟩, 1),$ $\begin{equation} \avg{fg}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} = \avg{f}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} \, \avg{g}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} + \frac{1}{48} \mathcal{C}^{(1)}_{\bs{i}\pm\frac{1}{2}\bs{e}^d,\perp}(\avg{f},\avg{g},1), \label{eq:face<fg>} \end{equation}$ (33)where $\begin{matrix} 𝒞_{i \pm \frac{1}{2} e^{d}, ⊥}^{(1)} (a,b,q) = \\ \begin{matrix} 3 \\ \sum \\ d^{'} = 1 d^{'} \neq d \end{matrix} (a_{i \pm \frac{1}{2} e^{d} + q e^{d^{'}}} - a_{i \pm \frac{1}{2} e^{d} - q e^{d^{'}}}) (b_{i \pm \frac{1}{2} e^{d} + q e^{d^{'}}} - b_{i \pm \frac{1}{2} e^{d} - q e^{d^{'}}}) . \end{matrix}$ $\begin{eqnarray} &&\mathcal{C}^{(1)}_{\bs{i}\pm\frac{1}{2}\bs{e}^d,\perp}(a,b,q) = \nonumber\\ &&\qquad \sum_{\substack{d^\prime=1\\d^\prime\ne d}}^3 (a_{\bs{i}\pm\frac{1}{2}\bs{e}^d+q\bs{e}^{d^\prime}} - a_{\bs{i}\pm\frac{1}{2}\bs{e}^d-q\bs{e}^{d^\prime}}) (b_{\bs{i}\pm\frac{1}{2}\bs{e}^d+q\bs{e}^{d^\prime}} - b_{\bs{i}\pm\frac{1}{2}\bs{e}^d-q\bs{e}^{d^\prime}}). \end{eqnarray}$ (34)Using Eq. (33) we can write the face average fluxes as $\begin{matrix} ⟨ N_{d} F^{s} ⟩_{i \pm \frac{1}{2} e^{d}} & = & ⟨ N_{d}^{1} ⟩_{i \pm \frac{1}{2} e^{d}} ⟨ F^{s} ⟩_{i \pm \frac{1}{2} e^{d}} + \frac{1}{48} 𝒞_{i \pm \frac{1}{2} e^{d}, ⊥}^{(1)} (⟨ N_{d}^{1} ⟩, ⟨ F^{s} ⟩, 1) \\ + ⟨ N_{d}^{2} ⟩_{i \pm \frac{1}{2} e^{d}} ⟨ G^{s} ⟩_{i \pm \frac{1}{2} e^{d}} + \frac{1}{48} 𝒞_{i \pm \frac{1}{2} e^{d}, ⊥}^{(1)} (⟨ N_{d}^{2} ⟩, ⟨ G^{s} ⟩, 1) \end{matrix}$ $\begin{eqnarray} \avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} & =\ & \avg{N_d^1}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} \avg{F^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} + \frac{1}{48} \mathcal{C}^{(1)}_{\bs{i}\pm\frac{1}{2}\bs{e}^d,\perp} (\avg{N_d^1},\avg{F^s},1)\nonumber\\ &&+ \avg{N_d^2}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} \avg{G^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} + \frac{1}{48}\mathcal{C}^{(1)}_{\bs{i}\pm\frac{1}{2}\bs{e}^d,\perp} (\avg{N_d^2},\avg{G^s},1)\nonumber\\ && + \avg{N_d^3}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} \avg{H^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} + \frac{1}{48}\mathcal{C}^{(1)}_{\bs{i}\pm\frac{1}{2}\bs{e}^d,\perp} (\avg{N_d^3},\avg{H^s},1), \label{eq:nf} \end{eqnarray}$ (35)where $N_{d}^{c}$ $\hbox{$N_d^c$}$ , with c = 1,2,3, is the cth component of the row matrix N_d.

The way one calculates $⟨ N_{d}^{c} ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{N_d^c}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ is an important step for a hydrodynamics code using a curvilinear grid, since it determines the freestream preservation property of the code. Consider a uniform initial state with ρ = p = 1 and u = v = w = 0 at every grid point. All conserved quantities do not change with time for such a state. Therefore, it follows from Eqs. (14) and (35) that $\sum_{\pm = +, -} \sum_{d = 1}^{3} \pm ⟨ N_{d}^{c} ⟩_{i \pm \frac{1}{2} e^{d}} = 0.$ $\begin{equation} \sum_{\pm=+,-}\sum_{d=1}^3\pm\avg{N_d^c}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}=0. \label{eq:freestream} \end{equation}$ (36)This condition is equivalent to the geometric identity that the sum of the surface normal vectors over all faces of a closed control volume must be zero. Violation of the freestream condition leads to errors that can destroy the numerical solution completely (see, e.g., Grimm-Strele et al. 2014).

If we define a matrix $\hbox{$\bs{\mathcal{N}}$}$ by (Colella et al. 2011) $𝒩 = \frac{1}{2} (\begin{matrix} y \frac{∂z}{∂ξ} - z \frac{∂y}{∂ξ} & z \frac{∂x}{∂ξ} - x \frac{∂z}{∂ξ} & x \frac{∂y}{∂ξ} - y \frac{∂x}{∂ξ} \\ y \frac{∂z}{∂η} - z \frac{∂y}{∂η} & z \frac{∂x}{∂η} - x \frac{∂z}{∂η} & x \frac{∂y}{∂η} - y \frac{∂x}{∂η} \\ y \frac{∂z}{∂ζ} - z \frac{∂y}{∂ζ} & z \frac{∂x}{∂ζ} - x \frac{∂z}{∂ζ} & x \frac{∂y}{∂ζ} - y \frac{∂x}{∂ζ} \end{matrix}),$ $\begin{equation} \bs{\mathcal{N}}=\frac{1}{2}\left( {\def\arraystretch{2.5} \begin{array}{ccc} y\dfrac{\partial z}{\partial\xi} -z\dfrac{\partial y}{\partial\xi} & z\dfrac{\partial x}{\partial\xi} -x\dfrac{\partial z}{\partial\xi} & x\dfrac{\partial y}{\partial\xi} -y\dfrac{\partial x}{\partial\xi} \\ y\dfrac{\partial z}{\partial\eta} -z\dfrac{\partial y}{\partial\eta} & z\dfrac{\partial x}{\partial\eta} -x\dfrac{\partial z}{\partial\eta} & x\dfrac{\partial y}{\partial\eta} -y\dfrac{\partial x}{\partial\eta} \\ y\dfrac{\partial z}{\partial\zeta}-z\dfrac{\partial y}{\partial\zeta} & z\dfrac{\partial x}{\partial\zeta}-x\dfrac{\partial z}{\partial\zeta} & x\dfrac{\partial y}{\partial\zeta}-y\dfrac{\partial x}{\partial\zeta} \\ \end{array} } \right), \label{eq:mathcaln-matrix} \end{equation}$ (37)it follows that $\nabla_{ξ} \times 𝒩^{c} = N^{c},$ $\begin{equation} \bs{\nabla_\xi} \times \bs{\mathcal{N}}^c = {\bf N}^c, \label{eq:crossmathcaln} \end{equation}$ (38)where $\hbox{$\bs{\mathcal{N}}^c$}$ and N^c are the cth column of the matrices $\hbox{$\bs{\mathcal{N}}$}$ and N, respectively. Using Eq. (38) and Stoke’s theorem, the face average $⟨ N_{d}^{c} ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{N_d^c}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ becomes $⟨ N_{d}^{c} ⟩_{i \pm \frac{1}{2} e^{d}} = \begin{matrix} 3 \\ \sum \\ d^{'} = 1 d^{'} \neq d \end{matrix} 􏽉 𝒩^{c} \cdot d E_{ξ_{d}},$ $\begin{equation} \avg{N_d^c}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} = \sum_{\substack{d^\prime=1\\d^\prime\ne d}}^3 \oint\bs{\mathcal{N}}^c \cdot\, \mrm{d}\mbf{E}_{\xi_d} \label{eq:intnde} , \end{equation}$ (39)where E_{ξ_d} is the right-handed tangent vector of the (hyper) edges of the cell face in the dth direction. For instance, $\begin{matrix} ⟨ N_{1}^{1} ⟩_{i \pm \frac{1}{2},j,k} = \\ \frac{1}{2} {\int_{η_{i \pm \frac{1}{2},j - \frac{1}{2},k - \frac{1}{2}}}^{η_{i \pm \frac{1}{2},j + \frac{1}{2},k - \frac{1}{2}}} (y \frac{∂z}{∂η} - z \frac{∂y}{∂η}) d η + \int_{ζ_{i \pm \frac{1}{2},j + \frac{1}{2},k - \frac{1}{2}}}^{ζ_{i \pm \frac{1}{2},j + \frac{1}{2},k + \frac{1}{2}}} (y \frac{∂z}{∂ζ} - z \frac{∂y}{∂ζ}) d ζ \end{matrix}$ $\begin{eqnarray} &&\avg{N_1^1}_{i\pm\frac{1}{2},j,k} = \nonumber\\ &&\quad \frac{1}{2}\left\{\;\int\limits_ {\eta_{i\pm\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2}}}^ {\eta_{i\pm\frac{1}{2},j+\frac{1}{2},k-\frac{1}{2}}}\, \left( y\frac{\partial z}{\partial\eta} - z\frac{\partial y}{\partial\eta} \right)\,\mrm{d}\eta + \int\limits_ {\zeta_{i\pm\frac{1}{2},j+\frac{1}{2},k-\frac{1}{2}}}^ {\zeta_{i\pm\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2}}}\, \left( y\frac{\partial z}{\partial\zeta} - z\frac{\partial y}{\partial\zeta} \right)\,\mrm{d}\zeta\right. \nonumber\\ &&\left. \quad - \int\limits_ {\eta_{i\pm\frac{1}{2},j-\frac{1}{2},k+\frac{1}{2}}}^ {\eta_{i\pm\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2}}}\, \left( y\frac{\partial z}{\partial\eta} - z\frac{\partial y}{\partial\eta} \right)\,\mrm{d}\eta - \int\limits_ {\zeta_{i\pm\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2}}}^ {\zeta_{i\pm\frac{1}{2},j-\frac{1}{2},k+\frac{1}{2}}}\, \left( y\frac{\partial z}{\partial\zeta} - z\frac{\partial y}{\partial\zeta} \right)\,\mrm{d}\zeta\,\right\}. \label{eq:<N11>} \end{eqnarray}$ (40)To calculate $⟨ N_{d}^{c} ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{N_d^c}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ , we apply Simpson’s rule to calculate the integrals in Eq. (39) and we approximate the first-order derivatives by second-order accurate finite differences. For example, the first integral on the right-hand side (RHS) of Eq. (40) is then given by $\begin{matrix} \int_{η_{i \pm \frac{1}{2},j - \frac{1}{2},k - \frac{1}{2}}}^{η_{i \pm \frac{1}{2},j + \frac{1}{2},k - \frac{1}{2}}} (y \frac{∂z}{∂η} - z \frac{∂y}{∂η}) d η = \\ \frac{1}{6} [y_{i \pm \frac{1}{2},j + \frac{1}{2},k - \frac{1}{2}} (3 z_{i \pm \frac{1}{2},j + \frac{1}{2},k - \frac{1}{2}} - 4 z_{i \pm \frac{1}{2},j,k - \frac{1}{2}} + z_{i \pm \frac{1}{2},j - \frac{1}{2},k - \frac{1}{2}}) \\ + 4 y_{i \pm \frac{1}{2},j,k - \frac{1}{2}} (z_{i \pm \frac{1}{2},j + \frac{1}{2},k - \frac{1}{2}} - z_{i \pm \frac{1}{2},j - \frac{1}{2},k - \frac{1}{2}}) \\ + y_{i \pm \frac{1}{2},j - \frac{1}{2},k - \frac{1}{2}} (- 3 z_{i \pm \frac{1}{2},j + \frac{1}{2},k - \frac{1}{2}} + 4 z_{i \pm \frac{1}{2},j,k - \frac{1}{2}} - z_{i \pm \frac{1}{2},j - \frac{1}{2},k - \frac{1}{2}}) \\ - z_{i \pm \frac{1}{2},j + \frac{1}{2},k - \frac{1}{2}} (3 y_{i \pm \frac{1}{2},j + \frac{1}{2},k - \frac{1}{2}} - 4 y_{i \pm \frac{1}{2},j,k - \frac{1}{2}} + y_{i \pm \frac{1}{2},j - \frac{1}{2},k - \frac{1}{2}}) \\ - 4 z_{i \pm \frac{1}{2},j,k - \frac{1}{2}} (y_{i \pm \frac{1}{2},j + \frac{1}{2},k - \frac{1}{2}} - y_{i \pm \frac{1}{2},j - \frac{1}{2},k - \frac{1}{2}}) \\ - z_{i \pm \frac{1}{2},j - \frac{1}{2},k - \frac{1}{2}} (- 3 y_{i \pm \frac{1}{2},j + \frac{1}{2},k - \frac{1}{2}} + 4 y_{i \pm \frac{1}{2},j,k - \frac{1}{2}} - y_{i \pm \frac{1}{2},j - \frac{1}{2},k - \frac{1}{2}})] . \end{matrix}$ $\begin{eqnarray} &&\int\limits_ {\eta_{i\pm\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2}}}^ {\eta_{i\pm\frac{1}{2},j+\frac{1}{2},k-\frac{1}{2}}}\, \left( y\frac{\partial z}{\partial\eta} - z\frac{\partial y}{\partial\eta} \right)\,\mrm{d}\eta = \nonumber\\ &&\quad\frac{1}{6}[y_{i\pm\frac{1}{2},j+\frac{1}{2},k-\frac{1}{2}} ( 3z_{i\pm\frac{1}{2},j+\frac{1}{2},k-\frac{1}{2}} -4z_{i\pm\frac{1}{2},j,k-\frac{1}{2}} + z_{i\pm\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2}}) \nonumber\\ && \quad +4y_{i\pm\frac{1}{2},j,k-\frac{1}{2}} ( z_{i\pm\frac{1}{2},j+\frac{1}{2},k-\frac{1}{2}} - z_{i\pm\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2}}) \nonumber\\ && \quad+ y_{i\pm\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2}} (-3z_{i\pm\frac{1}{2},j+\frac{1}{2},k-\frac{1}{2}} +4z_{i\pm\frac{1}{2},j,k-\frac{1}{2}} - z_{i\pm\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2}}) \nonumber\\ && \quad - z_{i\pm\frac{1}{2},j+\frac{1}{2},k-\frac{1}{2}} ( 3y_{i\pm\frac{1}{2},j+\frac{1}{2},k-\frac{1}{2}} -4y_{i\pm\frac{1}{2},j,k-\frac{1}{2}} + y_{i\pm\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2}}) \nonumber\\ && \quad -4z_{i\pm\frac{1}{2},j,k-\frac{1}{2}} ( y_{i\pm\frac{1}{2},j+\frac{1}{2},k-\frac{1}{2}} - y_{i\pm\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2}}) \nonumber\\ && \quad - z_{i\pm\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2}} (-3y_{i\pm\frac{1}{2},j+\frac{1}{2},k-\frac{1}{2}} +4y_{i\pm\frac{1}{2},j,k-\frac{1}{2}} - y_{i\pm\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2}})]. \end{eqnarray}$ (41)Using the relation ∇·x = 3, we compute ⟨ J ⟩ _i with the help of Eq. (9). Integrating over the cell volume, we obtain $\begin{matrix} ⟨ J ⟩_{i} & = & \frac{1}{3} \int_{Ω_{i}} \nabla_{ξ} \cdot (Nx) dΩ \\ = & \frac{1}{3} \sum_{d = 1}^{3} (⟨ N_{d} x ⟩_{i + \frac{1}{2} e^{d}} - ⟨ N_{d} x ⟩_{i - \frac{1}{2} e^{d}}) . \end{matrix}$ $\begin{eqnarray} \avg{J}_{\bs{i}}& =& \frac{1}{3}\int_{\Omega_{\bs{i}}}\bs{\nabla_\xi} \cdot (\mbf{Nx})\,\mrm{d}\Omega \nonumber\\ & =& \frac{1}{3}\sum_{d=1}^3 ( \avg{{\bf N}_d\mbf{x}}_{\bs{i}+\frac{1}{2}\bs{e}^d} -\avg{{\bf N}_d\mbf{x}}_{\bs{i}-\frac{1}{2}\bs{e}^d}). \end{eqnarray}$ (42)The face averages $⟨ N_{d} x ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{{\bf N}_d\mbf{x}}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ can be computed in a similar manner as the face average fluxes $⟨ N_{d} F^{s} ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{{\bf N}_d\bs{\mfk{F}}^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ in Eq. (35). With the help of Eq. (30) we write $\begin{matrix} ⟨ N_{d} x ⟩_{i \pm \frac{1}{2} e^{d}} & = & ⟨ N_{d}^{1} ⟩_{i \pm \frac{1}{2} e^{d}} ⟨ x ⟩_{i \pm \frac{1}{2} e^{d}} + \frac{1}{12} \begin{matrix} 3 \\ \sum \\ d^{'} = 1 d^{'} \neq d \end{matrix} \frac{\partial N_{d}^{1}}{\partial ξ_{d^{'}}} \frac{∂x}{\partial ξ_{d^{'}}} \\ + ⟨ N_{d}^{2} ⟩_{i \pm \frac{1}{2} e^{d}} ⟨ y ⟩_{i \pm \frac{1}{2} e^{d}} + \frac{1}{12} \begin{matrix} 3 \\ \sum \\ d^{'} = 1 d^{'} \neq d \end{matrix} \frac{\partial N_{d}^{2}}{\partial ξ_{d^{'}}} \frac{∂y}{\partial ξ_{d^{'}}} \end{matrix}$ $\begin{eqnarray} \avg{{\bf N}_d\mbf{x}}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} &=\ & \avg{N_d^1}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} \avg{x}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} + \frac{1}{12}\sum_{\substack{d^\prime=1\\d^\prime\ne d}}^3 \frac{\partial N_d^1}{\partial\xi_{d^\prime}} \frac{\partial x}{\partial\xi_{d^\prime}} \nonumber\\ &&+ \avg{N_d^2}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} \avg{y}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} + \frac{1}{12}\sum_{\substack{d^\prime=1\\d^\prime\ne d}}^3 \frac{\partial N_d^2}{\partial\xi_{d^\prime}} \frac{\partial y}{\partial\xi_{d^\prime}} \nonumber\\ &&+ \avg{N_d^3}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} \avg{z}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} + \frac{1}{12}\sum_{\substack{d^\prime=1\\d^\prime\ne d}}^3 \frac{\partial N_d^3}{\partial\xi_{d^\prime}} \frac{\partial z}{\partial\xi_{d^\prime}}, \label{eq:Nx} \end{eqnarray}$ (43)where the face averages $⟨ x ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{x}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ , $⟨ y ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{y}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ and $⟨ z ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{z}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ are given by the expressions $\begin{matrix} ⟨ y ⟩_{i \pm \frac{1}{2} e^{d}} = y_{i \pm \frac{1}{2} e^{d}} + \frac{1}{6} 𝒞_{i \pm \frac{1}{2}, ⊥}^{(2)} (y, \frac{1}{2}), \end{matrix}$ $\begin{eqnarray} &&\avg{x}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} = x_{\bs{i}\pm\frac{1}{2}\bs{e}^d}+ \frac{1}{6}\mathcal{C}^{(2)}_{\bs{i}\pm\frac{1}{2},\perp}\left(x,\frac{1}{2}\right), \label{eq:faceavgx} \\ && \avg{y}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} = y_{\bs{i}\pm\frac{1}{2}\bs{e}^d}+ \frac{1}{6}\mathcal{C}^{(2)}_{\bs{i}\pm\frac{1}{2},\perp}\left(y,\frac{1}{2}\right), \end{eqnarray}$ and $⟨ z ⟩_{i \pm \frac{1}{2} e^{d}} = z_{i \pm \frac{1}{2} e^{d}} + \frac{1}{6} 𝒞_{i \pm \frac{1}{2}, ⊥}^{(2)} (z, \frac{1}{2}),$ $\begin{equation} \avg{z}_{\bs{i}\pm\frac{1}{2}\bs{e}^d} = z_{\bs{i}\pm\frac{1}{2}\bs{e}^d}+ \frac{1}{6}\mathcal{C}^{(2)}_{\bs{i}\pm\frac{1}{2},\perp}\left(z,\frac{1}{2}\right), \end{equation}$ (46)respectively. The first-order derivatives of $N_{d}^{c}$ $\hbox{$N_d^c$}$ , x, y, and z appearing in the sums in Eq. (43) are computed as follows. To calculate $\frac{\partial N_{1}^{1}}{∂η}$ $\hbox{$\dfrac{\partial N_1^1}{\partial\eta}$}$ , for instance, we have $\frac{\partial N_{1}^{1}}{∂η} = \frac{∂y}{∂η} \frac{\partial}{∂η} (\frac{∂z}{∂ζ}) + \frac{∂z}{∂ζ} \frac{\partial^{2} y}{\partial η^{2}} - \frac{∂z}{∂η} \frac{\partial}{∂η} (\frac{∂y}{∂ζ}) - \frac{∂y}{∂ζ} \frac{\partial^{2} z}{\partial η^{2}} \cdot$ $\begin{equation} \frac{\partial N_1^1}{\partial\eta} = \frac{\partial y}{\partial\eta} \frac{\partial}{\partial\eta} \left(\frac{\partial z}{\partial\zeta} \right) + \frac{\partial z}{\partial\zeta}\frac{\partial^2 y} {\partial\eta^2} - \frac{\partial z}{\partial\eta} \frac{\partial}{\partial\eta} \left(\frac{\partial y}{\partial\zeta} \right) - \frac{\partial y}{\partial\zeta}\frac{\partial^2 z} {\partial\eta^2}\cdot \label{eq:dn11} \end{equation}$ (47)Approximating the first and second derivatives in this expression by second-order accurate finite differences, we obtain

$\begin{matrix} {(\frac{\partial N_{1}^{1}}{∂η})}_{i \pm \frac{1}{2},j,k} = \\ \times (y_{i \pm \frac{1}{2},j + \frac{1}{2},k} - y_{i \pm \frac{1}{2},j - \frac{1}{2},k}) \\ (z_{i \pm \frac{1}{2},j + \frac{1}{2},k + \frac{1}{2}} - z_{i \pm \frac{1}{2},j + \frac{1}{2},k - \frac{1}{2}} - z_{i \pm \frac{1}{2},j - \frac{1}{2},k + \frac{1}{2}} + z_{i \pm \frac{1}{2},j - \frac{1}{2},k - \frac{1}{2}}) \\ + 4 (z_{i \pm \frac{1}{2},j,k + \frac{1}{2}} - z_{i \pm \frac{1}{2},j,k - \frac{1}{2}}) (y_{i \pm \frac{1}{2},j + \frac{1}{2},k} - 2 y_{i \pm \frac{1}{2},j,k} + y_{i \pm \frac{1}{2},j - \frac{1}{2},k}) \\ - (z_{i \pm \frac{1}{2},j + \frac{1}{2},k} - z_{i \pm \frac{1}{2},j - \frac{1}{2},k}) \\ \times (y_{i \pm \frac{1}{2},j + \frac{1}{2},k + \frac{1}{2}} - y_{i \pm \frac{1}{2},j + \frac{1}{2},k - \frac{1}{2}} - y_{i \pm \frac{1}{2},j - \frac{1}{2},k + \frac{1}{2}} + y_{i \pm \frac{1}{2},j - \frac{1}{2},k - \frac{1}{2}}) \\ - 4 (y_{i \pm \frac{1}{2},j,k + \frac{1}{2}} - y_{i \pm \frac{1}{2},j,k - \frac{1}{2}}) (z_{i \pm \frac{1}{2},j + \frac{1}{2},k} - 2 z_{i \pm \frac{1}{2},j,k} + z_{i \pm \frac{1}{2},j - \frac{1}{2},k}) . \end{matrix}$ $\begin{eqnarray} &&\left(\frac{\partial N_1^1}{\partial\eta}\right)_{i\pm\frac{1}{2},j,k} = \nonumber \\ &&\qquad \times (y_{i\pm\frac{1}{2},j+\frac{1}{2},k} - y_{i\pm\frac{1}{2},j-\frac{1}{2},k}) \nonumber\\ &&\qquad ( z_{i\pm\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2}} - z_{i\pm\frac{1}{2},j+\frac{1}{2},k-\frac{1}{2}} - z_{i\pm\frac{1}{2},j-\frac{1}{2},k+\frac{1}{2}} + z_{i\pm\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2}} ) \nonumber\\ &&\qquad +4 ( z_{i\pm\frac{1}{2},j,k+\frac{1}{2}} - z_{i\pm\frac{1}{2},j,k-\frac{1}{2}} ) ( y_{i\pm\frac{1}{2},j+\frac{1}{2},k} -2y_{i\pm\frac{1}{2},j,k} + y_{i\pm\frac{1}{2},j-\frac{1}{2},k} ) \nonumber\\ &&\qquad - (z_{i\pm\frac{1}{2},j+\frac{1}{2},k} - z_{i\pm\frac{1}{2},j-\frac{1}{2},k}) \nonumber\\ && \qquad \times ( y_{i\pm\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2}} - y_{i\pm\frac{1}{2},j+\frac{1}{2},k-\frac{1}{2}} - y_{i\pm\frac{1}{2},j-\frac{1}{2},k+\frac{1}{2}} + y_{i\pm\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2}}) \nonumber\\ &&\qquad -4 (y_{i\pm\frac{1}{2},j,k+\frac{1}{2}} \!-\! y_{i\pm\frac{1}{2},j,k\!-\!\frac{1}{2}}) ( z_{i\pm\frac{1}{2},j\!+\!\frac{1}{2},k} -2z_{i\pm\frac{1}{2},j,k} \! +\! z_{i\pm\frac{1}{2},j-\frac{1}{2},k} ). \end{eqnarray}$ (48)The derivatives of x, y, and z are approximated by central differences, i.e. $, \frac{∂x}{∂η}$ $\hbox{$,\dfrac{\partial x}{\partial\eta}$}$ , for instance, is given by ${(\frac{∂x}{∂η})}_{i \pm \frac{1}{2},j,k} = x_{i \pm \frac{1}{2},j + \frac{1}{2},k} - x_{i \pm \frac{1}{2},j - \frac{1}{2},k} .$ $\begin{equation} \left( \frac{\partial x}{\partial\eta} \right)_{i\pm\frac{1}{2},j,k} = x_{i\pm\frac{1}{2},j+\frac{1}{2},k} - x_{i\pm\frac{1}{2},j-\frac{1}{2},k}. \label{eq:dxdeta} \end{equation}$ (49)When using Eqs. (44)–(49) to calculate the face averages, $⟨ N_{d} x ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{{\bf N}_d\mbf{x}}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ do not require any data from neighboring cells.

2.4. Evaluation of face average fluxes

To obtain the face average fluxes $⟨ F^{s} ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{F^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ , $⟨ G^{s} ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{G^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ , and $⟨ H^{s} ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{H^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ , we follow the procedure described in detail by (McCorquodale & Colella 2011), which we summarize in the following:

1.
Compute ⟨ U^s ⟩ _i from ⟨ JU^s ⟩ _i using Eq. (32) $⟨ U^{s} ⟩_{i} = \frac{1}{⟨ J ⟩_{i}} (⟨ J U^{s} ⟩_{i} - \frac{1}{48} 𝒞_{i}^{(1)} (⟨ J ⟩, ⟨ U^{s} ⟩, 1)) .$ $\begin{equation} \avg{U^s}_{\bs{i}}=\frac{1}{\avg{J}_{\bs{i}}}\left( \avg{JU^s}_{\bs{i}}-\frac{1}{48}\mathcal{C}^{(1)}_{\bs{i}} (\avg{J},\avg{U^s},1)\right). \label{eq:jutou1} \end{equation}$ (50)However, because ⟨ U^s ⟩ is not yet known, use instead $⟨ U^{s} ⟩_{i} = \frac{1}{⟨ J ⟩_{i}} (⟨ J U^{s} ⟩_{i} - \frac{1}{48} 𝒞_{i}^{(1)} (⟨ J ⟩, \frac{⟨ J U^{s} ⟩}{⟨ J ⟩}, 1)) .$ $\begin{equation} \avg{U^s}_{\bs{i}} = \frac{1}{\avg{J}_{\bs{i}}} \left( \avg{JU^s}_{\bs{i}}-\frac{1}{48}\mathcal{C}^{(1)}_{\bs{i}} (\avg{J},\frac{\avg{JU^s}}{\avg{J}},1) \right). \label{eq:jutou2} \end{equation}$ (51)
2.
Obtain the pointwise value $U_{i}^{s}$ $\hbox{$U^s_{\bs{i}}$}$ at the cell center using Eq. (25).
3.
Convert U_i to W_i using the EOS, where W = (ρ,u,v,w,p)^T is the vector of primitive variables.
4.
Calculate the cell average ⟨ W^s ⟩ _i from W_i using Eq. (24).
5.
Reconstruct the face average $⟨ W^{s} ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{W^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ from the cell average ⟨ W^s ⟩ _i using the expression $⟨ W^{s} ⟩_{i + \frac{1}{2} e^{d}} = \frac{7}{12} (⟨ W^{s} ⟩_{i} + ⟨ W^{s} ⟩_{i + e^{d}}) - \frac{1}{12} (⟨ W^{s} ⟩_{i - e^{d}} + ⟨ W^{s} ⟩_{i + 2 e^{d}}) .$ $\begin{equation} \avg{W^s}_{\bs{i}+\frac{1}{2}\bs{e}^d}= \frac{7}{12}(\avg{W^s}_{\bs{i}} +\avg{W^s}_{\bs{i}+\bs{e}^d}) -\frac{1}{12}(\avg{W^s}_{\bs{i}-\bs{e}^d} +\avg{W^s}_{\bs{i}+2\bs{e}^d}). \end{equation}$ (52)
6.
Apply a limiter to the interpolated value $⟨ W^{s} ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{W^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ as described in Sect. 2.4.1 in McCorquodale & Colella (2011). This limiter is a variant of the smooth extrema preserving limiter by Colella & Sekora (2008), which was proposed as an improvement to the limiter for the piecewise parabolic method (PPM) in Colella & Woodward (1984). The application of the limiter leads to two different values at cell faces, $⟨ W^{s} ⟩_{i \pm \frac{1}{2} e^{d}}^{L}$ $\hbox{$\avg{W^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}^L$}$ and $⟨ W^{s} ⟩_{i \pm \frac{1}{2} e^{d}}^{R}$ $\hbox{$\avg{W^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}^R$}$ , in regions where the interpolated profiles are not smooth.
7.
Solve a Riemann problem at each cell face with $⟨ W^{s} ⟩_{i \pm \frac{1}{2} e^{d}}^{L}$ $\hbox{$\avg{W^s}_{i\pm\frac{1}{2}e^d}^L$}$ and $⟨ W^{s} ⟩_{i \pm \frac{1}{2} e^{d}}^{R}$ $\hbox{$\avg{W^s}_{i\pm\frac{1}{2}e^d}^R$}$ as the left and right state, respectively. In Apsara we utilize the exact Riemann solver for real gases of Colella & Glaz (1985).
8.
Compute $W_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\mbf{W}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ at the face centers from the solution of the Riemann problem $⟨ W^{s} ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{W^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ using Eq. (29).
9.
Use $W_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\mbf{W}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ to compute $F_{i \pm \frac{1}{2} e^{d}}^{s}$ $\hbox{$\bs{\mfk{F}}^s_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ .
10.
Calculate $⟨ F^{s} ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{F^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ , $⟨ G^{s} ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{G^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ , and $⟨ H^{s} ⟩_{i \pm \frac{1}{2} e^{d}}$ $\hbox{$\avg{H^s}_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ from $F_{i \pm \frac{1}{2} e^{d}}^{s}$ $\hbox{$\bs{\mfk{F}}^s_{\bs{i}\pm\frac{1}{2}\bs{e}^d}$}$ using Eq. (28).

Fig. 1

Examples of 2D meshes generated using smooth mapping M₁ (left) and circular mapping M₂ (right).

3. Mapping functions

Apsara is capable of integrating the Euler equations on arbitrary curvilinear structured grids using the mapped grid technique. To perform the test problems that we present in Sect. 4 we used three different grid mappings. The computational space is discretized with an equidistant Cartesian mesh for each grid mapping, in which the grid spacing is unity in the three coordinate directions ξ, η, and ζ. The edge lengths of the computational domain are therefore $L_{ξ} = N_{ξ}, L_{η} = N_{η}, and L_{ζ} = N_{ζ},$ $\begin{equation} L_\xi=N_\xi,\; L_\eta=N_\eta,\; \mrm{and} \;L_\zeta=N_\zeta, \nonumber \end{equation}$ (53)where N_ξ,N_η, and N_ζ are the number of grid cells in the ξ, η, and ζ direction. We define normalized coordinates in computational space by $\begin{matrix} ˜ \\ ξ \end{matrix} = \frac{ξ}{N_{ξ}}, \begin{matrix} ˜ \\ η \end{matrix} = \frac{η}{N_{η}}, and \begin{matrix} ˜ \\ ζ \end{matrix} = \frac{ζ}{N_{ζ}},$ $\begin{equation} \tilde{\xi} = \frac{\xi }{N_\xi},\; \tilde{\eta} = \frac{\eta }{N_\eta},\;\mrm{and}\; \tilde{\zeta}= \frac{\zeta}{N_\zeta}, \nonumber \end{equation}$ (54)and we choose the coordinates of the inner (i.e., left) grid boundaries in computational space to be ξ_ib = η_ib = ζ_ib = 0.

3.1. Cartesian mapping M₀

In this simplest case, the coordinates (ξ,η,ζ) in computational space are mapped to physical space by $\begin{matrix} x = x_{ib} + (x_{ob} - x_{ib}) \cdot \begin{matrix} ˜ \\ ξ \end{matrix}, \\ y = y_{ib} + (y_{ob} - y_{ib}) \cdot \begin{matrix} ˜ \\ η \end{matrix}, \\ z = z_{ib} + (z_{ob} - z_{ib}) \cdot \begin{matrix} ˜ \\ ζ \end{matrix}, \end{matrix}$ $\begin{eqnarray} &&x=x_\mrm{ib}+(x_\mrm{ob}-x_\mrm{ib})\cdot\tilde{\xi},\\ &&y=y_\mrm{ib}+(y_\mrm{ob}-y_\mrm{ib})\cdot\tilde{\eta},\\ &&z=z_\mrm{ib}+(z_\mrm{ob}-z_\mrm{ib})\cdot\tilde{\zeta}, \end{eqnarray}$ where x_ib/ob,y_ib/ob, and z_ib/ob denote the coordinates of the inner/outer grid boundaries in physical space in the x,y, and z direction, respectively.

3.2. Smooth mapping M₁

For the test problems performed in a rectangular domain, we used a nonlinear mapping given in Colella et al. (2011). The mapping function is defined as $\begin{matrix} x = x_{ib} + (x_{ob} - x_{ib}) \cdot (\begin{matrix} ˜ \\ ξ \end{matrix} + c_{d} \cdot \sin 2 π \begin{matrix} ˜ \\ ξ \end{matrix} \sin 2 π \begin{matrix} ˜ \\ η \end{matrix} \sin 2 π \begin{matrix} ˜ \\ ζ \end{matrix}), \\ y = y_{ib} + (y_{ob} - y_{ib}) \cdot (\begin{matrix} ˜ \\ η \end{matrix} + c_{d} \cdot \sin 2 π \begin{matrix} ˜ \\ ξ \end{matrix} \sin 2 π \begin{matrix} ˜ \\ η \end{matrix} \sin 2 π \begin{matrix} ˜ \\ ζ \end{matrix}), \\ z = z_{ib} + (z_{ob} - z_{ib}) \cdot (\begin{matrix} ˜ \\ ζ \end{matrix} + c_{d} \cdot \sin 2 π \begin{matrix} ˜ \\ ξ \end{matrix} \sin 2 π \begin{matrix} ˜ \\ η \end{matrix} \sin 2 π \begin{matrix} ˜ \\ ζ \end{matrix}) . \end{matrix}$ $\begin{eqnarray} &&x=x_\mrm{ib}+(x_\mrm{ob}-x_\mrm{ib})\cdot(\tilde{\xi} +c_d\cdot\sin{2\pi\tilde{\xi}} \sin{2\pi\tilde{\eta}} \sin{2\pi\tilde{\zeta}}),\\ &&y=y_\mrm{ib}+(y_\mrm{ob}-y_\mrm{ib})\cdot(\tilde{\eta} +c_d\cdot\sin{2\pi\tilde{\xi}} \sin{2\pi\tilde{\eta}} \sin{2\pi\tilde{\zeta}}),\\ &&z=z_\mrm{ib}+(z_\mrm{ob}-z_\mrm{ib})\cdot(\tilde{\zeta} +c_d\cdot\sin{2\pi\tilde{\xi}} \sin{2\pi\tilde{\eta}} \sin{2\pi\tilde{\zeta}}). \end{eqnarray}$ This mapping function results in a smoothly deformed Cartesian mesh shown in the left panel of Fig. 1. The deformation is controlled by the parameter c_d. In all test problems, where we applied this smooth mapping, we set c_d equal to 0.1.

3.3. Circular mapping M₂

For the 2D test problems simulated in a circular domain, we used the singularity-free mapping function proposed by Calhoun et al. (2008), which is defined by $\begin{matrix} x = R_{D} \cdot \frac{D}{\begin{matrix} ˜ \\ r \end{matrix}} \cdot (2 \begin{matrix} ˜ \\ ξ \end{matrix} - 1), \\ y = R_{D} \cdot \frac{D}{\begin{matrix} ˜ \\ r \end{matrix}} \cdot (2 \begin{matrix} ˜ \\ η \end{matrix} - 1), \\ z = \begin{matrix} ˜ \\ ζ \end{matrix} \end{matrix}$ $\begin{eqnarray} &&x = R_D\cdot\frac{D}{\tilde{r}}\cdot(2\tilde{\xi}-1),\\ &&y = R_D\cdot\frac{D}{\tilde{r}}\cdot(2\tilde{\eta}-1),\\ &&z = \tilde{\zeta} \end{eqnarray}$ with $\begin{matrix} D = \begin{matrix} \max \end{matrix} (| 2 \begin{matrix} ˜ \\ ξ \end{matrix} - 1 |, | 2 \begin{matrix} ˜ \\ η \end{matrix} - 1 |), \\ \begin{matrix} ˜ \\ r \end{matrix} = \sqrt{(2 \begin{matrix} ˜ \\ ξ \end{matrix} - 1)^{2} + (2 \begin{matrix} ˜ \\ η \end{matrix} - 1)^{2}}, \end{matrix}$ $\begin{eqnarray} &&D = \max(\lvert2\tilde{\xi}-1\rvert,\lvert2\tilde{\eta}-1\rvert),\\ &&\tilde{r} = \sqrt{(2\tilde{\xi}-1)^2+(2\tilde{\eta}-1)^2}, \end{eqnarray}$ where R_D is the radius of the circular domain in physical space. This mapping function maps concentric squares in computational space to circular shells in physical space. As a result grid cells near the diagonals of the squares are severely distorted. An example of such a mesh is depicted in the right panel of Fig. 1.

4. Hydrodynamic tests

Table 1

L₁ and L_∞ errors of ⟨ ρ ⟩, and the convergence rates for the 1D linear advection test.

Table 2

Same as Table 1, but for the 2D linear advection test.

4.1. Linear advection

We performed one-dimensional (1D) and two-dimensional (2D) advection tests, where a Gaussian density profile $ρ (r) = 1 + e^{- 256 (r - \frac{1}{2})^{2}}$ $\begin{equation} \rho(r)=1+{\rm e}^{-256(r-\frac{1}{2})^2} \end{equation}$ (66)is advected with a constant velocity through a domain of size [0,1] (1D) and [0,1] × [0,1] (2D), respectively. In the 1D test r = x−x₀ and v = (1.0,0.0,0.0)^T, while r² = (x−x₀)² + (y−y₀)² and v = (1.0,0.5,0.0)^T in the 2D test. We placed the center of the Gaussian profile at (x₀,y₀) = (0.5,0.5) and imposed periodic boundary conditions in each coordinate direction.

While we employed a single grid with a constant spacing in the 1D test, we used two different meshes in the 2D tests: the Cartesian mesh M₀, and the smoothly deformed mesh M₁. We followed the advection of the Gaussian profile for two time units in 2D and for ten time units in 1D to compare our 1D results with those of McCorquodale & Colella (2011). The Gaussian profile returned to its original position at the end of each simulation. We set the CFL number to 0.2 and used the ideal gas EOS with $γ = \frac{5}{3}$ $\hbox{$\gamma=\frac{5}{3}$}$ .

We computed the L₁ and L_∞ norms of ⟨ ρ ⟩ from the simulation results, using ⟨ ρ ⟩ at time t = 0 as the reference solution. These norms, together with the corresponding convergence rates, are given in Tables 1 and 2 for the 1D and 2D simulations, respectively. In these tables we also give the results that we computed with the Prometheus code (Fryxell et al. 1991; Mueller et al. 1991).

For the 1D linear advection test as well as for the 2D test with M₀ and M₁, we find that Apsara shows 4th-order convergence according to both the L₁ and the L_∞ norm. The convergence rates for grid sizes below 256 zones in 1D and 256² zones in 2D are less than fourth order because at these grid resolutions the Gaussian profile is not well enough resolved. The L₁ and L_∞ errors are almost exactly the same in our 1D Apsara simulations as those obtained by McCorquodale & Colella (2011). The L₁ and L_∞ errors computed with Prometheus show convergence rates below third order. This results partly from the fact that the original PPM limiter implemented in Prometheus fails to preserve the extremum of the Gaussian profile (Colella & Sekora 2008). Because the limiter implemented in Apsara is similar to the extremum-preserving limiter proposed by Colella & Sekora (2008), clipping of extrema is considerably reduced in Apsara. We demonstrate this in Fig. 2, where the density profiles are plotted for 1D simulations that are performed with both codes on a grid of 128 grid cells.

Fig. 2

Density profiles for the 1D linear advection test obtained on a grid of 128 zones at t = 10. Results computed with Apsara and Prometheus are shown with blue crosses and red circles, respectively. The exact solution is shown in black.

4.2. Linear acoustic wave

In this test we simulated the propagation of a sound wave as described in Stone et al. (2008) in 1D, 2D, and 3D using the smooth mapping M₁ and the ideal gas EOS with $γ = \frac{5}{3}$ $\hbox{$\gamma=\frac{5}{3}$}$ . The background fluid has a density ρ₀ = 1, a pressure $p_{0} = \frac{3}{5}$ $\hbox{$p_0=\frac{3}{5}$}$ , and is initially at rest. We introduced a perturbation vector $δ U = A \sin (2 πx) \cdot (1, - 1, 1, 1, 1.5)^{T},$ $\begin{equation} \delta\mbf{U} = A \sin(2\pi x) \cdot (1,-1,1,1,1.5)^T, \end{equation}$ (67)which here only depends on the coordinate x, but in principle could depend on the other coordinates as well. We used an amplitude A = 10^-6 to perturb the conserved variables U₀ of the background fluid. The domain size was set equal to one wavelength of the pertubation, i.e., the size of the domain in each coordinate direction was unity.

Imposing periodic boundary conditions in all directions and using a CFL number of 1.3, we evolved the flow for one time unit so that the wave propagated once across the domain. In multi-D we also performed simulations using a second-order version of Apsara for comparison. This required an easy modification of the code, where one omits the second order correction terms when calculating face average fluxes. We used the RK4 scheme for time integration for both the fourth-order and second-order accurate code version. Reducing the order of the spatial discretization to second order was sufficient to downgrade the overall convergence from fourth order to second order, i.e., the accuracy of the code is dominated by spatial discretization errors in this test.

As described in Stone et al. (2008), we also computed the L₁ error for each component of the conserved variables U using the pointwise values at the cell centers. Fourth-order accurate approximations of these pointwise values of U^s were calculated as described in steps 1 and 2 in Sect. 2.4. Following these steps is essential when one wants to achieve fourth-order convergence. If one computes the L₁ errors without distinguishing between pointwise values and cell average values, one observes only second-order convergence even though the numerical scheme is of higher order. Table 3 gives the sum of the L₁ errors of all conserved variables $ϵ = \sum_{s} \sqrt{| L_{1} (U^{s}) |^{2}} .$ $\begin{equation} \epsilon = \sum_{\rm s}\sqrt{\lvert L_1({U}^s)\rvert^2}. \label{eq:linwav_epsiln} \end{equation}$ (68)This table also shows that the fourth-order implementation of Apsara achieves the design convergence rate only at low resolution, i.e., up to about 128 grid points per coordinate direction, while at higher resolutions the convergence rate drops dramatically. This is not unexpected because the reference solution used for computing the errors is only a solution of the linearized Euler equations. Hence, the convergence rate starts to degrade when the errors reach a level of A² ~ 10^-12 and nonlinear terms can no longer be neglected. On the other hand, second-order convergence is achieved with the second-order version of Apsara as expected, but here the convergence rate also degrades in 2D at high resolutions when ϵ ~ 10^-12.

To demonstrate that fourth-order convergence is indeed achieved with Apsara, we calculated the L₁ norm of the coarse-fine differences of the cell center values of ρ for different grid resolutions. More precisely, we first computed a fourth-order approximation of ρ at the cell centers of a coarser grid from the values given on the next finer grid, and then computed the differences between these cell center values and those obtained on the coarser grid. We refer to Appendix A for a more detailed explanation of this coarse-fine interpolation. The results shown in Table 4 indeed confirm the fourth-order convergence of the coarse-fine differences in ρ. In 2D, the convergence rate decreases slightly when comparing results from simulations with 512² and 1024² grid points owing to numerical round-off errors. As expected the second-order accurate implementation of Apsara yields second-order convergence.

Table 3

Sum of the mean errors of all conserved variables ϵ, defined in Eq. (68), and its convergence rate for the linear acoustic wave test.

Table 4

L₁-norm of the coarse-fine differences of the pointwise densities ρ as a function of grid resolution for the linear acoustic wave test.

4.3. Advection of a nonlinear vortex

In this 2D test problem we simulated the advection of an isentropic vortex with a uniform backgroud flow, following Yee et al. (2000). The background flow has a density, pressure, and temperature ρ₀ = p₀ = T₀ = 1, and moves with velocities u₀ = v₀ = 1. The isentropic vortex is added to the freestreaming background flow by introducing perturbations in velocity and temperature, which are given by $\begin{matrix} (δu,δv) = \frac{ϵ}{2 π} e^{\frac{1 - r^{2}}{2}} (- y̅, x̅), \\ δT = - \frac{(γ - 1) ϵ^{2}}{8 γ π^{2}} e^{\frac{1 - r^{2}}{2}}, \end{matrix}$ $\begin{eqnarray} &&(\delta u,\delta v) = \frac{\epsilon}{2\pi} {\rm e}^{\frac{1-r^2}{2}} (-\bar{y},\bar{x}),\\ &&\delta T = -\frac{(\gamma-1)\epsilon^2}{8\gamma\pi^2} {\rm e}^{\frac{1-r^2}{2}}, \end{eqnarray}$ where ϵ is a free parameter regulating the vortex strength, and $\hbox{$r^2=\bar{x}^2+\bar{y}^2$}$ . The coordinates $\hbox{$(\bar{x},\bar{y})$}$ are measured from the center of the vortex located at (x₀,y₀). We chose ϵ = 5, and used the ideal gas EOS with γ = 1.4 for this test. After adding the perturbation $\begin{matrix} (u,v) = (u_{0}, v_{0}) + (δu,δv), \\ T = T_{0} + δT, \end{matrix}$ $\begin{eqnarray} &&(u,v) = (u_0,v_0) + (\delta u,\delta v),\\ &&T = T_0 + \delta T, \end{eqnarray}$ all conserved quantities can be calculated using the fact that the flow is isentropic.

We performed this test on a 2D rectangular domain [−10,10] × [−10,10] using the mapping function M₁, and also on a 2D circular domain of radius R_D = 10 using the circular mapping function M₂. We assume a free outflow boundary in both coordinate directions. The vortex is initially placed at (x₀,y₀) = (−1,−1). We follow the advection of this vortex for two time units with a CFL number of 1.3. The analytic solution tells one that the center of the vortex is located at (x,y) = (1,1) at the end of the simulations. We computed the L₁ and L_∞ norms of the errors of the pointwise value of the density ρ at the cell center (see Table 5).

The results from the simulations performed on the M₁ grid show fourth-order convergence in both the L₁ and L_∞ norm. In contrast, the convergence order is reduced significantly on the M₂ grid. This illustrates that the accuracy and convergence rate of the fourth-order scheme implemented in Apsara is degraded when using a mapping function that is nonsmooth. Grimm-Strele et al. (2014) also performed the same test problem with the same M₂ mapping function, using the weighted essentially nonoscillatory (WENO) finite volume scheme. They also found a reduced convergence rate for the M₂ grid in comparison with the smooth mapping M₁. We speculate that the convergence rate on the grid M₂ is destroyed by the very poor mesh quality along the diagonals since the angle between two adjacent cell sides is nearly 180° for grid cells there.

Table 5

L₁ and L_∞ errors of ρ, and their convergence rates for the 2D advection of a nonlinear vortex test.

Table 6

L₁ and L_∞ errors and their convergence rates for the 2D advection of a nonlinear vortex on grids with four different smoothing widths b (see Sect. 4.4).

4.4. Influence of grid smoothness

In many astrophysical simulations in spherical geometry one usually employs logarithmic spacing for grid cells in the radial direction to cover large spatial length scales. For example, the grid spacing of a radial zone i, Δr_i, can be written as $Δ r_{i} = α Δ r_{i - 1} for 2 \leq i \leq N r$ $\begin{equation} \Delta r_i = \alpha \Delta r_{i-1}\;\;\textrm{for}\;\;2\le i\le N_r \end{equation}$ (73)where Δr_i−1 is the grid size of zone i−1, N_r is the number of radial grid zones, and α is the grid expansion ratio. If α is too large, the grid can be considered as nonsmooth and, as a result, the order of accuracy of the employed numerical scheme is degraded. Therefore, in this section we systematically study the influence of the smoothness of the mapping function on the performance of Apsara to determine approximately the maximum grid expansion ratio allowed in the case of a nonequidistant grid spacing without losing fourth-order convergence.

To this end, we considered the isentropic vortex test described in Sect. 4.3 on a 2D mesh with variable grid spacing in one coordinate direction. The 2D mesh used for this test is designed in the following way. Firstly, we defined two grid spacing parameters: δ₁ = 2L_x/N_ξ and δ₂ = 2L_x/ 3N_ξ, where L_x is the domain size in the x direction. The parameters δ₁ and δ₂ control the largest and smallest cell width, respectively. The cell width of the ith cell in the x direction is a linear combination of δ₁ and δ₂ , and is written as $δ x_{i} = (1 - w_{i}) \cdot δ_{1} + w_{i} \cdot δ 2$ $\begin{equation} \delta x_i = ( 1 - w_i ) \cdot \delta_1 + w_i \cdot \delta_2 \end{equation}$ (74)where w is a smoothing function of the form $w_{i} = {\begin{matrix} i \leq N_{ξ} / 2, \\ else . \end{matrix}$ $\begin{equation} w_i = \begin{cases} 0.5 \cdot \left( \tanh \left( 2 \pi\,b\,(i/N_\xi-1/8) \right) + 1 \right), & i \leq N_\xi/2, \\ 0.5 \cdot \left( \tanh \left( - 2 \pi\,b\,(i/N_\xi-7/8) \right) + 1 \right), & \text{else}. \end{cases} \label{eq:smoothfunc} \end{equation}$ (75)The smoothing width b regulates by how much the cell width varies in the x direction. In the limit of b → ∞, one obtains a nonsmooth mapping with a discontinuity in the grid spacing in the x direction. The cell spacing in the y direction is kept constant for simplicity and is δy = L_y/N_η where L_y is the domain size in the y direction.

Parameters for the numerical setup differ slightly from those described in Sect. 4.3. In this section, the computational domain is [0,15] × [0,15] , and the center of the vortex is placed at (x₀,y₀) = (7.5,7.5). The boundary condition is periodic in both coordinate directions. All other parameters have the same values as in Sect. 4.3. We simulated the advection of this vortex for 15 time units, i.e., the vortex is advected across the computational domain once and returns to its initial position at the end of the simulations.

In Table 6 we give the order of convergence in dependence of the smoothing width b computed both in the L₁ norm and the L_∞ norm. For small values of b, we find fourth-order convergence, whereas for large b the results only show second-order convergence. For intermediate values of b, the convergence order changes from 2 to 4 when the resolution is increased. For a more quantitative discussion of the grid smoothness, we define the maximum expansion ratio $c = \begin{matrix} \max \end{matrix} (\frac{Δ x_{i + 1}}{Δ x_{i}}, \frac{Δ x_{i - 1}}{Δ x_{i}}) \cdot$ $\begin{equation} c = \max \left( \frac{\Delta x_{i+1}}{\Delta x_i}, \frac{\Delta x_{i-1}}{\Delta x_i} \right)\cdot \end{equation}$ (76)As Table 7 shows there is a close correlation between the value of c and the L₁ convergence rate; the numerical solution converges with fourth order when c ≲ 1.5 and with second order otherwise.

Table 7

Dependence of the maximum expansion ratio of the grid c and the L₁ convergence rate on the grid resolution and smoothing width b.

In accordance with the conclusions of Visbal & Gaitonde (2002) and Bassi & Rebay (1997), we find that the actually achievable order of high-order hydrodynamics codes depends on the grid smoothness, which can be assessed, amongst others, by the expansion ratio c (see the 1D example above) or the cell orthogonality (in multi-d applications). The latter is very poor for the grids proposed by Calhoun et al. (2008) and, hence, explains the poor performance of Apsara on these grids.

4.5. Gresho vortex

We considered the Gresho vortex problem (Gresho & Chan 1990) to test the capability of Apsara to treat low-Mach number flows. The Gresho vortex is a rotating flow in which the centrifugal forces are balanced by the pressure gradients, resulting in a steady-state solution of a rotating vortex. Miczek et al. (2015) used the Gresho vortex test to demonstrate the inability of Godunov-type schemes relying on Roe’s flux to correctly simulate low-Mach number flows. We follow their numerical setup here.

We simulated the flow in a 2D periodic domain [0,1] × [0,1] using the mapping M₀. The fluid density ρ was set equal to 1 everywhere in the domain, and the vortex was placed at the center of the computational domain, i.e., the coordinates of its center were (x₀,y₀) = (0.5,0.5). The angular velocity of the rotating flow as a function of radius r measured from the center of the vortex is $\begin{matrix} u_{φ} (r) = 0.4 π {\begin{matrix} \end{matrix} \end{matrix}$ $\begin{eqnarray} u_\phi(r) = 0.4 \pi \begin{cases} 5r ;&\,0 \le r <0.2 \\ 2-5r ;&\,0.2\le r <0.4 \\ 0 ;&\,0.4\le r \end{cases} \end{eqnarray}$ (77)and the radial profile of the fluid pressure is $\begin{matrix} p (r) = {\begin{matrix} \end{matrix} \end{matrix}$ $\begin{eqnarray} p(r) = \begin{cases} p_0 + \frac{25}{2}r^2 ;&\,0\le r <0.2 \\ p_0 + \frac{25}{2}r^2 + 4(1-5r-\ln{0.2}+\ln{r}) ;&\,0.2\le r <0.4 \\ p_0 - 2 + 4\ln{2} ;&\,0.4\le r \end{cases} \label{eq:pgresho} \end{eqnarray}$ (78)with $p_{0} = \frac{(0.4 π)^{2}}{γ M_{\max}^{2}} - \frac{1}{2} \cdot$ $\begin{equation} p_0 = \frac{(0.4\pi)^2}{\gamma M_\mrm{max}^2}-\frac{1}{2} \cdot \end{equation}$ (79)The maximum Mach number M_max of the rotating flow is attained at a radius r = 0.2. We used the ideal gas EOS with $γ = \frac{5}{3}$ $\hbox{$\gamma = \frac{5}{3}$}$ and evolved the vortex until time t = 1, i.e., until the vortex completed one revolution.

Fig. 3

Mach number distributions normalized by the maximum Mach number M_max at t = 1 for the Gresho vortex test. The top and bottom rows show results computed by Apsara and Prometheus, respectively. Regardless of M_max the vortex remains intact when simulated with Apsara, while with Prometheus the vortex gradually dissolves because of numerical dissipation as M_max decreases.

Fig. 4

Total kinetic energy as a function of time E_kin(t) normalized by the total kinetic of the initial state E_kin(0) for the Gresho vortex problem. The solid and dashed lines show the results computed with Apsara and Prometheus, respectively. Curves of different colors correspond to different maximum Mach numbers. In the case of Apsara, these curves are visually indistinguishable.

Figure 3 shows the Mach number distributions at the final time t = 1 for decreasing values of M_max ranging from 10^-1 to 10^-4. The results obtained with Apsara for different values of M_max are almost identical, i.e., they are independent of M_max. This behavior differs drastically from that found with Prometheus, for which the vortex gradually dissolves as M_max decreases, and it eventually even disappears at M_max = 10^-4. The different performance of the two codes is also evident when one considers the time evolution of the total kinetic energy of the vortex, which is shown in Fig. 4 normalized by its initial value at t = 0. In the simulations performed with Apsara the total kinetic energy had decreased by less than 0.5% at t = 1 for all M_max, while with Prometheus the vortex lost an increasingly larger amount of kinetic energy as M_max was decreased. For M_max = 10^-4, more than half of the total kinetic energy was gone at the end of the simulation performed with Prometheus.

Fig. 5

Mean turbulent kinetic energy, mean dissipation rate of turbulent kinetic energy, mean enstrophy, and numerical Reynolds number (from top to bottom) plotted as function of the dimensionless time t^∗ for the Taylor-Green vortex test for the reference Mach number M_ref ~ 0.29.

Fig. 6

Same as Fig. 5, but for M_ref = 10^-2.

4.6. Taylor-Green vortex

We quantify the numerical viscosity of Apsara in the low-Mach number limit using the Taylor-Green vortex (Taylor & Green 1937) and compare it to Prometheus. Thereby, we follow the numerical setup described in Drikakis et al. (2007) and Edelmann (2014). The initial conditions were $\begin{matrix} ρ & = & ρ_{0}, \\ u & = & u_{0} \sin kx \cdot \cos ky \cdot \cos kz, \\ v & = & - u_{0} \cos kx \cdot \sin ky \cdot \cos kz, \\ w & = & 0, \\ p & = & p_{0} + (\frac{{u_{0}}^{2} ρ}{16}) (2 + \cos {2 kz}^{)} (\cos 2 kx + \cos 2 ky), \end{matrix}$ $\begin{eqnarray} \rho & =& \rho_0, \\ u & =& u_0 \sin{kx} \cdot \cos{ky} \cdot \cos{kz}, \\ v & =& - u_0 \cos{kx} \cdot \sin{ky} \cdot \cos{kz}, \\ w & =& 0, \\ p & =& p_0 + \left(\frac{{u_0}^2\rho}{16}\right) \left(2+\cos{2kz} \right) \left(\cos{2kx} + \cos{2ky} \right), \end{eqnarray}$ where ρ₀ = 1.178 × 10^-3 g cm^-3, p₀ = 10⁶ bar, and k = 0.01 cm^-1. The velocity u₀ was set by selecting a reference Mach number of the flow M_ref = u₀/c₀, where $c_{0} = \sqrt{γ p_{0} / ρ_{0}}$ $\hbox{$c_0=\sqrt{\gamma p_0/\rho_0}$}$ is the sound speed. We considered two cases: M_ref ~ 0.29 (i.e.u₀ = 10⁴ cm s^-1) and M_ref = 10^-2. For both cases we used an ideal gas EOS with γ = 1.4. The vortex was simulated in a periodic cubic domain with an edge length of 2π·100 cm using the Cartesian mapping M₀. The simulations were carried out with 64³ and 128³ grid cells until a dimensionless time t^∗ = ku₀t = 30.

Figures 5 and 6 show a comparison of the results obtained with Apsara and Prometheus for simulations with M_ref ~ 0.29 and M_ref = 10^-2, respectively. In both figures the top panel shows the time evolution of the nondimensionalized mean turbulent kinetic energy K^∗ defined by $K^{*} = \frac{K}{u_{o}^{2}} = \frac{1}{2 u_{0}^{2}} (u - u)^{2},$ $\begin{equation} K^*=\frac{K}{u_o^2}=\frac{1}{2u_0^2}\overline{(\bs{u}-\overline{\bs{u}})^2}, \end{equation}$ (85)where u = (u,v,w) is the velocity vector and the overbar denotes a volumetric average. As large-scale eddies decay to smaller ones, kinetic energy is dissipated by numerical viscosity. The corresponding energy dissipation rates are plotted in the second panels (from top) of Figs. 5 and 6, which show that the dissipation rate peaks roughly at the same time as the mean enstrophy, $Ω^{*} = \frac{Ω}{(k u_{0})^{2}} = \frac{1}{2 (k u_{o})^{2}} {| \nabla \times u |}^{2},$ $\begin{equation} \Omega^*=\frac{\Omega}{(ku_0)^2} =\frac{1}{2(ku_o)^2}\overline{\left|\bs{\nabla}\times\bs{u}\right|^2} , \end{equation}$ (86)given in the third panels of Figs. 5 and 6. From the enstrophy Ω^∗ and the dissipation rate dK^∗/ dt^∗ we calculate the numerical Reynolds number Re using the relation $Re = - \frac{Ω^{*}}{d K^{*} / d t^{*}}$ $\begin{equation} Re=-\frac{\Omega^*}{\mrm{d}K^*/\mrm{d}t^*} \end{equation}$ (87)for an incompressible Navier-Stokes fluid (Drikakis et al. 2007). The numerical Reynolds number are plotted in the bottom panels of Figs. 5 and 6.

Although the results computed with Apsara and Prometheus do not differ markedly for M_ref ~ 0.29 (Fig. 5), small differences can be recognized. The mean kinetic energy dissipates slightly faster near its peak for Prometheus than for Apsara for the same grid resolution. On the other hand, the mean enstrophy is slightly lower for Prometheus than for Apsara. As a consequence, the numerical Reynolds number is slightly higher for Apsara than for Prometheus, implying a lower numerical viscosity of the former code. This difference becomes more evident for M_ref = 10^-2 (Fig. 6); for a grid resolution of 128³ zones the Reynolds number is about 1.5 times larger for Apsara than for Prometheus.

Table 8

Parameters used in our simulations to generate a turbulent flow (for more details see text).

Fig. 7

Square root of the magnitude of the vorticity $\sqrt{| \nabla \times u |}$ $\hbox{$\sqrt{\left|\bs{\nabla}\times\bs{u}\right|}$}$ in the x-z plane along y = 0 from the turbulence simulation with 512³ grid cells at time t = 100.

5. Turbulence in the context of core collapse supernovae

Lately, the study of the properties of turbulent flows in CCSN has received great attention among supernova researchers because turbulent pressure is thought to provide additional support for the revival of the stalled supernova shock wave, thus aiding an explosion (Murphy et al. 2013; Couch & Ott 2015). Abdikamalov et al. (2015) investigated the resolving power of 3D CCSN simulations. They concluded that the grid resolutions presently used in 3D implicit large eddy simulations (ILES) of CCSN are still at least ~7−8 times too low to correctly model the turbulent energy cascade down to the inertial range. Instead, the excessive numerical viscosity resulting from a grid resolution that is too low causes inefficient transport of kinetic energy from large to small scales, creating a pile-up of energy between the energy injection scale and the dissipation range. Their conclusion concerning a grid resolution that is adequate for resolving the inertial range is based on earlier work on isotropic, stationary turbulence by Sytine et al. (2000). However, because turbulence in CCSN is driven by time-dependent neutrino heating of post-shock matter, the resulting turbulent flow is anisotropic and nonstationary.

To study turbulence in the regime relevant to CCSN Radice et al. (2015) performed simulations of anisotropic turbulence with a characteristic Mach number of ~0.3 in a 3D Cartesian box using the Flash code (Fryxell et al. 2000; Dubey et al. 2009). In their simulations, the stirring force is enhanced in the x direction, which they used to mimic the preferred radial direction in CCSN, such that the turbulent kinetic energy is twice larger in the x direction than in the other two coordinate directions. They considered five different HRSC methods varying the spatial reconstruction order, Riemann solver, and grid resolutions. Their findings agree with those of the previous study by Sytine et al. (2000) of isotropic turbulence. At low resolution, the energy cascade is strongly affected by the bottleneck effect. Radice et al. (2015) recovered the inertial range only partially even at the highest resolution (512 grid points per coordinate direction) considered. These authors estimated that at this resolution about 20% of the energy is still accumulated at intermediate scales using the least dissipative scheme considered in their work.

Motivated by the study of Radice et al. (2015) we performed a set of ILES of anisotropic turbulence in the context of CCSN with Apsara. Our numerical setup resembles that of Radice et al. (2015) and can be summarized as follows. Turbulence is driven by an external acceleration field by adding a stirring acceleration term to the RHS of the momentum equation. We closely followed the implementation of the stirring module in Flash and list the respective parameters in Table 8. We generated an acceleration field that was purely solenoidal (divergence-free), only containing power within a small range of Fourier modes between wavenumbers k_min and k_max. Six separate phases are evolved for each stirring mode after each time step Δt in Fourier space by an Ornstein-Uhlenbeck (OU) random process (Uhlenbeck & Ornstein 1930). The phases ψ are updated as $ψ (t + Δ t) = f ψ (t) + σ \sqrt{1 - f^{2}} ℛ_{g},$ $\begin{equation} \psi(t+\Delta t)=f\,\psi(t)+\sigma\sqrt{1-f^2}\,\mathcal{R}_g , \end{equation}$ (88)where f = exp(−Δt/τ) is a decay factor with correlation timescale τ, and σ is the variance of the OU process. The value ℛ_g is a Gaussian random variable drawn from a Gaussian distribution with unit variance. The power per mode is set by the variance σ and is chosen to be constant for all excited wavenumbers. To break isotropy, we multiplied the OU variance σ in x direction by a factor of four before carrying out the solenoidal projection in Fourier space. We tuned the OU variance σ such that the flow has a rms Mach number of ≈0.37.

We performed the simulations in a unit cube using the Cartesian mapping M₀ employing four grid resolutions of 64³, 128³, 256³, and 512³ zones, respectively. We simulated the evolution of the flow until t = 100. Figure 7 shows a snapshot at the final time of the square root of the magnitude of the vorticity $\sqrt{| \nabla \times u |}$ $\hbox{$\sqrt{\left|\bs{\nabla}\times\bs{u}\right|}$}$ in a 2D (x-z) slice through the middle of the simulation box (y = 0) for the simulation with 512³ zones. In Fig. 8 we show for all grid resolutions the compensated 1D energy spectra E(k) k^5/3, which we computed following Radice et al. (2015).

Fig. 8

Energy spectra (Eq. (89)) plotted vs. the dimensionless wavenumber 512kΔx for four different grid resolutions. The spectra are compensated by a k^5/3 spectrum, so that any region with Kolmogorov scaling appears flat.

Fig. 9

Effective numerical viscosity (Eq. (91)) as a function of wavenumber estimated from our ILES simulations for four different grid resolutions.

We computed first the 3D energy spectrum $E (k) = \frac{1}{2} \hat{u} \cdot \hat{u} *,$ $\begin{equation} E(\bs{k}) = \frac{1}{2} \hat{\bs{u}} \cdot \hat{\bs{u}}^*, \label{eq:energy-spectrum} \end{equation}$ (89)where $\hat{u}$ $\hbox{$\hat{\bs{u}}$}$ is the Fourier transform of the velocity field and $\hat{u}^{*}$ $\hbox{$\hat{\bs{u}}^*$}$ is its complex conjugate. To obtain the 1D energy spectrum E(k), we averaged the 3D spectrum over spherical shells in k space (Eswaran & Pope 1988) $E (k) = \frac{4 π k^{2}}{N_{k}} \sum_{k - \frac{1}{2} < | k | \leq k + \frac{1}{2}} E (k),$ $\begin{equation} E(k) = \frac{4\pi k^2}{N_k} \sum_{k-\frac{1}{2}<\left|\bs{k}\right|\le k+\frac{1}{2}} E(\bs{k}), \label{eq:kshellav} \end{equation}$ (90)where N_k denotes the number of discrete modes within the bin $k - \frac{1}{2} < | k | \leq k + \frac{1}{2}$ $\hbox{$k-\frac{1}{2} < \left|\bs{k}\right| \le k+\frac{1}{2}$}$ . We then averaged E(k) over time using 381 snapshots in the interval 5 ≤ t ≤ 100 spaced by 0.25 time units.

The resulting energy spectra are in very good agreement with those of the least dissipative numerical schemes considered by Radice et al. (2015), i.e., the third-order PPM reconstruction (Colella & Woodward 1984) and the improved fifth-order WENO reconstruction (WENO-Z; Borges et al. 2008) using the HLLC approximate Riemann solver of Toro et al. (1994). The compensated spectrum exhibits a flat region with Kolmogorov scaling (for 5 ≲ k ≲ 10) expected in the inertial range only at a resolution of 512³. At lower grid resolutions the compensated energy spectra are dominated by the bottleneck effect, which manifests itself with a ~k^-1 scaling.

On the other hand, the effective viscosity measured from our simulations differs from that reported in Radice et al. (2015). It is estimated according to $ν (k) = - \frac{1}{2} \frac{R (k)}{k^{2} E (k)},$ $\begin{equation} \nu(k )= -\frac{1}{2} \frac{R(k)}{k^2E(k)}, \label{eq:numvis} \end{equation}$ (91)where $R (k) = - T (k) - C (k) - ϵ (k)$ $\begin{equation} R(k) = -T(k) -C(k) -\epsilon(k) \end{equation}$ (92)is the residual of the 1D energy balance equation. The 1D energy transfer term T(k),C(k), and the energy injection rate ϵ(k) are averages of their 3D counterparts T(k),C(k), and ϵ(k), and are obtained in the same way we calculated the 1D energy spectrum E(k) (Eq. (90)). The 3D counterparts are defined as $\begin{matrix} ϵ (k) = ℜ [\hat{a} \cdot \hat{u} *], \end{matrix}$ $\begin{eqnarray} &&T(\bs{k}) = 2\pi\, \Re[(\hat{\bs{u}}*i\bs{k} \otimes \hat{\bs{u}})\cdot\hat{\bs{u}}^*], \label{eq:Tterm}\\ &&C(\bs{k}) = 2\pi\, \Re[(\frac{1}{\rho}*i\bs{k}p)\cdot\hat{\bs{u}}^*], \label{eq:Cterm}\\ &&\epsilon(\bs{k}) = \Re[\hat{\bs{a}}\cdot\hat{\bs{u}}^*], \end{eqnarray}$ where $\hat{a}$ $\hbox{$\hat{\bs{a}}$}$ is the Fourier transform of the acceleration field and ∗ denotes the convolution operator. In Fig. 9 we show the effective numerical viscosity ν(k) as a function of wavenumber estimated from our simulations. The effective viscosity exhibits fewer variations in k space than that estimated by Radice et al. (2015) from their simulations (see their Fig. 5).

6. Summary

We have developed a new numerical code for simulating astrophysical flows called Apsara. The code is based on the fourth-order implementation of a high-order, finite-volume method for mapped coordinates proposed by Colella et al. (2011). An extension of the method for solving nonlinear hyperbolic equations following McCorquodale & Colella (2011) and Guzik et al. (2012) is applied to solve the Euler equations of Newtonian gas dynamics. Apsara comprises great flexibility concerning grid geometry because of the implemented mapped-grid technique, which makes the code suitable for a wide range of applications.

By defining a mapping function, which describes the coordinate transformation between physical space and an abstract computational space, the governing equations written in Cartesian coordinates in physical space are transformed into equations of similar form in the computational space. The physical domain of interest can then be discretized as a general structured curvilinear mesh, while the computational space is still discretized by an equidistant Cartesian mesh that allows for easy and efficient implementation of numerical algorithms.

We strictly follow the implementation described in Colella et al. (2011) for computation of metric terms on faces of control volumes to ensure that Apsara preserves the freestream condition. The freestream property is crucial for grid-based codes in curvilinear coordinates because violation of this condition can introduce numerical artifacts as shown, for example, in Grimm-Strele et al. (2014) for the WENO-G finite difference scheme (Nonomura et al. 2010; Shu 2003).

The method of Colella et al. (2011) can be extended to a higher order of accuracy in time as described in Buchmüller & Helzel (2014). However, RK schemes of order higher than four are considerably more expensive because many more integration stages are needed in this case (Ketcheson 2008). Concerning fourth-order accurate RK schemes, the ten-stage algorithm SSPRK(10, 4) of Ketcheson (2008) has the advantage of possessing the strong stability-preserving (SSP) property, which is not the case for the classical RK4 scheme. In this work, however, we always used the RK4 scheme and the extension to the SSPRK(10, 4) scheme is subject of the future development of Apsara.

We validated our numerical code by simulating a set of hydrodynamic tests problems. In the case of smooth solutions, i.e., for the linear advection of a Gaussian profile, the propagation of a linear acoustic wave, and the advection of a nonlinear vortex, Apsara exhibits fourth-order accuracy both on a Cartesian mesh and the sinusoidally deformed mesh considered by Colella et al. (2011). However, our results from the advection of a nonlinear vortex test also revealed that the order of accuracy is severely degraded when using the circular mapping proposed by Calhoun et al. (2008), i.e., the grid smoothness directly impacts the achievable order of accuracy of the scheme. Our finding agrees with results obtained by Lemoine & Ou (2014), who used a modified version of grid mappings suggested by Calhoun et al. (2008). Lemoine & Ou also achieved roughly first-order convergence on these nonsmooth grid mappings. Thus, we conclude that the singularity-free mappings for circular and spherical domains of Calhoun et al. (2008) are not suitable for use in conjuction with high-order, finite-volume methods.

We quantified how the grid smoothness influences the order of accuracy of the high-order method of Colella et al. (2011) using the advection of a nonlinear vortex across a 2D mesh varying the nonequidistant grid spacing in one coordinate direction systematically. We found that the grid expansion ratio, i.e., the ratio of the cell size between two neighboring zones, must be less than ~1.5 in order to preserve fourth-order accuracy. This result can be regarded as a guideline for the grid setup in the case of a logarithmically spaced radial grid that is often employed in astrophysical simulations in spherical geometry.

Motivated by the work of Miczek et al. (2015), which demonstrated a failure of a Godunov-type scheme employing a Roe flux in the low-Mach number regime, we performed simulations of the Gresho vortex test varying the maximum Mach number of the vortex from 10^-1 to 10^-4. Our results show that the decrease of the kinetic energy of the vortex is independent of the maximum Mach number when computed with Apsara. In all our simulations the kinetic energy had decreased by less than ≈0.5% at the end of the simulations, i.e., after one complete rotation of the vortex. This contrasts with solutions calculated with Prometheus, which uses the dimensionally-split PPM method as do many other codes in astrophysics. The kinetic energy of the vortex decreased more and more rapidly when the maximum Mach number was decreased. In the lowest Mach number case, the Gresho vortex vanished completely and the kinetic energy was reduced by more than half of its initial value because of the numerical dissipation present in Prometheus.

The superior performance of the high-order method of McCorquodale & Colella (2011) for highly subsonic flows is further supported by our results obtained for the Taylor-Green vortex problem. We simulated this test for two reference Mach numbers and compared the results computed with Apsara and Prometheus. For a reference Mach number 0.29 we obtained similar results for both codes, whereas, for a reference Mach number 10^-2, Apsara achieved a larger numerical Reynolds number for the same resolution, implying a lower numerical viscosity than Prometheus.

To demonstrate Apsara’s performance on an astrophysical application, we performed ILES of anisotropic turbulence in a periodic Cartesian domain following the numerical setup used by Radice et al. (2015). The characteristic Mach number of the flow was 0.37, and the acceleration field driving the turbulence was enhanced along one coordinate direction to mimic the turbulent flow conditions in a CCSN. We obtained very similar results as those reported by Radice et al. (2015). At low resolution the energy spectra are dominated by the bottleneck effect at intermediate scales. We only began to recover the inertial range for a narrow range of wavenumbers at the highest grid resolution of 512 zones per coordinate direction. On the other hand, the effective viscosity as a function of wavenumber displayed fewer variations than that of Radice et al. (2015).

For problems to be simulated in a spherical domain, which we are particularly interested in, the singularity-free grid mappings of Calhoun et al. (2008) result in a loss of convergence order. Therefore, we plan to implement the mapped multiblock grid technique into Apsara following the strategy of McCorquodale et al. (2015), who extended the high-order, finite-volume method of Colella et al. (2011) to include this technique. They described an algorithm for data communication between blocks such that high-order accuracy is maintained. We also plan to implement additional physics into Apsara, for example, self-gravity, a nuclear EOS, and a nuclear reaction network, to tackle multiphysics astrophysical problems such as CCSN and stellar convection.

Acknowledgments

Most of the computations were performed at the Max Planck Computing & Data Facility (MPCDF). A.W. thanks Tomoya Takiwaki for fruitful discussions on turbulence in CCSN and H.G.S. acknowledges financial support from the Austrian Science fund (FWF), project P21742.

References

Abdikamalov, E., Ott, C. D., Radice, D., et al. 2015, ApJ, 808, 70 [NASA ADS] [CrossRef] [Google Scholar]
Appel, W. 2007, Mathematics for Physics and Physicists (Princeton University Press) [Google Scholar]
Bassi, F., & Rebay, S. 1997, J. Comput. Phys., 138, 251 [Google Scholar]
Borges, R., Carmona, M., Costa, B., & Don, W. S. 2008, J. Comput. Phys., 227, 3191 [Google Scholar]
Buchmüller, P., & Helzel, C. 2014, J. Sci. Comput., 61, 343 [CrossRef] [Google Scholar]
Calhoun, D. A., Helzel, C., & Leveque, R. J. 2008, SIAM Rev., 50, 723 [CrossRef] [Google Scholar]
Cardall, C. Y., Budiardja, R. D., Endeve, E., & Mezzacappa, A. 2014, ApJS, 210, 17 [NASA ADS] [CrossRef] [Google Scholar]
Colella, P., & Glaz, H. M. 1985, J. Comput. Phys., 59, 264 [NASA ADS] [CrossRef] [Google Scholar]
Colella, P., & Sekora, M. D. 2008, J. Comput. Phys., 227, 7069 [NASA ADS] [CrossRef] [Google Scholar]
Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
Colella, P., Dorr, M. R., Hittinger, J. A. F., & Martin, D. F. 2011, J. Comput. Phys., 230, 2952 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Couch, S. M., & Ott, C. D. 2015, ApJ, 799, 5 [NASA ADS] [CrossRef] [Google Scholar]
Drikakis, D., Fureby, C., Grinstein, F. F., & Youngs, D. 2007, J. Turbulence, 8, N20 [CrossRef] [Google Scholar]
Dubey, A., Antypas, K., Ganapathy, M. K., et al. 2009, Parallel Computing, 35, 512 [CrossRef] [Google Scholar]
Edelmann, P. 2014, Dissertation, Technische Universität München, München [Google Scholar]
Eswaran, V., & Pope, S. B. 1988, Computers and Fluids, 16, 257 [NASA ADS] [CrossRef] [Google Scholar]
Feng, X., Jiang, C., Xiang, C., Zhao, X., & Wu, S. T. 2012, ApJ, 758, 62 [Google Scholar]
Fragile, P. C., Lindner, C. C., Anninos, P., & Salmonson, J. D. 2009, ApJ, 691, 482 [NASA ADS] [CrossRef] [Google Scholar]
Fryxell, B., Arnett, D., & Mueller, E. 1991, ApJ, 367, 619 [NASA ADS] [CrossRef] [Google Scholar]
Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273 [NASA ADS] [CrossRef] [Google Scholar]
Gresho, P. M., & Chan, S. T. 1990, Int. J. Numer. Meth. Fluids, 11, 621 [Google Scholar]
Grimm-Strele, H., Kupka, F., & Muthsam, H. J. 2014, Comput. Phys. Comm., 185, 764 [Google Scholar]
Guzik, S., McCorquodale, P., & Colella, P. 2012, in 50th AIAA Aerospace Sciences Meetings (American Institute of Aeronautics and Astronautics, Nashville, TN), AIAA Paper 2012–0574 [Google Scholar]
Hayashi, H., & Kageyama, A. 2016, J. Comput. Phys., 305, 895 [NASA ADS] [CrossRef] [Google Scholar]
Hotta, H., Rempel, M., & Yokoyama, T. 2014, ApJ, 786, 24 [NASA ADS] [CrossRef] [Google Scholar]
Jiang, C., Feng, X., & Xiang, C. 2012, ApJ, 755, 62 [NASA ADS] [CrossRef] [Google Scholar]
Kageyama, A., & Sato, T. 2004, Geophys., Geosys., 5, 9005 [Google Scholar]
Ketcheson, D. I. 2008, SIAM J. Sci. Comput., 30, 2113 [CrossRef] [Google Scholar]
Kifonidis, K., & Müller, E. 2012, A&A, 544, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Koldoba, A. V., Romanova, M. M., Ustyugova, G. V., & Lovelace, R. V. E. 2002, ApJ, 576, L53 [NASA ADS] [CrossRef] [Google Scholar]
Lemoine, G. I., & Ou, M. Y. 2014, SIAM J. Sci. Comput., 36, B396 [CrossRef] [Google Scholar]
Lentz, E. J., Bruenn, S. W., Hix, W. R., et al. 2015, ApJ, 807, L31 [NASA ADS] [CrossRef] [Google Scholar]
LeVeque, R. 2002, Finite Volume Methods for Hyperbolic Problems, Cambridge Texts in Applied Mathematics (Cambridge University Press) [Google Scholar]
McCorquodale, P., & Colella, P. 2011, CAMCOS, 6, 1 [CrossRef] [Google Scholar]
McCorquodale, P., Dorr, M. R., Hittinger, J. A. F., & Colella, P. 2015, J. Comput. Phys., 288, 181 [NASA ADS] [CrossRef] [Google Scholar]
Melson, T., Janka, H.-T., Bollig, R., et al. 2015a, ApJ, 808, L42 [NASA ADS] [CrossRef] [Google Scholar]
Melson, T., Janka, H.-T., & Marek, A. 2015b, ApJ, 801, L24 [NASA ADS] [CrossRef] [Google Scholar]
Miczek, F., Röpke, F. K., & Edelmann, P. V. F. 2015, A&A, 576, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Mueller, E., Fryxell, B., & Arnett, D. 1991, A&A, 251, 505 [NASA ADS] [Google Scholar]
Müller, B. 2015, MNRAS, 453, 287 [NASA ADS] [CrossRef] [Google Scholar]
Murphy, J. W., Dolence, J. C., & Burrows, A. 2013, ApJ, 771, 52 [NASA ADS] [CrossRef] [Google Scholar]
Nonomura, T., Iizuka, N., & Fujii, K. 2010, Computers & Fluids, 39, 197 [CrossRef] [Google Scholar]
Peng, X., Xiao, F., & Takahashi, K. 2006, Quart. J. Roy. Meteor. Soc., 132, 979 [NASA ADS] [CrossRef] [Google Scholar]
Radice, D., Couch, S. M., & Ott, C. D. 2015, Comput. Astrophys. Cosmol., 2, 7 [NASA ADS] [CrossRef] [Google Scholar]
Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., & Lovelace, R. V. E. 2012, MNRAS, 421, 63 [NASA ADS] [Google Scholar]
Ronchi, C., Iacono, R., & Paolucci, P. S. 1996, J. Comput. Phys., 124, 93 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
Shiota, D., Kusano, K., Miyoshi, T., & Shibata, K. 2010, ApJ, 718, 1305 [NASA ADS] [CrossRef] [Google Scholar]
Shu, C.-W. 2003, Int. J. Comput. Fluid Dyn., 17, 107 [NASA ADS] [CrossRef] [Google Scholar]
Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137 [NASA ADS] [CrossRef] [Google Scholar]
Sytine, I. V., Porter, D. H., Woodward, P. R., Hodson, S. W., & Winkler, K.-H. 2000, J. Comput. Phys., 158, 225 [NASA ADS] [CrossRef] [Google Scholar]
Taylor, G. I., & Green, A. E. 1937, Proc. Roy. Soc. London A, 158, 499 [Google Scholar]
Toro, E. F., Spruce, M., & Speares, W. 1994, Shock Waves, 4, 25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Uhlenbeck, G. E., & Ornstein, L. S. 1930, Phys. Rev., 36, 823 [NASA ADS] [CrossRef] [Google Scholar]
Visbal, M. R., & Gaitonde, D. V. 2002, J. Comput. Phys., 181, 155 [NASA ADS] [CrossRef] [Google Scholar]
Warren, D. C., & Blondin, J. M. 2013, MNRAS, 429, 3099 [NASA ADS] [CrossRef] [Google Scholar]
Wilson, J. R., & Mayle, R. W. 1988, Phys. Rep., 163, 63 [NASA ADS] [CrossRef] [Google Scholar]
Wongwathanarat, A., Hammer, N. J., & Müller, E. 2010, A&A, 514, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Wongwathanarat, A., Müller, E., & Janka, H.-T. 2015, A&A, 577, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Yee, H. C., Vinokur, M., & Djomehri, M. J. 2000, J. Comput. Phys., 162, 33 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Coarse-fine interpolation

We consider a coarse grid and a fine grid in 1D, which is refined from the coarse grid by a factor of two, i.e., fine grid cells i and i + 1 are the subdivisions of the coarse grid cell i′. To compute a fourth-order accurate approximation of a scalar quantity $q_{i^{'}}^{c}$ $\hbox{$q^c_{i^\prime}$}$ at a cell center of a zone i′ on the coarse grid from cell averages ⟨ q^f ⟩ on the fine grid, we expand q with respect to the cell center of the coarse grid as $q (x) = a_{0} + a_{1} x + a_{2} x^{2} + a_{3} x^{3} .$ $\appendix \setcounter{section}{1} \begin{equation} q(x) = a_0+a_1x+a_2x^2+a_3x^3 \label{eq:polyexp1d} . \end{equation}$ (A.1)Integrating this equation to obtain cell averages ⟨ q^f ⟩ for the zones i−1, i, i + 1 and i + 2 yields a system of equations $\begin{matrix} ⟨ q^{f} ⟩_{i - 1} = a_{0} - \frac{3}{2} a_{1} h + \frac{7}{3} a_{2} h^{2} - \frac{15}{4} a_{3} h^{3}, \\ ⟨ q^{f} ⟩_{i} = a_{0} - \frac{1}{2} a_{1} h + \frac{1}{3} a_{2} h^{2} - \frac{1}{4} a_{3} h^{3}, \\ ⟨ q^{f} ⟩_{i + 1} = a_{0} + \frac{1}{2} a_{1} h + \frac{1}{3} a_{2} h^{2} + \frac{1}{4} a_{3} h^{3}, \\ ⟨ q^{f} ⟩_{i + 2} = a_{0} + \frac{3}{2} a_{1} h + \frac{7}{3} a_{2} h^{2} + \frac{15}{4} a_{3} h^{3}, \end{matrix}$ $\appendix \setcounter{section}{1} \begin{eqnarray} &&\avg{q^f}_{i-1}= a_0-\frac{3}{2}a_1h+\frac{7}{3}a_2h^2-\frac{15}{4}a_3h^3,\\ &&\avg{q^f}_{i\ph{-1}}= a_0-\frac{1}{2}a_1h+\frac{1}{3}a_2h^2-\frac{1}{4}a_3h^3,\\ &&\avg{q^f}_{i+1}= a_0+\frac{1}{2}a_1h+\frac{1}{3}a_2h^2+\frac{1}{4}a_3h^3,\\ &&\avg{q^f}_{i+2}= a_0+\frac{3}{2}a_1h+\frac{7}{3}a_2h^2+\frac{15}{4}a_3h^3, \end{eqnarray}$ where h is the cell spacing on the fine grid. Solving this system of equations for a₀ gives $q_{i^{'}}^{c} = a_{0} = \frac{7}{12} (⟨ q^{f} ⟩_{i} + ⟨ q^{f} ⟩_{i + 1}) - \frac{1}{12} (⟨ q^{f} ⟩_{i - 1} + ⟨ q^{f} ⟩_{i + 2}) .$ $\appendix \setcounter{section}{1} \begin{equation} q^c_{i^\prime} = a_0 = \frac{7}{12}(\avg{q^f}_{i}+\avg{q^f}_{i+1}) -\frac{1}{12}(\avg{q^f}_{i-1}+\avg{q^f}_{i+2}). \end{equation}$ (A.6)The same procedure can be applied in a straightforward manner in 2D and 3D. The respective expressions for the 2D case read $\begin{matrix} q_{i^{'}, j^{'}}^{c} & = & \frac{1}{3} (⟨ q^{f} ⟩_{i,j} + ⟨ q^{f} ⟩_{i + 1,j} + ⟨ q^{f} ⟩_{i,j + 1} + ⟨ q^{f} ⟩_{i + 1,j + 1}) \\ - \frac{1}{24} (⟨ q^{f} ⟩_{i - 1,j} + ⟨ q^{f} ⟩_{i - 1,j + 1} + ⟨ q^{f} ⟩_{i + 2,j} + ⟨ q^{f} ⟩_{i + 2,j + 1} \\ \times ⟨ q^{f} ⟩_{i,j - 1} + ⟨ q^{f} ⟩_{i,j + 2} + ⟨ q^{f} ⟩_{i + 1,j - 1} + ⟨ q^{f} ⟩_{i + 1,j + 2}), \end{matrix}$ $\appendix \setcounter{section}{1} \begin{eqnarray} q^c_{i^\prime,j^\prime}& =\ &\frac{1}{3} (\avg{q^f}_{i,j}+\avg{q^f}_{i+1,j}+\avg{q^f}_{i,j+1}+\avg{q^f}_{i+1,j+1})\nonumber\\ \quad &&-\frac{1}{24}(\avg{q^f}_{i-1,j}+\avg{q^f}_{i-1,j+1}+\avg{q^f}_{i+2,j}+ \avg{q^f}_{i+2,j+1}\nonumber\\ \quad&&\times \avg{q^f}_{i,j-1}+\avg{q^f}_{i,j+2}+\avg{q^f}_{i+1,j-1}+\avg{q^f}_{i+1,j+2}), \end{eqnarray}$ (A.7)

and $\begin{matrix} q_{i^{'}, j^{'}, k^{'}}^{c} & = & \frac{9}{48} (⟨ q^{f} ⟩_{i,j,k} + ⟨ q^{f} ⟩_{i + 1,j,k} + ⟨ q^{f} ⟩_{i,j + 1,k} + ⟨ q^{f} ⟩_{i,j,k + 1} \\ + ⟨ q^{f} ⟩_{i + 1,j + 1,k} + ⟨ q^{f} ⟩_{i + 1,j,k + 1} + ⟨ q^{f} ⟩_{i,j + 1,k + 1} + ⟨ q^{f} ⟩_{i + 1,j + 1,k + 1}) \\ - \frac{1}{48} (⟨ q^{f} ⟩_{i - 1,j,k} + ⟨ q^{f} ⟩_{i - 1,j + 1,k} + ⟨ q^{f} ⟩_{i - 1,j,k + 1} \\ + ⟨ q^{f} ⟩_{i - 1,j + 1,k + 1} + ⟨ q^{f} ⟩_{i + 2,j,k} + ⟨ q^{f} ⟩_{i + 2,j + 1,k} \\ + ⟨ q^{f} ⟩_{i + 2,j,k + 1} + ⟨ q^{f} ⟩_{i + 2,j + 1,k + 1} + ⟨ q^{f} ⟩_{i,j - 1,k} \\ + ⟨ q^{f} ⟩_{i,j - 1,k + 1} + ⟨ q^{f} ⟩_{i + 1,j - 1,k} + ⟨ q^{f} ⟩_{i + 1,j - 1,k + 1} \\ + ⟨ q^{f} ⟩_{i,j + 2,k} + ⟨ q^{f} ⟩_{i,j + 2,k + 1} + ⟨ q^{f} ⟩_{i + 1,j + 2,k} \\ + ⟨ q^{f} ⟩_{i + 1,j + 2,k + 1} + ⟨ q^{f} ⟩_{i,j,k - 1} + ⟨ q^{f} ⟩_{i,j + 1,k - 1} \\ + ⟨ q^{f} ⟩_{i + 1,j,k - 1} + ⟨ q^{f} ⟩_{i + 1,j + 1,k - 1} + ⟨ q^{f} ⟩_{i,j,k + 2} \\ + ⟨ q^{f} ⟩_{i,j + 1,k + 2} + ⟨ q^{f} ⟩_{i + 1,j,k + 2} + ⟨ q^{f} ⟩_{i + 1,j + 1,k + 2}) . \end{matrix}$ $\appendix \setcounter{section}{1} \begin{eqnarray} q^c_{i^\prime,j^\prime,k^\prime}& =\ &\frac{9}{48}( \avg{q^f}_{i,j,k}+\avg{q^f}_{i+1,j,k}+\avg{q^f}_{i,j+1,k}+\avg{q^f}_{i,j,k+1} \nonumber\\ &&+\avg{q^f}_{i+1,j+1,k}\!+\!\avg{q^f}_{i+1,j,k+1} \!+\!\avg{q^f}_{i,j+1,k+1}\!+\!\avg{q^f}_{i+1,j+1,k+1})\nonumber\\ &&-\frac{1}{48}( \avg{q^f}_{i-1,j,k}+\avg{q^f}_{i-1,j+1,k}+\avg{q^f}_{i-1,j,k+1} \nonumber\\ && +\avg{q^f}_{i-1,j+1,k+1}+\avg{q^f}_{i+2,j,k}+\avg{q^f}_{i+2,j+1,k} \nonumber\\ && +\avg{q^f}_{i+2,j,k+1}+\avg{q^f}_{i+2,j+1,k+1}+\avg{q^f}_{i,j-1,k} \nonumber\\ && +\avg{q^f}_{i,j-1,k+1}+\avg{q^f}_{i+1,j-1,k}+\avg{q^f}_{i+1,j-1,k+1} \nonumber\\ && +\avg{q^f}_{i,j+2,k}+\avg{q^f}_{i,j+2,k+1}+\avg{q^f}_{i+1,j+2,k} \nonumber\\ && +\avg{q^f}_{i+1,j+2,k+1}+\avg{q^f}_{i,j,k-1}+\avg{q^f}_{i,j+1,k-1} \nonumber\\ && +\avg{q^f}_{i+1,j,k-1}+\avg{q^f}_{i+1,j+1,k-1}+\avg{q^f}_{i,j,k+2} \nonumber\\ && +\avg{q^f}_{i,j+1,k+2}+\avg{q^f}_{i+1,j,k+2}+\avg{q^f}_{i+1,j+1,k+2}). \end{eqnarray}$ (A.8)for the 3D case.

All Tables

Table 1

L₁ and L_∞ errors of ⟨ ρ ⟩, and the convergence rates for the 1D linear advection test.

In the text

Table 2

Same as Table 1, but for the 2D linear advection test.

In the text

Table 3

Sum of the mean errors of all conserved variables ϵ, defined in Eq. (68), and its convergence rate for the linear acoustic wave test.

In the text

Table 4

L₁-norm of the coarse-fine differences of the pointwise densities ρ as a function of grid resolution for the linear acoustic wave test.

In the text

Table 5

L₁ and L_∞ errors of ρ, and their convergence rates for the 2D advection of a nonlinear vortex test.

In the text

Table 6

L₁ and L_∞ errors and their convergence rates for the 2D advection of a nonlinear vortex on grids with four different smoothing widths b (see Sect. 4.4).

In the text

Table 7

Dependence of the maximum expansion ratio of the grid c and the L₁ convergence rate on the grid resolution and smoothing width b.

In the text

Table 8

Parameters used in our simulations to generate a turbulent flow (for more details see text).

In the text

All Figures

	Fig. 1 Examples of 2D meshes generated using smooth mapping M₁ (left) and circular mapping M₂ (right).
In the text

	Fig. 2 Density profiles for the 1D linear advection test obtained on a grid of 128 zones at t = 10. Results computed with Apsara and Prometheus are shown with blue crosses and red circles, respectively. The exact solution is shown in black.
In the text

	Fig. 3 Mach number distributions normalized by the maximum Mach number M_max at t = 1 for the Gresho vortex test. The top and bottom rows show results computed by Apsara and Prometheus, respectively. Regardless of M_max the vortex remains intact when simulated with Apsara, while with Prometheus the vortex gradually dissolves because of numerical dissipation as M_max decreases.
In the text

	Fig. 4 Total kinetic energy as a function of time E_kin(t) normalized by the total kinetic of the initial state E_kin(0) for the Gresho vortex problem. The solid and dashed lines show the results computed with Apsara and Prometheus, respectively. Curves of different colors correspond to different maximum Mach numbers. In the case of Apsara, these curves are visually indistinguishable.
In the text

	Fig. 5 Mean turbulent kinetic energy, mean dissipation rate of turbulent kinetic energy, mean enstrophy, and numerical Reynolds number (from top to bottom) plotted as function of the dimensionless time t^∗ for the Taylor-Green vortex test for the reference Mach number M_ref ~ 0.29.
In the text

	Fig. 6 Same as Fig. 5, but for M_ref = 10^-2.
In the text

	Fig. 7 Square root of the magnitude of the vorticity $\sqrt{\| \nabla \times u \|}$ $\hbox{$\sqrt{\left\|\bs{\nabla}\times\bs{u}\right\|}$}$ in the x-z plane along y = 0 from the turbulence simulation with 512³ grid cells at time t = 100.
In the text

	Fig. 8 Energy spectra (Eq. (89)) plotted vs. the dimensionless wavenumber 512kΔx for four different grid resolutions. The spectra are compensated by a k^5/3 spectrum, so that any region with Kolmogorov scaling appears flat.
In the text

	Fig. 9 Effective numerical viscosity (Eq. (91)) as a function of wavenumber estimated from our ILES simulations for four different grid resolutions.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.

[1] Abdikamalov, E., Ott, C. D., Radice, D., et al. 2015, ApJ, 808, 70 [NASA ADS] [CrossRef] [Google Scholar]

[2] Appel, W. 2007, Mathematics for Physics and Physicists (Princeton University Press) [Google Scholar]

[3] Bassi, F., & Rebay, S. 1997, J. Comput. Phys., 138, 251 [Google Scholar]

[4] Borges, R., Carmona, M., Costa, B., & Don, W. S. 2008, J. Comput. Phys., 227, 3191 [Google Scholar]

[5] Buchmüller, P., & Helzel, C. 2014, J. Sci. Comput., 61, 343 [CrossRef] [Google Scholar]

[6] Calhoun, D. A., Helzel, C., & Leveque, R. J. 2008, SIAM Rev., 50, 723 [CrossRef] [Google Scholar]

[7] Cardall, C. Y., Budiardja, R. D., Endeve, E., & Mezzacappa, A. 2014, ApJS, 210, 17 [NASA ADS] [CrossRef] [Google Scholar]

[8] Colella, P., & Glaz, H. M. 1985, J. Comput. Phys., 59, 264 [NASA ADS] [CrossRef] [Google Scholar]

[9] Colella, P., & Sekora, M. D. 2008, J. Comput. Phys., 227, 7069 [NASA ADS] [CrossRef] [Google Scholar]

[10] Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]

[11] Colella, P., Dorr, M. R., Hittinger, J. A. F., & Martin, D. F. 2011, J. Comput. Phys., 230, 2952 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[12] Couch, S. M., & Ott, C. D. 2015, ApJ, 799, 5 [NASA ADS] [CrossRef] [Google Scholar]

[13] Drikakis, D., Fureby, C., Grinstein, F. F., & Youngs, D. 2007, J. Turbulence, 8, N20 [CrossRef] [Google Scholar]

[14] Dubey, A., Antypas, K., Ganapathy, M. K., et al. 2009, Parallel Computing, 35, 512 [CrossRef] [Google Scholar]

[15] Edelmann, P. 2014, Dissertation, Technische Universität München, München [Google Scholar]

[16] Eswaran, V., & Pope, S. B. 1988, Computers and Fluids, 16, 257 [NASA ADS] [CrossRef] [Google Scholar]

[17] Feng, X., Jiang, C., Xiang, C., Zhao, X., & Wu, S. T. 2012, ApJ, 758, 62 [Google Scholar]

[18] Fragile, P. C., Lindner, C. C., Anninos, P., & Salmonson, J. D. 2009, ApJ, 691, 482 [NASA ADS] [CrossRef] [Google Scholar]

[19] Fryxell, B., Arnett, D., & Mueller, E. 1991, ApJ, 367, 619 [NASA ADS] [CrossRef] [Google Scholar]

[20] Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273 [NASA ADS] [CrossRef] [Google Scholar]

[21] Gresho, P. M., & Chan, S. T. 1990, Int. J. Numer. Meth. Fluids, 11, 621 [Google Scholar]

[22] Grimm-Strele, H., Kupka, F., & Muthsam, H. J. 2014, Comput. Phys. Comm., 185, 764 [Google Scholar]

[23] Guzik, S., McCorquodale, P., & Colella, P. 2012, in 50th AIAA Aerospace Sciences Meetings (American Institute of Aeronautics and Astronautics, Nashville, TN), AIAA Paper 2012–0574 [Google Scholar]

[24] Hayashi, H., & Kageyama, A. 2016, J. Comput. Phys., 305, 895 [NASA ADS] [CrossRef] [Google Scholar]

[25] Hotta, H., Rempel, M., & Yokoyama, T. 2014, ApJ, 786, 24 [NASA ADS] [CrossRef] [Google Scholar]

[26] Jiang, C., Feng, X., & Xiang, C. 2012, ApJ, 755, 62 [NASA ADS] [CrossRef] [Google Scholar]

[27] Kageyama, A., & Sato, T. 2004, Geophys., Geosys., 5, 9005 [Google Scholar]

[28] Ketcheson, D. I. 2008, SIAM J. Sci. Comput., 30, 2113 [CrossRef] [Google Scholar]

[29] Kifonidis, K., & Müller, E. 2012, A&A, 544, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[30] Koldoba, A. V., Romanova, M. M., Ustyugova, G. V., & Lovelace, R. V. E. 2002, ApJ, 576, L53 [NASA ADS] [CrossRef] [Google Scholar]

[31] Lemoine, G. I., & Ou, M. Y. 2014, SIAM J. Sci. Comput., 36, B396 [CrossRef] [Google Scholar]

[32] Lentz, E. J., Bruenn, S. W., Hix, W. R., et al. 2015, ApJ, 807, L31 [NASA ADS] [CrossRef] [Google Scholar]

[33] LeVeque, R. 2002, Finite Volume Methods for Hyperbolic Problems, Cambridge Texts in Applied Mathematics (Cambridge University Press) [Google Scholar]

[34] McCorquodale, P., & Colella, P. 2011, CAMCOS, 6, 1 [CrossRef] [Google Scholar]

[35] McCorquodale, P., Dorr, M. R., Hittinger, J. A. F., & Colella, P. 2015, J. Comput. Phys., 288, 181 [NASA ADS] [CrossRef] [Google Scholar]

[36] Melson, T., Janka, H.-T., Bollig, R., et al. 2015a, ApJ, 808, L42 [NASA ADS] [CrossRef] [Google Scholar]

[37] Melson, T., Janka, H.-T., & Marek, A. 2015b, ApJ, 801, L24 [NASA ADS] [CrossRef] [Google Scholar]

[38] Miczek, F., Röpke, F. K., & Edelmann, P. V. F. 2015, A&A, 576, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[39] Mueller, E., Fryxell, B., & Arnett, D. 1991, A&A, 251, 505 [NASA ADS] [Google Scholar]

[40] Müller, B. 2015, MNRAS, 453, 287 [NASA ADS] [CrossRef] [Google Scholar]

[41] Murphy, J. W., Dolence, J. C., & Burrows, A. 2013, ApJ, 771, 52 [NASA ADS] [CrossRef] [Google Scholar]

[42] Nonomura, T., Iizuka, N., & Fujii, K. 2010, Computers & Fluids, 39, 197 [CrossRef] [Google Scholar]

[43] Peng, X., Xiao, F., & Takahashi, K. 2006, Quart. J. Roy. Meteor. Soc., 132, 979 [NASA ADS] [CrossRef] [Google Scholar]

[44] Radice, D., Couch, S. M., & Ott, C. D. 2015, Comput. Astrophys. Cosmol., 2, 7 [NASA ADS] [CrossRef] [Google Scholar]

[45] Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., & Lovelace, R. V. E. 2012, MNRAS, 421, 63 [NASA ADS] [Google Scholar]

[46] Ronchi, C., Iacono, R., & Paolucci, P. S. 1996, J. Comput. Phys., 124, 93 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]

[47] Shiota, D., Kusano, K., Miyoshi, T., & Shibata, K. 2010, ApJ, 718, 1305 [NASA ADS] [CrossRef] [Google Scholar]

[48] Shu, C.-W. 2003, Int. J. Comput. Fluid Dyn., 17, 107 [NASA ADS] [CrossRef] [Google Scholar]

[49] Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137 [NASA ADS] [CrossRef] [Google Scholar]

[50] Sytine, I. V., Porter, D. H., Woodward, P. R., Hodson, S. W., & Winkler, K.-H. 2000, J. Comput. Phys., 158, 225 [NASA ADS] [CrossRef] [Google Scholar]

[51] Taylor, G. I., & Green, A. E. 1937, Proc. Roy. Soc. London A, 158, 499 [Google Scholar]

[52] Toro, E. F., Spruce, M., & Speares, W. 1994, Shock Waves, 4, 25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[53] Uhlenbeck, G. E., & Ornstein, L. S. 1930, Phys. Rev., 36, 823 [NASA ADS] [CrossRef] [Google Scholar]

[54] Visbal, M. R., & Gaitonde, D. V. 2002, J. Comput. Phys., 181, 155 [NASA ADS] [CrossRef] [Google Scholar]

[55] Warren, D. C., & Blondin, J. M. 2013, MNRAS, 429, 3099 [NASA ADS] [CrossRef] [Google Scholar]

[56] Wilson, J. R., & Mayle, R. W. 1988, Phys. Rep., 163, 63 [NASA ADS] [CrossRef] [Google Scholar]

[57] Wongwathanarat, A., Hammer, N. J., & Müller, E. 2010, A&A, 514, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[58] Wongwathanarat, A., Müller, E., & Janka, H.-T. 2015, A&A, 577, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[59] Yee, H. C., Vinokur, M., & Djomehri, M. J. 2000, J. Comput. Phys., 162, 33 [NASA ADS] [CrossRef] [Google Scholar]

APSARA: A multi-dimensional unsplit fourth-order explicit Eulerian hydrodynamics code for arbitrary curvilinear grids

1. Introduction

2. Numerical methods

2.1. Governing equations

2.2. Temporal discretization

2.3. Cell and face averages

2.4. Evaluation of face average fluxes

3. Mapping functions

3.1. Cartesian mapping M0

3.2. Smooth mapping M1

3.3. Circular mapping M2

4. Hydrodynamic tests

4.1. Linear advection

4.2. Linear acoustic wave

4.3. Advection of a nonlinear vortex

4.4. Influence of grid smoothness

4.5. Gresho vortex

4.6. Taylor-Green vortex

5. Turbulence in the context of core collapse supernovae

6. Summary

Acknowledgments

References

Appendix A: Coarse-fine interpolation

All Tables

All Figures

3.1. Cartesian mapping M₀

3.2. Smooth mapping M₁

3.3. Circular mapping M₂