LAI Wencong, KHAN Abdul A.

Glenn Department of Civil Engineering, Clemson University, Clemson, S.C., USA, E-mail: wlai@g.clemson.edu

Time stepping in discontinuous Galerkin method*

LAI Wencong, KHAN Abdul A.

Glenn Department of Civil Engineering, Clemson University, Clemson, S.C., USA, E-mail: wlai@g.clemson.edu

(Received March 7, 2013, Revised April 15, 2013)

The time discretization in the Discontinuous Galerkin (DG) scheme has been traditionally based on the Total Variation Diminishing (TVD) second-order Runge-Kutta (RK2) scheme. Computational efficiency and accuracy with the Euler Forward (EF) and the TVD second-order RK2 time stepping schemes in the DG method are investigated in this work. Numerical tests are conducted with the scalar Burgers equation, 1-D and 2-D shallow water flow equations. The maximum Courant number or time step size required for stability for the EF scheme and RK2 scheme with different slope limiters are compared. Numerical results show that the slope limiters affect the stability requirement in the DG method. The RK2 scheme is generally more diffusive than the EF scheme, and the RK2 scheme allows larger time step sizes. The EF scheme is found to be more efficient and accurate than the RK2 scheme in the DG method in computation.

Discontinuous Galerkin (DG) method, Euler Forward (EF) scheme, second-order Runge-Kutta scheme

Introduction

In recent decades, the Discontinuous Galerkin (DG) method has gained popularity for solving numerically hyperbolic partial differential equations. The DG method is a combination of finite element method and finite volume method.The DG method was first introduced by Reed and Hill for the solution of the steady state neutron transport equation. Early numerical models with the DG method using the explicit time stepping scheme were unconditionally unstable unless a very restrictive time step was utilized. To overcome this drawback, slope limiter was introduced to obtain a Total Variation Bounded (TVB) scheme. The instability problem was further overcome by the introduction of the Runge-Kutta Discontinuous Galerkin (RKDG) method, where the Total Variation Diminishing (TVD) Runge-Kutta method and improved slope limiter were combined. The RKDG method was further developed in a series of papers. Since then, the RKDG method has been widely used in numerical simulation[1,2]. In addition to the popular TVD Runge-Kutta time discretization scheme, Qiu et al.[3]proposed an explicit onestep Lax-Wendroff time discretization in the DG method, and Gassner et al.[4]proposed an explicit onestep time discretization based on a predictor-corrector space-time Taylor expansion formulation. Local timestepping in the DG method has been also proposed[5,6]. The limiting procedures in the DG method brought about diffusive effects in the numerical solutions to depress oscillation[7,8]. Therefore, the use of the slope limiters in every intermediate stage of the Runge-Kutta time discretization will affect the accuracy of the solution. In this paper, the efficiency and accuracy of the first-order Euler Forward (EF) and second-order TVD Runge-Kutta (RK2) time stepping schemes for the DG method are investigated and compared. The two time-stepping schemes with different slope limiters are evaluated through the numerical tests involving 1-D scalar equation, 1-D and 2-D shallow water flow equations. The organization of the paper is as follows. The DG procedure is first outlined and discussed. Extensive numerical tests are presented to compare the DG method with different time stepping schemes and slope limiters.

1. DG procedure

The general procedure of the discontinuous Galerkin method for the hyperbolic conservation laws is briefly described in this section. A system ofm partial differential equations can be written in the conservation form as given in Eq.(1), where U is thevector of conserved variables,F is the flux vector, andS is the vector of source terms and they are defined in Eq.(2).

The problem domain Ωis divided into a collection of Nenon-overlapping elements as given in Eq.(3). The discontinuous Galerkin formulation is then applied element-wise. Inside an element, the variation of the conserved variableU can be approximated using the Lagrange polynomial, and the dependent variableFandS become the functions of the approximated variables. The approximations are given in Eq.(4), where Njis the basis or interpolating function.

The conservation laws are multiplied with the test function Ni, and the resulting equation is integrated over an element as given in Eq.(5). Substitute the approximation ofU and applying the divergence theorem results given in Eq.(6). Then the flux term (F(?)·n)at the boundary is replaced by the numerical flux)as shown in Eq.(7). In this paper, the HLL approximate Riemann solver is used to calculate the numerical flux for both 1-D and 2-D problems.

Equation (7) is integrated and the results can be written in a compact form as given in Eq.(8), where the matrixM and vectorR are given by Eq.(9). Finally, the solution of the unknown conserved variables U can be obtained with appropriate time integration method. In the hyperbolic system of equations, discontinuities and shocks may form even with smooth initial and boundary conditions. Numerical tests show that slope limiters are required to depress the spurious oscillations generated with higher order space approximation.

In the RKDG method, the TVD Runge-Kutta time integration is used, and the slope limiting procedure is applied in every intermediate step of the time integration procedure. In the TVD Runge-Kutta scheme, the time integration is one order higher than the spatial order. In this work, linear elements are used in both 1-D and 2-D cases, and the effects of the firstorder Euler forward and second-order Runge-Kutta scheme are investigated. The first-order EF and second-order Runge-Kutta (RK2) schemes are given by Eqs.(10) and (11), respectively. In the Euler forward method, the slope limiter is applied on Un+1after every time step, while in the TVD Runge-Kutta method, the slope limiter is applied on both U(1)and Un+1in every time step.

The slope limiter for the 1-D linear element is given in Eq.(12), whereis the average of the variable U( x)in elemente,Uel(x)is the limited value within the elemente,is the midpoint of the element,xeand xeare the start and end coordinates of se the element, and σeis the limited slope in an element e. The upwind slope a, downwind slope b, and central slope (a+ b )/2are given by Eqs.(13)-(15), respectively. Three TVD slope limiters[9]are used here, including the monotonized central (MC) slope limiter, the Minmod slope limiter and the Godunov method zero slope limiter, given respectively in Eqs.(16)-(18).

Table 1 Maximum Courant number in the shock wave and rarefaction wave

2. Numerical tests

2.1 Test 1: 1-D scalar equation

The 1-D Burgers equation written in its conservation form is given in Eq.(25). The computational domain is x=[-1,1]and the initial conditions corresponding to the shock and rarefaction waves are given in Eqs.(26) and (27) respectively. The domain is divided into 200 uniform size elements. The maximum Courant numbers (C r= uΔ t/Δx)that give stable results for both shock wave and rarefaction wave with the Euler forward and second-order TVD Runge-Kutta schemes are listed in Table 1.

Numerical results (Table 1) show that the stability requirements for the RK2 scheme are less restrictive than the EF scheme except for the rarefaction wave with zero slope limiter. However, the RK2 scheme requires twice the computational effort compared to the EF scheme in each time step. For the RK2 scheme, the maximum Courant numbers are not twice as large as those for the EF scheme. Thus, the RK2 scheme is not advantageous than the EF scheme in computational efficiency. The test shows that the maximum Courant number is slope limiter dependent. The MC limiter is the least efficient, while the zero limiter has the highest efficiency for the two time-marching schemes. The MC limiter provides the biggest difference between the Courant number based on the EF scheme and the RK2 scheme. The Courant number difference decreases with the use of the Minmod limiter and the least difference is achieved using the zero limiter.

Fig.1 Shock wave at t =0.6 s with MC slope limiter

Fig.2 Shock wave at t =0.6 s with Minmod slope limiter

Numerical results for the shock wave case with the Courant number of 0.5 using the three different slope limiters are shown in Figs.1-3. Numerical results show that the MC limiter is the least diffusive, while the zero limiter is the most diffusive. For each limiter, the RK2 scheme is more diffusive than the EF scheme, because the slope limiter is applied twice with the two-stage time stepping scheme. Further investigation show that the difference between the results of the RK2 and EF schemes reduces as the time step decreases. Numerical results for the rarefaction wave case with the Courant number of 0.5 using the RK2 scheme with the three different slope limiters are shown in Fig.4. The results using the MC slope limiter with large Courant number are poor, showing the results are slope limiter dependent. The MC limiter not only requires smaller time step to give stable results, it also need smaller time step to provide accurate results.This limitation is restricted to RK2 method only. For the EF method, the MC limiter provides stable and accurate results. Numerical results for the rarefaction wave are similar to the shock wave results. The results for the rarefaction wave for the Courant number of 0.1 are shown in Fig.5 with the MC slope limiters only. The MC limiter is the least diffusive, and the zero limiter is the most diffusive. The RK2 scheme is more diffusive than the EF scheme except for the zero slope limiter with small time step, for which the performance of the two time-stepping schemes are similar.

Fig.3 Shock wave at t =0.6 s with zero slope limiter

Fig.4 Rarefaction wave at t =0.6 s using Cr =0.5 with different slope limiters

Fig.5 Rarefaction wave at t =0.6 s using Cr =0.1 with MC slope limiter

2.2 Test 2: 1-D shallow water flows

The 1-D shallow water flows equations for arbitrary channel shape are given by Eq.(28), whereA is the flow area,Q is the flow rate,gis the gravitational acceleration,Zis the water surface elevation, andn is the Manning roughness coefficient. The friction slope Sfis given in Eq.(30), whereRis the hydraulic radius of the flow section. The details of the application of the DG scheme to the 1-D shallow water flow equations, including the wetting and drying scheme, are given by Lai and Khan[11]. The hydrostatic pressure term and the gravity force term are combined and treated as a source term. In order to have meaningful comparison, the water surface gradient term is treated in the same manner for both time stepping schemes.

The idealized dam-break problem with both wet bed and dry bed downstreams of the dam is simulated first, followed by the dry bed dam-break flow in a channel with friction. For the rectangular channel tests, the shallow water flow equations in terms of unit width discharge and flow depth are used[12]. For the idealized dam-break problem, the rectangular, horizontal channel is 1 000 m long with the dam located at 500 m, and the water depth upstream is 10 m. The depths downstream of the dam are 1 m and zero, respectively for the wet bed and dry bed cases. For the dam-break test with friction, the channel is 20 m long with the dam location at 10 m and dry bed downstream. Water depth upstream of the dam is 0.074 m. In all cases, the computational domain is discretized with 200 uniform size elements.

The maximum Courant numbers for the 1-D dam-break problems are listed in Table 2. The results for the dam-break problems are similar to the Burgers equation shock wave. The maximum time step size for the RK2 scheme is larger than the EF scheme for stable solutions. The EF scheme is computationally more efficient than the RK2 scheme except for the idealized dry bed case with minmod slope limiter, where the two schemes have similar efficiency. In the wet bed dam-break case and the dry bed dam-break case with friction, the maximum Courant numbers are similar. In addition, with the zero limiter, the EF scheme is almost twice as efficient as the RK2 scheme.

Table 2 Maximum Courant number in 1-D dam-break problems

Fig.6 Computed water depth in idealized wet bed dam-break with MC limiter

Fig.7 Computed flow rate in idealized wet bed dam-break with MC limiter

For the idealized wet bed dam-break problem,the numerical results at 30 s after dam-break with the Courant number of 0.1 are shown in Figs.6 and 7 for the MC slope limiters. For all three limiters, the EF scheme is more accurate than the RK2 scheme, as it is less diffusive. Out of the three slope limiters, the MC limiter is the least diffusive and the zero limiter is the most diffusive, similar to the results for the Burgers equation shock wave. In addition, the difference in results is the least using the two time-stepping schemes for the zero limiter. Since the MC limiter provides the most accurate results, only the results with the MC limiter are presented for the dry bed case. For the dry bed idealized dam-break, numerical result at 20 s is shown in Fig.8 with the MC limiter and the Courant number of 0.1. The EF scheme is more accurate in capturing the moving wave.

Fig.8 Computed water depth in idealized dry bed dam-break with MC limiter

Fig.9 Computed water depth in dry bed dam-break with friction using MC limiter

For the dam-break flow with friction, the measured water surface elevations after the dam-break are available. The rectangular, horizontal flume is 0.096 m wide, 0.08 m high. The Manning roughness coefficient is n=0.009 s/m1/3. The computed water surfaces are compared with the measured data at 3.75 s and 9.40 s after the removal of the dam in Fig.9. The Courant number of 0.2 is used in the simulation. The results show that the two time-stepping schemes provide similar results. A small difference in the rarefaction wave speed with the EF scheme is highlighted in the figure.

Table 3 Maximum time step (s) in the Toce River case

Fig.10 Simulated and measured hydrograph in Toce River (1-D test)

Fig.11 Computational domain and mesh for the symmetric channel contraction

To investigate the performance of the two timestepping schemes for the real world problems involving arbitrary channel shape and varying bed topography. The Toce River dam-break flow is investigated. The Toce River physical model was established at the ENEL-HYDRO Laboratory in Milan at a scale of 1:100. The physical model covered an area of 50 m× 12 m. The river was initially dry. Details of the modeling parameters and discussion about the governing equations can be found in Refs.[11,13]. In this case, 61 elements of non-uniform size are used to accommodate all the measure cross sections. The maximum allowable time steps for the three slope limiters are listed in Table 3. The maximum time step allowed for the EF scheme is slightly smaller than RK2, which issimilar to the dry bed dam-break flow with friction. This shows that the EF scheme is twice as efficient as the RK2 scheme. Computed (with the MC limiter) and measured hydrographs at different gauge points are presented in Fig.10. The two time-stepping schemes provide similar results, however, the EF scheme is clearly more efficient.

Table 4 Maximum time step (s)

2.3 Test 3: 2-D shallow water flows

The 2-D depth-averaged, shallow water flow equations are given by Eq.(30). The friction slopes in thexandy directions are given by Eq.(31). In these equations,h is the flow depth,qxand qyare the unit width flow rates in thexandydirections, respectively, and zbis the bed elevation, and nis the Manning roughness coefficient. To evaluate the two time-stepping schemes, several simulations are performed including steady state hydraulic jump in a channel contraction, dam-break flow in a horizontal channel with a 90obend, and dam-break in a natural channel. The application of discontinuous Galerkin method to 2-D shallow water flow equation is given by Lai and Khan[2,10]. In case of natural channels with varying bed topography, the pressure force term and the bed elevation term are combined into a water surface elevation term as in the cased of 1-D shallow water flows. The water surface gradient terms are treated in the same way for the two time stepping schemes. The discretization details are given by Lai and Khan[10].

For the steady hydraulic jump in channel contraction case, the plan view of geometry and computational mesh are shown in Fig.11. Lin et al.[14]provided detail description of the shock wave in channel contraction. The channel is 20 m wide at the upstream end and 10.548 m at the downstream end. The length of the contraction is 22.234 m, the angle of wall deflection is 12o. The water depth at the inflow boundary is 1 m, the longitudinal velocity is 8.4566 m/s, and the lateral velocity is zero. The inflow is supercritical with the Froude number of 2.7. In this steady hydraulic jump case, the maximum time step that give stable results are the same for the two different slopes with two different time stepping schemes as shown in Table 4. The EF time-stepping scheme is computationally more efficient for this case. The two slope limiters and the two time-stepping schemes provide similar numerical results. The computed water surface using the EF scheme and linear slope limiter are plotted in Fig.12. The computed and exact solutions along the dash line and the solid line (Fig.11) are shown in Figs.13 and 14 respectively.

Fig.12 Computed water surface profile in 3-D view

Fig.13 Computed and exact solutions along the dash line

Fig.14 Computed and exact solutions along the solid line

Fig.15 Computational domain and mesh for the channel with bend

Table 5 Gauge point locations

The geometry and computational domain for the horizontal channel with a 90obend is shown in Fig.15 and the gauge points locations are listed in Table 5. The dam is represented by the gate at the reservoir outlet. The initial water depth in the reservoir is 0.2 m and the downstream channel is dry. The Manning roughness coefficients are 0.0095 s/m1/3and 0.0195 s/m1/3for the bottom and wall respectively[15]. The computed hydrographs using the EF scheme at the gauge points G1 and G3 are compared with the measured data in Figs.16 and 17. Time step size of 0.0013 s is used for both zero slope limiter and linear slope limiter cases. As was expected, the zero slope limiter is more diffusive than the linear slope, especially at G1. For the linear slope limiter, the maximum time step size (Table 4) for the EF scheme is half that of the RK2 scheme and the two time-stepping schemes are equally efficient. However, for the zero slope limiter the EF scheme is more efficient.

Fig.16 Computed and measured hydrographs at G1

Fig.17 Computed and measured hydrographs at G3

Fig.18 Simulated and measured hydrograph in Toce River (2-D test)

The Toce River case, simulated earlier with the 1-D shallow water flow equations, is revisited with the 2-D shallow water flow equations. For the 2-D case, 13 316 triangular elements are used. The results are slope limiter dependent. For the zero slope limiter, the maximum time step in the EF is smaller than the RK2 as expected, while for the linear slope limiter, the maximum allowed time step in EF is slightly larger than RK2. Overall, the EF scheme is more efficient than the RK2 scheme. Numerical results of computed and measured hydrographs at different gauge pointsare presented in Fig.18. The linear slope limiter with the time step 0.002 s is used for comparison. The EF and RK2 schemes give similar accuracy.

3. Conclusions

The EF scheme and second-order total variation diminishing Runge-Kutta (RK2) scheme are investigated for time integration within the DG method have been investigated and discussed. The 1-D scalar equation, 1-D shallow water flow equations and 2-D shallow water flow equations are used to evaluate the two time-stepping schemes. The numerical tests include shock wave and rarefaction wave, steady and unsteady flows, wet bed and dry bed problems, idealized and real world channels. The effects of different slope limiters on time-stepping schemes in the discontinuous Galerkin method have also been investigated. The maximum Courant number or maximum time step required for stability in each test is evaluated with different time stepping schemes and slope limiters. Numerical results show that the EF scheme is computationally more efficient and accurate than the RK2 scheme.The RK2 scheme is more diffusive than the EF scheme, as the slope limiters are applied twice in the RK2 scheme. Even though the RK2 scheme can generally have larger time step size, and the two-step process makes it computationally less efficient than the EF scheme. In cases where bed friction is included, the EF and RK2 schemes give similar performance in accuracy. In addition, the more diffusive slope limiters (e.g., zero slope limiter) can have less restrictive Courant number (larger time step size). It would be important to generate slope limiters that have less diffusive effect and less restrictive Courant number stability requirement in future research.

[1] CHEN Er-yun, MA Da-wei and LE Gui-gao et al. Numerical simulation of highly underexpanded axisymmetric jet with Runge-Kutta discontinuous Galerkin finite element method[J]. Journal of Hydrodynamics, 2008, 20(5): 617-623.

[2] LAI Wencong, KHAN Abdul A. Modeling dam-break flood over natural rivers using discontinuous Galerkin method[J]. Journal of Hydrodynamics, 2012, 24(4): 467-478.

[3] QIU J., DUMBSER M. and SHU C. W. The discontinuous Galerkin method with Lax-Wendroff type time discretizations[J]. Computer Methods in Applied Mechanics and Engineering, 2005, 194(42-47): 4528- 4543.

[4] GASSNER G., DUMBSER M. and HINDENLANG F. et al. Explicit one-step time discretizations for discontinuous Galerkin and finite volume schemes based on local predictors[J]. Journal of Computational Physics, 2011, 230(11): 4232-4247.

[5] TRAHAN C. J., DAWSON C. Local time-stepping in Runge-Kutta discontinuous finite element methods applied to shallow-water equations[J]. Computer Methods in Applied Mechanics and Engineering, 2012, 217- 220: 139-152.

[6] DAWSON C., TRAHAN C. J. and KUBATKO E. J. et al. A parallel local timestepping Runge-Kutta discontinuous Galerkin method with applications to coastal ocean modeling[J]. Computer Methods in Applied Mechanics and Engineering, 2013,259(1): 154-165.

[7] ZHONG X., SHU C. W. A simple weighted essentially nonoscillatory limiter for Runge-Kutta discontinuous Galerkin methods[J]. Journal of Computational Physics, 2013, 232(1): 397-415.

[8] ZHAO J., TANG H. Runge-Kutta discontinuous Galerkin methods with WENO limiter for the special relativistic hydrodynamics[J]. Journal of Computa- tional Physics, 2013, 242(1): 138-168.

[9] LI B. Q. Discontinuous finite elements in fluid dynamics and heat transfer[M]. London, UK: Springer- Verlag, 2006.

[10] LAI W., KHAN A. A. A discontinuous Galerkin method for two-dimensional shallow water flows[J]. International Journal for Numerical Methods in Fluids, 2012, 70(8): 939-960.

[11] LAI W., KHAN A. A. Discontinuous Galerkin method for 1D shallow water flow in nonrectangular and nonprismatic channels[J]. Journal of Hydraulic Enginee- ring, ASCE, 2012, 138(3): 285-296.

[12] LAI W., KHAN A. A. Discontinuous Galerkin method for 1D shallow water surface flow with water surface slope limiter[J]. International Journal of Civil and Environmental Engineering, 2011, 3(3): 167-176.

[13] LAI W., KHAN A. A. Discontinuous Galerkin method for 1D shallow water flows in natural rivers[J]. Engineering Applications of Computational Fluid Me- chanics, 2012, 6(1): 74-86.

[14] LIN G., LAI J. and GUO W. High-resolution TVD schemes in finite volume method for hydraulic shock wave modeling[J]. Journal of Hydraulic Research, 2005, 43(4): 376-389.

[15] SOAREoS FRAZ?O S., ZECH Y. Dam-break in channel with 90 bend[J]. Journal of Hydraulic Engineering, ASCE, 2002, 128(11): 956-968.

10.1016/S1001-6058(11)60370-4

* Biography: LAI Wencong (1985-), Male, Ph. D.