Time Step Size Research Articles

Adaptive time step selection for spectral deferred correction

AbstractSpectral Deferred Correction (SDC) is an iterative method for the numerical solution of ordinary differential equations. It works by refining the numerical solution for an initial value problem by approximately solving differential equations for the error, and can be interpreted as a preconditioned fixed-point iteration for solving the fully implicit collocation problem. We adopt techniques from embedded Runge-Kutta Methods (RKM) to SDC in order to provide a mechanism for adaptive time step size selection and thus increase computational efficiency of SDC. We propose two SDC-specific estimates of the local error that are generic and do not rely on problem specific quantities. We demonstrate a gain in efficiency over standard SDC with fixed step size and compare efficiency favorably against state-of-the-art adaptive RKM.

Particle Virtual Element Method (PVEM): an agglomeration technique for mesh optimization in explicit Lagrangian free-surface fluid modelling

Explicit solvers are commonly used for simulating fast dynamic and highly nonlinear engineering problems. However, these solvers are only conditionally stable, requiring very small time-step increments determined by the characteristic length of the smallest, and often most distorted, element in the mesh. In the Lagrangian description of fluid motion, the computational mesh quickly deteriorates. To circumvent this problem, the Particle Finite Element Method (PFEM) creates a new mesh (e.g., through a Delaunay tessellation, based on node positions) when the current one becomes overly distorted. A fast and efficient remeshing technique is therefore of pivotal importance for an effective PFEM implementation in explicit dynamics. Unfortunately, the 3D Delaunay tessellation does not guarantee well-shaped elements, often generating zero- or near-zero-volume elements (slivers), which drastically reduce the stable time-step size. Available mesh optimization techniques have limited applicability due to their high computational cost when runtime remeshing is required. An innovative possibility to overcome this problem is the use of the Virtual Element Method (VEM), a variant of the finite element method that can make use of polyhedral elements of arbitrary shapes and number of nodes. This paper presents the formulation of a 3D first-order Particle Virtual Element Method (PVEM) for weakly compressible flows. Starting from a tetrahedral mesh, poorly shaped elements, such as slivers, are agglomerated to form polyhedral Virtual Elements (VEs) with a controlled characteristic length. This approach ensures full control over the minimum time-step size in explicit dynamics simulations, maintaining stability throughout the entire analysis.

Domain decomposition with local time discretization for the nonlinear Stokes–Biot system

This work presents a domain decomposition method for the fluid-poroelastic structure interaction (FPSI) system, which utilizes local time integration for subproblems. To derive the domain decomposition scheme, we introduce a Lagrange multiplier and define time-dependent Steklov–Poincaré-type operators based on the interface conditions. These operators are employed to transform the coupled system into an evolutionary nonlinear interface problem, which is then solved using an iterative algorithm. This approach provides the flexibility to use different time discretization schemes and step sizes in subdomains, making it an efficient method for simulating multiphysics systems. We present numerical tests for both non-physical and physical problems to demonstrate the accuracy and efficiency of this method.

Stability and computational analysis of Influenza-A epidemic model through double time delay

Delay factors demonstration has a significant role in controlling a strain of infectious disease instead of a pharmaceutical strategy. According to the World Health Organization (WHO), 3–5 million cases are reported annually and approximately 290,000 to 650,000 respiratory deaths annually. So, in the present study, we develop a delayed mathematical model based on delay differential equations (DDEs) for the influenza epidemic using a deterministic approach by introducing double delay parameters. The four distinct sub-populations are considered susceptible, exposed, infected, and recovered. For the rigorous analysis, the fundamental properties of the model like positivity, boundedness, existence, and uniqueness, were studied. The influenza-free equilibrium (IFE) and influenza-existing equilibrium (IEE), are the two nonnegative equilibrium points that the model demonstrates. Both locally and globally, the asymptotic stability of the equilibrium points of the model is established and shown under specific situations of reproduction number. Additionally, investigated the model's parameter sensitivity and determined the relative sensitivity of each parameter. Both standard and nonstandard methods—such as Euler, Runge-Kutta, and nonstandard finite difference with a delayed sense—are presented to make computational analysis support a dynamical analysis and the best visualization of results. The stability of the non-standard finite difference scheme is thoroughly analyzed around the steady states of the model. Additionally, the results show that the nonstandard finite difference approximation is an efficient, cost-effective method, independent of time step size, to solve such highly nonlinear and complex real-world problems.

A Novel Approach Integrating the Method of Characteristics with Large Time-Step Scheme (MOC-LTS) for Efficient Transient Flow Simulation in Liquid Pipelines

Water hammer incidents pose a significant risk to the safety and stability of pipeline operations. Therefore, rapid and precise transient flow simulation is essential for efficiently developing scientific water hammer control strategies. Nevertheless, the prevailing transient flow simulation methods for liquid pipelines predominantly employ explicit schemes to solve transient flow control equations, necessitating adherence to the stability criterion of Courant-Friedrichs-Lewy (CFL) ≤ 1. This results in limited time step sizes, which in turn constrains computational efficiency, particularly when simulating hydraulic behavior in large-scale pipeline networks. In this paper, a novel approach integrating the method of characteristics (MOC) with a large time-step scheme (LTS) is proposed to enable rapid and accurate simulation of transient flow in liquid pipelines. The proposed approach discretizes the computational domain into contiguous control volumes, ategorizing them as either boundary or internal control volumes depending on whether they are affected by boundary conditions within a large time step. For internal control volumes, an LTS scheme based on the first-order Godunov format is employed to improve computational efficiency. For boundary control volumes, MOC is applied iteratively to accurately capture boundary dynamics until the simulated time matches that of the internal control volumes. Quantitative analysis is conducted to verify the performance of the proposed approach. The results confirm that the proposed approach outperforms the MOC in computational efficiency while maintaining minimal accuracy loss.

Error Estimates for First- and Second-Order Lagrange–Galerkin Moving Mesh Schemes for the One-Dimensional Convection–Diffusion Equation

A new moving mesh scheme based on the Lagrange–Galerkin method for the approximation of the one-dimensional convection–diffusion equation is studied. The mesh movement is prescribed by a discretized dynamical system for the nodal points. This system is related to the velocity and diffusion coefficient in the convection–diffusion equation such that the nodal points follow the convective flow of the model. It is shown that under a restriction of the time step size the mesh movement cannot lead to an overlap of the elements and therefore an invalid mesh. Using a piecewise linear approximation, optimal error estimates in the ℓ∞(L2)∩ℓ2(H01)\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\\ell ^\\infty (L^2) \\cap \\ell ^2(H_0^1)$$\\end{document} norm are proved in case of both, a first-order backward Euler method and a second-order two-step method in time. These results are based on new estimates of the time dependent interpolation operator derived in this work. Preservation of the total mass is verified for both choices of the time discretization. Numerical experiments are presented that confirm the error estimates and demonstrate that the proposed moving mesh scheme can circumvent limitations that the Lagrange–Galerkin method on a fixed mesh exhibits.

Stability of step size control based on a posteriori error estimates

A posteriori error estimates based on residuals can be used for reliable error control of numerical methods. Here, we consider them in the context of ordinary differential equations and Runge-Kutta methods. In particular, we take the approach of Dedner & Giesselmann (2016) and investigate it when used to select the time step size. We focus on step size control stability when combined with explicit Runge-Kutta methods and demonstrate that a standard I controller is unstable while more advanced PI and PID controllers can be designed to be stable. We compare the stability properties of residual-based estimators and classical error estimators based on an embedded Runge-Kutta method both analytically and in numerical experiments.

An automatic simulation pipeline for coupled simulations of acoustic damping materials

AbstractFoamed materials are widely used to reduce noise due to their comparably good acoustic damping behavior. However, out of a large variety of these materials a suitable candidate has to be identified for each application. This is a challenging process that is typically guided by experiments and experience. While numerical simulations could support these experiments and reduce the effort to a great extent, no suitable discretization approach has yet been established that can fully capture the complex geometry of the foam. A fully resolved model is desirable in order to yield reliable predictions that can then be used to establish homogenized models. We established a monolithic coupling approach based on the finite cell method (FCM) that realizes a vibroacoustic simulation in this sense. The fluid and the structure domain are discretized by Cartesian grids and the geometry defined based on computed tomography scans is accounted for during the quadrature of the weak form. Our simulation in the time domain makes use of explicit time marching schemes and is therefore limited by a critical time step size. This is known to be arbitrarily low for discretizations with the FCM containing cells with arbitrarily small support. As a remedy against this we use the classical ‐stabilization technique and investigate its potentials and limitations.

A convection-diffusion-reaction system with discontinuous flux modelling biofilm growth in slow sand filters

A Slow Sand Filter (SSF) consists of biofilm attached to saturated packed sand with supernatant water on top, through which water flows by gravity. SSFs are used in the production of drinking water to improve water quality and remove pathogens, and, as they use biological filtration via the microbes present in the biofilm, they are simple and cheap to construct. A drawback of SSFs is that the difficult sampling and the complex dynamics associated with biofilm growth of microbes pose an obstacle for analysis and predictions. In order to overcome this issue, a mathematical model consisting of a non-linear convection-diffusion-reaction system with discontinuous flux is developed from first principles as a unified framework for the dynamics in time and depth of the biofilm and microbial community, where a variety of ecological models can be included in the reaction terms, which also contain attachment, detachment and transfer processes. The multi-phase approach considers three main physical volumes: a biofilm matrix, the water volume enclosed by it and the flowing suspension through the SSF surrounding the biofilm. The only input data needed are the inflow concentrations and rate, light irradiation and temperature. The model includes a Cahn–Hilliard-type equation for the biofilm velocity in the supernatant water. A numerical scheme that preserves physical properties of the solution is also presented. This is achieved with a time-splitting scheme using an upwind scheme with an adaptive time-step size that ensures positivity of the solution under a CFL condition and a reasonable assumption on the separate numerical solution of the Cahn–Hilliard-type equation. The latter is solved numerically with a semi-implicit finite-difference scheme using a mixed form of the equation. The flexible model framework with its numerical scheme allows for including a variety of ecological models, making the investigation and development of increasingly complex simulations possible. Simulations with different model parameters are shown to be qualitatively consistent with the literature.

Residual network with self-adaptive time step size

Residual Networks (ResNet) are pivotal in machine learning. The connection between ResNets and ordinary differential equations (ODEs) has inspired enhancements of ResNets using sophisticated numerical methods for ODE systems. Recent advancements in numerical self-adaptive schemes, which adjust time step sizes based on the feature maps or parameters of residual blocks, have demonstrated promising results in enhancing ResNet performance, surpassing those achieved with fixed time step methods. However, these self-adaptive time step constructions lack theoretical support and can limit performance improvements since the self-adaptive time steps should theoretically depend on both the feature maps and the parameters of the residual blocks. In this study, we conduct a rigorous theoretical analysis of the residual functions associated with fixed or self-adaptive time step methods to demonstrate the advantages and rational designs of the self-adaptive approach. Subsequently, we introduce a novel self-adaptive ResNet, AdaTS-ResNet, which effectively incorporates both feature maps and parameters in its time step adjustments. Experimental results reveal that AdaTS-ResNet surpasses TSCLSTM-ResNet (Yang et al., 2020) in prediction accuracy and computational efficiency, where TSCLSTM represents the latest advancements of the methods designed based on the self-adaptive scheme of ODE. Our findings highlight the potential of improving ResNet architectures through adaptive techniques of dynamical systems, offering insights for future enhancements in deep learning models.

An improved solution of the probability density evolution equation for analyzing stochastic structural responses with adaptive time steps and initial conditions

The probability density evolution method (PDEM) offers an innovative approach for analyzing the stochastic vibration responses of structures. Given the complexities inherent in real-world challenges, the probability density evolution equations (PDEE) often require numerical methods for effective resolution. Consequently, improving computational efficiency and solution accuracy is of paramount importance. Currently, the impulse function discretization method (IFDM) faces several challenges, including limited differentiability in the spatial direction for initial value conditions. Furthermore, the probability density function produced by this method frequently deviates significantly from the true value, undermining the accuracy of the PDEE. Additionally, when addressing vibration responses characterized by rapid evolutionary velocities, the use of a globally fixed uniform time step size (FUTSS) results in a marked decrease in computational efficiency. To tackle these challenges, this study introduces an improved probability density evolution method (IPDEM), which integrates the Gaussian kernel density estimation discretization method (GKDEDM) and an adaptive non-uniform time step size (ANUTSS). The advantages of the IPDEM are illustrated through two numerical simulation examples. The findings indicate that the GKDEDM not only enhances spatial differentiability but also provides a more accurate estimation of the probability density at the initial moment, thereby improving the accuracy of the PDEE. Furthermore, the ANUTSS significantly reduces the number of iterations required while preserving precision, leading to a substantial optimization of computation time.

A parallel coupling framework for DEM-MBD: Model verification and application

Bulk material handling equipment can be numerically studied using the DEM-MBD method, where the particulate dynamics are solved by DEM (Discrete Element Method) and the motion characteristics of the mechanical components are modeled by MBD (Multi-Body Dynamics). In this work, a two-way coupling framework is proposed for the co-simulation of our self-developed DEM solver (DEMSLab) with the commercial MBD software (RecurDyn) using a parallel staggered coupling approach. To verify the proposed coupling method, the simulation result of a torsional spring damper system is compared with the corresponding analytical solution. After that, to assess the feasibility of this method in industrial-scale applications, a series of simulations involving a wheel loader is carried out. Different simulation conditions are adopted to reveal their impact on computational efficiency and simulation accuracy. These conditions include the mesh representation of the equipment geometry, the MBD time step size, and the coupling interval. An improved mesh simplification algorithm is proposed and used to reasonably simplify the surface mesh of the wheel loader, thus lowering the computational intensity of particle-wall interactions in DEM. The comparison results indicate that the computational time can be largely saved while the simulation accuracy remains considerably high. The MBD step size and coupling interval are crucial for numerical stability and should be selected carefully.

Exponential time propagators for elastodynamics

We propose a computationally efficient and systematically convergent approach for elastodynamics simulations. We recast the second-order dynamical equation of elastodynamics into an equivalent first-order system of coupled equations, so as to express the solution in the form of a Magnus expansion. With any spatial discretization, it entails computing the exponential of a matrix acting upon a vector. We employ an adaptive Krylov subspace approach to inexpensively and accurately evaluate the action of the exponential matrix on a vector. In particular, we use an apriori error estimate to predict the optimal Krylov subspace size required for each time-step size. We show that the Magnus expansion truncated after its first term provides quadratic and superquadratic convergence in the time-step for nonlinear and linear elastodynamics, respectively. We demonstrate the accuracy and efficiency of the proposed method for one linear (linear cantilever beam) and three nonlinear (nonlinear cantilever beam, soft tissue elastomer, and hyperelastic rubber) benchmark systems. For a desired accuracy in energy, displacement, and velocity, our method allows for 10−100× larger time-steps than conventional time-marching schemes such as Newmark-β method. Computationally, it translates to a ∼1000× and ∼10−100× speed-up over conventional time-marching schemes for linear and nonlinear elastodynamics, respectively.

Exponential time difference methods with spatial exponential approximations for solving boundary layer problems

AbstractIn this work, an efficient and stable exponential time difference method is presented for solving boundary layer problems. By combining exponential time difference schemes with spatial direct discontinuous Galerkin discretization based on exponential boundary layer approximations, the proposed algorithm not only may admit large time step sizes but also could provide good spatial approximations even if on rather coarse spatial grids in the boundary layer. Some energy stabilities of the numerical scheme are rigorously derived. Numerical examples illustrate the accuracy, stability and efficiency of the algorithm.

Non-polynomial spline method for computational study of reaction diffusion system

This work addresses an efficient and new numerical technique utilizing non-polynomial splines to solve system of reaction diffusion equations (RDS). These system of equations arise in pattern formation of some special biological and chemical reactions. Different types of RDS are in the form of spirals, hexagons, stripes, and dissipative solitons. Chemical concentrations can travel as waves in reaction-diffusion systems, where wave like behaviour can be seen. The purpose of this research is to develop a stable, highly accurate and convergent scheme for the solution of aforementioned model. The method proposed in this paper utilizes forward difference for time discretization whereas for spatial discretization cubic non-polynomial spline is used to get approximate solution of the system under consideration. Furthermore, stability of the scheme is discussed via Von-Neumann criteria. Different orders of convergence is achieved for the scheme during a theoretical convergence test. Suggested method is tested for performance on various well known models such as, Brusselator, Schnakenberg, isothermal as well as linear models. Accuracy and efficiency of the scheme is checked in terms of relative error (E R ) and L ∞ norms for different time and space step sizes. The newly obtained results are analyzed and compared with those available in literature.

A second order difference method combined with time two-grid algorithm for two-dimensional time-fractional Fisher equation

In this paper, a finite difference (FD) scheme based on time two-grid algorithm is proposed for solving the two-dimensional time-fractional Fisher equation (2D-TFFE). Firstly, the Caputo fractional derivative and the spatial derivative are discretized by the L 2 − 1 σ formula and the central difference formula, respectively. In order to improve the efficiency of computation, the time two-grid algorithm is then constructed. Secondly, based on the cut-off function technic, stability and convergence of the time two-grid FD scheme are obtained by the energy method, and the global convergence order is O ( τ F 2 + τ C 4 + h x 2 + h y 2 ) , where τ F and τ C represent time-step sizes on the fine and coarse grid, respectively, while h x and h y represent the space-step sizes. Finally, numerical experiments are presented to show the feasibility and efficiency of the algorithm.

The impact of step sizes on pressure prediction in fracturing treatment via deep learning algorithms

Hydraulic fracturing relies on accurate pressure prediction for effective risk management and treatment evaluation. Deep learning models such as artificial neural networks, convolutional neural networks, recurrent neural networks, long short-term memory networks, gated recurrent units, and CNN-LSTM offer promising avenues for this prediction. However, the choice of time step size significantly influences the models’ predictive power and engineering relevance. Understanding how different models respond to varying time step sizes remains a challenge. This study systematically compares these models for both single-step and multistep pressure predictions in hydraulic fracturing scenarios. By examining the impact of time step size on prediction accuracy and lag, the research sheds light on the tradeoffs between model complexity, dataset size, and prediction performance. The integration of multiple models aims to leverage their individual strengths while mitigating weaknesses, potentially enhancing overall prediction accuracy. This comprehensive analysis aims to inform practitioners and researchers about the optimal selection and deployment of deep learning models for pressure prediction, thus advancing the effectiveness and reliability of hydraulic fracturing operations.

Efficient simulation of complex Ginzburg–Landau equations using high-order exponential-type methods

In this paper, we consider the task of efficiently computing the numerical solution of evolutionary complex Ginzburg–Landau equations on Cartesian product domains with homogeneous Dirichlet/Neumann or periodic boundary conditions. To this aim, we employ for the time integration high-order exponential methods of splitting and Lawson type with constant time step size. These schemes enjoy favorable stability properties and, in particular, do not show restrictions on the time step size due to the underlying stiffness of the models. The needed actions of matrix exponentials are efficiently realized by using a tensor-oriented approach that suitably employs the so-called μ-mode product (when the semidiscretization in space is performed with finite differences) or with pointwise operations in Fourier space (when the model is considered with periodic boundary conditions). The overall effectiveness of the approach is demonstrated by running simulations on a variety of two- and three-dimensional (systems of) complex Ginzburg–Landau equations with cubic or cubic-quintic nonlinearities, which are widely considered in literature to model relevant physical phenomena. In fact, we show that high-order exponential-type schemes may outperform standard techniques to integrate in time the models under consideration, i.e., the well-known second-order split-step method and the explicit fourth-order Runge–Kutta integrator, for stringent accuracies.

Energy-conserving SAV-Hermite–Galerkin spectral scheme with time adaptive method for coupled nonlinear Klein–Gordon system in multi-dimensional unbounded domains

In this article, we construct an efficient numerical scheme to solve the coupled nonlinear Klein–Gordon system in a multi-dimensional unbounded domain Rd(d=1,2,3). We first use the SAV method to equivalently deform the original system, then we treat the nonlinear part of the system explicitly in temporal discretization. Thanks to the SAV method, we ensure the unconditional energy conservation of the discrete system, as a cost, we only need to solve two completely linearized decoupled systems. To avoid the errors caused by domain truncation and the construction of artificial boundary conditions, we apply the Hermite–Galerkin spectral method with scaling factor for spatial discretization. In order to save computational costs, we utilize two time adaptive strategies to adjust the time step size, through which we can improve computational efficiency without sacrificing accuracy. Finally, we presented numerical experiments to demonstrate the accuracy, efficiency, and energy conservation properties of our algorithm, and applied it to simulate the interactions of two-dimensional and three-dimensional vector solitons.

Increasing numerical stability of mountain valley glacier simulations: implementation and testing of free-surface stabilization in Elmer/Ice

Abstract. This paper concerns a numerical stabilization method for free-surface ice flow called the free-surface stabilization algorithm (FSSA). In the current study, the FSSA is implemented into the numerical ice-flow software Elmer/Ice and tested on synthetic two-dimensional (2D) glaciers, as well as on the real-world glacier of Midtre Lovénbreen, Svalbard. For the synthetic 2D cases it is found that the FSSA method increases the largest stable time-step size at least by a factor of 5 for the case of a gently sloping ice surface (∼ 3°) and by at least a factor of 2 for cases of moderately to steeply inclined surfaces (∼ 6° to 12°) on a fine mesh. Compared with other means of stabilization, the FSSA is the only one in this study that increases largest stable time-step sizes when used alone. Furthermore, the FSSA method increases the overall accuracy for all surface slopes. The largest stable time-step size is found to be smallest for the case of a low sloping surface, despite having overall smaller velocities. For an Arctic-type glacier, Midtre Lovénbreen, the FSSA method doubles the largest stable time-step size; however, the accuracy is in this case slightly lowered in the deeper parts of the glacier, while it increases near edges. The implication is that the non-FSSA method might be more accurate at predicting glacier thinning, while the FSSA method is more suitable for predicting future glacier extent. A possible application of the larger time-step sizes allowed for by the FSSA is for spin-up simulations, where relatively fast-changing climate data can be incorporated on short timescales, while the slow-changing velocity field is updated over larger timescales.