Strong Stability Preserving Runge Kutta Research Articles

A fast explicit time-splitting spectral scheme for the viscous Cahn-Hilliard equation with nonlocal diffusion operator

The viscous Cahn-Hilliard equation (VCH) is a parameter-dependent model encapsulated Cahn-Hilliard (CH) model and constrained Allen-Cahn (CAC) model. In this paper, we propose a fast explicit time-splitting spectral scheme for the VCH with nonlocal diffusion operator. The numerical scheme is constructed based on the second-order symmetric time-splitting method, which results in a splitting of the original problem into nonlocal linear part and nonlinear part, respectively. The nonlocal linear subproblem can be solved exactly by applying the Fourier spectral discretization in space and thus has no stability restriction on the time-step size. The nonlinear subproblem is solved via a second-order strong stability preserving Runge–Kutta method. The new scheme is time symmetric, fully explicit and conserves the discrete mass exactly. Theoretical analysis shows that the stability of the scheme depends on both the viscosity coefficient and the range of nonlocal interactions. In addition, a rigorously convergence analysis of the scheme is established, and the convergence rate is proved to be second-order in time and spectral-order in space, respectively. Four numerical experiments are performed to demonstrate the effectiveness of the proposed scheme.

High order conservative Lagrangian schemes for two-dimensional radiation hydrodynamics equations in the equilibrium-diffusion limit

Radiation hydrodynamics equations (RHE) refer to the study of how interactions between radiation and matter influence thermodynamic states and dynamic flow, which has been widely applied to high temperature hydrodynamics, such as inertial confinement fusion (ICF) and astrophysical gaseous stars. Solving RHE accurately and robustly even under the equilibrium diffusion approximation is a challenging task. To address this, we develop two types of high order conservative Lagrangian schemes for RHE in the equilibrium-diffusion limit for the two dimensional case on the Lagrangian moving mesh. Based on the multi-resolution WENO reconstruction for the spatial discretization and strong stability preserving Runge-Kutta (SSP-RK) time discretization, we first develop an explicit Lagrangian scheme with the HLLC numerical flux to achieve high order accuracy in space and time. We also discuss the positivity-preserving property of the high order explicit Lagrangian scheme. To overcome the severe time step restriction arising from the nonlinear radiation diffusion term in the explicit scheme, we further present a high order explicit-implicit-null (EIN) Lagrangian scheme. By adding a sufficiently large linear diffusion term on both sides of the scheme, we treat the complicated nonlinear parts explicitly and efficiently, and treat the added linear diffusion term on the right-hand side implicitly with a relaxed time step restriction. According to our numerical experiments, these two types of Lagrangian schemes are high order accurate, conservative and can capture the interfaces automatically. Additionally, the explicit scheme is found to be non-oscillatory and can preserve positivity while maintaining the original high order accuracy.

CFL Optimized Forward–Backward Runge–Kutta Schemes for the Shallow-Water Equations

Abstract We present the formulation and optimization of a Runge–Kutta-type time-stepping scheme for solving the shallow-water equations, aimed at substantially increasing the effective allowable time step over that of comparable methods. This scheme, called FB-RK(3,2), uses weighted forward–backward averaging of thickness data to advance the momentum equation. The weights for this averaging are chosen with an optimization process that employs a von Neumann–type analysis, ensuring that the weights maximize the admittable Courant number. Through a simplified local truncation error analysis and numerical experiments, we show that the method is at least second-order in time for any choice of weights and exhibits low dispersion and dissipation errors for well-resolved waves. Further, we show that an optimized FB-RK(3,2) can take time steps up to 2.8 times as large as a popular three-stage, third-order strong stability-preserving Runge–Kutta method in a quasi-linear test case. In fully nonlinear shallow-water test cases relevant to oceanic and atmospheric flows, FB-RK(3,2) outperforms SSPRK3 in admittable time step by factors roughly between 1.6 and 2.2, making the scheme approximately twice as computationally efficient with little to no effect on solution quality. Significance Statement The purpose of this work is to develop and optimize time-stepping schemes for models relevant to oceanic and atmospheric flows. Specifically, for the shallow-water equations we optimize for schemes that can take time steps as large as possible while retaining solution quality. We find that our optimized schemes can take time steps between 1.6 and 2.2 times larger than schemes that cost the same number of floating point operations, translating directly to a corresponding speedup. Our ultimate goal is to use these schemes in climate-scale simulations.

BE-BDF2 Time Integration Scheme Equipped with Richardson Extrapolation for Unsteady Compressible Flows

In this work we investigate the effectiveness of the Backward Euler-Backward Differentiation Formula (BE-BDF2) in solving unsteady compressible inviscid and viscous flows. Furthermore, to improve its accuracy and its order of convergence, we have equipped this time integration method with the Richardson Extrapolation (RE) technique. The BE-BDF2 scheme is a second-order accurate, A-stable, L-stable and self-starting scheme. It has two stages: the first one is the simple Backward Euler (BE) and the second one is a second-order Backward Differentiation Formula (BDF2) that uses an intermediate and a past solution. The RE is a very simple and powerful technique that can be used to increase the order of accuracy of any approximation process by eliminating the lowest order error term(s) from its asymptotic error expansion. The spatial approximation of the governing Navier–Stokes equations is performed with a high-order accurate discontinuous Galerkin (dG) method. The presented numerical results for canonical test cases, i.e., the isentropic convecting vortex and the unsteady vortex shedding behind a circular cylinder, aim to assess the performance of the BE-BDF2 scheme, in its standard version and equipped with RE, by comparing it with the ones obtained by using more classical methods, like the BDF2, the second-order accurate Crank–Nicolson (CN2) and the explicit third-order accurate Strong Stability Preserving Runge–Kutta scheme (SSP-RK3).

Explicit strong stability preserving second derivative multistep methods for the Euler and Navier–Stokes equations

In this paper, we develop the explicit second derivative multistep methods (SDMMs) for the Euler and Navier–Stokes equations. The SDMMs use the calculated results from the previous steps, which can improve computational efficiency. The third-order and fourth-order SDMMs temporal discretization methods are constructed based on the order conditions. Additionally, the strong stability preserving (SSP) condition for the SDMMs is proposed, which provides the corresponding parameters for optimal SSP SDMMs. For numerical experiments, the SDMMs can achieve the order of accuracy as designed on the smooth regions, and the optimal SSP coefficient for the SDMMs is shown to be roughly the same in theory as in the numerical results. Moreover, the SDMMs have approximately twice the computational efficiency compared with the third-order and fourth-order two-derivative Runge–Kutta methods, and the SDMMs are also more efficient than the classical three-stage third-order SSP Runge–Kutta method.

An efficient scheme for solving nonlinear generalized kuramoto-sivashinksy processes

In this work, we investigate the numerical solution of generalized Kuramoto-Sivashinksy (GKS) problems based on the collocation of the quantic B-spline (QBS) and high-order strong stability-preserving Runge–Kutta (SSPRK54) scheme. When considering nonlinear parts that lose real features, we address the issue without resorting to any transformations or linearization. The efficiency and accuracy of our proposed technique are evaluated using a variety of illustrative examples. The numerical results show that our approach captured the natural behaviour of the problems well and consumed less storage space.

A new high-order maximum-principle-preserving explicit Runge–Kutta method for the nonlocal Allen–Cahn equation

We develop the high-order maximum principle preserving schemes based on combining the scalar variable with the explicit strong stability preserving Runge–Kutta (eSSPRK) methods, which can be applied to a class of gradient flows. In this work, we use the methods to deal with the nonlocal Allen–Cahn equation. There are two remarkable properties of the scheme. On the one hand, by using the introduced scalar variable to control the nonlinear term, and adopting the eSSPRK method to temporal discretization, we only need to solve the linear equations with constant coefficients at each time level. On the other hand, we introduce the first-order approximation of the energy balance equation instead of the dynamic equation for the scalar variable, and we demonstrate that the proposed schemes are unconditionally energy stable. In particular, the discrete schemes are verified to preserve the maximum bound principle under the time-step restriction. Moreover, error estimates and the asymptotic compatibility of the proposed methods are proved, respectively. Numerical examples are performed to illustrate the efficiency and accuracy of the numerical schemes.

HOW-MHD: A High-order WENO-based Magnetohydrodynamic Code with a High-order Constrained Transport Algorithm for Astrophysical Applications

Due to the prevalence of magnetic fields in astrophysical environments, magnetohydrodynamic (MHD) simulation has become a basic tool for studying astrophysical fluid dynamics. To further advance the precision of MHD simulations, we have developed a new simulation code that solves ideal adiabatic or isothermal MHD equations with high-order accuracy. The code is based on the finite-difference weighted essentially nonoscillatory (WENO) scheme and the strong stability-preserving Runge–Kutta (SSPRK) method. Most of all, the code implements a newly developed, high-order constrained transport (CT) algorithm for the divergence-free constraint of magnetic fields, completing its high-order competence. In this paper, we present the version in Cartesian coordinates, which includes a fifth-order WENO and a fourth-order five-stage SSPRK, along with extensive tests. With the new CT algorithm, fifth-order accuracy is achieved in convergence tests involving the damping of MHD waves in 3D space. And substantially improved results are obtained in magnetic loop advection and magnetic reconnection tests, indicating a reduction in numerical diffusivity. In addition, the reliability and robustness of the code, along with its high accuracy, are demonstrated through several tests involving shocks and complex flows. Furthermore, tests of turbulent flows reveal the advantages of high-order accuracy and show that the adiabatic and isothermal codes have similar accuracy. With its high-order accuracy, our new code would provide a valuable tool for studying a wide range of astrophysical phenomena that involve MHD processes.

A new type of improved third order WENO scheme with finite difference framework

In this paper, we present a novel variety of third order finite difference weighted essentially non-oscillatory (WENO) scheme based on WENO-ZQ scheme to approximate the computational solutions of one- and two-dimensional non-linear hyperbolic conservation laws. The new reconstruction process of fluxes is the convex combination of a second degree polynomial and two linear polynomials. We apply a new modified global smoothness indicator into the non-linear weights. We also use the power parameter ‘p’ to balance the accuracy of the proposed scheme. This novel third order WENO scheme applies three points data {xi−1,xi,xi+1} same as WENO-JS3 and WENO-Z3 schemes. In the proposed scheme, the main advantages are that the associated linear weights can be taken as any positive numbers on condition their summation is one. Also, it is easy to be used in multi-dimensional engineering problems. We implement the time discretization by using third order accurate strong stability preserving Runge–Kutta method. The extensive computational experiments of a collection of benchmark test problems are provided to demonstrate high resolution, robustness and absolute truncation errors in L1 and L∞ norms. We observe from the results that the new scheme contains better capability in performance than WENO-PZ3 scheme.

A well‐balanced and entropy stable scheme for a reduced blood flow model

AbstractA well‐known reduced model of the flow of blood in arteries can be formulated as a strictly hyperbolic system of two scalar balance laws in one space dimension where the unknowns are the cross‐sectional area of the artery and the average blood flow velocity as functions of the axial coordinate and time. This system is endowed with an entropy pair such that solutions of the balance equations satisfy an entropy inequality in the distributional sense. It is demonstrated that this property can be utilized to construct an entropy stable finite difference scheme for the blood flow model based on the general framework by Tadmor (E. Tadmor, Math. Comput., 49 (1987), 91–103). Furthermore, a fourth‐order extension of the resulting entropy conservative flux and a fourth‐order sign‐preserving reconstruction of the scaled entropy variables are employed as well as a second‐order strong stability preserving Runge–Kutta method for time discretization. The result is a computationally inexpensive and easy‐to‐implement explicit entropy stable scheme for the blood flow model. It is proven that the scheme is well‐balanced (i.e., preserves certain steady solutions of the model) and numerical examples are presented.

Strong Stability Preserving Runge–Kutta and Linear Multistep Methods

This paper reviews strong stability preserving discrete variable methods for differential systems. The strong stability preserving Runge–Kutta methods have been usually investigated in the literature on the subject, using the so-called Shu–Osher representation of these methods, as a convex combination of first-order steps by forward Euler method. In this paper, we revisit the analysis of strong stability preserving Runge–Kutta methods by reformulating these methods as a subclass of general linear methods for ordinary differential equations, and then using a characterization of monotone general linear methods, which was derived by Spijker in his seminal paper (SIAM J Numer Anal 45:1226–1245, 2007). Using this new approach, explicit and implicit strong stability preserving Runge–Kutta methods up to the order four are derived. These methods are equivalent to explicit and implicit RK methods obtained using Shu–Osher or generalized Shu–Osher representation. We also investigate strong stability preserving linear multistep methods using again monotonicity theory of Spijker.

A second-order Strang splitting scheme with exponential integrating factor for the Allen–Cahn equation with logarithmic Flory–Huggins potential

In this paper, we mainly consider the numerical approximation for the Allen-Cahn(AC) equation with logarithmic Flory–Huggins potential. It is well-known that the logarithmic Flory–Huggins potential is more physically significant than the generally employed fourth-order polynomial potential in various phase-field models. Nevertheless, owing to the singularity of the logarithmic term and the strong nonlinearity of the energy potential function, it is extremely challenging to develop an effective and easy-to-implement numerical scheme that preserves the maximum principle and energy dissipation law. To overcome these difficulties, we construct a second-order Strang splitting scheme with exponential integrating factor. Firstly, base on the idea of dimensional splitting, a new time dependent function, which is called exponential integrating factor is introduced. Then we propose the Strang splitting approach which is aim to decompose the original equation into linear part and nonlinear part. In particular, the logarithmic term is treated by explicit 2-stage strong stability preserving Runge–Kutta(SSP-RK2) method. Furthermore, we rigorously demonstrate the maximum principle, energy stability and convergence in details for the proposed scheme. A variety of numerical simulations in 2D and 3D are presented to confirm the validity of the proposed method.

The totally volume local discontinuous Galerkin method for viscous Burgers’ equations.

In this paper, the totally volume local discontinuous Galerkin TV-LDG method is proposed to solve the viscous Burgers’ equations with appropriate initial and boundary conditions. Time derivative is discretized by the third-order strong stability preserving Runge Kutta explicit SSP-RK (3, 3) method. Space derivatives discretization is performed by the totally volume local discontinuous Galerkin method. Finally, the validity of proposed scheme is demonstrated by numerical experiments and shows that the present scheme offers better accuracy in comparison with other existing numerical methods.

The positivity preserving property on the high order arbitrary Lagrangian-Eulerian discontinuous Galerkin method for Euler equations

This paper presents an almost arbitrary Lagrangian-Eulerian discontinuous Galerkin (ALE-DG) method to solve the compressible Euler equations in one and two space dimensions, and consider their positivity preserving property of states density and pressure. The ALE-DG method coupling with a modified strong stability-preserving Runge-Kutta method is developed to ensure the geometric conservation law and the positivity property of the scheme with the positivity preserving limiter developed by Zhang et al. (2010) [34] and Zhang et al. (2012) [38]. For the Lax-Friedrichs, HLL, and HLLC numerical fluxes, we prove that our proposed ALE-DG method can keep the positivity property of approximations of density and pressure from the first-order ALE-DG method and the cell averages of approximations of density and pressure from high order ALE-DG method under suitable time step. Numerical examples demonstrate that the proposed positivity-preserving ALE-DG method can keep the high order accuracy, numerical stability, and positivity on moving meshes for Euler equations.

A High-Order Velocity-Based Discontinuous Galerkin Scheme for the Shallow Water Equations: Local Conservation, Entropy Stability, Well-Balanced Property, and Positivity Preservation

The nonlinear shallow water equations (SWEs) are widely used to model the unsteady water flows in rivers and coastal areas. In this work, we present a novel class of locally conservative, entropy stable and well-balanced discontinuous Galerkin (DG) methods for the nonlinear shallow water equation with a non-flat bottom topography. The major novelty of our work is the use of velocity field as an independent solution unknown in the DG scheme, which is closely related to the entropy variable approach to entropy stable schemes for system of conservation laws proposed by Tadmor (in: Tezduyar, Hughes T (eds) Proceedings of the winter annual meeting of the American Society of Mechanical Engineering 1986) back in 1986, where recall that velocity is part of the entropy variable for the shallow water equations. Due to the use of velocity as an independent solution unknown, no specific numerical quadrature rules are needed to achieve entropy stability of our scheme on general unstructured meshes in two dimensions. The proposed DG semi-discretization is then carefully combined with the classical explicit strong stability preserving Runge–Kutta (SSP–RK) time integrators (Gottlieb et al. in SIAM Rev. 43, 89–112, 2001) to yield a locally conservative, well-balanced, and positivity preserving fully discrete scheme. Here the positivity preservation property is enforced with the help of a simple scaling limiter. In the fully discrete scheme, we re-introduce discharge as an auxiliary unknown variable. In doing so, standard slope limiting procedures can be applied on the conservative variables (water height and discharge) without violating the local conservation property. Here we apply a characteristic-wise TVB limiter (Cockburn and Shu in J Comput Phys 141:199–224, 1998) on the conservative variables using the Fu-Shu troubled cell indicator (Fu and Shu in J Comput Phys 347:305–327, 2017) in each inner stage of the Runge–Kutta time stepping to suppress numerical oscillations. This fully discrete can be readily applied to various SWEs simulations without dry areas where the water height is close to zero. The case with dry areas need further special attention, where the velocity approximation can be unphysically large near cells with a small water height, which may eventually crashes the simulation if no special treatment is used near these cells. Here we propose a simple wetting/drying treatment for the velocity update without violating the local conservation property to enhance the robustness of the overall scheme. One- and two-dimensional numerical experiments are presented to demonstrate the performance of the proposed methods.

Monolithic parabolic regularization of the MHD equations and entropy principles

We show at the PDE level that the monolithic parabolic regularization of the equations of ideal magnetohydrodynamics (MHD) is compatible with all the generalized entropies, fulfills the minimum entropy principle, and preserves the positivity of density and internal energy. We then numerically investigate this regularization for the MHD equations using continuous finite elements in space and explicit strong stability preserving Runge–Kutta methods in time. The artificial viscosity coefficient of the regularization term is constructed to be proportional to the entropy residual of MHD. It is shown that the method has a high order of accuracy for smooth problems and captures strong shocks and discontinuities accurately for non-smooth problems.

High-resolution ILW outflow boundary treatment for the Navier–Stokes equations

WENO schemes have high-resolution and great capability of capturing unsteady and non-linear phenomena. These phenomena occur in elaborated flows, such as cascade and external aerodynamic flows. Thus, high-resolution methods play an important role in modern CFD. Easiness and effectiveness are reached when implementing these methods in rectangular meshes. However, one has the non-boundary-conforming issue. To overcome it the ILW procedure can be employed. We showed that the ILW has a similar behavior to the far-field boundary treatment. Then, we proposed and tested a new ILW outflow boundary treatment and an WENO-type extrapolation. We employed the finite difference method with the positivity-preserving Lax–Friedrichs splitting, high-resolution viscous terms discretization, multi-resolution WENO scheme, and third-order strong stability preserving Runge–Kutta time discretization. We tested the proposed methods in a smooth test case, supersonic flow past a cylinder, NACA 0012 nitrogen flows with 0 and 12° angle of attack, and a NACA 9520 cascade airflow. The proposed methods have high order and resolution, allowed domain size reduction, and captured unsteady and non-linear phenomena. Although preliminary, the NACA 0012 and 9520 flows showed the methods are promising. Common issues with other methods are also reported and further improvements are needed. For example, local mesh refinement, high-order wall boundary treatment, and turbulence modeling.

An assessment of solvers for algebraically stabilized discretizations of convection–diffusion–reaction equations

Abstract We consider flux-corrected finite element discretizations of 3D convection-dominated transport problems and assess the computational efficiency of algorithms based on such approximations. The methods under investigation include flux-corrected transport schemes and monolithic limiters. We discretize in space using a continuous Galerkin method and ℙ1 or ℚ1 finite elements. Time integration is performed using the Crank–Nicolson method or an explicit strong stability preserving Runge–Kutta method. Nonlinear systems are solved using a fixed-point iteration method, which requires solution of large linear systems at each iteration or time step. The great variety of options in the choice of discretization methods and solver components calls for a dedicated comparative study of existing approaches. To perform such a study, we define new 3D test problems for time dependent and stationary convection–diffusion–reaction equations. The results of our numerical experiments illustrate how the limiting technique, time discretization and solver impact on the overall performance.

Embedded pairs for optimal explicit strong stability preserving Runge–Kutta methods

We construct a family of embedded pairs for optimal explicit strong stability preserving Runge–Kutta methods of order 2≤p≤4 to be used to obtain numerical solution of spatially discretized hyperbolic PDEs. In this construction, the goals include non-defective property, large stability region, and small error values as defined in Dekker and Verwer (1984) and Kennedy et al. (2000). The new family of embedded pairs offer the ability for strong stability preserving (SSP) methods to adapt by varying the step-size. Through several numerical experiments, we assess the overall effectiveness in terms of work versus precision while also taking into consideration accuracy and stability.

High-order multirate explicit time-stepping schemes for the baroclinic-barotropic split dynamics in primitive equations

In order to treat the multiple time scales of ocean dynamics in an efficient manner, the baroclinic-barotropic splitting technique has been widely used for solving the primitive equations for ocean modeling. Based on the framework of strong stability-preserving Runge-Kutta approach, we propose two high-order multirate explicit time-stepping schemes (SSPRK2-SE and SSPRK3-SE) for the resulting split system in this paper. The proposed schemes allow for a large time step to be used for the three-dimensional baroclinic (slow) mode and a small time step for the two-dimensional barotropic (fast) mode, in which each of the two mode solves just need to satisfy their respective CFL conditions for numerical stability. Specifically, at each time step, the baroclinic velocity is first computed by advancing the baroclinic mode and fluid thickness of the system with the large time-step and the assistance of some intermediate approximations of the barotropic mode obtained by substepping with the small time step; then the barotropic velocity is corrected by using the small time step to re-advance the barotropic mode under an improved barotropic forcing produced by interpolation of the forcing terms from the preceding baroclinic mode solves; lastly, the fluid thickness is updated by coupling the baroclinic and barotropic velocities. Additionally, numerical inconsistencies on the discretized sea surface height caused by the mode splitting are relieved via a reconciliation process with carefully calculated flux deficits. Two benchmark tests from the “MPAS-Ocean” platform are carried out to numerically demonstrate the performance and parallel scalability of the proposed SSPRK-SE schemes.