Maximum Stable Time Step Research Articles

A time splitting scheme for the nonhydrostatic atmospheric model

A semi-implicit finite difference scheme with time splitting corresponding to the multi-scale nature of atmospheric dynamics is developed for the nonhydrostatic atmospheric model. The fast acoustic and gravity waves are approximated implicitly, while relatively slow advective terms and Rossby modes are treated explicitly. The used level of implicitness requires a solution of the three-dimensional linear elliptic equation at each time step, which is accomplished efficiently by vertical decoupling and subsequent application of fast and accurate multigrid solvers for the most significant vertical modes, and fast and simpler solvers for the secondary ones. Analysis of the linear stability of the scheme reveals that the Courant–Friedrichs–Lewy condition on the time step associated with gravity and acoustic waves can be circumvented and the maximum stable time step can be brought into closer agreement with the time step limitations arising from accuracy considerations. The performed numerical experiments show the computational efficiency of the designed scheme and accuracy of the predicted atmospheric fields.

Numerical methods for nonlinear hillslope transport laws

The numerical methods used to solve nonlinear sediment transport equations often set very restrictive limits on the stability and accuracy of landscape evolution models. This is especially true for hillslope transport laws in which sediment flux increases nonlinearly as the surface slope approaches a limiting value. Explicit‐time finite difference methods applied to such laws are subject to fundamental limits on numerical stability that require time steps much shorter than the timescales over which landscapes evolve, creating a heavy computational burden. I present an implicit method for nonlinear hillslope transport that builds on a previously proposed approach to modeling alluvial sediment transport and improves stability and accuracy by avoiding the direct calculation of sediment flux. This method can be adapted to any transport law in which the expression for sediment flux is differentiable. Comparisons of numerical solutions with analytic solutions in one and two dimensions show that the implicit method retains the accuracy of a standard explicit method while permitting time steps several orders of magnitude longer than the maximum stable time step for the explicit method. The ability to take long time steps affords a substantial savings in overall computation time, despite the implicit method's higher per‐iteration computational cost. Implicit models for hillslope evolution also offer a distinct advantage when modeling the response of hillslopes to incising channels.

Implementation of an implicit shear Alfvén operator in the M3D code

A new pair of linear operators is introduced into the partially implicit but largely explicit time advance algorithm for the MHD equations in the M3D code. These operators jointly render the algorithm implicit with respect to propagation of the shear Alfvén wave to leading order. As a result, the maximum stable time step is increased by up to a factor of three for many practical applications. The additional computational cost is trivial for the linear mode of operation and adds only approximately 30% additional execution time per step for nonlinear calculations.

Application of Dey–Mittra conformal boundary algorithm to 3D electromagnetic modeling

The Dey–Mittra conformal boundary conditions have been implemented for the finite-difference time-domain (FDTD) electromagnetic solver of the VORPAL plasma simulation framework and studied in the context of three-dimensional, large-scale computations. The maximum stable time step when using these boundary conditions can be arbitrarily small, due to the presence of small fractional cells inside the vacuum region. Use of the Gershgorin Circle theorem allows the determination of a rigorous criterion for exclusion of small cells in order to have numerical stability for particular values of the ratio f DM ≡ Δ t / Δ t CFL of the time step to the Courant–Friedrichs–Lewy value for the infinite system. Application to a spherical cavity shows that these boundary conditions allow computation of frequencies with second-order error for sufficiently small f DM . However, for sufficiently fine resolution, dependent on f DM , the error becomes first order, just like the error for stair-step boundary conditions. Nevertheless, provided one does use a sufficiently small value of f DM , one can obtain third-order accuracy through Richardson extrapolation. Computations for the TESLA superconducting RF cavity design compare favorably with experimental measurements.

Ridging, strength, and stability in high‐resolution sea ice models

In multicategory sea ice models the compressive strength of the ice pack is often assumed to be a function of the potential energy of pressure ridges. This assumption, combined with other standard features of ridging schemes, allows the ice strength to change dramatically on short timescales. In high‐resolution (∼10 km) sea ice models with a typical time step (∼1 hour), abrupt strength changes can lead to large internal stress gradients that destabilize the flow. The unstable flow is characterized by large oscillations in ice concentration, thickness, strength, velocity, and strain rates. Straightforward, physically motivated changes in the ridging scheme can reduce the likelihood of abrupt strength changes and improve stability. In simple test problems with flow toward and around topography, stability is significantly enhanced by eliminating the threshold fraction G* in the ridging participation function. Use of an exponential participation function increases the maximum stable time step at 10‐km resolution from less than 30 min to about 2 hours. Modifying the redistribution function to build thinner ridges modestly improves stability and also gives better agreement between modeled and observed thickness distributions. Allowing the ice strength to increase linearly with the mean ice thickness improves stability but probably underestimates the maximum stresses.

Conformal FDTD-methods to avoid time step reduction with and without cell enlargement

During the last decades there have been considerable efforts to develop accurate and yet simple conformal methods for modelling curved boundaries within the finite difference time domain (FDTD) algorithm. In an earlier publication we proposed the uniformly stable conformal (USC) approach as a general three-dimensional extension of FDTD without the need to reduce the maximum stable time step. The main idea of USC is the usage of virtually enlarged cells near to the boundary, leading to an increased implementation effort. In this paper we review the USC method and introduce a new simple and accurate conformal scheme which does not use such enlarged cells. This simplified conformal (SC) scheme has the same number of operations and algorithmic logic as the standard “staircase” method, and thus is easily realizable in existing FDTD codes. Like USC, it leads to accurate results without time step reduction, showing a nearly second order convergence in practice. The method is verified and compared to other approaches by means of several numerical 2D and 3D examples.

Density-contrast instabilities in finite element modelling of slow viscous flow

Finite element methods are often used to model Earth processes involving slow viscous or viscoelastic flow. Inertial terms of the Navier-Stokes equations are neglected in very slow flows, so timestep size is not limited by the Courant instability. However, where there is advection of density contrasts in a gravitational field, over-advection can lead to numerically induced flow oscillations. We derive analytic results for the maximum stable timestep size in two cases: a free surface over a fluid of uniform density, and a free surface kept level by sedimentation/erosion, but with a density gradient in the underlying medium. Using parameters appropriate to the Earth's crust we show that the density-contrast instability occurs for timesteps larger than 3000 years for the constant-density case. For a fluid with a density gradient of 10 kg/m3 per km the solution is stable for timesteps up to about 200,000 years if full erosion/sedimentation is implemented.

Transient Modal Analysis of Quasi-Implicit FDTD Schemes

Time domain simulations for high-frequency applications are widely dominated by the leapfrog time-integration scheme, especially in combination with finite-difference methods (FDTD) or the finite integration technique (FIT). As an explicit method, however, the leapfrog algorithm is restricted to a maximum stable time step, and recently some alternative, unconditionally stable techniques have been proposed to overcome this limitation. We analyze such schemes using a transient modal decomposition of the electric fields. It is shown that stability alone is not sufficient to guarantee correct results, but additionally important conservation properties have to be met.

Finite‐difference modeling of SH‐wave propagation in nonwelded contact media

Many geological structures of interest are known to exhibit fracturing. Fracturing directly affects seismic wave propagation because, depending on its scale, fracturing may give rise to scattering and/or anisotropy. A fracture may be described mathematically as an interface in nonwelded contact (i.e., as a displacement discontinuity). This poses a difficulty for finite‐difference modeling of seismic wave propagation in fractured media, because the standard heterogeneous approach assumes welded contact. In the past, this difficulty has been circumvented by incorporating nonwelded contact into the medium parameters using equivalent medium theory. We present an alternate method based on the homogeneous approach to finite differencing, whereby nonwelded contact boundary conditions are imposed explicitly. For simplicity, we develop the method in the SH‐wave case.In the homogeneous approach, nonwelded contact boundary conditions are discretized by introducing auxiliary, so‐called fictitious, grid points. Wavefield values at fictitious grid points are then used in the discrete equation of motion, so that the time‐evolved wavefield satisfies the correct boundary conditions. Although not as general as the heterogeneous approach, the homogeneous approach has the advantage of being relatively simple and manifestly satisfying nonwelded contact boundary conditions. For fractures aligned with the numerical grid, the homogeneous and heterogeneous approaches yield identical results. In particular, in both approaches nonwelded contact results in a larger maximum stable time step size than in the welded contact case.

Implicit solution of uncertain volatility/transaction cost option pricing models with discretely observed barriers

Option pricing models with uncertain volatility/transaction costs give rise to a nonlinear PDE. Previous work has focused on explicit methods. However, pricing discretely observed barrier options requires a very small grid spacing near the barrier, and as a result, the maximum stable timestep for an explicit method is impractically small. A fully implicit method is developed for nonlinear option pricing models, and applied to arithmetic step options, where the option loses a fraction of its value for every day over the barrier.

Instability in Runge-Kutta schemes for simulation of oil recovery

Runge-Kutta timestep schemes have been used recently in an attempt to increase the maximum stable timestep for IMPES (implicit pressure, explicit saturation) simulation of oil recovery. It has been claimed that anm-stage method will increase the real stability boundary bym2. A rigorous stability analysis for a Buckley-Leverett problem shows that this claim is false, and that previous stabilized IMPES schemes are never more efficient than ordinary IMPES. Test calculations support these results.