Implicit Runge-Kutta Schemes Research Articles

On convergence of implicit Runge-Kutta methods for the incompressible Navier-Stokes equations with unsteady inflow

This study investigates the convergence properties of implicit Runge-Kutta (IRK) methods when applied to the temporal solution of incompressible Navier-Stokes (N-S) equations with unsteady inflow. Owing to the differential-algebraic nature of spatially discretized N-S equations, conventional IRK methods may experience a significant order reduction while requiring exact satisfaction of the divergence-free constraint on the velocity field. Notably, the enhanced performance achieved through modified IRK techniques, such as projected Runge-Kutta methods and specialized Runge-Kutta methods, is confined to Runge-Kutta coefficients with specific attributes. In response to these limitations, this paper proposes a perturbed IRK scheme, modifying the intermediate stage equations of standard IRK methods by incorporating predefined perturbations, aiming to enhance convergence properties for the incompressible N-S equations accompanied by unsteady inflow. These perturbations within the scheme not only alleviate order reduction but also ensure exact enforcement of the divergence-free constraint. Moreover, the proposed scheme remains applicable in cases where the unsteady inflow is only available as discrete-time fields, rather than explicit functions of time. To demonstrate the efficiency of the proposed enhancement, an extensive analysis of the convergence properties for all considered IRK methods through a series of numerical experiments, is conducted.

The Matrix‐Free Macro‐Element Hybridized Discontinuous Galerkin Method for Steady and Unsteady Compressible Flows

ABSTRACTThe macro‐element variant of the hybridized discontinuous Galerkin (HDG) method combines advantages of continuous and discontinuous finite element discretization. In this paper, we investigate the performance of the macro‐element HDG method for the analysis of compressible flow problems at moderate Reynolds numbers. To efficiently handle the corresponding large systems of equations, we explore several strategies at the solver level. On the one hand, we utilize a second‐layer static condensation approach that reduces the size of the local system matrix in each macro‐element and hence the factorization time of the local solver. On the other hand, we employ a multi‐level preconditioner based on the FGMRES solver for the global system that integrates well within a matrix‐free implementation. In addition, we integrate a standard diagonally implicit Runge–Kutta scheme for time integration. We test the matrix‐free macro‐element HDG method for compressible flow benchmarks, including Couette flow, flow past a sphere, and the Taylor–Green vortex. Our results show that unlike standard HDG, the macro‐element HDG method can operate efficiently for moderate polynomial degrees, as the local computational load can be flexibly increased via mesh refinement within a macro‐element. Our results also show that due to the balance of local and global operations, the reduction in degrees of freedom, and the reduction of the global problem size and the number of iterations for its solution, the macro‐element HDG method can be a competitive option for the analysis of compressible flow problems.

Krylov-based adaptive-rank implicit time integrators for stiff problems with application to nonlinear Fokker-Planck kinetic models

We propose a high order adaptive-rank implicit integrators for stiff time-dependent PDEs, leveraging extended Krylov subspaces to efficiently and adaptively populate low-rank solution bases. This allows for the accurate representation of solutions with significantly reduced computational costs. We further introduce an efficient mechanism for residual evaluation and an adaptive rank-seeking strategy that optimizes low-rank settings based on a comparison between the residual size and the local truncation errors of the time-stepping discretization. We demonstrate our approach with the challenging Lenard-Bernstein Fokker-Planck (LBFP) nonlinear equation, which describes collisional processes in a fully ionized plasma. The preservation of the equilibrium state is achieved through the Chang-Cooper discretization, and strict conservation of mass, momentum and energy via a Locally Macroscopic Conservative (LoMaC) procedure. The development of implicit adaptive-rank integrators, demonstrated here up to third-order temporal accuracy via diagonally implicit Runge-Kutta schemes, showcases superior performance in terms of accuracy, computational efficiency, equilibrium preservation, and conservation of macroscopic moments. This study offers a starting point for developing scalable, efficient, and accurate methods for high-dimensional time-dependent problems.

Mass and energy conservative high-order diagonally implicit Runge–Kutta schemes for nonlinear Schrödinger equation

We present and analyze a series of structure-preserving diagonally implicit Runge–Kutta schemes for the nonlinear Schrödinger equation. These schemes possess not only high accuracy, high order convergence (up to fifth order) and efficiency due to the diagonally implicity but also mass and energy conservative properties. Theoretical analysis and numerical experiments are conducted to verify the accuracy, invariant conservative properties and longtime simulation stability.

Diagonally implicit Runge–Kutta schemes: Discrete energy-balance laws and compactness properties

Abstract We study diagonally implicit Runge–Kutta (DIRK) schemes when applied to abstract evolution problems that fit into the Gelfand-triple framework. We introduce novel stability notions that are well-suited to this setting and provide simple, necessary and sufficient, conditions to verify that a DIRK scheme is stable in our sense and in Bochner-type norms. We use several popular DIRK schemes in order to illustrate cases that satisfy the required structural stability properties and cases that do not. In addition, under some mild structural conditions on the problem we can guarantee compactness of families of discrete solutions with respect to time discretization.

On the accuracy, stability and computational efficiency of explicit last-stage diagonally implicit Runge–Kutta methods (ELDIRK) for the adaptive solution of phase-field problems

Runge–Kutta algorithms with adaptive step size control provide reliable tools for the solution of initial value problems. Most widespread within the vast family of Runge–Kutta methods are the diagonally implicit Runge–Kutta (DIRK) methods. In particular, the new low-order explicit last-stage diagonally implicit Runge–Kutta (ELDIRK) methods are investigated, combining implicit schemes with an additional explicit evaluation as an explicit last stage. This results in Butcher tableaus with two solutions of different convergence orders suitable for embedded methods, where the higher-order solution is achieved by additional explicit evaluations. Thus, the iterative solution of non-linear systems is omitted for the additional stage, presenting a major reduction in computational cost for the determination of a local error estimate. The key contribution is the application of the novel Butcher tableaus to phase-field problems, solved with the finite-element method, leading to substantial numerical investigations with an efficient approach for diagonally implicit Runge–Kutta schemes. The most important aspects are the extension towards fourth-order methods, the study of the convergence orders for the ELDIRK schemes, their respective stability regions and computational efficiency. A local error estimator is presented capturing the evolution of the phase-field problems, such that adaptive step size control for the new low-order embedded schemes based on an empirical approach for error estimation is achieved. A suitable parallel algorithm is presented with conclusive two-dimensional phase-field simulations based on a Kobayashi–Warren–Carter model including benchmarks for computational efficiency. The higher-order convergence suggested by the novel schemes is confirmed, and their effective results are demonstrated, resulting in a valuable semi-explicit addition to the family of Runge–Kutta time integration schemes.

Numerical investigations of new low‐order explicit last stage diagonal implicit Runge–Kutta schemes with the finite‐element method

AbstractInitial value problems can be solved efficiently by means of Runge–Kutta algorithms with adaptive step size control. Diagonally implicit Runge–Kutta (DIRK) methods are the most popular class among the diverse family of Runge–Kutta algorithms. In this paper, the novel class of low‐order explicit last‐stage diagonally implicit Runge–Kutta (ELDIRK) methods are explored, which combine implicit schemes with an additional explicit evaluation as an explicit last stage. ELDIRK Butcher tableaus are used to control embedded RK methods to obtain solutions of different orders. The lower‐order solution is obtained by classical implicit RK stages and the higher‐order solution is obtained by additional explicit evaluation. As a result, a significant reduction in computational cost is achieved by skipping the iterative solution of nonlinear systems for the additional step. The examination of the heat problem and the use of the innovative Butcher tableau in the finite‐element method are the main contributions of this work. Thus, it is possible to establish adaptive step size control for the new low‐order embedded methods based on an empirical method for error estimation. Two‐dimensional simulations are used to show an appropriate algorithm for the ELDIRK schemes. The new Runge–Kutta schemes' predictions of higher‐order convergence are confirmed, and their successful outcomes are illustrated.

Implicit Runge–Kutta Schemes for Optimal Control Problems with Evolution Equations

Abstract In this paper we discuss the use of implicit Runge–Kutta schemes for the time discretization of optimal control problems with evolution equations. The specialty of the considered discretizations is that the discretizations schemes for the state and adjoint state are chosen such that discretization and optimization commute. It is well known that for Runge–Kutta schemes with this property additional order conditions are necessary. We give sufficient conditions for which class of schemes these additional order condition are automatically fulfilled. The focus is especially on implicit Runge–Kutta schemes of Gauss, Radau IA, Radau IIA, Lobatto IIIA, Lobatto IIIB and Lobatto IIIC collocation type up to order six. Furthermore, we also use a SDIRK (singly diagonally implicit Runge–Kutta) method to demonstrate, that for general implicit Runge–Kutta methods the additional order conditions are not automatically fulfilled. Numerical examples illustrate the predicted convergence rates.

A high-order linearly implicit energy-preserving Partitioned Runge-Kutta scheme for a class of nonlinear dispersive equations

<abstract><p>In this paper, we design a novel class of arbitrarily high-order, linearly implicit and energy-preserving numerical schemes for solving the nonlinear dispersive equations. Based on the idea of the energy quadratization technique, the original system is firstly rewritten as an equivalent system with a quadratization energy. The prediction-correction strategy, together with the Partitioned Runge-Kutta method, is then employed to discretize the reformulated system in time. The resulting semi-discrete system is high-order, linearly implicit and can preserve the quadratic energy of the reformulated system exactly. Finally, we take the Camassa-Holm equation as a benchmark to show the efficiency and accuracy of the proposed schemes.</p></abstract>

Numerical evolution of the resistive relativistic magnetohydrodynamic equations: A minimally implicit Runge-Kutta scheme

We present the Minimally-Implicit Runge-Kutta (MIRK) methods for the numerical evolution of the resistive relativistic magnetohydrodynamic (RRMHD) equations, following the approach proposed by Komissarov (2007) of an augmented system of evolution equations to numerically deal with constraints. Previous approaches rely on Implicit-Explicit (IMEX) Runge-Kutta schemes; in general, compared to explicit schemes, IMEX methods need to apply the recovery (which can be very expensive computationally) of the primitive variables from the conserved ones in numerous additional times. Moreover, the use of an iterative process for the recovery could have potential convergence problems, increased by the additional number of required loops. In addition, the computational cost of the previous IMEX approach in comparison with the standard explicit methods is much higher. The MIRK methods are able to deal with stiff terms producing stable numerical evolutions, minimize the number of recoveries needed in comparison with IMEX methods, their computational cost is similar to the standard explicit methods and can actually be easily implemented in numerical codes which previously used explicit schemes. Two standard numerical tests are shown in the manuscript.

Spiral Bevel Gears: nonlinear dynamic model based on accurate static stiffness evaluation

In the present paper non-linear dynamics of a spiral bevel gear pair with backlash are investigated in order to clarify the internal excitations of major importance from the vibration point of view: manufacturing errors in the teeth profile, teeth spacing errors, and elastic deformation of the teeth. In some conditions, like in the case of backside contact, the destructive effect of internal excitations can be intensified leading to complex dynamics; for such reasons here backside contacts and reverse rotation are investigated in detail using a nonlinear time-varying model. The effect of damping is investigated as well. A one-DOF model is developed in order to study the dynamic behavior; the resulting a nonlinear differential equation with time-varying mesh stiffness is solved via numerical integration based on an adaptive step-size implicit Runge-Kutta scheme. The dynamic response of the system is analyzed through time histories, phase portraits, bifurcation diagrams, and Poincaré maps. Results show that for small backlash values, the possibility of backside contact increases. Meanwhile, by increasing the backlash value, the amplitude vibration of the gear rotation rises as well. By comparing the dynamic response of the system with different damping ratios, the results show that higher damping effectively reduces gear vibration resonance, although the probability of unsteady response still exists.

A fully coupled, fully implicit simulation method for unsteady flames using Jacobian approximation and clustering

Jacobian clustering is proposed to reduce the computational cost associated with matrix operations encountered in the Newton iteration in fully coupled, fully implicit schemes for unsteady reactive flow simulations with detailed chemistry. The iterative solver is based on the Lower-Upper Symmetric Gauss-Seidel (LUSGS) algorithm and sparse matrix technique. The evaluation and sparse Lower-Upper (LU) factorization of the diagonal block of the system Jacobian are performed within clusters rather than individual cells for Computational Fluid Dynamics (CFD) simulations. Cells that are close in the state space are clustered to provide the averaged states for calculating the Jacobians of chemical source terms and transport fluxes. The cells retrieve the factorized sparse matrices from the belonging clusters to perform the necessary iterations. For the purpose of clustering, the spatial dependency of transport Jacobian in the diagonal block is eliminated. To further reduce the computational cost, the sparsity of chemical Jacobian is augmented by removing the insignificant matrix elements. The method is tested in various one-dimensional hydrocarbon flames with both the second order Crank-Nicolson scheme and a third order implicit Runge-Kutta scheme. Various chemical mechanisms with 9 to 111 species are used to test the performance of the iterative solver. Fast convergence of Newton iteration is achieved, and the formal order of accuracy is demonstrated with Jacobian clustering. The overall costs of evaluation and factorization of the block diagonal Jacobian are negligible compared to the cost of calculating transport fluxes and chemical source terms. The averaged costs of Jacobian evaluation, LU factorization and Newton iteration, all increase only linearly with the number of chemical species. The fully coupled, fully implicit Crank-Nicolson scheme with Jacobian clustering shows 4 to 42 times speedup in computational time compared to the decoupled implicit scheme with Strang operator splitting. Jacobian clustering is promising to increase the computational efficiency of high order fully coupled, fully implicit schemes for unsteady reactive flow simulations with detailed chemistry.

Discrete analysis of Schwarz waveform relaxation for a diffusion reaction problem with discontinuous coefficients

In this paper, we investigate the effect of the space and time discretisation on the convergence properties of Schwarz Waveform Relaxation (SWR) algorithms. We consider a reaction-diffusion problem with discontinuous coefficients discretised on two non-overlapping domains with several numerical schemes (in space and time). A methodology to determine the rate of convergence of the classical SWR method with standard interface conditions (Dirichlet-Neumann or Robin-Robin) accounting for discretisation errors is presented. We discuss how such convergence rates differ from the ones derived at a continuous level (i.e. assuming an exact discrete representation of the continuous problem). In this work we consider a second-order finite difference scheme and a finite volume scheme based on quadratic spline reconstruction in space, combined with either a simple backward Euler scheme or a two-step “Padé” scheme (resembling a Diagonally Implicit Runge Kutta scheme) in time. We prove those combinations of space-time schemes to be unconditionally stable on bounded domains. We illustrate the relevance of our analysis with specifically designed numerical experiments.

Large-Eddy Simulation of Transonic Buffet Using Matrix-Free Discontinuous Galerkin Method

We present an implicit large-eddy simulation of transonic buffet over the OAT15A supercritical airfoil at Mach number 0.73, angle of attack 3.5 deg, and Reynolds number . The simulation is performed using a matrix-free discontinuous Galerkin (DG) method and a diagonally implicit Runge–Kutta scheme on graphics processor units. We propose a Jacobian-free Newton–Krylov method to solve nonlinear systems arising from the discretization of the Navier–Stokes equations. The method successfully predicts the buffet onset, the buffet frequency, and turbulence statistics owing to the high-order DG discretization and an efficient mesh refinement for the laminar and turbulent boundary layers. A number of physical phenomena present in the experiment are captured in our simulation, including periodical low-frequency oscillations of shock wave in the streamwise direction, strong shear layer detached from the shock wave due to shock-wave/boundary-layer interaction and small-scale structures broken down by the shear-layer instability in the transition region, and shock-induced flow separation. The pressure coefficient, the root mean square of the fluctuating pressure, and the streamwise range of the shock wave oscillation agree well with the experimental data. The results suggest that the proposed method can accurately predict the onset of turbulence and buffet phenomena at high Reynolds numbers without a subgrid scale model or a wall model.

Efficient high-order radial basis-function-based differential quadrature-finite volume method for incompressible flows on unstructured grids.

This paper presents an efficient high-order radial basis-function-based differential quadrature-finite volume method for incompressible flows on unstructured grids. In this method, a high-order polynomial based on the Taylor series expansion is applied within each control cell to approximate the solution. The derivatives in the Taylor series expansion are approximated by the mesh-free radial basis-function-based differential quadrature method. The recently proposed lattice Boltzmann flux solver is applied to simultaneously evaluate the inviscid and viscous fluxes at the cell interface by the local solution of the lattice Boltzmann equation. In the present high-order method, a premultiplied coefficient matrix appears in the time-dependent term, reflecting the implicit nature. The implicit time-marching techniques, i.e., the lower-upper symmetric Gauss-Seidel and the explicit first stage, singly diagonally implicit Runge-Kutta schemes, are incorporated to efficiently solve the resultant ordinary differential equations. Several numerical examples are tested to validate the accuracy, efficiency, and robustness of the present method on unstructured grids. Compared with the k-exact method, the present method enjoys higher accuracy and better computational efficiency.

High-order solution of Generalized Burgers–Fisher Equation using compact finite difference and DIRK methods

The main goal of this paper is to developed a high-order and accurate method for the solution of one-dimensional of generalized Burgers-Fisher with Numman boundary conditions. We combined between a fourth-order compact finite difference scheme for spatial part with diagonal implicit Runge Kutta scheme in temporal part. In addition, we discretized boundary points by using a compact finite difference scheme in terms of fourth order accuracy. This combine leads to ordinary differential equation which will take full advantage of method of line (MOL). Some numerical experiments presented to show that the combination give an accurate and reliable for solving the generalized Burgers-Fisher problems.

Hybrid High-Order Methods for the Acoustic Wave Equation in the Time Domain

We devise hybrid high-order (HHO) methods for the acoustic wave equation in the time domain. We first consider the second-order formulation in time. Using the Newmark scheme for the temporal discretization, we show that the resulting HHO-Newmark scheme is energy-conservative, and this scheme is also amenable to static condensation at each time step. We then consider the formulation of the acoustic wave equation as a first-order system together with singly-diagonally implicit and explicit Runge-Kutta (SDIRK and ERK) schemes. HHO-SDIRK schemes are amenable to static condensation at each time step. For HHO-ERK schemes, the use of the mixed-order formulation, where the polynomial degree of the cell unknowns is one order higher than that of the face unknowns, is key to benefit from the explicit structure of the scheme. Numerical results on test cases with analytical solutions show that the methods can deliver optimal convergence rates for smooth solutions of order $$\mathcal{O}(h^{k+1})$$ in the $$H^1$$ -norm and of order $$\mathcal{O}(h^{k+2})$$ in the $$L^2$$ -norm. Moreover, test cases on wave propagation in heterogeneous media indicate the benefits of using high-order methods.

Spiral Bevel Gears Nonlinear Vibration Having Radial and Axial Misalignments Effects

In gear transmissions, vibration causes noise and malfunction. In actual applications, misalignments contribute to intensifying the destructive effect of vibrations. In this paper, the nonlinear dynamics of a spiral bevel gear pair, with small helix angle, considering different misalignments, are deeply investigated. Axial misalignment, radial misalignment, and the combination of these two types are considered in this study. The governing equation is numerically solved through an implicit Runge–Kutta scheme. Since the main goal of this study is the analysis of the dynamic scenario, the mesh stiffness of the gear pair is obtained from the literature. The dynamical system is nonlinear and time-varying; it is analyzed through time responses, phase portraits, Poincaré maps, and bifurcation diagrams. Results show that, among the considered three cases with different types of misalignments, the spiral bevel gear with axial misalignment is the worst destructive case; aperiodic, subharmonic, and multiperiod responses are observable for this case. It is interesting that the chaotic responses for the case, having both types of misalignments, are less likely for the case with axial misalignment, only.

Low Mach preconditioned density-based methods with implicit Runge–Kutta schemes in physical-time

Low Mach preconditioned density-based methods with implicit Runge–Kutta schemes in physical-time

A symmetry-preserving second-order time-accurate PISO-based method

A new conservative symmetry-preserving second-order time-accurate PISO-based pressure-velocity coupling for solving the incompressible Navier-Stokes equations on unstructured collocated grids is presented in this paper. This new method for implicit time stepping is an extension of the conservative symmetry-preserving incremental-pressure projection method for explicit time stepping and unstructured collocated meshes of Trias et al. [35]. In order to assess and compare both methods, we have implemented them within one unified solver in the open source code OpenFOAM where we use a Butcher array to prescribe the Runge-Kutta method. Thus, by changing the entries of the Butcher array, explicit and diagonally implicit Runge-Kutta schemes can be combined into one solver. We assess the energy conservation properties of the implemented discretisation methods and the temporal consistency of the selected Runge-Kutta schemes using Taylor-Green vortex and lid-driven cavity flow test cases. Finally, we use a more complex turbulent channel flow test case in order to further assess the performance of the presented new conservative symmetry-preserving incremental-pressure PISO-based method.Although both implemented methods are based on a symmetry-preserving discretisation, we show they still produce a small amount of numerical dissipation when the total pressure is directly solved from a Poisson equation. When an incremental-pressure approach is used, where a pressure correction is solved from a Poisson equation, both methods are effectively fully-conservative. For high-fidelity simulations of incompressible turbulent flows, it is highly desirable to use fully-conservative methods. For such simulations, the presented numerical methods are therefore expected to have large added value, since they pave the way for the execution of truly energy-conservative high-fidelity simulations in complex geometries. Furthermore, both methods are implemented in OpenFOAM, which is widely used within the CFD community, so that a large part of this community can benefit from the developed and implemented numerical methods.