A Gas Dynamics Method Based on the Spectral Deferred Corrections ( SDC ) Time Integration Technique and the Piecewise Parabolic Method ( PPM )

We present a computational gas dynamics method based on the Spectral Deferred Corrections (SDC) time integration technique and the Piecewise Parabolic Method (PPM) finite volume method. The PPM framework is used to define edge-averaged quantities, which are then used to evaluate numerical flux functions. The SDC technique is used to integrate solution in time. This kind of approach was first taken by Anita et al in [1]. However, [1] is problematic when it is implemented to certain shock problems. Here we propose significant improvements to [1]. The method is fourth order (both in space and time) for smooth flows, and provides highly resolved discontinuous solutions. We tested the method by solving variety of problems. Results indicate that the fourth order of accuracy in both space and time has been achieved when the flow is smooth. Results also demonstrate the shock capturing ability of the method.


Introduction
In this paper, we present a conservative scheme based on the SDC and PPM methods.The integration of the SDC method to the PPM method was first carried out by Anita et al in [1].However, [1] is problematic in the sense that it is oscillatory when it is applied to certain shock problems unless complicated extra repair steps are introduced.The oscillations are developed around the shocks at earlier times, then spread into entire computational region.These oscillations are neutral oscillations (extraneous wiggles) that pollute the solution and never disappear.The main reason for having such wiggles is that the PPM fluxes (without the time averaging in SDC framework) are evaluated in a way to obtain sharper discontinuous profiles, therefore lack of necessary numerical diffusion mechanism.Here, we introduce a strategy that eliminates the oscillatory behavior of [1].
The SDC-PPM method falls in to the class of higher -order high-resolution-schemes.Here, we shall provide a short historical perspective to higher-order high-resolution-schemes. High resolution schemes are designed to solve gas dynamics equations (Gas dynamics equations are also referred to as the inviscid Euler equations or the conservation laws.Hereafter, we will use all of three terminologies interchangeably).In the last several decades, many numerical schemes were introduced for this purpose.Godunov [2] initiated a novel approach that is now accepted as one of the main building blocks for construction of a high resolution scheme.Godunov supposed that the initial data could be replaced by a piecewise constant set of states (i.e, cell averages of the initial solution) with discontinuities located at computational cell edges (cell faces in 3-D).He then found exact solutions to this simplified problem by locally performing the well-known Riemann solution theory.Finally, he replaced exact solutions by a set of piecewise constant approximations (i.e, exact solutions are averaged down to the local cells).Godunov's method revolutionized the field of computational gas dynamics, by overcoming many of the difficulties that have been persistent for many years.However, Godunov's method was only first order accurate.The first major improvement to Godunov's work was made by van Leer [3] who approximated the initial data and solutions at each subsequent time level by piecewise linear segments allowing discontinuities be-tween the segments.Van Leer's work, also known as the MUSCL (Monotonic Upwind Scheme for Conservation Laws) scheme, raised the order of accuracy of the Godunov's method to two.Later, van Leer's MUSCL scheme was reconsidered or revised by others [4][5][6].For instance, Colella's MUSCL scheme [5] defines the slopes for the linear reconstruction based on the average values as oppose to van Leer [3] treats them as separate variables.When the flow is smooth, Colella's slope definition corresponds to a fourth order finite difference approximation to the derivative of a given state variable.Colella's work [5] was another step forward to increase the accuracy of the Gudonov's method.So far the Godunov type methods used constant or linear piecewise reconstructions.In [7], Colella and Woodward introduced a new method which is based on a piecewise parabolic reconstruction of solutions.Their method, famously known as the Piecewise Parabolic Method (PPM), achieves more accurate (fourth order in space) solution representations for smooth flows as well as it captures steeper shock discontinuities.
The MUSCL schemes [4][5][6] or the PPM method [7] are the few examples of higher order high-resolutionschemes.There exists number of other high-resolution schemes such as ENO [8], WENO [9], TVD [10], FCT [11], PHM [12], and LLR [13,14] methods.In several instances, comprehensive studies have been carried out to compare the performance of these schemes [15][16][17][18][19].We remark that these above cited studies mostly favor the PPM method.A particularly attractive feature of PPM method is that it can produce fourth order accurate calculations as long as the solutions stay smooth.We exploited this aspect when we studied zero Mach number flow problems (our collaborated work with M. Minion in [20]).The drawbacks of the PPM method reveal mostly when it is applied to shock problems.For instance, the PPM method [7] needs several external repairs in order to produce discontinuous solutions without spurious oscillations.The oscillatory behavior is associated with the method not having enough numerical diffusion (dissipation).In order to add appropriate numerical diffusion, the fluxing steps are modified in a way that the so-called smoothening and flattening algorithms have to be performed.We note that the implementation of these additional steps can be complicated.In our work, we avoid these external fixes by using a rather simple approach which will be explained next.
The SDC method is based on a number of deferred corrections of a low order provisional solution, which is predicted by forward or backward (depending on the stiffness of the problem) Euler method in order to achieve higher order of accuracy in time.In our case, we perform three deferred corrections to obtain fourth order accurate solutions.The SDC-PPM method of [1] uses the PPM fluxing procedure to evaluate the numerical flux functions and employs the SDC method to integrate solutions in time.We observed that the SDC-PPM method develops neutrally stable oscillations at the prediction step and maintains these oscillations during the deferred correction iterations.The main reason for this behavior (as mentioned above) is that the higher order flux evaluation at the prediction step lacks of necessary numerical diffusion so as the correction steps leading to unwanted oscillations around discontinuities.We fix this behavior with the following strategy.We know from our observation that one has to use more diffusive fluxes at the prediction step to avoid potential oscillations.Thus instead of full PPM fluxing, we employ an up-winding procedure that is naturally more diffusive.On the other hand, we let the correction iterations include the higher order PPM fluxing.With this strategy, we found out that solutions from the prediction step have enough numerical diffusion so that potential oscillations that might come from the correction steps are killed off.One question about this strategy is that does the up-winding procedure excessively smear the discontinuities?Our finding indicates that although the shocks are smeared at the prediction step, the correction iterations sharpen them back.In fact, we have obtained the same sharpness (and same shock amplitudes in that matter) as the full PPM would predict.We remark that the full PPM method [7] as well as the SDC-PPM method [1] has to make use of the smoothening, flattening, and artificial diffusion steps in order to keep solutions oscillation free.We avoid all of these extra steps in our approach.
The organization of this paper is as follows.In Section 2, the governing equations are presented.In Section 3, some notational conventions and the numerical algorithm are described.In Section 4, the computational results are presented.In Section 5, our concluding remarks are given.

Governing Equations
The in viscid Euler partial differential equations are widely used to model gas dynamics problems.These equations describe the physical evolution of conserved quantities such as mass, momentum, and total energy in space and time.Therefore, they are often referred to as the conservation laws.The conservation laws fall into class of the hyperbolic partial equations that admit discontinuous solutions such as shock waves.Therefore these equations are the best model for a gas problem in which shock waves frequently occur.
We consider the following two dimensional conservation equations, where U is the vector of conserved quantities, and F(U) and G(U) are the flux functions in x-and y-directions.For instance, where and E denote density, velocity, pressure, and the total energy.E satisfies This equation is called the equation of state for an ideal polytropic gas.An ideal gas obeying (2) is also referred to as the gamma-law gas with  denoting the specific heat ratio.Under ordinary circumstances air is composed primarily of 2 and 2 , so . We will use 1.4   in all of our computations.

Notations and Formulation
In this section, we briefly review several notational conventions that we have introduced in [20].For ease of presentation, we assume that the physical domain is two-dimensional and divided into a uniform array of cells of width and height h.Let the cell with center at x y be denoted by the pair (i, j), and let the half integer subscripts i+ 1/2 and j + 1/2 denote a shift by distance h/2 in the x-and y-direction respectively.The extension to rectangular cells and three dimensions is straightforward.
The PPM method or many other high resolution schemes rely on the so called finite-volume approach.The finite-volume approach is based on the evolution of cell averages.A cell average for some quantity   and Figure 1 gives a graphical illustration of the above definitions.To specify the finite volume formulation of the conservation law where   , , U x y t is the vector of conserved quantities and 6) is equivalent to Equation ( 1)), we integrate Equation ( 6) over the computational cells and use the divergence theorem to attain , , 0 where the flux integral above is defined as Using the definitions of the edge average in Equation ( 4)-(5) Therefore, Furthermore, if we let , , 0.

The Fourth Order Spectral Deferred Corrections (SDC) Time Integration Algorithm
The construction of stable and accurate numerical methods for solving initial value problems is a well studied subject [21][22][23].Runga-Kutta methods or linear multi-step methods are highorder discretization schemes by construction.On the other hand, Richardson extrapolation or deferred correction methods are based on increasing the accuracy of a low order method through iterations.The second group is also known as the predictor-corrector techniques.High order direct discretization schemes can be expensive due to the stability constraint.In this case, most practitioners recommend the predictor-corrector techniques.
In this paper, we consider the spectral deferred corrections method.The spectral deferred corrections (SDC) method, introduced in [24], is a known stable numerical technique with, in principle, arbitrarily high order of accuracy.The SDC method proceeds by first using a simple low order numerical method to compute a provisional (predicted) solution at a set of sub-steps within a given time step.Then, an equation for the correction to this provisional solution is constructed and also approximated on the sub-steps with the same numerical scheme.Each iteration at the correction step increases the accuracy of the provisional solution by one.
Consider the initial value problem, The SDC method is based on the observations concerning the integral form of the solution to Equations ( 12)-( 13) which then can be combined with Equation ( 15) to give This equation is referred to as the correction equation.Now given Equations ( 12) and ( 13) and time interval   is computed by approximating Equation ( 17) to provide an increasingly accurate approximation to the solution To achieve this, the function 17) must be calculated with spectral accuracy, i.e, by a Gaussian quadrature.Our choice of quadrature nodes is the Gauss-Lobatto nodes.This is because the Gauss-Lobatto nodes already include the end points, therefore, one does not have to do extrapolations of the final solution to those points.See [25] for different choices of quadrature nodes.
A discretization for Equation (17) based on the forward Euler method is given by where m Each iteration for Equation (19) raises the formal temporal order of accuracy of the scheme by one, therefore for a fourth order method, three iterations (e.g, k = 3) of the correction equation are needed.Also, for a fourth order method, four Gauss-Lobatto points are necessary, i.e, p = 3 in Equation (18).Now reconsider Equation (11) in the following form where , . Here

The Fourth Order PPM Oriented Flux Calculations
In this section, we define the cell edge averaged quantities that will be used to calculate the numerical flux functions, i.e, in (11) or in (20).The procedure is based on Colella and Woodward [7].In [7], cell edge averages are calculated by interpolating quartic polynomials which are constructed through cell averages.Formulae, in [7], are given for one dimensional non-equally spaced mesh.The formulae below are the two dimensional versions of [7] in a uniformly spaced mesh.Below, we systematically describe the steps to calculation of 1 2, i j  and ler method, i.e, . First, evaluating the interpolation functions ( the construction of interpolation functions and more detail can be found in [7]) at the  1 2, th i j  Correct the provisi ction equation for three iterations, for k = 1, 3  and , 1 2 th i j   edges of cell gives the following edge formulae: 6 where Equations ( 21) and ( 22) are the slope limited edge quantities.Notice that if .
Equations ( 27) and ( 28) are fourth order extrapolations to the edge quantities.Consequently, this mak scheme fourth order accurate in space in regions where th es the PPM e solution is smooth.To prevent possible overshoots or undershoots, the slope limited edge quantities are further monotonized with the monotonicity algorithm.
Monotonicity algorithm: Initialize We note that the aim is to construct the Riemann variables, with which the numerical flux functions are evaluated.To define the point-wise Riemann variables, first we set , .
Then the slope limited and monotonized corner points at each direction can be calculated through following Equ ations ( 21 Now we can define the Riemann variables along cyell edges by taking the two corner values and a cell edge averaged quantity.First we construct a quad nomial along each cell edges.For instance, along the ratic poly- U  , we determine the coefficients of the quadratic polynomial (i.e.,

U y ay by c
structing the quadrat efine o left point-wise Riemann variables at the two Gauss points, i.e.,

After con
ic polynomial, we d the tw for , in , y y .
The right point-wise Riemann variables, R y U can be defined similarly.
Using these Riemann variables, we calc mann solutions at the two Gauss points, i.e. ulate the Rie-, , , Then the Riemann solution is giv The Riemann solution in y-dir milarly.
e employ a words, we define the left and right values in (36) by using the left and right state of the solution vector, e.g., en by  ection is calculated si-Remark: W n up-winding flux procedure at the SDC's prediction step.In other x evaluation for Equations ( 43) and (44).
. This provides more diffusive numerical flu

Solving the Linear Advection Equation
We consider, 0, As a first test, we use a smooth initial solution; In this test, the solution stays smooth i u n time.We monstrate the fourth order of periodic boundary otice sim-) and (28).Also one doesn't have to apply ity algorithm r this specific test case.Figure 2 shows the orm of the absolu di e indicating the stability of the decreases by the factor of en for a finer me in solve this test problem to de accuracy of our method.We apply conditions in both x and y-coordinate directions.N that the edge average formulae for the state variable plifies into (27 the monotonic fo n te error for 2 L sixte fferent meshes.The error grows linearly in tim numerical scheme, and sh indicatg the fourth order convergence.Tables 1-3 show the convergence rates with different norms.The convergence rate  is estimated using the two finest grids according to the following formula

15] and follows
The following test case is conas setting the initial solution to a smooth parabolic profile.For instance, where y = 5x − 2.5 and x belongs to [0, 1].Periodic boundary conditions are applied at both x = 0 and x = 1.At t = 0.2, shock waves start forming in the solution nea x = x = r 0.3 and 0.7.After this time, the shock waves propagate in the solution with To be consistent with [15], we use the same number of grid points (e.g., 40) to produce o sults.The time stepping criteria cept omitting the c, since the only characteristic speed is lution u itself.W ning the code with 200 grid points.Comparing our results to the results reported in [15], we observed that our r better or comparable to the most of the high res schemes (e.g., ENO, WENO, FCT, PPM, TVD-MUSCL urgers equation; ur computational reuses Equation (59) exthe so e set cfl = 1.0 in (59).Figure 4 represents the computed solutions at t = 0.5 and t = 2.4.The reference solutions are calculated by run esults are olution etc).Notice that there is no over/under shoots around the shock profiles.

2-D Burgers Test Problem
We consider the two dimensional B We solve (57 We apply periodic boundary conditions in both coordinate directions.Again the time steps are calculated by (59) with cfl = 0.8.The smooth initial solution turns into a diagonally propagating shock wave around t = 0.05.
Figure 5 shows the time history of the numerical solutions.The solutions are calculated with 80 grid points in each coordinate directions.From Figure 5, we observe that the shock discontinuities are captured sharply without excessive oscillations.The long time run (e.g.
,0 , ,0 0.125,0,0.1 if 0.5 This problem is well studied over the years, i.e., in some cases an exact solution can be constructed.Thus it became a benchmark problem in scientific community to test new numerical schemes.We solve this problem in [0,1] interval with inflow boundary conditions.Our computational results are produced using 400 grid points.The time step is calculated by , max icates the stability of our scheme.Tables the solution is still hen cate the fourth order convergence of our scheme.

4.3.
4-6 are created at t = 0.05 while smooth.We note that the limiting procedure was off w producing the Tables 4-6.These tables clearly indi where c represents the sound speed and cfl = 0.3.Figure 6 shows the numerical solution at t = 0.15. Figure 6 clearly indicates highly resolved discontinuous solutions (e.g, sharp shock and contact discontinuities) with correct shock locations.

Sod's Shock Tube Problem
The shock tube problem considers a long, thin, cylindrical tube containing a gas separated by a thin membrane.The gas is assumed to be at rest on both sides of the membrane, but it has different constant pressures an

Interaction of Shock-Acoustic Waves
This problem is first studied in [8] and describes the interaction of a shock wave with a smooth acoustic entropy wave.The problem solves the one dimensional Euler equations with d densities on each side.At time t = 0, the membrane is ruptured, and the problem is to determine the ensuing motion of the gas.The solution to this problem consists of a shock wave moving into the low pressure region, a rare-faction wave that expands into the high pressure region, and a contact discontinuity which represents the interface.
The shock tube problem of Sod [28] solves the one dimensional Euler equations with the following initial The shock wave interacts with the acoustic wave.As a result, the post-shock acoustic profile is amplified (Figure 7).We solve this problem in [0, 1] interval with inflow boundary conditions.We used 400 grid points and   59) to generate our results.The reference solution is produced with 1600 grid points.Figure 7 shows our computational result.We obtained highly resolved solution comparable to [8] and [12,14].This is an important problem to study, since the solution has some structure in the smooth region.As noted in [8], higher order methods perform better in such regions.Our results confirm this.In other words, our results are comparable to higher order ENO [8], higher order PHM (Piecewise Hyperbolic Methods) [12], or higher order LLR (Local  Logarithmic Reconstructions Method) yet better than the second order MUSCL-type schemes.

Interaction of Blast Waves
This problem is introduced by Woodward and Colella [17] and is often referred to as the Woodward-Colella blast-wave problem.The problem solves the one dimensional Euler equations with  The two discontinuities in the initial data each have the form of a shock-tube problem and yield strong shock waves and contact discontinuities going inwards and rare-faction waves going outwards.The boundary conditions at x = 0 and x = 1 are reflective solid wall conditions, i.e., We used 400 grid points in the [0, 1] interval and cfl = 0.3 in (59) to generate our results.The reference solution is produced with 1600 grid points.Figure 8 shows the numerical results at t = 0.38.It is clear from Figure 8 that we obtained sharp discontinuities without spurious oscillations.

Smooth Euler Solution
We solve this problem to verify the fourth order convergence of our scheme for the Euler equations.The problem consists of an initial value problem with a smooth solution [29].The problem considers the gas initially at rest and 1,0,10 if 0 0.1 1,0,10 if 0.1 0.9 1,0,10 if 0.9 1 , , ,0 , , ,0 1 0.1e ,   7 shows convergence rates for the density field.Clearly, it is fourth order.Sin xact solution for this problem, the convergence analysis , again no limiting procedure is performed.rred Correction (SDC) time integration technique and the Piecewise Parabolic Method (PPM) finite volume method.Our method introduces significant improvements to [1] such as the removal of unphysical oscillations due to flux treatment of these oscillations is that [1] uses the same higher order less diffusive fluxes at all SDC steps (prediction plus deferred correction steps).In order to cure the oscillatory be ce we don't have an e is done similar to the one described in Section 4. We also note that since the solution is smooth

Conclusions
We have presented an improved gas dynamics method based on the Spectral Defe s.The source havior, [1] has to employ couple of external fixes such as the smoothening, flattening, and artificial diffusion algorithms which are not straightforward to implement in a computer code.We avoid all of these extra steps by making use of an up-winding flux procedure at the prediction step.This brings sufficient amount of numerical dif-s t = 0.15. Figure 9.A smooth solution of the Euler equation fusion so that the deferred correction steps can use less diffusive higher order flux evaluations.With this strategy, we have obtained oscillation-free discontinuous solutions with the quality comparable to the original PPM have also demonstrated the fourth order of accurac [7].We y of this method in both space and time for smooth flow problems.

Appendix
Here we define stencils for multiplication and division of the cell averaged (also cell edge averaged in 2-D) quantities in order to compute higher-order accurate numerical flux functions.The primary reason for this i eraged value of a product (or quotient) is not equal to the product (or quotient) of averaged values.The goal here is to express averaged values of quotients or products as the quotient or produ rrection term which must be computed by applying a stencil to neighboring averaged values.
If we were to derive formula for an averaged value of a product, i.e.,   i u  on the th i cell, we first con the following local series expansions sider . 24 Derivation of an averaged value for a quotient and several more details can be found in [20].

Figure 1 .
Figure 1.Cell and cell-edge averaged representation of a quantity  .


then the wave (either ock or rare-factio e) moves from right to left.

Figure 2 .
Figure 2. Demonstrating the fourth order of accuracy of the method with the L 2 -norm when solving Problem (4.1).

Figure 3 .
Figure 3. Propagation of the initial square pulse after one uta-period (t = 1.0).80 grid points are used for this comp tion.

Figure 4 .Figure 5 .
Figure 4. Solution to the one dimensionalBurgers equations t = 0.5 for the left figure and t = 2.4 for the right figure.Numerical solutions are produced with 40 grid points.Solid lines represent reference solutions.,
r  x y z  The solution of this problem

Figure 8 .
Figure 8. Colella's double blast wave problem.t = 0.38.remainsspherically symmetric.Therefore, the three dimensional Euler equations reduce to the one dimensional conservation laws with a source term[29].The problem is solved in a [0, 1] domain.On the left end, the reflective wall boundary conditions are used as in Section 4.3.3.On the right end, the outflow boundary conditions are imposed.The numerical solution in Fig- ure9, produced at t = 0.15 with 160 mesh points, is superimposed on a reference solution produced with 1280 points.Table7shows convergence rates for the density field.Clearly, it is fourth order.Sin xact solution for this problem, the convergence analysis Considering series expansions (71), (72), (73), and (74) for xx  terms in (76) and carrying out the similar algebra, we arrive the following fourth order expression )-(35).At this point we have the right/left corner values, i.e, ,  .Thus we have to consider the following four cases in both coordinate directions.x-direction:Case 1: If L u and , then the wave (either shock herefore n solution becomes

Table 2 . Error table in terms of norm for the linear ad- vection problem (4.1). The final is set to . Here the solution uses the smooth initial Function (52).
quation2.1.1-DBurgers Test ProblemWe consider the one dimensional Burgers equation;