WENO–WOMBAT: Scalable Fifth-order Constrained-transport Magnetohydrodynamics for Astrophysical Applications

Following Jiang & Shu (1996), we use a fourth-order, four-stage Runge–Kutta (RK4) time integrator. It has been shown that such a scheme cannot have the strong-stability-preserving property (e.g., Spiteri & Ruuth 2002). In fact, the scheme used here is not even total variation diminishing (TVD; Shu 1988; Shu & Osher 1988). Nonetheless, it allows us to use $\mathrm{CFL}=0.8$ everywhere, even in three dimensions, and is sufficiently robust in practice. Modern SSP RK schemes are better behaved mathematically, but have to involve five or more stages at fourth order and thus are computationally 25% more expensive than the fourth-order scheme at the same time-step size (e.g., Gottlieb et al. 2001; Spiteri & Ruuth 2002; Gottlieb et al. 2009). As we aim for computational efficiency per solution quality, we leave tests of SSP integrators to future work.

In N_d dimensions, the time step is set as

$$\begin{eqnarray}&&{\rm{\Delta }}t\lt \displaystyle \frac{\mathrm{CFL}{\rm{\Delta }}x}{{\displaystyle \sum }_{d=1}^{{N}_{d}}\max (| {v}_{d}| +{c}_{\mathrm{fast}}^{d}))}.\end{eqnarray} \tag{ 17 }$$

Omitting CT for now, the RK4 scheme from Jiang & Shu (1996) and Jiang & Wu (1999) can be written as

$$\begin{eqnarray*}\begin{array}{rcl}{\boldsymbol{q}} & = & {{\boldsymbol{q}}}_{{\boldsymbol{0}}}={{\boldsymbol{q}}}_{i,j,k}^{n}\\ {{\boldsymbol{q}}}_{\mathrm{save}} & = & -\displaystyle \frac{4}{3}{{\boldsymbol{q}}}_{{\boldsymbol{0}}}\\ {\boldsymbol{q}} & = & {{\boldsymbol{q}}}_{{\boldsymbol{0}}}-\displaystyle \frac{1}{2}\displaystyle \frac{{\rm{\Delta }}t}{{\rm{\Delta }}x}\left({{\boldsymbol{F}}}_{i+1/2,j,k}({\boldsymbol{q}})-{{\boldsymbol{F}}}_{i-1/2,j,k}({\boldsymbol{q}})\right)\\ {{\boldsymbol{q}}}_{\mathrm{save}} & = & {{\boldsymbol{q}}}_{\mathrm{save}}+\displaystyle \frac{1}{3}{\boldsymbol{q}}\\ {\boldsymbol{q}} & = & {{\boldsymbol{q}}}_{0}-\displaystyle \frac{1}{2}\displaystyle \frac{{\rm{\Delta }}t}{{\rm{\Delta }}x}\left({{\boldsymbol{F}}}_{i+1/2,j,k}({\boldsymbol{q}})-{{\boldsymbol{F}}}_{i-1/2,j,k}({\boldsymbol{q}})\right)\\ {{\boldsymbol{q}}}_{\mathrm{save}} & = & {{\boldsymbol{q}}}_{\mathrm{save}}+\displaystyle \frac{2}{3}{\boldsymbol{q}}\\ {\boldsymbol{q}} & = & {{\boldsymbol{q}}}_{0}-\displaystyle \frac{{\rm{\Delta }}t}{{\rm{\Delta }}x}\left({{\boldsymbol{F}}}_{i+1/2,j,k}({\boldsymbol{q}})-{{\boldsymbol{F}}}_{i-1/2,j,k}({\boldsymbol{q}})\right)\\ {{\boldsymbol{q}}}_{\mathrm{save}} & = & {{\boldsymbol{q}}}_{\mathrm{save}}+\displaystyle \frac{1}{3}{\boldsymbol{q}}\\ {\boldsymbol{q}} & = & {{\boldsymbol{q}}}_{0}-\displaystyle \frac{1}{6}\displaystyle \frac{{\rm{\Delta }}t}{{\rm{\Delta }}x}\left({{\boldsymbol{F}}}_{i+1/2,j,k}({\boldsymbol{q}})-{{\boldsymbol{F}}}_{i-1/2,j,k}({\boldsymbol{q}})\right)\\ {{\boldsymbol{q}}}_{\mathrm{save}} & = & {{\boldsymbol{q}}}_{\mathrm{save}}+{\boldsymbol{q}}\\ {{\boldsymbol{q}}}_{i,j,k}^{n+1} & = & {{\boldsymbol{q}}}_{\mathrm{save}}.\end{array}\end{eqnarray*}$$

The CT time integration scheme is the same, but updates the face-centered magnetic fields after the WENO step (see Section 2.3). In three dimensions, it adds the three-component magnetic field on the boundary and the three-component corner flux to global storage.

For reference, we outline the WENO5 spatial discretization to obtain fluxes at the zone boundaries in Equation (16). The classical scheme from Jiang & Wu (1999) is an FD approach that efficiently computes boundary fluxes from point-valued cell-centered states and fluxes. This is in contrast to modern schemes that often use a more flexible, but also more computationally expensive, FV approach. For more details, we kindly ask the reader to refer to the extensive literature on the subject (e.g., Jiang & Shu 1996; Shu 1998; Jiang & Wu 1999; Shu 2009; H. Jang et al. 2019, in preparation).

The fifth-order WENO scheme uses a weighted average of three third-order stencils to approximate the solution at the boundary. Thus, the WENO averaging itself requires two boundary zones. Including an additional zone from the flux difference in Equation (16) gives a total of three boundary zones per RK substep. To avoid Gibbs phenomena (ringing) near shocks, the weights "switch off" one of the three polynomials adjacent to discontinuities and inside a discontinuity, where the method becomes first order.

We note that in two or three dimensions, all fluxes are applied at the same time, starting from the same initial state vector. Strang splitting is not used (Shu 2003). For the sake of efficiency, the scheme does not include dimensional cross-terms found in many recent FV schemes (see, e.g., Balsara et al. 2009; Verma & Müller 2018). This strategy is common in the FD WENO literature that commonly aims at cheap high-order implementations (e.g., Shu 2003). The general argument is that the time-integration scheme implicitly averages over all three directions four times and thus maintains the high order as long as cross-dimensional fluxes are not too large. The rotated three-dimensional linear wave advection test (Section 4.1) confirms that the WENO5 scheme is effectively unsplit at fifth order in the linear regime. Cross-terms might be more relevant for strongly nonlinear problems, where the influence of the grid can be seen, e.g., in the Sedov–Taylor blast wave test, Section 4.16. As we will show, the algorithmic fidelity of the FD WENO5 scheme is then comparable to a third-order discontinuous Galerkin scheme, however still at a fraction of the computational cost.

Following Jiang & Wu (1999), we denote quantities in the decoupled system with superscript s, denote the individual components with $m\in [1,7]$, and drop the vector notation. Thus, the fluxes, state vectors, and their differences in the decoupled system are

$$\begin{eqnarray}&&{F}_{k}^{s}(m)={{\boldsymbol{L}}}_{i+\tfrac{1}{2}}(m)\cdot {{\boldsymbol{F}}}_{k},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{q}_{k}^{s}(m)={{\boldsymbol{L}}}_{i+\tfrac{1}{2}}(m)\cdot {{\boldsymbol{q}}}_{k},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{F}_{k+\tfrac{1}{2}}^{s}(m)={F}_{k+1}^{s}(m)-{F}_{k}^{s}(m),\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{q}_{k+\tfrac{1}{2}}^{s}(m)={q}_{k+1}^{s}(m)-{q}_{k}^{s}(m),\end{eqnarray} \tag{ 21 }$$

respectively. Note that the left-hand eigenvector is taken at the zone boundary using simple arithmetic averaging. In the decoupled system, WENO uses Lax–Friedrichs-type flux splitting to upwind the fluxes (Lax 1954). Jiang & Wu (1999) and Shu (2003) have argued that this usually very diffusive approach is sufficiently accurate and computationally cheap in the FD WENO context:

$$\begin{eqnarray}&&{F}_{i}^{s\pm }(m)=\displaystyle \frac{1}{2}\left({F}_{i}^{s}(m)\pm \alpha (m){q}_{i}^{s}(m)\right)\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{F}_{k+\tfrac{1}{2}}^{s\pm }(m)=\displaystyle \frac{1}{2}\left({\rm{\Delta }}{F}_{k+\tfrac{1}{2}}^{s}(m)\pm \alpha (m){\rm{\Delta }}{q}_{k+\tfrac{1}{2}}^{s}(m)\right),\end{eqnarray} \tag{ 23 }$$

where $\alpha (m)=\max (| {\lambda }_{i}(m)| ,| {\lambda }_{i+1}(m)| )$ is the larger of the two eigenvalues $\lambda (m)$ between zone i and $i+1$.

The WENO-interpolated fluxes at the zone boundary for Equation (16) are given by ${{\boldsymbol{F}}}_{i+\tfrac{1}{2}}={F}_{i+\tfrac{1}{2}}^{s}(m)\cdot {{\boldsymbol{R}}}_{i+\tfrac{1}{2}}(m)$, where again the right-hand eigenvector is taken at the zone boundary. In the decoupled system, the flux of each component m is (Jiang & Wu 1999)

$$\begin{eqnarray}\begin{array}{rcl}{F}_{i+\tfrac{1}{2}}^{s} & = & \displaystyle \frac{1}{12}\left(-{F}_{i-1}^{s}+7{F}_{i}^{s}+7{F}_{i+1}^{s}-{F}_{i+2}^{s}\right)\\ & & -\,{\varphi }_{N}\left({\rm{\Delta }}{F}_{i-\tfrac{3}{2}}^{s+},{\rm{\Delta }}{F}_{i-\tfrac{1}{2}}^{s+},{\rm{\Delta }}{F}_{i+\tfrac{1}{2}}^{s+},{\rm{\Delta }}{F}_{i+\tfrac{3}{2}}^{s+}\right)\\ & & +\,{\varphi }_{N}\left({\rm{\Delta }}{F}_{i+\tfrac{5}{2}}^{s-},{\rm{\Delta }}{F}_{i+\tfrac{3}{2}}^{s-},{\rm{\Delta }}{F}_{i+\tfrac{1}{2}}^{s-},{\rm{\Delta }}{F}_{i-\tfrac{1}{2}}^{s-}\right),\end{array}\end{eqnarray} \tag{ 24 }$$

where we have dropped the (m) notation for clarity. The WENO interpolant ${\varphi }_{N}(a,b,c,d)$ is defined as

$$\begin{eqnarray}&&{\varphi }_{N}=\displaystyle \frac{1}{3}{\omega }_{0}\left(a-2b+c\right)+\displaystyle \frac{1}{6}\left({\omega }_{2}-\displaystyle \frac{1}{2}\right)\left(b-2\,c+d\right).\end{eqnarray} \tag{ 25 }$$

The nonlinear weights are

$$\begin{eqnarray}&&{\omega }_{0}=\displaystyle \frac{{\alpha }_{0}}{{\alpha }_{0}+{\alpha }_{1}+{\alpha }_{2}},{\omega }_{2}=\displaystyle \frac{{\alpha }_{2}}{{\alpha }_{0}+{\alpha }_{1}+{\alpha }_{2}}\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{\alpha }_{0}=\displaystyle \frac{1}{{\left(\epsilon +{{IC}}_{0}\right)}^{2}},{\alpha }_{1}=\displaystyle \frac{6}{{\left(\epsilon +{{IC}}_{1}\right)}^{2}},\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{\alpha }_{2}=\displaystyle \frac{3}{{\left(\epsilon +{{IC}}_{2}\right)}^{2}},\end{eqnarray} \tag{ 28 }$$

with  = 10⁻⁶ and

$$\begin{eqnarray}&&{{IS}}_{0}=13{\left(a-b\right)}^{2}+3{\left(a-3b\right)}^{2},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{{IS}}_{1}=13{\left(b-c\right)}^{2}+3{\left(b+c\right)}^{2},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{{IS}}_{2}=13{\left(c-d\right)}^{2}+3{\left(3c-d\right)}^{2}.\end{eqnarray} \tag{ 31 }$$

This description of the weights is generally fifth-order accurate in smooth flows, but only third-order accurate near critical points (extrema and saddle points). Borges et al. (2008) proposed the WENO-Z weights to make the scheme truly fifth order everywhere:

$$\begin{eqnarray}&&{\tau }_{5}=| {{IS}}_{0}-{{IS}}_{2}| \end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{\alpha }_{0}=1+\displaystyle \frac{{\tau }_{5}}{{\left({\epsilon }_{Z}+{{IS}}_{0}\right)}^{2}},\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{\alpha }_{1}=6\left(1+\displaystyle \frac{{\tau }_{5}}{{\left({\epsilon }_{Z}+{{IS}}_{1}\right)}^{2}}\right),\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{\alpha }_{2}=3\left(1+\displaystyle \frac{{\tau }_{5}}{{\left({\epsilon }_{Z}+{{IS}}_{2}\right)}^{2}}\right),\end{eqnarray} \tag{ 35 }$$

and ${\epsilon }_{Z}={10}^{-40}$. As we will see, this is useful in advection-dominated problems, but reduces the general robustness of the scheme. For complex problems, WENO-Z will fall back into protection fluxes more often, which might not aid the solution for all problems.

Equation (5) requires maintenance of a zero magnetic field divergence at all times. In CT, the magnetic field is defined on a staggered mesh, i.e., on the faces of the computational zones (Evans & Hawley 1988; Dai & Woodward 1998; Balsara & Spicer 1999b; Tóth 2000). This way, magnetic field divergence can be maintained to machine precision. For computational efficiency, we follow (Ryu et al. 1998) in our implementation of CT, which uses second-order accurate averages to obtain the corner fluxes. As we will see, this choice introduces additional dispersion to the solution, but not additional diffusion.

We denote the magnetic field on the forward zone faces (i.e., $i+\tfrac{1}{2}$ for the x-component) as ${\boldsymbol{b}}={({b}_{x},{b}_{y},{b}_{z})}^{T}$. From the Riemann solver, we obtain the fifth-order accurate magnetic field fluxes across the zone faces f_y, f_z during the sweep in the x-direction (components six and seven); g_z, g_x for the y-direction; and h_x, h_y for the z-direction:

$$\begin{eqnarray}&&{f}_{y}^{* }\,=\,{f}_{y}+\displaystyle \frac{1}{2}\left({B}_{x,i,j,k}{v}_{y,i,j,k}+{B}_{x,i+1,j,k}{v}_{y,i+1,j,k}\right)\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&{f}_{z}^{* }\,=\,{f}_{z}+\displaystyle \frac{1}{2}\left({B}_{x,i,j,k}{v}_{z,i,j,k}+{B}_{x,i+1,j,k}{v}_{z,i+1,j,k}\right)\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{g}_{x}^{* }\,=\,{g}_{x}+\displaystyle \frac{1}{2}\left({B}_{y,i,j,k}{v}_{z,i,j,k}+{B}_{y,i,j+1,k}{v}_{z,i,j+1,k}\right)\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&{g}_{z}^{* }\,=\,{g}_{z}+\displaystyle \frac{1}{2}\left({B}_{y,i,j,k}{v}_{x,i,j,k}+{B}_{y,i,j+1,k}{v}_{x,i,j+1,k}\right)\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&{h}_{x}^{* }={h}_{x}+\displaystyle \frac{1}{2}\left({B}_{z,i,j,k}{v}_{x,i,j,k}+{B}_{z,i,j,k+1}{v}_{x,i,j,k+1}\right)\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&{h}_{y}^{* }={h}_{y}+\displaystyle \frac{1}{2}\left({B}_{z,i,j,k}{v}_{y,i,j,k}+{B}_{z,i,j,k+1}{v}_{u,i,j,k+1}\right).\end{eqnarray} \tag{ 41 }$$

Due to our implementation involving the rotation of the grid, the components actually remain the same for every direction. However, the fluxes have to be rotated backward (${g}_{x},{g}_{z}$) or forward (${h}_{x},{h}_{y}$) into the original frame. The corner fluxes ${{\rm{\Omega }}}_{x},{{\rm{\Omega }}}_{y},{{\rm{\Omega }}}_{z}$ are then

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{x}=\displaystyle \frac{1}{2}\left({g}_{x,i,j,k}^{* }+{g}_{x,i+1,j,k}^{* }\right)-\displaystyle \frac{1}{2}\left({f}_{y,i,j,k}^{* }+{f}_{y,i,j+1,k}^{* }\right)\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{y}=\displaystyle \frac{1}{2}\left({h}_{y,i,j,k}^{* }+{h}_{y,i,j+1,k}^{* }\right)-\displaystyle \frac{1}{2}\left({g}_{z,i,j,k}^{* }+{g}_{z,i,j,k+1}^{* }\right)\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{z}=\displaystyle \frac{1}{2}\left({f}_{z,i,j,k}^{* }+{f}_{z,i,j,k+1}^{* }\right)-\displaystyle \frac{1}{2}\left({h}_{x,i,j,k}^{* }+{h}_{x,i+1,j,k}^{* }\right).\end{eqnarray} \tag{ 44 }$$

The face-centered magnetic fields are then updated analogously to Equation (16):

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{{db}}_{x}}{{dt}}+\displaystyle \frac{1}{{\rm{\Delta }}x}\left({{\rm{\Omega }}}_{x,i,j,k}-{{\rm{\Omega }}}_{x,i,j-1,k}\right)\\ \qquad -\,\displaystyle \frac{1}{{\rm{\Delta }}x}\left({{\rm{\Omega }}}_{z,i,j,k}-{{\rm{\Omega }}}_{z,i,j,k-1}\right)=0\end{array}\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{{db}}_{y}}{{dt}}+\displaystyle \frac{1}{{\rm{\Delta }}x}\left({{\rm{\Omega }}}_{y,i,j,k}-{{\rm{\Omega }}}_{y,i,j,k-1}\right)\\ \qquad -\,\displaystyle \frac{1}{{\rm{\Delta }}x}\left({{\rm{\Omega }}}_{x,i,j,k}-{{\rm{\Omega }}}_{x,i-1,j,k}\right)=0\end{array}\end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{{db}}_{z}}{{dt}}+\displaystyle \frac{1}{{\rm{\Delta }}x}\left({{\rm{\Omega }}}_{z,i,j,k}-{{\rm{\Omega }}}_{z,i-1,j,k}\right)\\ \qquad -\,\displaystyle \frac{1}{{\rm{\Delta }}x}\left({{\rm{\Omega }}}_{y,i,j,k}-{{\rm{\Omega }}}_{y,i,j-1,k}\right)=0.\end{array}\end{eqnarray} \tag{ 47 }$$

The zone-centered magnetic field is interpolated at fourth order from face values:

$$\begin{eqnarray}\begin{array}{rcl}{B}_{x,i,j,k} & = & \displaystyle \frac{1}{16}\left(-{b}_{x,i-2,j,k}+9{b}_{x,i-1,j,k}\right.\\ & & \left.+\,9{b}_{x,i,j,k}-{b}_{x,i+1,j,k}\right),\end{array}\end{eqnarray} \tag{ 48 }$$

and the other components follow accordingly along the j and k indices.

We use the same fourth-order Runge–Kutta scheme presented in Section 2.1. The CT update is interleaved with the WENO update, with boundary communication of Ω between a WENO and a CT step. As we will see, the CT scheme converges to second order. However, it can be shown that magnetic field dissipation converges to fifth order, so the additional error is in the shape of the field, not its energy (see also H. Jang et al. 2019, in preparation).

In WOMBAT, every solver works on a single boundary communicated patch, where the grid data reside in rank global memory. At this point, no communication is required any longer; work and communication are completely separate in the implementation. This eases maintenance of the code considerably.

In the WENO solver, we first copy the grid data into OpenMP thread private buffers to minimize Non-Uniform Memory Access effects. This introduces memory overhead of about 25 MB per thread for patches with 18³ zones. In multiple dimensions, the grid is explicitly flattened as well, i.e., the number of indices of all arrays is reduced to one along spatial dimensions and all loops run over the flattened array, ignoring the boundaries between dimensions. Thus, memory and code-level layout of the data are identical. This increases vector length and decreases cache-blocked patch size, which in turn eases MPI load balancing. It also improves code readability, as all lower level routines (fluxes, eigenvectors) remain inherently one-dimensional. This layout naturally leads to an explicit separation of computation and data movement in the code. In two and three dimensions, the sweeps in the y-direction and z-direction are implemented as rotations of the plane/cube in flattened memory. The resulting WENO fluxes are later rotated back into the original coordinate system to compute the state-vector updates (Equation (16) and Section 2.1). Using the Fortran CONTIGUOUS keyword and implied array sizes, all relevant loops in the implementation auto-vectorize with current Cray and Intel compilers. This is a necessary prerequisite to achieve a significant fraction of peak performance on current CPUs. We anticipate that this implementation will be very convenient to port to accelerators (GPUs) as it naturally exposes long vectors and makes memory movement explicit.

Due to WOMBAT's parallelization strategy, the RK4 scheme presented in Section 2.1 is implemented using a running state vector q, a state vector buffer q_save to accumulate intermediate results, and a copy of the initial state ${{\boldsymbol{q}}}_{0}={{\boldsymbol{q}}}_{i}^{n}$. Note that we communicate boundaries twice (once for WENO, once for CT update) per RK4 substep. This allows us to retain three boundary zones on the patch, which is optimal regarding flops and memory consumption given that communication in WOMBAT is local, asynchronous, and thus extremely cheap.

3.1.1. WENO–WOMBAT in One Dimension

3.1.2. WENO–WOMBAT in Two Dimensions

3.1.3. WENO–WOMBAT in Three Dimensions

The convergence order of the scheme can be exposed by advecting a perturbation corresponding to one of the (M)HD eigenvectors through a periodic box with 2N × N × N zones and a domain of 3 × 1.5 × 1.5 in three dimensions. This test evaluates the code in the linear regime and is very valuable for debugging. As the setup is not trivial, especially in three dimensions, we describe it in greater detail here, but note that the WOMBAT or ATHENA source code is likely indispensable to fully reproduce the results. A very instructive discussion of advection comparing WENO5 with TVD in one dimension can be found in the introduction of Shu (2003).

Following Gardiner & Stone (2005, 2008) and Stone et al. (2008), we set γ = 5/3 and $\bar{{\boldsymbol{U}}}={(1,{v}_{x},0,0,{B}_{x},{B}_{y},{B}_{z},1/\gamma )}^{T}$, where v_x = 1 for shear and entropy modes and v_x = 0 otherwise. For hydrodynamic waves, B = 0; for MHD waves, ${\boldsymbol{B}}\,={(1,\sqrt{2},1/2)}^{T}$. We add a perturbation on the conserved variables: ${\boldsymbol{Q}}=\bar{{\boldsymbol{Q}}}+{A}_{0}{{\boldsymbol{R}}}_{k}\sin (2\pi x)$. Here, R_k is the right-hand eigenvector of mode k and ${A}_{0}={10}^{-6}$. We note that the perturbation is applied only once per component, i.e., the momentum fluctuation uses the unperturbed background density. The right eigenvectors for the (M)HD waves are

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{\mathrm{Alfven}}=\displaystyle \frac{1}{6\sqrt{5}}{\left(0,0,1,-2\sqrt{2},0,-\mathrm{1,\; 2}\sqrt{2},0\right)}^{T}\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{\mathrm{fast}}=\displaystyle \frac{1}{6\sqrt{5}}{\left(\mathrm{6,12},-4\sqrt{2},-\mathrm{2,\; 0,\; 8}\sqrt{2},\mathrm{4,27}\right)}^{T}\end{eqnarray} \tag{ 51 }$$

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{\mathrm{slow}}=\displaystyle \frac{1}{6\sqrt{5}}{\left(\mathrm{12,\; 6,\; 8}\sqrt{2},\mathrm{4,\; 0},-4\sqrt{2},-\mathrm{2,9}\right)}^{T}\end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{\mathrm{entropy}}={\left(\mathrm{1,1},\mathrm{0,\; 0},\mathrm{0,\; 0,\; 0},\displaystyle \frac{1}{2}\right)}^{T}.\end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{\mathrm{sound}}={\left(\mathrm{1,1},\mathrm{0,\; 0},\mathrm{0,\; 0},\mathrm{0,1}/(\gamma -1)\right)}^{T}\end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{\mathrm{shear},{\rm{y}}}={\left(\mathrm{0,\; 0},\mathrm{1,\; 0},\mathrm{0,\; 0},\mathrm{0,0}\right)}^{T}.\end{eqnarray} \tag{ 55 }$$

In three dimensions, we rotate the zone positions parallel to the wave vector following Gardiner & Stone (2008), with the rotation matrix

$$\begin{eqnarray}M(\alpha ,\beta )=\left[\begin{array}{ccc}\cos \alpha \cos \beta & -\sin \beta & \sin \alpha \cos \beta \\ \cos \alpha \sin \beta & \cos \beta & -\sin \alpha \sin \beta \\ \sin \alpha & 0 & \cos \alpha \end{array}\right],\end{eqnarray} \tag{ 56 }$$

where $\sin \alpha =2/3$ and $\cos \beta =2/\sqrt{5}$. We use its inverse/transpose to rotate back into the laboratory frame. We set a magnetic vector potential:

$$\begin{eqnarray}&&{A}_{x}=0\end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray}&&{A}_{y}=\displaystyle \frac{{A}_{0}}{2\pi }{R}_{k,7}\sin (2\pi x)+{B}_{z}{x}_{k}\end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray}&&{A}_{z}=-\displaystyle \frac{{A}_{0}}{2\pi }{R}_{k,6}\sin (2\pi x)-{B}_{y}{x}_{k}+{B}_{x}{y}_{k},\end{eqnarray} \tag{ 59 }$$

where x_k and y_k are the x and y components of the zone position rotated with the inverse rotation matrix. We note that the vector potential has to be evaluated on the edges of the zone in three dimensions, so that the resulting magnetic field is field-centered. An example in our case Ax is defined at $\left(i,j+\tfrac{1}{2},k+\tfrac{1}{2}\right)$. The vector potential is then multiplied with the forward rotation matrix, and the magnetic field on the interfaces is obtained using a standard first-order FD rotation operator. We observe ${\rm{\nabla }}\cdot {\boldsymbol{B}}\lt {10}^{-12}$ in all tests at all times.

In Figure 1, we compare the profile of the density contrast at t = 10 in one dimension with the analytic solution. Simulations are done with 8, 16, 32, and 64 zones and shown as violet diamonds, green square, blue crosses, and red circles, respectively. Dispersion is visible only at the lowest resolution, eight zones. This compares well with results shown in Wongwathanarat et al. (2016). In Figure 2, we plot the resulting L₁ error for resolutions between 8 and 128 zones run with WENO5-Z. From top left to bottom right, we show 1D hydrodynamic waves, 3D hydrodynamic waves, 1D MHD waves, and 3D MHD waves. In the 1D case, the simulation converges as ${N}^{-5}$, as expected, down to 10⁻¹¹, where the effect of wave steepening begins to limit further convergence of the compressive modes. We note that for $N\gt 128$, the round-off error from the 8 byte variables leads to an increase in L₁ again (see also Donnert et al. 2018). This suggests that the scheme needs to be run with at least 8 byte precision to take advantage of the high-order convergence properties.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In three dimensions, all hydrodynamic waves converge to fifth order. The entropy MHD wave converges to fifth order in L₁ as well. This mode does not feature a fluctuation in the magnetic field (Equation (53)); the background magnetic field is only advected. All other MHD waves converge to second order in L₁, because our CT scheme is geometrically only second order. However, the fluxes used are of course fifth-order accurate, thus this error should be dominated by dispersion, not dissipation (LeVeque 2002). The L₁ norm measures both errors, dispersion and diffusion.

To further investigate this, we consider the decay rate of the magnetic field in the wave that is sensitive only to the dissipation error of the scheme (numerical diffusion), not the dispersion of the wave. The decay rate ${\tau }_{{\rm{D}}}$ of B_y is defined as

$$\begin{eqnarray}&&{\tau }_{{\rm{D}}}({B}_{y})=-\displaystyle \frac{1}{{t}_{\mathrm{end}}}\mathrm{log}\left(\displaystyle \frac{\delta {B}_{y}(t={t}_{\mathrm{end}})}{\delta {B}_{y}(t=0)}\right),\end{eqnarray} \tag{ 60 }$$

where δB_y(t) is the root-mean-square magnetic field fluctuation in the frame parallel to the wave vector. We plot the decay rate of the fast, slow, and Alfvén waves in Figure 3 between the initial conditions and time t = 2. The time was chosen so that all wave families have been advected through the domain at least once. It is important to evaluate the decay time at a late enough time, as the rms values fluctuate with time, indicative of additional waves being present in the test setup. Fifth-order convergence of the decay rate is observed for all three remaining MHD waves.

Zoom In Zoom Out Reset image size

This confirms that the error introduced by the second-order CT scheme is indeed in the form of wave dispersion only, not dissipation, i.e., the price paid for the computationally cheap CT method from Ryu et al. (1998) is in the form of the magnetic field shape, but not magnetic energy conservation. This view is supported by results in the Alfvén wave and magnetic field loop advection tests shown later in the paper. Nonetheless, a comparison of the absolute value of the L₁ error with Stone et al. (2008 their Figure 32) yields that CT-WENO5 increases the effective resolution compared to ATHENA in three dimensions by roughly a factor of 2, even for linear MHD waves propagating with an oblique angle to the grid in three dimensions. The L₁ error is significantly smaller even at N = 16, so the statement is true despite the second-order convergence for MHD waves.

Polarized Alfvén waves are important in heating the solar corona due to instabilities and have been studied for many decades (e.g., Goldstein 1978; Del Zanna et al. 2001; Ruderman & Simpson 2005; De Pontieu et al. 2007). As a numerical test, polarized Alfvén waves can be used to evaluate the fidelity of the CT scheme (Tóth 2000).

Following Gardiner & Stone (2005), we set ${\boldsymbol{U}}=(1,0,{v}_{y},{{v}_{z},1,{B}_{y},{B}_{z},0.1)}^{T}$, with ${v}_{y}={B}_{y}\,=0.1\sin (2\pi x)$, ${v}_{z}={B}_{z}\,=0.1\cos (2\pi x)$. In two dimensions, we rotate U by the second Euler angle $\theta ={\tan }^{-1}(2)\approx 63\buildrel{\circ}\over{.} 5$. The computational domain has $2N\times N$ zones and covers ${L}_{\mathrm{Box},x},{L}_{\mathrm{Box},y}=\sqrt{5},\sqrt{5}/2$. In three dimensions, the wave is rotated similarly to Section 4.1. The computational domain of 2N × N × N zones covers $({L}_{x},{L}_{y},{L}_{z})=(3,3/2,3/2)$. Gardiner & Stone (2005) did not find an instability using these parameters. We note that the analytical rms at t = 0 is given as $0.1/\sqrt{2}$ for fluctuating quantities, while for all other quantities the initial rms is zero. Thus, the former errors dominate the simulation, while the latter quantities are subject to truncation errors from the rotation operation.

On the top left of Figure 4, we show the L₁ error norm in 2D (blue) and 3D (green) over the number of zones after the wave has traveled for a five-wave period (${c}_{{\rm{A}}}=1$) in a frame parallel to the wave vector. The scheme converges at second order in L₁. This is in line with the result from the previous linear wave convergence test.

Zoom In Zoom Out Reset image size

We also show the decay rate convergence as dashed lines in Figure 4. However, this time, oscillations in the rms of the state vector components can lead to negative decay rates even at rather late times. In Figure 4, bottom, we show these rms fluctuations over time for $64\times {32}^{2}$ zones in three dimensions for the components that carry a fluctuation in the initial conditions. Dissipation in the solution is observed as asymmetry in the rms at late times. To obtain a robust measure for the wave decay, one may consider the rms of the transverse magnetic flux in the wave (red line):

$$\begin{eqnarray}&&\delta {F}_{{\rm{B}},\mathrm{tv}}=\mathrm{rms}\left(\sqrt{{\left({B}_{y}{v}_{y}\right)}^{2}+{\left({B}_{z}{v}_{z}\right)}^{2}}\right)\end{eqnarray} \tag{ 61 }$$

For this quantity, fifth-order convergence is observed at all times in Figure 4, top left.

In Figure 4, top right, we show the z-component of the magnetic field in two dimensions after advecting for five wave periods for resolutions of $N\in [8,16,32,64,128]$ in red, blue, green, purple, and orange, respectively. Only the lowest resolution shows some dissipation; dispersion is observed for the two lowest resolutions. This is a significant improvement over the results shown in Tóth (2000) for the same CT scheme. It also compares favorably with WOMBAT's TVD+CTU implementation (see Mendygral et al. 2017) and matches results from the ATHENA solver (Gardiner & Stone 2008, their Figure 4), despite the simpler CT scheme.

We compute shock tube number 2a from Ryu & Jones (1995) with WOMBAT's TVD+CTU and WENO5 solver. The shock tube features all four MHD modes. The left-hand state vector is ${{\boldsymbol{U}}}_{{\rm{L}}}={(1.08,1.2,0.01,2/\sqrt{4\pi },3.6/\sqrt{4\pi },2/\sqrt{4\pi },0.95)}^{T}$. The right-hand state vector is ${{\boldsymbol{U}}}_{{\rm{R}}}\,=(1,0,0,0,2/\sqrt{4\pi },{4/\sqrt{4\pi },2/\sqrt{4\pi },1)}^{T}$. We plot the TVD+CTU result with 10⁴ zones using the red line and the WENO5 result with 256 zones as black circles. The evolved state vector components at t = 0.2 for the one-dimensional computation are found in Figure 5; those in three dimensions are in Figure 6. For three dimensions, we rotated the shock normal along k = (1, 1, 1)^T in a 128³ zone simulation. The primitive variables alongside the magnetic field angle $\phi ={\arctan }^{-1}({B}_{y}/{B}_{z})$ in degrees are shown in Figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We compute shock tube number 4d from Ryu & Jones (1995) with WOMBAT's TVD+CTU and WENO5 solver. The left-hand state vector is ${{\boldsymbol{U}}}_{{\rm{L}}}={(1,0,0,0,0.7,0,0,1)}^{T}$. The right-hand state vector is ${{\boldsymbol{U}}}_{{\rm{R}}}={(0.3,\mathrm{0.0.1},0.7,1,0,0.2)}^{T}$. We plot the TVD+CTU result with 10⁴ zones using the red line and the WENO5 result with 256 zones as black circles. The evolved state vector components at t = 0.16 are shown in Figure 7.

Zoom In Zoom Out Reset image size

The MHD shock tube from Brio & Wu (1988a) has the initial conditions ${{\boldsymbol{U}}}_{{\rm{L}}}={(0.125,0,0,0,0.75,1,0,1)}^{T}$, ${{\boldsymbol{U}}}_{{\rm{R}}}=(0.125,{0,0,0,0.75,-1,0,0.1)}^{T}$, with $\gamma =2$. In Figure 8, we show the solution at t = 0.08 with 400 zones as black squares and the solution from TVD+CTU with 10⁴ zones as the red line.

Zoom In Zoom Out Reset image size

We note that this problem starts with a degeneracy in the magnetic field when computing the eigenvectors of the initial conditions in the single zone at x = 0.5. Following Brio & Wu (1988a), we resolve this degeneracy by setting ${B}_{y}={B}_{z}=1/\sqrt{2}$. However, in this particular problem, B_z = 0 at all times. Thus, the test exposes this choice in the eigenvectors. Setting B_z = 0 for this problem only would remove the oscillations in Figure 8. Balbás & Tadmor (2006) have shown that the issue can also be fixed with global smoothness indicators. In real world applications in two or more dimensions, this is not an issue. H. Jang et al. (2019, in preparation) found that the oscillations disappear with WENO3, i.e., in lower order codes, increased dissipation leads to a constant solution.

The Shu–Osher shock tube simulates the interaction of a shock with a smooth flow (Shu & Osher 1989). Greenough & Rider (2004) used it to argue that second-order PPM methods give results comparable to third-order WENO methods and are thus computationally more efficient. It exposes that the standard WENO5 weights are at best third-order accurate in critical points (Borges et al. 2008).

The initial conditions contain a shock tube with left and right states ${{\boldsymbol{U}}}_{{\rm{L}}}={(3.857143,2.629369,0,0,0,0,0,31/3)}^{T}$ and ${{\boldsymbol{U}}}_{{\rm{R}}}={(1+0.2\sin (5x),0,0,0,0,0,0,1)}^{T}$ on a domain of size L_x = 10. The shock is located at x = −4. The shock tube is evolved until t = 2.

In Figure 9, we show the result from the WENO5 (green circles) and WENO5-Z simulations (black squares) with 300 zones and TVD+CTU with 10,000 zones (red line). Clearly, the WENO-Z result traces the complex flow behind the shock more accurately than the WENO5 result.

Zoom In Zoom Out Reset image size

Another variant of the test, run on a domain of [−1, 1] and with a factor π in the sine function of the initial conditions, is shown in Figure 9, bottom. The result can be directly compared to Felker & Stone (2018), their Figure 5. WENO5 roughly performs as well as their PPM+RK4 scheme, while WENO5-Z traces the complex shock wave structure better than PPM+RK4. This is in line with expectations from the WENO literature (e.g., Borges et al. 2008).

Following McCorquodale & Colella (2011), we simulate the time evolution of a Gaussian acoustic pulse in two dimensions. In a periodic domain of size one, we set

$$\begin{eqnarray}\rho (r)=\left\{\begin{array}{ll}\rho 0+(\delta {\rho }_{0}){e}^{-16{r}^{2}}\cos (\pi r) & \mathrm{if}\quad r\leqslant \tfrac{1}{2}\\ \rho 0 & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 62 }$$

$$\begin{eqnarray}&&p={\left(\displaystyle \frac{\rho }{\rho 0}\right)}^{\gamma }\end{eqnarray} \tag{ 63 }$$

with ${\rho }_{0}=1.4$, $\delta {\rho }_{0}=0.14$, γ = 1.4. The simulation ran until t = 0.24.

We show the resulting density in Figure 10 at t = 0 (left) and at t = 0.24 (right). The color scale ranges from 0.136 to 1.54; we add contours at 1.375, 1.4, 1.425, 1.45, 1.475, 1.5, and 1.525.

Zoom In Zoom Out Reset image size

To compare directly with McCorquodale & Colella (2011), we run a convergence study of the density differences at time 0.24. To this end, we compute the density error using cubic spline interpolation to obtain a reference solution from a run with 1024² zones. We obtain second-order convergence.

The noise inherent in an implementation can be exposed in the two-dimensional implosion test (Liska & Wendroff 2003). The setup leads a series of shock waves, which are intersecting many times in the simulation. The problem is highly nonlinear and thus not very useful to compare algorithmic fidelity. There is simply no convergence to a universal solution that results could be compared to. However, the nonlinearity ensures that computational noise gets amplified very quickly and solutions diverge strongly from correct, less noisy results. Using this test, we found noise injected by the compiler at the $\sqrt{{10}^{-16}}$ level in the eigenvector calculation (see the Appendix) by comparing the solution to TVD+CTU. This may be a word of warning for running Eulerian methods with 4 byte floating point precision, where numerical noise itself resides at 10⁻⁸ and thus will influence this test significantly. Many real world applications are strongly nonlinear and thus may not be correct with 4 byte variables.

We set ${\boldsymbol{U}}={(1,0,0,0,0,0,0,1)}^{T}$ everywhere in a periodic domain with size ${L}_{x}={L}_{y}=0.3$ and 200² zones. When x + y < 0.15, we set ${\boldsymbol{U}}=(0.125,0,0,0,0,0,0,{0.14)}^{T}$. We show the resulting density in Figure 11 at t = 0.2, 0.3, 0.5, and 0.7 (top left to bottom right) on a color scale ranging from 0.4 to 1.2. We notice that the symmetry in the simulation is not visibly broken, and the results are reasonably similar to the tests shown in Liska & Wendroff (2003). We show the vorticity of the test at t = 0.3 in Figure 12, where no noise is visible. In particular, the result is symmetric along the x = y axis down to single zones.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The implosion test was used by Sijacki et al. (2012) to compare results between a traditional SPH and a Lagrangian finite volume method. Both show more computational noise than our Eulerian method due to particle motion. The same seems to be true for the low-order DG scheme presented in Mocz et al. (2014). There are sizable differences in the shape of the "bird" in the lower left corner between all these codes and our results. In particular, interfaces in the results of the Lagrangian FV and DG schemes differ from our WENO5 simulation. Numerical noise seems to trigger Raleigh–Taylor instabilities in the Lagrangian simulation that are absent in the Eulerian simulations. This noise is clearly visible in the vorticity maps (their Figure 3). Runs at much higher resolution show that indeed these interfaces do become unstable eventually. However, due to the nonlinearity of the flow, the high-resolution Eulerian result differs significantly from both low-resolution runs in that all instabilities grow faster in the high-resolution Eulerian runs. As mentioned above it is unclear which result is closer to the "true solution" in this test, and we conclude that numerical noise affects the growth of instabilities in highly nonlinear solutions of the hydrodynamic equations, which is of course not a new result at all (e.g., McNally et al. 2012).

The Orszag–Tang vortex is a standard MHD test that mimics the transition of a smooth flow into 2D turbulence (Orzang & Tang 1979; Tóth 2000; Stone et al. 2008). The initial conditions in a unit domain are ${\boldsymbol{U}}=(25/(36\pi ),{-\sin (2\pi y),\sin (2\pi y),0,{B}_{x},{B}_{y},{B}_{z}5/(12\pi ))}^{T},$ where the magnetic field components are obtained from a vector potential with the z-component ${A}_{z}=({B}_{0}/4\pi )\cos (4\pi x)\,+({B}_{0}/2\pi )\cos (2\pi y)$, with ${B}_{0}=1/\sqrt{4\pi }$.

We show in Figure 13, top left to bottom right, the density, pressure, velocity magnitude, and magnetic field strength at time t = 0.5 at a resolution of 192² zones with WENO5-Z. In Figure 14, we show the vortex at t = 1 run with 500² zones. Some schemes produce negative pressures at these late times. The results compare very well to the high-order code shown in Felker & Stone (2018). In particular, we see more structure in the central vortex and no break in symmetry. In Figure 15, we show the pressure in two slices through the domain at $t=0.5,y\,=-0.1875$ (top) and t = 0.5, y = 0.073 (bottom). A high-resolution run with TVD+CTU is included as a reference solution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 16, we compare the pressure at y = −0.1875, −0.5 < x < 0, and t = 0.5 of WENO5-Z with 128² zones (green squares) with TVD+CTU (blue circles) with 256² zones. The reference solution with 1024² zones is plotted by the red line. Apparently, WENO5-Z has fidelity comparable to or better than that of TVD+CTU at double the resolution in this test. In particular, the blip at $x=-0.45$ is better resolved in the WENO5-Z solution. A comparison with Figure 15 shows that the blip becomes fully resolved at 192².

Zoom In Zoom Out Reset image size

In Figure 17, we compare Schlieren images (Hooke 1665; Helzel et al. 2011) of the density in TVD+CTU runs at 512² and 256² zones (left, middle), with the WENO5-Z run at 128². Following Guo & Zhang (2004), we obtained a synthetic brightness value by adding the red and yellow channels with a Gladstone–Dale constant of ${K}_{\mathrm{GD}}=1$, L = 1, and a contrast parameter of C = 0.075. Schlieren images visualize the gradient of the solution, which makes the Orszag–Tang vortex test more discriminating between schemes. We observe that all main features of the low-resolution TVD+CTU simulation are reproduced in the WENO5-Z simulation at half that resolution. We conclude that WENO5-Z doubles the effective resolution compared to TVD+CTU in this test.

Zoom In Zoom Out Reset image size

The MHD rotor test simulates a rapidly rotating disk in two dimensions, which produces shocks and rarefaction waves as the disk expands due to centrifugal forces (Stone & Norman 1992; Balsara & Spicer 1999a; Stone et al. 2008). The solution needs to remain symmetric, and no spurious oscillations may occur. We set initial conditions with γ = 5/3 in a unit domain of P = 1, ${\boldsymbol{B}}={(5/\sqrt{4\pi },0)}^{T}$, ${v}_{0}=2$, ${r}_{0}=0.1$, ${r}_{1}=0.115$, $f(r)=\tfrac{{r}_{1}-r}{{r}_{1}-{r}_{0}}$, and (Tóth 2000)

$$\begin{eqnarray}\rho =\left\{\begin{array}{ll}10 & \iff \ r\lt {r}_{0}\\ 1+9f(r) & \iff \ {r}_{0}\leqslant r\leqslant {r}_{1}\\ 1 & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 64 }$$

The velocity components are

$$\begin{eqnarray}{v}_{x}=\left\{\begin{array}{ll}-y\tfrac{{v}_{0}}{{r}_{0}} & \iff \ r\lt {r}_{0}\\ -f(r)y\tfrac{{v}_{0}}{{r}_{0}} & \iff \ {r}_{0}\leqslant r\leqslant {r}_{1}\\ 0 & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 65 }$$

$$\begin{eqnarray}{v}_{y}=\left\{\begin{array}{ll}-x\tfrac{{v}_{0}}{{r}_{0}} & \iff \ r\lt {r}_{0}\\ -f(r)x\tfrac{{v}_{0}}{{r}_{0}} & \iff \ {r}_{0}\leqslant r\leqslant {r}_{1}\\ 0 & \mathrm{otherwise},\end{array}\right.\end{eqnarray} \tag{ 66 }$$

We show the test result with WENO5 and 200² zones at time t = 0.15 in Figure 18. A slice through the rotor at t = 0.5 and $y=-0.1875$ (top) and y = 0.073 are shown in Figure 19. We show the WENO5 solution with black squares, and TVD+CTU with 1024² zones with the red line and with 400² zones with blue dots. No asymmetries or oscillations are visible; the WENO5 solution at 200² resolves the small features comparably to or better than TVD+CTU. A by-eye comparison with published results from ATHENA (Stone et al. 2008, their Figure 26), yields that WENO5 at half resolution is not quite as resolved as ATHENA at full resolution, likely because our CT scheme is only second order in space. There is some ringing visible in the density maps, which could likely be alleviated using global weights. We leave the improvement to future work.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The double Mach reflection test simulates the evolution of a Mach 10 shock at an oblique 60° angle to the grid with reflecting boundaries on the bottom x-axis, continuous upper x-boundary, and open y-boundaries (Woodward & Colella 1984; Stone et al. 2008). Here, the adiabatic index γ = 1.4 for air, the state ahead of the shock is ${\boldsymbol{U}}={(1.4,0,0,0,0,0,0,1)}^{T}$, and the domain size is 3.25 × 1. Because the shock is reflected off the bottom wall, a second weak shock and a jet form. In the upper boundary, the evolution of the shock is followed analytically, i.e., ${x}_{\mathrm{shock}}=1/6+(1+20t)$. As the shock in the setup is only one cell wide, the initial conditions contain an error, which is also continuously injected at the upper boundary. A one-cell-wide shock is not a "natural" solution to the numerical scheme, similar to the common shock tube setup used in SPH code tests. The test foremost exposes the robustness of the numerical scheme that has to handle these errors by injecting numerical diffusion. This differs from SPH shock tube tests, where dissipation and conduction has to be explicitly added. This test produces negative pressures when run with the WENO5-Z scheme without protection floors or protection fluxes.

In Figure 20, we show the density from the test run with the standard WENO5 scheme at times t = 0.05, 0.1, and 0.2 (top to bottom). We also overlay 30 contours between 1.73 and 20.92, which can be directly compared to Woodward & Colella (1984), their Figure 4. In the top panel, the error of the initial shock is visible in the color scheme as two stripes with slightly lower density left of and parallel to the shock. Over time, another such dip forms, starting at the location of the shock at the upper y-boundary traveling toward the lower left.

Zoom In Zoom Out Reset image size

A zoom into a simulation with 1024 × 256 zones at t = 2 is shown in Figure 21 and can be directly compared to McCorquodale & Colella (2011), their Figure 7. In some codes, the reflection of the shock off the bottom y-boundary causes a carbuncle instability, because the Mach number is very high here. We do observe this instability with the standard WENO5 scheme, although only at high resolution.

Zoom In Zoom Out Reset image size

The growth of instabilities plays a crucial role for mixing and the injection of vorticity in astrophysical fluid flows (e.g., Sheardown et al. 2018). We evaluate the performance of the code using the Kelvin–Helmholtz (KH) test, where the kinetic energy of a slightly perturbed shear flow generates vorticity at the interface (e.g., Parker 1958; Miura & Pritchett 1982; Frank et al. 1996). This test has been used extensively in the literature to test hydrodynamic codes (e.g., Springel 2010; McNally et al. 2012; Hopkins 2015; Schaal et al. 2015; Beck et al. 2016). In particular, the KH problem was used to expose the artificial surface tension of traditional SPH methods (e.g., Hopkins 2013). It also exposes the velocity-dependent truncation error in Eulerian grid methods (Robertson et al. 2010). Our setup follows Lecoanet et al. (2016):

$$\begin{eqnarray}&&\rho =1+\displaystyle \frac{{\rm{\Delta }}\rho }{{\rho }_{0}}\displaystyle \frac{1}{2}\left[\tanh \left(\displaystyle \frac{z-{z}_{1}}{a}\right)-\tanh \left(\displaystyle \frac{z-{z}_{2}}{a}\right)\right]\end{eqnarray} \tag{ 67 }$$

$$\begin{eqnarray}&&{v}_{x}={v}_{\mathrm{flow}}\left[\tanh \left(\displaystyle \frac{z-{z}_{1}}{a}\right)-\tanh \left(\displaystyle \frac{z-{z}_{2}}{a}\right)-1\right]\end{eqnarray} \tag{ 68 }$$

$$\begin{eqnarray}&&{v}_{z}=A\sin (2\pi x)\left[\exp \left(-\displaystyle \frac{{\left(z-{z}_{1}\right)}^{2}}{{\sigma }^{2}}\right)+\exp \left(-\displaystyle \frac{{\left(z-{z}_{2}\right)}^{2}}{{\sigma }^{2}}\right)\right]\end{eqnarray} \tag{ 69 }$$

$$\begin{eqnarray}&&P={P}_{0},\end{eqnarray} \tag{ 70 }$$

with the parameters A = 0.01, ${P}_{0}=10$, a = 0.05, $\sigma =0.4$, ${z}_{1}=0.5$, z₂ = 1.5, ${u}_{\mathrm{flow}}=1$, and ${\rm{\Delta }}\rho /{\rho }_{0}={\rho }_{0}=1$. We note that the instability is not resolved in our runs with this choice of parameters, i.e., the instability grows too slowly due to numerical effects. Simulations are run in a periodic domain with ${L}_{y}=2,{L}_{x}=1$.

In Figure 22, we show the density of four simulations at t = 2. The top-left to bottom-right panels show WENO5-Z without boost with 512 × 1024 zone resolution, WENO5 with 32 × 64 zones without a velocity boost, WENO5 with 32 × 64 zones boosted by v_x = 50, and TVD+CTU with 128 × 256 zones without boost.

Zoom In Zoom Out Reset image size

The TVD+CTU run does not develop a vortex roll due to numerical diffusivity, even at twice the resolution from the WENO5 runs. In contrast, the low-resolution WENO5 runs develop a vortex similar to the high-resolution simulation (see Lecoanet et al. 2016). This shows that WENO5 resolves instabilities much better than the second-order TVD+CTU, even at half the grid resolution. WENO5 is also much less susceptible to advection errors. The boost of v_x = 50 (!) barely changes the result. Some differences are visible at the tip of the roll, because the truncation error is not Galilean invariant. Nonetheless, we expect significant improvements in the growth of instabilities in cosmological simulations, where advection is commonplace and may be a major source of diffusivity for second-order Eulerian methods (Springel 2010).

The Gresho vortex (Gresho & Chan 1990; Liska & Wendroff 2003; Springel 2010) is a two-dimensional rotating structure in hydrodynamic equilibrium that evaluates the angular momentum and vorticity conservation of the code. The test becomes very challenging for Eulerian codes, when the vortex is advected to the grid. It exposes how the increase in truncation error affects angular momentum conservation. Instructive discussions can be found in Miczek et al. (2015) and Wongwathanarat et al. (2016).

In a two-dimensional unit domain, we set

$$\begin{eqnarray}&&{\boldsymbol{U}}={\left(1,{v}_{x}(r),{v}_{y}(r),\mathrm{0,\; 0},\mathrm{0,\; 0},P(r)\right)}^{T},\end{eqnarray} \tag{ 71 }$$

with the radius r and¹⁰

$$\begin{eqnarray}&&{v}_{x}(r)=-v(r)\cos (\phi )+{v}_{\mathrm{boost}},\end{eqnarray} \tag{ 72 }$$

$$\begin{eqnarray}&&{v}_{y}(r)=-v(r)\sin (\phi ),\end{eqnarray} \tag{ 73 }$$

$$\begin{eqnarray}&&\phi ={\tan }^{-1}\left(\displaystyle \frac{y}{x}\right)\end{eqnarray} \tag{ 74 }$$

$$\begin{eqnarray}v(r)=0.4\pi \left\{\begin{array}{ll}5r & \iff \,r\leqslant 0.2\\ 2-5r & \iff \,r\lt 0.4\\ 0 & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 75 }$$

$$\begin{eqnarray}P(r)=\left\{\begin{array}{ll}{p}_{0}+\tfrac{25}{2}{r}^{2} & \iff \,r\leqslant 0.2\\ {p}_{0}+\tfrac{25}{2}{r}^{2}+4\left(1-5r+\mathrm{log}\left(\tfrac{r}{0.2}\right)\right) & \iff \,r\lt 0.4\\ {p}_{0}-2+4\mathrm{log}(2) & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 76 }$$

$$\begin{eqnarray}&&{p}_{0}=\displaystyle \frac{{\left(0.4\pi \right)}^{2}}{\gamma {M}_{\max }^{2}}-\displaystyle \frac{1}{2}.\end{eqnarray} \tag{ 77 }$$

The tests shown in Springel (2010) correspond to ${p}_{0}=5$, a relatively high maximum Mach number of ${M}_{\max }\approx 0.41$.

We run a convergence study of the test with the WENO5-Z algorithm using three different boost velocities in the x-direction: $0,1,3$. The results are shown in Figure 23 where we show the L₁ error over resolution. Blue, green, and purple lines correspond to boost velocities of $0,1,3$, respectively. We add a dotted line scaling with ${N}^{-1.4}$, which corresponds to the scaling obtained with AREPO in Springel (2010). In this regime, our results compare quite favorably to the results obtained with lower order FV codes. In particular, the error remains small even with the velocity boost at most resolutions and is smaller than the errors observed with ATHENA at all times.

Zoom In Zoom Out Reset image size

The boost mostly affects resolutions below 100², where the vortex becomes underresolved. Below 20², the L₁ error with velocity boost is reduced again. We note that the error changes depending on how the center of the vortex is placed on the grid, likely because the initial conditions contain a discontinuity at r = 0.2. As discussed in Muñoz et al. (2013), the sampling of the discontinuity affects the convergence properties of the simulations.

To illustrate this further, we show in Figure 24 (left to right) the angular velocity v_ϕ, pressure, and vorticity $\omega (r)={\rm{\nabla }}\times {\boldsymbol{v}}$ over radius. We compute

$$\begin{eqnarray}&&v(\phi )=r\displaystyle \frac{{{xv}}_{y}-{{yv}}_{x}}{{x}^{2}+{y}^{2}},\end{eqnarray} \tag{ 78 }$$

and the vorticity from the standard FD operator. Analytic profiles are shown as red lines, the initial conditions as red dots, the simulation at t = 3 without boost as black dots, and the simulation with velocity boost of v_x = 3 as blue dots. From top to bottom, we show resolutions of ${16}^{2},{32}^{2}$ and 64² zones. The results for the static compare well with the results of the low-order DG scheme on a rectangular grid shown in Mocz et al. (2014), their Figure 5. Low-order DG codes on polar grids, however, show better conservation properties in vorticity (e.g., Velasco Romero et al. 2018, their Figure 3). The run with 32² zones and velocity boost shows the largest scatter, in line with the increase in L₁ error observed before. At 16² zones, the profiles show less scatter, but are only sampled with 15 unique values below r = 0.4 in the ICs. At this low resolution, the problem is underresolved, and it stands to reason that the reduced scatter in the evolved boosted profiles is a result of the poor sampling. In particular, the discontinuous peak of the v_ϕ profile is not sampled directly, which should lead to additional error (Springel 2010; Muñoz et al. 2013). We speculate that the pronounced increase in L₁ error around 32³ zones is a feature of the FD WENO5 reconstruction without cross-terms. Indeed, Liska & Wendroff (2003) found a relatively high L₁ error in total kinetic energy for their WENO5 scheme at low resolutions.

Zoom In Zoom Out Reset image size

To fully expose the fidelity in momentum conservation, we show the normalized velocity at t = 1 with 40² zones in Figure 25. From left to right, we show runs with decreasing ${M}_{\max }={10}^{-1},{10}^{-2},{10}^{-3}$. The vortex clearly loses integrity at low Mach numbers. This result can be compared directly to Wongwathanarat et al. (2016), their Figure 3. Following the discussion in Miczek et al. (2015), this result is likely due to the scaling of the Roe fluxes, and the missing cross-terms in our FD WENO scheme.

Zoom In Zoom Out Reset image size

We conclude that the stationary solution with WENO5 is competitive with the AREPO and ATHENA solutions presented in Springel (2010), but not with high-order DG codes (Velasco Romero et al. 2018) or FV codes on polar grids (Wongwathanarat et al. 2016). At relatively high Mach numbers, the error in the boosted solutions is comparable with second-order static Cartesian mesh solutions at intermediate resolutions and approaches the Lagrangian error at resolutions above 100² zones.

The velocity-dependent truncation error of Eulerian grid methods can be clearly exposed in the density advection test (Robertson et al. 2010; Hopkins 2015). In two dimensions, a density square in hydrodynamic equilibrium with the surrounding medium is advected supersonically through a unit domain: ${\boldsymbol{U}}={(1,100,50,0,0,2.4)}^{T}$, except $\rho =4$, when $-0.25\,\lt x\lt 0.25$ and $-0.25\lt y\lt 0.25$. Here, γ = 5/3, thus c_s = 2, so the advection is supersonic ($M\gt 100$) with respect to the reference frame of the mesh. Although this test is trivial with Lagrangian schemes such as SPH, meshless finite volume, and moving finite volume schemes, it presents a very serious challenge to Eulerian grid methods.

We show the result of the test run with TVD+CTU, WENO5, and WENO5-Z in Figure 26 at time t = 10. We also include a WENO5-Z run with CFL = 0.2. The convergence rate is shown in Figure 27.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Due to the high velocity relative to the grid, the solution is smeared out significantly by the TVD+CTU algorithm, as expected from a second-order method. WENO5 performs a bit better, although still with significant dispersion. This is because the edges of the density square are critical points in the flow, where the classical WENO5 weights are only third-order accurate. WENO5-Z performs significantly better, with the density square being preserved, but strongly distorted. Quantitatively, WENO5-Z shows an improvement in L₁ error by roughly a factor of 2 at late times.

These results can be directly compared to the results shown with the discontinuous Galerkin (DG) code TENET in Schaal et al. (2015). In terms of L₁ error, WENO5 performs roughly like a second-order DG code. WENO5-Z performs better than second-order DG, but worse than third-order DG. We note that increasing the convergence order in a DG scheme reduces the time step by a factor of 2. Thus, the time step of the DG3 scheme is 20 times smaller than in our WENO5-Z run (compare their Equation (24), CFL = 0.2, with our Equation (17), CFL = 0.8, in three dimensions). Thus, these DG schemes are more computationally expensive than WENO5, where the time step does not change with the order of the scheme.¹¹ DG schemes also require more memory, because the full polynomial has to be stored in every zone in one way or another. We ran the WENO5-Z simulation with a reduced CFL number of 0.2 to approach the computational cost of a DG3 algorithm. The distortion of the square is reduced further, but the L₁ error remains roughly the same.

This result suggests running strongly advection-dominated problems with WENO5-Z, not with the classical WENO5 method. As the method is significantly less robust, this results in a trade-off between more protection fluxes and reduced advection error for complex problems. It will depend on the problem under consideration on which approach delivers better results.

The field loop advection test is instrumental for the development and testing of CT schemes for the magnetic field evolution in Eulerian codes (Gardiner & Stone 2008). This is a difficult test for the CT scheme of the code, as only the electric field and magnetic field are moving through the domain. In applications, this will rarely be the case.

In two dimensions, we set up a periodic domain of size ${L}_{x}\times {L}_{y}=2\times 1$ with ${\boldsymbol{U}}={(1,2,1,0,{B}_{x},{B}_{y},0,1)}^{T}$ and a vector potential on the interfaces with zero x and y components, but

$$\begin{eqnarray}&&{A}_{z}={A}_{0}\left({r}_{c}-r\right),\end{eqnarray} \tag{ 79 }$$

where ${A}_{0}={10}^{-3}$, r_c = 0.3. The magnetic field on the interfaces is then obtained as the curl of the vector potential on the faces. In Figure 28, we show the resulting magnetic field energy density of the initial conditions (left) and the simulation at time t = 2 with TVD+CTU (middle) and WENO5 (right) with 128 × 64 zones. The WENO5 scheme keeps the shape of the loop well while it is advected through the box. The shape of the loop is comparable to second-order CT schemes (Stone et al. 2008; Lee 2013). In contrast, the TVD+CTU loop is barely visible; most of the magnetic energy has dissipated into heat. To quantify this, we show the time evolution of the normalized magnetic energy in three runs with 32, 64, and 128 zone base resolution in Figure 29. We mark 98% by the dotted horizontal line.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

WENO5 retains more than 90% of the magnetic energy of the field loop, even at the lowest resolution run, while TVD+CTU loses more than half of the magnetic energy until that time, even at the highest resolution. The conservation of magnetic energy also compares favorably with the results shown in Lee (2013) and Felker & Stone (2018). Even though our CT scheme is much less elaborate than theirs, the magnetic field fluxes are fifth-order accurate and thus improve energy conservation. Again, WENO5 roughly doubles the effective resolution of the simulation, i.e., our run with 64² zones conserves magnetic energy roughly like ATHENA with 128² zones. These results are in line with our findings from the wave convergence tests, where wave diffusion is fifth-order accurate, but wave dispersion is second-order accurate. The structural integrity of our field loop is comparable to the results in Felker & Stone (2018), which is third order in space and fourth order in time. As expected, WENO5 with CT conserves magnetic energy better than Felker & Stone (2018), due to the fifth-order fluxes. The scheme also compares well with the FDs schemes shown in Mignone et al. (2010). This shows the advantage of the CT approach over divergence-cleaning methods, despite the second-order convergence in the linear wave advection tests. The fourth-order FV CWENO code presented in Verma et al. (2019) shows even better conservation properties. The CWENO code retains 98% of magnetic energy at t = 100 for 128² zones. These results place our FD WENO5 scheme in between high-order FV WENO codes and low-order FV codes or high-order codes without CT. Considering the usual trade-off in efficiency and fidelity, our simple but fast CT approach sacrifices some magnetic field structure for computational efficiency, but not magnetic energy conservation. In three dimensions, the loop is rotated around the y-axis by $\phi =\arctan (0.5)$ (Gardiner & Stone 2008). We set up the magnetic vector potential as

$$\begin{eqnarray}&&{A}_{x}=-{A}_{0}\displaystyle \frac{\left({r}_{c}-{r}^{\star }\right)}{\sqrt{5}}\end{eqnarray} \tag{ 80 }$$

$$\begin{eqnarray}&&{A}_{y}=0\end{eqnarray} \tag{ 81 }$$

$$\begin{eqnarray}&&{A}_{z}=-{A}_{0}\displaystyle \frac{2\left({r}_{c}-{r}^{\star }\right)}{\sqrt{5}},\end{eqnarray} \tag{ 82 }$$

where r^⋆ is the radius obtained from the coordinates

$$\begin{eqnarray}{x}^{\star }=\left\{\begin{array}{ll}\tfrac{(2x+z)+2}{\sqrt{5}} & \iff \ x\lt -\sqrt{5}\\ \tfrac{(2x+z)-2}{\sqrt{5}} & \iff \ x\gt -\sqrt{5}\end{array}\right.\end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray}&&{y}^{\star }=y\end{eqnarray} \tag{ 84 }$$

$$\begin{eqnarray}&&{z}^{\star }=\displaystyle \frac{-x+2z}{\sqrt{5}}.\end{eqnarray} \tag{ 85 }$$

A rendering of the resulting magnetic energy at time t = 2 is shown in Figure 30, i.e., after advecting the loop twice through the grid. Following Gardiner & Stone (2008), we impose a threshold in the projection at ${B}^{2}={10}^{-7}$ to show the shape of the loop edge. Again, the structure of the loop is comparable with the rendering shown in Gardiner & Stone (2008).

Zoom In Zoom Out Reset image size

This test models the self-similar evolution of a point-like thermal detonation (von Neumann 1942; Sedov 1946; Taylor 1950). Despite its unpleasant historical context, it also is a model for early supernova explosions and is useful to expose the effects of the regular grid on the evolution of a strong spherical shock wave propagating into a thin medium with negligible pressure. We set ${\boldsymbol{U}}=(1,0,0,0,{10}^{-8},0,0,{10}^{-5})$ everywhere in a three-dimensional computational domain with ${L}_{x}={L}_{y}={L}_{z}=1$. We inject a unit energy ${E}_{1}=1$ in the center of the box and distribute it using a Gaussian with FWHM $\sigma =0.5\,\mathrm{dx}$, where dx is the zone size. Our simulation uses γ = 5/3 and N = 64³ zones. The self-similar analytic solution of the problem can be found in Landau & Lifshitz (1966), with β = 1.15.

In the left panel of Figure 31, we plot the radial profile of the shock at t = 0.06: the analytic solution in the red line and every 20th zone of the simulation by a black dot. On the right of Figure 31, we show a slice through the simulation at $z=0,t=0.06$ with colors on a log scale from 10⁻³ to 2 and contours at 0.01, 0.05, 0.1, 0.2, 0.5, and 2. The imprint of the grid on the explosion is clearly visible at radii below 0.3 and comparable to the DG3 method reported in Schaal et al. (2015), but better than the lower order WENO method presented in Feng et al. (2004), likely due to our higher order time integrator. Although these results seem disappointing at first, it is worth remembering that the DG3 scheme results used a 20 times smaller time step than ours and are an order of magnitude more expensive to compute.

Zoom In Zoom Out Reset image size

In this test, Lagrangian schemes perform much better than our Eulerian scheme, due to the adaptive sampling (e.g., Springel 2010).

The magnetized blast wave simulation tests code performance in the low-β regime, i.e., when shocks are evolved in the presence of strong magnetic fields. Following Londrillo & Del Zanna (2000), we set ${\boldsymbol{U}}={(1,0,0,0,{B}_{0}/\sqrt{2},{B}_{0}/\sqrt{2},0,P)}^{T}$, with $P(r\lt {r}_{c})=1$, $P(r\gt {r}_{c})=100$, ${B}_{0}=10$, and r_c = 0.125. It follows that $\beta =0.02$ in the medium ahead of the shock. The domain is ${L}_{x}\times {L}_{y}\times {L}_{z}=1\times 1.5\times 1$. Slices through density, pressure, v², and B² from a WENO5 simulation at time t = 0.2 and z = 0 with resolution of 200 × 300 × 200 zones are shown in Figure 32. The imprint of the magnetic field that is oriented along ${(1,1,0)}^{T}$ on the shock is clearly visible. We note that this test evaluates the robustness of the algorithm, and indeed, WENO5 required protection floors to complete the simulation. We leave a more elaborate implementation of protection fluxes to future work.

Zoom In Zoom Out Reset image size

All modern HPC systems feature a cache hierarchy to increase effective memory bandwidth to the CPUs (level 1, level 2, ...), with the fastest cache usually being the smallest due to silicon cost. To achieve good performance, it is desirable to divide the computational problem into subproblems that fit into the cache hierarchy, so that memory bandwidth is increased. This technique is called cache blocking. The optimal size will depend on the architecture of choice and the overhead in the algorithm, thus it has to be determined empirically. In WOMBAT, cache blocking is mitigated by tuning the patch size for a given architecture. Small patch sizes are desirable, because they enable more fine-grained load balancing.

We run a patch size optimization study with 27 nodes to saturate MPI communication on the Aries interconnect. In Figure 33, we show the performance of the WENO5 algorithm on the Broadwell architecture with 2 × 18 cores per node in million zones per second per node as a function of number of zones per dimension per patch. The problem size is at least 512³ zones, but varies by about 30% between runs, as the number of patches must be evenly divisible by the number of threads to not induce load imbalance in the threading. Colors correspond to runs with a varying number of MPI processes per node—blue: 2 MPI ranks/12 OpenMP threads; green: 4 ranks/9 threads; purple: 12 ranks/3 threads; orange: 18 ranks/2 threads; yellow: 36 ranks/1 thread.

Zoom In Zoom Out Reset image size

For all decompositions, the performance is maximized at 18³ zone patches, then drops slightly and increases again toward 70³ zone patches. This behavior differs from the TVD+CTU solver that peaks at 32³ and shows a strong drop in performance for a larger patch size (see green dashed line from Mendygral et al. (2017). The flattening of the data arrays increases the vector length of all loops in the algorithm and shifts the optimal cache block to smaller patch sizes in the WENO5 solver. It also leads to effective use of the hardware prefetcher at large patch sizes, so the cache blocking is effectively done in hardware. This is not possible in the TVD+CTU solver that is not flattened. With the side constraint of small patch sizes, we conclude that WENO5 performs at 0.6 million zones per node with 18³ zones per patch on Broadwell, with four MPI ranks per node. WENO5 is thus a factor of 7.5 slower than the TVD+CTU solver at the same resolution, but doubles the effective resolution. It follows that the WENO5 solver is more efficient than the second-order scheme, because increasing the resolution by a factor of 2 per dimension increases the runtime of the TVD+CTU solver by roughly a factor of 16: 2× per dimension and roughly 2× because the time step halves due to the factor 2 in ${\rm{\Delta }}x$ in Equation (17) (TVD+CTU uses a 1D criterion, but with CFL = 0.4 in 3D).

Our implementation uses vector-processing registers (SIMD) in modern CPUs extensively. Flattening the data arrays leads to large trip counts in the vector loops that are well cache-blocked. For every memory access, the virtual address has to be translated into a physical address. In the translation look-aside buffer (TLB), the memory management unit buffers page addresses for this translation from virtual to physical memory. On a standard system with 4 kb memory pages, our optimized vector loops can lead to a large miss rate in the TLB; the resulting page walk reduces performance by a factor of 2 (TLB thrashing). Thus, it is important to run WOMBAT with some form of huge pages enabled.

To demonstrate the performance gain from vectorization, we show the performance of the WENO5 solver on a 3D problem in GFlops on the Intel Broadwell architecture in Figure 34: in red, a single core without threads; in blue, a single core with whole node active (36 cores as 4 ranks with 9 threads). The performance of the run was measured by Cray Perftools.

Zoom In Zoom Out Reset image size

In Figure 35, we compare the fractional number of instructions generated by the Cray compiler for the Intel Broadwell and Intel Skylake architecture, given a largest vector length using the -vector0 and -preferred-vector-length compiler flags. The top panel shows the instruction for the Intel Broadwell architecture. Clearly, the compiler vectorizes virtually all of the code at the highest vector length. The bottom panel shows the same graph, but for the Intel Skylake architecture. Here, the compiler only vectorizes the complete code at 128 bit vector length, but chooses a mix of vectorization and scalar instructions at broader vector lengths. In particular, at 512 bit vectors, scalar instructions make up half of the program. We note that the version with the widest vectors is the fastest for both architectures, reaching about 20% of double precision peak performance. Thus, Skylake is still faster than Broadwell by about 500 MFlops.

Zoom In Zoom Out Reset image size

This shows that on some architectures, constraints other than vector length influence performance of the code. This is a good argument against using intrinsics or compiler pragmas to vectorize the whole program manually. Aside from the additional maintenance effort required to port the program to new architectures, the approach can actually increase execution time on architectures like Skylake. The better strategy is to expose instruction-level parallelism to the compiler and let the auto-vectorizer choose the instructions depending on its internal metrics. Given the increasing variety of HPC architectures competing for new exascale systems in the next years, the compiler remains the central tool for the programmer to achieve optimal performance on a given architecture.

We test the weak scaling of WOMBAT on the Cray XC40 "Hazel Hen" at HLRS Stuttgart. The machine features a Cray Aries interconnect, which uses a Dragonfly network with adaptive routers that are well suited for the high rate of small MPI messages in WOMBAT's communication pattern. Large pages of 16 MB are enabled by default. We note that the problem is well balanced, i.e., that the workload per node is the same on all ranks. Thus, WOMBAT's load-balancing capabilities are not tested here. For large problems, communication is only a small fraction of the total workload, which effectively hides communication inefficiencies. Hence, we choose a particularly small workload to clearly expose overhead from the MPI communication in the run. Unfortunately, this is not true for all weak scaling tests in the literature, and direct comparison with other codes is not always straightforward. We also note that for large-enough machines (millions of cores), communication overhead will eventually always be exposed. Thus, the argument that enough work is available per time step to hide communication inefficiencies is just another way of saying that an implementation does not scale beyond a certain point.

Here we use 3 × 4 × 4 patches per MPI rank with 18³ zones per rank, which results in a step time of about 2.5 s and a throughput of about 500,000 zones per second per node. We evolve the 3D problem for 100 steps, which means a total runtime of about 4 minutes. Every run was a separate submission on a different set of nodes on the system. The resulting mean step time over the number of nodes and cores is shown in Figure 36, from 1 node (24 cores) to 4096 nodes (98,304 cores). The minimum and maximum time per step is shown by the error bars. Ideal scaling at 2.45 s per update is shown by the dotted line. We also run the same test with the TVD+CTU solver, but with 32³ patches, which is its optimal cache-blocked patch size on Broadwell. The difference stems from the flattened array implementation, which TVD does not share with WENO5. The resulting scaling is shown by the blue line. The step time is reduced by a factor of 1.2 relative to the WENO solver, and the problem is a factor of 5.6 larger, so the TVD+CTU solver is a factor of 6.7 faster at scale than the WENO solver. But the WENO solver doubles the effective resolution, which would increase the runtime of TVD+CTU by a factor of 16. TVD+CTU runs with CFL = 0.4, half the step size of WENO. Thus, the WENO solver is more efficient ("faster") than TVD+CTU by a factor of about 2 × 2.3 = 4.6 at scale.

Zoom In Zoom Out Reset image size

We expect this difference in throughput (zones per second per code) to hold across architectures. Although TVD does not use flattened arrays, both implementations reach similar performance on the Broadwell architecture (14% versus 20% of peak), if properly cache-blocked. WOMBAT is written in vanilla Fortran 2008, which makes the code as portable as possible and gives access to the most mature compilers in the industry. Thus, the higher efficiency of WENO5 cannot be explained simply by implementation details or architecture and will, due its complete vectorization, very likely translate to accelerators as well. The cache-blocking properties will change of course, but the program design allows adapting to this side constraint and achieving optimal performance independent of cache size and properties.

We note that the scaling is quite flat at about 2.45 s per update with about 10% scatter around the mean. Differences arise from network contention on the system, as runtime is dominated by the slowest node during the run. They are much smaller than seen previously on the Blue Waters system with the Gemini interconnect (Mendygral et al. 2017). This view is supported by the decrease in scatter in step times for the largest runs, where a significant fraction of the machine is used. The largest run uses more than half of the Hazel Hen system and performs at about 150 TFlops, or about 5% of the peak performance of the system. We leave it to the reader to compare with other community codes (e.g., Nordlund et al. 2018).

Strong scaling refers to the speed-up of an algorithm as the problem size is kept fixed and the compute power is increased. Ideally, execution time should half as the computing power doubles. According to Amdahl's law (Amdahl 1967), the execution time of any parallel algorithm will eventually be dominated by the nonparallelizable part of the work, which can be in the form of branching, communication, or load imbalance.

We run a strong scaling test on the Cray XC40 "Hazel Hen" at HLRS Stuttgart 4 MPI ranks and 12 OpenMP threads per node. We use 18³ patches with an initial patch distribution of 96 × 12 × 12 patches, so the world grid has $1728\times {216}^{2}$ zones. For the final point at 512 nodes, we reduced the patch size to 9 × 18² zones.

The result is shown in Figure 37 as seconds per time step over the number of nodes and cores. WENO–Wombat follows the ideal scaling closely down to 128 nodes. The point at 512 nodes deviates from the ideal scaling, due to the decreased patch size. There was simply not enough work available anymore with 512 nodes, which limits WOMBAT's ability to further strong scale.

Zoom In Zoom Out Reset image size

WENO–WOMBAT already includes a treatment of cold supersonic flows using an elegant entropy scheme based on flux splitting (H. Jang et al. 2019, in preparation). A few improvements to the implementation presented here come to mind. A high-order CT scheme that matches the accuracy of the spatial interpolation would reduce magnetic field dispersion, albeit at considerable computational cost. High-order CT schemes exist for FV approaches (e.g., Londrillo & del Zanna 2004; Felker & Stone 2018; Verma et al. 2019), but need to be adapted to our FD algorithm. Global smoothness indicators could be used to improve the robustness of the scheme (e.g., Balbás & Tadmor 2006). Worth another look is likely a low-storage fourth-order Runge–Kutta integrator. A wide variety of strong-stability-preserving schemes exist, usually with five of more stages and CFL numbers >1 (Spiteri & Ruuth 2002). Furthermore, a ninth-order WENO implementation might become affordable once supercomputers reach exaflop performance, further doubling the effective resolution of the simulation and reducing storage and memory requirements. An extension of our WENO5 implementation would be straightforward. On the technical side, the flattened array implementation will make it trivial to port the WENO solver to accelerators. The OpenMP 4 standard would be the natural choice, providing a portable approach for a wide variety of accelerator technologies.

It has been shown that computational fluid dynamics without code-level transparency is not reproducible, even by the authors of the original publication: Mesnard & Barba (2017) found that open source code, logs, and data are a bare minimum to reproduce results years later. WOMBAT follows their reproducibility policy: sources and logs used to produce this document are available at https://wombatcode.org/publications/. Most data are made available, but have to adhere to certain space limitations. This data file also contains the complete source code tree of the WOMBAT version used in this paper, for reference purposes only. We did not yet include WENO5 in the public version of WOMBAT as we do not yet consider the scheme robust enough for applications.