Severe Stability Constraints Research Articles

Stabilized leapfrog based local time-stepping method for the wave equation

Local time-stepping methods permit to overcome the severe stability constraint on explicit methods caused by local mesh refinement without sacrificing explicitness. Diaz and Grote [SIAM J. Sci. Comput. 31 (2009), pp. 1985–2014] proposed a leapfrog based explicit local time-stepping (LF-LTS) method for the time integration of second-order wave equations. Recently, optimal convergence rates were proved for a conforming FEM discretization, albeit under a CFL stability condition where the global time-step, Δt, depends on the smallest elements in the mesh (see M. J. Grote, M. Mehlin, and S. A. Sauter [SIAM J. Numer. Anal. 56 (2018), pp. 994–1021]). In general one cannot improve upon that stability constraint, as the LF-LTS method may become unstable at certain discrete values of Δt. To remove those critical values of Δt, we apply a slight modification (as in recent work on LF-Chebyshev methods by Carle, Hochbruck, and Sturm [SIAM J. Numer. Anal. 58 (2020), pp. 2404–2433]) to the original LF-LTS method which nonetheless preserves its desirable properties: it is fully explicit, second-order accurate, satisfies a three-term (leapfrog like) recurrence relation, and conserves the energy. The new stabilized LF-LTS method also yields optimal convergence rates for a standard conforming FE discretization, yet under a CFL condition where Δt no longer depends on the mesh size inside the locally refined region.

High order integration factor methods for systems with inhomogeneous boundary conditions

Abstract Due to the high order spatial derivatives and stiff reactions, severe temporal stability constraints on the time step are generally required when developing numerical methods for solving high order partial differential equations. Implicit integration method (IIF) method along with its compact form (cIIF), which treats spatial derivatives exactly and reaction terms implicitly, provides excellent stability properties with good efficiency by decoupling the treatment of reaction and spatial derivatives. One major challenge for IIF is storage and calculation of the potential dense exponential matrices of the sparse discretization matrices resulted from the linear differential operators. The compact representation for IIF (cIIF) was introduced to save the computational cost and storage for this purpose. Another challenge is finding the matrix of high order space discretization, especially near the boundaries. In this paper, we extend IIF method to high order discretization for spatial derivatives through an example of reaction diffusion equation with fourth order accuracy, while the computational cost and storage are similar to the general second order cIIF method. The method can also be efficiently applied to deal with other types of partial differential equations with both homogeneous and inhomogeneous boundary conditions. Direct numerical simulations demonstrate the efficiency and accuracy of the approach.

Compact implicit integration factor methods for a family of semilinear fourth-order parabolic equations

When developing efficient numerical methods for solving parabolictypes of equations, severe temporal stability constraints on thetime step are often required due to the high-order spatialderivatives and/or stiff reactions. The implicit integration factor(IIF) method, which treats spatial derivative terms explicitly andreaction terms implicitly, can provide excellent stabilityproperties in time with nice accuracy. One major challenge for the IIFis the storage and calculation of the dense exponentials of the sparse discretization matrices resulted from thelinear differential operators. The compact representation of the IIF(cIIF) can overcome this shortcoming and greatly savecomputational cost and storage. On the other hand, the cIIF is oftenhard to be directly applied to deal with problems involving cross derivatives. In this paper, bytreating the discretization matrices in diagonalized forms, wedevelop an efficient cIIF method for solving a family of semilinear fourth-orderparabolic equations, in which the bi-Laplace operator is explicitlyhandled and the computational cost and storage remain the same as to theclassic cIIF for second-order problems. In particular, the proposed method can deal with not onlystiff nonlinear reaction terms but also various types ofhomogeneous or inhomogeneous boundary conditions. Numericalexperiments are finally presented to demonstrate effectiveness and accuracyof the proposed method.

Array-representation integration factor method for high-dimensional systems

High order spatial derivatives and stiff reactions often introduce severe temporal stability constraints on the time step in numerical methods. Implicit integration method (IIF) method, which treats diffusion exactly and reaction implicitly, provides excellent stability properties with good efficiency by decoupling the treatment of reactions and diffusions. One major challenge for IIF is storage and calculation of the potential dense exponential matrices of the sparse discretization matrices resulted from the linear differential operators. Motivated by a compact representation for IIF (cIIF) for Laplacian operators in two and three dimensions, we introduce an array-representation technique for efficient handling of exponential matrices from a general linear differential operator that may include cross-derivatives and non-constant diffusion coefficients. In this approach, exponentials are only needed for matrices of small size that depend only on the order of derivatives and number of discretization points, independent of the size of spatial dimensions. This method is particularly advantageous for high-dimensional systems, and it can be easily incorporated with IIF to preserve the excellent stability of IIF. Implementation and direct simulations of the array-representation compact IIF (AcIIF) on systems, such as Fokker–Planck equations in three and four dimensions and chemical master equations, in addition to reaction–diffusion equations, show efficiency, accuracy, and robustness of the new method. Such array-presentation based on methods may have broad applications for simulating other complex systems involving high-dimensional data.

A small-scale decomposition for 3D boundary integral computations with surface tension

An efficient, non-stiff boundary integral method for the initial value problem for interfacial Darcy flow (which is a model of porous media flow) in three space dimensions is presented. We consider a ‘doubly-periodic’ interface separating two fluids, with surface tension present at the boundary. Surface tension introduces high order (i.e., high derivative) terms in the governing equation, and this imposes a severe stability constraint on explicit time-integration methods. Furthermore, the high order terms appear in a nonlocal operator, which makes it difficult to design an efficient implicit method. The stiffness is removed by developing a small-scale decomposition in the spirit of prior work in the two-dimensional problem by Hou, Lowengrub, and Shelley. In order to develop this small-scale decomposition, we formulate the problem using a generalized isothermal parameterization of the free surface. An additional difficulty is the efficient calculation of the Birkhoff–Rott integral for the velocity of the interface. We present a new algorithm, based on Ewald summation, to compute this in O(NlogN) operations, where N is the number of interface grid points. Our non-stiff method is expected to apply widely to problems for doubly-periodic interfacial flow with surface tension which have a boundary integral formulation.

Locally Implicit Time Integration Strategies in a Discontinuous Galerkin Method for Maxwell’s Equations

An attractive feature of discontinuous Galerkin (DG) spatial discretization is the possibility of using locally refined space grids to handle geometrical details. However, locally refined meshes lead to severe stability constraints on explicit integration methods to numerically solve a time-dependent partial differential equation. If the region of refinement is small relative to the computational domain, the time step size restriction can be overcome by blending an implicit and an explicit scheme where only the solution variables living at fine elements are treated implicitly. The downside of this approach is having to solve a linear system per time step. But due to the assumed small region of refinement relative to the computational domain, the overhead will also be small while the solution can be advanced in time with step sizes determined by the coarse elements. In this paper, we present two locally implicit time integration methods for solving the time-domain Maxwell equations spatially discretized with a DG method. Numerical experiments for two-dimensional problems illustrate the theory and the usefulness of the implicit–explicit approaches in presence of local refinements.

High-order explicit local time-stepping methods for damped wave equations

Locally refined meshes impose severe stability constraints on explicit time-stepping methods for the numerical simulation of time dependent wave phenomena. Local time-stepping methods overcome that bottleneck by using smaller time-steps precisely where the smallest elements in the mesh are located. Starting from classical Adams–Bashforth multi-step methods, local time-stepping methods of arbitrarily high order of accuracy are derived for damped wave equations. When combined with a finite element discretization in space with an essentially diagonal mass matrix, the resulting time-marching schemes are fully explicit and thus inherently parallel. Numerical experiments with continuous and discontinuous Galerkin finite element discretizations corroborate the expected rates of convergence and illustrate the usefulness of these local time-stepping methods.

On the acceleration of explicit finite difference methods for option pricing

Implicit finite difference methods are conventionally preferred over their explicit counterparts for the numerical valuation of options. In large part the reason for this is a severe stability constraint known as the Courant–Friedrichs–Lewy (CFL) condition which limits the latter class's efficiency. Implicit methods, however, are difficult to implement for all but the most simple of pricing models, whereas explicit techniques are easily adapted to complex problems. For the first time in a financial context, we present an acceleration technique, applicable to explicit finite difference schemes describing diffusive processes with symmetric evolution operators, called Super-Time-Stepping. We show that this method can be implemented as part of a more general approach for non-symmetric operators. Formal stability is thereby deduced for the exemplar cases of European and American put options priced under the Black–Scholes equation. Furthermore, we introduce a novel approach to describing the efficiencies of finite difference schemes as semi-empirical power laws relating the minimal real time required to carry out the numerical integration to a solution with a specified accuracy. Tests are described in which the method is shown to significantly ameliorate the severity of the CFL constraint whilst retaining the simplicity of the underlying explicit method. Degrees of acceleration are achieved yielding comparable, or superior, efficiencies to a set of benchmark implicit schemes. We infer that the described method is a powerful tool, the explicit nature of which makes it ideally suited to the treatment of symmetric and non-symmetric diffusion operators describing complex financial instruments including multi-dimensional systems requiring representation on decomposed and/or adaptive meshes.

All Speed Scheme for the Low Mach Number Limit of the Isentropic Euler Equations

AbstractAn all speed scheme for the Isentropic Euler equations is presented in this paper. When the Mach number tends to zero, the compressible Euler equations converge to their incompressible counterpart, in which the density becomes a constant. Increasing approximation errors and severe stability constraints are the main difficulty in the low Mach regime. The key idea of our all speed scheme is the special semi-implicit time discretization, in which the low Mach number stiff term is divided into two parts, one being treated explicitly and the other one implicitly. Moreover, the flux of the density equation is also treated implicitly and an elliptic type equation is derived to obtain the density. In this way, the correct limit can be captured without requesting the mesh size and time step to be smaller than the Mach number. Compared with previous semi-implicit methods, firstly, nonphysical oscillations can be suppressed by choosing proper parameter, besides, only a linear elliptic equation needs to be solved implicitly which reduces much computational cost. We develop this semi-implicit time discretization in the framework of a first order Local Lax-Friedrichs (or Rusanov) scheme and numerical tests are displayed to demonstrate its performances.

A non-stiff boundary integral method for 3D porous media flow with surface tension

We present an efficient, non-stiff boundary integral method for 3D porous media flow with surface tension. Surface tension introduces high order (i.e., high derivative) terms into the evolution equations, and this leads to severe stability constraints for explicit time-integration methods. Furthermore, the high order terms appear in non-local operators, making the application of implicit methods difficult. Our method uses the fundamental coefficients of the surface as dynamical variables, and employs a special isothermal parameterization of the interface which enables efficient application of implicit or linear propagator time-integration methods via a small-scale decomposition. The method is tested by computing the relaxation of an interface to a flat surface under the action of surface tension. These calculations employ an approximate interface velocity to test the stiffness reduction of the method. The approximate velocity has the same mathematical form as the exact velocity, but avoids the numerically intensive computation of the full Birkhoff–Rott integral. The algorithm is found to be effective at eliminating the severe time-step constraint that plagues explicit time-integration methods.

Dynamic simulation of locally inextensible vesicles suspended in an arbitrary two-dimensional domain, a boundary integral method

We consider numerical algorithms for the simulation of hydrodynamics of two-dimensional vesicles suspended in a viscous Stokesian fluid. The motion of vesicles is governed by the interplay between hydrodynamic and elastic forces. Continuum models of vesicles use a two-phase fluid system with interfacial forces that include tension (to maintain local “surface” inextensibility) and bending. Vesicle flows are challenging to simulate. On the one hand, explicit time-stepping schemes suffer from a severe stability constraint due to the stiffness related to high-order spatial derivatives in the bending term. On the other hand, implicit time-stepping schemes can be expensive because they require the solution of a set of nonlinear equations at each time step. Our method is an extension of the work of Veerapaneni et al. [S.K. Veerapaneni, D. Gueyffier, D. Zorin, G. Biros, A boundary integral method for simulating the dynamics of inextensible vesicles suspended in a viscous fluid in 2D, Journal of Computational Physics 228(7) (2009) 2334–2353], in which a semi-implicit time-marching scheme based on a boundary integral formulation of the Stokes problem for vesicles in an unbounded medium was proposed. In this paper, we consider two important generalizations: (i) confined flows within arbitrary-shaped stationary/moving geometries; and (ii) flows in which the interior (to the vesicle) and exterior fluids have different viscosity. In the rest of the paper, we will refer to this as the “ viscosity contrast”. These two problems require solving additional integral equations and cause nontrivial modifications to the previous numerical scheme. Our method does not have severe time-step stability constraints and its computational cost-per-time-step is comparable to that of an explicit scheme. The discretization is pseudo-spectral in space, and multistep BDF in time. We conduct numerical experiments to investigate the stability, accuracy and the computational cost of the algorithm. Overall, our method achieves several orders of magnitude speed-up compared to standard explicit schemes. As a preliminary validation of our scheme, we study the dependence of the inclination angle of a single vesicle in shear flow on the viscosity contrast and the reduced area of the vesicle, the lateral migration of vesicles in shear flow, the dispersion of two vesicles, and the effective viscosity of a dilute suspension of vesicles.

A boundary integral method for simulating the dynamics of inextensible vesicles suspended in a viscous fluid in 2D

We present a new method for the evolution of inextensible vesicles immersed in a Stokesian fluid. We use a boundary integral formulation for the fluid that results in a set of nonlinear integro-differential equations for the vesicle dynamics. The motion of the vesicles is determined by balancing the non-local hydrodynamic forces with the elastic forces due to bending and tension. Numerical simulations of such vesicle motions are quite challenging. On one hand, explicit time-stepping schemes suffer from a severe stability constraint due to the stiffness related to high-order spatial derivatives and a milder constraint due to a transport-like stability condition. On the other hand, an implicit scheme can be expensive because it requires the solution of a set of nonlinear equations at each time step. We present two semi-implicit schemes that circumvent the severe stability constraints on the time step and whose computational cost per time step is comparable to that of an explicit scheme. We discretize the equations by using a spectral method in space, and a multistep third-order accurate scheme in time. We use the fast multipole method (FMM) to efficiently compute vesicle–vesicle interaction forces in a suspension with a large number of vesicles. We report results from numerical experiments that demonstrate the convergence and algorithmic complexity properties of our scheme.

On the Acceleration of Explicit Finite Difference Methods for Option Pricing

Implicit finite difference methods are conventionally preferred over their explicit counterparts for the valuation of options. In large part the reason for this is a severe stability constraint known as the Courant-Friedrichs-Lewy (CFL) condition which limits the latters' efficiencies. Implicit methods, however, are difficult to implement for all but the most simple of pricing models whereas explicit techniques are easily adapted to complex problems. In this work we present an acceleration technique for explicit finite difference schemes called Super-Time-Stepping (STS) for the first time in a financial context. Furthermore, we introduce a novel method for describing the efficiencies of finite difference schemes as semi-empirical power laws relating the minimal walltime W required to attain a solution with an error of magnitude E. For European and American put option test cases we demonstrate degrees of acceleration over standard explicit methods resulting in efficiencies comparable, or superior, to a set of implicit scheme benchmarks. We conclude that STS is a powerful tool for the numerical pricing of options and propose it as the method-of-choice for exotic financial intruments such as those requiring multi-dimensional descriptions on adaptive meshes.

Energy Conserving Explicit Local Time Stepping for Second-Order Wave Equations

Locally refined meshes impose severe stability constraints on explicit time-stepping methods for the numerical simulation of time dependent wave phenomena. To overcome that stability restriction, local time-stepping methods are developed, which allow arbitrarily small time steps precisely where small elements in the mesh are located. When combined with a symmetric finite element discretization in space with an essentially diagonal mass matrix, the resulting discrete numerical scheme is explicit, is inherently parallel, and exactly conserves a discrete energy. Starting from the standard second-order “leap-frog” scheme, time-stepping methods of arbitrary order of accuracy are derived. Numerical experiments illustrate the efficiency and usefulness of these methods and validate the theory.

Pseudo-time stepping methods for space–time discontinuous Galerkin discretizations of the compressible Navier–Stokes equations

The space–time discontinuous Galerkin discretization of the compressible Navier–Stokes equations results in a non-linear system of algebraic equations, which we solve with pseudo-time stepping methods. We show that explicit Runge–Kutta methods developed for the Euler equations suffer from a severe stability constraint linked to the viscous part of the equations and propose an alternative to relieve this constraint while preserving locality. To evaluate its effectiveness, we compare with an implicit–explicit Runge–Kutta method which does not suffer from the viscous stability constraint. We analyze the stability of the methods and illustrate their performance by computing the flow around a 2D airfoil and a 3D delta wing at low and moderate Reynolds numbers.

Simulating the dynamics and interactions of flexible fibers in Stokes flows

The dynamics of slender filaments or fibers suspended in Stokesian fluids are fundamental to understanding many flows arising in physics, biology and engineering. Such filaments can have aspect ratios of length to radius ranging from a few tens to several thousands. Full discretizations of such 3D flows are very costly. Instead, we employ a non-local slender body theory that yields an integral equation, along the filament centerline, relating the force exerted on the body to the filament velocity. This hydrodynamical description takes into account the effect of the filament on the fluid, and is extended to capture the interaction of multiple filaments as mediated by the intervening fluid. We consider filaments that are inextensible and elastic. Replacing the force in the slender body integral equation by an explicit expression that uses Euler–Bernoulli theory to model bending and tensile forces yields an integral expression for the velocity of the filament centerlines, coupled to auxiliary integro-differential equations for the filament tensions. Based on a regularized version of these slender body equations that is asymptotically equivalent to the original formulation, we construct a numerical method which uses a combination of finite differences, implicit time-stepping to avoid severe stability constraints, and special quadrature methods for nearly singular integrals. We present simulations of single flexible filaments, as well as multiple interacting filaments, evolving in a background shear flow. These simulations show shear induced buckling and relaxation of the filaments, leading to the storage and release of elastic energy. These dynamics are responsible for the development of positive first normal stress differences, commonly associated with visco-elastic fluids that are suspensions of microscopic elastic fibers.

The Nonlinear Evolution of Vortex Sheets with Surface Tension in Axisymmetric Flows

The presence of surface tension for interfacial flows usually leads to severe stability constraints for explicit time integration methods. Moreover, the nonlocality and nonlinearity of the high-order terms make the application of implicit methods difficult. In this paper, a computational strategy is presented for computing the motion of fluid interfaces with surface tension in axisymmetric flows using boundary integral techniques. This method is based on adaptive quadratures for the principal-value integrals and a small-scale decomposition for the treatment of surface tension through a vector-potential formulation. A study of the method is conducted in the context of vortex sheet evolution with surface tension in axisymmetric flows. The method is found to be accurate, efficient, and robust. Numerical simulations indicate that the dynamics of vortex sheets with surface tension frequently result in topological singularities (i.e., self-intersection). Away from the axis of symmetry, these singularities are similar to those found in the two-dimensional flows. Singularities occurring near the axis of symmetry take a different form.

Numerical Solutions to Free Boundary Problems

Many physically interesting problems involve propagation of free surfaces. Vortex-sheet roll-up in hydrodynamic instability, wave interactions on the ocean's free surface, the solidification problem for crystal growth and Hele-Shaw cells for pattern formation are some of the significant examples. These problems present a great challenge to physicists and applied mathematicians because the underlying problem is very singular. The physical solution is sensitive to small perturbations. Naïve discretisations may lead to numerical instabilities. Other numerical difficulties include singularity formation and possible change of topology in the moving free surfaces, and the severe time-stepping stability constraint due to the stiffness of high-order regularisation effects, such as surface tension.This paper reviews some of the recent advances in developing stable and efficient numerical algorithms for solving free boundary-value problems arising from fluid dynamics and materials science. In particular, we will consider boundary integral methods and the level-set approach for water waves, general multi-fluid interfaces, Hele–Shaw cells, crystal growth and solidification. We will also consider the stabilising effect of surface tension and curvature regularisation. The issue of numerical stability and convergence will be discussed, and the related theoretical results for the continuum equations will be addressed. This paper is not intended to be a detailed survey and the discussion is limited by both the taste and expertise of the author.

Removing the stiffness from interfacial flows with surface tension

A new formulation and new methods are presented for computing the motion of fluid interfaces with surface tension in two-dimensional, irrotational, and incompressible fluids. Through the Laplace-Young condition at the interface, surface tension introduces high-order terms, both nonlinear and nonlocal, into the dynamics. This leads to severe stability constraints for explicit time integration methods and makes the application of implicit methods difficult. This new formulation has all the nice properties for time integration methods that are associated with having a linear, constant coefficient, highest order term. That is, using this formulation, we give implicit time integration methods that have no high order time step stability constraint associated with surface tension and are explicit in Fourier space. The approach is based on a boundary integral formulation and applies more generally, even to problems beyond the fluid mechanical context. Here they are applied to computing with high resolution the motion of interfaces in Hele-Shaw flows and the motion of free surfaces in inviscid flows governed by the Euler equations. One Hele-Shaw computation shows the behavior of an expanding gas bubble over long-time as the interface undergoes successive tip-splittings and finger competition. A second computation shows the formation of a very ramified interface through the interaction of surface tension with an unstable density stratification. In Euler flows, the computation of a vortex sheet shows its roll-up through the Kelvin-Helmholtz instability. This motion culminates in the late time self-intersection of the interface, creating trapped bubbles of fluid. This is, we believe, a type of singularity formation previously unobserved for such flows in 2D. Finally, computations of falling plumes in an unstably stratified Boussinesq fluid show a very similar behavior.

Advection‐diffusion modelling using the modified QUICK scheme

AbstractIn recent years the QUICK finite difference scheme has been increasingly used in solving the advection‐diffusion equation, particularly for water quality modelling studies relating to coastal and estuarine flows. This scheme has the benefits of mass conservation, reasonably high accuracy and computational efficiency in comparison with many other higher‐order‐accurate schemes reported in the recent literature. A von Neumann stability analysis showed that the explicit QUICK scheme has a severe stability constraint which depends upon the diffusion coefficient. It can be proved that this scheme is numerically unstable for the case of pure advection. Various modified forms of the implicit QUICK scheme have been formulated and their numerical stability properties have been studied and analysed. The modified QUICK schemes considered have been tested for transient simulations for the cases of pure advection and of advection and diffusion in an idealized one‐dimensional basin using three different initial boundary conditions: (a) a sharp front concentration gradient, (b) a Gaussian concentration distribution and (c) a plug source. Details of the comparisons between these modified schemes and with other typical second‐order‐accurate difference schemes are given, together with comparisons with the analytical solutions for each case. A two‐dimensional version of the semi‐time‐centred QUICK scheme (ADI‐QUICK), has also been applied to a two‐dimensional test case using the standard ADI technique and has been shown to be attractive in comparison with other comparable second‐order schemes.