Implicit-explicit Runge-Kutta Schemes Research Articles

Stability analysis and error estimates of implicit-explicit Runge-Kutta least squares RBF-FD method for time-dependent parabolic equation

In this paper, for the time-dependent parabolic equations defined on complex geometries domain, we develop and analyze the least-squares radial basis function finite difference method (RBF-FD) coupled with the implicit-explicit Runge-Kutta (IMEX-RK) time discretization up to third order accuracy, which improves stability and accuracy. We derive the absolute stability region and the optimal time-step constraint for four kinds of IMEX-RK schemes. Compared to the traditional explicit or implicit time discretization, these are not trivial. Under a wide time-step constraint, the stability and the error estimates in l2-norm are established. Finally, several numerical experiments on the regular domain and non-convex domain are performed to validate the theoretical analysis.

Uniform accuracy of implicit-explicit Runge-Kutta (IMEX-RK) schemes for hyperbolic systems with relaxation

Implicit-explicit Runge-Kutta (IMEX-RK) schemes are popular methods to treat multiscale equations that contain a stiff part and a non-stiff part, where the stiff part is characterized by a small parameter ε \varepsilon . In this work, we prove rigorously the uniform stability and uniform accuracy of a class of IMEX-RK schemes for a linear hyperbolic system with stiff relaxation. The result we obtain is optimal in the sense that it holds regardless of the value of ε \varepsilon and the order of accuracy is the same as the design order of the original scheme, i.e., there is no order reduction.

General-relativistic Radiation Transport Scheme in Gmunu. I. Implementation of Two-moment-based Multifrequency Radiative Transfer and Code Tests

We present the implementation of a two-moment-based general-relativistic multigroup radiation transport module in the General-relativistic multigrid numerical (Gmunu) code. On top of solving the general-relativistic magnetohydrodynamics and the Einstein equations with conformally flat approximations, the code solves the evolution equations of the zeroth- and first-order moments of the radiations in the Eulerian-frame. An analytic closure relation is used to obtain the higher order moments and close the system. The finite-volume discretization has been adopted for the radiation moments. The advection in spatial space and frequency-space are handled explicitly. In addition, the radiation–matter interaction terms, which are very stiff in the optically thick region, are solved implicitly. The implicit–explicit Runge–Kutta schemes are adopted for time integration. We test the implementation with a number of numerical benchmarks from frequency-integrated to frequency-dependent cases. Furthermore, we also illustrate the astrophysical applications in hot neutron star and core-collapse supernovae modelings, and compare with other neutrino transport codes.

High order asymptotic preserving finite difference WENO schemes with constrained transport for MHD equations in all sonic Mach numbers

In this paper, a high-order semi-implicit (SI) asymptotic preserving (AP) and divergence-free finite difference weighted essentially nonoscillatory (WENO) scheme is proposed for magnetohydrodynamic (MHD) equations. We consider the sonic Mach number ε ranging from 0 to O(1). High-order accuracy in time is obtained by SI implicit-explicit Runge–Kutta (IMEX-RK) time discretization. High-order accuracy in space is achieved by finite difference WENO schemes with characteristic-wise reconstructions. A constrained transport method is applied to maintain a discrete divergence-free condition. We formally prove that the scheme is AP. Asymptotic accuracy (AA) in the incompressible MHD limit is obtained if the implicit part of the SI IMEX-RK scheme is stiffly accurate. Numerical experiments are provided to validate the AP, AA, and divergence-free properties of our proposed approach. Besides, the scheme can well capture discontinuities such as shocks in an essentially non-oscillatory fashion in the compressible regime, while it is also a good incompressible solver with uniform large-time step conditions in the low sonic Mach limit.

Unified preserving properties of kinetic schemes.

The kinetic theory provides a physical basis for developing multiscale methods for gas flows covering a wide range of flow regimes. A particular challenge for kinetic schemes is whether they can capture the correct hydrodynamic behaviors of the system in the continuum regime (i.e., as the Knudsen number ε≪1) without enforcing kinetic scale resolution. At the current stage, the main approach to analyze such a property is the asymptotic preserving (AP) concept, which aims to show whether a kinetic scheme reduces to a solver for the hydrodynamic equationsas ε→0, such as the shock capturing scheme for the Euler equations. However, the detailed asymptotic properties of the kinetic scheme are indistinguishable when ε is small but finite under the AP framework. To distinguish different characteristics of kinetic schemes, in this paper we introduce the concept of unified preserving (UP) aiming at assessing asymptotic orders of a kinetic scheme by employing the modified equationapproach and Chapman-Enskon analysis. It is shown that the UP properties of a kinetic scheme generally depend on the spatial and temporal accuracy and closely on the interconnections among three scales (kinetic scale, numerical scale, and hydrodynamic scale) and their corresponding coupled dynamics. Specifically, the numerical resolution and specific discretization of particle transport and collision determine the flow physics of the scheme in different regimes, especially in the near continuum limit. As two examples, the UP methodology is applied to analyze the discrete unified gas-kinetic scheme and a second-order implicit-explicit Runge-Kutta scheme in their asymptotic behaviors in the continuum limit.

An adaptive nonhydrostatic dynamical core using a multimoment finite‐volume method on a cubed sphere

Abstract: An adaptive nonhydrostatic atmospheric dynamical core was developed on a Cartesian grid and extended to spherical geometry with the application of a cubed‐sphere grid. To assure a practical Courant‐Friedrichs‐Lewy number for stability, a horizontally explicit and vertically implicit algorithm was applied in this model through a third‐order implicit–explicit Runge–Kutta scheme. A two‐dimensional adaptive grid was applied in the horizontal directions to generate the computational meshes, with variable resolutions according to the evolution of the predicted variables during the simulations, while the vertical grid is fixed considering the characteristics of atmospheric flows and the computational efficiency in practical applications. The Berger–Oliger block‐structured adaptive algorithm was adopted in this study. Blocks with different resolutions can be constructed straightforwardly on each patch of the cubed sphere and an algorithm to exchange both solution and block information between adjacent patches was designed in order to implement a global model. The nonhydrostatic governing equations are solved by a three‐point multimoment constrained finite‐volume scheme. Using the compact stencil for spatial reconstructions, the interpolation operations between coarse–fine blocks can be implemented efficiently and the local‐based scheme is also helpful to suppress the computational modes around the coarse–fine interfaces due to the abrupt change of grid resolution. Additionally, flux corrections were conducted along the block boundaries and the resulting dynamical core is rigorously conservative. The proposed model was evaluated by calculating several idealized benchmark tests, and the effectiveness of the adaptive model in saving computational costs was verified in this study.

An Extension of Gmunu: General-relativistic Resistive Magnetohydrodynamics Based on Staggered-meshed Constrained Transport with Elliptic Cleaning

We present the implementation of general-relativistic resistive magnetohydrodynamics solvers and three divergence-free handling approaches adopted in the General-relativistic multigrid numerical (Gmunu) code. In particular, implicit–explicit Runge–Kutta schemes are used to deal with the stiff terms in the evolution equations for small resistivity. The three divergence-free handling methods are (i) hyperbolic divergence cleaning (also known as the generalized Lagrange multiplier), (ii) staggered-meshed constrained transport schemes, and (iii) elliptic cleaning through a multigrid solver, which is applicable in both cell-centered and face-centered (stagger grid) magnetic fields. The implementation has been tested with a number of numerical benchmarks from special-relativistic to general-relativistic cases. We demonstrate that our code can robustly recover from the ideal magnetohydrodynamics limit to a highly resistive limit. We also illustrate the applications in modeling magnetized neutron stars, and compare how different divergence-free handling methods affect the evolution of the stars. Furthermore, we show that the preservation of the divergence-free condition of the magnetic field when using staggered-meshed constrained transport schemes can be significantly improved by applying elliptic cleaning.

TVD-MOOD schemes based on implicit-explicit time integration

The context of this work is the development of first order total variation diminishing (TVD) implicit-explicit (IMEX) Runge-Kutta (RK) schemes as a basis of a Multidimensional Optimal Order detection (MOOD) approach to approximate the solution of hyperbolic multi-scale equations. A key feature of our newly proposed TVD schemes is that the resulting CFL condition does not depend on the fast waves of the considered model, as long as they are integrated implicitly. However, a result from Gottlieb et al. [1] gives a first order barrier for unconditionally stable implicit TVD-RK schemes and TVD-IMEX-RK schemes with scale-independent CFL conditions. Therefore, the goal of this work is to consistently improve the resolution of a first-order IMEX-RK scheme, while retaining its L∞ stability and TVD properties. In this work we present a novel approach based on a convex combination between a first-order TVD IMEX Euler scheme and a potentially oscillatory high-order IMEX-RK scheme. We derive and analyse the TVD property for a scalar multi-scale equation and numerically assess the performance of our TVD schemes compared to standard L-stable and SSP IMEX RK schemes from the literature. Finally, the resulting TVD-MOOD schemes are applied to the isentropic Euler equations.

A nearly-conservative, high-order, forward Lagrange–Galerkin method for the resolution of compressible flows on unstructured triangular meshes

In this work, we present a novel Lagrange–Galerkin method for the resolution of compressible and inviscid flows. The scheme considers: (i) high-order continuous space discretizations on unstructured triangular meshes, (ii) high-order implicit–explicit Runge-Kutta schemes for the time discretization, (iii) conservation of mass, momentum and total energy, as long as some integrals in the formulation are computed exactly, and (iv) subgrid-stabilization and discontinuity-capturing operators based on Brenner's model [51] (2006) for viscous flows. The method has been tested on several benchmark problems using a fourth-order time-marching formula and up to fifth-order continuous finite elements, yielding the expected results both for smooth and discontinuous solutions.

On the Construction of Conservative Semi-Lagrangian IMEX Advection Schemes for Multiscale Time Dependent PDEs

This article is devoted to the construction of a new class of semi-Lagrangian (SL) schemes with implicit-explicit (IMEX) Runge-Kutta (RK) time stepping for PDEs involving multiple space-time scales. The semi-Lagrangian (SL) approach fully couples the space and time discretization, thus making the use of RK strategies particularly difficult to be combined with. First, a simple scalar advection-diffusion equation is considered as a prototype PDE for the development of a high order formulation of the semi-Lagrangian IMEX algorithms. The advection part of the PDE is discretized explicitly at the aid of a SL technique, while an implicit discretization is employed for the diffusion terms. In this way, an unconditionally stable numerical scheme is obtained, that does not suffer any CFL-type stability restriction on the maximum admissible time step. Second, the SL-IMEX approach is extended to deal with hyperbolic systems with multiple scales, including balance laws, that involve shock waves and other discontinuities. A conservative scheme is then crucial to properly capture the wave propagation speed and thus to locate the discontinuity and the plateau exhibited by the solution. A novel SL technique is proposed, which is based on the integration of the governing equations over the space-time control volume which arises from the motion of each grid point. High order of accuracy is ensured by the usage of IMEX RK schemes combined with a Cauchy–Kowalevskaya procedure that provides a predictor solution within each space-time element. The one-dimensional shallow water equations (SWE) are chosen to validate the new conservative SL-IMEX schemes, where convection and pressure fluxes are treated explicitly and implicitly, respectively. The asymptotic-preserving (AP) property of the novel schemes is also studied considering a relaxation PDE system for the SWE. A large suite of convergence studies for both the non-conservative and the conservative version of the novel class of methods demonstrates that the formal order of accuracy is achieved and numerical evidences about the conservation property are shown. The AP property for the corresponding relaxation system is also investigated.

A coupled high-order implicit-explicit flux reconstruction lattice Boltzmann method for nearly incompressible thermal flows

In this paper, based on the double distribution function (DDF) model and the Boussinesq approximation, a high-order implicit-explicit flux reconstruction thermal lattice Boltzmann method (FRTLBM) is proposed for simulating nearly incompressible thermal flows. Two discrete velocity Boltzmann equations (DVBEs), one governing the flow field and the other governing the temperature field, are solved by using a high-order flux reconstruction scheme for spatial discretization and a second-order implicit-explicit Runge-Kutta scheme (IMEX RK) for time discretization in the generalized curvilinear coordinates with uniform and non-uniform grids. Numerical validations of the proposed method are implemented by simulating (a) porous plate problem, (b) natural convection in a square cavity, (c) Rayleigh-Bénard convection, (d) natural convection in a concentric annulus and (e) mixed heat transfer from a heated circular cylinder. Numerical results demonstrate that the present method can achieve the high-order accuracy and it is stable and accurate even at the relatively high Rayleigh number (Ra = 107 and 108). The present method is efficient due to the small amount of mesh, the simple algorithm and the time step unlimited by the relaxation time. Since the rescontruction is conducted in a single cell, the present method is compact for parallel computing. Moreover, the present method has a good adaptability of complicated geometries.

A high-order implicit-explicit flux reconstruction lattice Boltzmann method for viscous incompressible flows

In this paper, a high-order implicit-explicit flux reconstruction lattice Boltzmann method (FRLBM) is proposed for simulating viscous incompressible flows. The discrete velocity Boltzmann equation (DVBE) is solved by using a high-order flux reconstruction scheme (FR) for spatial discretization and a high-order implicit-explicit Runge-Kutta scheme (IMEX RK) for time discretization in the generalized curvilinear coordinates with uniform and non-uniform grids to obtain high simulation efficiency. Numerical validations of the proposed method are implemented by simulating (a) Taylor-Green vortex problem, (b) doubly periodic shear layer flow, (c) lid-driven cavity flow, (d) cylindrical Couette flow and (e) flow over a circular cylinder. Numerical results demonstrate that the present method is stable and accurate even at high Reynolds number. Compared with the existing high-order off-lattice Boltzmann methods, the present method is more efficient. In addition, the present method holds the compactness of the standard LBM for parallel computing, but greatly improves the adaptability of complicated geometries.

A unified asymptotic preserving and well-balanced scheme for the Euler system with multiscale relaxation

The design and analysis of a unified asymptotic preserving (AP) and well-balanced scheme for the Euler Equations with gravitational and frictional source terms is presented in this paper. The asymptotic behaviour of the Euler system in the limit of zero Mach and Froude numbers, and large friction is characterised by an additional scaling parameter. Depending on the values of this parameter, the Euler system relaxes towards a hyperbolic or a parabolic limit equation. Standard Implicit–Explicit Runge–Kutta schemes are incapable of switching between these asymptotic regimes. We propose a time semi-discretisation to obtain a unified scheme which is AP for the two different limits. A further reformulation of the semi-implicit scheme can be recast as a fully-explicit method in which the mass update contains both hyperbolic and parabolic fluxes. A space–time fully-discrete scheme is derived using a finite volume framework. A hydrostatic reconstruction strategy, an upwinding of the sources at the interfaces, and a careful choice of the central discretisation of the parabolic fluxes are used to achieve the well-balancing property for hydrostatic steady states. Results of several numerical case studies are presented to substantiate the theoretical claims and to verify the robustness of the scheme.

High order modal Discontinuous Galerkin Implicit–Explicit Runge Kutta and Linear Multistep schemes for the Boltzmann model on general polygonal meshes

Deterministic solutions of the Boltzmann equation represent a real challenge due to the enormous computational effort which is required to produce such simulations and often stochastic methods such as Direct Simulation Monte Carlo (DSMC) are used instead due to their lower computational cost. In this work, we show that combining different technologies for the discretization of the velocity space and of the physical space coupled with suitable time integration techniques, it is possible to compute very precise deterministic approximate solutions of the Boltzmann model in different regimes, from extremely rarefied to dense fluids, with CFL conditions only driven by the hyperbolic transport term. To that aim, we develop modal Discontinuous Galerkin (DG) Implicit–Explicit Runge Kutta schemes (DG-IMEX-RK) and Implicit–Explicit Linear Multistep Methods based on Backward-Finite-Differences (DG-IMEX-BDF) for solving the Boltzmann model on multidimensional unstructured meshes. The solution of the Boltzmann collision operator is obtained through fast spectral methods, while the transport term in the governing equations is discretized relying on an explicit shock-capturing DG method on polygonal tessellations in the physical space. A novel class of WENO-type limiters, based on a shifting of the moments of inertia for each zone of the mesh, is used to control spurious oscillations of the DG solution across discontinuities. The use of Linear Multistep Methods (LMM) allows the Boltzmann solutions to be consistent not only with the compressible Euler limit but also with the Navier–Stokes asymptotic regime. In addition, as numerically proven, they also permit to strongly reduce the computational effort compared to Runge–Kutta approaches while maintaining the same or even larger accuracy. The performances of these different time discretization techniques are measured comparing both precision and efficiency. At the same time, comparisons against simpler relaxation type kinetic models such as the BGK model are proposed. The order of convergence is numerically measured for different regimes and found to agree with the theoretical findings. The new methods are validated considering two-dimensional benchmark test cases typically used in the fluid dynamics community. A prototype engineering problem consisting of a supersonic flow around a NACA 0012 airfoil with space–time-dependent boundary conditions is also presented for which the pressure coefficients are measured.

A compact Fourth-Order Implicit-Explicit Runge-Kutta Type Method for Solving Diffusive Lotka–Volterra System

This paper aims to developed a high-order and accurate method for the solution of one-dimensional Lotka-Volterra-diffusion with Numman boundary conditions. A fourth-order compact finite difference scheme for spatial part combined with implicit-explicit Runge Kutta scheme in temporal are proposed. Furthermore, boundary points are discretized by using a compact finite difference scheme in terms of fourth order accuracy. A key idea for proposed scheme is to take full advantage of method of line (MOL), this is consequently enabling us to use implicit-explicit Runge Kutta method, that are of fourth order in time. We constructed fourth order accuracy in both space and time and is unconditionally stable. This is consequently leading to a reduction in the computational cost of the scheme. Numerical experiments show that the combination of the compact finite difference with IMEX- RK methods give an accurate and reliable for solving the Lotka-Volterra-diffusion.

On the preserving of the maximum principle and energy stability of high-order implicit-explicit Runge-Kutta schemes for the space-fractional Allen-Cahn equation

We put forward and analyze the high-order (up to fourth) strong stability-preserving implicit-explicit Runge-Kutta schemes for the time integration of the space-fractional Allen-Cahn equation, which inherits the maximum principle preserving and energy stability. The space-fractional Allen-Cahn equation with homogeneous Dirichlet boundary condition is first discretized in the spatial direction by using a second-order fractional centered difference scheme that preserves the semi-discrete maximum principle. It is subsequently integrated in the temporal direction by a class of strong stability-preserving implicit-explicit Runge-Kutta schemes that are specifically designed to preserve the maximum principle to the optimal time step size. The convergence order in the discrete $L^{\infty }$ norm and energy boundedness are provided by using the established maximum principle. Finally, a series of numerical experiments are carried out to demonstrate the high-order convergence, maximum principle preserving, and energy stability of the proposed schemes.

Accelerated implicit-explicit Runge-Kutta schemes for locally stiff systems

In this paper we introduce a family of accelerated implicit-explicit (AIMEX) schemes for the solution of stiff systems of equations. Similar to conventional IMEX schemes, AIMEX schemes allow a problem to be split into its stiff and non-stiff constituent parts. These are then advanced in time using an implicit and explicit scheme, respectively. By design, AIMEX schemes have an arbitrarily large number of stages. This allows for optimization of the explicit part to improve its stability properties, increasing the allowable time step size. Importantly, only two implicit stages are required regardless of the total number of stages, meaning a larger global time step can be taken with significantly fewer implicit stages for a given simulation time. Numerical results demonstrate that AIMEX schemes achieve their designed order of accuracy for linear and non-linear problems involving mesh induced stiffness. Simulations of unsteady flow over an SD7003 airfoil using the compressible Navier-Stokes equations demonstrate that AIMEX schemes can significantly outperform classical explicit Runge-Kutta schemes by over a factor of 20, and conventional IMEX schemes by over a factor of two, with negligible impact on quantitative results. Finally, a demonstration case of Implicit Large Eddy Simulation (ILES) of flow over a stalled NACA0020 airfoil shows the utility of AIMEX schemes for wall-bounded turbulent flows.

Design of IMEXRK time integration schemes via Delaunay-based derivative-free optimization with nonconvex constraints and grid-based acceleration

This paper develops a powerful new variant, dubbed $$\varDelta $$ -DOGS( $$\varOmega _Z$$ ), of our Delaunay-based Derivative-free Optimization via Global Surrogates family of algorithms, and uses it to identify a new, low-storage, high-accuracy, Implicit/Explicit Runge–Kutta (IMEXRK) time integration scheme for the stiff ODEs arising in high performance computing applications, like the simulation of turbulence. The $$\varDelta $$ -DOGS( $$\varOmega _Z$$ ) algorithm, which we prove to be globally convergent under the appropriate assumptions, combines (a) the essential ideas of our $$\varDelta $$ -DOGS( $$\varOmega $$ ) algorithm, which is designed to efficiently optimize a nonconvex objective function f(x) within a nonconvex feasible domain $$\varOmega $$ defined by a number of constraint functions $$c_\kappa (x)$$ , with (b) our $$\varDelta $$ -DOGS(Z) algorithm, which reduces the number of function evaluations on the boundary of the search domain via the restriction that all function evaluations lie on a Cartesian grid, which is successively refined as the iterations proceed. The optimization of the parameters of low-storage IMEXRK schemes involves a complicated set of nonconvex constraints, which leads to a challenging disconnected feasible domain, and a highly nonconvex objective function; our simulations indicate significantly faster convergence using $$\varDelta $$ -DOGS( $$\varOmega _Z$$ ) as compared with the original $$\varDelta $$ -DOGS( $$\varOmega $$ ) optimization algorithm on the problem of tuning the parameters of such schemes. A low-storage third-order IMEXRK scheme with remarkably good stability and accuracy properties is ultimately identified using this approach, and is briefly tested on Burgers’ equation.

Modeling blood flow in viscoelastic vessels: the 1D augmented fluid–structure interaction system

Nowadays mathematical models and numerical simulations are widely used in the field of hemodynamics, representing a valuable resource to better understand physiological and pathological processes in different medical sectors. The theory behind blood flow modeling is closely related to the study of incompressible flow through compliant thin-walled tubes, starting from the incompressible Navier–Stokes equations. Furthermore, the mechanical interaction between blood flow and vessels wall must be properly described by the model. Recent works showed the benefits of characterizing the rheology of the vessel wall through a viscoelastic law. Taking into account the viscous contribution of the wall material and not simply the elastic one leads to a more realistic representation of the vessel behavior, which manifests not only an instantaneous elastic strain but also a viscous damping effect on pulse pressure waves, coupled to energy losses. In this context, the aim of this work is to propose an easily extensible one-dimensional mathematical model able to accurately capture fluid–structure interactions. The originality of the model lies in the introduction of a viscoelastic tube law in PDE form, valid for both arterial and venous networks, leading to an augmented fluid–structure interaction system. In contrast to well established mathematical models, the proposed one is natively hyperbolic. The model is solved with an efficient and robust second-order numerical scheme; the time integration is based on an Implicit–Explicit Runge–Kutta scheme conceived for applications to hyperbolic systems with stiff relaxation terms. The validation of the proposed model is performed on several different test cases. Results obtained in Riemann problems, adopting a simple elastic tube law for the characterization of the vessel wall, are compared with available exact solutions. To validate the contribution given by the viscoelastic term, the Method of Manufactured Solutions has been applied. Specific tests have been designed to verify the well-balancing with respect to fluid-at-rest condition and the accuracy-preserving property of the scheme. Finally, a specific test case with an inlet pulse pressure wave has been designed to assess the effects of viscoelasticity with respect to a simple elastic behavior of the vessel wall. The complete code, written in MATLAB (MathWorks Inc.) language, with the implemented test cases, is made available in Mendeley Data repository.

L2-stability analysis of IMEX-(σ,μ)DG schemes for linear advection-diffusion equations

A fully discrete L2-stability analysis is carried out for linear advection-diffusion problems in one space dimension, discretized in space by discontinuous Galerkin methods based on the (σ,μ)-family of diffusion schemes and implicit-explicit Runge-Kutta schemes in time. Advection terms are discretized explicitly while diffusion terms are solved implicitly. Conditions on the parameters σ and μ are derived which guarantee L2-stability for time steps Δt=O(d/a2), where a and d denote the advection and diffusion coefficient, respectively, i.e. the allowable time step size does not decrease under grid refinement. These conditions are fulfilled in particular by the BR2 scheme as well as the recent (14,94)-recovery scheme. However, the BR1 scheme and the Baumann-Oden method do not possess a similar stability property which is also confirmed by numerical experiments.