The Spectral Difference Scheme

The spectral difference method has quite recently been proposed by Liu et al. [8] and further developed by Wang et al. [13] and by the present authors [9]. Consider the conservation law

$$\dd{u}{t} + \nabla \cdot F = 0 \ , \qquad (x,t) \in \Omega \times [0,T] \ ,$$ (1)

subject to suitable initial and boundary conditions, where $\Omega\subset \mathbb{R}^d$, and $u=(u_1,\dots,u_p)$, and the divergence operation is to be taken componentwise for $p>1$. The Spectral Difference method uses a pseudo-spectral collocation-based reconstruction for both the dependent variables $u(x)$ and the flux function $F(u)$ inside a mesh element, $T_i$, say. Assume a triangulation of $\mathbb{R}^d$ consisting of simplexes. The reconstruction for the dependent variables can be written

$$u_{i}(x) = \mathcal I^m(u)(x) = \sum_{j=0}^{N(m)} L_j(x) u_{ij} \ , \qquad x \in T_i \ ,$$ (2)

where $u_{ij}= u(x_{ij})$, and $x_{ij}$ is the $j^{th}$ solution collocation node in the $i^{th}$ mesh element (here and in the following the double index notation refers to a node inside a cell). The interpolation operator $\mathcal I^m$ denotes a collocation using polynomials of total degree $m$. The $L_j(x)$ are the cardinal basis functions for the chosen set of collocation nodes $x_{j}$, where $j=0,\dots,N(m)$, and $N(m)$ is determined by the chosen order of accuracy. The reconstruction of the flux function in $T_i$ reads

$$F_i(u(x)) = \mathcal I^{m+1}(F(u))(x) = \sum_{k=0}^{N(m+1)} M_k(x) F_{ik} \ , \qquad x \in T_i \ ,$$ (3)

where the $M_k$ are the cardinal basis functions corresponding to the collocation nodes $x_{ik}\ , k=0\dots N(m+1)$, and $F_{ik}=F(u(x_{ik}))$. If the solution is reconstructed to order $n$, the flux nodes are interpolated to order $n+1$, because of the differentiation operation in eq. (1). If the numerical flux is defined as the union of all local interpolations, it will be discontinuous at element boundaries. Instead we define the numerical flux function $F_h$ on the triangle $T_i$ as

$$F_h = \left\{\begin{array}{ll} F_i & , \quad x \in T_i\\ G (u_i(x),u^{e,1},\dots,u^{e,k}\dots )& , \quad x \in \partial T_i \end{array}\right. \ ,$$ (4)

where the $u^{e,1},\dots,u^{e,k}\dots$ are ``external'' solutions on triangles $T_j$ with $j(k)\neq i$ such that $x \in \partial T_j$, and the value is given by $u^{e,k} = \lim_{y\to x} u_{j}(y)$. For nodes on edges of triangles there is only one external solution $u^e$. It is necessary for discrete conservation that the normal flux component be continuous across the edge, which suggests one uses numerical flux functions, standard in finite-volume formulations, such that the normal flux component $F_n = F_i\cdot n$, for the node $x_k$, say, is replaced by the numerical flux $h(u_i(x_k),u^{e},n)$, where $n$ is the edge normal, and $h$ is the numerical flux function approximating $F_n$. For flux nodes on corners of elements, one may compute the flux using the numerical flux functions associated with both incident edges, see [13]. For the numerical flux we have used central average with both CUSP construction of artificial diffusion [7] and simple scalar dissipation, as well as Roe's approximate Riemann solver with good results.

Any combination of collocation nodes may be used, provided that the nodes for $u$ support a quadrature of the order of the interpolation $n$, and the restriction of the flux nodes to the boundaries supports a $d-1$-dimensional quadrature of order $n+1$. This ensures discrete conservation in the sense that is satisfied exactly for the solution and reconstructed flux function [9]. For the solution nodes one can choose Gauss quadrature points. Hesthaven proposed nodes based on the solution of an electrostatics problem for simplexes [6] , which support both a volume and a surface integration to the required degree of accuracy. These nodes can be used for both flux and solution collocation. The baseline scheme is now readily defined in ODE form as $u_h = L_h(u_h)$, where the degrees of freedom are given by the values of the solution at the collocation nodes, $u_{ij}= u_h(x_{ij})$, where again $x_{ij}$ is the $j^{th}$ solution collocation node in the $i^{th}$ mesh element, and the right-hand side solution operator is given by the exact differentiation of the reconstructed flux function:

$$\frac{du_h(x_{ij})}{dt} = \left( \nabla \cdot F_h\right)(x_{ij} )$$ (5)