Search

Fixed-point fast sweeping WENO methods are a class of efficient high-order numerical methods to solve steady-state solutions of hyperbolic partial differential equations (PDEs). The Gauss-Seidel iterations and alternating sweeping strategy are used to cover characteristics of hyperbolic PDEs in each sweeping order to achieve fast convergence rate to steady-state solutions. A nice property of fixed-point fast sweeping WENO methods which distinguishes them from other fast sweeping methods is that they are explicit and do not require inverse operation of nonlinear local systems. Hence, they are easy to be applied to a general hyperbolic system. To deal with the difficulties associated with numerical boundary treatment when high-order finite difference methods on a Cartesian mesh are used to solve hyperbolic PDEs on complex domains, inverse Lax-Wendroff (ILW) procedures were developed as a very effective approach in the literature. In this paper, we combine a fifth-order fixed-point fast sweeping WENO method with an ILW procedure to solve steady-state solution of hyperbolic conservation laws on complex computing regions. Numerical experiments are performed to test the method in solving various problems including the cases with the physical boundary not aligned with the grids. Numerical results show high-order accuracy and good performance of the method. Furthermore, the method is compared with the popular third-order total variation diminishing Runge-Kutta (TVD-RK3) time-marching method for steady-state computations. Numerical examples show that for most of examples, the fixed-point fast sweeping method saves more than half CPU time costs than TVD-RK3 to converge to steady-state solutions.

The discontinuous Galerkin (DG) method is widely used in numerical solution of partial differential equations, especially for hyperbolic equations. However, for problems containing strong shocks, the DG method often needs to be supplemented by a limiter to control spurious oscillations and to ensure nonlinear stability. The total variation bounded (TVB) limiter is a popular choice and can maintain the original high order accuracy of the DG scheme in smooth regions and keep a sharp and non-oscillatory discontinuity transition, when a certain TVB constant M is chosen adequately. For scalar conservation laws, suitable choice of this constant M can be based on solid mathematical analysis. However, for nonlinear hyperbolic systems, there is no rigorous mathematical guiding principle for the determination of this constant, and numerical experiments often use ad hoc choices based on experience and through trial and error. In this paper, we develop a TVB constant artificial neural network (ANN) based estimator by constructing a multi-layer perceptron (MLP) model. We generate the training data set by constructing piecewise smooth functions containing local maxima, local minima, and discontinuities. By using the supervised learning strategy, the MLP model is trained offline. The proposed method gives the TVB constant M with robust performance to capture sharp and non-oscillatory shock transitions while maintaining the original high order accuracy in smooth regions. Numerical results using this new estimator in the TVB limiter for DG methods in one and two dimensions are given, and its performance is compared with the classical ad hoc choices of this TVB constant.

Robust numerical methods for CFD applications, such as WENO schemes, quickly evolved in the past few decades. Together with the Inverse Lax–Wendroff (ILW) procedure, WENO ideas were also applied in the boundary treatment. Those methods are known for their high-resolution property, i.e., good representation of nonlinear phenomena, which is an important property in solving challenging engineering problems. In light of that, the objective of this work is to present a review of well-established high-resolution numerical methods to solve the Euler equations and adapt the Navier–Stokes viscous terms discretization and boundary treatment. To test the modifications, we employed the positivity-preserving Lax–Friedrichs splitting, multi-resolution WENO scheme, third-order strong stability preserving Runge–Kutta time discretization, and ILW boundary treatment. The first problems were simple flows with analytical solutions for accuracy tests. We also tested the accuracy with nontrivial phenomena in the vortex flow. Oblique shock and complicated flow structures were captured in the Rayleigh–Taylor instability and flow past a cylinder. We showed the discretization and boundary treatment can handle non-constant viscosity, are high-order, high-resolution, and behave similarly to the well-established numerical methods. Furthermore, the methods discussed here can preserve symmetry and no approximations regarding the boundary layer were made. Therefore, the discretization and boundary treatment can be considered when solving direct numerical simulations.

In this article, we introduce a new WENO algorithm that aims to calculate an approximation to derivative values of a function in a non-regular grid. We adapt the ideas presented in [Amat et al., SIAM J. Numer. Anal. (2020)] to design the nonlinear weights in a manner such that the order of accuracy is maximum in the intervals close to the discontinuities. Some proofs, remarks on the choice of the stencils and explicit formulas for the weights and smoothness indicators are given. We also present some numerical experiments to confirm the theoretical results.

We propose and analyze a class of robust, uniformly high-order accurate discontinuous Galerkin (DG) schemes for multidimensional relativistic magnetohydrodynamics (RMHD) on general meshes. A distinct feature of the schemes is their physical-constraint-preserving (PCP) property, i.e., they are proven to preserve the subluminal constraint on the fluid velocity and the positivity of density, pressure, and internal energy. This is the first time that provably PCP high-order schemes are achieved for multidimensional RMHD. Developing PCP high-order schemes for RMHD is highly desirable but remains a challenging task, especially in the multidimensional cases, due to the inherent strong nonlinearity in the constraints and the effect of the magnetic divergence-free condition. Inspired by some crucial observations at the PDE level, we construct the provably PCP schemes by using the locally divergence-free DG schemes of the recently proposed symmetrizable RMHD equations as the base schemes, a limiting technique to enforce the PCP property of the DG solutions, and the strong-stability-preserving methods for time discretization. We rigorously prove the PCP property by using a novel “quasi-linearization” approach to handle the highly nonlinear physical constraints, technical splitting to offset the influence of divergence error, and sophisticated estimates to analyze the beneficial effect of the additional source term in the symmetrizable RMHD system. Several two-dimensional numerical examples are provided to further confirm the PCP property and to demonstrate the accuracy, effectiveness and robustness of the proposed PCP schemes.

When solving CFD problems, the solver, or the numerical code, plays an important role. Depending on the phenomena and problem domain, designing such numerical codes can be hard work. One strategy is to start with simple problems and construct the code as building blocks. The purpose of this work is to provide a detailed review of the theory to compute analytical and exact solutions, and recent numerical methods to construct a code to solve compressible and inviscid fluid flows with high-resolution, arbitrary domains, non-linear phenomena, and on rectangular meshes. We also propose a modification to the inverse Lax–Wendroff procedure solid wall treatment and two-dimensional WENO-type extrapolation stencil selection and weights to handle more generic situations. To test our modifications, we use the finite difference method, Lax–Friedrichs splitting, WENO-Z+ scheme, and third-order strong stability preserving Runge-Kutta time discretization. Our first problem is a simple one-dimensional transient problem with periodic boundary conditions, which is useful for constructing the core solver. Then, we move to the one-dimensional Rayleigh flow, which can handle flows with heat exchange and requires more detailed boundary treatment. The next problem is the quasi-one-dimensional nozzle flow with and without shock, where the boundary treatment needs a few adjustments. The first two-dimensional problem is the Ringleb flow, and despite being smooth, it has a curved wall as the left boundary. Finally, the last problem is a two-dimensional conical flow, which presents an oblique shock and an inclined straight line wall being the cone surface. We show that the designed accuracy is being reached for smooth problems, that high-resolution is being attained for non-smooth problems, and that our modifications produce similar results while providing a more generic way to treat solid walls.

Local discontinuous Galerkin methods with generalized alternating numerical fluxes coupled with implicit–explicit time marching for solving convection–diffusion problems is analyzed in this paper, where the explicit part is treated by a strong-stability-preserving Runge–Kutta scheme, and the implicit part is treated by an L-stable diagonally implicit Runge–Kutta method. Based on the generalized alternating numerical flux, we establish the important relationship between the gradient and interface jump of the numerical solution with the independent numerical solution of the gradient, which plays a key role in obtaining the unconditional stability of the proposed schemes. Also by the aid of the generalized Gauss–Radau projection, optimal error estimates can be shown. Numerical experiments are given to verify the stability and accuracy of the proposed schemes with different numerical fluxes.

In this paper, we propose, analyze, and implement a new reduced basis method (RBM) tailored for the linear (ordinary and partial) differential equations driven by arbitrary (i.e. not necessarily Gaussian) types of noise. There are four main ingredients of our algorithm. First, we propose a new space-time-like treatment of time in the numerical schemes for ODEs and PDEs. The second ingredient is an accurate yet efficient compression technique for the spatial component of the space-time snapshots that the RBM is adopting as bases. The third ingredient is a non-conventional “parameterization” of a non-parametric problem. The last is a RBM that is free of any dedicated offline procedure yet is still efficient online. The numerical experiments verify the effectiveness and robustness of our algorithms for both types of differential equations.

This paper proposes and analyzes arbitrarily high-order discontinuous Galerkin (DG) and finite volume methods which provably preserve the positivity of density and pressure for the ideal MHD on general meshes. Unified auxiliary theories are built for rigorously analyzing the positivity-preserving (PP) property of MHD schemes with a HLL type flux on polytopal meshes in any space dimension. The main challenges overcome here include establishing relation between the PP property and discrete divergence of magnetic field on general meshes, and estimating proper wave speeds in the HLL flux to ensure the PP property. In 1D case, we prove that the standard DG and finite volume methods with the proposed HLL flux are PP, under condition accessible by a PP limiter. For multidimensional conservative MHD system, standard DG methods with a PP limiter are not PP in general, due to the effect of unavoidable divergence-error. We construct provably PP high-order DG and finite volume schemes by proper discretization of symmetrizable MHD system, with two divergence-controlling techniques: locally divergence-free elements and a penalty term. The former leads to zero divergence within each cell, while the latter controls the divergence error across cell interfaces. Our analysis reveals that a coupling of them is important for positivity preservation, as they exactly contribute the discrete divergence-terms absent in standard DG schemes but crucial for ensuring the PP property. Numerical tests confirm the PP property and the effectiveness of proposed PP schemes. Unlike conservative MHD system, the exact smooth solutions of symmetrizable MHD system are proved to retain the positivity even if the divergence-free condition is not satisfied. Our analysis and findings further the understanding, at both discrete and continuous levels, of the relation between the PP property and the divergence-free constraint., Comment: 49 pages, 11 figures

In this paper, we study the superconvergence properties of the energy-conserving discontinuous Galerkin (DG) method in [18] for one-dimensional linear hyperbolic equations. We prove the approximate solution superconverges to a particular projection of the exact solution. The order of this superconvergence is proved to be $$k+2$$ when piecewise $$\mathbb {P}^k$$ polynomials with $$k \ge 1$$ are used. The proof is valid for arbitrary non-uniform regular meshes and for piecewise $$\mathbb {P}^k$$ polynomials with arbitrary $$k \ge 1$$ . Furthermore, we find that the derivative and function value approximations of the DG solution are superconvergent at a class of special points, with an order of $$k+1$$ and $$k+2$$ , respectively. We also prove, under suitable choice of initial discretization, a ( $$2k+1$$ )-th order superconvergence rate of the DG solution for the numerical fluxes and the cell averages. Numerical experiments are given to demonstrate these theoretical results.

In this paper, we shall establish the superconvergence property of the Runge–Kutta discontinuous Galerkin (RKDG) method for solving a linear constant-coefficient hyperbolic equation. The RKDG method is made of the discontinuous Galerkin (DG) scheme with upwind-biased numerical fluxes coupled with the explicit Runge–Kutta algorithm of arbitrary orders and stages. Superconvergence results for the numerical flux, cell averages as well as the solution and derivative at some special points are shown, which are based on a systematical study of the $$\hbox {L}^2$$ -norm stability for the RKDG method and the incomplete correction techniques for the well-defined reference functions at each time stage. The result demonstrates that the superconvergence property of the semi-discrete DG method is preserved, and the optimal order in time is provided under the smoothness assumption that is independent of the number of stages. As a byproduct of the above superconvergence study, the expected order of the post-processed solution is obtained when a special initial solution is used. Some numerical experiments are also given.

In this paper, an ultra-weak discontinuous Galerkin (DG) method is developed to solve the generalized stochastic Korteweg–de Vries (KdV) equations driven by a multiplicative temporal noise. This method is an extension of the DG method for purely hyperbolic equations and shares the advantage and flexibility of the DG method. Stability is analyzed for the general nonlinear equations. The ultra-weak DG method is shown to admit the optimal error of order $$k+1$$ in the sense of the spatial $$L^2(0,2\pi )$$-norm for semi-linear stochastic equations, when polynomials of degree $$k\ge 2$$ are used in the spatial discretization. A second order implicit–explicit derivative-free time discretization scheme is also proposed for the matrix-valued stochastic ordinary differential equations derived from the spatial discretization. Numerical examples using Monte Carlo simulation are provided to illustrate the theoretical results.

In this paper, we extend our previous work in Huang and Shu (J Sci Comput, 2018. https://doi.org/10.1007/s10915-018-0852-1 ) and develop a third-order unconditionally positivity-preserving modified Patankar Runge–Kutta method for production–destruction equations. The necessary and sufficient conditions for the method to be of third-order accuracy are derived. With the same approach as Huang and Shu (2018), this time integration method is then generalized to solve a class of ODEs arising from semi-discrete schemes for PDEs and coupled with the positivity-preserving finite difference weighted essentially non-oscillatory schemes for non-equilibrium flows. Numerical experiments are provided to demonstrate the performance of our proposed scheme.

We present this special issue of the Journal of Scientific Computing to celebrate Bernardo Cockburn’s sixtieth birthday. The theme of this issue is discontinuous Galerkin methods, a hallmark of Bernardo’s distinguished professional career. This foreword provides an informal but rigorous account of what enabled Bernardo’s achievements, based on the concluding presentation he gave at the the IMA workshop “Recent Advances and Challenges in Discontinuous Galerkin Methods and Related Approaches” on July 1, 2017 which was widely deemed as the best lecture of his career so far.

In this paper, we construct a family of modified Patankar Runge–Kutta methods, which is conservative and unconditionally positivity-preserving, for production–destruction equations, and derive necessary and sufficient conditions to obtain second-order accuracy. This ordinary differential equation solver is then extended to solve a class of semi-discrete schemes for PDEs. Combining this time integration method with the positivity-preserving finite difference weighted essentially non-oscillatory (WENO) schemes, we successfully obtain a positivity-preserving WENO scheme for non-equilibrium flows. Various numerical tests are reported to demonstrate the effectiveness of the methods.

In this paper, we continue our research on the numerical study of convergence to steady-state solutions for a new class of finite volume weighted essentially non-oscillatory (WENO) schemes in Zhu and Shu (J Comput Phys 349:80–96, 2017), from tensor product meshes to triangular meshes. For the case of triangular meshes, this new class of finite volume WENO schemes was designed for time-dependent conservation laws in Zhu and Qiu (SIAM J Sci Comput 40(2):A903–A928, 2018) for the third- and fourth-order versions. In this paper, we extend the design to a new fifth-order version in the same framework to keep the essentially non-oscillatory property near discontinuities. Similar to the case of tensor product meshes in Zhu and Shu (2017), by performing such spatial reconstruction procedures together with a TVD Runge–Kutta time discretization, these WENO schemes do not suffer from slight post-shock oscillations that are responsible for the phenomenon wherein the residues of classical WENO schemes hang at a truncation error level instead of converging to machine zero. The third-, fourth-, and fifth-order finite volume WENO schemes in this paper can suppress the slight post-shock oscillations and have their residues settling down to a tiny number close to machine zero in steady-state simulations in our extensive numerical experiments.

So far the theory and numerical practice of stochastic partial differential equations (SPDEs) have dealt almost exclusively with Gaussian noise or Levy noise. Recently, Mikulevicius and Rozovskii (Stoch Partial Differ Equ Anal Comput 4:319–360, 2016) proposed a distribution-free Skorokhod–Malliavin calculus framework that is based on generalized stochastic polynomial chaos expansion, and is compatible with arbitrary driving noise. In this paper, we conduct systematic investigation on numerical results of these newly developed distribution-free SPDEs, exhibiting the efficiency of truncated polynomial chaos solutions in approximating moments and distributions. We obtain an estimate for the mean square truncation error in the linear case. The theoretical convergence rate, also verified by numerical experiments, is exponential with respect to polynomial order and cubic with respect to number of random variables included.

We further investigate the high order positivity-preserving discontinuous Galerkin (DG) methods for linear hyperbolic and radiative transfer equations developed in Yuan et al. (SIAM J Sci Comput 38:A2987–A3019, 2016). The DG methods in Yuan et al. (2016) can maintain positivity and high order accuracy, but they rely both on the scaling limiter in Zhang and Shu (J Comput Phys 229:8918–8934, 2010) and a rotational limiter, the latter may alter cell averages of the unmodulated DG scheme, thereby affecting conservation. Even though a Lax–Wendroff type theorem is proved in Yuan et al. (2016), guaranteeing convergence to weak solutions with correct shock speed when such rotational limiter is applied, it would still be desirable if a conservative DG method without changing the cell averages can be obtained which has both high order accuracy and positivity-preserving capability. In this paper, we develop and analyze such a DG method for both linear hyperbolic equations and radiative transfer equations. In the one-dimensional case, the method uses traditional DG space $$P^k$$ of piecewise polynomials of degree at most k. A key result is proved that the unmodulated DG solver in this case can maintain positivity of the cell average if the inflow boundary value and the source term are both positive, therefore the positivity-preserving framework in Zhang and Shu (2010) can be used to obtain a high order conservative positivity-preserving DG scheme. Unfortunately, in two-dimensions this is no longer the case. We show that the unmodulated DG solver based either on $$P^k$$ or $$Q^k$$ spaces (piecewise kth degree polynomials or piecewise tensor-product kth degree polynomials) could generate negative cell averages. We augment the DG space with additional functions so that the positivity of cell averages from the unmodulated DG solver can be restored, thereby leading to high order conservative positivity-preserving DG scheme based on these augmented DG spaces following the framework in Zhang and Shu (2010). Computational results are provided to demonstrate the good performance of our DG schemes.

In this article, we present a second-order in time implicit–explicit (IMEX) local discontinuous Galerkin (LDG) method for computing the Cahn–Hilliard equation, which describes the phase separation phenomenon. It is well-known that the Cahn–Hilliard equation has a nonlinear stability property, i.e., the free-energy functional decreases with respect to time. The discretized Cahn–Hilliard system modeled by the IMEX LDG method can inherit the nonlinear stability of the continuous model. We apply a stabilization technique and prove unconditional energy stability of our scheme. Numerical experiments are performed to validate the analysis. Computational efficiency can be significantly enhanced by using this IMEX LDG method with a large time step.

In this paper, we propose a simple bound-preserving sweeping procedure for conservative numerical approximations. Conservative schemes are of importance in many applications, yet for high order methods, the numerical solutions do not necessarily satisfy maximum principle. This paper constructs a simple sweeping algorithm to enforce the bound of the solutions. It has a very general framework acting as a postprocessing step accommodating many point-based or cell average-based discretizations. The method is proven to preserve the bounds of the numerical solution while conserving a prescribed quantity designated as a weighted average of values from all points. The technique is demonstrated to work well with a spectral method, high order finite difference and finite volume methods for scalar conservation laws and incompressible flows. Extensive numerical tests in 1D and 2D are provided to verify the accuracy of the sweeping procedure.

In this paper, we apply the local discontinuous Galerkin (LDG) method to 2D Keller---Segel (KS) chemotaxis model. We improve the results upon (Epshteyn and Kurganov in SIAM J Numer Anal, 47:368---408, 2008) and give optimal rate of convergence under special finite element spaces before the blow-up occurs (the exact solutions are smooth). Moreover, to construct physically relevant numerical approximations, we consider $$P^1$$P1 LDG scheme and develop a positivity-preserving limiter to the scheme, extending the idea in Zhang and Shu (J Comput Phys, 229:8918---8934, 2010). With this limiter, we can prove the $$L^1$$L1-stability of the numerical scheme. Numerical experiments are performed to demonstrate the good performance of the positivity-preserving LDG scheme. Moreover, it is known that the chemotaxis model will yield blow-up solutions under certain initial conditions. We numerically demonstrate how to find the approximate blow-up time by using the $$L^2$$L2-norm of the $$L^1$$L1-stable numerical solution.

In this paper, high order central finite difference schemes in a finite interval are analyzed for the diffusion equation. Boundary conditions of the initial-boundary value problem are treated by the simplified inverse Lax–Wendroff procedure. For the fully discrete case, a third order explicit Runge–Kutta method is used as an example for the analysis. Stability is analyzed by both the Gustafsson, Kreiss and Sundstrom theory and the eigenvalue visualization method on both semi-discrete and fully discrete schemes. The two different analysis techniques yield consistent results. Numerical tests are performed to demonstrate and validate the analysis results.

We consider the drift-diffusion (DD) model of one dimensional semiconductor devices, which is a system involving not only first derivative convection terms but also second derivative diffusion terms and a coupled Poisson potential equation. Optimal error estimates are obtained for both the semi-discrete and fully discrete local discontinuous Galerkin (LDG) schemes with smooth solutions. In the fully discrete scheme, we couple the implicit-explicit (IMEX) time discretization with the LDG spatial discretization, in order to allow larger time steps and to save computational cost. The main technical difficulty in the analysis is to treat the inter-element jump terms which arise from the discontinuous nature of the numerical method and the nonlinearity and coupling of the models. A simulation is also performed to validate the analysis.

In this paper, we develop and analyze a new multiscale discontinuous Galerkin (DG) method for one-dimensional stationary Schrodinger equations with open boundary conditions which have highly oscillating solutions. Our method uses a smaller finite element space than the WKB local DG method proposed in Wang and Shu (J Comput Phys 218:295---323, 2006) while achieving the same order of accuracy with no resonance errors. We prove that the DG approximation converges optimally with respect to the mesh size $$h$$h in $$L^2$$L2 norm without the typical constraint that $$h$$h has to be smaller than the wave length. Numerical experiments were carried out to verify the second order optimal convergence rate of the method and to demonstrate its ability to capture oscillating solutions on coarse meshes in the applications to Schrodinger equations.

Fourier spectral methods achieve exponential accuracy both on the approximation level and for solving partial differential equations, if the solution is analytic. If the solution is discontinuous but piecewise analytic up to the discontinuities, Fourier spectral methods produce poor pointwise accuracy, but still maintain exponential accuracy after post-processing (Gottlieb and Shu in SIAM Rev 30:644---668, 1997) . In Chen and Shu (J Comput Appl Math 265:83---95, 2014), an extended technique is provided to recover exponential accuracy for functions which have end-point singularities, from the knowledge of point values on standard collocation points. In this paper, we develop a technique to recover exponential accuracy from the first $$N$$N Fourier coefficients of functions which are analytic in the open interval but have unbounded derivative singularities at end points. With this post-processing method, we are able to obtain exponential accuracy of spectral methods applied to linear transport equations involving such functions.

In this paper, we propose a multi-scale discontinuous Galerkin (DG) method for second-order elliptic problems with curvilinear unidirectional rough coefficients by choosing a special non-polynomial approximation space. The key ingredient of the method lies in the incorporation of the local oscillatory features of the differential operators into the approximation space so as to capture the multi-scale solutions without having to resolve the finest scales. The unidirectional feature of the rough coefficients allows us to construct the basis functions of the DG non-polynomial approximation space explicitly, thereby greatly increasing the algorithm efficiency. Detailed error estimates for two-dimensional second-order DG methods are derived, and a general guidance on how to construct such non-polynomial basis is discussed. Numerical examples are also presented to validate and demonstrate the effectiveness of the algorithm.

In this paper we present an error estimate for the explicit Runge-Kutta discontinuous Galerkin method to solve a linear hyperbolic equation in one dimension with discontinuous but piecewise smooth initial data. The discontinuous finite element space is made up of piecewise polynomials of arbitrary degree $$k\ge 1$$ k ? 1 , and time is advanced by the third order explicit total variation diminishing Runge-Kutta method under the standard CFL temporal-spatial condition. The $$L^2(\mathbb R \backslash \mathcal R _T)$$ L 2 ( R ? R T ) -norm error at the final time $$T$$ T is optimal in both space and time, where $$\mathcal R _T$$ R T is the pollution region due to the initial discontinuity with the width $$\mathcal O (\sqrt{T\beta }h^{1/2}\log (1/h))$$ O ( T β h 1 / 2 log ( 1 / h ) ) . Here $$h$$ h is the maximum cell length and $$\beta $$ β is the flowing speed. These results are independent of the time step and hold also for the semi-discrete discontinuous Galerkin method.

The convergence to steady state solutions of the Euler equations for weighted compact nonlinear schemes (WCNS) [Deng X. and Zhang H. (2000), J. Comput. Phys. 165, 22–44 and Zhang S., Jiang S. and Shu C.-W. (2008), J. Comput. Phys. 227, 7294–7321] is studied through numerical tests. Like most other shock capturing schemes, WCNS also suffers from the problem that the residue can not settle down to machine zero for the computation of the steady state solution which contains shock waves but hangs at the truncation error level. In this paper, the techniques studied in [Zhang S. and Shu. C.-W. (2007), J. Sci. Comput. 31, 273–305 and Zhang S., Jiang S and Shu. C.-W. (2011), J. Sci. Comput. 47, 216–238], to improve the convergence to steady state solutions for WENO schemes, are generalized to the WCNS. Detailed numerical studies in one and two dimensional cases are performed. Numerical tests demonstrate the effectiveness of these techniques when applied to WCNS. The residue of various order WCNS can settle down to machine zero for typical cases while the small post-shock oscillations can be removed.

The gas dynamics equations, coupled with a static gravitational field, admit the hydrostatic balance where the flux produced by the pressure is exactly canceled by the gravitational source term. Many astrophysical problems involve the hydrodynamical evolution in a gravitational field, therefore it is essential to correctly capture the effect of gravitational force in the simulations. Improper treatment of the gravitational force can lead to a solution which either oscillates around the equilibrium, or deviates from the equilibrium after a long time run. In this paper we design high order well-balanced finite difference WENO schemes to this system, which can preserve the hydrostatic balance state exactly and at the same time can maintain genuine high order accuracy. Numerical tests are performed to verify high order accuracy, well-balanced property, and good resolution for smooth and discontinuous solutions. The main purpose of the well-balanced schemes designed in this paper is to well resolve small perturbations of the hydrostatic balance state on coarse meshes. The more difficult issue of convergence towards such hydrostatic balance state from an arbitrary initial condition is not addressed in this paper.

High order path-conservative schemes have been developed for solving nonconservative hyperbolic systems in (Pares, SIAM J. Numer. Anal. 44:300---321, 2006; Castro et al., Math. Comput. 75:1103---1134, 2006; J. Sci. Comput. 39:67---114, 2009). Recently, it has been observed in (Abgrall and Karni, J. Comput. Phys. 229:2759---2763, 2010) that this approach may have some computational issues and shortcomings. In this paper, a modification to the high order path-conservative scheme in (Castro et al., Math. Comput. 75:1103---1134, 2006) is proposed to improve its computational performance and to overcome some of the shortcomings. This modification is based on the high order finite volume WENO scheme with subcell resolution and it uses an exact Riemann solver to catch the right paths at the discontinuities. An application to one-dimensional compressible two-medium flows of nonconservative or primitive Euler equations is carried out to show the effectiveness of this new approach.

The entropy solutions of the compressible Euler equations satisfy a minimum principle for the specific entropy (Tadmor in Appl Numer Math 2:211–219, 1986). First order schemes such as Godunov-type and Lax-Friedrichs schemes and the second order kinetic schemes (Khobalatte and Perthame in Math Comput 62:119–131, 1994) also satisfy a discrete minimum entropy principle. In this paper, we show an extension of the positivity-preserving high order schemes for the compressible Euler equations in Zhang and Shu (J Comput Phys 229:8918–8934, 2010) and Zhang et al. (J Scientific Comput, in press), to enforce the minimum entropy principle for high order finite volume and discontinuous Galerkin (DG) schemes.

We are interested in the deterministic computation of the transients for the Boltzmann-Poisson system describing electron transport in semiconductor devices. The main difficulty of such computation arises from the very high dimensions of the model, making it necessary to use relatively coarse meshes and hence requiring the numerical solver to be stable and to have good resolution under coarse meshes. In this paper we give a brief survey of the discontinuous Galerkin (DG) method, which is a finite element method using discontinuous piecewise polynomials as basis functions and numerical fluxes based on upwinding for stability, for solving the Boltzmann-Poisson system. In many situations, the deterministic DG solver can produce accurate solutions with equal or less CPU time than the traditional DSMC (Direct Simulation Monte Carlo) solvers. In order to make the presentation more concise and to highlight the main ideas of the algorithm, we use a simplified model to describe the details of the DG method. Sample simulation results on the full Boltzmann-Poisson system are also given.

In Zhang and Shu (J. Comput. Phys. 229:3091---3120, 2010), two of the authors constructed uniformly high order accurate finite volume and discontinuous Galerkin (DG) schemes satisfying a strict maximum principle for scalar conservation laws on rectangular meshes. The technique is generalized to positivity preserving (of density and pressure) high order DG or finite volume schemes for compressible Euler equations in Zhang and Shu (J. Comput. Phys. 229:8918---8934, 2010). The extension of these schemes to triangular meshes is conceptually plausible but highly nontrivial. In this paper, we first introduce a special quadrature rule which is exact for two-variable polynomials over a triangle of a given degree and satisfy a few other conditions, by which we can construct high order maximum principle satisfying finite volume schemes (e.g. essentially non-oscillatory (ENO) or weighted ENO (WENO) schemes) or DG method solving two dimensional scalar conservation laws on triangular meshes. The same method can preserve the maximum principle for DG or finite volume schemes solving two-dimensional incompressible Euler equations in the vorticity stream-function formulation, or any passive convection equation with an incompressible velocity field. We also obtain positivity preserving (for density and pressure) high order DG or finite volume schemes solving compressible Euler equations on triangular meshes. Numerical tests for the third order Runge-Kutta DG (RKDG) method on unstructured meshes are reported.

The convergence to steady state solutions of the Euler equations for high order weighted essentially non-oscillatory (WENO) finite difference schemes with the Lax-Friedrichs flux splitting (Jiang and Shu, in J. Comput. Phys. 126:202---228, 1996) is investigated. Numerical evidence in Zhang and Shu (J. Sci. Comput. 31:273---305, 2007) indicates that there exist slight post-shock oscillations when we use high order WENO schemes to solve problems containing shock waves. Even though these oscillations are small in their magnitude and do not affect the "essentially non-oscillatory" property of the WENO schemes, they are indeed responsible for the numerical residue to hang at the truncation error level of the scheme instead of settling down to machine zero. Differently from the strategy adopted in Zhang and Shu (J. Sci. Comput. 31:273---305, 2007), in which a new smoothness indicator was introduced to facilitate convergence to steady states, in this paper we study the effect of the local characteristic decomposition on steady state convergence. Numerical tests indicate that the slight post-shock oscillation has a close relationship with the local characteristic decomposition process. When this process is based on an average Jacobian at the cell interface using the Roe average, as is the standard procedure for WENO schemes, such post-shock oscillation appears. If we instead use upwind-biased interpolation to approximate the physical variables including the velocity and enthalpy on the cell interface to compute the left and right eigenvectors of the Jacobian for the local characteristic decomposition, the slight post-shock oscillation can be removed or reduced significantly and the numerical residue settles down to lower values than other WENO schemes and can reach machine zero for many test cases. This new procedure is also effective for higher order WENO schemes and for WENO schemes with different smoothness indicators.

In this paper we continue our effort in Liu-Shu (2004) and Liu-Shu (2007) for developing local discontinuous Galerkin (LDG) finite element methods to discretize moment models in semiconductor device simulations. We consider drift-diffusion (DD) and high-field (HF) models of one-dimensional devices, which involve not only first derivative convection terms but also second derivative diffusion terms, as well as a coupled Poisson potential equation. Error estimates are obtained for both models with smooth solutions. The main technical difficulties in the analysis include the treatment of the inter-element jump terms which arise from the discontinuous nature of the numerical method, the nonlinearity, and the coupling of the models. A simulation is also performed to validate the analysis.

This note aims at demonstrating the advantage of moving-water well-balanced schemes over still-water well-balanced schemes for the shallow water equations. We concentrate on numerical examples with solutions near a moving-water equilibrium. For such examples, still-water well-balanced methods are not capable of capturing the small perturbations of the moving-water equilibrium and may generate significant spurious oscillations, unless an extremely refined mesh is used. On the other hand, moving-water well-balanced methods perform well in these tests. The numerical examples in this note clearly demonstrate the importance of utilizing moving-water well-balanced methods for solutions near a moving-water equilibrium.

A fifth order weighted essentially non-oscillatory (WENO) fast sweeping method is designed in this paper, extending the result of the third order WENO fast sweeping method in J. Sci. Comput. 29, 25---56 (2006) and utilizing the two approaches of accurate inflow boundary condition treatment in J. Comput. Math. 26, 1---11 (2008), which allows the usage of Cartesian meshes regardless of the domain boundary shape. The resulting method is tested on a variety of problems to demonstrate its good performance and CPU time efficiency when compared with lower order fast sweeping methods.

High order accurate weighted essentially non-oscillatory (WENO) schemes have been used extensively in numerical solutions of hyperbolic partial differential equations and other convection dominated problems. However the WENO procedure can not be applied directly to obtain a stable scheme when negative linear weights are present. In this paper, we first briefly review the WENO framework and the role of linear weights, and then present a detailed study on the positivity of linear weights in a few typical WENO procedures, including WENO interpolation, WENO reconstruction and WENO approximation to first and second derivatives, and WENO integration. Explicit formulae for the linear weights are also given for these WENO procedures. The results of this paper should be useful for future design of WENO schemes involving interpolation, reconstruction, approximation to first and second derivatives, and integration procedures.

In this paper, we develop a local discontinuous Galerkin (LDG) finite element method for surface diffusion and Willmore flow of graphs. We prove L 2 stability for the equation of surface diffusion of graphs and energy stability for the equation of Willmore flow of graphs. We provide numerical simulation results for different types of solutions of these two types of the equations to illustrate the accuracy and capability of the LDG method.

Strong stability preserving (SSP) high order time discretizations were developed to ensure nonlinear stability properties necessary in the numerical solution of hyperbolic partial differential equations with discontinuous solutions. SSP methods preserve the strong stability properties--in any norm, seminorm or convex functional--of the spatial discretization coupled with first order Euler time stepping. This paper describes the development of SSP methods and the connections between the timestep restrictions for strong stability preservation and contractivity. Numerical examples demonstrate that common linearly stable but not strong stability preserving time discretizations may lead to violation of important boundedness properties, whereas SSP methods guarantee the desired properties provided only that these properties are satisfied with forward Euler timestepping. We review optimal explicit and implicit SSP Runge---Kutta and multistep methods, for linear and nonlinear problems. We also discuss the SSP properties of spectral deferred correction methods.

In this paper, we develop a multiscale local discontinuous Galerkin (LDG) method to simulate the one-dimensional stationary Schrodinger-Poisson problem. The stationary Schrodinger equation is discretized by the WKB local discontinuous Galerkin (WKB-LDG) method, and the Poisson potential equation is discretized by the minimal dissipation LDG (MD-LDG) method. The WKB-LDG method we propose provides a significant reduction of both the computational cost and memory in solving the Schrodinger equation. Compared with traditional continuous finite element Galerkin methodology, the WKB-LDG method has the advantages of the DG methods including their flexibility in h-p adaptivity and allowance of complete discontinuity at element interfaces. Although not addressed in this paper, a major advantage of the WKB-LDG method is its feasibility for two-dimensional devices.

In this paper we develop a high order explicit finite difference weighted essentially non-oscillatory (WENO) scheme for solving a hierarchical size-structured population model with nonlinear growth, mortality and reproduction rates. The main technical complication is the existence of global terms in the coefficient and boundary condition for this model. We carefully design approximations to these global terms and boundary conditions to ensure high order accuracy. Comparing with the first order monotone and second order total variation bounded schemes for the same model, the high order WENO scheme is more efficient and can produce accurate results with far fewer grid points. Numerical examples including one in computational biology for the evolution of the population of Gambussia affinis, are presented to illustrate the good performance of the high order WENO scheme.

The convergence to steady state solutions of the Euler equations for the fifth-order weighted essentially non-oscillatory (WENO) finite difference scheme with the Lax---Friedrichs flux splitting [7, (1996) J. Comput. Phys. 126, 202---228.] is studied through systematic numerical tests. Numerical evidence indicates that this type of WENO scheme suffers from slight post-shock oscillations. Even though these oscillations are small in magnitude and do not affect the "essentially non-oscillatory" property of WENO schemes, they are indeed responsible for the numerical residue to hang at the truncation error level of the scheme instead of settling down to machine zero. We propose a new smoothness indicator for the WENO schemes in steady state calculations, which performs better near the steady shock region than the original smoothness indicator in [7, (1996) J. Comput. Phys. 126, 202---228.]. With our new smoothness indicator, the slight post-shock oscillations are either removed or significantly reduced and convergence is improved significantly. Numerical experiments show that the residue for the WENO scheme with this new smoothness indicator can converge to machine zero for one and two dimensional (2D) steady problems with strong shock waves when there are no shocks passing through the domain boundaries.

In this note we consider the application of the WENO scheme to simulations of steady-state flow in a converging diverging nozzle. We demonstrate the recovery of design accuracy through Gegenbauer postprocessing, despite the degradation of the order of accuracy for the numerical solution of the Euler equations to first-order in regions where the characteristics passed through the shock. We have shown a case in which the Gegenbauer postprocessing can recover the order of accuracy right up to the shock location. This suggests that high-order accurate information which crosses through the shock may not be irretrievably lost, and we can strive to recover it through various types of postprocessing.