Time Step Size Research Articles

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.

Harmonic Oscillator Staging Coordinates for Efficient Path Integral Simulations of Quantum Oscillators and Crystals.

Imaginary-time path integral (PI) is a rigorous quantum mechanical tool to compute static properties at finite temperatures. However, the stiff nature of the internal PI modes poses a sampling challenge. This is commonly tackled using staging coordinates, in which the free particle (FP) contribution of the PI action is diagonalized. We introduce novel and simple staging coordinates that diagonalize the entire action of the harmonic oscillator (HO) model, rendering it efficiently applicable to (exclusively) systems with harmonic character, such as quantum oscillators and crystals. The method is not applicable to fluids or systems with imaginary modes. Unlike FP staging, the HO staging provides a unique treatment of the centroid mode. We provide implementation schemes for PIMC and PIMD simulations in NVT ensemble. Sampling efficiency is assessed in terms of the precision and accuracy of estimating the energy and heat capacity of a one-dimensional HO and an asymmetric anharmonic oscillator (AO). In PIMC, the HO coordinates propose collective moves that perfectly sample the HO contribution, then (for AO) the residual anharmonic term is sampled using standard Metropolis method. This results in a high acceptance rate and, hence, high precision, in comparison to the FP staging. In PIMD, the HO coordinates naturally prescribe definitions for the fictitious masses, yielding equal frequencies of all modes when applied to the HO model. This allows for a substantially larger time step sizes relative to standard staging, without affecting accuracy or integrator stability. For completeness, we also present results using normal mode (NM) coordinates, based on both HO and FP models. While staging and NM coordinates show similar performance (for FP or HO), staging is computationally preferable due to its cheaper scaling with the number of beads. The simplicity and the enhanced sampling gained by the HO coordinates open avenues for efficient estimation of nuclear quantum effects in more complex systems with harmonic character, such as real molecular bonds and quantum crystals.

[formula omitted]-adaptive discontinuous Galerkin solution of transonic viscous flows with variable time step-size

The discontinuous Galerkin methods demonstrated to be well suited for scale resolving simulations of complex configurations, characterized by different fluid dynamics phenomena, such as transition, separation, shock/boundary layer interaction. Unfortunately, solvers based on these methods cannot yet reach the computational efficiency of well-established standard solvers, e.g., based on finite volume methods. To reduce the computational cost and the memory footprint, while not spoiling the accuracy, different approaches were investigated, which act on the spatial or the temporal discretization, and on the linear algebra. The above approaches have been extensively investigated, but rarely considering their mutual interactions. In this work, a p-adaptation algorithm is coupled with a time step-size adaptation algorithm to mitigate robustness issues arising after each adaptation cycle, probably related to the projection of the old solution on the new polynomial orders distribution. Moreover, this strategy is able to control the transient phase before each adaptation cycle, and can handle also the presence of initial steady solution, when the mesh is too coarse and the polynomial degree has not yet compensated for the lack of spatial accuracy. The effectiveness and accuracy of the proposed algorithm are demonstrated on transonic viscous flows with shock-wave/boundary-layer interaction.

Surface Green's functions for an anisotropic viscoelastic half-plane and their application to contact problems

In this paper, the elastic-like surface Green's functions for an anisotropic viscoelastic half-plane are derived using the time-stepping method. Using the elastic-like surface Green's functions as the core analytical solutions, we develop semi-analytical models (SAMs) and apply them to solve two different contact problems with anisotropic viscoelastic materials. As new modeling approaches, the SAMs developed here can provide fast and efficient approaches to solving contact problems. These methods enable us to consider contact problems with generally anisotropic viscoelastic solids, in which the contact surface is frictional and either smooth or rough, and the applied loads and boundaries can be time-variant. The correctness of the derived surface Green's functions is demonstrated by comparing the numerical results obtained by SAMs and those achieved from the analytical solutions or boundary element methods. Using the obtained numerical results, the impacts of time step size, anisotropy, frictional coefficient, roughness, and applied loads on the contact responses are further analyzed and discussed.

Second-order error analysis of a corrected average finite difference scheme for time-fractional Cable equations with nonsmooth solutions

In this paper, we first formulate an average finite difference scheme with appropriate correction terms for the time-fractional initial-value problem of y′(t)+0RLDt1−αy(t)=g(t), where 0RLDt1−α denotes the Riemann–Liouville fractional derivative of order 1−α(α∈(0,1)). It is shown that, under some suitable conditions on the data, the numerical solution converges to the exact solution with order O(τ2), where τ is the time step size. The method is then extended to solve the time-fractional Cable equation, combined with a standard discretisation of the spatial derivatives on a uniform mesh. The second-order time convergence rate is proved. Several numerical examples are given to illustrate the good agreement with the theoretical analysis of the presented methods.

On the convergence of a linearly implicit finite element method for the nonlinear Schrödinger equation

AbstractWe consider a model initial‐ and Dirichlet boundary–value problem for a nonlinear Schrödinger equation in two and three space dimensions. The solution to the problem is approximated by a conservative numerical method consisting of a standard conforming finite element space discretization and a second‐order, linearly implicit time stepping, yielding approximations at the nodes and at the midpoints of a nonuniform partition of the time interval. We investigate the convergence of the method by deriving optimal‐order error estimates in the and the norm, under certain assumptions on the partition of the time interval and avoiding the enforcement of a Courant‐Friedrichs‐Lewy (CFL) condition between the space mesh size and the time step sizes.

A hybrid shape functions based temporal finite element method to solve transient heat conduction problems

A novel temporal finite element method is presented to solve Fourier and non-Fourier heat conduction problems. Two kinds of temporal finite element models are developed using Gurtin variational principle and weighted residual technique, respectively. A kind of hybrid shape function with polynomial and trigonometric basis is stressed to give more flexible and appropriate description of temporal behavior of temperature. A recursive algorithm is developed by which the temperature solution at a specific time can be obtained only via a matrix power product with the initial condition, and a criterion of stability analysis is derived, which is can numerically be conducted when the constitution of shape function and time step size are prescribed. The proposed approach is available for Fourier and non-Fourier heat transfer problems, and can conveniently be combined with well-developed numerical algorithms for the boundary value problem, such as FE, SBFEM, etc. Various numerical examples, including those with the nonlinear thermal conductivity, radiative boundary condition, etc. are provided to illustrate the efficiency of proposed approaches, and impacts of the temporal FE model, constitution of shape functions, and step size, etc. are taken into account. Satisfactory results are achieved in comparison with the analytical and the ABAQUS based numerical solutions.

Numerical Computation of 2D Domain Integrals in Boundary Element Method by (α, β) Distance Transformation for Transient Heat Conduction Problems

When the time-dependent boundary element method, also termed the pseudo-initial condition method, is employed for solving transient heat conduction problems, the numerical evaluation of domain integrals is necessitated. Consequently, the accurate calculation of the domain integrals is of crucial importance for analyzing transient heat conduction. However, as the time step decreases progressively and approaches zero, the integrand of the domain integrals is close to singular, resulting in large errors when employing standard Gaussian quadrature directly. To solve the problem and further improve the calculation accuracy of the domain integrals, an (α, β) distance transformation is presented. Distance transformation is a simple and efficient method for eliminating near-singularity, typically applied to nearly singular integrals. Firstly, the (α, β) coordinate transformation is introduced. Then, a new distance transformation for the domain integrals is constructed by replacing the shortest distance with the time step. With the new method, the integrand of the domain integrals is substantially smoothed, and the singularity arising from small time steps in the domain integrals is effectively eliminated. Thus, more accurate results can be obtained by the (α, β) distance transformation. Different sizes of time steps, positions of source point, and shapes of integration elements are considered in numerical examples. Comparative studies of the numerical results for the domain integrals using various methods demonstrate that higher accuracy and efficiency are achieved by the proposed method.

Coupled modeling and simulation of tank self-pressurization and thermal stratification

Numerical analysis of thermo-fluid dynamics for cryogenic propellant storage primarily consists of nodal and computational fluid dynamics (CFD) modeling approaches. While nodal approaches prioritize faster computation over fidelity, CFD approaches promise accuracy and fidelity at the expense of significant computational resources. This paper presents an efficient coupling methodology developed to facilitate the simulation of a long-term self-pressurization process of a cryogenic propellant tank as well as associated thermal stratification phenomenon under both normal and reduced gravity conditions. Key highlight of the methodology is the development of a coupling scheme based on a domain decomposition approach that separates the tank into two computational regions at the liquid-vapor interface. SINDA/FLUINT, the nodal code, is utilized to simulate the liquid region, while ANSYS Fluent, the CFD code, handles the vapor region. A data exchange algorithm was developed to compute required boundary conditions for each region using the local thermodynamic properties of assumed infinitely thin interface. The coupling between the nodal and CFD codes was demonstrated by simulating a self-pressurization process in a small size tank using hydrogen as the working fluid. A sensitivity study on the grid sizing and coupling time step was conducted to determine appropriate spatial and ensure a divergent-free explicit data exchange time step size. Using the coupling approach, two numerical case studies were performed to study tank self-pressurization by considering two initial hydrogen fill levels. The effectiveness and efficiency of the coupling methodology was assessed by comparing the temperature and pressure evolution results of the coupled simulation with those obtained solely from the CFD simulation. Results showed that the coupling scheme and data exchange algorithm were successfully implemented, and the coupled simulation agreed well with the CFD simulations. In addition, a stratification number was used to characterize the transient dynamic development of the thermal stratification. Lastly, the computational costs associated with both approaches were compared. Significant reductions in simulation time were achieved for both cases.

Unconditional optimal [formula omitted]-norm error estimate and superconvergence analysis of a linearized nonconforming finite element variable-time-step BDF2 method for the nonlinear complex Ginzburg–Landau equation

In this paper, a linearized fully discrete scheme combining the variable-time-step two-step backward differentiation formula (VSBDF2) in time and the nonconforming finite element methods (FEMs) in space is constructed and analyzed for the nonlinear complex Ginzburg–Landau equation. A novel convergence analysis approach is proposed, which shows that the H1-norm error of this scheme can reach optimal order in space and sharp second-order convergence in time. The key ingredients of our analysis are temporal–spatial error splitting approach, a new positive definiteness property of discrete orthogonal convolution (DOC) kernels with respect to complex sequences, an energy projection operator and a discrete Laplace operator Δh related to the nonconforming element space, and a relationship between Δh and classical L2 projection operator in complex space. Furthermore, through a specific treatment of the error equation and the combination technique of interpolating operator and projection operator, we prove that the proposed scheme has the global superconvergence result of order O(τ2+h2) in H1-norm for the first time, which further improves the efficiency of numerical computation. Here τ and h are the maximum time and spatial step sizes, respectively. It should be pointed out that all these error results are rigorously established under mild assumptions of the time step, and get rid of the usual grid ratio constraints between temporal and spatial step sizes in linearized scheme. At last, numerical experiments are provided to verify the theoretical analysis.

Analytical and numerical study of diffusion propelled surface growth phenomena

To understand the complex problem of surface growth phenomena, we developed a model in which the regular diffusion equation is coupled to the Kardar–Parisi–Zhang (KPZ) equation. The fundamental or Gaussian solution of the regular diffusion equation is considered as an external noise or source term in the KPZ equation. The obtained system of partial differential equations is analytically solved by a self-similar Ansatz and expressed the solution as a combination of elementary and special functions. Using this solution, the effects of the physical parameters were explicitly investigated. Our recent explicit numerical method, the leapfrog-hopscotch algorithm, is also tested for this problem. It is shown that this method can be safely used with orders of magnitude larger time step sizes than the usual explicit (Euler) scheme as well if we are far from the cusp-like solutions. We pointed out that the cusps themselves cannot be properly simulated by any method that we know.

Temporal consistency loss for physics-informed neural networks

Physics-informed neural networks (PINNs) have been widely used to solve partial differential equations (PDEs) in a forward and inverse manner using neural networks. However, balancing individual loss terms can be challenging, mainly when training these networks for stiff PDEs and scenarios requiring enforcement of numerous constraints. Even though statistical methods can be applied to assign relative weights to the regression loss for data, assigning relative weights to equation-based loss terms remains a formidable task. This paper proposes a method for assigning relative weights to the mean squared loss terms in the objective function used to train PINNs. Due to the presence of temporal gradients in the governing equation, the physics-informed loss can be recast using numerical integration through backward Euler discretization. The physics-uninformed and physics-informed networks should yield identical predictions when assessed at corresponding spatiotemporal positions. We refer to this consistency as “temporal consistency.” This approach introduces a unique method for training physics-informed neural networks (PINNs), redefining the loss function to allow for assigning relative weights with statistical properties of the observed data. In this work, we consider the two- and three-dimensional Navier–Stokes equations and determine the kinematic viscosity using the spatiotemporal data on the velocity and pressure fields. We consider numerical datasets to test our method. We test the sensitivity of our method to the timestep size, the number of timesteps, noise in the data, and spatial resolution. Finally, we use the velocity field obtained using particle image velocimetry experiments to generate a reference pressure field and test our framework using the velocity and pressure fields.

Asymptotically compatible energy of variable-step fractional BDF2 scheme for the time-fractional Cahn–Hilliard model

Abstract A novel discrete gradient structure of the variable-step fractional BDF2 formula approximating the Caputo fractional derivative of order $\alpha \in (0,1)$ is constructed by a local-nonlocal splitting technique, that is, the fractional BDF2 formula is split into a local part analogue to the two-step backward differentiation formula (BDF2) of the first derivative and a nonlocal part analogue to the L1-type formula of the Caputo derivative. Then a local discrete energy dissipation law of the variable-step fractional BDF2 implicit scheme is established for the time-fractional Cahn–Hilliard model under a weak step-ratio constraint $0.3960\le \tau _{k}/\tau _{k-1}<r^{*}(\alpha )$, where $\tau _{k}$ is the $k$th time-step size and $r^{*}(\alpha )\ge 4.660$ for $\alpha \in (0,1)$. The present result provides a practical answer to the open problem in [SINUM, 57: 218-237, Remark 6] and significantly relaxes the severe step-ratio restriction [Math. Comp., 90: 19–40, Theorem 3.2]. More interestingly, the discrete energy and the corresponding energy dissipation law are asymptotically compatible with the associated discrete energy and the energy dissipation law of the variable-step BDF2 method for the classical Cahn–Hilliard equation, respectively. To the best of our knowledge, such type energy dissipation law is established at the first time for the variable-step L2 type formula of Caputo’s derivative. Numerical examples with an adaptive stepping procedure are provided to demonstrate the accuracy and the effectiveness of our proposed method.

Spatiotemporal Patterns in a Space–Time Discrete SIRS Epidemic Model with Self- and Cross-Diffusion

In this work, a two-dimensional Coupled Map Lattice (CML) method is applied to a continuous Reaction–Diffusion (RD) SIRS epidemic model. The model involves a complex, nonlinear incidence rate and the spatially heterogeneous self- and cross-diffusion coefficients. The CML method applied here results in a space–time discrete counterpart of the model. The basic reproduction number is derived, and the local stability of the fixed points of the discrete model is established. The two typical bifurcations of codimension 1, such as the flip and the Neimark–Sacker (NS) bifurcations, are investigated by utilizing the center manifold theory and the normal form theorem. Furthermore, the codimension-2 fold–flip bifurcation about the endemic fixed point is discussed in detail. Turing instability criterion of the cross-diffusion epidemic model is established by utilizing the CML method under periodic boundary conditions. Finally, the numerical simulations clarify both the effectiveness of all theoretical findings and the effects of time step size and cross-diffusion on the spatiotemporal patterns. The main contribution of this work is the discovery of some interesting Turing and non-Turing patterns which may be induced by flip–Turing instability, NS–Turing instability or fold–flip–Turing instability. Concretely speaking, spiral patterns are observed near the flip bifurcation threshold value and ultimately evolve into an irregular defect pattern. Stripe pattern, labyrinth pattern and circular pattern are observed near the NS bifurcation threshold value. The disorderly interwoven ribbon pattern, mosaic pattern and some irregular oscillatory patterns are observed near the fold–flip bifurcation point. Our results greatly help in predicting the long-term behavior of the space–time discrete SIRS epidemic model under study.

Analysis of the dynamic response for Kirchhoff plates by the element-free Galerkin method

Dynamic response analysis involves examining structural behavior under dynamic loading conditions. Transient dynamic analysis emerges as a method employed to assess the responses of deformable bodies. Due to the depth analysis of transient responses for thin plates, the associated error analysis work becomes particularly important. In this work, we conduct stability and convergence analysis of the element-free Galerkin (EFG) method for the dynamic response of Kirchhoff plate model. We start by discretizing the thin plate dynamics using the EFG method for the spatial domain and the central difference scheme for the temporal domain. Stability and convergence analysis for transient responses, including deflection and moment responses, are derived similarly to those in two-dimensional transient structural dynamic analysis. The key to the analysis is transforming the governing equations and clamped boundary conditions into a weak formulation with penalty terms using the penalty method. The main results demonstrate a significant correlation between the error estimate and both nodal spacing and time step size. Additionally, the continuity of the approximation functions and the selection of penalty factors significantly affect numerical accuracy. To validate the theoretical results, we conduct numerical experiments involving clamped square and L-shaped plates with uniform and nonuniform node distributions, as well as a perforated plate model. Furthermore, we investigate the impact of node influence domains and penalty factors on numerical accuracy, facilitating parameter selection for practical engineering applications.

DynGMA: A robust approach for learning stochastic differential equations from data

Learning unknown stochastic differential equations (SDEs) from observed data is a significant and challenging task with applications in various fields. Current approaches often use neural networks to represent drift and diffusion functions, and construct likelihood-based loss by approximating the transition density to train these networks. However, these methods often rely on one-step stochastic numerical schemes, necessitating data with sufficiently high time resolution. In this paper, we introduce novel approximations to the transition density of the parameterized SDE: a Gaussian density approximation inspired by the random perturbation theory of dynamical systems, and its extension, the dynamical Gaussian mixture approximation (DynGMA). Benefiting from the robust density approximation, our method exhibits superior accuracy compared to baseline methods in learning the fully unknown drift and diffusion functions and computing the invariant distribution from trajectory data. And it is capable of handling trajectory data with low time resolution and variable, even uncontrollable, time step sizes, such as data generated from Gillespie's stochastic simulations. We then conduct several experiments across various scenarios to verify the advantages and robustness of the proposed method.

An Integral-like Numerical Approach for Solving Burgers’ Equation

The Burgers’ equation, commonly appeared in the study of turbulence, fluid dynamics, shock waves, and continuum mechanics, is a crucial part of the dynamical core of any numerical weather model, influencing simulated meteorological phenomena. While past studies have suggested several robust numerical approaches for solving the equation, many are too complicated for practical adaptation and too computationally expensive for operational deployment. This paper introduces an unconventional approach based on spline polynomial interpolations and the Hopf-Cole transformation. Using Taylor expansion to approximate the exponential term in the Hopf-Cole transformation, the analytical solution of the simplified equation is discretized to form our proposed scheme. The scheme is explicit and adaptable for parallel computing, although certain types of boundary conditions need to be employed implicitly. Three distinct test cases were utilized to evaluate its accuracy, parallel scalability, and numerical stability. In the aspect of accuracy, the schemes employed cubic and quintic spline interpolation perform equally well, managing to reduce the <i>ӏ</i><sub>1</sub>, <i>ӏ</i><sub>2</sub>, and <i>ӏ</i><sub>∞</sub> error norms down to the order of 10<sup>−4</sup>. Parallel scalability observed in the weak-scaling experiment depends on time step size but is generally as good as any explicit scheme. The stability condition is <i>ν</i>∆<i>t</i>/∆<i>x</i><sup>2</sup> > 0.02, including the viscosity coefficient <i>ν</i> due to the Hopf-Cole transformation step. From the stability condition, the schemes can run at a large time step size ∆<i>t</i> even when using a small grid spacing ∆<i>x</i>, emphasizing its suitability for practical applications such as numerical weather prediction.

Precise integration boundary element method for solving nonlinear transient heat conduction problems with temperature dependent thermal conductivity and heat capacity

The precise integration boundary element method (PIBEM) is used to solve nonlinear transient heat conduction problems with temperature dependent thermal conductivity and heat capacity in this article. At first, a boundary-domain integral equation can be obtained by using the fundamental solution for steady linear heat conduction problems. Secondly, the domain integrals are converted into boundary integrals by using the radial integration method to get a pure boundary integral equation, which can retain the advantage of dimension reduction of the boundary element method. Then, after discretization, a first-order nonlinear differential equation system in time is obtained, and the precise integration method with predictor-corrector technique is used to solve the system of nonlinear equations. In the end, several numerical examples are presented to verify the accuracy and stability of the presented method. The results obtained by the presented method are less affected by the time step size and satisfactory results can be obtained even for a large time step size.

An optimized time-adaptive aerothermal coupling calculation method for aerothermal analysis of disc cavity system

In order to further achieve the balance between the calculation accuracy and efficiency of the transient analysis of the aero-engine disc cavity system, an Optimized Time-adaptive Aerothermal Coupling calculation (OTAC) method has been proposed. It combines one-dimensional transient calculation of air system, Conventional Sequence Staggered (CSS) method, Time-adaptive Aerothermal Coupling calculation (TAC) method and differential evolution optimization algorithm to obtain an efficient and high-precision aerothermal coupling calculation method of air system. Considering both the heat conduction in the solid domain and the flow in the fluid domain as unsteady states in the OTAC, the interaction of fluid–solid information within a single coupling time step size was implemented based on the CSS method. Furthermore, the coupling time step size was automatically adjusted with the number of iterations by using the Proportional-Integral-Derivative (PID) controller. Results show that when compared with the traditional loosely coupling method with a fixed time step size, the computational accuracy and efficiency of the OTAC method are improved by 8.9% and 30%, respectively. Compared with the tight coupling calculation, the OTAC method can achieve a speedup of 1 to 2 orders of magnitude, while the calculation error is maintained within 6.1%.

Numerical Calculations of Adiabatic Invariants From MHD‐Driven Magnetic Fields

AbstractThe adiabatic invariants (M, J, Φ) and the related invariants (M, K, L*) have been established as effective coordinate systems for describing radiation belt dynamics at a theoretical level, and through numerical techniques, can be paired with in situ observations to order phase‐space density. To date, methods for numerical techniques to calculate adiabatic invariants have focused on empirical models such the Tsyganenko models TS05, T96, and T89. In this work, we develop methods based on numerical integration and variable step size iteration for the calculation of adiabatic invariants, applying the method to the Lyon‐Fedder‐Mobarry (LFM) global magnetohydrodynamics (MHD) simulation code, with optional coupling to the Rice Convection Model (RCM). By opening the door to adiabatic invariant modeling with MHD magnetic fields, the opportunity for exploratory modeling work of radiation belt dynamics is enabled. Calculations performed using LFM are cross‐referenced with the same code applied to the T96 and TS05 Tsyganenko models evaluated on the LFM grid. Important aspects of the adiabatic invariant calculation are reviewed and discussed, including (a) sensitivity to magnetic field model used, (b) differences in the problem between quiet and disturbed geomagnetic states, and (c) the selection of key parameters, such as the magnetic local time step size for drift shell determination. The rigorous development and documentation of this algorithm additionally acts as preliminary step for future thorough reassessment of in situ phase‐space density results using alternative magnetic field models.