Runge-Kutta Time Discretization Research Articles

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.

Double-pole unsplit complex-frequency-shifted multiaxial perfectly matched layer combined with strong-stability-preserved Runge-Kutta time discretization for seismic wave equation based on the discontinuous Galerkin method

The classical complex-frequency-shifted perfectly matched layer (CFS-PML) technique has attracted widespread attention for seismic wave simulations. However, few studies have addressed the double-pole variant of the CFS-PML scheme. The double-pole CFS-PML has a stronger capacity to absorb near-grazing incident waves and evanescent waves than the classical CFS-PML. Using the discontinuous Galerkin (DG) method, we derive a double-pole unsplit auxiliary ordinary differential equation CFS-multiaxial PML (AODE CFS-MPML) formulation, which combines a fourth-order strong-stability-preserved Runge-Kutta time discretization for wavefield simulation on an unstructured grid. The double-pole unsplit CFS-MPML formulations are obtained by introducing auxiliary memory variables and AODEs. The original stress-velocity equations and the double-pole unsplit AODE CFS-MPML equations are all first-order hyperbolic systems and suitable for the DG method. The attenuative variables are added directly to the original seismic wave equations without changing their formats. In contrast to the split perfectly matched layer (PML), we avoid reformulating PML equations in the nonattenuative modeling region. The original seismic wave equation is solved in the nonattenuative modeling domain, whereas the double-pole unsplit AODE CFS-MPML equation is implemented in the PML absorbing region. Three numerical examples validate the performance of the double-pole unsplit AODE CFS-MPML technique. The isotropic and anisotropic experiments demonstrate that our developed double-pole unsplit AODE CFS-MPML is more stable and obtains more accurate solutions than the classical CFS-PML. The second example indicates the flexibility of the combination of the DG method with the double-pole CFS-MPML on undulating topography. The final example displays the applicability and effectiveness of our method in a 3D situation.

Uniform stability for local discontinuous Galerkin methods with implicit-explicit Runge-Kutta time discretizations for linear convection-diffusion equation

In this paper, we consider the linear convection-diffusion equation in one dimension with periodic boundary conditions, and analyze the stability of fully discrete methods that are defined with local discontinuous Galerkin (LDG) methods in space and several implicit-explicit (IMEX) Runge-Kutta methods in time. By using the forward temporal differences and backward temporal differences, respectively, we establish two general frameworks of the energy-method based stability analysis. From here, the fully discrete schemes being considered are shown to have monotonicity stability, i.e. the L 2 L^2 norm of the numerical solution does not increase in time, under the time step condition τ ≤ F ( h / c , d / c 2 ) \tau \le \mathcal {F}(h/c, d/c^2) , with the convection coefficient c c , the diffusion coefficient d d , and the mesh size h h . The function F \mathcal {F} depends on the specific IMEX temporal method, the polynomial degree k k of the discrete space, and the mesh regularity parameter. Moreover, the time step condition becomes τ ≲ h / c \tau \lesssim h/c in the convection-dominated regime and it becomes τ ≲ d / c 2 \tau \lesssim d/c^2 in the diffusion-dominated regime. The result is improved for a first order IMEX-LDG method. To complement the theoretical analysis, numerical experiments are further carried out, leading to slightly stricter time step conditions that can be used by practitioners. Uniform stability with respect to the strength of the convection and diffusion effects can especially be relevant to guide the choice of time step sizes in practice, e.g. when the convection-diffusion equations are convection-dominated in some sub-regions.

On the conservative phase-field method with the N-component incompressible flows

This paper presents a conservative Allen–Cahn model coupled with the incompressible Navier–Stokes equation for tracking the interface with the N-component immiscible fluids system. The proposed conservative phase-field model can track the interface with large deformation in divergence-free velocity fields. The erroneous estimation of the normal vector is a challenging numerical issue for the interface capturing due to the appearance of spurious oscillations. The improved phase-field-based method combines the nonlinear preprocessing operation guided by the level-set method with local artificial viscosity stabilization to improve the computation of the discrete normal vector. The interfaces between different immiscible components are replaced by the transition region with finite thickness in the continuous phase field. The surface tension effects are represented with the continuous surface tension force in the system, which is not limited by the number of components. The third-order Runge-Kutta time discretization and second-order spatial discretization are applied for the multi-component system. To eliminate the spurious oscillations caused by discontinuous and steep gradient for capturing the shocks and sharp interfaces, we apply the third-order weighted essentially non-oscillatory method for the advection term. Several quantitative results of numerical tests, such as error estimation with the proposed method, comparative tests with different methods, and convergence rate for classical benchmark test, have been performed to illustrate that our method works well for the interface tracing issue with high numerical accuracy. In addition, various representative qualitative tests have been presented to demonstrate the applicability of our method.

Fully scalable solver for frequency-domain visco-elastic wave equations in 3D heterogeneous media: A controllability approach

We develop a controllability strategy for the computation of frequency-domain solutions of the 3D visco-elastic wave equation, in the perspective of seismic imaging applications. We generalize the controllability results for such equations beyond the sound-soft scattering (obstacle) problem. We detail the conjugate gradient implementation and show how an inner elliptic problem needs to be solved to compute the Riesz representative of the gradient at each iteration. We select a spectral-element spatial discretization and a fourth-order Runge-Kutta time discretization. We implement the controllability method in the framework of the SEM46 full waveform modeling and inversion software, to inherit for its excellent scalability which relies on an efficient domain decomposition algorithm. We perform a series of numerical experiments to validate the approach and illustrate its scalability up to more than fifteen hundred cores. In this case, with an elapsed time of less than 50 minutes, we solve a problem on a cubic domain containing up to 160 wavelengths in each direction, involving more than 1.7 billion unknowns.

Energy Plus Maximum Bound Preserving Runge–Kutta Methods for the Allen–Cahn Equation

It is difficult to design high order numerical schemes which could preserve both the maximum bound property (MBP) and energy dissipation law for certain phase field equations. Strong stability preserving (SSP) Runge–Kutta methods have been developed for numerical solution of hyperbolic partial differential equations in the past few decades, where strong stability means the non-increasing of the maximum bound of the underlying solutions. However, existing framework of SSP RK methods can not handle nonlinear stabilities like energy dissipation law. The aim of this work is to extend this SSP theory to deal with the nonlinear phase field equation of the Allen–Cahn type which typically satisfies both maximum bound preserving (MBP) and energy dissipation law. More precisely, for Runge–Kutta time discretizations, we first derive a general necessary and sufficient condition under which MBP is satisfied; and we further provide a necessary condition under which the MBP scheme satisfies energy dissipation.

An implicit monolithic AFC stabilization method for the CG finite element discretization of the fully-ionized ideal multifluid electromagnetic plasma system

This study considers an implicit finite element formulation for an ideal fully-ionized multifluid electromagnetic plasma system. The formulation is based on fully-implicit Runge-Kutta time discretizations and a monolithic discrete algebraic flux corrected (AFC) continuous Galerkin (CG) spatial discretization of the coupled system. The AFC approach adds scalar artificial diffusion to the high-order, semi-discrete Galerkin method and uses mass lumping in the time derivative term. The result is a low-order method that attempts to enforce local-extremum-diminishing properties for the hyperbolic system. An element-based iterative limiter is applied to reduce the amount of artificial diffusion that is used in regions where the solution is smooth and the additional stabilization is not required. Two models are considered for the electromagnetics portion of the system: an electrostatic model, and a full Maxwell system with a parabolic divergence cleaning approach that enforces the required involutions on the electric and magnetic fields. Results are presented that demonstrate the accuracy and robustness of the formulation for smooth and discontinuous solutions to challenging plasma physics problems. This includes a demonstration that the solution of the full multifluid system yields the expected behavior in the ideal shock-MHD limit.

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.

Energy conserving and well-balanced discontinuous Galerkin methods for the Euler–Poisson equations in spherical symmetry

ABSTRACT This paper presents high-order Runge–Kutta (RK) discontinuous Galerkin methods for the Euler–Poisson equations in spherical symmetry. The scheme can preserve a general polytropic equilibrium state and achieve total energy conservation up to machine precision with carefully designed spatial and temporal discretizations. To achieve the well-balanced property, the numerical solutions are decomposed into equilibrium and fluctuation components that are treated differently in the source term approximation. One non-trivial challenge encountered in the procedure is the complexity of the equilibrium state, which is governed by the Lane–Emden equation. For total energy conservation, we present second- and third-order RK time discretization, where different source term approximations are introduced in each stage of the RK method to ensure the conservation of total energy. A carefully designed slope limiter for spherical symmetry is also introduced to eliminate oscillations near discontinuities while maintaining the well-balanced and total-energy-conserving properties. Extensive numerical examples – including a toy model of stellar core collapse with a phenomenological equation of state that results in core bounce and shock formation – are provided to demonstrate the desired properties of the proposed methods, including the well-balanced property, high-order accuracy, shock-capturing capability, and total energy conservation.

High order well-balanced asymptotic preserving finite difference WENO schemes for the shallow water equations in all Froude numbers

In this paper, high order semi-implicit well-balanced and asymptotic preserving finite difference WENO schemes are proposed for the shallow water equations with a non-flat bottom topography. We consider the Froude number ranging from O(1) to 0, which in the zero Froude limit becomes the “lake equations” for balanced flow without gravity waves. We apply a well-balanced finite difference WENO reconstruction, coupled with a stiffly accurate implicit-explicit (IMEX) Runge-Kutta time discretization. The resulting semi-implicit scheme can be shown to be well-balanced, asymptotic preserving (AP) and asymptotically accurate (AA) at the same time. Both one- and two-dimensional numerical results are provided to demonstrate the high order accuracy, AP property and good performance of the proposed methods in capturing small perturbations of steady state solutions.

A High–Order WENO Scheme Based on Different Numerical Fluxes for the Savage–Hutter Equations

The study of rapid free surface granular avalanche flows has attracted much attention in recent years, which is widely modeled using the Savage–Hutter equations. The model is closely related to shallow water equations. We employ a high-order shock-capturing numerical model based on the weighted essential non-oscillatory (WENO) reconstruction method for solving Savage–Hutter equations. Three numerical fluxes, i.e., Lax–Friedrichs (LF), Harten–Lax–van Leer (HLL), and HLL contact (HLLC) numerical fluxes, are considered with the WENO finite volume method and TVD Runge–Kutta time discretization for the Savage–Hutter equations. Numerical examples in 1D and 2D space are presented to compare the resolution of shock waves and free surface capture. The numerical results show that the method proposed provides excellent performance with high accuracy and robustness.

An Efficient Discontinuous Galerkin Method Using a Tetrahedral Mesh for 3D Seismic Wave Modeling

ABSTRACTThe discontinuous Galerkin (DG) method is a numerical algorithm that is widely used in various fields. It has high accuracy and low numerical dispersion with advantages of easy handling boundary conditions and good parallel performance. In this study, we develop an efficient parallel weighted Runge–Kutta discontinuous Galerkin (WRKDG) method on unstructured meshes for solving 3D seismic wave equations. The DG method we use is based on the first-order formulation of a hyperbolic system with an explicit weighted Runge–Kutta time discretization. We describe the numerical scheme and parallel implementation in detail, and analyze the stability condition and numerical dispersion and dissipation. We carry out a convergence test on unstructured meshes and investigate the parallel efficiency of the implementation of the WRKDG method. The speedup curve shows that the method has good parallel performance. Finally, we present several numerical simulation examples, including acoustic and elastic wave propagations in isotropic and anisotropic media. Numerical results further verify the effectiveness of the WRKDG method in solving 3D wave propagation problems.

Bound-preserving Flux Limiting for High-Order Explicit Runge–Kutta Time Discretizations of Hyperbolic Conservation Laws

We introduce a general framework for enforcing local or global maximum principles in high-order space-time discretizations of a scalar hyperbolic conservation law. We begin with sufficient conditions for a space discretization to be bound preserving (BP) and satisfy a semi-discrete maximum principle. Next, we propose a global monolithic convex (GMC) flux limiter which has the structure of a flux-corrected transport (FCT) algorithm but is applicable to spatial semi-discretizations and ensures the BP property of the fully discrete scheme for strong stability preserving (SSP) Runge–Kutta time discretizations. To circumvent the order barrier for SSP time integrators, we constrain the intermediate stages and/or the final stage of a general high-order RK method using GMC-type limiters. In this work, our theoretical and numerical studies are restricted to explicit schemes which are provably BP for sufficiently small time steps. The new GMC limiting framework offers the possibility of relaxing the bounds of inequality constraints to achieve higher accuracy at the cost of more stringent time step restrictions. The ability of the presented limiters to recognize undershoots/overshoots, as well as smooth solutions, is verified numerically for three representative RK methods combined with weighted essentially nonoscillatory (WENO) finite volume space discretizations of linear and nonlinear test problems in 1D. In this context, we enforce global bounds and prove preservation of accuracy for the linear advection equation.

High order semi-implicit weighted compact nonlinear scheme for the full compressible Euler system at all Mach numbers

A high order semi-implicit weighted compact nonlinear scheme (WCNS) is presented to solve the full compressible Euler equations at all Mach numbers. To avoid the stringent acoustic CFL restriction for an explicit time discretization method, the pressure splitting methodology is introduced to split the Euler equations into nonstiff and stiff terms. The nonstiff and stiff terms are treated explicitly and implicitly, respectively, based on the high order IMEX Runge-Kutta time discretization method that is asymptotic preserving and asymptotically accurate in the zero Mach number limit. A fifth-order WCNS and fourth-order centered finite difference schemes with zero numerical viscosity are used for the spatial discretization. Numerical tests in one and two space dimensions are displayed to demonstrate the performance of the high order semi-implicit WCNS in compressible and incompressible regimes. Comparisons with the third-order explicit weighted essentially non-oscillatory (WENO) scheme are also made in order to better assess the proposed scheme.

A novel hybrid method based on discontinuous Galerkin method and staggered‐grid method for scalar wavefield modelling with rough topography

ABSTRACTIn recent years, the discontinuous Galerkin method and the finite‐difference method have been widely applied to simulate seismic wave propagation. However, few studies have focused on the combination of the finite‐difference method with the discontinuous Galerkin method. We develop a new hybrid algorithm based on domain decomposition that combines the high efficiency of the finite‐difference method with the flexibility of the discontinuous Galerkin method. A computational domain is decomposed into two types of subregions where the finite‐difference method is applied in the main part of the model, whereas the discontinuous Galerkin method is employed to handle the complex structure with a free surface or caves. Total variation diminishing Runge–Kutta time discretization is used to enhance the stability of implementation. The approach inherits the advantages of the discontinuous Galerkin method and the finite‐difference method and does not need a transition zone, which reduces the computational cost caused by the additional conversion. In addition, some mesh patterns are provided to make the discontinuous Galerkin method domain discretized as needed without changing the grid size of the finite‐difference method region, which makes our hybrid algorithm more flexible. Numerical experiments prove that our hybrid scheme is capable of dealing with complex structures and maintains moderate computational efficiency.

High-Order Godunov-Type Scheme for Conservation Laws with Discontinuous Flux Function in Space

The numerical method for the conservation laws with discontinuous flux in space is considered. The difficulty in designing a high-order accurate and maintaining an efficient well-balanced scheme for this system lies in the fact that the flux is discontinuous across the stationary discontinuity, which results in the jump of the unknown function. In order to overcome this difficulty, the Godunov-type numerical flux is written in the positive and negative fluxes for the first time, which is embedded between two discontinuous flux functions. With the Godunov flux-splitting, the high-order accurate scheme is proposed through the WENO reconstruction of flux and the third-order TVD Runge–Kutta time discretization. Some tests are presented to demonstrate that our new high-order scheme is capable of exactly preserving steady-state solutions.

Internal resonance and nonlinear dynamics of a dielectric elastomer circular membrane

Dielectric elastomers (DE) are widely used in several applications such as artificial muscles, sensors or energy harvesters. Their hyperelastic nature complicates their design because of the obtained rich dynamics. The electrical field applied in the material, couples the dynamics of these devices to the applied voltage through compliant electrodes. In this work, a circular DE membrane under mechanical and electrical excitations is analyzed in order to seek internal resonances capable of increasing its performance for future applications. The hyperelastic properties of the material and the Maxwell stresses are taken into account to derive the nonlinear governing equation. Analytical and numerical approaches are used to obtain the dynamic solution of the device under different prestretch and electrical conditions. The effect of the prestretch on the first transverse and radial natural frequencies is first calculated analytically and validated by a finite element software. One-to-one, two-to-one and three-to-one internal resonances are possible as the prestretch applied to the membrane is varied. Using a Galerkin procedure and a multiple scale technique, the frequency response curves are obtained for different applied voltages. The results are validated using a Runge–Kutta time discretization technique. Finally, the effect of the applied voltage is analyzed for the frequency response curves of the radial and transverse displacements.

A cell-centered implicit-explicit Lagrangian scheme for a unified model of nonlinear continuum mechanics on unstructured meshes

A cell-centered implicit-explicit updated Lagrangian finite volume scheme on unstructured grids is proposed for a unified first-order hyperbolic formulation of continuum fluid and solid mechanics, namely the Godunov-Peshkov-Romenski (GPR) model. The scheme provably respects the stiff relaxation limits of the continuous model at the fully discrete level, thus it is asymptotic preserving. Furthermore, the GCL is satisfied by a compatible discretization that makes use of a nodal solver to compute vertex-based fluxes that are used both for the motion of the computational mesh as well as for the time evolution of the governing PDEs. Second order of accuracy in space is achieved using a TVD piecewise linear reconstruction, while an implicit-explicit (IMEX) Runge-Kutta time discretization allows the scheme to obtain higher accuracy also in time. Particular care is devoted to the design of a stiff ODE solver, based on approximate analytical solutions of the governing equations, that plays a crucial role when the visco-plastic limit of the model is approached. We demonstrate the accuracy and robustness of the scheme on a wide spectrum of material responses covered by the unified continuum model that includes inviscid hydrodynamics, viscous heat conducting fluids, elastic and elasto-plastic solids in multidimensional settings.

A Low Mach Number IMEX Flux Splitting for the Level Set Ghost Fluid Method

Considering droplet phenomena at low Mach numbers, large differences in the magnitude of the occurring characteristic waves are presented. As acoustic phenomena often play a minor role in such applications, classical explicit schemes which resolve these waves suffer from a very restrictive timestep restriction. In this work, a novel scheme based on a specific level set ghost fluid method and an implicit-explicit (IMEX) flux splitting is proposed to overcome this timestep restriction. A fully implicit narrow band around the sharp phase interface is combined with a splitting of the convective and acoustic phenomena away from the interface. In this part of the domain, the IMEX Runge-Kutta time discretization and the high order discontinuous Galerkin spectral element method are applied to achieve high accuracies in the bulk phases. It is shown that for low Mach numbers a significant gain in computational time can be achieved compared to a fully explicit method. Applications to typical droplet dynamic phenomena validate the proposed method and illustrate its capabilities.

Positivity preserving and entropy consistent approximate Riemann solvers dedicated to the high-order MOOD-based Finite Volume discretization of Lagrangian and Eulerian gas dynamics

In this paper we propose to revisit the notion of simple Riemann solver both in Lagrangian and Eulerian coordinates following the seminal work of Gallice ”Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates” in Numer. Math., 94, 2003. We provide in this work the relation between the Eulerian and Lagrangian forms of systems of conservation laws in 1D. Then an approximate (simple) Lagrangian Riemann solver for gas dynamics is derived based on the notions of positivity preservation and entropy control. Its Eulerian counterpart is further deduced. Next we build the associated 1D first-order accurate cell-centered Lagrangian Godunov-type Finite Volume scheme and show numerically its behaviors on classical test cases. Then using the Lagrangian–Eulerian relationships, we derive and test the Eulerian Godunov-type Finite Volume scheme, which inherits by construction the properties of the Lagrangian solver in terms of positivity preservation and well-defined CFL condition. At last we extend this Eulerian scheme to arbitrary orders of accuracy using a Runge–Kutta time discretization, polynomial reconstruction and an a posteriori MOOD limiting strategy. Numerical tests are carried out to assess the robustness, accuracy, and essentially non-oscillatory properties of the numerical methods.