Runge-Kutta Time Stepping Research Articles

Finite Elements with Switch Detection for direct optimal control of nonsmooth systems

This paper introduces Finite Elements with Switch Detection (FESD), a numerical discretization method for nonsmooth differential equations. We consider the Filippov convexification of these systems and a transformation into dynamic complementarity systems introduced by Stewart (Numer Math 58(1):299–328, 1990). FESD is based on solving nonlinear complementarity problems and can automatically detect nonsmooth events in time. If standard time-stepping Runge–Kutta (RK) methods are naively applied to a nonsmooth ODE, the accuracy is at best of order one. In FESD, we let the integrator step size be a degree of freedom. Additional complementarity conditions, which we call cross complementarities, enable exact switch detection, hence FESD can recover the high order accuracy that the RK methods enjoy for smooth ODE. Additional conditions called step equilibration allow the step size to change only when switches occur and thus avoid spurious degrees of freedom. Convergence results for the FESD method are derived, local uniqueness of the solution and convergence of numerical sensitivities are proven. The efficacy of FESD is demonstrated in several simulation and optimal control examples. In an optimal control problem benchmark with FESD, we achieve up to five orders of magnitude more accurate solutions than a standard time-stepping approach for the same computational time.

On the Convergence of Monolithic Multigrid for Implicit Runge–Kutta Time Stepping of Finite Element Problems

On the Convergence of Monolithic Multigrid for Implicit Runge–Kutta Time Stepping of Finite Element Problems

A Wetting and Drying Approach for a Mode-Nonsplit Discontinuous Galerkin Hydrodynamic Model with Application to Laizhou Bay

A wetting and drying treatment for a three-dimensional discontinuous Galerkin hydrodynamic model without mode splitting (external and internal modes) was developed. In this approach, computing elements are classified into wet, dry, and semidry elements, which are treated differently. In a Runge–Kutta time step, the reconstruction of the semidry elements and the combined utilization of two- and three-dimensional limiters help the model maintain stability. Numerical results show that the wetting and drying method can achieve a well-balanced property under the condition of still-water equilibrium and can reasonably describe the variation process of wetting and drying regions during a long wave run-up on a uniform slope and a tidal cycle in a basin with a variable slope. Analysis of the role of the limiters in the model indicated that the robustness of the three-dimensional hydrodynamic model can be effectively maintained when the two- and three-dimensional limiters are jointly applied for wetting and drying process simulation. A three-dimensional discontinuous Galerkin hydrodynamic model was applied with the presented wetting and drying method to simulate the tidal current evolution of a spring tidal cycle in southwestern Laizhou Bay in the Bohai Sea, in November 2003, and the simulated results of the water surface elevation and vertical layered current velocities agreed well with the measured data.

Convergence of a second-order scheme for non-local conservation laws

In this article, we present the convergence analysis of a second-order numerical scheme for traffic flow models that incorporate non-local conservation laws. We combine a MUSCL-type spatial reconstruction with strong stability preserving Runge-Kutta time-stepping to devise a fully discrete second-order scheme. The resulting scheme is shown to converge to a weak solution by establishing the maximum principle, bounded variation estimates and L1Lipschitz continuity in time. Further, using a space-step dependent slope limiter, we prove its convergence to the entropy solution. We also propose a MUSCL-Hancock type second-order scheme which requires only one intermediate stage unlike the Runge-Kutta schemes and is easier to implement. The performance of the proposed second-order schemes in comparison to a first-order scheme is demonstrated through several numerical experiments.

Solitary-Wave Solutions of Benjamin–Ono and Other Systems for Internal Waves: II. Dynamics

Considered here are two systems of equations modeling the two-way propagation of long-crested, long-wavelength internal waves along the interface of a two-layer system of fluids in the Benjamin–Ono and the Intermediate Long-Wave regime, respectively. These systems were previously shown to have solitary-wave solutions, decaying to zero algebraically for the Benjamin–Ono system, and exponentially in the Intermediate Long-Wave regime. Several methods to approximate solitary-wave profiles were introduced and analyzed by the authors in Part I of this project. A natural continuation of this previous work, pursued here, is to study the dynamics of the solitary-wave solutions of these systems. This will be done by computational means using a discretization of the periodic initial-value problem. The numerical method used here is a Fourier spectral method for the spatial approximation coupled with a fourth-order, explicit Runge–Kutta time stepping. The resulting, fully discrete scheme is used to study computationally the stability of the solitary waves under small and large perturbations, the collisions of solitary waves, the resolution of initial data into trains of solitary waves, and the formation of dispersive shock waves. Comparisons with related unidirectional models are also undertaken.

Verification and accuracy check of simulations with PoPe and iPoPe

The theoretical background of the PoPe and iPoPe verification scheme is presented. Verification is performed using the simulation output of production runs. The computing overhead is estimated to be at most 10%. PoPe or iPoPe calculations can be done offline provided the necessary data is stored, for example additional time slices, or online where iPoPe is more effective. The computing overhead is mostly that of storing the necessary data. The numerical error is determined and split into a part proportional to the operators, which are combined to form the equations to be solved, thus modifying their control parameters, completed by a residual error orthogonal to these operators. The accuracy of the numerical solution is determined by this modification of the control parameters. The PoPe and iPoPe methods are illustrated in this paper with simulations of a simple mechanical system with chaotic trajectories evolving into a strange attractor with sensitivity to initial conditions. We show that the accuracy depends on the particular simulation both because the properties of the numerical solution depend on the values of the control parameter, and because the target accuracy will depend on the problem that is addressed. One shows that for a case close to bifurcations between different states, the accuracy is determined by the level of detail of the bifurcation phenomena one aims at describing. A unique verification index, the PoPe index, is proposed to characterise the accuracy, and consequently the verification, of each production run. The PoPe output allows one to step beyond verification and analyse for example the numerical scheme efficiency. For the chosen example at fixed PoPe index, therefore at fixed numerical error, one finds that the higher order integration scheme, comparing order 4 to order 2 Runge-Kutta time stepping, reduces the computation cost by a factor 4.

Discontinuous curvilinear collocated grid combined with nonuniform time step Runge-Kutta scheme for poroelastic finite-difference modeling

For poroelastic media, the existence of a slow P-wave mode, next to the standard fast P and S waves, hinders efficient numerical implementations to propagate poroelastic waves through arbitrary seismic models. The slow P-wave speed can be an order of magnitude smaller than the fast P-wave speed. Hence, a stable and accurate simulation that can capture the slow P wave requires a fine grid and a small time step, which increases the overall computation cost greatly. To decrease the computation cost, we propose a poroelastic finite-difference simulation method that combines a discontinuous curvilinear collocated-grid method with a nonuniform time step Runge-Kutta (NUTS-RK) scheme. The fine grid and small time step are only used for areas near interfaces, where the contribution of the slow P wave is nonnegligible. The NUTS-RK scheme is derived from a Taylor expansion and it can circumvent the need for interpolations or extrapolations otherwise required by communications between different time levels. The accuracy and efficiency of the proposed method are verified by numerical tests. Compared with the curvilinear collocated-grid finite-difference method that uses a globally uniform space grid as well as a uniform time step, the proposed method requires fewer computing resources and can reduce the computing time greatly.

High-order bicompact schemes for the quasilinear multidimensional diffusion equation

High-order bicompact schemes are designed for the quasilinear multidimensional diffusion equation. They are constructed using the method of lines, the finite volume method, the bicubic (or tricubic) Hermite interpolation. Implicit-explicit and diagonally-implicit Runge-Kutta methods are applied to time integration. The resulting schemes are stable for any ratios between grid steps, are conservative, have approximation of the fourth order in space and the third order in time. To implement the implicit-explicit schemes, an iteration method is proposed based on the approximate factorization of their multidimensional difference operators. This method is modified to implement the schemes with diagonally-implicit Runge-Kutta time stepping. High-order grid convergence and implementation efficiency of the new bicompact schemes are demonstrated on numerical examples.

Irksome: Automating Runge–Kutta Time-stepping for Finite Element Methods

While implicit Runge–Kutta (RK) methods possess high order accuracy and important stability properties, implementation difficulties and the high expense of solving the coupled algebraic system at each time step are frequently cited as impediments. We present Irksome , a high-level library for manipulating UFL (Unified Form Language) expressions of semidiscrete variational forms to obtain UFL expressions for the coupled Runge–Kutta stage equations at each time step. Irksome works with the Firedrake package to enable the efficient solution of the resulting coupled algebraic systems. Numerical examples confirm the efficacy of the software and our solver techniques for various problems.

An alternative SPH formulation: ADER-WENO-SPH

We present a new class of fully-discrete one-step SPH schemes based on a mesh-free ADER (Arbitrary DERivatives in space and time) reconstruction on moving particles in multiple space dimensions. In particular, the new SPH scheme computes mesh-free and local high order accurate polynomials in space and time to evaluate numerical fluxes at the midpoint between two interaction particles with a proper Riemann solver within the general SPH framework of Vila (1999) for nonlinear systems of hyperbolic conservation laws. The new scheme has been carefully tested against reference solutions for both the compressible Euler and the magneto-hydrodynamics (MHD) equations. The capability of the proposed scheme to accurately capture shocks and rarefaction waves for 1D and 2D problems with minimal amount of diffusion has been demonstrated. Via numerical evidence it has been shown that the new fully-discrete one-step ADER-WENO-SPH method is computationally more efficient than WENO-SPH schemes based on classical Runge–Kutta time-stepping. This is mainly due to the fact that with ADER timestepping the expensive stencil and neighbor search needs to be done only once per time step, while with Runge–Kutta time integrators the neighbor and stencil search is needed in each Runge–Kutta stage again.

The validation of a semi-recursive vehicle dynamics model for a real-time simulation

Semi-recursive formulations and their various versions have made it possible to describe complex nonlinear systems such as vehicles precisely while still solving the relevant equations of motion in real time. An optimal combination of an efficient multibody formulation and a fast numerical time integration scheme are needed to accurately simulate complex systems in real time. This paper introduces a double-step semi-recursive multibody formulation and analyzes its performance with high-order numerical time-integration algorithms for real-time simulation. The Runge-Kutta, Gill, Runge-Kutta-Fehlberg, Adams-Bashforth-Moulton, and adaptive time step Runge-Kutta numerical methods are explained and compared. Results are verified against a commercial multibody software solution. A 15-degree-of-freedom sedan vehicle model serves as a benchmark to verify the theoretical results. The results highlight the differences between the numerical algorithms and suggest appropriate approaches for a nonlinear vehicle dynamics model, particularly for cases where simulation times are long.

Optimal, globally constraint-preserving, DG(TD)2 schemes for computational electrodynamics based on two-derivative Runge-Kutta timestepping and multidimensional generalized Riemann problem solvers – A von Neumann stability analysis

Discontinuous Galerkin (DG) methods have become mainstays in the accurate solution of hyperbolic systems, which suggests that they should also be important for computational electrodynamics (CED). Typically DG schemes are coupled with Runge-Kutta timestepping, resulting in RKDG schemes, which are also sometimes called DGTD schemes in the CED community. However, Maxwell's equations, which are solved in CED codes, have global mimetic constraints. In Balsara and Käppeli [von Neumann Stability Analysis of Globally Constraint-Preserving DGTD and PNPM Schemes for the Maxwell Equations using Multidimensional Riemann Solvers, Journal of Computational Physics, 376 (2019) 1108-1137] the authors presented globally constraint-preserving DGTD schemes for CED. The resulting schemes had excellent low dissipation and low dispersion properties. Their one deficiency was that the maximal permissible CFL of DGTD schemes decreased with increasing order of accuracy. The goal of this paper is to show how this deficiency is overcome. Because CED entails the propagation of electromagnetic waves, we would also like to obtain DG schemes for CED that minimize dissipation and dispersion errors even more than the prior generation of DGTD schemes.Two recent advances make this possible. The first advance, which has been reported elsewhere, is the development of a multidimensional Generalized Riemann Problem (GRP) solver. The second advance relates to the use of Two Derivative Runge Kutta (TDRK) timestepping. This timestepping uses not just the solution of the multidimensional Riemann problem, it also uses the solution of the multidimensional GRP. When these two advances are melded together, we arrive at DG(TD)2 schemes for CED, where the first “TD” stands for time-derivative and the second “TD” stands for the TDRK timestepping. The first goal of this paper is to show how DG(TD)2 schemes for CED can be formulated with the help of the multidimensional GRP and TDRK timestepping. The second goal of this paper is to utilize the free parameters in TDRK timestepping to arrive at DG(TD)2 schemes for CED that offer a uniformly large CFL with increasing order of accuracy while minimizing the dissipation and dispersion errors to exceptionally low values. The third goal of this paper is to document a von Neumann stability analysis of DG(TD)2 schemes so that their dissipation and dispersion properties can be quantified and studied in detail.At second order we find a DG(TD)2 scheme with CFL of 0.25 and improved dissipation and dispersion properties; for a second order scheme. At third order we present a novel DG(TD)2 scheme with CFL of 0.2571 and improved dissipation and dispersion properties; for a third order scheme. At fourth order we find a DG(TD)2 scheme with CFL of 0.2322 and improved dissipation and dispersion properties. As an extra benefit, the resulting DG(TD)2 schemes for CED require fewer synchronization steps on parallel supercomputers than comparable DGTD schemes for CED. We also document some test problems to show that the methods achieve their design accuracy.

Power management control strategy using physics-based battery models in standalone PV-battery hybrid systems

Abstract The rechargeable lithium-ion batteries (LiBs) are becoming a technology of choice for energy storage applications, due to its high energy and power density. In this paper, a power management control strategy is proposed for a standalone PV-Battery Energy Storage (BES) hybrid system. The proposed control algorithm tracks the Maximum Power Point (MPP) of the solar-cells while avoiding overcharging of LiBs under different solar radiation and load conditions. The controller has two regulation modes; (a) MPP Tracking mode, (b) battery state of charge limit mode. The standalone PV-BES hybrid system is represented as a system of differential algebraic equations and enables effective implementation of control algorithms. A physics-based single particle model accounting for internal states of a battery is implemented in the simulation and control algorithm. A mixed-order finite difference method with optimal node spacing and strong stability-preserving time-stepping Runge–Kutta scheme is used to solve the battery model. A case study is provided to demonstrate the effectiveness of the proposed strategy using real world data. The study reveals that the proposed power control strategy is robust and meets multiple objectives of standalone PV-BES hybrid systems such as no overcharging, 0% excess output power production, and ensuring no energy is transferred to the dump load.

Numerical Solutions for Unsteady Flows of a Magnetohydrodynamic Jeffrey Fluid Between Parallel Plates Through a Porous Medium

In the present article, the numerical solutions for three fundamental unsteady flows (namely Couette, Poiseuille, and generalized Couette flows) of an incompressible magnetohydrodynamic Jeffrey fluid between two parallel plates through a porous medium are presented using differential quadrature method. The equations governing the flow of Jeffrey fluid are modeled in Cartesian coordinate system. The resulting non-dimensional differential equations are approximated by using a new scheme that is trigonometric B-spline differential quadrature method. The scheme is based on the differential quadrature method in which the weighting coefficients are obtained by using trigonometric B-splines as a set of basis functions. This scheme reduces the equation into the system of first-order ordinary differential equation which is solved by adopting strong stability-preserving time-stepping Runge–Kutta scheme. The effects of the sundry parameters of interest on the velocity profiles are studied and the results are presented through graphs. It is observed that, the velocity increases from the horizontal channel to vertical channel. The velocity is a decreasing function of magnetic parameter. With an increase in time, the velocity increases.

Stationary and transient analyses of a pulley-belt system based on an Eulerian approach

In some engineering industries such as the transportation domain, the belt drives are widely used to transmit the power and the motion. They have a complex dynamic behavior because of the tension forces, the dimensions and the utilization purposes. In this paper, extant models are summarized to present a better understanding of the belt mechanics and the dynamic behavior of pulley-belt systems. Then, a new approach based on the Eulerian description as well as on the normal and shear stresses is introduced in order to discuss the belt mechanics, to build an analytical model of a pulley-belt system and to investigate its dynamic behavior. It can describe the zones of contact between the belt and the pulleys in the stationary and transient regimes. Analytical solutions and formulations for the tensions and the relative displacements of the belt over the pulleys in the stationary regime are established, while differential equations for the belt tensions and the driven pulley velocity in the transient regime are derived. These differential equations are solved using the explicit Runge–Kutta time-step integration scheme. The study is also conducted to observe slip phenomena due to the contact between elastic and rigid surfaces and multiple zones of slip can be found by means of an iterative Newton–Raphson method.

Stability and nonlinear self-excited friction-induced vibrations for a minimal model subjected to multiple coalescence patterns

In certain industrial applications with frictional interfaces such as brake systems, the friction-induced vibrations created by coupled modes can lead to a dynamic instability and thus to an important deterioration in the operating condition. As a result, they are considered as a source of critical engineering problem. In addition, the presence of the nonlinearity makes necessary the consideration of the nonlinear dynamic analysis in order to explain clearly the complexity of the contribution of different frequency components due to unstable modes in the self-excited friction-induced vibrations, to get a design as reliable as possible and to avoid catastrophic failure during the operation phase of the mechanism. The present paper is based on previous works of Sinou and Jézéquel and extends them to include a developed damped four-degree-of-freedom system with frictional contact and spring cubic nonlinearities. Its essential goals are to analyze numerically the mode-coupling instability of the four-degree-of-freedom system owing to the friction between the surfaces of contact and to predict its nonlinear dynamic behavior. The numerical study of stability for the static solution of the mechanical system is performed by applying the complex eigenvalue analysis of the linearized differential equations of motion and by identifying the Hopf bifurcation points as a function of the coefficient of kinetic friction. Depending on the Runge-Kutta time-step integration scheme and the fast Fourier transforms, quantitative and qualitative nonlinear phenomena related to self-excited friction-induced oscillations and limit cycle evolutions are observed and discussed for various friction coefficients.

Von Neumann stability analysis of globally divergence-free RKDG schemes for the induction equation using multidimensional Riemann solvers

In this paper we focus on the numerical solution of the induction equation using Runge–Kutta Discontinuous Galerkin (RKDG)-like schemes that are globally divergence-free. The induction equation plays a role in numerical MHD and other systems like it. It ensures that the magnetic field evolves in a divergence-free fashion; and that same property is shared by the numerical schemes presented here. The algorithms presented here are based on a novel DG-like method as it applies to the magnetic field components in the faces of a mesh. (I.e., this is not a conventional DG algorithm for conservation laws.) The other two novel building blocks of the method include divergence-free reconstruction of the magnetic field and multidimensional Riemann solvers; both of which have been developed in recent years by the first author.Since the method is linear, a von Neumann stability analysis is carried out in two-dimensions to understand its stability properties. The von Neumann stability analysis that we develop in this paper relies on transcribing from a modal to a nodal DG formulation in order to develop discrete evolutionary equations for the nodal values. These are then coupled to a suitable Runge–Kutta timestepping strategy so that one can analyze the stability of the entire scheme which is suitably high order in space and time.We show that our scheme permits CFL numbers that are comparable to those of traditional RKDG schemes. We also analyze the wave propagation characteristics of the method and show that with increasing order of accuracy the wave propagation becomes more isotropic and free of dissipation for a larger range of long wavelength modes. This makes a strong case for investing in higher order methods. We also use the von Neumann stability analysis to show that the divergence-free reconstruction and multidimensional Riemann solvers are essential algorithmic ingredients of a globally divergence-free RKDG-like scheme.Numerical accuracy analyses of the RKDG-like schemes are presented and compared with the accuracy of PNPM schemes. It is found that PNPM retrieve much of the accuracy of the RKDG-like schemes while permitting a larger CFL number.

A note on solving the fourth-order Kuramoto-Sivashinsky equation by the compact finite difference scheme

The present article is concerned with the implementation of the compact finite difference scheme, in the space and the optimal four-stage, order three strong stability-preserving time-stepping Runge-Kutta (SSP-RK43) scheme, in time for computation of one dimensional Kuramoto-Sivashinsky equation (KSE), arises in the study of flame front propagation, phase turbulence in reaction-diffusion system and in many other biological and chemical processes. The efficiency of proposed scheme is confirmed by six test problems with known exact solutions. The numerical results demonstrate the reliability and efficiency of the algorithm developed.

Extended modified cubic B-spline algorithm for nonlinear Burgers' equation

In this paper, an extended modified cubic B-Spline differential quadrature method is proposed to approximate the solution of the nonlinear Burgers' equation. The proposed method is used in space and a five-stage and four order strong stability-preserving time-stepping Runge–Kutta (SSP-RK54) method is used in time. The accuracy and efficiency of the method is illustrated by considering four numerical problems. The numerical results of the method are compared with some existing methods and it was found that the proposed numerical method produces acceptable results and even more accurate results in comparison with some existing methods. The stability analysis of the scheme is also carried out and was found to be unconditionally stable.

Hp discontinuous Galerkin methods for the vertical extent of the water column in coastal settings part I: Barotropic forcing

In this article, we present novel, high-order, discontinuous Galerkin (DG) methods for the vertical extent of the water column in coastal settings. We examine the shallow water equations (SWE) in the context of DG spatial discretizations coupled with explicit Runge-Kutta (RK) time stepping. All the primary variables, including the free surface elevation, are discretized using discontinuous polynomial spaces of arbitrary order. The difficulty of mismatched lateral boundary faces that accompanies the use of a discontinuous free surface is overcome through the use of a so-called sigma-coordinate system in the vertical, which transforms the bottom boundary and free surface into coordinate surfaces. We develop high-order methods for the SWE that exhibit optimal orders of convergence for all the primary variables via two distinct paths: the first involves the use of a convolution kernel made up of B-splines to filter out errors in the DG discretization of the surface elevation and the corresponding pressure flux. The second involves a method that evaluates the discrete depth-integrated velocity exactly, eliminating the need to solve the depth-integrated momentum equation altogether. The result is a simple and efficient high-order scheme that can be extended to the full three-dimensional SWE.