Spectral difference method for unstructured grids I: Basic formulation

Abstract

A new, high-order, conservative, and efficient method for conservation laws on unstructured grids is developed. It combines the best features of structured and unstructured grid methods to attain computational efficiency and geometric flexibility; it utilizes the concept of discontinuous and high-order local representations to achieve conservation and high accuracy; and it is based on the finite-difference formulation for simplicity. Universal reconstructions are obtained by distributing unknown and flux points in a geometrically similar manner for all unstructured cells. Placements of these points with various orders of accuracy are given for the triangular elements. Accuracy studies of the method are carried out with the two-dimensional linear wave equation and Burgers’ equation, and each order of accuracy is verified numerically. Numerical solutions of plane electromagnetic waves incident on perfectly conducting circular cylinders are presented and compared with the exact solutions to demonstrate the capability of the method. Excellent agreement has been found. The method is much simpler than the discontinuous Galerkin and spectral volume methods for unstructured grids.

Introduction

A current problem of great interest is to develop a numerical method for conservation laws with the following properties: that it be (a) conservative, (b) high-order accurate, (c) geometrically flexible, (d) computationally efficient, and (e) simply formulated. Simplicity and computational efficiency can be achieved using structured grids. The earliest and most widely used method is the finite-difference (FD) method employing a body-fitted curvilinear coordinate system [3], [37], with the equations written in strong conservation law form [39]. The spatial differencing is essentially one-dimensional, and carried out along coordinate directions. Thus a large number of data points are ignored in high-order stencils. Near boundaries, the stencil has to be modified with one-sided formulas. Since numerical grid generators are mostly only second-order accurate, the numerical differencing of grid point coordinates in evaluating metric terms can severely degrade the accuracy of the solution if the grid is not sufficiently smooth. The unknowns are solution values at grid points. Therefore the true integral conservation laws can only be satisfied to second-order accuracy. When a single structured grid is not feasible for very complex geometries, multi-block patched or overlapping grids are employed [4]. At interface boundaries between patches, or in overlapped regions, the high accuracy is generally degraded or sophisticated interface algorithms are needed.

In order to satisfy the integral conservation laws, finite-volume (FV) methods were developed, e.g. [29], [19]. The unknowns are now cell averages over quadrilaterals (2D) or hexahedra (3D). A high order reconstruction in terms of neighboring unknowns is used to calculate flux integrals over cell boundaries, using Riemann solvers and appropriate limiters. In practice, the conventional finite-volume method for structured grids does not overcome the limitations of the finite-difference method. The reconstruction is still one-dimensional along coordinate directions. While geometric quantities such as surface area vectors or cell volumes can be precisely calculated, flux integrals and volume integrals are usually evaluated by one-point quadratures, and are only second-order accurate. Both methods suffer significant loss of accuracy for very unsmooth, highly curved grids.

In order to achieve geometric flexibility along with high accuracy, we normally use an unstructured grid consisting of triangles in 2D and tetrahedra in 3D. The most commonly used conservative unstructured method is the finite-volume (UFV) method, applied to the integral form of the conservation law with cell averages of the conservative variables as the unknowns [2], [13], [26]. A polynomial reconstruction of any desired order of accuracy for each cell is obtained in terms of unknowns at neighboring cells. The flux integral for each face is evaluated using the reconstructed solutions from the two cells sharing the face and an approximate Riemann solver. A quadrature approximation is employed for non-linear flux functions. Thus, conservation is satisfied locally for each cell. However, due to the unstructured nature of the grid, it is difficult to obtain a non-singular stencil [13]. This necessitates a least-squares inversion in general. For very high order of accuracy, the number of cells, and thus the number of operations to carry out the numerical procedure, can become very large in three dimensions. This would hamper the efficiency of the method. Furthermore, since each unknown employs a different stencil, one must repeat the least-squares inversion for every cell at each time step, or must store the inversion coefficients. In a high-order, three-dimensional computation, the former would involve impractically large CPU times, while for the latter the memory requirement becomes prohibitive. In addition, the data from neighboring cells required for the computation can be far apart in memory. This further degrades the efficiency of the method due to data gathering and scattering. As a result of these deficiencies, the UFV method is limited to second-order accuracy in most applications.

An alternate method for unstructured grids is the finite-difference (UFD) method, applied to the differential form of the conservation law with values of the conservative variables at grid nodes as the unknowns [8], [1], [36]. Actually, one only needs an arbitrary set of nodal points, without the connections that define a grid. Here one employs a local polynomial reconstruction of the fluxes in terms of neighboring values determined by the unknowns. The method is simpler than the UFV method, since one only needs to differentiate the reconstructed solution. However, the UFD method has a major disadvantage of not being locally or globally conservative.

Finite-element (FE) methods [18] have long been used for unstructured grids because of their geometric flexibility. As in the unstructured FV method, the global domain is subdivided into triangular or tetrahedral cells, called elements. A major difference between the FE and FV or FD methods is that in the former we employ reconstruction data from within the element, while in the latter the reconstruction data comes from outside the cell. In the FE formulation, the unknowns are nodal values at nodes which are placed at geometrically similar points in each element. As a result, the local reconstructions become universal for all elements in terms of the same set of cardinal basis functions or shape functions. By assigning a single value to the nodes on element boundaries, a global reconstruction is now piecewise continuous. The conservation equations are then satisfied in a weak form via the method weighted of residuals (MWR), by multiplying the equations with the requisite number of test functions, integrating over the global domain, and using integration by parts. Usually the Galerkin approach is used, in which the test functions are the same as the basis functions. This results in a set of coupled equations of all unknowns. Their solution involves a very large, sparse matrix, whose entries depend on the element geometries. For non-linear equations, quadrature approximations are necessary to evaluate the matrix entries. While the integral conservation law is satisfied for the global domain, it is not satisfied for each element.

In practice, all the above methods normally employ relatively low-order approximations in their formulations. We now turn to a class of methods called spectral methods which have the properties of very high accuracy and spectral (or exponential) convergence [14], [5]. In traditional spectral methods, the unknown variable is expressed as a truncated series expansion in terms of some basis functions (trial functions) and solved using the MWR. The trial functions are infinitely differentiable global functions, and the most frequently used ones are trigonometric functions or Chebyshev and Legendre polynomials. In the spectral Galerkin method the test functions are the same as the trial functions, while in the spectral collocation method they are the translated Dirac delta functions centered at so-called collocation points. There are two types of formulations. In the modal formulation, the unknowns are the expansion coefficients. In contrast, for the nodal formulation, the unknowns are the nodal values of the unknown variables at the collocation points. For the Galerkin method, one can use either formulation. However, since the test functions and the trial functions are in general orthogonal to each other only in the modal space, the modal formulation will result in an uncoupled system, but not in the nodal formulation. For the collocation method, the nodal formulation is the more natural choice, and it always results in an uncoupled system since the delta functions are used as the test functions. For nonlinear fluxes, the pseudo-spectral method [30] is commonly used for computational efficiency, in which fluxes are calculated at nodal points using the nodal values of unknowns. For the modal formulation, this requires that unknowns and possibly their derivatives be first transferred to the nodal space, and then fluxes are transferred back to the modal space to compute their derivatives. While the original formulations were carried out in one dimension, their tensor products can be employed for problems in a simple multi-dimensional rectangle or box.

One of the major shortcomings of the traditional spectral methods is their restriction to problems in simple domains. Recent developments have extended these methods to multiple domains, including the spectral element (SE) method [31], [25] based on the Galerkin approach and the multi-domain spectral method [21], [22], [23], [24], [16] based on the collocation approach. Domains containing general quadrilateral and hexahedral elements can be handled by mapping the general elements to the standard elements using the conventional iso-parameterization. The SE method has also been extended to triangular and tetrahedral elements [20].

The SE method can be viewed as a high-order FE method with the nodal points placed at proper locations so that the spectral convergence can be obtained. In order to achieve local conservation for the FE or SE methods, the discontinuous Galerkin (DG) method was developed [32], [9], [10], [11], [12], [17]. Nodes on element boundaries are allowed to have multiple values, so that the local reconstruction in each element is in general discontinuous with that of its neighbors. The Galerkin MWR method is now applied locally to each element, using the local shape functions. As in the unstructured FV method, a Riemann solver is employed at element boundaries to compute the numerical fluxes. The integral conservation law is now satisfied for each element. While we must still solve a large set of coupled equations, each set involves only the unknowns in a few neighboring elements. Some of the integrals in the matrix entries involve quadratic terms. For non-linear flux functions, the required quadrature formulas must have twice the degree of precision as the precision of the reconstruction. In order to obtain stable and spectral convergence, unknowns for the DG method are normally placed at points where the reconstruction matrix is optimized. One choice involves the Fekete points where the determinant of the reconstruction matrix is maximized [38]. Other choices include points where the maximum Lebesgue constant of the reconstruction matrix is minimized [6], [7] or multivariate point sets through the electrostatic analogy [15]. In general, these points may not provide the necessary precisions of quadrature approximations for the surface and volume integrals, and one must therefore obtain solutions at other quadrature points through interpolations.

The universal local reconstruction concept inherent in the FE method can be utilized to overcome the computational inefficiencies of the more direct unstructured FV method. In the spectral volume (SV) method [40], [41], [42], [43], [28], each triangular or tetrahedral cell, here called a spectral volume (SV), is partitioned into structured subcells called control volumes (CV). These are polygons in 2D, and polyhedra in 3D. The latter can have non-planar faces, which must be subdivided into planar facets in order to perform flux integrations. The unknowns are now cell averages over the CV’s. If the SV’s are partitioned in a geometrically similar manner, a single, universal reconstruction results. Thus only a few coefficients need to be stored in advance to evaluate all flux integrals. For high orders of accuracy in 3D, the partitioning requires the introduction of a large number of parameters, whose optimization to achieve spectral convergence becomes increasingly more difficult. The growth in the number of interior facets and the increase in the number of quadrature points for each facet, raises the computational cost greatly. The computational cost of the SV method, and the difficulties in determining the parameters for spectral convergence, can both be significantly reduced if one were to apply the universal local reconstruction concept to the simpler unstructured FD method using nodal unknowns.

In this paper, we introduce a new, high-order, conservative, and efficient method, named the spectral difference (SD) method, for conservation laws on unstructured grids. The method combines the best features of structured and unstructured grid methods to obtain computational efficiency and geometric flexibility. It utilizes the concept of discontinuous and high-order local representations to achieve conservation and high accuracy in a manner similar to the DG and SV methods, but the new method is based on the finite-difference formulation to attain a simpler form and higher efficiency. Specifically, the differential form of the conservation laws is satisfied at nodal unknown points, with flux derivatives expressed in terms of values at flux points.

The paper is organized as follows. In the next section, we first describe the basic formulation of the method. Some representative placements of the unknown and flux points with various orders of accuracy for triangular elements are then presented in Section 3. Accuracy studies of the method are carried out with the two-dimensional linear wave equation and Burgers’ equation in Section 4. Each order of accuracy is numerically verified with five unstructured grids of consecutive refinement. Numerical solutions of plane electromagnetic waves incident on perfectly conducting circular cylinders are presented and compared with the exact solutions also in Section 4. Finally, some concluding remarks and suggestions for future study are given in Section 5.