Multidimensional upwind hydrodynamics on unstructured meshes using Graphics Processing Units I. Two-dimensional uniform meshes

TL;DR

This paper introduces a multidimensional hydrodynamics method on unstructured meshes optimized for GPUs, outperforming traditional schemes in accuracy and recognizing complex stationary states, with significant speedups.

Contribution

It presents a novel multidimensional Roe solver on unstructured meshes tailored for GPUs, enhancing accuracy and efficiency over traditional methods.

Findings

Significantly outperforms traditional Roe solver in test problems

Achieves 100-250x speedup on GPUs compared to desktop CPU

Recognizes multidimensional stationary states beyond hydrostatic solutions

Abstract

We present a new method for numerical hydrodynamics which uses a multidimensional generalisation of the Roe solver and operates on an unstructured triangular mesh. The main advantage over traditional methods based on Riemann solvers, which commonly use one-dimensional flux estimates as building blocks for a multidimensional integration, is its inherently multidimensional nature, and as a consequence its ability to recognise multidimensional stationary states that are not hydrostatic. A second novelty is the focus on Graphics Processing Units (GPUs). By tailoring the algorithms specifically to GPUs we are able to get speedups of 100-250 compared to a desktop machine. We compare the multidimensional upwind scheme to a traditional, dimensionally split implementation of the Roe solver on several test problems, and we find that the new method significantly outperforms the Roe solver in almost all cases. This comes with increased computational costs per time step, which makes the new method approximately a factor of 2 slower than a dimensionally split scheme acting on a structured grid.

Figures

Tables

Table 1. Table 1: Speed comparison of various codes on a hydrostatic test problem (constant density, constant pressure, zero velocity) on the unit square at a resolution equivalent to 128 × 128 128 128 128\times 128 . First column lists the name of the code, the second column lists relevant parameters, and the third column gives the computing time spent per cell per time step. Only a single CPU was used in all cases. Tests were performed on an Intel 2.2 2.2 2.2 GHz Core i7 with 8 GB DDR3 memory.

Code	Parameters	Time/cell/step
		($\mu s$)
astrix	$N$1	0.888
astrix	$B$2	2.09
rodeo	sb1	0.938
pluto	RK2, Roe, linear	0.864
gizmo	MFM	9.42
gizmo	MFV	10.2
rich	HLLC, Eulerian	19.6

Equations

$$x_i$$

$$y_j$$

$$r=\frac{l_{a}l_{b}l_{c}}{2\left|\begin{array}[]{cc}a_{x}-c_{x}&a_{y}-c_{y}\\ b_{x}-c_{x}&b_{y}-c_{y}\end{array}\right|},$$

$$\textsc{Orient2D}(a,b,c)=\left|\begin{array}[]{cc}a_{x}-c_{x}&a_{y}-c_{y}\\ b_{x}-c_{x}&b_{y}-c_{y}\end{array}\right|.$$

$$\left({\mathbf{r}}_a-{\mathbf{r}}_v\right)\cdot \left({\mathbf{r}}_b-{\mathbf{r}}_v\right)<0$$

$$\text{textsc}InCircle2D(a,b,c,d)=$$

$$\displaystyle\left|\begin{array}[]{ccc}a_{x}-d_{x}&a_{y}-d_{y}&(a_{x}-d_{x})^{2}+(a_{y}-d_{y})^{2}\\ b_{x}-d_{x}&b_{y}-d_{y}&(b_{x}-d_{x})^{2}+(b_{y}-d_{y})^{2}\\ c_{x}-d_{x}&c_{y}-d_{y}&(c_{x}-d_{x})^{2}+(c_{y}-d_{y})^{2}\end{array}\right|.$$

$$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x}(\rho u)+\frac{\partial }{\partial y}(\rho v)$$

$$\frac{\partial }{\partial t}(\rho u)+\frac{\partial }{\partial x}(\rho u^2+p)+\frac{\partial }{\partial y}(\rho uv)$$

$$\frac{\partial }{\partial t}(\rho v)+\frac{\partial }{\partial x}(\rho uv)+\frac{\partial }{\partial y}(\rho v^2+p)$$

$$\frac{\partial e}{\partial t}+\frac{\partial }{\partial x}(\rho hu)+\frac{\partial }{\partial y}(\rho hv)$$

$$e=\frac{1}{2}\rho \left(u^2+v^2\right)+\frac{p}{\gamma -1},$$

$$\frac{\partial \mathbf{W}}{\partial t}+\nabla \cdot \mathcal{F}=0,$$

$$\varphi ={\oint }_{\partial V}\mathcal{F}\cdot \mathbf{n}dS,$$

$$\frac{\partial \mathbf{W}}{\partial t}+\frac{\partial \mathbf{F}}{\partial x}=0,$$

$$\frac{\partial \mathbf{W}}{\partial t}+\mathcal{A}(\mathbf{W})\frac{\partial \mathbf{W}}{\partial x}=0,$$

$${\mathbf{W}}_i-{\mathbf{W}}_{i-1}=\sum _{p=1}^q{\alpha }_{i-1/2,p}{\mathbf{r}}_{i-1/2,p},$$

$${\mathbf{\alpha }}_{i-1/2}={\mathcal{R}}_{i-1/2}^{-1}\left({\mathbf{W}}_i-{\mathbf{W}}_{i-1}\right),$$

$${\mathbf{W}}_i(t+\Delta t)=\frac{1}{\Delta x}{\int }_{x_{i-1/2}}^{x_i+1/2}{\mathbf{W}}^{\ast }(x,t+\Delta t)dx.$$

$${\mathbf{W}}_i(t+\Delta t)={\mathbf{W}}_i(t)-\frac{\Delta t}{\Delta x}\sum _{p=1}^q{\left({\lambda }^p\right)}^{+}{\alpha }^p{\mathbf{r}}^p,$$

$$\mathbf{W}=(\rho ,\rho u,e{)}^T,$$

$$\mathbf{Z}=\sqrt{\rho }(1,u,h{)}^T,$$

$$\mathcal{A}(\bar{\mathbf{Z}})({\mathbf{W}}_i-{\mathbf{W}}_{i-1})={\mathbf{F}}_i-{\mathbf{F}}_{i-1};$$

$$\mathbf{W}(t+\Delta t,x)$$

$$\frac{d\mathbf{W}}{dt}=\mathcal{L}(\mathbf{W}),$$

$$\frac{\partial \mathbf{W}}{\partial t}+\frac{\partial \mathbf{F}}{\partial x}+\frac{\partial \mathbf{G}}{\partial y}=0,$$

$$\frac{\partial \mathbf{W}}{\partial t}+\mathcal{A}\frac{\partial \mathbf{W}}{\partial x}+\mathcal{B}\frac{\partial \mathbf{W}}{\partial y}=0,$$

$$\mathbf{W}(t+\Delta t,x)$$

$$\varphi ={\mathbf{F}}_i-{\mathbf{F}}_{i-1}=\bar{\mathcal{A}}\left({\mathbf{W}}_i-{\mathbf{W}}_{i-1}\right).$$

$$\sum _{p=1}^q{\left({\lambda }^p\right)}^{+}{\alpha }^p{\mathbf{r}}^p$$

$$\sum _{p=1}^q{\left({\lambda }^p\right)}^{-}{\alpha }^p{\mathbf{r}}^p={\bar{\mathcal{A}}}^{-}\left({\mathbf{W}}_i-{\mathbf{W}}_{i-1}\right).$$

$${\varphi }_i$$

$${\varphi }_{i-1}$$

Full text

††pagerange: Multidimensional upwind hydrodynamics on unstructured meshes using Graphics Processing Units I. Two-dimensional uniform meshes–173††pubyear: 2016

1 Numerical gas dynamics

The observation that $99\%$ of the visible matter in the Milky Way is in gaseous form (Draine et al., 2007) makes gas dynamics an important part of the study of various systems in astrophysics, from stars and supernovae to accretion discs and gas giant planets. The non-linear nature of the governing equations makes numerical simulations an essential tool to make progress in our understanding of these systems. A wide variety of methods exists for solving the equations of gas dynamics, all of which perform better at some problems than others. Below, we give a very brief overview of the numerical landscape, which serves to put our new method in context.

Numerical methods for gas dynamics solve discretised versions of the governing equations. A first choice when selecting a method to tackle a particular problem in hydrodynamics is whether the resolution elements move with the gas flow (the Lagrangian approach) or are fixed in space (the Eulerian approach). Eulerian methods use a computational mesh to discretise space into small elements (usually squares in two dimensions, cubes in three dimensions). Popular methods include methods based on finite difference approximations such as Pencil (Brandenburg & Dobler, 2002) and zeus (Stone & Norman, 1992) or methods based on Riemann solvers such as flash (Fryxell et al., 2000) and athena (Stone et al., 2008). These two classes of Eulerian methods take in some sense opposite viewpoints of the underlying solution. Finite difference methods assume the flow to be smooth, and in order to prevent unphysical oscillations near discontinuities add artificial viscosity in order to smooth these out. Godunov methods based on Riemann solvers view the underlying solution as a set of discontinuities, for which the time evolution can be computed by solving Riemann problems. While this means discontinuities in the flow can be handled in a natural way, it also restricts the method to be first order accurate in space and time. In regions of smooth flow, higher-order methods can be used safely. In order to avoid unphysical oscillations near discontinuities, the contribution of the high-order method should be limited. If the chosen limiter function is Total Variation Diminishing (TVD, van Leer, 1974) oscillations near discontinuities can be avoided.

In addition to these two classes of solvers, spectral methods are available (e.g Lesur & Longaretti, 2005; Burns et al., 2016; Lecoanet et al., 2014), that solve the Navier-Stokes equations in spectral (usually Fourier) space. Advantages of Eulerian methods in general are the low intrinsic numerical dissipation, and, in the case of Riemann solvers (Falle, 2002), the automatic addition of the correct amount of dissipation for shock waves. Main disadvantages include that dissipation is highly non-linear, making it more difficult to control , and for example velocity-dependent (Springel, 2010, hereafter S10), and that it is not trivial to vary the resolution within the computational domain while maintaining low dissipation. For example, adaptive mesh refinement (AMR, Berger & Oliger, 1984) leads to locations in the mesh where dissipation is especially high (where jumps in resolution occur and the update is usually only correct up to first order, see however Schaal et al. (2015)). Slowly varying spatial resolution can be achieved by choosing an appropriate coordinate transformation, for which second-order methods exist (e.g. Eulderink & Mellema, 1995). However, in that case it has to be known in advance where high resolution is needed, and, furthermore, if orthogonal coordinates are desired, this limits the ability to achieve high resolution locally.

While staggered mesh Lagrange plus remap methods do exist, in which the flow is remapped onto the mesh every time step (e.g. Woodward & Colella, 1984; Pember & Anderson, 2000), the most well-known Lagrangian method in astrophysics is the mesh-free method of Smoothed Particle Hydrodynamics (SPH, Lucy, 1977; Gingold & Monaghan, 1977), where the gas is represented by a set of particles that move with the flow. The Lagrangian nature of SPH gives it two main advantages over traditional grid-based methods: errors associated with large-scale bulk motion of the fluid are virtually non-existent (S10), and resolution automatically follows the concentration of mass. This makes SPH competitive in for example collapse problems in cosmology (e.g. Schaye et al., 2015) and star and planet formation (e.g. Bate et al., 2003; Mayer et al., 2002). Disadvantages of SPH include its relatively large numerical dissipation for certain problems, especially when low-density regions are dynamically important (e.g. de Val-Borro et al., 2006) or in shear flows (Agertz et al., 2007). A meshless method that can have high resolution in arbitrary locations was presented in Maron et al. (2012).

Recent studies have focused on bridging the gap between SPH and grid-based methods, in an effort to get the best of both worlds. Examples include arepo (S10), rich (Yalinewich et al., 2015) and gizmo (Hopkins, 2014, 2015). One can arrive at the class of moving mesh codes by starting with SPH, but seeing the particles as mesh generation points, and subsequently solving the Euler equations on this (necessarily unstructured) mesh. If the mesh points are fixed, an Eulerian method on an unstructured grid is obtained. If the points are allowed to move, a moving mesh code results. The use of these mesh-generating points and their Voronoi tessellation makes the mesh evolve in a continuous manner, and one can avoid mesh-tangling problems of traditional Arbitraly Lagrangian-Eulerian (ALE) methods (e.g. Vachal et al., 2004). While in principle one is free to move the mesh points with any velocity, if the mesh velocity is set equal to the local gas velocity, one obtains a Lagrangian method. A cylindrical moving mesh code was recently presented by Duffell (2016). One problem encountered with moving meshes is grid noise (Bauer & Springel, 2012; Hopkins, 2015), caused by changes in topology as the mesh evolves which lead to volume inconsistency errors (Yalinewich et al., 2015). Several fixes have been proposed, from smoothing the velocities of mesh points (Duffell & MacFadyen, 2015), to regularising the mesh (Mocz et al., 2015), to directly attacking the volume inconsistency (Steinberg et al., 2016).

Given an Eulerian method, for example based on a Riemann solver, acting on a structured grid, adapting it to work on an unstructured grid is a non-trivial undertaking. Indeed, most multi-dimensional Eulerian methods are built using one-dimensional solvers, combined in such a way, making use of the structured nature of the mesh (often logically rectangular), to yield an accurate approximation to the multidimensional problem (e.g. Strang, 1968; Colella, 1990; Balsara, 2012). Even more, Fourier spectral methods require a regular mesh from the very beginning. On the other hand, unstructured meshes offer desirable properties such as more isotropic numerical diffusion, and the ability to refine the mesh in arbitrary ways without the need for jumps in resolution.

Unstructured meshes usually come in the form of Delaunay triangulations (Delaunay, 1934) or Voronoi tesselations (Dirichlet, 1850; Voronoi, 1907), for reasons that we will go into in section 2. Triangular grids, structured or unstructured, have an additional advantage that they allow for a multidimensional analog of Roe’s approximate Riemann solver (Struijs, 1994). That is, it is possible to design a multidimensional solver without having to rely on one-dimensional building blocks. This clearly has advantages over traditional Eulerian methods, especially for flows that are not aligned with any axis of the grid. While these multidimensional upwind, or residual distribution, methods have been around for several decades, initially they were designed to solve for steady flows only. The time-dependent case turned out to be relatively tricky to work out with various different formulations (e.g. Ferrante & Deconinck, 1997; De Palma et al., 2005). Moreover, these were all implicit time integration schemes and therefore relatively expensive.

More recently, an explicit formulation was derived (Ricchiuto & Abgrall, 2010), making multidimensional upwind methods potentially competitive for time dependent astrophysical flows, which are often very compressible but also multidimensional. In this paper, we present and test a two-dimensional version of a multidimensional upwind method in an astrophysical context in the form of astrix111Freely available as an open-source project at https://github.com/SijmeJan/Astrix/ (AStrophysical fluid dynamics on unstructured TRiangular eXtreme grids). Another advantage of multidimensional upwind methods compared to the more common Riemann solvers is that they employ a very compact stencil: the vast majority of all calculations (in particular calculating the fluxes between cells) are done using data from one triangle only. In comparison, a dimensionally split scheme based on a Riemann solver needs four cells in each direction in order to compute a second order accurate interface flux (see e.g. LeVeque, 2002). The number of bytes that need to be read in order to compute a flux is therefore much larger in traditional Eulerian codes, which makes multidimensional upwind methods excellent candidates to port to Graphics Processing Units (GPUs).

The rest of this paper is structured as follows. In section 2 we introduce unstructured grids and how to generate them, while in section 3 we describe the residual distribution methods that are part of astrix. In section 4 we discuss the GPU implementation of both mesh generation and residual distribution . In section 5 these methods are tested on one and two dimensional problems and we give a discussion in section 6. We conclude in section 7.

2 Unstructured grids

2.1 Basic definitions

A grid, or mesh, is defined by a set of nodes, or vertices, together with a recipe for getting from one vertex to its neighbours. A mesh can be said to be structured if the location of the vertices and their interconnections follow a simple pattern, which usually leads to a high degree of symmetry. In two dimensions, an $m\times n$ structured Cartesian mesh on the unit square can be defined as a collection of vertices $(i,j)$ with coordinates

[TABLE]

where $0\leq i<m$ and $0\leq j<n$, together with the connectivity rules that vertex $(i,j)$ is connected to $(i-1,j)$ and $(i+1,j)$ in the $x$-direction, and to $(i,j-1)$ and $(i,j+1)$ in the $y$-direction. Such a mesh consists of square cells. The same collection of vertices with different connectivity rules could for example lead to triangular cells (for an example see Fig. 1).

On the other hand, for an unstructured mesh no simple rules exist for the location of the vertices. An example is shown in Fig. 2. Obviously this makes the implementation more complicated and more memory intensive, since the location of all vertices has to be stored explicitly, rather than using simple rules to work out the coordinates. However, unstructured grids provide absolute freedom on where to place the vertices, which has three main advantages:

• •

It is possible to handle complex geometries. When studying the flow across an aircraft wing, for example, it is necessary to get as close to the real shape of the wing as possible . Embedding a shape that is not rectangular in a regular Cartesian mesh is hopeless.

• •

It is possible to construct meshes with no preferred directions , which leads to more uniform numerical dissipation. For a 2D structured Cartesian mesh, numerical dissipation will strongly depend on the direction of the flow with respect to the coordinate axes. This leads to the famous carbuncle instability (Peery & Imlay, 1988; Quirk, 1994).

• •

One has much more freedom to vary the resolution of the mesh from place to place in a smooth way. For structured meshes the options are limited. For a structured cylindrical mesh for example, it is possible to increase the resolution towards small radii by choosing a logarithmic radial coordinate. However, this affects all cells in the inner parts, while the region where high resolution is required may be very limited in azimuthal extent. Of course, for structured meshes there exists the powerful technique of AMR to obtain high resolution locally, but this technique leads to boundaries between regions of coarse and fine resolution (Berger & Oliger, 1984) where additional interpolation errors occur, which in some cases may be unacceptable.

While it is possible to generate quadrilateral unstructured meshes, the more common cell choice is the triangle. For a given set of points , there exist many ways of interconnecting them using triangles. This means we can choose a triangulation that is in some sense optimal. A useful criterion of the quality of the grid is the minimum opening angle of any triangle. Meshes with small angles often lead to numerical problems for simulations, since numerical diffusion will be very non-isotropic (e.g. Babuska & Aziz, 1976) . This former problem happens in Cartesian structured meshes when the cells have a very large aspect ratio ; for example $m\ll n$ in equation (1). For Cartesian structured meshes, we would usually like the cells to be as square as possible, i.e. $m=n$ in equation (1). For a triangular mesh this translates into having triangles that are as close to equilateral as possible. While it is not possible to achieve this limit in practice, for example because of boundary constraints, we would still like to maximise the minimum angle for any triangle in the grid. For a given set of vertices, this fixes the triangulation, since it is the Delaunay triangulation (Delaunay, 1934) that achieves this.

2.2 Delaunay triangulation

2.2.1 Circumcircles

The circumcircle of a set of two or three vertices is a circle that passes through all vertices in the set. For a set of three vertices, this circle is unique. For a set of two vertices, there are infinitely many circumcircles. A triangle $T$ is called Delaunay if its circumcircle $C$ is empty, i.e. there are no mesh-generating points inside $C$ . An edge $E$, consisting of vertices $u$ and $v$, is called Delaunay if there exists a circumcircle of $u$ and $v$ that is empty. Note that if a triangle is Delaunay, all of its three edges are automatically Delaunay, since there exists a circumcircle $C$ of any of the edges of $T$ that is empty (take $C$ to be the circumcircle of $T$). The reverse is also true, but less easy to show formally.

A triangulation of a two-dimensional space where all triangles and all edges are Delaunay is guaranteed to exist, and in addition it is unique if no four vertices lie on the same circle (and, more trivially, not all vertices are collinear). If this condition is violated, the triangulation is no longer unique: edges that are Delaunay can be crossing, and in order to obtain a valid triangulation a selection of edges has to be made. In practice, this happens automatically during mesh construction (see below), but it is a first indication that numerical roundoff errors will play a prominent role in mesh construction. If we bring four vertices closer and closer to being on the same circle, at some point the triangulation will become degenerate because of roundoff errors. This situation has to be handled with care, since if some parts of the algorithm detect a degeneracy while other parts do not, which can easily happen when working close to roundoff limits, the whole algorithm will break down. There is therefore a need for exact geometric predicates (see section 4.1.3). In Fig. 3 a Delaunay triangulation is shown together with each triangle’s circumcircle.

2.2.2 Edge flipping

Consider an edge $E$ in a triangulation (not necessarily Delaunay), together with only its two neighbouring triangles (see Fig. 4). The four vertices define a containing quadrilateral. Define $E$ to be locally Delaunay if empty circumcircle of $E$ exists. Obviously, if $E$ is not locally Delaunay it is not Delaunay. The converse is not true: $E$ may be locally Delaunay but not Delaunay. However, if all edges in the triangulation are locally Delaunay this means that all edges are Delaunay.

Define $E^{\prime}$ as the edge that would exist if we connected the two vertices of the containing quadrilateral not part of $E$ ($a$ and $d$ in Fig. 4). If $E$ is not locally Delaunay, it follows that $E^{\prime}$ would be locally Delaunay. This means that we can remove an edge that is not locally Delaunay by flipping $E$ (i.e. replacing it with $E^{\prime}$). Not all edges are flippable, for example if the containing quadrilateral is not convex. However, it can be shown that all edges that are not locally Delaunay can be flipped.

It follows that if a triangulation is not Delaunay, at least one edge is not locally Delaunay and can be flipped. It is intuitively clear, but slightly more difficult to show, that each flip moves the triangulation closer to the Delaunay triangulation222This is only true in 2D; in 3D, convergence is not guaranteed.. When there are no more edges to flip, the resulting triangulation is the Delaunay triangulation.

Note that some edges can not be flipped because it would make the triangulation invalid. In Fig. 3, this is true for the almost horizontal edge near the top. The reason is that the containing quadrilateral is not convex. However, it can be shown that if an edge can not be flipped, it has to be Delaunay.

2.2.3 Triangle quality

The main reason for the popularity of Delaunay triangulations for mesh generation is that the resulting triangles are of high quality in the sense that small angles can be avoided as long as the boundary of the domain does not require them . In fact, among all possible triangulations of a set of vertices, the Delaunay triangulation maximises the minimum angle present in the mesh333Unfortunately, this result only holds in two dimensions (Lawson, 1977). This can be shown in a straightforward way by first noting that flipping an edge to make it locally Delaunay always increases the minimum angle present in the containing quadrilateral. Now any valid triangulation can be transformed into a Delaunay triangulation by a sequence of edge flips. Since any edge flip increases the minimum angle, this means that the Delaunay triangulation maximises the minimum angle.

2.2.4 Incremental insertion algorithm

The discussion on edge flipping suggests the following algorithm for finding the Delaunay triangulation of a set of vertices: start with any valid triangulation, and perform edge flips until all edges are locally Delaunay. The resulting triangulation is the Delaunay triangulation.

Finding an initial triangulation for a large set of input vertices is cumbersome. Therefore, in practice vertices are inserted one by one, starting from an initial triangle large enough to contain all subsequent vertices (Lawson, 1977). Such algorithms are usually called incremental insertion algorithms. When implemented in its simplest form, incremental insertion can be slow, because in principle the addition of a single vertex may lead to the flipping of all edges in the mesh. However, if the vertices are inserted in random order (Guibas et al., 1992), incremental insertion becomes a competitive method (Su & Drysdale, 1997). Other algorithms include divide and conquer (Guibas & Stolfi, 1985; Dwyer, 1987), sweepline (Fortune, 1986), gift-wrapping (Dwyer, 1991) and algorithms based on the convex hull (Barber et al., 1996). The reason for choosing incremental insertion over any of the other methods is that by its very nature is very well suited to deal with cases where we do not know in advance the position of all vertices in the mesh. This is the situation we are in when generating unstructured meshes.

2.3 Delaunay refinement

Even though the Delaunay triangulation is optimal in the sense that it maximises the minimum angle, low-quality triangles can still be seen in Fig. 3. This is because the location of the vertices were chosen randomly. In practice, we are free to choose the majority of the positions of the vertices. Possible constraints on the location of the vertices include any fixed boundaries (a wall, an aircraft wing, a planet) and any resolution constraints on the local density of vertices. Within these constraints, there is still a lot of freedom in choosing vertex locations. This can be done in an optimal way as to guarantee a mesh of a certain quality, that is, a minimum opening angle that is larger than a certain value. This way of choosing vertices, often called Delaunay refinement, is discussed next. For a more detailed discussion, see Shewchuk (2002).

2.3.1 Removing low-quality triangles

First of all, the idea of a triangle of low quality can be made quantitative. Define $\beta$ to be the ratio of the circumradius to the shortest edge of a triangle. A straightforward calculation shows that $2\beta=1/\sin\alpha$, where $\alpha$ is the smallest angle of the triangle. Therefore, in order to avoid small angles, we need to avoid triangles with large values of $\beta$, say we require $\beta\leq B$ for some appropriate bound $B$.

A low-quality triangle $T$ can be removed by adding a new vertex $v$ to the mesh, exactly at the circumcentre of $T$. It is clear that $T$ can not be part of the new mesh, since its circumcircle is no longer empty. However, it is also clear that any new triangle created by inserting $v$, will have a shortest edge that is at least the radius of the circumcircle of $T$. Since $\beta>B$ for $T$, we have that the shortest edge of any new triangle is at least $B$ times the shortest edge of $T$. For $B\geq 1$, this means that new edges will have at least the length of the minimum edge length in the initial triangulation $|E|_{\mathrm{min}}$. Therefore, an algorithm that removes triangles with $\beta>B$ must terminate eventually if $B\geq 1$, since at some point the vertex density becomes so high that no new triangles can be created with minimum edge length larger than $|E|_{\mathrm{min}}$. The three most well-known Delaunay refinement algorithms use $B=\sqrt{2}$ (Ruppert, 1995) and $B=1$ (Chew, 1989, 1993).

In addition to the quality constraint $\beta\leq B$, we can impose a size constraint, by removing triangles that are too large in the same way as removing low-quality triangles. The maximum size $|T|_{\mathrm{max}}$ can be a function of space, allowing for non-uniform meshes. As long as $|T|_{\mathrm{max}}({\bf x})>M$ for some positive constant $M$, the algorithm is still guaranteed to terminate.

2.3.2 Splitting segments

The main distinction between the algorithms of Chew (1989), Chew (1993) and Ruppert (1995) lies in the way mesh boundaries are treated. These boundaries are part of the input of the Delaunay refinement algorithm: it takes as input the vertices that make up the boundary together with their interconnection. The simplest case would be if the computational domain is a rectangle, which can be specified with four vertices (the corners) and four edges (the sides). More complicated cases would include an outer boundary that is a circle, or a hole in the computational domain (an inner boundary) in the form of an aircraft wing. These edges are special in the sense that they will have to be part of the final mesh. Such edges will be referred to as segments444Segments do not have to be part of the boundary: in principle the term refers to any input edge that needs to be part of the final triangulation.

Ruppert (1995) provides an elegant way of dealing with segments. Define the diametric circle of a segment to be the smallest circle that encloses the segment. A segment is defined to be encroached if any vertex lies inside the diametric circle. During Delaunay refinement, if a new vertex would lead to encroached segments, the vertex is not inserted: instead the segments in question are split by inserting vertices at their midpoints. Note that a non-encroached segment is an edge that is Delaunay since there exists a circumcircle, namely the diametric circle, that is empty. This way, the algorithm ensures that boundary segments will be part of the final triangulation. An additional advantage is that it is also guaranteed that a new vertex will be inside the existing triangulation: if the circumcentre of a triangle happens to be outside the triangulation, this must mean that somewhere a segment is encroached.

2.4 Grid generation

We now go over the different steps in our Delaunay refinement algorithm. Every iteration consists of six steps, which are described below and are repeated until no more vertices need to be inserted. It is important to realise that the resulting mesh, defined by the input boundary vertices and the desired quality, is not unique. In particular, the locations of the vertices depend on the order in which they were inserted.

2.4.1 Finding low-quality triangles

First of all, for every triangle $T$ with vertices $(a,b,c)$ in the mesh we compute its circumradius $r$:

[TABLE]

where $l_{a}$, $l_{b}$ and $l_{c}$ are the length of the edges opposite vertices $a$, $b$ and $c$, respectively, and its circumradius-to-shortest-edge ratio $\beta$. If either $r>R$ (triangle too big) or $\beta>B$ (triangle too low quality), a vertex will be inserted at its circumcentre. Here, $R$ and $B$ are user input quantities specifying the desired triangle size ($R$) and the desired triangle quality ($B$). This step leads to a list of vertices to add into the mesh.

2.4.2 Finding triangles containing the new vertices

Next, for every new vertex we find the triangle containing this vertex (often this is not the original triangle, especially if this triangle has a large value of $\beta$). A triangle $T$ with vertices $a$, $b$ and $c$ contains vertex $v$ if $v$ lies on the ‘correct’ side of all three edges of $T$. If $(a,b,c)$ are in counterclockwise order, then $v$ lies in $T$ if both $(a,b,v)$, $(b,c,v)$ and $(c,a,v)$ are in counterclockwise order. Following Shewchuk (1997), we will refer to the counterclockwise test as Orient2D$(a,b,c)$, a function that returns a positive value if $(a,b,c)$ are in counterclockwise order, negative if they are in clockwise order, and zero if they are collinear. This function can be implemented as a matrix determinant555Note that the same determinant appears in the denominator in equation (2):

[TABLE]

It should be clear that in order to maintain a valid triangulation, it is extremely important that Orient2D$(a,b,c)$ gives the correct result, even in the presence of round-off errors. This is difficult when $(a,b,c)$ are close to collinear. If we try and place a vertex $v$ in triangle $T$ but extremely close to the edge that is shared by $T$ and $T^{\prime}$, due to round-off errors a naive implementation of Orient2D may decide that $v$ lies in both $T$ and $T^{\prime}$, or not in either of them. Such mutually contradictory results inevitably lead to nonsensical results and invalid triangulations. One might hope that for straightforward computational domains (no complicated boundaries) such cases never show up, but in practice they always do. Therefore, we evaluate Orient2D using exact geometric predicates (see section 4.1.3).

While the triangle containing $v$ may not be the original triangle $T_{0}$ leading to the insertion of $v$, often it is close to it. Therefore, we start searching in $T_{0}$, and move to a neighbouring triangle in the direction of $v$ until we have found a triangle containing $v$. In some cases, $v$ will be located exactly on an edge and those have to be dealt with separately.

2.4.3 Testing if any new vertex encroaches upon a segment

Vertices can not be inserted if they lead to an encroached segment. A vertex $v$ with coordinates ${\bf r}_{v}$ encroaches upon a segment with vertices $a$ and $b$ with coordinates ${\bf r}_{a}$ and ${\bf r}_{b}$ if

[TABLE]

There is usually no need for exact predicates in this computation.

If vertex $v$ is to be inserted in triangle $T$ with vertices $(a,b,c)$ we check all triangles that have either $a$, $b$ or $c$ as a vertex. Should $v$ be inserted onto an existing edge $E$, we check all triangles that share a vertex with the two triangles sharing $E$ as an edge. If any of these triangles has a segment for an edge upon which $v$ encroaches, $v$ is not inserted at its original location but on the middle of this particular segment, thereby splitting the segment. Of course, inserting $v$ might have lead to multiple encroached segments, but this approach is computationally advantageous as the number of vertices to be inserted does not change. Any difficulties arising from this approach are dealt with in section 2.4.5.

Note that by checking all triangles that have either $a$, $b$ or $c$ as a vertex, we make sure that $v$ does not encroach upon any segment. For if $v$ encroaches upon any other segment, it automatically follows that either $a$, $b$ or $c$ encroaches upon this segment as well, which can not be the case if the original mesh had no encroached segments.

2.4.4 Insert new vertices into mesh

Next, we actually insert the vertices into the mesh. If vertex $v$ is to be inserted in triangle $T$, $T$ is split into three triangles that all have $v$ as a vertex. This operation therefore adds two triangles to the mesh. If $v$ is inserted on an edge $E$, the two triangles sharing $E$ are each split into two, thereby again adding two triangles to the mesh. When splitting a segment, the single triangle of which the segment is part is split into two, thereby adding one triangle to the mesh. Note that at this point, we do not maintain a Delaunay triangulation.

2.4.5 Split any encroached segments

Occasionally, splitting a segment leads to another encroached segment. In order to deal with this, we check any vertex that was inserted on a segment if it encroaches any other segment. If this is the case, the segment in question is split by adding an extra vertex to the mesh. In practice this only very rarely happens.

2.4.6 Maintain Delaunay triangulation

With the insertion of new vertices and therefore the addition of new triangles the current triangulation will usually no longer be Delaunay. The final step in the iteration is to transform the current triangulation into a Delaunay triangulation, which is achieved by edge-flipping.

First of all, all relevant edges in the mesh are checked for Delaunay-hood. Consider edge $E$ and its containing quadrilateral, which consists of triangles $T$ and $T^{\prime}$ (see Fig. 4). Let the vertices of $T$ be $(a,b,c)$, and let $d$ be the one vertex of $T^{\prime}$ that is not part of $T$. Then $E$ is Delaunay if $d$ does not lie inside the circumcircle of $T$. Following Shewchuk (1997), define the InCircle2D$(a,b,c,d)$ test to return a positive value when $d$ lies inside the circle defined by $(a,b,c)$ (assuming $(a,b,c)$ are in counterclockwise order), a negative value when $d$ lies outside this circle, and zero if $d$ lies exactly on this circle. As with Orient2D, InCircle2D can be expressed as a determinant:

[TABLE]

Again, we need the InCircle2D test to give the exact result, even in the case of finite precision arithmetic. Otherwise, not only could the flip algorithm get stuck, continuously flipping the same edge, but it could also flip an edge that is not flippable because the containing quadrilateral is concave, yielding an invalid triangulation and therefore breaking the Delaunay refinement algorithm.

Flipping an edge may create new edges that are not Delaunay. Therefore, the process of checking and flipping is iterated until all edges are Delaunay. When the triangulation is Delaunay again, the next iteration can be started by testing all triangles if they match the quality constraints. When there are no more vertices to add, the algorithm exits.

3 Residual distribution schemes

We now turn to the hydrodynamic solver acting on our unstructured grid. The equations to be solved are conservation of mass, momentum and energy in two spatial dimensions:

[TABLE]

where $\rho$ is the density, ${\bf v}=(u,v)^{T}$ is the velocity vector, $p$ is the pressure and $e$ is the total energy for an ideal gas equation of state :

[TABLE]

in which the last term denotes the internal energy, which is specified by the pressure under the assumption of a perfect gas, with ratio of specific heats $\gamma$. Finally, $h=(e+p)/\rho$ is the fluid enthalpy.

Above conservation laws can be written concisely as:

[TABLE]

where ${\bf W}$ is the state vector and $\mathcal{F}=({\bf F},{\bf G})$ is the flux term, containing the flux in the $x$ direction ${\bf F}$ and the flux in the $y$ direction ${\bf G}$. Integrating over a volume $V$ and using Gauss’ theorem shows that the time evolution of ${\bf W}$ is governed by the residual $\phi$:

[TABLE]

where ${\bf n}$ is the outward-pointing unit normal vector of surface element $dS$. If $\phi=0$, the state is stationary. A residual distribution scheme, as the name suggests, is a numerical scheme that solves a discrete version of (14) by taking $\phi$ and distributing it over neighbouring cells.

3.1 Roe solver

The Roe solver is a well-known approximate Riemann solver that is part of many astrophysical fluid dynamics packages (e.g. Mignone et al., 2007; Stone et al., 2008). However, it can also be viewed as a residual distribution scheme in one spatial dimension, and thereby making a connection between more traditional methods and the framework presented in the next sections.

3.1.1 Brief derivation

Consider a system of $q$ hyperbolic conservation laws in one spatial dimension:

[TABLE]

where ${\bf W}$ is the state vector and ${\bf F}$ the flux vector. Write the conservation laws in quasi-linear form:

[TABLE]

where $\mathcal{A}=\partial{\bf F}/\partial{\bf W}$ is the Jacobian.

Now consider a uniform grid with cell centres $x_{i}$ and spacing $\Delta x$. Consider two neighbouring grid cells $i-1$ and $i$, and corresponding states ${\bf W}_{i-1}$ and ${\bf W}_{i}$. The interaction of these two grid cells can be seen to arise from a state jump at the cell interface. The setup of two constant states separated by a discontinuity is known as a Riemann problem and has an analytic solution (e.g. Lax, 1957). Every cell interface has its own Riemann problem defined by the neighbouring states, and for small enough time steps these Riemann problems will be independent. Solving the Riemann problems then yield for example interface fluxes that can then be used to update the state in the cells (see e.g. LeVeque, 2002).

For a linear system, i.e. $\mathcal{A}$ does not depend on ${\bf W}$, the solution to the Riemann problem consists of a set of discontinuities travelling at speeds given by the eigenvalues of $\mathcal{A}$ (see Fig. 5). The strength of each discontinuity can be found by projecting the initial state jump onto the right eigenvectors of $\mathcal{A}$:

[TABLE]

where ${\bf r}_{i-1/2,p}$ is the $p$th eigenvector of $\mathcal{A}$ and $\alpha_{i-1/2,p}$ is the corresponding projection coefficient, which are found by

[TABLE]

where $\mathcal{R}$ is the matrix containing the right eigenvectors of $\mathcal{A}$ as columns.

This completes the solution ${\bf W}^{*}(t,x)$, which we can use to do a time step from $t$ to $t+\Delta t$ by averaging the solution over a grid cell:

[TABLE]

Since ${\bf W}^{*}$ is piecewise constant, the integration is easily done, yielding

[TABLE]

where $\left(\lambda^{p}\right)^{+}$ stands for $\mathrm{max}(0,\lambda^{p})$, with $\lambda^{p}$ the $p$th eigenvalue of $\mathcal{A}$. Note that we have only considered the cell interface between $i$ and $i-1$: there will be a similar contribution to the update of ${\bf W}_{i}$ from the interface with cell $i+1$.

Of course, the governing equations of gas dynamics are non-linear. The Roe solver (Roe, 1981) is defined by a suitable linearisation of $\mathcal{A}({\bf W})$. For the one-dimensional Euler equations, with

[TABLE]

where $\rho$ is the density, $u$ the velocity and $e$ the total energy, Roe (1981) found that by using a parameter vector

[TABLE]

where $h=(e+p)/\rho$ is the fluid enthalpy, the matrix $\mathcal{A}$ evaluated at $\bar{\bf Z}=({\bf Z}_{i}+{\bf Z}_{i-1})/2$ provides a linearisation with desirable properties. In particular

[TABLE]

a property necessary for a conservative scheme.

3.1.2 Second-order accuracy

The scheme presented above is only first-order accurate, i.e. the largest error term is proportional to $\Delta x$. It is possible to increase the order of accuracy by considering higher-order terms in the Taylor expansion of the solution:

[TABLE]

The term proportional to $\Delta t$ is dealt with by the first order scheme above, while the term proportional to $\Delta t^{2}$ must come from considering a linear reconstruction of the solution (e.g. LeVeque, 2002). The contribution of this last term involves information from cells further away from the interface under consideration (i.e. $i-2$ or $i+1$) and has to be limited in order to avoid spurious oscillations near shocks. TVD limiter functions come in many flavours, from the least compressive minmod limiter to the very compressive superbee limiter (e.g. Sweby, 1984). In regions of smooth flow, the update is done using the Lax-Wendroff scheme (Lax & Wendroff, 1960) , and is second-order in both space and time.

A different approach to improve the order of accuracy is not based on (25), but in stead separates discretization in space and time. Dealing with space first leads to an ordinary differential equation (ODE)

[TABLE]

where $\mathcal{L}$ is the operator governing the spatial discretization. Above equation can now be solved by a second-order ODE solver, leading to a second-order accurate method if the spatial discretization is second order as well. Again, a limiter has to be applied to avoid oscillations near discontinuities, and care must be taken that the ODE integrator does not introduce oscillations. This method of lines has the advantage that, unlike the Taylor series approach, it is straightforward to extend the method to higher than second order. The main disadvantage is that it needs to solve more than one Riemann problem per cell interface per time step, again unlike the Taylor series approach.

3.1.3 Approaches for more than one spatial dimension

The one-dimensional Roe solver, or any other Riemann solver, can be used to build a multidimensional numerical method. Here we highlight some of the problems that arise, since they are pertinent to our discussion later.

Consider the two-dimensional hyperbolic system

[TABLE]

with quasi-linear form

[TABLE]

with $\mathcal{B}=\partial{\bf G}/\partial{\bf W}$. In this section, we can afford to deal only with the linear problem and therefore take $\mathcal{A}$ and $\mathcal{B}$ to be constant. A Taylor expansion of the solution is given by

[TABLE]

The first-order Roe solver takes care of the terms proportional to $\Delta t$. If only a first order method is required, all interactions between grid cells can be taken into account simultaneously through applying (21) four times for every cell (one for every neighbour).

Terms proportional to $\Delta t^{2}$ come in two flavours: those containing $\mathcal{A}^{2}$ and $\mathcal{B}^{2}$, which in the one-dimensional case are taken care of by correction fluxes or slope limiters, and cross-derivative terms proportional to $\mathcal{A}\mathcal{B}$ and $\mathcal{B}\mathcal{A}$. The latter come about because flow at an angle to the grid may take material from cell $(i,j)$ directly to for example cell $(i+1,j+1)$: there can be transport across the corners of the grid.

The simplest way of dealing with these cross-derivatives is by using dimensional splitting: treat each dimension separately, varying the order in such a way as to minimize the splitting error (e.g. Strang, 1968). Alternatively, one can adopt an unsplit method such as the Corner Transport Upwind method (CTU, Colella, 1990), which then has to take care of the cross-derivatives directly, an approach that is taken for example in ATHENA (Stone et al., 2008), while PLUTO (Mignone et al., 2007) offers both options. Note that neither of these options are available for unstructured meshes, which makes a code such as arepo (Springel, 2010) formally only first-order accurate, although in practice it shows second-order convergence (e.g. Pakmor et al., 2016).

While both dimensional splitting and CTU offer formal second-order accuracy for smooth flows, in some cases it can be advantageous to adopt a fully multidimensional approach, in particular in regions where the flow is not extremely well-resolved (Balsara, 2010). One concern is that in many implementations, higher-order corrections are basically one-dimensional: when deciding if the interaction between cell $(i,j)$ and $(i-1,j)$ can be second order, information from cells $(i-2,j)$ and $(i+1,j)$ is used, all at the same value of $j$ (e.g. LeVeque, 2002). It is therefore interesting to look at alternatives, such as fully multidimensional Riemann solvers (Balsara, 2010). Here, we explore multidimensional upwind methods in the framework of residual distribution.

3.1.4 Residual distribution formulation

Before diving into the residual distribution framework, we first show how the familiar one-dimensional Roe solver can be formulated as a residual distribution scheme. This sets the scene for the next few sections.

A different way of looking at the one-dimensional linear problem ($\mathcal{A}$ does not depend on ${\bf W}$) involves defining the residual associated with cell interface $i-1/2$:

[TABLE]

If the residual is zero, the flux is constant and there should be no evolution of the state. Note that, in Fig. 5, the left state ${\bf W}_{i-1}$ is only modified by the jump associated with a negative propagation speed, and therefore a negative eigenvalue of $\bar{\mathcal{A}}$, while the right state is modified by the two jumps associated with positive eigenvalues of $\bar{\mathcal{A}}$. The residual is therefore split between the two neighbouring cells according to the sign of the eigenvalues of $\bar{\mathcal{A}}$666LeVeque (2002) uses the term fluctuation splitting rather than residual distribution.. This can be appreciated even more when writing the update (21) in matrix form, introducing $\Lambda^{+}$ as the diagonal matrix with $\lambda^{+}$ as entries on the diagonal:

[TABLE]

where $\mathcal{R}$ denotes the matrix of right eigenvectors of $\bar{\mathcal{A}}$. The contribution of the Riemann problem at interface $i-1/2$ to cell $i-1$ is, using the same notation:

[TABLE]

Note that the contributions to the two neighbouring cells (31) and (32) sum up to $\phi$. The residual $\phi$ is split, or redistributed, amongst the neighbouring cells, according to the recipe:

[TABLE]

This formulation of the Roe solver is interesting because it has higher dimension counterparts, which we introduce next.

3.2 Residual distribution basics

Consider a system of $q$ hyperbolic conservation laws in $d$ spatial dimensions:

[TABLE]

Write in quasi-linear form, introducing a parameter vector ${\bf Z}$, to be specified later:

[TABLE]

Now we introduce our triangular mesh. Assume we know ${\bf Z}$ at the vertices (nodes), and since the nodes are connected by triangles there exists a piecewise linear interpolation of the nodal values:

[TABLE]

where $N$ is the total number of nodes in the mesh and $\omega_{i}^{h}$ is the piecewise linear shape function equal to unity at node $i$ and vanishing outside of the triangles sharing $i$ as a vertex.

The triangle residual is given by (cf. the one-dimensional case (30)):

[TABLE]

where $\mathcal{A}_{j}=\partial{\bf F}_{j}/\partial{\bf Z}$. For the piecewise linear interpolation (36) we have that

[TABLE]

where $|T|$ denotes the area of triangle $T$, ${\bf n}_{i}$ is the inward pointing normal to the edge opposite node $i$ of triangle $T$ and ${\bf\hat{x}}_{j}$ is the unit vector in direction $j$. Plugging this into (37) yields

[TABLE]

with

[TABLE]

with

[TABLE]

In order for the resulting scheme to be conservative, the average matrix $\bar{\mathcal{A}}_{j}$ must obey above relation given ${\bf Z}^{h}$. While this is difficult in general, if the fluxes ${\bf F}$ are at most quadratic functions of ${\bf Z}$, this means that the entries of $\mathcal{A}_{j}$ are at most linear in ${\bf Z}$, making the integrals trivial to evaluate so that the average matrix $\bar{\mathcal{A}}_{j}$ is just $\mathcal{A}_{j}$ evaluated at the nodal average of ${\bf Z}$:

[TABLE]

For the Euler equations, a parameter vector can be found that leads to a quadratic flux function, which is basically a multidimensional analogue of Roe’s original parameter vector (Deconinck et al., 1993).

The hyperbolic nature of equation (34) guarantees that $\mathcal{K}_{i}$ has $q$ real eigenvalues and a complete set of linearly independent eigenvectors. Diagonalization yields

[TABLE]

where $\mathcal{R}_{i}$ is a matrix whose columns are the right eigenvectors of $\mathcal{K}_{i}$, $\mathcal{L}_{i}=\mathcal{R}_{i}^{-1}$, and $\Lambda_{i}$ is a diagonal matrix containing the eigenvalues of $\mathcal{K}_{i}$. We can now define the multidimensional upwind parameter as

[TABLE]

with

[TABLE]

Explicit expressions for the matrix elements are provided in appendix A.

Residual distribution schemes take the cell residual, and, as the name suggests, redistributes it over neighbouring nodes (cf. the one-dimensional case (33)). They take into account the hyperbolic nature of the equations by looking at upwind directions, and do so in a multidimensional way by looking at the eigenvalues of $\mathcal{K}$ rather than $\mathcal{A}$ and $\mathcal{B}$ separately. All nodes gather parts of the residuals of all triangles that have that specific node as a vertex. This process is illustrated in Fig. 6. Triangle $T$ has residual $\phi^{T}$, which is redistributed over the three vertices of $T$. The total update at vertex $a$ is the sum of the contributions of the seven triangles that have $a$ as a vertex. Defining $\phi^{T}_{i}$ as the part of the residual of triangle $T$ to be sent to node $i$, an update would look like

[TABLE]

where $n$ denotes the number of time steps taken so far, and $V_{i}$ is the volume associated with node $i$. For a Delaunay triangulation, $V_{i}$ is the volume of the Voronoi cell centred on $i$ (see Fig. 6). For stability, the time step is limited by a CFL condition:

[TABLE]

where the minimum is taken over all vertices in the triangulation, $l^{T}_{\mathrm{max}}$ is the length of the longest edge of $T$ and $\lambda^{T}_{\mathrm{max}}$ is the maximum possible signal speed. For the Euler equations,

[TABLE]

where ${\bf v}_{j}$ and $c_{j}$ are the velocity and sound speed at vertex $j$.

The update (46) is only first order accurate in time. Higher order temporal accuracy is possible but will depend on the exact scheme used and will be discussed in section 3.4. A residual distribution scheme is defined by how it defines $\phi_{i}$, or, in other words, how the residual is distributed over the neighbouring nodes. Below, we discuss two possible choices for the distribution function.

3.3 Distribution schemes

Traditionally, residual distribution schemes have been used to find steady solutions to the Euler equations, with (46) used only to reach the required steady state. In this case, temporal accuracy is not an issue. The final steady state will depend on the distribution coefficients $\phi_{i}$. Several design criteria have been identified (e.g. van der Weide, 1998), of which we have encountered two already: conservation and multi-dimensional upwinding. The former sets the linearization (42), while the latter was introduced through the use of $\mathcal{K}$ rather than using $\mathcal{A}$ and $\mathcal{B}$ separately. All schemes discussed below are both conservative and multidimensional upwind.

Two other important considerations are positivity, or monotinicity, and linearity preservation. A scheme is said to be positive when no new extrema are introduced in the solution when going from one time step to the next. It is therefore related to the concept of total variation diminishing (LeVeque, 2002), and is especially important in compressible flows since a positive scheme does not introduce oscillations near discontinuities. If a scheme is linearity preserving it means that exact linear solutions are recovered by the scheme. In a steady state, this means such a scheme is second-order accurate in space (Abgrall, 2001).

A scheme is said to be linear if, when applied to a linear partial differential equation such as

[TABLE]

the solution update can be expressed as

[TABLE]

where $N$ is the total number of grid points and the coefficients $c_{j}$ are independent of $u$. For example, the one-dimensional first-order Roe scheme is linear. Unfortunately, as a consequence of Godunov’s theorem, a linear scheme can not be both positive (monotone) and linearity preserving (second-order) (Struijs, 1994). Just as in the one-dimensional case, non-linear schemes have to be designed in order to get the best of both worlds.

3.3.1 Linear $N$ scheme

The $N$ scheme ($N$ for narrow, Struijs, 1994; van der Weide, 1998) is a monotonic scheme that is at most first-order accurate in space. The distribution function is given by

[TABLE]

where

[TABLE]

There are certain cases for which the inverse matrix $\mathcal{\hat{N}}$ does not exist, for example at stagnation points. However, the product $\mathcal{K}^{+}_{i}\mathcal{\hat{N}}$ always has meaning, making the $N$ scheme always well-defined.

3.3.2 Linear $LDA$ scheme

While the $N$ scheme can deal with shocks in a stable and satisfactory way, its first-order nature makes it too diffusive in smooth flows to be of practical use. A popular second-order scheme is $LDA$ (Low Diffusion A, Struijs, 1994; van der Weide, 1998), for which the distribution function is given by:

[TABLE]

with distribution coefficients

[TABLE]

Using the $LDA$ scheme in the presence of discontinuities leads to unphysical oscillations, as expected.

3.3.3 Non-linear blended schemes

In order to get the best of both worlds, second-order accuracy in regions of smooth flow while remaining monotone in the presence of discontinuities, schemes that blend the $N$ and $LDA$ residue have been designed:

[TABLE]

where $\Theta^{E}$ is a diagonal non-linear blending matrix

[TABLE]

where subscript $k$ indicates the $k$th equation of the system (e.g. Csík et al., 2002). This scheme will be referred to as the $B$ (for blended) scheme. Two useful variants can be obtained by setting all diagonal values of $\Theta^{E}$ to $\theta$, there $\theta$ can be taken to be either the maximum ($B$max) or the minimum ($B$min) over all $\Theta_{k,k}$. $B$max favours the $N$ scheme when in doubt, and is therefore a good choice when strong shocks are present in the solution, while $B$min favours the $LDA$ scheme, and is therefore a better choice for solutions that are relatively smooth.

A different blending scheme was proposed in Dobes & Deconinck (2008), where the blending coefficient is taken to be to be a scalar $\theta$:

[TABLE]

which is based on a shock sensor $s$:

[TABLE]

where $h$ is a measure of the size of the element and

[TABLE]

where $L$ is the domain size, is a sensor that is non-zero only in regions of compression and is $O(1)$ in regions of smooth flow, which makes $\theta$ $O(h)$ and the scheme second-order accurate (Dobes & Deconinck, 2008). This scheme will be referred to as the $Bx$ scheme.

3.4 Second order temporal accuracy

When time-dependent problems are considered, temporal accuracy becomes an issue. While there exists a residual distribution scheme that is second-order in space and time, the Lax-Wendroff scheme (cf. section 3.1.2), this has to be mixed with a first-order scheme whenever discontinuities are present, and such schemes have been only moderately successful (März & Degrez, 1996; Ferrante & Deconinck, 1997; Hubbard & Roe, 2000).

The second method for improving the order of accuracy of section 3.1.2, the method of lines, has seen much more progress in recent years. Unfortunately, a straightforward implementation of the method of lines is impossible due to the fact that residual distribution schemes in their basic formulation suffer from an inconsistent spatial discretization (März & Degrez, 1996). For linearity preserving schemes such as LDA, a consistent update can be obtained through a Petrov-Galerkin formulation well known from finite element analysis , but for positive schemes another approach is needed (Abgrall & Mezine, 2003).

In the time-dependent case, the total residual $\Phi^{T}$has both a space and a time component:

[TABLE]

As for the steady case, the question is how to redistribute this residual over the nodes. In order for the scheme not to suffer from an inconsistent discretization we introduce a mass matrix $m_{ij}$, well known from finite element analysis:

[TABLE]

The particular form of $m_{ij}$ is discussed below. A time step can now be taken by requiring that for every node

[TABLE]

For a general choice of $m_{ij}$, this leads to an implicit solver. Explicit schemes were derived only very recently (Rossiello et al., 2009; Ricchiuto & Abgrall, 2010). In the following, we present the second-order scheme of Ricchiuto & Abgrall (2010), which is a two-stage Runge-Kutta method:

[TABLE]

where we have indicated specifically that the spatial residue in the first step is based on the interpolation ${\bf W}^{h}$ at time level $n$. The total residue $\Phi_{i}$ appearing in the second step is for the $N$ scheme given by

[TABLE]

while for the $LDA$ scheme it reads:

[TABLE]

Note that this integration scheme is closely related to total-variation-diminishing time integration schemes (Shu, 1988). For the $N$ residual, the mass matrix is simply

[TABLE]

which is not consistent but since the $N$ scheme is only first order anyway this is not a problem. Several choices can be made for the mass matrix in the $LDA$ residue (Ricchiuto & Abgrall, 2010), from which we choose777While different versions of $m_{ij}$ can formally be ranked according to their dissipative nature, much less is known about their stability. The choice of (68) is based on simplicity and apparent stability.

[TABLE]

where $\beta_{i}$ is given by (54). A slightly more complicated update, called selective lumping (Ricchiuto & Abgrall, 2010) adds an anti-diffusive term to the residual:

[TABLE]

where $\Phi_{i}$ can be either $\Phi^{LDA}_{i}$ or $\Phi^{N}_{i}$ and $m^{G}_{ij}$ is the Galerkin mass matrix

[TABLE]

For blended schemes, the $N$ and $LDA$ total residue are mixed in the usual way (see section 3.3.3). The resulting scheme is second order accurate in space and time wherever the blending procedure favours the $LDA$ scheme (Ricchiuto & Abgrall, 2010).

To summarise, a complete integration scheme is defined by the residual distribution scheme ($N$, $LDA$, $B$, etc), a mass matrix, and a choice of lumping. In all test problems discussed below, we stick to the mass matrix of equation (68) together with global lumping. The only choice left is the residual distribution scheme, with the annotation that when using the first order $N$ scheme we also use a first order time integration.

3.5 Boundary conditions

On structured grids, boundary conditions are often imposed by adding a layer of ghost cells to the computational domain, whose states are set in such a way to achieve the desired boundary condition. Periodic boundary conditions, for example, are simply achieved by copying the relevant states from the other side of the computational domain into the ghost cells. A reflecting boundary can be achieved by copying the states next to the boundary from the computational domain into the ghost cells but reversing the velocity normal to the boundary. One reason why this approach is very effective in the case of structured grids is that the boundaries always align with one of the coordinate axis. A second reason is that all cells have the same shape and volume so that copying the state is trivial888Exceptions include for example structured grids in curvilinear coordinates, where cells at different spherical radii $r$ will have different volumes $\propto r^{2}$..

For an unstructured grid, we do not necessarily have the boundary aligned with one of the coordinate axis, and, since all computational cells are slightly different, having a layer of ghost cells would mean copying part of the grid structure, which is expensive both in terms of computational effort and memory requirement, and should therefore be avoided. An exception to this rule are non-reflecting boundaries, where the boundaries are taken to be so far away from the region of interest that their exact shape does not matter. In this case, we promote the boundary vertices to ghost vertices, whose states never changes from the initial conditions, which are taken to be a stationary state. All vertices connected to the ghost vertices, because of multidimensional upwinding, ‘see’ waves through the usual characteristic decomposition. A wave trying to leave the computational domain can do so, but if information needs to be drawn from outside the computational domain because one of the characteristics points inward, this information is drawn from the ghost cells containing no wave. Therefore, no waves enter the computational domain and therefore we call these boundary conditions non-reflecting, and they are relatively trivial to implement: at the start of a time step, all boundary cells have to be set to the initial condition. This procedure can also be used to specify an inflow boundary.

Periodic boundary conditions are completely handled by the mesh. If the mesh is periodic in both $x$ and $y$, all vertices and triangles have neighbours in all directions (see Fig. 2). Therefore, for periodic boundaries, no extra work is required.

The final boundary discussed here is a reflecting wall. Such a boundary is most easily enforced in a weak sense, which means specifying the flux across the boundary rather than the state at the boundary. For the implementation we follow van der Weide (1998). Consider a triangle of which one edge $e$ with vertices $a$ and $b$, belongs to a reflecting wall, and let the normalized normal of this edge be ${\bf n}=(n_{x},n_{y})^{T}$. The desired flux at a reflecting wall is such that the velocity normal to the wall vanishes: $v_{n}={\bf v}\cdot{\bf n}=0$. Usually the flux as computed as if there was no wall does not obey this condition, and therefore a correction flux has to be applied:

[TABLE]

The correction residual is then given by

[TABLE]

where the trapezium rule was used to approximate the integral. This residual is distributed over the nodes $a$ and $b$ using a parameter $\alpha\in[0,1]$:

[TABLE]

Following van der Weide (1998), we choose $\alpha=0.75$.

4 GPU implementation

GPUs have emerged relatively recently as viable alternatives to large distributed-memory machines. Modern single GPU cards are capable of Teraflops performance, comparable in speed to a CPU cluster of a few $100$ cores but at a fraction of the cost. However, getting close to peak performance of a GPU is not straightforward, even though in recent years is has become much easier.

The computational intensity of numerical fluid dynamics made it a prime candidate to be ported to GPUs. Both SPH (Hérault et al., 2010) and structured grid-based methods (Hagen et al., 2006) were successfully run on GPUs, with unstructured grid methods for magneto-hydrodynamics following later (Lani et al., 2014). The latter was based on an implementation of the Rusanov flux (Rusanov, 1961). A method specific for turbulence calculations on hybrid grids was presented in Asouti et al. (2011). To the best of our knowledge, astrix is the first implementation of an explicit residual distribution method on GPUs.

astrix is written using NVIDIA’s CUDA (Compute Unified Device Architecture) programming model. CUDA-capable cards come in different generations or compute capabilities. The higher the compute capability, the newer the card and this means more features may be available. The first generation of CUDA-capable cards were built to comply with single precision IEEE requirements999The first generation CUDA cards with compute capability $1.$x does not support denormal numbers, and the precision of division and square root operations are slightly below IEEE 754 standards (Whitehead & Fit-Florea, 2011). Cards with compute capability $2.$x and higher do not suffer from these issues., which, as we saw turns out to be important for generating unstructured meshes as at several stages we need exact geometric predicates, while newer cards (compute capability $2.$x and higher) fully support double precision arithmetic.

Designing algorithms for the GPU is fundamentally different from designing traditional parallel CPU algorithms, simply because the GPU is a very different beast. A useful analogy is that the CPU is a single genius, while the GPU represents millions of unskilled labourers. The GPU gets its performance not by doing single computations fast, but by taking on a lot of single computations at the same time and switching between them if necessary. Of course, this can only be done if the computations are independent. Therefore, we need to expose enough parallellism: we need to load the GPU with millions of independent relatively small tasks in order to get close to peak performance. This is done by launching many threads (typically millions) that each work independently. These threads are organised into warps that can be moved between computing units very efficiently, thereby hiding instruction and memory latency. Exposing enough parallellism is straightforward in the residual distribution schemes discussed in section 3: the bulk of the computations can be done independently for all triangles, and one expects a two-dimensional mesh to contain typically millions of triangles. The situation is more tricky for mesh generation as will be explained in section 4.1 below.

A second condition for good GPU performance is adequate use of GPU memory. Data transfer from CPU to GPU is slow since it has to go through the PCI bus (typical speed $6$ GB/s). This has to be kept in mind when porting only part of an application to GPU. While the GPU may speed up a particular computation 100 times, if this computation is preceded and followed by a few seconds of data transfer the overall speedup may be negligible or even negative. The converse is also true: if a particular part of an algorithm is not very well suited to run on a GPU, the benefit of running this part on the CPU may be outweighed by the increased memory traffic. astrix is designed with minimal CPU-GPU memory transfers: all data reside on the GPU and stay on the GPU, unless explicit output is required.

Finally, this brings us to performance metrics. From a user’s point of view, it is important to know how much the code speeds up when running on the GPU compared to the CPU. Unfortunately, this is extremely dependent on the GPU/CPU combination. Moreover, one should compare an algorithm optimised for the GPU to an algorithm optimised for the CPU, and usually these are very different algorithms because the GPU works differently from a CPU. A fair comparison therefore requires designing and optimising two different algorithms performing the same task. While for completeness we do mention GPU/CPU speedups when measuring performance of astrix, above considerations should be kept in mind.

4.1 Mesh generation

As discussed above, efficient use of a GPU is non-trivial. In order to expose as much parallelism as possible, most steps in the Delaunay refinement algorithm are done launching one CUDA thread per element. For example, when finding low-quality triangles, each thread will check a single triangle if the quality and size constraints are met. This leads to a list of vertices to add to the mesh, for which we can find their containing triangles independently, and also check independently whether they encroach upon any segment.

4.1.1 Parallel insertion set

Inserting new vertices requires more attention, since not all vertices can be inserted independently. For example, two new vertices may find themselves in the same triangle, in which case only one of them can be inserted at a time. If a new vertex is to be inserted on an edge, no vertices can be inserted in the neighbouring triangles. A more stringent constraint arises from the demand that it must be possible to generate the mesh by inserting the vertices one at a time. This is important because all proofs of quality and termination of the algorithm rely on inserting new vertices sequentially. Therefore, we have to select a subset of circumcentres that can be inserted independent from each other.

It is obvious that a single new vertex $v$ will only affect those triangles in the mesh whose circumcircles contain $v$. All other triangles have empty circumcircles and will therefore be part of the new triangulation. This leads to the definition of the cavity of $v$: the set of triangles whose circumcircles contain $v$. Then two vertices can be inserted independently from each other if their cavities do not overlap. An example of a cavity is shown in Fig. 7.

Finding all triangles belonging to the cavity of $v$ is relatively straightforward. We already have found the insertion triangle (see section 2.4.2). Starting from this triangle, which is of course part of the cavity, move into one of the neighbouring triangles if it also belongs to the cavity. By always checking the neighbours in anti-clockwise order, it is possible to find all triangles belonging to the cavity in relatively few steps.

A simple algorithm for selecting non-overlapping cavities is the following. For every new vertex $v$, flag all triangles belonging to its cavity with a unique integer $i(v)$. If a triangle has already been flagged with $i^{\prime}(v^{\prime})$, take $\mathrm{max}(i,i^{\prime})$ as the flag so that a triangle will always be associated with at most one new vertex, with higher values of $i(v)$ being given priority. Taking the maximum as mentioned above involves three steps: 1) reading the current value of the flag, 2) comparing it with the new value, 3) write back the new value if it is higher than the old value. In order to prevent a race condition, these three operations have to be done without interference from other GPU threads, which can be done within CUDA through so-called atomic functions, in this case atomicMax. Once all new vertices are processed in this way, we walk through the cavities a second time, and check for every new vertex $v$ if all triangles belonging to its cavity are flagged with $i(v)$. If so, the vertex can be inserted.

Of course, we want to insert as many vertices in one parallel step as possible. There exists what is called a maximal independent set (e.g. Luby, 1986); a maximum number of vertices that can be inserted in a single parallel step. While Spielman et al. (2007) present an algorithm for calculating the maximal independent set for this specific problem, the simple algorithm presented above performs remarkably well if the integers $i(v)$ are chosen to be random. That is, if there are $N$ bad triangles in the mesh, and therefore potentially $N$ vertices to insert $(v_{1}...v_{N})$, give each vertex a unique random integer as its $i(v)$. This leads to independent cavities covering the whole domain relatively uniformly.

The algorithm presented in Spielman et al. (2007) takes $O(\log^{2}(L/s))$ iterations, where $L$ is the domain size and $s$ the smallest circumradius in the final mesh. In Fig. 8, we compare our simple algorithm to this theoretical limit, both for regular integer assignment (i.e. $i(v_{n})=n$) and random integer assignment. Note that because of Morton ordering (see section 4.1.4 below), vertices with large values of $i$ are located very close to each other, making the parallel selection method very inefficient. However, for random integer assignment, the number of iterations necessary follows the theoretical limit rather nicely. In addition, finding the cavities associated with the new vertices has the additional advantage that we know which edges may need flipping: only edges part of a cavity of a new vertex plus any newly created edges need to be checked for the Delaunay property. This saves a lot of redundant checking of edges.

4.1.2 Data structure

The mesh contains nVertex vertices, nTriangle triangles and nEdge edges. The basic structure of the grid is stored in four arrays, using CUDA intrinsics:

• •

vertexCoordinates; a float2 array of size nVertex containing the $x$ and $y$ coordinates of all vertices in the mesh.

• •

triangleVertices; an int3 array of size nTriangle containing for every triangle the three vertices that make up the triangle.

• •

triangleEdges; an int3 array of size nTriangle containing for every triangle its three edges.

• •

edgeTriangles; an int2 array of size nEdge containing for every edge the two neighbouring triangles (or just one if the edge is part of the boundary and therefore a segment).

4.1.3 Exact geometric predicates

As indicated in the previous sections, at several stages (finding insertion triangles and testing edges for the Delaunay property) we need exact geometric predicates, i.e. the exact sign of the determinants Orient2D (3) and InCircle2D (8). While this can be done in principle using exact arithmetic, the price is quite high: up to two orders of magnitude reduction in speed. Fortunately, an adaptive method was designed by Shewchuk (1997), based on earlier work by Priest (1991). The key insight is that the exact determinant is not needed: all we are interested in is the sign. If we can be sure that a calculation at finite precision gives the correct sign, there is no need to make it more precise. These algorithms work on most processors, in particular those complying to the IEEE 754 standard, and can therefore be ported in a straightforward way to modern GPUs.

4.1.4 Morton ordering

Data locality has always been critical for efficient use of GPUs. On older cards (compute capability 1.x), when reading an array from global device memory, it was critical for neighbouring threads to read neighbouring data: if thread 0 reads $\texttt{array}[0]$, it was necessary for thread 1 to read $\texttt{array}[1]$ and so on; any other order would incur a speed penalty of up to 2 orders of magnitude. More recent GPUs have relaxed these requirements by introducing on-chip cache, but this still means that data locality is highly desirable: if neighbouring threads read data that is close together in memory, chances are it can be found in the cache, which means a read from global device memory is unnecessary.

Unstructured meshes pose a challenge for maintaining data locality due to the non-trivial interconnections between vertices. Moreover, in the process of creating the mesh new vertices, edges and triangles are added, quickly destroying data locality even if it was present at some stage. In order to mitigate this, after every parallel insertion step, we reorder vertices, edges and triangles as to maintain as much data locality as possible. This is done by assigning a Morton value (Morton, 1966) to for example each vertex, and then sorting the vertices according to their Morton value. The same for edges and triangles. While sorting itself is non-local and not trivial to implement of a GPU, CUDA has fast built-in sorting algorithms so that the overall effect on execution speed of Morton ordering is positive.

4.1.5 Delaunay triangulator

For efficient use of the GPU in maintaining a Delaunay triangulation, we want to flip as many edges as possible in parallel. First, we use the robust InCircle2D (8) test to generate a list of edges that do not satisfy the Delaunay requirement. From this list, edges can be flipped independently if they are not part of the same triangle. We select an independent set in much the same way as done in section 4.1.1, where the ‘cavity’ of an edge is now defined by the two neighbouring triangles. Randomization was not found to be necessary in this case, since because the ‘cavities’ are so small a large independent set can always be found.

Unfortunately, edge flipping can corrupt the data structure, in particular edgeTriangles (triangleVertices and triangleEdges are updated during the flip). While the flipped edge still has the same neighbouring triangles $t_{1}$ and $t_{2}$, other edges belonging to for example the original $t_{1}$ may suddenly have $t_{2}$ rather than $t_{1}$ as a neighbour. Fortunately, this is straightforward to correct. For every edge $e$, we look at its neighbouring triangles $t_{1}$ and $t_{2}$ through edgeTriangles. If none of the edges of $t_{1}$, as per triangleEdges, is equal to $e$ then this means that $e$ now neighbours $t_{2}$ rather than $t_{1}$. Therefore, before flipping edges, we create a triangle substitution array triangleSub so that for every edge to be flipped $\texttt{triangleSub}[t_{1}]=t_{2}$ and vice versa. Note that no conflicts can arise since any triangle can only be associated with one edge that will be flipped (otherwise these edges can not be flipped in parallel). After a parallel step of edge flipping, we can then replace $t_{1}$ or $t_{2}$ with $\texttt{triangleSub}[t_{1}]$ or $\texttt{triangleSub}[t_{2}]$, respectively, where necessary. See Navarro et al. (2011) for more details.

4.1.6 GPU performance

As a test case, we consider the generation of a uniform unstructured mesh, periodic in both $x$ and $y$, with $1.3$ million vertices. We compare a GPU version to a CPU version, using exactly the same algorithms using single precision floating points. The test was run on a system consisting of an Intel Xeon $1.8$ GHz CPU and a NVIDIA Tesla K20m GPU, which has CUDA compute capability $3.5$.

In Fig. 9 the total time spent on the CPU in each ‘kernel’101010When running on the CPU a kernel is replaced by a for-loop performing exactly the same task is shown. Most time is spent checking if edges are Delaunay, followed by the quality check of triangles. It should be noted that the kernel TestQuality is called once every refine cycle, while the kernel CheckEdge is called multiple times in the same cycle until the mesh is Delaunay. The time spent in individual instances of CheckEdge is actually smaller than that for CheckTriangle, but the number of kernel calls make CheckEdge the most time-consuming kernel.

The corresponding timings for the GPU are shown in Fig. 10, in the same order from left to right as Fig. 9. Checking edges for Delaunay-hood is still the most expensive operation, but the costs have been reduced by a factor of $\sim 80$ compared to the CPU version. While this may seem as a healthy speedup, it is nowhere near the maximum capability of the Tesla K20m GPU, as we will se below. This is partly due to the fact that many calls to CheckEdge involve only very few edges to be processed, hence limiting the parallelization. The maximum number of edges to be checked in a single kernel call is $\sim 650000$ for this particular mesh, giving a speedup of $\sim 100$ compared to the CPU. The brute force approach of checking all edges in every Delaunay iteration therefore gives a bigger speedup, but the overall computational costs would still increase. A second reason for the relatively poor performance is that the kernel requires a lot of memory traffic. For every edge checked, we need to know all coordinates and all edges of the two neighbouring triangles, and because of the unstructured nature of the grid the memory access involved is not ideal for the GPU. The third reason is that the kernel has to make use of exact geometric predicates, which first of all makes the algorithm more complicated, which increases the number of registers used and therefore limits the amount of blocks that can be run simultaneously, and at the same time leads to warp divergence: the amount of computation performed can differ significantly for different edges.

On the other hand, the kernel TestQuality show a much better speedup, from $\sim 220$ on average to $\sim 270$ maximum. The achieved bandwidth of $120$ Gb/s comes reasonably close to the theoretical maximum of the Tesla K20m GPU of $200$ Gb/s, given the unfavourable memory access pattern due to the unstructured nature of the mesh. Nevertheless, even here there is room for improvement, although the focus should of course be on the CheckEdge kernel.

The other kernels worth mentioning are TestEncroach, LockTriangle, IndependentCavities and FlagEdges. These have in common that they walk through the grid around a certain point, visiting an unknown number of triangles, for example the insertion cavity in the case of LockTriangle. These kernels show the worst speedup on the GPU ($\sim 50$), first of all for similar reasons as CheckEdge mentioned above. In addition, there is the extra complication of insertion cavities having different sizes, which leads to different work loads for different GPU threads. Moreover, the size of the cavity is unknown beforehand, making optimisations more difficult for the compiler.

Overall, the creation of the $1.3$ million vertex mesh has sped up by roughly a factor of $100$ compared to the CPU, on a graphics card that costs only a fraction of a CPU compute cluster, making the effort of specialising to the GPU worthwhile.

4.2 Hydrodynamics

4.2.1 Residual distribution

The two-stage Runge Kutta update (63)-(64) consists of four steps:

• •

Calculate $\phi({\bf W}_{h}^{n})$, $\phi_{i}^{N}({\bf W}_{h}^{n})$ and $\phi_{i}^{LDA}({\bf W}_{h}^{n})$

• •

Blend into $\phi_{i}({\bf W}_{h}^{n})$ and calculate ${\bf W}^{*}$

• •

Calculate $\phi({\bf W}_{h}^{*})$, $\Phi_{i}^{N}$ and $\Phi_{i}^{LDA}$

• •

Blend into $\Phi_{i}$ and calculate ${\bf W}^{n+1}$

These steps are distributed over the following GPU kernels:

• •

CalcResidual: calculate $\phi({\bf W}_{h}^{n})$, $\phi_{i}^{N}({\bf W}_{h}^{n})$, $\phi_{i}^{LDA}({\bf W}_{h}^{n})$

• •

AddResidual: blend into $\phi_{i}({\bf W}_{h}^{n})$ and calculate ${\bf W}^{*}$

• •

CalcTotalResNtot: calculate $\phi({\bf W}_{h}^{*})$, $\Phi_{i}^{N}$

• •

CalcTotalResLDA: calculate $\Phi_{i}^{LDA}$

• •

AddResidual: blend into $\Phi_{i}$ and calculate ${\bf W}^{n+1}$

In addition, there are kernels for calculating the allowed time step, the parameter vector and to set the boundary conditions, but the vast majority of the computational time is spent in the kernels mentioned above. Note that the kernel AddResidual performs exactly the same task twice but with different residuals.

Calculations of the residuals are independent for each triangle and can therefore be parallelized very efficiently. The node updates involve the contributions from all triangles sharing a particular node. While this could be parallelized over the nodes, we do not have direct information on which triangles share for example node $i$ from the mesh data structure. This would be difficult to achieve, since the number of triangles per node can vary quite a lot. It would be possible to assign one triangle to every node, and walk around the node collecting the contribution from all triangles sharing the node, but since all nodes have to be updated we found it more efficient to again paralellize over triangles and update the nodal values using atomic operations.

4.2.2 GPU performance

As a test case, we consider the mesh generated in section 4.1.6 and consider the cost of a single time step, using the same CPU/GPU combination as in section 4.1.6. The results for the CPU are shown in Fig. 11. It is clear that the bulk of the computational time is spent calculating the residuals. Comparing with Fig. 9, we see that the cost of setting up the mesh is roughly $100$ time steps. While in this paper we are considering static meshes only, this high cost of generating the mesh should be kept in mind when contemplating dynamic meshes. We will see below that this issue is even more important on the GPU.

In Fig. 12 we show the corresponding timings for the GPU. Now most time is spent in adding the residuals to the vertices. This kernel shows a modest speedup with respect to the CPU of roughly $40$.This is not because the kernel makes inefficient use of the GPU: the achieved bandwidth is $\sim 150$ Gb/s, which is good compared to the theoretical maximum of $200$ Gb/s considering that all additions have to be atomic. Rather, it is the relatively low amount of computations compared to the memory traffic in this kernel that limits the speedup.

The situation is much better for the computationally intensive kernels calculating the residuals. The kernels CalcResidual and CalcTotalResLDA show a speedup of $500$ at a bandwidth of $110$ Gb/s and $130$ Gb/s, respectively, while CalcTotalResNtot shows a speedup of $320$ at a bandwidth of $135$ Gb/s. This shows the true power of the GPU as these kernels do a lot of computations per memory element. Overall, this makes a GPU time step $250$ times faster compared to the CPU. Since this skews the performance on the GPU towards the hydrodynamics compared to the generation of the mesh, a dynamic mesh will hurt performance more on the GPU than on the CPU.

5 Test problems

In this section we discuss the performance of astrix in various standard test problems. Since the schemes implemented in astrix can be seen as multidimensional variants of the Roe solver, we will use a dimensionally split version of the Roe solver, as implemented in rodeo (see e.g. Paardekooper & Mellema, 2006; Paardekooper, 2012), as a benchmark. In addition, we also show results obtained with a different approximate Riemann solver (HLLC, Toro et al., 1992). In all cases, except where explicitly mentioned, the minmod limiter was used as a flux limiter.

5.1 One-dimensional tests

We begin by discussing some standard one-dimensional test problems. Since astrix is a purely multidimensional method, we have to choose a two-dimensional mesh that does not allow for variations in the state in one of the coordinate directions. This can be achieved by making a structured mesh periodic in $y$ with one cell covering the whole domain, see Fig. 13. The ‘one-dimensional’ domain is located at $y=0$, and the domain is periodic in $y$ with period $0.1$. All vertices are located either at $y=0$ or at $y=\pm 0.1$, which means all vertices map onto $y=0$, making the calculation effectively one-dimensional.

5.1.1 Linear sound wave

The first problem consists of a linear sound wave of amplitude $10^{-4}$, in a uniform medium of $\rho=1$, ${\bf v}={\bf 0}$ and adiabatic sound speed $c=1$. The ratio of specific heats is set to $\gamma=1.4$. The domain size in $x$ as well as the wavelength of the perturbation is set to unity. Boundary conditions are periodic, and the wave is evolved to $t=1$. We compute the $L_{1}$ error norm in the density as

[TABLE]

where $\rho_{i}$ is the density in the $i$th cell after the final time step, $\rho_{0,i}$ is the exact solution at that particular location, $A_{i}$ is the volume of the $i$th cell, and $A=\sum_{i}A_{i}$ is the total volume of the computational domain. Note that the scaling with volume effectively makes this an estimate for the mass error. In the case of a 1D mesh, all cells have the same volume and the scaling with volume has no effect.

The results for three first-order schemes, the $N$ scheme, the first order Roe scheme and the first order HLLC scheme are shown by the solid curves in Fig. 14. The thin solid lines indicate errors $\propto h$ and $\propto h^{2}$. The differences between the first order schemes are so small they are hardly visible, and in particular the Roe scheme and the HLLC scheme give the same $L_{1}$ error up to three significant digits. For high enough resolution, all schemes show linear convergence as expected. Second order Roe and HLLC, shown in Fig. 14 by dashed curves are again indistinguishable. The same is true for the two second order astrix schemes $B$ and $Bx$, although they do show smaller errors than both Roe and HLLC. Second order convergence is obtained for all second order schemes around $h\sim 0.01$. Of course, the linear sound wave itself is only an approximation to the true solution, and towards higher resolution, the errors become dominated by the departure from linearity, which removes the second order convergence for all schemes. We note that results obtained with the $LDA$ scheme are indistinguishable from those obtained with both $B$ and $Bx$.

5.1.2 Sod shock tube

This one-dimensional Riemann problem is a well-known test case for gas dynamics codes (Sod, 1978). The initial conditions consist of two constant states, separated by a membrane at $x=0.5$. The left state has density $\rho_{L}=1$, pressure $p_{L}=1$, while the right state has density $\rho_{R}=0.125$ and pressure $p_{R}=0.1$. The velocity is zero everywhere initially, and the ratio of specific heats is set to $\gamma=1.4$. The domain $0<x<1$ is covered by $100$ grid cells, and the solution is evolved until $t=0.2$.

The results are shown in Fig. 15. The analytic solution, shown with the dashed curve, consists of a shockwave and a contact discontinuity traveling to the right, and a left-traveling rarefaction wave. First thing to note is that the results obtained with astrix are almost indistinguishable from those obtained with the Roe solver. The second-order astrix scheme produces very minor overshoots near discontinuities, but the first-order $N$ scheme is monotone as expected. The first order schemes show more numerical diffusion, but in all cases the correct signal speeds are recovered. In the case of the Roe scheme, the sharpness of the shocks mildly depends on the choice of the flux limiter, while for astrix, they depend on the exact blending procedure. The $Bx$ scheme gives slightly sharper shocks compared to the $B$ scheme, at the expense of slightly stronger overshoots.

The $L_{1}$ density error is shown in Fig. 16 for different resolutions. The convergence obtained is less than first order for all algorithms. While this may be surprising at first, this is due to the presence of a contact discontinuity for which it is well-known that shock-capturing schemes show sublinear convergence rates111111In fact, the theoretical convergence rate is $1/2$ (Hedstrom, 1968). Figure 16 shows a slightly higher convergence rate because of the additional presence of both a shock and a rarefaction wave. (e.g. Hedstrom, 1968; Orszag & Jayne, 1974). It should be noted in particular that this behaviour is not limited to the Roe solver: the HLLC scheme shows similar sublinear convergence and gives results almost indistinguishable from the results obtained with the Roe solver. We note that this behaviour for second order methods is sensitive to the choice of flux limiter, and, in the case of astrix, one of the few cases where both the chosen mass matrix and the lumping strategy matter. Both a more compressive flux limiter than the minmod limiter used in Fig. 16 for the Roe solver (Orszag & Jayne, 1974) and a less diffusive mass matrix combined with selective lumping in the case of astrix give slightly better results than given in Fig. 16.

5.1.3 Blast wave

Similarity solutions to strong blast waves have been known for a long time (e.g. Bethe et al., 1947; Sedov, 1946; Taylor, 1950). While in the astrophysical community they are already well-known because they are useful models for supernova explosions, they can be used as test problems for numerical hydrodynamical codes as well. Here, we discuss a one dimensional variant in planar geometry, for which the method of obtaining a reference solution is discussed in appendix B.

The initial conditions consist of uniform density $\rho_{0}=1$, zero velocity and negligible pressure $p_{0}=10^{-6}$. To this we add a Gaussian pressure perturbation of the form

[TABLE]

with ratio of specific heats $\gamma=5/3$ and cell size $h$. The computational domain is $|x|<2$ and we evolve the system until $t=10$. Note that the total energy injected is constant for different resolutions and equals $E=1/50$, which means that for the lowest resolution we consider (200 grid cells) this amounts to an energy perturbation of unity within one grid cell.

A strong shock emerges that propagates away from $x=0$. The shock position and velocity as obtained with the $N$ scheme with 200 grid cells are displayed in Fig. 17. The position of the shock was taken to be the cell where the density passes through $\rho=2.5$ (halfway between the expected post- and pre-shock values) and decreases with $|x|$. Since we can not determine the position of the shock within one grid cell, the shock position shows a stair step pattern. This leads to a lot of noise when we differentiate with respect to time in order to get the shock velocity, hence these were smoothed with a window of $\Delta t=0.05$ in order to make the results readable. Results for the $B$ and $Bx$ schemes as well as the Roe solver are almost indistinguishable from the results with the $N$ scheme and are therefore not shown.

The $L_{1}$ density errors are shown in Fig. 18. All tested methods show linear convergence, and the errors are comparable in magnitude. As in the case of Sod’s shock tube, the results with the HLLC solver are indistinguishable form those obtained with the Roe solver. Results for the $Bx$ scheme (not shown) are very similar to those obtained with the $B$ scheme.

5.2 Two-dimensional tests

In two dimensions, we can make full use of an unstructured grid. An example mesh is shown in Fig. 2, for a computational domain that is periodic in both directions; as a consequence, the jagged outer right edge slots in the left edge, and similarly the bottom edge slots in the top edge. The number of vertices is $\sim 400$, so the resolution is equivalent so a structured mesh of $20\times 20$. Such a structured triangular mesh is shown in Fig. 1. Note that the structured mesh has clear symmetries, while the unstructured mesh is locally isotropic.

5.2.1 Isentropic vortex

As a first test case, we consider an isentropic stationary vortex located at $x=x_{c}$, $y=y_{c}$. The velocity profile is given by

[TABLE]

where $r=\sqrt{(x-x_{c})^{2}+(y-y_{c})^{2}}$ and the angular velocity is given by

[TABLE]

Equilibrium requires

[TABLE]

which, together with the isentropic assumption $p=K\rho^{\gamma}$, leads to

[TABLE]

Solving for $T=p/\rho$ we find

[TABLE]

where $p_{\infty}$ and $\rho_{\infty}$ are the pressure and density far away from the vortex. Pressure and density distributions follow from (85):

[TABLE]

We take $p_{\infty}=\rho_{\infty}=K=1$ and a vortex with $\beta=5$ placed at $x_{c}=5$, $y_{c}=5$ in a computational domain $0<(x,y)<10$. We let the vortex evolve until $t=10$, and measure the $L_{1}$ error in the density.

The results are shown in Fig. 19. It is clear that for this problem, the multidimensional upwind methods really come into their own, outperforming the Roe solver by more than an order of magnitude in the second order case. Note that since this is a stationary problem, the second-order time integration adds nothing and similar results can be obtained using the $Bx$ scheme with only first-order time integration. Comparing first and second order results, we see that part of this big difference in error between Roe and $Bx$ is due to the multidimensional characteristic decomposition, which is the only major difference between first order Roe and the $N$ scheme. The results obtained with a first order HLLC scheme are indistinguishable from those obtained with the first order Roe scheme.

All second order schemes show second order convergence, as expected. The second order Roe scheme does slightly better ($10-20\%$) than the second order HLLC scheme, which is probably linked to the way the flux limiter is applied when the characteristic waves are not aligned with with eigenvectors of the linearised Jacobian (LeVeque, 1997). The blended $B$ scheme significantly outperforms both HLLC and Roe, and the $Bx$ scheme does in fact more than an order of magnitude better than the second order Roe scheme. In this case, since there are no shocks present, the $Bx$ scheme yields exactly the same results as the $LDA$ scheme.

We now make the problem time-dependent by giving the whole domain an $x$ velocity boost of $v_{\mathrm{advect}}=1$. We enlarge the computational domain in the $x$ direction to $0<x<20$, and again run the simulation until $t=10$. This should place the vortex at $x=15$, with velocity and density structure unchanged from $t=0$. The resulting $L_{1}$ errors in the density are displayed in Fig. 20. Again, we see that astrix outperforms the Roe solver for both first and second order updates, but not by as much as in the previous case of a stationary vortex. The difference in error has gone down from roughly a factor of $3$ to a factor of $2$ for the first order schemes. Again, results obtained by HLLC and Roe are exactly the same at first order, while the Roe solver shows slightly better results at second order. The $Bx$ scheme, which again gives the same results as the $LDA$ scheme, outperforms Roe and HLLC by roughly a factor of $3$. The $B$ scheme shows a departure from second-order convergence towards high resolution. The cause of this remains to be investigated.

The reason that astrix performs so well on a stationary vortex lies in the linearity preserving nature of the $LDA$ scheme: if a stationary solution is linear over an element, it is recovered exactly by the integration scheme. Of course, the isentropic vortex is not a linear solution, but it can up to some precision be represented by a linear solution, which is the solution the $LDA$ scheme will evolve towards. Note that the solution only needs to be linear over every single triangle; it does not have to be linear globally. The difference between the $LDA$ solution and the true solution scales with $h^{2}$ (Abgrall, 2001). In other words: the scheme tries to find a solution where the residual $\phi^{T}=0$ everywhere, and since the distribution coefficients satisfy $\phi_{i}^{LDA}=\beta_{i}\phi^{T}=0$ no evolution takes place. It is not obvious, especially with a dimensionally split scheme, whether such a multidimensional numerically stationary state exists, unless it is a hydrostatic solution. Even in one spatial dimension in the presence of source terms one needs to be careful not to evolve away from stationary states (e.g. Eulderink & Mellema, 1995). What is happening in the case of the Roe solver, is that the solution continues to evolve to the only state that is numerically stationary for this scheme, which is the hydrostatic solution, in which case there is of course no vortex. Note that the $N$ scheme is not linearity-preserving and will therefore behave more like the Roe solver. Nevertheless, there is still an advantage coming from the multidimensional upwinding.

5.2.2 Noh problem

Here we discuss a demanding problem originally solved by Noh (1987), of which we consider the two-dimensional version. The initial conditions are uniform density $\rho=1$, and a velocity of magnitude 1 everywhere in the direction of the origin. If we take the initial pressure to be zero, the problem has an analytic solution consisting of a shock of formally infinite Mach number moving radially outward. Defining $r=\sqrt{x^{2}+y^{2}}$, and taking the ratio of specific heats $\gamma=5/3$, the solution reads:

[TABLE]

[TABLE]

[TABLE]

We take the computational domain to be $-1\leq x\leq 1$ and $-1\leq y\leq 1$, and as boundary conditions we impose the analytic solution. Since both astrix and the Roe solver are based on a characteristic decomposition, they can not handle zero pressure, and therefore, following Liska & Wendroff (2003), we set it to $10^{-6}$ initially. The solution is evolved up to $t=2$, at which point the shock is located at $r=2/3$.

Results are shown in Fig. 21. The top panels show results obtained with the first order Roe scheme, which clearly shows the carbuncle instability (Peery & Imlay, 1988; Quirk, 1994). This well-known instability that occurs for strong shocks that are aligned with the grid is usually associated with multidimensional flows, but a one-dimensional version does exist as well (e.g. Dumbser et al., 2004). In this case, the instability may be associated with the nonlinearity of the Hugoniot locus (Zaide, 2012). In multidimensional flows, it is thought to be associated with using one-dimensional fluxes (Stone et al., 2008), and the instability can be corrected efficiently by adding extra multidimensional dissipation (Quirk, 1994; Stone et al., 2008), a feature not implemented in the Roe scheme used here. Going to second order makes the results worse rather than better and therefore these are not shown. The HLLC solver suffers from similar problems.

Results obtained with astrix do not show the carbuncle instability. This is largely due to the use of an unstructured grid: there is no place where the shock can be said to be aligned with the grid. However, the multidimensional upwind nature of astrix appears to play a role as well. Results for structured grids do show some artifacts along the coordinate directions, but they remain confined to the postshock region as in the top left panel of Fig. 21, with no real carbuncles developing even at high resolution.

A numerical artifact that can be seen especially in the lower left panel is “wall heating" (Noh, 1987), resulting in an underdense region in the centre while the pressure is constant. This rise in temperature is seen in most Riemann solver codes (Liska & Wendroff, 2003; Stone et al., 2008), and even though a complete understanding of this phenomenon is still lacking, it may be due in part by an ambiguity in the position of a shock within a cell (Zaide, 2012). We get a minimum density in the central region of $14$, which makes it comparable to the schemes tested in Liska & Wendroff (2003) (better than PPM, slightly worse than the WENO scheme). While we only show results for the first order $N$ scheme, since the solution is dominated by a strong shock a blended second order scheme give exactly the same result as the $N$ scheme in the case where we take the blending coefficient as the maximum over all equations. Less diffusive blended schemes either fail on this problem, or need backup fluxes produced by the $N$ scheme so that again the results look exactly the same as in Fig. 21.

In Fig. 22 we show the $L_{1}$ density error for those algorithms that did not show any carbuncle instabilities. As a comparison to the astrix results we also show results obtained with the original Harten-Lax-van Leer solver (HLL, Harten et al., 1983), which is known to be carbuncle-free (Pandolfi & D’Ambrosio, 2001). Even the HLL method does show some artefacts around the coordinate axes; however, these have very limited impact on the total error. The errors are very comparable between the three schemes and show linear convergence at low resolution. At high resolution, a departure from linearity can be seen, which is probably due to the wall heating phenomenon.

5.2.3 Cylindrical blast wave

We now consider a blast wave in cylindrical geometry. The initial conditions consist of uniform density $\rho_{0}=1$, zero velocity and negligible pressure $p_{0}=10^{-6}$. To this we add a Gaussian pressure perturbation of the form

[TABLE]

with ratio of specific heats $\gamma=5/3$ and $h$ is a linear measure of the cell size. In the case of an unstructured grid this is taken to be the average over the whole domain of the square root of the cell volume. The total energy injected is equal to an order unity perturbation in a disc with radius $1/100$. A reference solution is calculated as outlined in appendix B. The computational domain is taken to be $|x|,|y|<0.2$ and the solution is evolved until $t=1$.

As in the planar case (see section 5.1.3), a strong shock emerges from the location of the energy perturbation. As in the case of the Noh problem (see section 5.2.2), both the Roe solver and the HLLC solver suffer from a carbuncle instability, while the astrix results are carbuncle-free because of the use of an unstructured grid.

In Fig. 23 we show the shock position (top panel) and shock velocity (bottom panel) as a function of time as obtained with the $N$ scheme at an equivalent resolution of $100\times 100$, together with the analytical expectation. As in the planar case, good agreement is found. Since the shock position is now obtained by an angular average, the result is less noisy and no smoothing is necessary to obtain the shock velocity, which again shows good agreement with analytical expectations.

In Fig. 24 we show the resulting $L_{1}$ error in the density for three algorithms that did not suffer from numerical instabilities. The $N$ scheme produces results that are almost indistinguishable from those obtained with the HLL scheme. Unlike for the Noh problem, the HLL solution looks perfectly smooth and with the largest part of the error originating from the strong shock it is hard to do any better on this problem. The $B$ scheme shows larger errors but linear convergence as do the $N$ and HLL schemes. Note that in this problem wall heating is not an issue, and the errors show a linear decrease for all resolutions considered.

5.2.4 Riemann problem

Here we study the two-dimensional Riemann problem originally introduced by Schulz-Rinne et al. (1993). The problem is defined on a square $0\leq x<\leq 1$ and $0\leq y\leq 1$, with initial conditions

[TABLE]

where ${\bf U}=(\rho,v_{x},v_{y},p)^{T}$ is the vector of primitive variables and

[TABLE]

The ratio of specific heats is taken to be $\gamma=1.4$. The solution consists of four shocks travelling along the wall with a complex interaction region with shear flow, susceptible to Kelvin-Helmholtz instabilities (e.g. San & Kara, 2014). The solution along the walls are moving single shocks, for which the speeds can be computed easily from the Rankine-Hugoniot jump conditions. We use a one-dimensional shock solution with the computed shock speeds as boundary conditions, and compute the solution until $t=0.8$,

Results are shown in Fig. 25. Each panel contains a result obtained with astrix for $y>x$ and a result obtained with the Roe solver for $y<x$. All schemes agree on the position of the shocks, as is to be expected for conservative schemes. The jet along the diagonal gets longer with increasing resolution and decreasing numerical dissipation. The left panels show results obtained with first-order schemes at resolutions $128\times 128$ and $1024\times 1024$. It is clear from the position of the jet that the $N$ scheme has lower dissipation than the first order Roe scheme. Moreover, the Roe scheme shows an instability at $x=0.8$ at the position of the shock in the lower left panel. The seeds of this instability can be seen in most panels as light horizontal streaks near $y=0.8$ and light vertical streaks near $x=0.8$. These streaks are present in almost all codes, regardless of Riemann solver (Liska & Wendroff, 2003), and are present even if the solution between for example ${\bf U}_{2}$ and ${\bf U}_{4}$ is computed in one dimension. Their source is purely numerical and is a consequence of a “startup error" (Jin & Liu, 1996; Zaide, 2012): the initial conditions are not a solution to the modified equation (including numerical dissipation) the code is solving. The adjustment to the numerical solution leads to small artifacts that in this case stay put as the post-shock velocity is close to zero. Remarkably, this leads to a numerical instability for the first-order Roe scheme. The second-order Roe schemes, possibly because it increases the stencil, are able to correct this instability very efficiently. The results obtained with the minmod limiter (middle panels in Fig. 25) still show the streaks at a resolution of $128\times 128$, albeit at lower amplitude compared to the first-order Roe scheme, but at a resolution of $1024\times 1024$ they have disappeared. Similar features show up with the HLLC solver, with some subtle differences. No instability is observed in the first order HLLC scheme, but on the other hand the streaks remain far more pronounced in the second order case.

A second interesting feature is the appearance of Kelvin-Helmholtz vortices in the shear flow past the shocks. Except for the first order Roe scheme, which becomes dominated by the numerical instability due to the start-up error, these vortices show up in all schemes at sufficient resolution. It has been argued that these are physical rather than numerical (e.g. San & Kara, 2014), and the fact that they show up in all schemes lends support to this conclusion. However, this suggests that for higher and higher resolutions, the schemes will not converge to a well-defined solution, making the problem less suited as a numerical test (a similar problem haunts the Kelvin-Helmholtz instability, see section 5.2.5). This is especially apparent from the lower right panel of Fig. 25, where the solutions look very noisy, especially in the case of the Roe solver with superbee flux limiter. The same holds for results obtained with the HLLC solver.

5.2.5 Kelvin-Helmholtz instability

The Kelvin-Helmholtz instability (KHI) is a classical hydrodynamical instability (e.g. Chandrasekhar, 1961) that operates in systems where a velocity shear is present. It has attracted significant attention as a test problem ever since Agertz et al. (2007) showed that SPH has difficulties resolving the KHI. In addition, Springel (2010) claimed that grid-based methods violated Galilean invariance using the KHI with an added bulk velocity. This makes it an interesting test problem to consider.

The simplest two-dimensional setup demonstrating the KHI consists of two fluids shearing past each other on a square $0<(x,y)\leq 1$, periodic in both $x$ and $y$, so that the unperturbed velocity ${\bf v}=(u,v)^{T}=(u_{0}(y),0)^{T}$, with

[TABLE]

where $U_{0}$ is a constant. While strictly not necessary for the KHI, we take a density profile of similar profile in order to highlight the mixing properties of the schemes:

[TABLE]

where $\rho_{\mathrm{H}}$ and $\rho_{\mathrm{L}}$ are a constant high and low density, respectively. Taking a constant pressure $p_{0}$, we have an equilibrium solution that is however unstable to the KHI.

Analysis of the incompressible problem (e.g. Chandrasekhar, 1961) reveals that perturbations grow $\propto\exp(st)$ with growth rate

[TABLE]

where $k_{x}$ is the horizontal wave number. Interestingly, all wavenumbers are unstable, and the highest wave numbers have the largest growth rates. On an unstructured grid, initial perturbations are present on the grid scale, and because they grow fastest they come to dominate the solution. Fortunately, in the case of the KHI there is a simple way to regularize the problem by considering smooth profiles of velocity and density in stead of step functions. As shown in Chandrasekhar (1961), if we take the transition between $U_{0}$ and $-U_{0}$ to be linear in $y$ over a distance $d$, unstable wave numbers must satisfy

[TABLE]

where $\kappa_{0}\approx 1.27846$ is the solution of $\kappa=1+\exp(-\kappa)$. Roughly speaking, wavelengths smaller than $d$ are stable. Smooth initial conditions therefore lead to a better-posed problem, at least at early times (McNally et al., 2012), but in order to make the analysis more quantitative we need linear growth rates to compare against, which are computed using the method described in appendix C.

In order to obtain a smooth profile, we first define a function

[TABLE]

so that

[TABLE]

describes a smooth transition from [math] to $1$ over an interval $0\leq t\leq 1$. A smooth version of the discontinuous velocity profile (97) is then given by

[TABLE]

where the discontinuity has been spread over a distance $d$. Similarly for the initial density profile:

[TABLE]

Growth rates were computed for $U_{0}=1/2$, $\rho_{\mathrm{L}}=1$, $\rho_{\mathrm{H}}=2$, $p_{0}=2.5$ and $\gamma=1.4$ for various values of $d$ and $k_{x}$. The results are shown in Fig. 26. It is clear that, as expected from the incompressible analysis, for finite $d$ there is a maximum wave number for which the flow is unstable. The effect of compressibility is to reduce the growth rate, but with $p_{0}=2.5$ the velocities are all subsonic and for $d\rightarrow 0$ the growth rates come close to the incompressible result.

We now choose a relatively large value of $d=0.25$, which makes for only one unstable mode with $k_{x}=2\pi$. The linear growth rate was computed to be $s=1.551$. Following McNally et al. (2012), we monitor the kinetic energy using the vertical velocity only, but integrated over the entire computational domain:

[TABLE]

This quantity is expected to grow at a rate $2s$. As initial conditions, on top of the background flow as specified above, we put in the eigenvector belonging to the growth rate $s$ with a velocity amplitude of $\sim 10^{-3}$.

Results are shown in Fig. 27 for first order schemes at a resolution of $128\times 128$ and second order schemes at a resolution of $32\times 32$. These resolutions were chosen as to highlight differences between the various schemes; at higher resolutions they all converge to the linear result, which is shown by the thick solid line. The linear phase lasts until roughly $t=3.5$, after which the instability saturates. As an average growth rate between $t=0$ and $t=2$, we measure $s=0.6$ for the first order Roe scheme and $s=0.9$ for the first order $N$ scheme. This trend of the $N$ scheme showing higher growth rates continues towards higher resolution until both first order schemes approach the theoretical growth rates. In order for the first order Roe scheme to produce similar growth rates, we need roughly a factor of 2 higher resolution compared to the $N$ scheme. A similar story, but less dramatic, holds for the second order schemes: at a resolution of $32\times 32$, we measure $s=1.4$ for the $B$min scheme and $s=1.3$ for the second order Roe scheme. For resolutions $64\times 64$ and higher, the results of the two second order schemes become indistinguishable, at least in the linear phase.

The KHI saturates by perturbations rolling up into large vortices. The onset of this phase at $t=3.5$ is depicted in Fig. 28. While at this early non-linear stage, all schemes still give fairly similar results, two things are worth pointing out. First of all, it is clear from the top row, which shows results for the first order schemes, that the $N$ scheme outperforms the first order Roe scheme at a resolution of $128\times 128$. This was clear from the growth rates in Fig. 27, but the top left panels really drive this point home. The first order Roe scheme at $256\times 256$ gives comparable results to the $N$ scheme at resolution $128\times 128$. All second order schemes show very similar amplitudes at this stage at the resolutions considered. The second point is the small artifact around $(x,y)=(0.6,0.7)$ seen in the first-order Roe scheme as well as the astrix results. This kink occurs at the density interface where the $y$-velocity is close to zero, which means the eigenvalue corresponding to a contact discontinuity in the $y$-direction is close to zero, which leads to low numerical dissipation. In this particular case, the numerical dissipation is in fact too low for this kink to disappear. Numerical experiments show that this artifact disappears when adding a constant $y$-velocity to the whole domain.

It is a well-known problem of linearised Riemann solvers like the Roe solver that eigenvalues close to zero can lead to numerical artifacts. Most notably, the Roe solver can not deal properly with transsonic rarefaction waves, where the eigenvalue corresponding to the slow sonic wave passes through zero, which requires some form of entropy fix (Harten & Hyman, 1983). While entropy fixes have been designed for multidimensional upwind methods (Sermeus & H., 2005), these have not been implemented in astrix at present and in any case they would not remove the kinks seen in Fig. 28. In fact, results obtained with the HLLC solver, which in general does not require an entropy fix to deal with transsonic rarefactions, show exactly the same feature. Fortunately, its effect is much less severe than errors in transsonic rarefactions, as it appears to be limited to a few wiggles in the density profile and does not grow with time.

The KHI for the same setup but at at later time $t=5$ is shown in Fig. 29. From the first-order schemes in the top row, it is clear that indeed the kinks visible in Fig. 28 do not grow. The most obvious differences between the various panels in Fig. 29 are the secondary instabilities seen in the second-order astrix results as well as in the results obtained with the Roe solver using the superbee flux limiter. The same holds for results obtained with the HLLC solver. There has been a lot of discussion in the literature on these instabilities (McNally et al., 2012; Lecoanet et al., 2016). While it is clear that there is a relation between the occurrence of these secondary instabilities and numerical diffusion (for the Roe and HLLC solver, they only show up when using the least-diffusive flux limiter), if there is no convergence with resolution to a well-defined solution this relation is meaningless. And this appears to be the case with the KHI without any physical dissipation (i.e. viscosity), even with smooth initial conditions as applied here. While the Roe solver with the minmod flux limiter does not show small-scale structure at the resolutions shown in Fig. 29, they do show up at later times at higher resolution. Therefore, while smoothing the initial conditions leads to well-posed problem in the linear phase, where growth rates can be compared to results from linear calculations, at later times physical dissipation is needed to regularise the solution (Lecoanet et al., 2016).

Interestingly, there is a variant of the KHI that appears to be well-posed even at later times. The KHI as discussed until now has velocities that are subsonic: $U_{0}^{2}<\gamma p_{0}/\rho_{\mathrm{L}}$. It is known that for supersonic velocities, the KHI changes character drastically (Karimi & Girimaji, 2016), with no roll-up into large vortices. In order to study the supersonic KHI, we take the same initial conditions but lower the initial pressure to $p_{0}=0.1$. While again linear growth rates can be computed, it is more difficult to focus on a single growing mode, since even at $d=0.25$ growth rates are positive up to at least $k_{x}=20\pi$, although the growth rates decrease rapidly with $k_{x}$.

The non-linear phase for the supersonic KHI at $t=8.5$ is shown in Fig. 30. The top row shows results for the first-order schemes, and while the $N$-scheme shows stronger growth for $k_{x}=2\pi$, it is also clear that other modes are present, unlike in the first order Roe scheme. This is again due to the unstructured grid: all modes are present initially, and they all grow. Only if the scheme gets the relative growth rates correct will the $k_{x}=2\pi$ mode stand out, which happens at much higher resolution than depicted in Fig. 30 in both the first Roe scheme and the $N$ scheme.

The results improve drastically for the second-order schemes, for which all results look remarkably similar. The few small-scale features that can be seen are agreed upon by all schemes, suggesting that the solution to the supersonic KHI is well-behaved even in the absence of physical dissipation. The maximum density increases slightly when going from $128\times 128$ to $256\times 256$ for all schemes, but the results appear to be almost converged at this resolution. Note that the astrix results show a higher maximum density compared to the Roe results. Results obtained with the HLLC solver are almost indistinguishable from the Roe results. It is interesting to note that the $LDA$ scheme performs equally well on this problem despite sharp gradients present in the solution. The same can not be said for the Roe and HLLC schemes if we force a second-order update everywhere (not shown in Fig. 30), in which case spurious oscillations appear in the solution.

6 Discussion

In this paper, we have presented an implementation of a residual distribution method in an astrophysical fluid dynamics package astrix. A key difference between astrix and other grid based methods is its inherently multidimensional nature. While multidimensional integration schemes, usually in the form suggested by Colella (1990), have become a standard part of codes working on structured meshes (e.g. Fromang et al., 2006; Mignone et al., 2007; Stone et al., 2008), for methods employing unstructured grids these are rarely mentioned. This is probably partly due to the inherent complexity of unstructured meshes, but also due to the fact that methods based on one-dimensional flux estimates perform very well even for a naive implementation of multidimensional integration. For example, the Roe solver method as tested against astrix in this paper performs equally well when using formally second order Strang splitting (Strang, 1968) as when just alternating the direction of integration. However, especially in the case of the isentropic vortex problem, there are clear advantages for a truly multidimensional update as employed in astrix. In this section, we discuss the limitations of the current version and pathways to improvement.

From the test problems presented in section 5 it is clear that astrix performs at least as well as the dimensionally split Roe solver in all cases, and performs significantly better at some, notably the isentropic vortex and the Noh problem. This performance comes at a price of increased complexity of the method and increased computational time. In order to assess this quantatively we ran a simple speed test comparing astrix to various publicly available codes. Besides rodeo, which was used to compare to astrix in terms of accuracy, we took the structured grid code pluto (Mignone et al., 2007)121212See http://plutocode.ph.unito.it/, the meshless code gizmo (Hopkins, 2014, 2015)131313See http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html and the moving mesh code rich (Yalinewich et al., 2015)141414see https://github.com/bolverk/huji-rich. We chose a particularly simple test problem for the speed comparison: constant density, constant pressure and zero velocity on the unit square. The main reason for this choice it that it makes the Lagrangian and Eulerian approaches equal, allowing for a fairer comparison of the different codes. For example, the public version of gizmo only allows for Lagrangian meshless integration, either Meshless Finite Mass (MFM) or Meshless Finite Volume (MFV), see Hopkins (2015) for details. In this case, a hydrostatic calculation should eliminate overhead for example from finding new neighbours, or, in the case of rich, determining the new Voronoi tesselation. A second reason for this simple problem is that all codes should be able to do this problem without any special effort. For example, rich can only do the Noh problem when the "cold flow" option is activated, which makes the code three times slower.

The results of the speed test are shown in Table 1. All codes except astrix with the $N$ scheme were run in order to achieve second-order accuracy. The third column lists the CPU time per cell per time step, as measured on an Intel Core i7 2.2 GHz CPU with 8 GB DDR3 memory. These timings represent the speed of the algorithms themselves, and ignore for example the fact that astrix can run with larger time steps than a dimensionally split scheme such as rodeo. They also ignore the initial construction of the mesh. While this takes virtually no time in the case of the Roe solver, for astrix it is a more significant effort. As a rule of thumb, the creation of the mesh takes about as much time as taking ten time steps with second order $LDA$. In all of the test problems performed here, this means mesh creation takes roughly $1\%$ of the total time spent. While this is clearly not an issue for a static mesh, if the mesh is to be updated frequently this can quickly become a bottleneck.

The results in Table 1 show that the CPU version of astrix is roughly $2$-$2.5$ times slower per cell per time step than structured grid codes such as rodeo and pluto. This is due to the increased complexity of a multidimensional upwind update. On the other hand, astrix compares favourably to methods that do not employ a structured grid: it is roughly $5$ times faster per cell per time step than gizmo, and roughly $10$ times faster per cell per time step than rich.

We have focused on static uniform meshes, to simplify the error analysis and to make a fair comparison to the Roe solver. One of the beautiful features of unstructured mesh solvers is that no extra effort is required to run on a mesh with varying resolution, unlike for example in traditional AMR, where interpolation at resolution jumps is necessary. As explained in section 2.3.1, it is straightforward to generate a mesh with spatially varying resolution, which therefore makes running on non-uniform meshes a trivial matter. However, this is only useful if it is known in advance the location where high resolution is going to be necessary, and if that location is fixed in time. While such problems do exist (these can be tackled in traditional methods using static mesh refinement or nested grids), in general it is necessary to be able to dynamically update the mesh. This can be done by adding and removing vertices, but since these operations change the triangulation and therefore the connectivity of the mesh these tasks are not trivial and will be the subject of future work. Here we also note that a moving mesh approach is also possible (Michler et al., 2003).

The current implementation works in two spatial dimensions only. This is the simplest setup in which to demonstrate the power of multidimensional upwind methods, but the approach of section 3 can be generalised in a straightforward way to three dimensions, at the expense of adding a row and a column to all matrices, and the resulting increase in complexity of the matrix element computations. Most of the work in going to three dimensions will go into the generation of the mesh. Some of the nice properties of Delaunay triangulations do not generalise to three dimensions, and in particular some low-quality tetrahedra can survive Delaunay refinement (see e.g. Shewchuk, 2002). Fortunately there are ways to remove these (Cheng et al., 2000).

We have considered Cartesian coordinates only. For many astrophysical problems, spherical or cylindrical coordinates are better choices to represent the flow. Astrophysical discs, for example, in traditional Eulerian methods are best described in cylindrical coordinates, so that the flow is mostly aligned with the mesh. A Cartesian frame leads to excessive diffusion as angular momentum is not conserved, which is particularly a problem in the field of disc-planet interactions (de Val-Borro et al., 2006). Even when using cylindrical coordinates, care must be taken to conserve angular momentum to machine precision (Kley, 1998). The test problem of the isentropic vortex clearly shows that angular momentum is much better conserved using a multidimensional upwind method compared to the dimensionally split Roe solver on a Cartesian grid. However, when exact angular momentum conservation is required, this is possible in two steps. First, there is no reason why a Delaunay triangulation can not be used to tesselate a region in cylindrical coordinates $(r,\varphi)$: just replace $x$ with $r$ and $y$ with $\varphi$ in Fig. 2 (and remove the periodicity in $x$). Second, solve the Euler equations in cylindrical coordinates. The main problem here is that while a conservative linearisation in general coordinates does exist (Eulderink & Mellema, 1995), this does not conserve angular momentum to machine precision (Paardekooper & Mellema, 2006). Fortunately, it is possible to formulate the residual distribution schemes presented here even when a conservative linearisation is not available (Csík et al., 2002). This formulation can bring the power of multidimensional upwind methods and unstructured grids to simulations of astrophysical discs and this will be considered in a future work.

While we have considered hydrodynamics only, it is possible to add additional physics to the residual distribution methods presented here. For example, Csík et al. (2002) use their residual distribution formulation without a conservative linearisation to study the case of ideal magnetohydrodynamics, which included a source term to clean the divergence of the magnetic field. Constrained transport algorithms on unstructured meshes were recently presented in Mocz et al. (2014); Mocz et al. (2016). A general diffusion solver in residual distribution form was presented in Nishikawa (2007), which can for example be used to include self-gravity or radiative transport in the diffusion limit. A more accurate radiative transfer method on unstructured grid can for example be found in Ritzerveld (2007) and Paardekooper et al. (2010). For examples of (general) relativistic solvers on unstructured grids see Anninos et al. (2005); Duffell & MacFadyen (2011). Therefore, while unstructured meshes call for different implementations, they impose no limit of the amount of physics that can be included in simulations.

We have only considered schemes that are at most second order accurate in space and time. While higher-order explicit residual distribution schemes are not available as yet, it is expected that by considering higher order elements the extension of the current method to higher order would be straightforward (Ricchiuto & Abgrall (2010), see also Cohen et al. (2001); Jund & Salmon (2007)).

7 Conclusions

We have presented a GPU implementation of a multidimensional upwind, or residual distribution, method acting on an unstructured Delaunay mesh in two spatial dimensions. While more expensive than than traditional grid-based methods, it has clear advantages when dealing with multidimensional flows over methods that use one-dimensional flux estimates as building blocks for a multidimensional integration. This is exemplified in several test problems, most notably the problem of a stationary isentropic vortex. Our GPU implementation shows speedups of $\sim 100$ for mesh generation and $\sim 250$ for the hydrodynamics, with the most expensive kernels coming close to the theoretical performance limit of our Tesla K20m GPU.

Acknowledgements

SJP is supported by a Royal Society University Research Fellowship.

Appendix A Multidimensional upwind matrices

Consider the Euler equations in two spatial dimensions in the form:

[TABLE]

then the matrix $\mathcal{K}=\mathcal{A}n_{x}+\mathcal{B}n_{y}$ is given by

[TABLE]

with $w=un_{x}+vn_{y}$, $\alpha=\gamma_{1}(u^{2}+v^{2})/2$ and $\gamma_{1}=\gamma-1$, and has eigenvalues $\lambda_{1}=w+c$, $\lambda_{2}=w-c$, $\lambda_{3}=\lambda_{4}=w$, and right eigenvectors

[TABLE]

Putting the eigenvectors as columns in a matrix R:

[TABLE]

then

[TABLE]

and $\mathcal{K}$ can be diagonalised with $\mathcal{K}=\mathcal{R}^{-1}\Lambda\mathcal{R}$, where $\Lambda$ is a diagonal matrix with the eigenvalues on the diagonal. If we put generic values $l_{1}$, $l_{2}$, $l_{3}=l_{4}$ in $\Lambda$, we get

[TABLE]

and matrix entries

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with $\alpha_{c}=\alpha/c$, $u_{c}=u/c$, $v_{c}=v/c$, $\gamma_{2}=\gamma-2$, $l_{123}=(l_{1}+l_{2}-2l_{3})/2$ and $l_{12}=(l_{1}-l_{2})/2$. If we set $l_{k}=\lambda_{k}$ for $k=1,2,3,4$, we recover $\mathcal{K}$, while if we set $l_{k}=\lambda_{k}^{\pm}$ we get $\mathcal{K}^{\pm}$.

Appendix B Blast wave solutions

Analytical blast wave solutions emerged shortly after the Second World War (Bethe et al., 1947; Sedov, 1946; Taylor, 1950) and make useful test problems for numerical gas dynamics codes. Most well known in the astrophysical community is the spherical blast wave, which is a useful model for a supernova explosion. However, in this work we are concerned with at most two spatial dimensions, and therefore we need lower-dimensional analogs of the spherical blast wave. Fortunately, these exist (e.g. Sedov, 1959; Kamm & Timmes, 2007), and below we briefly describe the procedure of obtaining reference blast wave solutions to compare to hydrodynamical calculations.

Consider a point explosion in spherical geometry, a line explosion in cylindrical geometry and a plane explosion in Cartesian geometry. Assume a strong shock, so that the preshock pressure $P_{1}$ plays no role. The only parameters in the problem are the preshock density $\rho_{1}$ and a measure of the total input energy $E$, which in the cylindrical case is an energy per unit length, and in the Cartesian case an energy per unit area. This leads to the definition of a similarity variable

[TABLE]

where $r$ denotes the distance to the initial energy release, and $a=0$ for Cartesian, $a=1$ for cylindrical, and $a=2$ for spherical coordinates.

Jump conditions for a strong shock give

[TABLE]

where $v_{s}$ is the velocity of the shock and $\gamma$ is the ratio of specific heats. In the region behind the shock, use the ansatz

[TABLE]

Note that $\hat{\rho}(\eta_{s})=\hat{P}(\eta_{s})=\hat{v}(\eta_{s})=1$, where $\eta_{s}$ denotes the position of the shock.

The continuity equation reads

[TABLE]

where again $a=0$ for Cartesian, $a=1$ for cylindrical, and $a=2$ for spherical coordinates. Transforming to $(\eta,t)$ gives:

[TABLE]

The momentum equation reads

[TABLE]

Transforming to $(\eta,t)$ gives:

[TABLE]

The energy equation reads

[TABLE]

Transforming to $(\eta,t)$ gives:

[TABLE]

The three equations for $\hat{\rho}$, $\hat{P}$ and $\hat{v}$ need to be solved subject to the constraint that the total energy of the solution equals the initial energy input:

[TABLE]

where $b=4\pi$ for spherical coordinates, $2\pi$ in cylindrical coordinates, and $1/2$ in Cartesian coordinates. In dimensionless form:

[TABLE]

In order to facilitate numerical integration, we rewrite the equations of mass and momentum conservation as:

[TABLE]

Use these two equations in the energy equation to obtain

[TABLE]

The resulting expression for $d_{\eta}\hat{v}$ can be used in equations (146) and (147) so that equations for $d_{\eta}\hat{\rho}$, $d_{\eta}\hat{P}$ and $d_{\eta}\hat{v}$ in terms of $\hat{\rho}$, $\hat{P}$ and $\hat{v}$ result. These can be solved numerically using standard techniques. First, one needs to guess a value of $\eta_{s}$, integrate the equations from $\eta=\eta_{s}$ to $\eta=0$, and check energy conservation (145). This process is repeated for different $\eta_{s}$ until the energy constraint is met. The resulting solutions for $a=0,1,2$ are shown in Fig. 31.

Appendix C Compressible Kelvin-Helmholtz growth rates

Consider a 2D domain, periodic in $x$ and $y$ with periods $L_{x}$ and $L_{y}$ and no gravity:

[TABLE]

Give the fluid a smooth profile in $x$-velocity and density that depends only on $y$: $\rho_{0}=\rho_{0}(y)$, $u_{0}=u_{0}(y)$, $v_{0}=0$, and constant pressure $p_{0}$. The linear perturbation equations read:

[TABLE]

together with adiabatic pressure perturbations $p_{1}=c_{0}^{2}\rho_{1}$ with $c_{0}^{2}=\gamma p_{0}/\rho_{0}$ the square of the unperturbed sound speed. Take perturbations $\propto\exp(ik_{x}x-i\omega t)$:

[TABLE]

Write in terms of momenta $a_{1}=\rho_{0}u_{1}$ and $b_{1}=\rho_{0}v_{1}$:

[TABLE]

Write all quantities, both perturbed and unperturbed, as a Fourier series, e.g.:

[TABLE]

The equations for component $n$ read:

[TABLE]

Fourier components of products are of course convolutions, so that for example

[TABLE]

where ${\bf d}_{1}$ is a vector of length $N$ with entries $\rho_{1m}$ and ${\bf u}_{n}$ is a vector with entries $u_{0}(n-m)$, with $0\leq m<N$. We can therefore write all $N$ equations for the components in the form of matrices:

[TABLE]

where matrix ${\bf U}$ has entries

[TABLE]

and matrix ${\bf U^{\prime}}$ has entries:

[TABLE]

and matrix ${\bf N}$ has entries:

[TABLE]

and matrix ${\bf C}$ has entries:

[TABLE]

Combining ${\bf d}_{1}$, ${\bf a}_{1}$ and ${\bf b}_{1}$ into a single vector ${\bf e}$ of length $3N$:

[TABLE]

It therefore remains to find the eigenvalues and eigenvectors of the $3N$ by $3N$ matrix

[TABLE]

If the maximum of the imaginary parts of the eigenvalues is larger than zero, the wavenumber $k_{x}=2\pi m/L_{x}$ for some integer $m$ is unstable, with growth rate equal to this imaginary part.