High Order Accuracy Research Papers

Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws

1998, Advanced numerical approximation of nonlinear …

In these lecture notes we describe the constraction, analysis, and applica-tion of ENO (Essentially Non-Oscillatory) and WENO (Weighted Essen-tially Non-Oscillatory) schemes for hyperbolic conservation laws and related Hamilton-Jacobi... more

Fic. 2.2. First order (top) and fourth order (bottom) ENO schemes for the Burgers equation, with the Godunov flux (left, and the Lax-Friedrichs flux (right). Solid lines: exact solution; Circles: the computed solution at t = 4.

is the cell average. We approximate (3.22) by the following conservative scheme describe the ENO and WENO schemes for 2D conservation laws:

Although there are many zero terms in the inner sum of (2.20) when x;, 1 is a node in the interpolation, we will keep this general form so that it applies also to the case where 2, 2 is not an interpolation point.

For easier manipulation we subtract a constant V(z,;_,_1) from (2.17), and use the fact that :— Evaluating the expression (2.19) at « = 2,,1, we finally obtain

Fic. 2.1. Fixed central stencil cubic interpolation (left) and ENO cubic interpolation (right) for the step function. Solid: exact function; Dashed: interpolant piecewise cubic polynomials.

with a CFL coefficient c = 0.936 in (4.14). Notice that two L’s must be computed. The effective CFL coefficient, comparing with an ideal case without L’s, is 0.936 x 4 = 0.624. Since it is difficult to solve the global optimization problem, we do not claim that (4.15) is the optimal 4 stage, 4th order TVD Runge-Kutta method.

Fic. 4.1. Second order TVD MUSCL spatial discretization. Solution after 500 time steps. Left: TVD time discretization (4.10); Right: non-TVD time discretization (4.22).

Fic. 5.1. Shock tube, Lax problem, density. (a): third order WENO; (b): fifth order WENO; (c): fifth order WENO with cheaper characteristic decomposition; (d): fifth order WENO with artificial compression. applied to these cases [43]. We show the results of WENO (third order and fifth order) schemes for the Lax problem, in Fig. 5.1. Notice that “PS” in the pictures means a way of treating the system cheaper than the local characteristic decompositions (for details, see [43]). “A” stands for Yang’s artificial compression [83] applied to these cases [43].

Fic. 5.3. 2D shock entropy wave interaction. Amplitude of amplified entropy waves. 800 points (about 20 points per entropy wave length).

Fic. 5.4. Density. Steady quasi-1D nozzle flow. 34 points. Left: fourth order WENO; right: fifth order WENO.

Fic. 5.5. Flow past a forward facing step. Density: 242 x 79 grid points. Top: fifth order WENO; bottom: fourth order ARNO.

Fic. 5.6. Double Mach reflection. Density: 480 x 119 grid points. Top: fifth order WENO; bottom: fourth order ENO.

Fic. 5.8. Flow past a cylinder. (a) Physical grid. (b) Pressure. WENO-LF-5. 20 contours

Fic. 5.9. Vortex evolution. Cut at x = 5. Solid: exact solution; circles: computed solution. Top: t = 50 (after 5 time periods); bottom: t = 100 (after 10 time periods). Left: second order TVD scheme; right: fifth order WENO scheme. incompressible Navier-Stokes or Euler equations

We then use the third order TVD (total variation diminishing) Runge-Kutta method (4.11) to discretize the resulting ODE:

Fic. 5.10. Double shear layer. Contours of vorticity. t = 8. Top left: spectral with 512? points; top right: ENO with 64? points; bottom left: ENO with 128? points; bottom right: the cut at x = 7 of v, spectral method with 512? points, ENO method with 647 and with 128? points.

Fic. 5.11. Two dimensional vortex sheet simulation. t = 4. Top ce ENO with 128? points, 6 function width « = 12Ax = 2; Top right: ENO with 256? points, 6 function width « = 12Ax = 3; Bottom left: ENO with 512? points, 6 function width e = 24Ax = 3; Bottom right: ENO with 1024? points, 6 function width « = 48Ax = =. The 6 function is approximated as in [61],[77] by

Fic. 5.12. Three dimensional vortex sheet simulation. t = 5. ENO with 64° points. 5 function width e = 6Ax. Left: three dimensional level surface; Right: z = 0 plane cut. where n is the electron concentration, p = (0,7) is the momenta, and W is the total energy. The equations, in two dimensions, take the form

Fic. 5.13. The one dimensional n+-n-n+ channel. Left: the doping ng; Right: the velocity v, comparison of the HD model and the reduced HD model. where nq is the doping (a given function which is typically discontinuous).

Accuracy on uz + uz =0 with uo(x) = sin4 (rz). TABLE 2.3

Accuracy on ut + Ur = 0 with uo(x) = sin(ra).

3.1. Reconstruction and Approximation in Multi Dimensions. In this section we describe how and assume that Az and Ay are of the same magnitude (their ratio is bounded from above and below during refinement). Finally,

Problem 3.2. Two dimensional conservative approximation.

Prominent examples are vortex sheets and vortex filaments in three dimensions, and vortex sheets, vortex dipole sheets and point vortices in two dimensions. In three dimensions, the equations are written in the form

descriptionView Paper arrow_downwardDownload

Application of the generalized differential quadrature method to the study of pull-in phenomena of MEMS switches

by Ghader Rezazadeh

2007, … Systems, Journal of

This paper reports on the pull-in behavior of nonlinear microelectromechanical coupled systems. The generalized differential quadrature method has been used as a high-order approximation to discretize the governing nonlinear... more

descriptionView Paper arrow_downwardDownload

On maximum-principle-satisfying high order schemes for scalar conservation laws

by Chi-wang Shu

2010, Journal of Computational Physics

Keywords: Hyperbolic conservation laws Finite volume scheme Discontinuous Galerkin method Essentially non-oscillatory scheme Weighted essentially non-oscillatory scheme Maximum principle High order accuracy Strong stability preserving... more

5 Third order SSP time discretization and fifth order finite volume WENO scheme with the limiter, for the Burgers’ equation. t = 0.15, Ax = 2 and At = 4 Axi, Table 5.2

Third order SSP time discretization and fifth order finite volume WENO scheme with the limiter, for the two-dimensional linear equation with initial data sin(2m(x +y)).t=0.1, At =4Ax3 and Ax = Ay =}. Table 5.3 Example 5.3. We use the nonconvex Buckley-Leverett flux

then the coefficients can be given explicitly as In summary, the algorithm for the limiter on the fifth order finite volume WENO with the Euler forward time discretization is:

Remark 4.4. We can prove the same results for the upwind flux defined in [15] following the same lines Remark 4.5. The results of Theorem 4.3 holds also for any passive convection linear equations with divergence-free velocity coefficients, namely Eq. (4.1) in which u and vare given functions satisfying u, + v, = 0, as long as the quadratures are exact for the integrands in the scheme. This can be easily achieved if we pre-process the divergence-free velocity field so that it is piecewise polynomial of the right degree for accuracy, continuous in the normal component across cell boundaries, and divergence-free.

Fig. 5.1. Example 5.2. N = 80, At = LAX. Solid line: exact solution; symbols: numerical solution (cell averages). Left: t = 0.6366; right: t = 2

‘ig. 5.2. Example 5.3. t= 0.4, N= 160, Ax =2 and “a = 4. Solid line: exact solution; symbols: numerical solution (cell averages

Fig. 5.3. Example 5.5. The numerical solution of the two-dimensional Burgers’ equation. t = 0.6. At = LAX and Ax = Ay = Third order SSP time discretization and fifth order finite volume WENO scheme with the limiter, for the two-dimensional Burgers’ equation with the initial data 0.5+sin(z(x+y)).t=0.05, At =1Ax3 and Ax = Ay =2.

Fig. 5.4. Example 5.5. Cut along the diagonal. Solid line: the exact solution; symbols: numerical solution. Left: t = 0.23, right: t = 0.6. Ax = Ay = 2.

Fig. 5.5. Example 5.6 (traffic flow problem). N = 800, Ax = *? and #4a = 4 where a = 100. Solid line: exact solution; circles: numerical solution (cell averages). 5.2. Standard two-dimensional test cases

Fig. 5.6. Left: the graph of q(p); right: the graph of q’(p).

Fig. 5.7. Example 5.7 (traffic flow problem). t = 18 min, N = 1600, Ax = 4? and 4a = ;4, where a ~ 100. Solid line: exact solution; circles: numerical solution cell averages). The region inside the rectangle on the right is magnified in Fig. 5.8.

Fig. 5.8. Magnified graph of the region inside the rectangle in Fig. 5.7.

Fig. 6.1. Example 6.1. N= 160, t= 100, Ax =2 and At = 4. Ax. Solid line: exact solution; squares: numerical solution.

Incompressible Euler equations. P? for vorticity, t = 0.5.

Fig. 6.3. Vorticity at t = 8, P’, 30 equally spaced contours from —4.9 to 4.9.

Fig. 6.4. Vorticity at t = 5, P’. Top: 30 equally spaced contours from —1.1 to 1.1; bottom: cut along the diagonal. comparison, we also show the result of the same DG scheme without the limiter. The overshoot near the discontinuity is then apparent.

Fig. 6.5. Vorticity at t = 10, P’. Top: 30 equally spaced contours from —1.1 to 1.1; bottom: cut along the diagonal.

The CFL number (2.6) for 2 < k <5 and the Gauss-Lobatto quadrature points on [—4,}). Table 2.1

Third order SSP time discretization and fifth order finite volume WENO scheme with the limiter, for the linear equation with initial data uo(x). Ax = 4 and At = 1 Ax3. Table 5.1

Third order SSP time discretization and the DG scheme with the limiter, for the linear equation with initial data sin(27x), Ax =}, t=0.1 Table 6.1 with suitable initial and boundary conditions. Here p € (0, p,,,,) is the density, p,,,, is the maximum (jam) density, and q(p) is the traffic flow on a homogeneous highway, which is assumed to be a function of the density only in the LWR model. The flow q, the density p and the equilibrium speed u are related by

Third order SSP time discretization and the DG scheme with the limiter, for the Burgers’ equation with initial data 0.25 + 0.5 sin(x), Ax=2, t=0.15 Table 6.2

Third order SSP time discretization and the DG scheme with the limiter, for the two-dimensional linear equation with initial data sin(7z(x +y)), Ax =, t= 0.1. 2IN, Table 6.3 Table 6.4

Third order SSP time discretization and the DG scheme with the limiter, for the two-dimensional Burgers’ equation with initial data sin(z(x+y)), Ax =}, t = 0.05.

descriptionView Paper arrow_downwardDownload

Runge–Kutta discontinuous Galerkin method using a new type of WENO limiters on unstructured meshes

by Chi-wang Shu and

2013, Journal of Computational Physics

In this paper we generalize a new type of limiters based on the weighted essentially non-oscillatory (WENO) finite volume methodology for the Runge-Kutta discontinuous Galerkin (RKDG) methods solving nonlinear hyperbolic conservation... more

descriptionView Paper arrow_downwardDownload

Central WENO schemes for hyperbolic systems of conservation laws

by Giovanni Russo

1999, Mathematical Modelling and Numerical Analysis

We present a family of high-order, essentially non-oscillatory, central schemes for approximating solutions of hyperbolic systems of conservation laws. These schemes are based on a new centered version of the Weighed Essentially... more

descriptionView Paper arrow_downwardDownload

High order finite difference WENO schemes with the exact conservation property for the shallow water equations

by Chi-wang Shu

2005, Journal of Computational Physics

Shallow water equations with nonflat bottom have steady state solutions in which the flux gradients are nonzero but exactly balanced by the source term. It is a challenge to design genuinely high order accurate numerical schemes which... more

descriptionView Paper arrow_downwardDownload

Numerical methods of high-order accuracy for nonlinear boundary value problems

by Martin Schultz

1968, Numerische Mathematik

descriptionView Paper arrow_downwardDownload

A locally conservative LDG method for the incompressible Navier-Stokes equations

by Bernardo Cockburn

2004, Mathematics of Computation

In this paper a new local discontinuous Galerkin method for the incompressible stationary Navier-Stokes equations is proposed and analyzed. Four important features render this method unique: its stability, its local conservativity, its... more

descriptionView Paper arrow_downwardDownload

High-order well-balanced finite volume WENO schemes for shallow water equation with moving water

by Chi-wang Shu and

2007, Journal of Computational Physics

Hyperbolic balance laws have steady state solutions in which the flux gradients are nonzero but are exactly balanced by the source term. In our earlier work [J. Comput. Phys. 208 (2005) 206-227; J. Sci. Comput., accepted], we designed a... more

Fig. 21. Comparison of the convergence history in L' error. Left: FV WENO schemes; right: RKDG schemes.

We integrate (2.9) over the interval J;; to obtain: The monotone flux F(a,b) is the same as defined above (for example, Formula (2.3)). We 9 4 p,Ay, are the (2k — 1)th order accurate reconstructed values obtained by a WENO reconstruction procedure. In this proce- dure, for rectangular meshes, if we use the tensor products of one dimensional polynomials, i.e., polynomials in QX~', things can proceed as in one dimension. A practical way to perform the reconstruction in two dimen- sions is given as follows. We first perform a one dimensional reconstruction in one of the directions (e.g., the y-direction), obtaining one dimensional cell averages of the function u in the other direction (e.g., the x-direc- tion). We then perform a reconstruction in the other direction to obtain the approximated point values, see 197.901

We need to decompose it further in the following way in order to obtain a well-balanced scheme where (7), 44 is a high order approximation to ¢;(x;,1), whose definition will be described below, and g;; is any high order approximation to the integral

Fig. 1. The initial surface level h + b and the bottom 5 for a small perturbation of a steady state water. Left: a big pulse « = 0.2; right: a small pulse « = 0.001. the RKDG scheme. The results for the small pulse « = 0.001 are shown in Figs. 3 and 5. For this small pulse problem, we take ¢ = 10~? in the WENO weight formula (2.5), such that it is smaller than the square of the perturbation. At this time, the downstream-traveling water pulse has already passed the bump. We can clearly see that there are no spurious numerical oscillations.

Fig. 2. FV scheme: small perturbation of a steady state water with a big pulse. t = 0.2 s. Left: surface level A + 5; right: the discharge hu.

Fig. 4. RKDG scheme: small perturbation of a steady state water with a big pulse. t = 0.2 s. Left: surface level h + 6; right: the dischargs hu. Fig. 3. FV scheme: small perturbation of a steady state water with a small pulse. t = 0.2 s. Left: surface level h + 5; right: the discharge hu.

Fig. 6. FV scheme: The surface level / + b for the dam breaking problem at time t = 15 s. Left: the numerical solution using 400 grid cells plotted with the initial condition and the bottom topography; right: the numerical solution using 400 and 4000 grid cells. Fig. 5. RKDG scheme: small perturbation of a steady state water with a small pulse. ¢ = 0.2 s. Left: surface level h + b; right: the dischargs hu.

Fig. 7. FV scheme: The surface level i + 6 for the dam breaking problem at time ¢ = 60 s. Left: the numerical solution using 400 grid cells, plotted with the initial condition and the bottom topography; right: the numerical solution using 400 and 4000 grid cells.

Fig. 8. RKDG scheme: The surface level 4 + b for the dam breaking problem at time ¢ = 15 s. Left: the numerical solution using 400 grid cells, plotted with the initial condition and the bottom topography; right: the numerical solution using 400 and 4000 grid cells.

Fig. 9. RKDG scheme: The surface level / + b for the dam breaking problem at time ¢ = 60 s. Left: the numerical solution using 400 gric cells, plotted with the initial condition and the bottom topography; right: the numerical solution using 400 and 4000 grid cells.

Fig. 10. FV scheme: steady transcritical flow over a bump without a shock. Left: the surface level 4+ b; right: the discharge hu as th numerical flux for the water height h.

Fig. 11. RKDG scheme: steady transcritical flow over a bump without a shock. Left: the surface level + b; right: the discharge hu as the numerical flux for the water height h.

Fig. 12. FV scheme: steady transcritical flow over a bump with a shock. Left: the surface level h+ 5; right: the discharge hu as the numerical flux for the water height h. 6.2. Two dimensional shallow water equations

Fig. 13. FV scheme: steady transcritical flow over a bump with a shock. Pointwise error comparison between numerical solutions using 200 and 400 cells. Left: the surface level h + b; right: the discharge hu as the numerical flux for the water height h.

Fig. 14. RKDG scheme: steady transcritical flow over a bump with a shock. Left: the surface level h + b; right: the discharge hu as th numerical flux for the water height h.

Fig. 15. RKDG scheme: steady transcritical flow over a bump with a shock. Pointwise error comparison between numerical solutions using 200 and 400 cells. Left: the surface level + b; right: the discharge hu as the numerical flux for the water height h.

Fig. 16. FV scheme: steady subcritical flow over a bump. Left: the surface level 4 + 5; right: the discharge hu as the numerical flux for the water height h.

Fig. 17. RKDG scheme: steady subcritical flow over a bump. Left: the surface level h + 5; right: the discharge hu as the numerical flux fo: the water height h.

Fig. 18. FV scheme: The contours of the surface level / + 6 for the problem in Section 6.2.3. 30 uniformly spaced contour lines. From top to bottom: at time t = 0.12 from 0.99942 to 1.00656; at time t = 0.24 from 0.99318 to 1.01659; at time ¢ = 0.36 from 0.98814 to 1.01161; at time ¢ = 0.48 from 0.99023 to 1.00508; and at time t = 0.6 from 0.99514 to 1.00629. Left: results with a 200 x 100 uniform mesh. Right: results with a 600 x 300 uniform mesh.

Fig. 19. RKDG scheme: The contours of the surface level + 6 for the problem in Section 6.2.3. 30 uniformly spaced contour lines. From top to bottom: at time t= 0.12 from 0.99942 to 1.00656; at time t= 0.24 from 0.99318 to 1.01659; at time t= 0.36 from 0.98814 to 1.01161; at time t= 0.48 from 0.99023 to 1.00508; and at time t= 0.6 from 0.99514 to 1.00629. Left: results with a 200 x 100 uniform mesh. Right: results with a 600 x 300 uniform mesh.

Fig. 20. The numerical (symbols) and the “exact” reference (solid line) stress o(x) at time t=0.4s. Left: FV schemes; right: RKDG schemes.

L' and L® errors for different precisions for the stationary solution with a smooth bottom (6.2

L! errors and numerical orders of accuracy for the example in Section 6.1.2 Table 3 where € is a non-zero perturbation constant. Two cases have been run: ¢ = 0.2 (big pulse) and « = 0.001 (small pulse). Theoretically, for small ¢, this disturbance should split into two waves, propagating left and right at the characteristic speeds +,/gh. Many numerical methods have difficulty with the calculations involving such small perturbations of the water surface [21]. Both sets of initial conditions are shown in Fig. |. The solution at time ¢ = 0.2 s for the big pulse « = 0.2, obtained on a 200 cell uniform grid with simple transmissive bound- ary conditions, and compared with a 3000 cell solution, is shown in Fig. 2 for the FV scheme and in Fig. 4 for

L' and L™ errors for different precisions for the stationary solution with a nonsmooth bottom (6.3) Table 2 with periodic boundary conditions, see [35]. Since the exact solution is not known explicitly for this case, we use the fifth order finite volume WENO scheme with N = 12,800 cells to compute a reference solution, and treat this reference solution as the exact solution in computing the numerical errors. We compute up to t= 0.1 when the solution is still smooth (shocks develop later in time for this problem). Table 3 contains the L' errors for the cell averages and numerical orders of accuracy for the finite volume and RKDG schemes, respectively. We can clearly see that fifth order accuracy is achieved for the WENO scheme, and third order accuracy is achieved for the RKDG scheme. For the RKDG scheme, the TVB constant M is taken as 32. Notice that the CFL number we have used for the finite volume scheme decreases with the mesh size and is recorded in Table 3. For the RKDG method, the CFL number is fixed at 0.18.

L' errors for different precisions for the stationary solution in Section 6.2.1 Table 4 defined over a unit square, with periodic boundary conditions. The terminal time is taken as t = 0.05 to avoid the appearance of shocks in the solution. Since the exact solution is also not known explicitly for this case, we use the same fifth order WENO scheme with an extremely refined mesh consisting of 1600 x 1600 cells to

RKDG scheme: L! errors and numerical orders of accuracy for the example in Section 6.2.2 Table 6

FV scheme: L! errors and numerical orders of accuracy for the example in Section 6.2.2 Table 5

L' errors and numerical orders of accuracy for the example in Section 6.4 Table 7 satisfying our requirement. For the RKDG scheme, we define

L' and L® errors for different precisions for the steady state (6.17)-(6.18)

descriptionView Paper arrow_downwardDownload

Resolution of high order WENO schemes for complicated flow structures

by Chi-wang Shu and

2003, Journal of Computational Physics

In this short note we address the issue of numerical resolution and efficiency of high order weighted essentially nonoscillatory (WENO) schemes for computing solutions containing both discontinuities and complex solution features, through... more

Fig. 1. Double Mach reflection problem. Blown-up region around the double Mach stems. Dynsity ps 30 equally spaced contour lines from p = 1.5 to p = 22.9705. Lek from top to bottom: fifth order WENO results withh =5L, ., |; right from top to bottom: ninth an aap) 480° sea order WENO results with 4 = a aan. 560? J. Shi et al. | Journal of Computational Physics 186 (2003) 690-696

For these two test cases with both discontinuities and complex solution structure, it seems that we can get comparable numerical resolution hence comparably small numerical viscosity using WENO9 with half Fig. 2. Rayleigh-Taylor Instability. Density p; 15 equally spaced contour lines from p = 0.952269 to p = 2.14589. Top from left to right: fifth order WENO results with h=35, gg sp qa bottom from left to right: ninth order WENO results with oe — 540? 480? 960°? 1920°

Parallel CPU time for WENOS and WENOS, unit: second, 48 processors, Rayleigh-Taylor instability problem, t = 1.95 Parallel efficiency, ninth order WENO scheme, Rayleigh-Taylor instability problem, t = 1.95

of the mesh points in each direction compared with the results using WENOS. This will more than offset the additional CPU cost per grid point of WENO9 versus WENOS. In Table 1, we list the CPU time for WENOS and WENO9 using 48 processors. We can see that, when comparable resolution is obtained with WENOS9 using half the number of mesh points in each direction than WENOS, the CPU time needed by WENOS is only about 20-34% of that needed by WENOS by our implementation. It is also important to note that this saving is more significant for more refined meshes.

descriptionView Paper arrow_downwardDownload

High Order Fast Sweeping Methods for Static Hamilton–Jacobi Equations

by Yong-Tao Zhang

2006, Journal of Scientific Computing

We construct high order fast sweeping numerical methods for computing viscosity solutions of static Hamilton-Jacobi equations on rectangular grids. These methods combine high order weighted essentially non-oscillatory (WENO)... more

descriptionView Paper arrow_downwardDownload

An Object-Oriented Parallel Particle-in-Cell Code for Beam Dynamics Simulation in Linear Accelerators

by Salman Habib

2000, Journal of Computational Physics

We present an object-oriented three-dimensional parallel particle-in-cell (PIC) code for simulation of beam dynamics in linear accelerators (linacs). An important feature of this code is the use of split-operator methods to integrate... more

descriptionView Paper arrow_downwardDownload

Numerical simulation of low Mach number reacting flows

by Marcus Day

2007, Journal of Physics: Conference Series

Using examples from active research areas in combustion and astrophysics, we demonstrate a computationally efficient numerical approach for simulating multiscale low Mach number reacting flows. The method enables simulations that... more

descriptionView Paper arrow_downwardDownload

A High-Order Nodal Discontinuous Galerkin Method for 1D Morphodynamic Modelling

by Mohammed Seaid

2012, International Journal of Computer Applications

In this paper we present a numerical solution of the sediment transport equations in one horizontal dimension, based on a discontinuous Galerkin finite-element method. The continuous equations are discretized using nodal polynomial basis... more

descriptionView Paper arrow_downwardDownload

A WENO-solver for the transients of Boltzmann–Poisson system for semiconductor devices: performance and comparisons with Monte Carlo methods

by Chi-wang Shu and

2003, Journal of Computational Physics

In this paper we develop a deterministic high order accurate finite-difference WENO solver to the solution of the 1-D Boltzmann-Poisson system for semiconductor devices. We follow the work in Fatemi and Odeh [9] and in Majorana and... more

descriptionView Paper arrow_downwardDownload

Anti-diffusive flux corrections for high order finite difference WENO schemes

by Chi-wang Shu

2005, Journal of Computational Physics

In this paper, we generalize a technique of anti-diffusive flux corrections, recently introduced by Després and Lagoutière [Journal of Scientific Computing 16 (2001) 479-524] for first-order schemes, to high order finite difference... more

descriptionView Paper arrow_downwardDownload

A numerical study for the performance of the Runge–Kutta discontinuous Galerkin method based on different numerical fluxes

by Chi-wang Shu

2006, Journal of Computational Physics

Runge-Kutta discontinuous Galerkin (RKDG) method is a high order finite element method for solving hyperbolic conservation laws employing useful features from high resolution finite volume schemes, such as the exact or approximate Riemann... more

descriptionView Paper arrow_downwardDownload

The discontinuous Galerkin method with Lax–Wendroff type time discretizations

by Chi-wang Shu

2005, Computer Methods in Applied Mechanics and Engineering

In this paper we develop a Lax-Wendroff time discretization procedure for the discontinuous Galerkin method (LWDG) to solve hyperbolic conservation laws. This is an alternative method for time discretization to the popular total variation... more

descriptionView Paper arrow_downwardDownload

Two barriers on strong-stability-preserving time discretization methods

by Raymond Spiteri

2002

with equality holding if and only if A,(a, 8) = 0. But A,(a,Z) £0 for a method of order p so which is a contradiction. We therefore conclude that the CFL coefficient must be less than or equal to (s — p+ 1). We note that the order conditions that are relevant for linear, constant- coefficient problems are the so-called sub-quadrature or tall-tree [5] condi- tions given in Butcher notation by:

descriptionView Paper arrow_downwardDownload

Finite Difference WENO Schemes with Lax--Wendroff-Type Time Discretizations

by Chi-wang Shu

In this paper we develop a Lax-Wendroff time discretization procedure for high order finite difference weighted essentially nonoscillatory schemes to solve hyperbolic conservation laws. This is an alternative method for time... more

Fig. 1. The Burgers equation. u(z,0) = 0.54 sin(ra). t = 1.5/7. N = 80 points. Left WENO5-LW3 and WENOS-RK3; right: WENOS-LW4 and WENOS-RK4. TABLE 6 figure compares WENO5-LW4 with WENO5-RK4. N = 80 grid points are used. We can see that both schemes give equally good nonoscillatory shock transitions for this problem. Notice that the central difference approximations for the higher order time derivatives in the WENO5-LW schemes do not lead to spurious oscillations. THe SE EE OR Re eee bee cw een De ere ern eon eens bce ESia cla pees OT xeeereeddd seen Tk bec

Fic. 2. The Buckley—Leverett problem. t = 0.4. N = 160 points. Left: WENO5-LW3 and WENO5-RK8; right: WENO5-LW4 and WENOS-RK4.

Fic. 3. WENO5-LW4 comparing with WENO5-RK4. N = 200 points. Density. Left: the Lax problem at t = 0.16; right: the Sod problem at t = 0.1644.

Fic. 4. The shock density wave interaction problem. t = 1.8. by WENO5-LW4. Left: N = 200 jrid points; right: N = 400 grid points.

Fic. 5. The two-dimensional Burgers equation. u(x,0) = 0.5 + sin(m(a + y)/2). t = 1.5/7. 80 x 80 grid points, by WENO5S-LWA. Left: a cut of the solution at x = y, where the solid line is the exact solution and the squares are the computed solution; right: the surface of the solution.

Fic. 6. Double Mach reflection problem by WENO5-LW4 with 960 x 239 (top) and 1920 x 479 (bottom) grid points. 30 equally spaced contours from 1.5 to 22.7 for the density p.

Fic. 7. Zoomed-in figures. Double Mach reflection problem by WENO5-LWA4 with 960 x 239 (left) and 1920 x 479 (right) grid points. 30 equally spaced contours from 1.5 to 22.7 for the density p.

CPU time (in seconds) for the WENO5-LW and WENO5-RK schemes to compute the double Mach reflection problem in Example 3.9. the local characteristic decomposition is performed only for the first inner stage, the result becomes oscillatory, according to our numerical experiments. We note that the second and higher order time derivatives, when converted to spatial derivatives as before, involve expressions like f’(u), which is a matrix (the Jacobian), f”(w), which is a three-dimensional “matrix” (tensor), etc., which could become very complicated. A symbolic manipulator such as MAPLE would be helpful to avoid mistakes. The code is also quite long and messy compared with codes using Runge-Kutta time dis- cretizations. However, we will see in the next section that one can save CPU time by this approach for certain problems. TABLE 1

utt (3 a = = 0. u(x,0) = 0.5+sin(zxz). WENO5-LW4 using N equally spaced points. t = 0.5/7 v1 and Loo errors. equation, double Mach reflection test case in Example 3.9. The computer we used is an IBM PC with a Pentium-4 processor of 1.4G cpu and a 128M ram. We can see that the WENO5-LW schemes cost about half the CPU time compared with the WENO5-RK schemes for this example. Of course, this is closely related to the fact that the costly local characteristic decomposition and WENO approximation are used only once per time step in computing the first order time derivative, and the much less costly central difference approximations are used for all higher order time deriva- tives. However, in all the numerical experiments we have performed, nonoscillatory results can be obtained by this approach. We should remark, however, that there are also possibilities for reducing costs in Runge-Kutta methods, such as the possibility of freezing the nonlinear smoothness indicators for the inner stages, etc., which could potentially change the comparison in Table 1 but are not addressed in this paper.

urt(% a =0. u(x, 0) = 0.5+sin(ra). WENOS- -RKA using N equally spaced points. t = 0.5/7. Ly and Loo errors. TABLE 3

by the WENO5-LW4 scheme are shown in Table 4. We can see that the scheme achieves close to its designed order of accuracy. Example 3.3. We solve the following nonlinear svstem of Euler equations:

descriptionView Paper arrow_downwardDownload

A regional hydrodynamic model of upwelling in the Southern Benguela

by Pierre Fréon

2001

Hydrodynamic modelling is an active field in oceanography. The increase in available computing power as well as the incorporation of advanced features in numerical codes allow hydrodynamic models efficiently and robustly to simulate the... more

Fig. 1. The low-resolution spatial grid (left) and the vertical grid (right) used in the regional configuration of ROMS in the Southern Benguela. ROMS is a community model shared by a large user group around the world, with applications ranging from the study of entire ocean basins to coastal subregions. ROMS incorporates advanced features and high-order numerics, allowing efficient and robust resolution of mesoscale dynamics in the oceanic and coastal domains. The model solves the free surface, hydrostatic, primitive equations of the fluid dynamics over variable topogra- phy using stretched, terrain-following coordinates in the verti- cal and orthogonal, curvilinear coordinates in the horizontal. This allows enhancement of the spatial resolution in the regions of interest. Active open boundaries, connecting the regional model with the open ocean, are also incorporated.’ A pie-shaped grid that follows the southwestern corner of the African conti- nent from 40°S to 28°S and from 10°E to 24° E was developed. The topography is derived from the ETOPO2 dataset. Both a low- (18 km at the coast) and a high-resolution (9 km at the coast) grid are configured (Fig. 1). Along the vertical, stretched terrain- following coordinates and 20 levels provide enhanced resolu- tion at the surface while preserving an adequate resolution in the deeper layers (Fig. 1). The third-order, upstream-biased advection scheme implemented in ROMS allows the generation of steep gradients, enhancing the effective resolution of the solu- tion for a given grid size.° The model is forced with winds, heat and salinity fluxes extracted from the COADS ocean surface monthly climatology.’ At the three lateral boundaries facing the

Fig. 2. Twenty-year time series of volume-averaged temperature (°C) using the low-resolution configuration.

Fig. 3. Top: volume-averaged kinetic energy (cm’s-’) for the low-resolution experi- ment; bottom: surface-averaged kinetic energy (cm’s-’) for the low-resolution experiment.

Fig. 4. An example of the surface structure of the temperature and currents (arrows) in the Southern Benguela region simulated by the high-resolution configu- ration. This snapshot of the surface structure on 1 September of the second year of the simulation shows the main oceanographic features of the region: 1) the warm and intense Agulhas C urrent flowing westward in the southern part of the domain, and 2) the coastal upwelling off the west coast. Upwelling filaments off the west coast, a shear edge eddy on the eastern part of the Agulhas Bank as well as eddy shedding in the vicinity of the Agulhas Current, of the Cape Peninsula and of Cape Columbine are examples of the mesoscale activity simulated by the model. low-resolution configuration. Because the bathymetry is of primary importance in the control of the circulation dynamics over the continental shelf and slopes, it is expected that this configuration of the model would better resolve the topographi- cally controlled features observed over the continental shelf and adjacent waters. Whereas very small-scale processes are still not solved by the model, the 9-km resolution grid is fine enough to reproduce most of the mesoscale dynamics of the Benguela upwelling system. Using the same forcing conditions as in the low-resolution experiment, a 10-year simulation was conducted. The output data were stored every two simulated days in order to analyse the outputs with a high level of spatial and temporal accuracy as well as to provide the basis for providing inputs for biological models with high-resolution environmental data. An example of simulated sea-surface temperature (SST) over the entire domain is given in Fig. 4. A high level of mesoscale activity

Fig. 5. Comparison of the seasonal mean spatial structures of sea-surface temperature (°C) from the model and the PATHFINDER sSatelliteSST database.

descriptionView Paper arrow_downwardDownload

Finite-element solution of the coupled-channel Schrödinger equation using high-order accuracy approximations

by Dmitri Abrashkevich

1995, Computer Physics Communications

The finite element method (FEM) is applied to solve the bound state (Sturm-Liouville) problem for systems of ordinary linear second-order differential equations. The convergence, accuracy and the range of applicability of the high-order... more

descriptionView Paper arrow_downwardDownload

Fourth-order balanced source term treatment in central WENO schemes for shallow water equations

by Alessandro Valiani

2006, Journal of Computational Physics

The aim of this work is to develop a well-balanced Central Weighted Essentially Non Oscillatory (CWENO) method, fourth-order accurate in space and time, for shallow water system of balance laws with bed slope source term. Time accuracy is... more

descriptionView Paper arrow_downwardDownload

Integral deferred correction methods constructed with high order Runge-Kutta integrators

by Andrew Christlieb

2010, Mathematics of Computation

Spectral deferred correction (SDC) methods for solving ordinary differential equations (ODEs) were introduced by Dutt, Greengard and Rokhlin (5). It was shown in (5) that SDC methods can achieve arbitrary high order accuracy and possess... more

descriptionView Paper arrow_downwardDownload

High order conservative Lagrangian schemes with Lax–Wendroff type time discretization for the compressible Euler equations

by Juan Cheng and

2009, Journal of Computational Physics

In this paper, we explore the Lax-Wendroff (LW) type time discretization as an alternative procedure to the high order Runge-Kutta time discretization adopted for the high order essentially non-oscillatory (ENO) Lagrangian schemes... more

Fig. 1. The density of the Lax problem using the Dukowicz flux with 100 cells at t = 1. Left: second-order schemes; right: third-order schemes A reflective boundary condition is applied at both x = 0 and x = 1. We use an initially uniform mesh with 400 cells to compute this problem. Fig. 3 is the density at time t = 0.038. The left plot is the second-order result and the right one is

Fig. 2. The density of Noh problem using the Dukowicz flux with 200 cells at t = 0.6. Left: second-order schemes; right: third-order schemes

Fig. 3. The density of the blast waves problem using the HLLC flux with 400 cells at t = 0.038. Left: second-order schemes; right: third-order schemes. the third-order result. For this problem, it seems that the Lax-Wendroff time discretization does not produce any spurious overshoots near the contact discontinuity, which is better than the performance of the Runge-Kutta time discretization.

Fig. 4. Second-order results: internal energy of the Leblanc problem using 1000 cells at t = 6. Left: with the Dukowicz flux; right: with the HLLC flux.

Fig. 5. Third-order results: internal energy of the Leblanc problem using 1000 cells at t = 6. Left: with the Dukowicz flux; right: with the HLLC flux.

"ig. 6. Second-order results: density of the shock entropy wave interactions problem using 200 cells at t = 1.8. Left: with the Dukowicz flux; right: with th LLC flux.

Fig. 7. Third-order results: density of the shock entropy wave interactions problem using 200 cells at t = 1.8. Left: with the Dukowicz flux; right: with th HLLC flux. 3. Two-dimensional case 3.1. Construction of the scheme in the 2D case

The Taylor expansion tells us that in order to get a rth order accurate scheme, we should find a reconstruction to satisfy the following equations Since we only have cell average variables at t", we would need to use a Taylor expansion to get the value at the Gaussian quadrature points Similar to the 1D case, if we use a consistent and Lipschitz continuous numerical flux, combined with accurate numerical integration, we only need to find a reconstruction of U which is at least rth order accurate on the Gaussian quadrature points for a rth order accurate scheme

Again, this is automatically achieved by a rth order ENO reconstruction. Therefore, we can use an ENO reconstruction to get our second-order scheme on quadrilateral meshes and third-order scheme on curvilinear meshes.

Fig. 8. The initial mesh of the Saltzman problem.

Fig. 9. Second-order LW-ENO results for the Saltzman problem. Left of the top two rows: meshes; right of the top two rows: contours of density; bottom row: surfaces of the density. First row and left of the last row: results with the Dukowicz flux; Second row and right of the last row: results with the HLLC flux. The Dukowicz problem is a two-dimensional shock refraction problem on an uneven mesh designed by Dukowicz and Meltz [7].

Fig. 10. Third-order LW-ENO results for the Saltzman problem with the Dukowicz flux. Left of the top two rows: meshes; right of the top two rows contours of density; bottom row: surfaces of the density. First row and left of the last row: results on a quadrilateral mesh; Second row and right of the las row: results on a curvilinear mesh.

Fig. 11. LW-WENO results for the Sedov problem with the Dukowicz flux. Top row: meshes; middle row: contours of the density; bottom row: surfaces of the density. Left: second-order scheme; right: third-order scheme.

Fig. 12. The initial mesh for the Dukowicz problem.

Fig. 13. LW-WENO results for the Dukowicz problem with the HLLC flux. Top row: contours of the density; bottom row: surfaces of the density. Lef second-order results; right: third-order results.

Fig. 14. LW-WENDO results for the Double Mach problem with the Dukowicz flux. Top row: contours of the density; bottom row: surfaces of the densit: Left: second-order results; right: third-order results.

Second-order scheme. L; error. Runge-Kutta and Lax-Wendroff (boldface) time discretizations. Table 1

Third-order scheme. L; error. Runge-Kutta and Lax-Wendroff (boldface) time discretizations. Table 2 discretization of the same Lagrangian ENO scheme [3] are also included. We can see that the Lax-Wendroff time discretiz« tion produces smaller errors and costs less CPU time.

In order to compute p,,Mi, Ex, first Compute Uy, Uyy, Uxy, Vex, Yyy, Uxys Dawes Pyys Pays Pres yes Mxe, Mye, Nee, Nye, Exe, Eye, Uxe, Uye, Uxe; Yyt Pxt> Pyt

Second-order schemes. L, error. Runge-Kutta and Lax-Wendroff (boldface).

Third-order schemes on quadrilateral meshes with straight-line edges. L, error. Runge-Kutta and Lax-Wendroff (boldface).

Third-order schemes on curved quadrilateral meshes with quadratic curved edges. L; error. Runge-Kutta and Lax-Wendroff (boldface)

descriptionView Paper arrow_downwardDownload

Computation of free-surface flows and fluid–object interactions with the CIP method based on adaptive meshless soroban grids

by 弘宇江

2007, Computational Mechanics

are combined for computation of threedimensional fluid-object and fluid-structure interactions, while maintaining high-order accuracy. For the robust computation of free-surface and multi-fluid flows, we adopt the CCUP method [Phys Soc... more

descriptionView Paper arrow_downwardDownload

Numerical methods for reactive transport on rectangular and streamline-oriented grids

by Olaf Cirpka

1999, Advances in Water Resources

Coupling advection-dominated transport to reactive processes leads to additional requirements and limitations for numerical simulation beyond those for non-reactive transport. Particularly, both monotonicity avoiding the occurence of... more

The main differences for numerical simulation of transport on streamline-oriented grids compared to ar- bitrarily oriented ones may be explained by the trans- port equation of a reacting compound which is given for two dimensions in streamline coordinates by eqn (3) and in arbitrarily oriented coordinates by eqn (4): where R; is the retardation coefficient for compound i, 7; is a source-sink term due to reactive processes which

Fig. 1. Differences in the approximation of advective transport on rectangular versus streamline-oriented grids. As will be shown, adopting streamline-oriented grids does not only simplify the discretization of advective— dispersive transport but also avoids adding artificial transverse diffusion which may be necessary in order to achieve monotonicity on an arbitrarily oriented grid. This is of particular interest for the transport simulation of interacting compounds in which reaction rates are dependent on mixing of the compounds. On the other hand, the streamline-oriented grid-generation scheme? is restricted to steady-state flowfields, and the extension to

Fig. 2. Differences in the approximation of dispersive transport on rectangular versus streamline-oriented grids.

Fig. 4. Evaluation of fluxes related to the off-diagonal entry D*” at a vertex. (a) Diffusive and anti-diffusive flux at the vertex i+1/2,7+1/2. (b) Differentiation stencil for cell i,7 and choice of the points of evaluation: black dots: fluxes related to D»”, circles: fluxes related to D™, gray dots: fluxes related to D”®.

Fig. 3. Influence of cell (i, /) on cell (i+ 1,7 + 1) and vice versa by cross-diffusion terms. Black dots: interfaces for cross-dif- fusive fluxes related to cell (i, 7), grey dots: interfaces for cross- diffusive fluxes related to cell (i+ 1,7 + 1).

with clockwise numbering of the cells around the vertex and the bilinear shape function N. In the FEM, the approximation of the related temporal concentration change is of course different to the cell-centered FVM. Only for the case of regular spacing in both directions the related volumes are identical. For this particular case eqn (9) can also be retrieved by the FVM.' Nevertheless, in the present study eqn (9) is also applied to grids with irregular spacing. Equation (9) may diagonal diffusive fluxes across the vertex in which the product of the diffusion coefficient times the interfacial area divided by t ne d between the lower —1/2@D” for the flux between the upper left and lower right cell, respec eqn (5) it is clear ative. It vanishes grid and the dis ive be inter istance eft and on. preted as summation of two equals 1/2¢D” for the flux the upper right cell and y. From the definition of @¢D” in that its va ue may be positive or neg- y if the principal directions of the persion tensor coincide. For all other cases one relation between diagonally positioned cells will be of a diffusive nature but with a negative diffusion coefficient. This is refered to as anti-diffusive relation. It is important to notice that anti-diffusion does not occur in nature, leads to unphysical sharpening of concentra- tion distributions and may cause oscillations in the vi- cinity of discontinuities.

Fig. 5. Principle of the Corner Transport Upwind method (CTU) '°. Fluxes originated from cell i, j. Numerical methods for reactive transport on rectangular and streamline-oriented grids with the correction factor C;,’, ; of cell (i+ 1,/) related to the outward-directed flux g;,, ,,,;. as defined by:

mentioned that this approximation is slightly different to the one obtained by the Lax—Wendroff scheme’’. However, both authors claim second-order accuracy if only the directional second-order derivatives are added to the CTU solution. This results in the following high- order fluxes: LeVeque**** suggested that the transverse propaga- tion of the low-order flux may be applied as well to the second-order correction terms. He did this by modifying cross-diffusion terms. These terms may again be refor- mulated so that diagonal fluxes occur directly. This yields to the following additional modifications related to the correction term B*,,,, ;:

Fig. 6. Approximation of gradients oriented into the direction of flow and orthogonal to it. (a) Perfectly streamline-oriented ele- ment; (b) Approximation by a bilinear quadrilateral; (c) Quadrilateral in local coordinates. C: Center of gravity with respect to the local coordinates €, 7.

Fig. 7. Apparent widths and lengths of element k. Direction of flow: cells k —i— m. from element *& into the direction of flow, whereas the elements A and m are directed into the direction of pseodopotential lines. For illustration see Fig. 7.

Fig. 8. Definition of slopes for element k. Open circles: aver- aged concentrations in elements p (upstream), k and q (downstream).

Fig. 10. Five-spot configuration in a heterogeneous aquifer. Streamline-oriented grid. Numerical methods for reactive transport on rectangular and streamline-oriented grids

Fig. 11. Five-spot configuration in a heterogeneous aquifer. Concentration distribution of a conservative tracer for the streamline-oriented grid three days after the start of injection.

Fig. 12. Five-spot configuration in a heterogeneous aquifer. Concentration distribution of a conservative tracer for the rectangular grid three days after the start of injection.

Fig. 14. Reactive transport in a perfect well couple. Stream- line-oriented grid adopted for the test case. Bold line: dividing streamline.

Fig. 13. Five-spot configuration in a heterogeneous aquifer. Breakthrough curves obtained at the outflow boundary for the rectangular and the streamline-oriented grid.

Fig. 15. Reactive transport in a perfect well couple. Concen- tration distribution 50 d after start of the injection adopting the streamline-oriented grid. % = 1 m, « =0.1 m.

Fig. 16. Reactive transport in a perfect well couple. Total mass of compound C 50 d after start of the injection as function of the transverse dispersivity. Comparison of grids adopted.

Fig. 17. Reactive transport in a perfect well couple. Total mass of compound C 50d after start of the injection calculated on the rectangular grid under assumption of strictly advective transport as function of grid refinement.

descriptionView Paper arrow_downwardDownload

HSTERM — A program to calculate potential curves and radial matrix elements for two-electron systems within the hyperspherical adiabatic approach

by Dmitri Abrashkevich

1995, Computer Physics Communications

A FORTRAN 77 program is pw.sented which calculates potential curves and matrix elements of radial coupling for twoelectron systems using the hyperspherical coordinate method. The adiabatic and diabatic-by-scctor close-coupling approaches... more

descriptionView Paper arrow_downwardDownload

KANTBP: A program for computing energy levels, reaction matrix and radial wave functions in the coupled-channel hyperspherical adiabatic approach

by Alexander Abrashkevich

2007, Computer Physics Communications

A FORTRAN 77 program is presented which calculates energy values, reaction matrix and corresponding radial wave functions in a coupledchannel approximation of the hyperspherical adiabatic approach. In this approach, a multi-dimensional... more

where O. Chuluunbaatar et al. / Computer Physics Communications 177 (2007) 649-675 where I7(®, E, Pmax) is the symmetric functional

Fig. 1. Flow diagram of the KANTBP program. O. Chuluunbaatar et al. / Computer Physics Communications 177 (2007) 649-675

In case of the repulsive interaction, we have used asymptotic solutions (84) with and indexes i, j start from 1. In order to compare with results [31], 60 = z/4, and 2E = 0.01 (k = 0.1) have been chosen. The matrix-solutions (84), (86) are included in the subroutine ASYMSC by default for c > 0. The following values of numerical parameters and characters have been used in the test run via the supplied input file 3DDSCP.INP

5. Description of subprogram units O. Chuluunbaatar et al. / Computer Physics Communications 177 (2007) 649-675

SUBDIVISION OF RHO-REGION ON THE FINITE-ELEMENT GROUPS: KKEKKKKKKKKKKKKKKKKKKEKKK KKK KKK KK KKKKKKKKKKKKKKK KKK KK KKK

O. Chuluunbaatar et al. / Computer Physics Communications 177 (2007) 649-675 SUBDIVISION OF RHO-REGION ON THE FINITE-ELEMENT GROUPS: KKEKKKKKKKKKK KKK KKK KKK KKK KKK KKK KEK KKK KKKEKKEKKKK KKK KK KKK KK

descriptionView Paper arrow_downwardDownload

High-Order Well-Balanced Finite Difference WENO Schemes for a Class of Hyperbolic Systems with Source Terms

by Chi-wang Shu

2006, Journal of Scientific Computing

In this paper, we generalize the high order well-balanced finite difference weighted essentially non-oscillatory (WENO) scheme, designed earlier by us in Xing and Shu (2005, J. Comput. phys. 208, 206-227) for the shallow water equations,... more

Fig. 1. The impedance Z(x) and the sound speed c(x) for the smooth media. 486 It is a test case where the impedance Z(x) and hence the eigenvectors are both spatially varying. We perform the computation with 200 uniform cells, with the ending time t=0.4s. An “exact” reference solution is com- puted with the same scheme over a 2000 grid point uniform cells. The sim- ulation results are shown in Fig. 2. The numerical resolution shows very good agreement with the “exact” reference solution.

Fig. 2. The numerical (symbols) and the “exact” reference (solid line) stress o(x) at time t=0.4s. High-Order Well-Balanced Finite Difference WENO Schemes

Fig. 3. The shape of a convergent-divergent nozzle. High-Order Well-Balanced Finite Difference WENO Schemes

Fig. 4. Steady state pressure for the nozzle flow.

fable I. ZL! Errors and Numerical Orders of Accuracy for the Example in Sec. 3.2 485 with the density o(x) and bulk modulus of compressibility K(x) given by:

Table II. ZL! Errors for Different Pre- cisions for the Stationary Solution in Sec. 3.3 10 Sj =Sj(Qj)=e es” “6%, $9 =S72(a2) =ONU, xX) =—er* vv”, ANG 12(X) =X. We now show two numerical examples to demonstrate the fifth order well balanced finite difference WENO scheme for (3.3). The first example is to test the well balancedness property of the scheme. We take the initial conditions as

Table II. ZL! Errors and Numerical Orders of Accuracy for the Example ir Sec. 3.3. with a periodic boundary condition. Since the exact solution is not known explicitly for this problem, we use the same fifth-order WENO scheme with N =5120 points to compute a reference solution and treat it as the exact solution when computing the numerical errors. Final time tf =0.5s is used to avoid the development of shocks. Table III contains the L! errors and numerical orders of accuracy. We can clearly see that fifth-order accu- racy is achieved for this example.

Table IV. L!Errors and Numerical Orders of Accuracy for the Example in Sec. 3.4

descriptionView Paper arrow_downwardDownload

Hermite WENO schemes for Hamilton–Jacobi equations on unstructured meshes

by Chi-wang Shu

2013, Journal of Computational Physics

In this paper, a class of weighted essentially non-oscillatory (WENO) schemes based on Hermite polynomials, termed HWENO (Hermite WENO) schemes, for solving Hamilton-Jacobi equations is presented. The idea of the reconstruction in the... more

descriptionView Paper arrow_downwardDownload

Filtering in Legendre spectral methods

by Jan S Hesthaven

2008, Mathematics of Computation

We discuss the impact of modal filtering in Legendre spectral methods, both on accuracy and stability. For the former, we derive sufficient conditions on the filter to recover high order accuracy away from points of discontinuity.... more

descriptionView Paper arrow_downwardDownload

On the Nonlinear Hydrodynamic Stability of Thin Keplerian Disks

by Mario Livio

1999, The Astrophysical Journal

descriptionView Paper arrow_downwardDownload

Robust HLLC Riemann solver with weighted average flux scheme for strong shock

by Hyoung Jin Lee

2009, Journal of Computational Physics

Many researchers have reported failures of the approximate Riemann solvers in the presence of strong shock. This is believed to be due to perturbation transfer in the transverse direction of shock waves. We propose a simple and clear... more

The HLL intercell flux is written as [3,4] where S, is the smallest wave speed and Sp is the largest wave speed. The subscripts L and R imply the left and right values at the cell interface, respectively. The subscripts “Land R imply the left and right values at the contact discontinuity between two nonlinear waves, respectively. The HLL Riemann solver assumes a single constant state between two nonlinear waves (shock or rarefaction) [3,4]. This assumption averages the spatial variations across the contact discontinuity resulting in a smeared solution of contact and shear waves [4]. The HLLC scheme is a modification of the HLL scheme wherein the missing contacts and shear waves are restored [4,6].

where the subscripts n and t represent the normal and tangential velocity components, and S, is the middle wave speed. To compute wave speeds S,, Sr, and S., the pressure-velocity based wave estimations presented by Toro [4] are used to estimate the shock and the rarefaction waves accurately. More detailed expressions are presented in [4]. Einfeldt proposed wave speeds motivated by the Roe-averaged values i, and 4 [5]. Using the modified wave speeds, the HLL scheme is often called the HLLE scheme [4,5,13].

Fig. 1. Odd-even grid perturbation problem: density contours (second-order WAF with MINBEE (Minmod) limiter).

Fig. 2. Steady normal shock test problem: density contours (second-order WAF with MINBEE (Minmod) limiter).

Fig. 3. Supersonic circular blunt body: density and switching function distributions (second-order WAF with MINBEE (Minmod) limiter)

Fig. 4. Kinked Mach stem problem: density and switching function distributions (second-order WAF with MINBEE (Minmod) limite

Fig. 5. Supersonic corner flow: density and switching function distributions (second-order WAF with MINBEE (Minmod) limiter).

descriptionView Paper arrow_downwardDownload

Modelling of rupture propagation using high-order mimetic finite differences

by Jose Castillo

2008, Geophysical Journal International

We present a split-node finite difference (FD) method for modelling shear ruptures that is consistently fourth-order accurate in its spatial discretization, both in the interior of the model and at the fault. The method, called mimetic... more

descriptionView Paper arrow_downwardDownload

High order conservative differencing for viscous terms and the application to vortex-induced vibration flows

by Ge-cheng Zha

2009, Journal of Computational Physics

A set of conservative 4th-order central differencing schemes for the viscous terms of Navier-Stokes equations are proved in this paper. These schemes are used with the 5th order WENO schemes for inviscid flux. The stencil width of the... more

descriptionView Paper arrow_downwardDownload

Semi discrete discontinuous Galerkin methods and stage-exceeding-order, strong-stability-preserving Runge–Kutta time discretizations

by Clint Dawson

2007, Journal of Computational Physics

This paper investigates the use of a special class of strong-stability-preserving (SSP) Runge-Kutta time discretization methods in conjunction with discontinuous Galerkin (DG) finite element spatial discretizatons. The class of SSP... more

descriptionView Paper arrow_downwardDownload

High order residual distribution conservative finite difference WENO schemes for convection–diffusion steady state problems on non-smooth meshes

by Chi-wang Shu

2007, Journal of Computational Physics

In this paper, we propose a high order residual distribution conservative finite difference scheme for solving steady state hyperbolic conservation laws on non-smooth Cartesian or other structured curvilinear meshes. WENO (weighted... more

descriptionView Paper arrow_downwardDownload

A stable and high-order accurate conjugate heat transfer problem

by Jens Lindström

2010, Journal of Computational Physics

This paper analyzes well-posedness and stability of a conjugate heat transfer problem in one space dimension. We study a model problem for heat transfer between a fluid and a solid. The energy method is used to derive boundary and... more

descriptionView Paper arrow_downwardDownload

A staggered grid, high-order accurate method for the incompressible Navier–Stokes equations

by Nikolaos A. Kampanis

2006, Journal of Computational Physics

A high-order accurate, finite-difference method for the numerical solution of the incompressible Navier–Stokes equations is presented. Fourth-order accurate discretizations of the convective and viscous fluxes are obtained on staggered... more

Fig. 1. Schematic of staggered grid. and using functional values at integer points x; as

Fig. 2. Comparison of the u-velocity computed by the fourth- and second-order accurate methods with a reference solution by Ghia et a [9]; driven cavity flow (Re = 100).

Fig. 3. Comparison of the v-velocity computed by the fourth- and second-order accurate methods with a reference solution by Ghia et al. [9]; driven cavity flow (Re = 100).

Fig. 4. Velocity field and velocity magnitude contours computed by the fourth-order method on the 221 x 221 point grid; driven cavity flow (Re = 100).

Fig. 5. Comparison of the L” norm of the residuals R of the numerical solution for the fourth- (on the 221 x 221 points grid) and the second-order (on the 221 x 221, 301 x 301 point grids) accurate methods; driven cavity flow (Re = 100).

Fig. 6. Grid convergence shown by the L* norms of the error obtained by the fourth- and second-order accurate methods; driven cavity flow (Re = 100). The decay of an ideal vortex is an unsteady flow problem with an exact solution [26]. This problem is of interest to numerical simulations of trailing vortices, (LES), and (DES) simulations. The initial velocity dis- tribution of the Oseen vortex is given by

Fig. 7. Comparison of u-velocity computed by the fourth- and second-order accurate methods on the 51 x 51 point grid with the exact solution at T = 4; Oseen vortex decay. MT EE TO AS ee, Se A Bel ee ee ene: In Fig. 12, the kinetic energy defined as E(t) = f,, w’(r, 0,t)rdrd0 is plotted for Re = 1, for the exact solu- tion and the fourth-order accurate numerical nate, raped on the 101 x 101 point ori using the explicit finite-differences of Section 3.2.1 to discretize in space the momentum equations. The kinetic energy is plotted in the same figure for the inviscid limit Re = 10'°, for the exact solution and the fourth-order accurate numer- ical solution on the 1001 x 1001 point grid. In both cases, the energy conservation of the proposed fourth- order accurate numerical method is evident. Fourth-order accurate numerical solutions computed on the same grid and the compact schemes of Section 3.2.2 for the space discretization, exhibited the same behavior. For the high Reynolds number solutions the numerical mesh is shifted to avoid having a grid point at the origin

Fig. 8. Comparison of u-velocity computed by the fourth- and second-order accurate methods on the 101 x 101 point grid with the exact solution at T = 4; Oseen vortex decay.

Fig. 9. Comparison of the L” norm of the spatial error obtained by the fourth- and second-order accurate methods at T = 1; Oseen vortex decay.

Fig. 10. Comparison of the L” norm of the temporal error obtained by the fourth-order accurate in space method at T = 1; Oseen vorter decay.

Fig. 11. Comparison of u-velocity computed by the fourth- and second-order accurate methods on the 1001 x 1001 point grid with exact solution for Re = 10!°, at T= 4; Oseen vortex decay.

Fig. 12. Comparison of the kinetic energy E as a function of time obtained by the w-velocity computed by the fourth-order accurate method for Re = 1 (on the 101 x 101 point grid) and Re = 10!° (on the 1001 x 1001 point grid), and the corresponding exact solution; Oseen vortex decay.

Fig. 13. Comparison of u-velocity computed by the fourth- and second-order accurate methods on the 6 x 51 point grid with the exact solution at T= 1/4; Stokes oscillating plate. Owl Yue. Comparisons of the computed solution, from both numerical methods, (at T = 1/4 after a time periodic solution was achieved) with the exact solution (48), are shown in Fig. 13 for the 6 x 101 point grid. The agree- ment with the exact solution is very good and the better resolution of the fourth-order method is evident. The variation of the L* norm of the error obtained by both methods with grid refinement are plotted in Fig. 14. It is evident that grid convergence is achieved and the proposed fourth-order accurate method retains fourth order of accuracy. It was found that the variant of the fourth-order method with compact schemes for the

Fig. 14. Grid convergence shown by the L? norm of the errors obtained by the fourth- and second-order accurate methods at T= 1/4 Stokes oscillating plate.

Fig. 15. Comparison of v-velocity component computed by the fourth- and second-order accurate methods on the 41 x 101 point grid with the exact solution at T= 1/4; pulsatile flow with « = 6.3245 and w= 4.

Fig. 16. Comparison of the computed and exact axial velocity variation over half period at r= R/2 and r = 3R/4; pulsatile flow with a = 6.3245 and wm =4.

descriptionView Paper arrow_downwardDownload

High-Order Finite-Volume Methods on Locally-Structured Grids

by Phillip Colella

2009

We describe recent developments in high-order (greater than second-order) accurate finite-volume discretizations of partial differential equations (PDE) in divergence form. We will address a number of algorithmic issues that arise in... more

descriptionView Paper arrow_downwardDownload

Higher-order mimetic methods for unstructured meshes

by Blair Perot

2006, Journal of Computational Physics

A higher-order mimetic method for the solution of partial differential equations on unstructured meshes is developed and demonstrated on the problem of conductive heat transfer. Mimetic discretization methods create discrete versions of... more

descriptionView Paper arrow_downwardDownload

Application of the discontinuous spectral Galerkin method to groundwater flow

by Sergio Fagherazzi

2004, Advances in Water Resources

The discontinuous spectral Galerkin method uses a finite-element discretization of the groundwater flow domain with basis functions of arbitrary order in each element. The independent choice of the basis functions in each element permits... more

descriptionView Paper arrow_downwardDownload

Construction of low dissipative high-order well-balanced filter schemes for non-equilibrium flows

by Wei Wang and

2011, Journal of Computational Physics

Recent progress in the development of a class of low-dissipative high-order filter schemes for multiscale Navier-Stokes and magnetohydrodynamics (MHD) systems by , and shows good performance in multiscale shock/turbulence simulations. The... more

descriptionView Paper arrow_downwardDownload

ODPEVP: A program for computing eigenvalues and eigenfunctions and their first derivatives with respect to the parameter of the parametric self-adjoined Sturm–Liouville problem

by Alexander Abrashkevich

2009, Computer Physics Communications

A FORTRAN 77 program is presented for calculating with the given accuracy eigenvalues, eigenfunctions and their first derivatives with respect to the parameter of the parametric self-adjoined Sturm-Liouville problem with the parametric... more

This article appeared in a journal published by Elsevier. The attachec copy is furnished to the author for internal non-commercial research and education use, including for instruction at the authors institution and sharing with colleagues. Provided for non-commercial research and education use. Not for reproduction, distribution or commercial use.

The functions {NP (Z)}Los L=np, form a basis in the space of polynomials of the pth order. Now, each function w(z; o) is approximated by a finite sum of local functions NP (Z) are determined. By means of the Lagrange elements o? ,(Z), we define a set of local functions N;(z) as follows After substituting expansion (20) into the variational functional (15) and minimizing it [10,11] we obtain the generalized eigenvalue problem

» NMESH1 is the dimension of the DOUBLE PRECISION array RMESH containing the information about the subdivision of the interval [Zmin, Zmax] on subintervals and the number of elements on each one of them. NMESH1 is always odd and > 3. » NROOT1 is the number of eigenvalues and eigenvectors required, and also the dimension of the DOUBLE PRECISION arrays HH, QQ, EIGV.

Thus, problem (6)-(8) with additional conditions (11), (12) has now a unique solution. In the most of applications the following formulas The calculated eigenvalues and potential matrix elements can be used for solving the bound state and multi-channel scattering problems for a system of the coupled second-order ordinary differential equations with the help of the KANTBP programs [8,9]. Taking a derivative of the boundary problem (1)-(3) with respect to parameter p, we get that dywj(z; 0)/do can be obtained as a solution of the following boundary problem

descriptionView Paper arrow_downwardDownload

Klein–Gordon equation with advection on unbounded domains using spectral elements and high-order non-reflecting boundary conditions

by Beny Neta

2010, Applied Mathematics and Computation

A reduced shallow water model under constant, non-zero advection in the infinite channel is considered. High-order (Givoli-Neta) non-reflecting boundary conditions are introduced in various configurations to create a finite computational... more

descriptionView Paper arrow_downwardDownload

High-order accurate modeling of electromagnetic wave propagation across media – Grid conforming bodies

by Eli Turkel

2006, Journal of Computational Physics

The Maxwell equations contain a dielectric permittivity e that describes the particular media. For homogeneous materials at low temperatures this coefficient is constant within a material. However, it jumps at the interface between... more

Fig. 1. Approximation of the spatial derivatives at the nodes/half-nodes. The total spatial portion of the scheme, in matrix form, is given in Fig. 1

At each interface L, and L» we supplement (4) by the jump conditions: E. Kashdan, E. Turkel | Journal of Computational Physics xxx (2006) xxx-xxx

We then define in the interval x,_12 < x < Xx41/2 the following parameters:

Fig. 3. Approximation of relative permittivity by continuous function for model problem. The initial conditions are given by (13) with t = 0. This solves Maxwell’s equations (4) together with the jump conditions (5). We solve numerically the dimensionless equations (c = 1) and E and H are normalized to have the same magnitude. In the following plots we show the L, error of the approximation of FE as function of time. The error in H behaves similarly. We locate the interface at a node so only the electric field component is located on interface. The “regularization” error QO; between the continuous solution based on variable ¢ and piecewise continuous ¢é (denoted as éjm) Is given by

Fig. 4. Arithmetic (a) and harmonic (h) averaging of és on interface, Yee scheme.

Fig. 5. Arithmetic (a) and harmonic (h) averaging of ¢ on interface, fourth order scheme.

Fig. 6. Regularization of ¢ by Hermite cubic spline w/wo monotonicity restraint (m).

Fig. 7. Regularization of dielectric permittivity on interface by high-order polynomial with fixed (f) and variable (v) length of d-interval.

Fig. 8. Global regularization of ¢ and 1/e using Padé interpolation.

Fig. 9. Global regularization of ¢ and 1/e using Padé interpolation for Yee scheme.

Any coordinate system that is not aligned with the bodies has the disadvantage that the body cannot be represented exactly. A general body in a Cartesian coordinate system gives rise to stair-casing and its resultant errors (see, for instance, [4,14]). In this paper we only consider bodies aligned with the coordinate system and so staircasing does not occur. Consider a PEC sphere surrounded by two different homogeneous media separated with an interface in the radial direction. Each medium has its own dielectric permittivity ¢ (the outermost medium is free space). This problem shown in Fig. 10 has no explicit solution. The Maxwell eqrnatinne in cnherical eanrdinatee (r A m) are oiven hv: In addition to the time dependent equations we have Gauss’ law. In the absence of sources both the divergence of £ and H are zero:

Fig. 10. Propagation of electromagnetic pulse in inhomogeneous media.

To prevent reflections from the outer artificial boundary in the radial direction we implement a perfectly natched layer (PML). We consider the generalization of the uniaxial PML, derived by Gedney in [8], tc pherical coordinates. We convert the Maxwell equations to Fourier space as Teixeira and Chew [23] have uggested:

and ae for H. We replace E by E* and add two variables P, and Q,. This is converted to the time-domain using iw — 2. This yields:

Fig. 12. Converted to Cartesian coordinates. E,. component of the total electric field over entire physical domain. GU dlisOU LU LOU CULUIEU OUsSeiaulOls Li fry Tb Ue VIRIILY OL UI UIeLlAcr, In Fig. 13 we show a slice of the Ey component in the 6 direction (taken at the source). This shows that some of the waves are reflected by the interface. The physical explanation of this phenomena that refraction causes the waves to change direction and to move from the selected slice, needs do be combined with the numerical problem of under-resolution since the Gaussian pulse includes a full spectrum of wavelengths. When the waves travel through a medium with a smaller speed of light then some of the waves can be only be represented on a finer erid.

Fig. 13. The slice of E, component in 0-direction, taken at source point. 7. Concluding remarks

Relative permittivities for different materials Table 1

Regularization (Q)) error for various é Table 2

Total error using global regularization for various ¢ Table 4 Total error using local regularization for various &

descriptionView Paper arrow_downwardDownload

The discontinuous Galerkin method for fractal conservation laws

by Espen Jakobsen and

2011, IMA Journal of Numerical Analysis

We propose, analyze, and demonstrate a discontinuous Galerkin method for fractal conservation laws. Various stability estimates are established along with error estimates for regular solutions of linear equations. Moreover, in the... more

Note that, since both e,e € H*/?(R), each term in the above expression is wel defined (cf. the discussion in the proof of Theorem |3.3). One can argue as in [14] Theorem 2.1] to bound the local terms by cep Aa® . Hence,

Proof. The plan is to estimate —A, s/t, u], and, then, use Lemma [4.3] to conclude. Proof for the implicit-explicit method. Let us introduce the notation a A 6b = min{a, b}, aVb = max{a, b}, n? = |UP — ul and q? = f(UP Vu) — f(UP Au), where u=u(y,s). Note that —A..s[t, u] can be rewritten as literature up to now (cf. Droniou [19] for an alternative convergence proof, without convergence rate, for the implicit-explicit case).

Next, the right-hand side of the above inequality needs to be estimated. To this end, let us point out that, as proved in 23) Example 3.14],

Proceeding as done in the proof for the implicit-explicit method, we obtain the cell entropy inequality Proof. Let us give the proof using the implicit-explicit method say, the fully explicit method 4.2 would also do. 4.1 . Needless

FIGURE 1. Initial data (piecewise linear) and solutions of the pure fractional equation (A = 0.5) with k = 0 and Az = 1/160.

FIGURE 2. Initial data (piecewise linear) and solutions of the frac- tional transport equation (A = 0.5) with k = 0 and Az = 1/160.

FIGURE 3. Initial data and solutions of the fractional Burgers equation (A = 0.5) using & = 0; T = 0.5 and Ax = 1/160.

FIGURE 4. Initial data and solutions of the fractional Burgers’ equation for different values of A using k = 0; T = 0.5, Ax = 1/200, and uo(#) = —arctan(15x)/90.

FIGURE 5. Initial data (piecewise linear) and solutions of the frac- tional Burgers’ equation (A = 0.5) at different times T using k = 0; Ax = 1/200.

We now prove the remaining J; (i = 2,3,4) to be of order Ax?*+?->. First note that since v is smooth on each interval J;, the fundamental theorem of calculus First of all, let us bound from above the H*/?-norm of v as

descriptionView Paper arrow_downwardDownload

A spectral element solution of the Klein-Gordon equation with high-order treatment of time and non-reflecting boundary

by Beny Neta

2010, Wave Motion

A spectral element (SE) implementation of the Givoli-Neta non-reflecting boundary condition (NRBC) is considered for the solution of the Klein-Gordon equation. The infinite domain is truncated via an artificial boundary B, and a... more

descriptionView Paper arrow_downwardDownload

A high order accurate conservative remapping method on staggered meshes

by Juan Cheng and

2008, Applied Numerical Mathematics

An accurate remapping algorithm is an essential component of the Arbitrary Lagrangian-Eulerian (ALE) methods. Most ALE codes applied to high speed flow problems use a staggered mesh, i.e., all the solution variables except the velocities... more

descriptionView Paper arrow_downwardDownload

Log In

High Order Accuracy

Related Topics