The Piecewise Cubic Method (PCM) for Computational Fluid Dynamics

# The Piecewise Cubic Method (PCM) for Computational Fluid Dynamics

Dongwook Lee Hugues Faller Adam Reyes Applied Mathematics and Statistics, University of California, Santa Cruz, CA, U.S.A Department of Physics, University of California, Santa Cruz, CA, U.S.A Département de Physique, École Normale Supérieure, Paris, France
###### Abstract

We present a new high-order finite volume reconstruction method for hyperbolic conservation laws. The method is based on a piecewise cubic polynomial which provides its solutions a fifth-order accuracy in space. The spatially reconstructed solutions are evolved in time with a fourth-order 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 fourth-order Runge-Kutta method. We demonstrate that the solutions of PCM converges in fifth-order 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:
High-order methods; piecewise cubic method; finite volume method; gas dynamics; magnetohydrodynamics; Godunov’s method.
journal: Journal of Computational Physics

## 1 Introduction

In this paper we are interested in solving multidimensional conservation laws of the Euler equations and the ideal MHD equations, written as

 ∂U∂t+∇⋅F(U)=0, (1)

where is the vector of the conservative variables, and

 F(U)=[F(U),G(U),H(U)]T=[F,G,H]T (2)

is the flux vector.

We present a new high-order 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 high-order accuracy in smooth flows and well-resolved solutions in shock/discontinuous flows. With the advent of high-performance 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 high-order property with great shock-capturing 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 convergence-to-solution rates are the key advantages in high-order schemes. The advantage of using high-order 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 high-performance 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 high-order 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 “dimension-by-dimension” approach for simplicity. Although this approach has an advantage in its simplicity, it unfortunately fails to retain the high-order accurate property of the 1D baseline algorithm in multidimensional problems. Instead, it provides only a second-order accuracy in multidimensional nonlinear advection in finite volume method due to the lack of accuracy in approximating a face-averaged flux function as a result of mis-using 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 high-order 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 semi-discrete form,

 ∂¯¯¯¯¯Ui∂t=−1Δx(Fi+1/2−Fi−1/2) (3)

to get an equation for the evolution of the volume averaged variables, . Typically to achieve high-order accuracy in time the temporal update is done using a TVD Runge-Kutta scheme in method-of-lines form shu1988tvd (); mccorquodale2011high (). In this approach the high-order 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 predictor-corrector type temporal update formulation in which we take the time-average of Eq. (3)

 ¯¯¯¯¯Un+1i=¯¯¯¯¯Uni−ΔtΔx(Fn+1/2i+1/2−Fn+1/2i−1/2). (4)

Here is the volume averaged quantity at , and is the time average of the interface flux from to . In this way high-order 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 fifth-order 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 fourth-order accurate temporal updating scheme of PCM using a predictor-corrector type characteristic tracing, which is much simpler than the typical high-order Runge-Kutta ODE updates. In Section 5 we discuss how to extend the 1D scheme in Section 2 to multiple spatial dimensions following the dimension-by-dimension approach This approach is the same as the Class A approach in zhang2011order (), and should not be confused with the so-called 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 One-Dimensional Piecewise Cubic Method (PCM) Spatial Reconstruction

In this section we describe a new PCM scheme in a one-dimensional finite volume formulation for solving hyperbolic conservation laws of hydrodynamics and magnetohydrodynamics. The new PCM scheme is a higher-order extension of Godunov’s method godunov1959difference (), bearing its key components in the reconstruction algorithm on the relevant ideas of its high-order predecessors, the PPM scheme colella1984piecewise (), the WENO schemes jiang1996efficient (); shi2002technique (); borges2008improved (); castro2011high (), and Hermite-WENO schemes qiu2004hermite (); qiu2005hermite (); zhu2009hermite (); balsara2007sub ().

For the purpose of this section, we take the hyperbolic system of conservation laws of the 1D Euler equations

 ∂U∂t+∂F(U)∂x=0. (5)

The notations used are the vector of the conservative variables and fluxes , respectively, defined as

 U=⎡⎢⎣ρρuE⎤⎥⎦,F(U)=⎡⎢⎣ρuρu2+pu(E+p)⎤⎥⎦. (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,

 E=pγ−1+ρu22, (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,

 ∂W∂t+Λ∂W∂x=0. (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 third-order especially with better non-oscillatory controls, in particular when considering wave-by-wave 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

 pi(x)=c0+c1(x−xi)+c2(x−xi)2+c3(x−xi)3. (9)

The goal is now to determine the four coefficients , , , which can be achieved by imposing the following four conditions:

 1Δx∫Iipi(x)dx = ¯¯¯qi, (10) pi(xi−1/2) = qL,i, (11) pi(xi+1/2) = qR,i, (12) p′i(xi) = q′C,i, (13)

where

 ¯¯¯qi=1Δx∫Iiq(x,tn)dx (14)

is the cell-averaged quantity at on which is given as an initial condition;

 qL,i=q(xi−1/2,tn)+O(Δxp),qR,i=q(xi+1/2,tn)+O(Δxp) (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

 q′C,i=q′(xi,tn)+O(Δxr) (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:

 c0+c2Δx212 = ¯¯¯qi, (17) c0−c1Δx2+c2Δx24−c3Δx38 = qL,i, (18) c0+c1Δx2+c2Δx24+c3Δx38 = qR,i, (19) c1 = q′C,i, (20)

which, in turn, can be solved for all four , . The final expressions of the coefficients in terms of , and are given as:

 c0 = 14(−qR,i−qL,i+6¯¯¯qi), (21) c1 = q′C,i, (22) c2 = 3Δx2(qR,i+qL,i−2¯¯¯qi), (23) c3 = 4Δx3(qR,i−qL,i−Δxq′C,i). (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 .

We now devote the following sections to describe how to determine , and so that the resulting PCM approximation to the variable lend its accuracy a fifth-order in space (Sections 2.2 and 2.3) and a fourth-order in time (Section 3).

### 2.2 Reconstruction of the Riemann States qL,i and qR,i

We follow the fifth-order finite volume WENO approach, either of the classical WENO-JS jiang1996efficient () or WENO-Z borges2008improved (); castro2011high (), in order to reconstruct the left and right Riemann states, and , on each cell . For the sake of providing a full self-contained 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 sub-stencils, , , each of which consisting three cells , . Let us first define

 S1 = {Ii−2,Ii−1,Ii}, (25) S2 = {Ii−1,Ii,Ii+1}, (26) S3 = {Ii,Ii+1,Ii+2}. (27)

Formulating the WENO reconstruction consists of the following three steps:

##### Step 1: ENO-Build

We begin with building a second degree polynomial for each ,

 pℓ(x)=2∑k=0aℓ,k(x−xi)k, (28)

each of which is defined on , satisfying

 1Δx∫Ikpℓ(x)dx=¯qk, (29)

for . After a bit of algebra, we obtain the coefficients that determine in Eq. (28).

For ,

 a1,0=(−124¯qi−2+112¯qi−1+2324¯qi), (30) a1,1=(12¯qi−2−2¯qi−1+32¯qi)1Δx, (31) a1,2=(12¯qi−2−¯qi−1+12¯qi)1Δx2, (32)

and for ,

 a2,0=(−124¯qi−1+1312¯qi−124¯qi+1), (33) a2,1=(−12¯qi−1+12¯qi+1)1Δx, (34) a2,2=(12¯qi−1−¯qi+12¯qi+1)1Δx2. (35)

Lastly, for , we get

 a3,0=(2324¯qi+112¯qi+1−124¯qi+2), (36) a3,1=(−32¯qi+2¯qi+1−12¯qi+2)1Δx, (37) a3,2=(12¯qi−¯qi+1+12¯qi+2)1Δx2. (38)

Then the three sets of left and right states follow as

 {p1(xi−1/2),p2(xi−1/2),p3(xi−1/2)}, and {p1(xi+1/2),p2(xi+1/2),p3(xi+1/2)}, (39)

where each of is the ENO approximation and is given by, first for ,

 p1(xi−1/2)=−16¯qi−2+56¯qi−1+13¯qi, (40) p1(xi+1/2)=13¯qi−2−76¯qi−1+116¯qi, (41)

and for ,

 p2(xi−1/2)=13¯qi−1+56¯qi−16¯qi+1, (42) p2(xi+1/2)=−16¯qi−1+56¯qi+13¯qi+1, (43)

and finally for ,

 p3(xi−1/2)=116¯qi−76¯qi+1+13¯qi+2, (44) p3(xi+1/2)=13¯qi+56¯qi+1−16¯qi+2. (45)

These left and right states respectively approximate the pointwise values at the interfaces with third-order accuracy, i.e., (see jiang1996efficient ()) by using the given cell-averaged quantities .

##### Step 2: Linear Constant Weights

The next step is to construct a fourth-degree polynomial

 ϕ(x)=4∑k=0bk(x−xi)k (46)

over the entire stencil which satisfies

 1Δx∫Ikϕ(x)dx=¯qk, (47)

for . We can show that the coefficients are given as

 b0=3640¯qi−2−29480¯qi−1+1067960¯qi−29480¯qi+1+3640¯qi+2, (48) b1=(548¯qi−2−1724¯qi−1+1724¯qi+1−548¯qi+2)1Δx, (49) b2=(−116¯qi−2+34¯qi−1−118¯qi+34¯qi+1−116¯qi+2)1Δx2, (50) b3=(−112¯qi−2+16¯qi−1−16¯qi+1+112¯qi+2)1Δx3, (51) (52)

WENO uses to determine three linear constant weights , , with , such that

 ϕ(xi±1/2)=3∑ℓ=1γ±ℓpℓ(xi±1/2). (53)

The values on the left-hand side become

 ϕ(xi−1/2)=−120¯qi−2+920¯qi−1+4760¯qi−1360¯qi+1+130¯qi+2, (54)

and

 ϕ(xi+1/2)=130¯qi−2−1360¯qi−1+4760¯qi+920¯qi+1−120¯qi+2. (55)

Now, by inspection, we obtain a set of linear weights for the left state,

 γ−1=310,γ−2=610,γ−3=110, (56)

and for the right state,

 γ+1=110,γ+2=610,γ+3=310. (57)
##### Step 3: Nonlinear Weights

The last step that imposes the non-oscillatory feature in the WENO approximations is to measure how smoothly the three polynomials vary on . This is done by determining non-constant, nonlinear weights (three of them for each state) that rely on the so-called smoothness indicator , defined by

 βℓ=2∑s=1(Δx2s−1∫Ii[dsdxspℓ(x)]2dx). (58)

With this definition becomes small for smooth flows, and large for discontinuous flows.

For explicit expressions, we attain

 β1=1312(¯qi−2−2¯qi−1+¯qi)2+14(¯qi−2−4¯qi−1+3¯qi)2, (59) β2=1312(¯qi−1−2¯qi+¯qi+1)2+14(¯qi−1−¯qi+1)2, (60) β3=1312(¯qi−2¯qi+1+¯qi+2)2+14(3¯qi−4¯qi+1+¯qi+2)2. (61)

Equipped with these , the nonlinear weights are defined as:

• For WENO-JS:

 ω±ℓ=~ω±ℓ∑s~ω±s, where ~ω±ℓ=γ±ℓ(ϵ+βℓ)m, (62)
• For WENO-Z:

 ω±ℓ=~ω±ℓ∑s~ω±s, where ~ω±ℓ=γ±ℓ(1+(|β0−β2|ϵ+βℓ)m). (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 fifth-order spatially accurate, non-oscillatorily reconstructed values at each cell interface at each time step jiang1996efficient (); shu2009high (),

 qL;R,i=3∑ℓ=1ω±ℓpℓ(xi±1/2). (64)

### 2.3 Reconstruction of the Derivative q′C,i

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 WENO-type 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 ENO-build strategy in Section 2.2 in which the three second degree ENO polynomials in Eq. (28) are constructed over the five-point stencil . However, this setup will provide only a third-order 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,

 βℓ=Δx3∫Ii[p′′ℓ(x)]2dx,ℓ=1,2,3. (65)

Through a Taylor expansion analysis on Eq. (65) we see

 βℓ=D(1+O(Δx)), (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

 ωℓ=γℓ+O(Δx), (67)

where the linear constant weights are assumed to exist, when is smooth in , such that

 q′C,i=3∑ℓ=1γℓp′ℓ(xi)=q′(xi,tn)+O(Δx3). (68)

This finally implies the accuracy of is found out to be third-order,

 q′C,i=3∑ℓ=1ωℓp′ℓ(xi)=q′(xi,tn)+O(Δx3), (69)

because

 3∑ℓ=1ωℓp′ℓ(xi)−3∑ℓ=1γℓp′ℓ(xi) = 3∑ℓ=1(ωℓ−γℓ)(p′ℓ(xi)−q′(xi,tn)) (70) = 3∑ℓ=1O(Δx)O(Δx2)=O(Δx3).

In the last equality, we used the fact that, for each , is only a first degree polynomial which is accurate up to second-order when approximating .

For this reason, we want a better strategy to obtain an approximation at least fourth-order accurate in order that the overall nominal accuracy of the 1D PCM scheme achieves at least fourth-order accurate in both space and time.

##### Step 1: PPM-Build

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 (),

 ϕ±(x)=3∑k=0a±k(x−xi±1/2)k. (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 :

 1Δx∫Ikϕ−(x)dx=¯qnk, % for i−2≤k≤i+1, (72)

and

 1Δx∫Ikϕ+(x)dx=¯qnk, % for i−1≤k≤i+2. (73)

After a bit of algebra we obtain the coefficients , with for , while for :

 a±0=112(−¯qi−2+s+7¯qi−1+s+7¯qi+s−¯qi+1+s), (74)
 a±1=112Δx(¯qi−2+s−15¯qi−1+s+15¯qi+s−¯qi+1+s), (75)
 a±2=14Δx2(¯qi−2+s−¯qi−1+s−¯qi+s+¯qi+1+s), (76)
 a±3=16Δx3(−¯qi−2+s+3¯qi−1+s−3¯qi+s+¯qi+1+s). (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 ,

 q′C,i=γ−ϕ′−(xi)+γ+ϕ′+(xi). (78)

The two linear weights can be determined by comparing Eq. (78) with in Eq. (46),

 γ−ϕ′−(xi)+γ+ϕ′+(xi)=ϕ′(xi) (79)

This gives us

 γ−(a−1+a−2Δx+3a−3Δx24)+γ+(a+1−a+2Δx+3a+3Δx24)=b1 (80)

where is defined in Eq. (49). By inspection, we obtain

 γ−=γ+=12. (81)
##### Step 3: Nonlinear Weights

The smoothness indicators are now constructed using as

 β±=3∑s=2(Δx2s−1∫Ii[dsdxsϕ±(x)]2dx). (82)

They can be written explicitly as

 β− = 4(a−2)2Δx4+12(a−2)(a−3)Δx5+48(a−3)2Δx6 = 14(¯qi−2−¯qi−1−¯qi+¯qi+1)2 + 12(¯qi−2−¯qi−1−¯qi+¯qi+1)(−¯qi−2+3¯qi−1−3¯qi+¯qi+1) + 43(−¯qi−2+3¯qi−1−3¯qi+¯qi+1)2, (84)

and

 β+ = 4(a+2)2Δx4−12(a+2)(a+3)Δx5+48(a+3)2Δx6 = 14(¯qi−1−¯qi−¯qi+1+¯qi+2)2 + 12(¯qi−1−¯qi−¯qi+1+¯qi+2)(−¯qi−1+3¯qi−3¯qi+1+¯qi+2) + 43(−¯qi−1+3¯qi−3¯qi+1+¯qi+2)2. (86)

Upon conducting Taylor series expansion analysis on , we can see that

 β±=D(1+O(Δx)), (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 WENO-JS:

 ω±=~ω±~ω−+~ω+, where ~ω±=γ±(ϵ+β±)m, (88)
• For WENO-Z:

 ω±=~ω±~ω−+~ω+, where ~ω±=γ±(1+(|β+−β−|ϵ+β±)m). (89)

The final representation of the approximation becomes

 q′C,i=ω−ϕ′−(xi)+ω+ϕ′+(xi). (90)

Let us now verify that this approximation is fourth-order accurate after all, that is,

 q′C,i=ω−ϕ−(xi)+ω+ϕ+(xi)=q′(xi,tn)+O(Δx4). (91)

Similarly as before, using Eqs. (87), (88), and (89), we can see that, with the help of the binomial series expansion,

 ω±=γ±+O(Δx). (92)

Therefore, the desired accuracy claimed in Eq. (91) is readily verified by repeating the similar relationship in Eq. (70):

 ∑ℓ=−,+ωℓϕ′ℓ(xi)−∑ℓ=−,+γℓϕ′ℓ(xi) = ∑ℓ=−,+(ωℓ−γℓ)(ϕ′ℓ(xi)−q′(xi,tn)) (93) = ∑ℓ=−,+O(Δx)O(Δx3)=O(Δx4).

Comparing Eq. (93) with Eq. (70), we now see that it is fourth-order accurate due to the improved third-order accuracy in calculating