The Piecewise Cubic Method (PCM) for Computational Fluid Dynamics
Abstract
We present a new highorder finite volume reconstruction method for hyperbolic conservation laws. The method is based on a piecewise cubic polynomial which provides its solutions a fifthorder accuracy in space. The spatially reconstructed solutions are evolved in time with a fourthorder accuracy by tracing the characteristics of the cubic polynomials. As a result, our temporal update scheme provides a significantly simpler and computationally more efficient approach in achieving fourth order accuracy in time, relative to the comparable fourthorder RungeKutta method. We demonstrate that the solutions of PCM converges in fifthorder in solving 1D smooth flows described by hyperbolic conservation laws. We test the new scheme in a range of numerical experiments, including both gas dynamics and magnetohydrodynamics applications in multiple spatial dimensions.
keywords:
Highorder methods; piecewise cubic method; finite volume method; gas dynamics; magnetohydrodynamics; Godunov’s method.1 Introduction
In this paper we are interested in solving multidimensional conservation laws of the Euler equations and the ideal MHD equations, written as
(1) 
where is the vector of the conservative variables, and
(2) 
is the flux vector.
We present a new highorder piecewise cubic method (PCM) algorithm that is extended from the classical PPM and WENO schemes colella1984piecewise (); jiang1996efficient (). These two algorithms, by far, have been extremely successful in various scientific fields where there are challenging computational needs for both highorder accuracy in smooth flows and wellresolved solutions in shock/discontinuous flows. With the advent of highperformance computing (HPC) in recent years, such needs have been more and more desired, and have become a necessary requirement in conducting large scale, cutting edge simulations of gas dynamics and magnetohydrodynamics (MHD) Dongarra2010 (); Dongarra2012 (); Subcommittee2014 (); Keyes2013 ().
As observed in the success stories of the PPM and WENO methods, discrete algorithms of data interpolation and reconstruction play a key role in numerical methods for PDE approximations LeVeque2002 (); leveque2007finite (); Toro2009 () within the broad framework of finite difference and finite volume discretization methods. In view of this, computational improvements of such interpolation and reconstruction schemes, particularly focused on the highorder property with great shockcapturing capability, take their positions at the center of HPC in modern computational fluid dynamics.
The properties of enhanced solution accuracy with lower numerical errors on a given grid resolution and faster convergencetosolution rates are the key advantages in highorder schemes. The advantage of using highorder methods in HPC is therefore clear: one can obtain reproducible, admissible, and highly accurate numerical solutions in a faster computational time at the expense of increased rate of floating point operations, while at the same time, with the use of smaller size of grid resolutions. This is by no means exceedingly efficient in highperformance computing (HPC), in view of the fact that the increase of grid resolutions has a direct impact to an increase of memory footprints which are bounded in all modern computing architecture.
In this regards, our goal in this paper is to lay down a mathematical foundation in designing a new highorder method using piecewise cubic polynomials. We mainly focus on describing the detailed PCM algorithm in 1D finite volume framework for the scope of the current paper. For multidimensional problems, we adopt the classical “dimensionbydimension” approach for simplicity. Although this approach has an advantage in its simplicity, it unfortunately fails to retain the highorder accurate property of the 1D baseline algorithm in multidimensional problems. Instead, it provides only a secondorder accuracy in multidimensional nonlinear advection in finite volume method due to the lack of accuracy in approximating a faceaveraged flux function as a result of misusing an averaged quantity in place of a pointwise quantity, or vice versa shu2009high (); buchmuller2014improved (); zhang2011order (); mccorquodale2011high (). Although the baseline 1D PCM scheme can be extended to multiple spatial dimensions preserving its highorder accuracy by following more sophisticated treatments shu2009high (); buchmuller2014improved (); zhang2011order (); mccorquodale2011high (), more careful work is needed to carry out the detailed design, and this will be considered in our future research.
For a finite volume scheme in 1D we take the spatial average of Eq. (1) over the cell , yielding a semidiscrete form,
(3) 
to get an equation for the evolution of the volume averaged variables, . Typically to achieve highorder accuracy in time the temporal update is done using a TVD RungeKutta scheme in methodoflines form shu1988tvd (); mccorquodale2011high (). In this approach the highorder accuracy comes from taking the multiple Euler stages of the RK time discretizations, which require repeated reconstructions in a single time step, increasing the computational cost.
Instead, as will be fully described in Section 2, one of the novel ideas in PCM is to employ the simple single stage predictorcorrector type temporal update formulation in which we take the timeaverage of Eq. (3)
(4) 
Here is the volume averaged quantity at , and is the time average of the interface flux from to . In this way highorder in space and time is accomplished with a single reconstruction in contrast to the multiple Euler stages of the RK time discretizations, providing significant benefits in computational efficiency per solution accuracy.
The organization of the paper is as follows: Section 2 describes the fifthorder accurate spatial reconstruction algorithm of PCM in 1D. We highlight several desirable properties of the PCM scheme in terms of computational efficiency and solution accuracy. Section 3 introduces the fourthorder accurate temporal updating scheme of PCM using a predictorcorrector type characteristic tracing, which is much simpler than the typical highorder RungeKutta ODE updates. In Section 5 we discuss how to extend the 1D scheme in Section 2 to multiple spatial dimensions following the dimensionbydimension approach ^{}^{}This approach is the same as the Class A approach in zhang2011order (), and should not be confused with the socalled dimensionally split approach. Our spatial integration scheme in this paper is directionally unsplit. buchmuller2014improved (); zhang2011order ().
In Section 6 we test the PCM scheme on a wide spectrum of benchmark problems in 1D, 2D and 3D, both for hydrodynamics and magnetohydrodynamics (MHD) applications. We also compare the PCM solutions with PPM and WENO solutions in order to examine numerical accuracy, capability and efficiency in both smooth and shock flow regimes. We conclude our paper in Section 7 with a brief summary.
2 The OneDimensional Piecewise Cubic Method (PCM) Spatial Reconstruction
In this section we describe a new PCM scheme in a onedimensional finite volume formulation for solving hyperbolic conservation laws of hydrodynamics and magnetohydrodynamics. The new PCM scheme is a higherorder extension of Godunov’s method godunov1959difference (), bearing its key components in the reconstruction algorithm on the relevant ideas of its highorder predecessors, the PPM scheme colella1984piecewise (), the WENO schemes jiang1996efficient (); shi2002technique (); borges2008improved (); castro2011high (), and HermiteWENO schemes qiu2004hermite (); qiu2005hermite (); zhu2009hermite (); balsara2007sub ().
For the purpose of this section, we take the hyperbolic system of conservation laws of the 1D Euler equations
(5) 
The notations used are the vector of the conservative variables and fluxes , respectively, defined as
(6) 
Here is the fluid density, is the fluid velocity in direction, and is the total energy as the sum of the internal energy and the kinetic energy obeying the ideal gas law,
(7) 
where is the gas pressure, with the ratio of specific heats denoted as . We denote the cells in direction by . We assume our grid is configured on an equidistant uniform grid for simplicity.
In addition to the system of the Euler equations in the conserved variables as given in Eq. (5), we often use the two other equivalent system of equations each of which can be written either in the primitive variables or in the characteristic variables . The characteristic variable is readily obtained from or by multiplying the left eigenvectors corresponding to either or , for instance, . In the latter ) represents the matrix obtained from diagonalizing the coefficient matrix , whose rows are the th left eigenvectors , . The representation of the system in furnishes a completely linearly decoupled 1D system of equations,
(8) 
The above system in the characteristic variables is therefore very handy for analyses, and also is a preferred choice of variable in order to furnish numerical solutions more accurate than thirdorder especially with better nonoscillatory controls, in particular when considering wavebywave propagations in a system of equations shu2009high (). In this reason the characteristic variable is taken as our default variable choice in the 1D PCM reconstruction steps via characteristic decompositions, albeit with an increased computational cost, among the other two choices of the primitive or the conservative variables .
The methodology presented below can be similarly applied to the 1D ideal MHD equations (see for instance, brio1988upwind ()).
2.1 Piecewise Cubic Profile
To begin with we first define a cubic polynomial to approximate a th characteristic variable on each interval by
(9) 
The goal is now to determine the four coefficients , , , which can be achieved by imposing the following four conditions:
(10)  
(11)  
(12)  
(13) 
where
(14) 
is the cellaveraged quantity at on which is given as an initial condition;
(15) 
are respectively the th order accurate pointwise left and the right Riemann states at on the cell that are unknown yet but are to be determined as described below; and lastly
(16) 
is the th order accurate approximation to the slope of at evaluated at , which is again unknown at this point but is to be determined as below.
For the moment let us assume that all four quantities and are known. It can be shown that the system of relations in Eqs. (10) (13) is equivalent to a system given as:
(17)  
(18)  
(19)  
(20) 
which, in turn, can be solved for all four , . The final expressions of the coefficients in terms of , and are given as:
(21)  
(22)  
(23)  
(24) 
Therefore once we figure out the three unknowns, , and , the cubic profile in Eq. (9) can be completely determined and is ready to approximate on each .
2.2 Reconstruction of the Riemann States and
We follow the fifthorder finite volume WENO approach, either of the classical WENOJS jiang1996efficient () or WENOZ borges2008improved (); castro2011high (), in order to reconstruct the left and right Riemann states, and , on each cell . For the sake of providing a full selfcontained description of the PCM scheme, we briefly present the two WENO Riemann state reconstruction strategies here.
The main idea in WENO is to employ its reconstruction procedure according to the nonlinear smoothness measurements on three ENO substencils, , , each of which consisting three cells , . Let us first define
(25)  
(26)  
(27) 
Formulating the WENO reconstruction consists of the following three steps:
Step 1: ENOBuild
We begin with building a second degree polynomial for each ,
(28) 
each of which is defined on , satisfying
(29) 
for . After a bit of algebra, we obtain the coefficients that determine in Eq. (28).
For ,
(30)  
(31)  
(32) 
and for ,
(33)  
(34)  
(35) 
Lastly, for , we get
(36)  
(37)  
(38) 
Then the three sets of left and right states follow as
(39) 
where each of is the ENO approximation and is given by, first for ,
(40)  
(41) 
and for ,
(42)  
(43) 
and finally for ,
(44)  
(45) 
These left and right states respectively approximate the pointwise values at the interfaces with thirdorder accuracy, i.e., (see jiang1996efficient ()) by using the given cellaveraged quantities .
Step 2: Linear Constant Weights
The next step is to construct a fourthdegree polynomial
(46) 
over the entire stencil which satisfies
(47) 
for . We can show that the coefficients are given as
(48)  
(49)  
(50)  
(51)  
(52) 
WENO uses to determine three linear constant weights , , with , such that
(53) 
The values on the lefthand side become
(54) 
and
(55) 
Now, by inspection, we obtain a set of linear weights for the left state,
(56) 
and for the right state,
(57) 
Step 3: Nonlinear Weights
The last step that imposes the nonoscillatory feature in the WENO approximations is to measure how smoothly the three polynomials vary on . This is done by determining nonconstant, nonlinear weights (three of them for each state) that rely on the socalled smoothness indicator , defined by
(58) 
With this definition becomes small for smooth flows, and large for discontinuous flows.
For explicit expressions, we attain
(59)  
(60)  
(61) 
Equipped with these , the nonlinear weights are defined as:

For WENOJS:
(62) 
For WENOZ:
(63)
Here is any arbitrarily small positive number that prevents division by zero, for which we choose . One of the classical choice of in many WENO literatures is found to be jiang1996efficient (); shi2002technique (); however, it was suggested in borges2008improved () that should be chosen to be much smaller in order to force this parameter to play only its original role of avoiding division by zero in the definitions of the weights, Eqs. (62) and (63).
Another closely related point of discussion is with the value of , the power in the denominators in Eqs. (62) and (63). The parameter determines the rate of changes in , and most of the WENO literatures use . However, we observe that using resolves discontinuities sharper in most of our numerical simulations without exhibiting any numerical instability, so became the default value in our implementation. For more detailed discussions on the choices of and , see borges2008improved (); castro2011high ().
Using these nonlinear weights, we complete the WENO reconstruction procedure of producing the fifthorder spatially accurate, nonoscillatorily reconstructed values at each cell interface at each time step jiang1996efficient (); shu2009high (),
(64) 
2.3 Reconstruction of the Derivative
The spatial reconstruction part of PCM proceeds to the next final step to obtain the derivative in Eq. (13). The approach is again to take the WENOtype reconstruction as before, but this time, to approximate a first derivative of a function shu2009high (), i.e., .
For this, we might consider using the same ENObuild strategy in Section 2.2 in which the three second degree ENO polynomials in Eq. (28) are constructed over the fivepoint stencil . However, this setup will provide only a thirdorder accurate approximation to the exact derivative . To see this, we first observe that the smoothness indicators with this setup will be including only a single term,
(65) 
Through a Taylor expansion analysis on Eq. (65) we see
(66) 
where is a nonzero quantity independent of but may depend on , assuming on . This results in a set of three nonlinear weights , , obtained either by Eq. (62) or Eq. (63), satisfying
(67) 
where the linear constant weights are assumed to exist, when is smooth in , such that
(68) 
This finally implies the accuracy of is found out to be thirdorder,
(69) 
because
(70)  
In the last equality, we used the fact that, for each , is only a first degree polynomial which is accurate up to secondorder when approximating .
For this reason, we want a better strategy to obtain an approximation at least fourthorder accurate in order that the overall nominal accuracy of the 1D PCM scheme achieves at least fourthorder accurate in both space and time.
Step 1: PPMBuild
An alternate strategy for this goal therefore would be to use a set of third degree polynomials instead. This can be designed using the two third degree polynomials, , from the PPM algorithm colella1984piecewise (),
(71) 
Following the description of PPM, we carry out to determine the coefficients by imposing the following constraints on that are essential to keeping the volume averages on each cell :
(72) 
and
(73) 
After a bit of algebra we obtain the coefficients , with for , while for :
(74) 
(75) 
(76) 
(77) 
Step 2: Linear Constant Weights
Now that the polynomials are determined over the stencil , we use their first derivatives to obtain a convex combination with two linear weights and ,
(78) 
The two linear weights can be determined by comparing Eq. (78) with in Eq. (46),
(79) 
This gives us
(80) 
where is defined in Eq. (49). By inspection, we obtain
(81) 
Step 3: Nonlinear Weights
The smoothness indicators are now constructed using as
(82) 
They can be written explicitly as
(84) 
and
(86) 
Upon conducting Taylor series expansion analysis on , we can see that
(87) 
where is a nonzero quantity independent of but might depend on , assuming on .
The remaining procedure is to obtain the two nonlinear weights in the similar way done in the edge reconstructions in Eqs. (62) – (63),

For WENOJS:
(88) 
For WENOZ:
(89)
The final representation of the approximation becomes
(90) 
Let us now verify that this approximation is fourthorder accurate after all, that is,
(91) 
Similarly as before, using Eqs. (87), (88), and (89), we can see that, with the help of the binomial series expansion,
(92) 
Therefore, the desired accuracy claimed in Eq. (91) is readily verified by repeating the similar relationship in Eq. (70):
(93)  
Comparing Eq. (93) with Eq. (70), we now see that it is fourthorder accurate due to the improved thirdorder accuracy in calculating