The wave equation in second-order form is an important model for many applications in science and engineering involving wave propagation. Example applications include acoustics, electromagnetics and elasticity; such problems are often posed mathematically as partial differential equations with appropriate initial and boundary conditions. Wave propagation problems are often solved most efficiently using high-order accurate explicit time-stepping schemes. Explicit schemes can be fast and memory efficient. The time-step for such schemes is limited by the usual CFL stability condition involving the size of grid cells and the wave speed. Thus, there are some situations, such as with a locally fine mesh or a locally large wave speed, when an explicit scheme with a global time-step is inefficient since a small time-step would be required everywhere. For such situations, we say the problem is geometrically stiff or materially stiff. An example of a geometrically stiff problem is the diffraction of an incident wave from a knife-edge as shown in Fig. 1. The solution of this problem is computed using an overset grid for which there are small grid cells near the tip of the knife-edge. These small cells force the time-step of an explicit method to be reduced by a factor of 20 from that required by the Cartesian background grid. (More information concerning overset grids and our numerical schemes for such grids is given later.) There are two common approaches to overcome this stiffness, local time-stepping (LTS) and locally implicit methods (LIM). LTS methods use a local time-step dictated by the local time-step restriction. LIMs use an implicit method on only part of the domain, usually where the grid cells are smallest.
In this article we develop new locally implicit time-stepping schemes for the wave equation in second-order form on overset grids based on the modified equation (ME) approach. These schemes are high-order accurate single-step schemes that use three time-levels and a compact spatial stencil. The schemes depend on one or more parameters that determine the degree of implicitness; the secondorder accurate scheme depends on one parameter while the fourth-order accurate scheme depends on two parameters. For certain ranges of these parameters the schemes are unconditionally stable in time. A small amount of upwind dissipation is added to the schemes for stability on overset grids. The upwind dissipation can be added in several ways, for example, in a fully implicit manner or in a predictor-corrector fashion where the upwinding is added in a separate explicit step.
Our implicit time-stepping ME scheme, denoted by IME, is combined with a ME-based explicit time-stepping scheme, denoted by EME, in a spatially partitioned manner. The EME schemes we use have a compact stencil and have a time-step restriction that is independent of the order of accuracy. 2 We say that these compact EME schemes are able to take a CFL-one time-step. This is in Local Time Stepping LIM:
Locally Implicit Method contrast to typical linear multi-step methods where explicit higher-order schemes tend to have smaller time-step restrictions, or to popular Discontinuous Galerkin (DG) methods where the time-step restriction typically scales as 1∕(2 + 1) with being the degree of the polynomial basis [1]. For an overset grid, a locally implicit scheme can be used, for example, on boundary-fitted component grids that resolve small geometrical features. The EME scheme can then be used on Cartesian background grids or on curvilinear grids that have similar grid spacings to the background grids. In this way the time-step for the EME scheme is not restricted by the small grid cells used to resolve small geometrical features. In typical applications, the majority of the grid points belong to background Cartesian grids, and the solution on these grid points can be advanced very efficiently with a CFL-one time-step. This can make the hybrid IME-EME scheme much more efficient than using the EME scheme everywhere with a small (global) CFL time-step. We refer to this hybrid scheme as a Spatially Partitioned Implicit-Explicit (SPIE) scheme. Note that the EME scheme is more accurate since ME schemes are most accurate for the CFL number close to one: unlike method-of-lines schemes, the accuracy of ME schemes is degraded for small CFL numbers. The implicit matrix formed by the IME schemes is definite, and it is well suited to a solution by modern Krylov-based methods or multigrid. Normally there is no benefit in using implicit time-stepping and taking a large (greater than one) global CFL time-step for wave propagation problems as the accuracy of the solution is usually degraded. However, there are applications where implicit timestepping methods for the wave equation using a large CFL number can be useful. For example, implicit methods for the wave equation are an attractive option for each iteration step of the WaveHoltz algorithm [2][3][4] which solves for time periodic (Helmholtz) solutions. 3The WaveHoltz algorithm can solve Helmholtz problems for frequencies anywhere in the spectrum without the need to invert an indefinite matrix as is common with many approaches. Each iteration of the WaveHoltz algorithm requires a solution of a wave equation over a given period, and just a few implicit time-steps per period (e.g. 5-10) are needed which leads to very large CFL numbers on fine grids. This is one of our motivations for developing stable IME schemes for overset grids.
There is a large literature on ME, LTS, and LIM schemes for the wave equation. Here we provide a brief synopsis, for further references, see [5,6], for example. Explicit ME schemes for the wave equation go back, at least, to the work of Dablain [7] and of Shubin and Bell [8]. Local time-stepping has most often been used for PDEs that are written as first-order systems in time. LTS has been used for decades with adaptive mesh refinement (AMR) since the pioneering work of Berger and Oliger [9]. Local time-stepping has also been developed, for example, for Runge-Kutta time-stepping [5,6,10], leap-frog time-stepping [11], and arbitrary high-order ADER schemes [1]. Of note is the ME-based LTS method for the wave equation of Diaz and Grote [12], where it was found necessary to have a small overlap of one or two cells between the coarse and fine cells in order to retain the time-step dictated by the coarse mesh. Also of note is the locally time-stepped Runge-Kutta finite difference scheme of Liu, Li, and Hu [10] for the wave equation on block-structured grids. Another reason for using local time-stepping is to couple two difference schemes together. For example, Beznsov and Appelö [13] use local time-stepping for the wave equation to couple a DG scheme (which has a small time-step) with an Hermite scheme (which has a CFL-one time-step). The DG scheme is used on boundary-fitted grids of an overset grid for accurate treatment of the boundary conditions.
Compact implicit difference approximations lead to globally implicit systems (although sometimes with a time-step restriction) and these have been used for wave equations by a number of authors [14][15][16]. To overcome the cost of the global implicit solve, it is common to use locally one-dimensional approximate factorizations such as the alternating-direction-implicit (ADI) scheme [17,18]. A variety of locally implicit methods for wave equations have been developed, for example [19][20][21][22]. For some formulations care is required in coupling the implicit and explicit schemes to avoid an order of accuracy reduction in time. Of particular note is the fourthorder accurate locally implicit method for the wave equation of Chabassier and Imperiale [21]. They use a Finite Element Method (FEM) discretization (with mass lumping to form a diagonal mass matrix) and a mortar element method with Lagrange multipliers to couple the implicit and explicit methods. The implicit ME scheme in [21] is similar to our implicit scheme except that we use finite difference approximations and a more compact approximation (which leads to different stability restrictions). Our approach uses a simple coupling between implicit and explicit regions based on overset grid interpolation. The price for this simpler coupling is that upwind dissipation is needed to ensure stability.
We have been developing high-order accurate algorithms for a variety of wave propagation problems on overset grids. These include the solution to Maxwell's equations of electromagnetics for linear and nonlinear dispersive materials [23][24][25][26][27], the solution of linear and non-linear compressible elasticity [28,29] and the solution of incompressible elasticity [30]. A fourth-order accurate ME scheme for Maxwell's equations in second-order form on overset grids was developed in [31]. Extensions of the implicit and implicitexplicit time-stepping methods developed in this article will be very useful for these other applications, both to treat geometric stiffness and for solving Helmholtz problems using the WaveHoltz algorithm.
The work presented in the remaining sections of this article is organized as follows. Explicit and implicit ME schemes for the wave equation are introduced in Section 2, which also serves to establish some notation. Details of the second and fourth-order accurate IME schemes for Cartesian grids are given in Section 3 where a von Neumann stability analysis is also performed. Methods of upwind dissipation for IME schemes are described in Section 4, and this is followed in Section 5 by a formulation and a GKS stability analysis of our new SPIE schemes. Section 6 discusses the implementation of the new ME schemes for overset grids, and Section 7 provides results of a matrix stability analysis the ME schemes for one-dimensional overset grids. Numerical results are discussed in Section 8 and concluding remarks are offered in Section 9.
We are interested in solving an initial-boundary-value problem (IBVP) for the wave equation in second-order form for a function ( , ) on a domain Ω, with boundary Γ, in space dimensions,
Here, = [ 1 , ..., ] ∈ ℝ is the vector of spatial coordinates, is time, and is the spatial part of the wave operator,
with wave speed > 0. The initial conditions on and are specified by the given functions 0 ( ) and 1 ( ), respectively, and the boundary conditions, denoted by the boundary condition operator , may be of Dirichlet, Neumann, or Robin type with given boundary data ( , ).
We begin with a description of the three-level ME schemes that ignores the specifics of the spatial discretizations. Details of the grids and spatial discretization are left to later sections. The explicit and implicit ME schemes are both based on the standard second-order accurate central difference approximation to the second time-derivative of ,
where Δ is the time-step and + and -are forward and backward divided difference operators in time given by + ( ,
respectively. Expanding the terms in (3) using Taylor series gives the following expansion, 4 + Δ 4 360 6 + ⋯ .
(5)
A.M. Carson, J.W. Banks, W.D. Henshaw et al.
Even time derivatives on the right-hand side of (5) are replaced by space derivatives using the governing equation (1a) to give
To form a th order accurate in time scheme, the expansion ( 6) is truncated to ∕2 terms, and the spatial operators in the resulting truncated expansion are discretized with a time-weighted average of three time levels. For example, second-order accurate explicit or implicit three-level ME schemes take the form
where 2,ℎ is a second-order accurate approximation of on a grid with representative grid spacing ℎ. The coefficients ( 2 , 2 , 2 ) are the weights in the time-weighted average of ( , ). The fourth-order accurate scheme has the form
where 4,ℎ is a fourth-order accurate approximation of and ( 4 , 4 , 4 ) are coefficients involved in the time-average of the correction term. Higher-order accurate schemes for = 6, 8, … can be defined in a similar way but for this article we only consider schemes for = 2 and 4. The explicit ME schemes we use have 2 = 2 = 0 and 2 ≠ 0 for = 1 and 2, while the implicit schemes have 2 ≠ 0 for = 1 or 2.
Truncation error analysis can be used to determine the constraints on the parameters 2 , 2 and 2 for ℎ order accuracy. The truncation error of the = 2 scheme in (7), denoted by 2 ( , ), is
For second-order accuracy in Δ and ℎ we take
These two conditions involving the three parameters ( 2 , 2 , 2 ) give a single-parameter family of second-order accurate schemes discussed further in Section 3. Note that the condition 2 = 2 implies that the schemes are symmetric in time which implies the schemes are time reversible. A similar analysis of the = 4 scheme in (8) leads to the conditions
The four constraints in (11) involving the six parameters ( 2 , 2 , 2 ), = 1, 2, implies a two-parameter family of fourth-order accurate schemes as discussed further in Section 3. Note that these = 4 schemes are also symmetric in time. Choices of the parameters that lead to stable schemes for = 2 and 4 are discussed in Section 3.2.
In order to analyze the proposed IME schemes in more detail we introduce a spatial approximation for Cartesian grids in Section 3.1. This allows us to show the form of the fully discrete schemes and to perform a von Neumann stability analysis in Section 3.2.
Let the domain Ω = [0, 2 ] be a box in dimensions discretized with a Cartesian grid with grid points in each direction. Denote the grid points as = [ 1 ℎ 1 , ..., ℎ ] , (12) for multi-index ∈ ℤ and grid spacings ℎ = 2 ∕ . Let ≈ ( , ) be elements of a grid function at time = Δ . Define the usual divided difference operators to be
where ∈ ℝ is the unit vector in direction (e.g. 2 = [0, 1, 0]).
A.M. Carson, J.W. Banks, W.D. Henshaw et al.
The compact ℎ order accurate discretization of the operator can be written in the form
where, for example, 0 = 1, 1 = -1∕12, 2 = 1∕90 and 3 = -1∕560. The compact second-order accurate approximation to 2 is just the square of the second-order accurate approximation 2,ℎ
Although not used here, note that the compact fourth-order accurate approximation to 2 is not the square of 4,ℎ as ( 4,ℎ ) 2 has a wider stencil than is needed [32].
Given the accuracy requirements (10), we write the fully discrete second-order accurate ME scheme (denoted by IME2) in terms of a single free parameter 2 ,
Note that larger values of 2 correspond to schemes that are more implicit with 2 = 0 being the explicit EME2 scheme. Similarly the fully discrete fourth-order accurate ME scheme (denoted by IME4) involves two free parameters 2 and 4 ,
Larger values of 2 and 4 correspond to schemes that are more implicit with 2 = 4 = 0 being the explicit EME4 scheme.
The stability of the IME schemes ( 16) and ( 17) is now studied using von Neumann analysis, assuming solutions that are periodic in space. Von Neumann analysis expands the solutions in a Fourier series and determines conditions so that all Fourier modes remain stable. There are numerous definitions for stability, but for our purposes here we make the following definition: Definition 1 (Stability). A numerical scheme for the wave equation is stable if there are no Fourier modes with non-zero wavenumber, ≠ , whose magnitude grow in time. For the zero wave-number, = , case the linear in time mode, given by = 0 + 1 for constants 0 and 1 , is permitted since this is an exact solution to the wave equation.
The explicit ME schemes (with 2 = 4 = 0) are known to be CFL-one stable (at least for = 2, 4, 6), meaning stable for
For implicit ME schemes we are generally interested in unconditional stability, that is stability for any Δ > 0. The constraint on 2 for unconditional (von Neumann) stability of the second-order accurate IME2 scheme is summarized by the following theorem.
Theorem 1 (IME2 Stability). The IME2 scheme ( 16) is unconditionally stable on a periodic domain provided
The constraints on 2 and 4 for unconditional stability of the fourth-order accurate IME4 scheme are summarized by the following theorem.
Theorem 2 (IME4 Stability). The IME4 scheme (17) is unconditionally stable on a periodic domain provided
The proofs of Theorems 1 and 2 are given in Appendix A.1 and Appendix A.2, respectively.
The EME and IME schemes described in previous sections have no dissipation and are neutrally stable. As a result, perturbations to the schemes, such as with variable coefficients or when the schemes are applied on overset grids may lead to instabilities. For single curvilinear grids, stable schemes can be defined using special discretizations such as summation by parts (SBP) [33][34][35][36][37] methods or the schemes described in [31]. Overset grids are a greater challenge and in this case we add dissipation for stability. Upwind dissipation for the wave equation in second-order form was first developed in [38] and applied to Maxwell's equations in [23]. An optimized version of the upwind dissipation was developed in [27] and it is this optimized version that we use here as a template for the IME scheme.
We first consider upwind dissipation for the explicit ME scheme to establish our basic approach and to introduce some notation. For the explicit scheme, dissipation is added in a predictor-corrector fashion,
where denotes the (full) spatial operator for the th -order accurate scheme, is an upwind dissipation parameter, and is a dissipation operator, which on a Cartesian grid takes the form
where Δ ± are undivided difference operators corresponding to the divided difference operators defined in (13). Note that the dissipation operator has a stencil of width + 3 compared to the stencil width of + 1 for . The wider stencil for the dissipation reflects the upwind character of the operator [38]. Also note that the addition of the upwind dissipation does not change the order of accuracy of the scheme. The dissipation operator in (21b) acts on an approximation of ( , ). The treatment of this approximation in the predictorcorrector scheme in (21) ensures that the scheme, with dissipation, remains explicit and p th -order accurate. For implicit ME schemes, there is more flexibility in the treatment of this approximation. Two approaches are described in the next subsections.
Upwind dissipation for the implicit ME schemes can be added directly into the single step update (denoted as the IME-UW scheme) as
Here denotes the spatial part of the IME scheme as given in ( 16) and ( 17) for some choice of the parameters 2 and 4 . A von Neumann stability analysis for a Cartesian grid leads to the following result.
Theorem 3. The IME-UW schemes (23) for = 2, 4 on a periodic Cartesian grid are unconditionally stable for any > 0 provided 2 satisfies the conditions of Theorem 1 for = 2, or 2 and 4 satisfy the conditions of Theorem 2 for = 4.
The proof of this theorem is given in Appendix A.3. The monolithic upwind dissipation allows for any > 0 and there are various possible strategies for choosing this value [43].
One disadvantage of the upwind scheme (23) is that the dissipation operator changes the implicit matrix, increasing the stencil size. This may increase the cost of the implicit solve and be undesirable, if for example, one wants to use an existing multigrid solver not designed for this special matrix. Dissipation can be added to the IME scheme is a separate explicit correction step as in (21). Allowing for multiple corrections leads to the implicit-predictor, explicit-corrector upwind scheme (IME-UW-PC)
A.M. Carson, J.W. Banks, W.D. Henshaw et al.
where u denotes the number of upwind correction steps. The sequence of corrections in ( 24) can be combined and written succinctly as
where
The conditions on for stability are specified in the following theorem. The theorem covers the cases of using an implicit or an explicit ME predictor.
Theorem 4. The upwind predictor-corrector scheme (24) is stable on the periodic domain provided the non-dissipative predictor scheme (explicit or implicit) is stable and provided
where u = 2 for u even and u = 1 for u odd, and where is the CFL parameter in coordinate direction ,
The proof of Theorem 4 is given in Appendix A. 4. In practice a reasonable choice might be
where ∈ (0, 1) is a safety factor.
Note from ( 27) that the coefficient of dissipation, , decreases as the CFL parameter increases, and thus less dissipation is added as the CFL number increases. Thus, for large CFL it may become necessary to use more than one correction step.
The IME and EME schemes can be combined in a spatially partitioned manner. For overset grids, the IME scheme is used on certain components grids while the EME scheme is applied on all other component grids. A typical strategy is to employ the EME scheme on background Cartesian grids and any curvilinear grids with grids spacings close to a nominal Cartesian value, and then use the IME scheme on any curvilinear component grids that have a minimum grid spacing that is relatively small as compared to the nominal value.
Let e denote the set of grids that use explicit time-stepping (e.g. Cartesian grids), and let i denote the set of grids that use implicit time-stepping (e.g. curvilinear grids). The SPIE algorithm consists of the following three stages: Stage 1. Update explicit grids,
Stage 2. Update implicit grids, interpolating from the solution on explicit grids from Stage 1,
Stage 3. Add dissipation to all grids. For example, if all grids use a predictor-corrector upwind formulation, then use
Multiple upwind correction steps can also be used as discussed in Section 4.2.
Fig. 2. One-dimensional overset grid used to assess the stability of the SPIE scheme. The explicit scheme is used on the left grid and the implicit scheme is used on the right. Interpolation points are marked as circles.
This section investigates the stability of the SPIE scheme for a one-dimensional overset grid. Normal mode (GKS) analysis [39,40] is used to show that the second-order accurate SPIE scheme is stable when solving the wave equation on an overset grid for a onedimensional infinite domain. Matrix stability analysis on a finite domain is presented later in Section 7 and the result of the matrix analysis supports the conclusions of the GKS analysis discussed in this section.
Consider the one-dimensional overlapping grid for Ω = (-∞, ∞) shown in Fig. 2. Let , = ( + 1)ℎ and , = ℎ denote the grid points for the left and right grids, respectively. Let , , for = , , denote the discrete solutions on the two grids. The grids overlap by a distance ℎ and the solution is interpolated at the interpolation points shown in Fig. 2. The second-order accurate SPIE scheme is used. The explicit (EME2) scheme is applied on the left grid and the implicit (IME2) scheme is applied on the right grid. The discrete equations are
,0 = ,-
where we have assumed that the solution is bounded as | | → ∞. Each individual scheme is assumed to be stable and so we take = Δ ∕ℎ < 1 and 2 ⩾ 0. Note that the IME2 scheme is unconditionally stable for 2 ⩾ 1∕4, but 2 ⩾ 0 is sufficient when < 1.
Before proceeding with the stability analysis, it is useful to first state the following lemma related to the stability of the EME2 and IME2 schemes for the Cauchy problem. Lemma 5.1. Suppose is a root of quadratic equation,
where is defined in terms of some ∈ ℂ by
) .
(32b)
The proof of Lemma 5.1 is given in Appendix A.5. We are now ready to prove the main theorem of this section.
Theorem 5. The SPIE scheme in (31), for the one-dimensional infinite domain overset grid, has no unstable solutions provided < 1 and 2 ⩾ 0.
Proof. We look for unstable mode solutions of the form
for some ∈ ℂ with | | > 1. Note that the same amplification factor must appear in both the left and right grid functions in order to match the interpolation equations (31e) and (31f). Substituting the ansatz (33) into (31a) and (31b) implies,
(34b) Fig. 3. Left: composite of a background grid ( 1 , blue) and a boundary-fitted grid ( 2 , green) in physical space for the domain defined by the interior of the red boundary. The grid points on 1 with green dots interpolate from 2 and the grid points on 2 with blue dots interpolate from 1 . Middle: Plot of 1 showing interpolation points, ghost points (grid points which exist outside the physical boundary), and unused points (grid points which do not affect the computation). Right:
The green boundary fitted grid, 2 , is mapped to a unit square. The plot shows interpolation points and ghost points.
The equations ( 34) can be written as quadratics for ,
for some , = , , with roots denoted by ,± . The general solutions for the left and right sides then take the form (using -= -1
where ± and ± are constants. Note that both equations in (34) are of the form (32) of Lemma 5.1. Since we have assumed | | > 1, it follows from Lemma 5.1 that the roots ,± cannot have magnitude equal to one. Since the product of the roots ,± is one, we can therefore, without loss of generality, take | ,+ | < 1, = , . The boundedness conditions (31c) and (31d) at infinity imply + = 0 and -= 0, reducing the solutions to
Applying the interpolation conditions (31e) and (31f) gives
which implies, assuming -≠ 0 and + ≠ 0, that
This last equation cannot hold since | + + | < 1. Therefore, only the trivial solution remains, thus yielding no unstable solutions with | | > 1. □
The new ME schemes have been implemented for complex geometry using overset grids, which are also known as composite overlapping grids or Chimera grids. As shown in Fig. 3, an overset grid, denoted as , consists of a set of component grids { }, = 1, … , , that cover the entire domain Ω. Solutions on the component grids are matched by interpolation [41]. Overset grids enable the use of efficient finite difference schemes on structured grids, while simultaneously treating complex geometry with highorder accuracy up to and including boundaries. Each component grid, , is a logically rectangular, curvilinear grid defined by a smooth mapping from a unit cube in dimensions (called the parameter space with coordinates ) to physical space , = ( ),
All grid points in are classified as discretization, interpolation or unused points [41]. The overlapping grid generator Ogen [42] from the Overture framework is used to construct the overlapping grid information.
Approximations to derivatives on a curvilinear grid can be formed using the mapping method. Given a mapping = ( ) and its
Derivatives of with respect to are then approximated with standard finite differences. Let denote grid points on the unit cube, where = 0, 1, … , . Let Δ = 1∕ denote the grid spacing on the unit cube with = ( 1 Δ 1 , 2 Δ 2 , 3 Δ 3 ). Let ≈ ( ) and define the difference operators,
where is the unit vector in direction . Second-order accurate approximations to the first derivatives, for example, are
where we assume the metric terms ∕ are known at grid points from the mapping. We do not, however, assume the second derivatives of the mapping are known (to avoid the extra storage) and these are computed using finite differences of the metrics. Using the chain rule, the second derivatives are
Fourth and higher-order accurate approximations are straightforward to form using similar techniques.
Careful attention to the discrete boundary conditions is important for accuracy and stability, especially for the wave equation which has no natural dissipation. We use compatibility boundary conditions (CBCs) which are generally more accurate and stable than one-sided approximations [32]. A simple CBC uses the governing equation on the boundary. More generally, CBCs for the case of the wave equation are formed by taking even time-derivatives of the boundary condition and then using the governing equation to replace 2 by . For flat boundaries with homogeneous Dirichlet or Neumann boundary conditions, CBCs lead to odd or even reflection conditions, respectively. For more details on CBCs see [31,32] for example.
The upwind dissipation operator was introduced in Section 4 and defined in (22) for the case of a Cartesian grid. More generally for a curvilinear grid the upwind dissipation operator is taken as
where
Here, Δ ± are undivided difference operators in the coordinate direction of the parameter space corresponding to the divided difference operators defined in (42).
Implicit three-level ME schemes require two time levels to initiate the time stepping. The solution at = 0 is found directly from the initial condition (1b). The solution at the first time-step = Δ could be found from a Taylor series in time using both initial conditions (1b), and (1c), together with the governing equation (1a). It would be convenient to use an explicit version of this Taylor series approximation to obtain the solution at the first time-step; formally this would not change the stability of the scheme. In practice, however, very large errors can be introduced in the first explicit time-step at = Δ when the CFL number is large, often rendering the full computation useless. Thus the first-time step should be taken implicitly when the CFL number is large. Further, it would be convenient if this implicit first time-step utilized the same implicit matrix as subsequent time steps of the full three-level scheme. In this section we show how this can be accomplished.
Consider the second-order accurate implicit IME2 scheme, re-written here for clarity,
which has the implicit operator
Given the initial conditions,
approximate (51b) to second-order accuracy using
with = 0. Solving for -1 gives
Substituting ( 53) into (49) and dividing by 2 gives the following implicit scheme for the first step ( = 0),
This gives the following update for the first time-step ( = 0)
that uses the same implicit operator 2 as the later time steps.
Consider the fourth-order accurate implicit IME4 scheme, re-written here for clarity,
(56) Scheme (56) has the implicit operator
To approximate the second initial condition (51b) to fourth-order accuracy, we use
To treat the 3 term in (58), take the first time derivative of the governing equation and write it in the form
Using this expression for 3 in (58) gives the approximation
Solving (60) for -1 gives
Substituting (61b) into (56) and dividing by 2 gives the following fourth-order accurate implicit scheme for the first step ( = 0),
(62)
Rearranging (62) gives
and substituting for leads to
The terms in blue 4 are dropped, based on accuracy, and this keeps the stencil compact. Note that scheme (64) has the same implicit operator 4 as scheme (56). If the red term 2,ℎ 1, in (64) is replaced by 4,ℎ 1, , then the form of ( 64) is similar to the usual interior update and the same code can be used for both the first step and later steps with the appropriate choice of coefficients.
For the overset grid results presented in this article the implicit equations are solved either with a direct sparse solver (for problems that are not too large) or with a Krylov space method (bi-CG-stab with an ILU preconditioner). All two-dimensional examples presented in this article use a direct sparse solver. For the three-dimensional sphere problem, a fixed (relative) residual tolerance of 10 -10 is used with the bi-CG-stab method. When iterative methods are used, the initial guess for the iterative solver at each time step is chosen by linear extrapolation in time. In general, the residual tolerance could be chosen based on the expected truncation error in the discrete solution, thus reducing the cost of the iterative solve. However, doing this could affect the stability properties of the scheme and thus this strategy requires some study. We leave this to future work.
The approach we currently use to solve the implicit systems is not efficient since an implicit matrix is formed for all grids points, on both explicit and implicit grids. The implicit matrix entries corresponding to the equations at points treated explicitly are simply the identity. This approach was done so that existing software could be used. A more efficient approach would be to only form a system for the implicit points. Moreover, it is often the case that the implicit points on different component grids are not coupled and in this case multiple smaller implicit systems could be formed.
In this section, matrix stability analysis is used to study the behavior of the new schemes on a collection of one-dimensional overset grids. A large number of overset grids with different grid spacings are considered to determine the stability behavior for a wide range of grid configurations. A scaled upwind dissipation coefficient = is used with ∈ [0, 1] to show how the stability of the scheme depends on the amount of dissipation ranging from no dissipation = 0 to full dissipation = 1. In particular, the results for = 0 show the necessity of upwind dissipation. For some cases, the number of upwind corrections, u , must also be chosen greater than one for stability.
The time-stepping update on an overset grid is written in the form of a single vector update for the active unknowns (i.e. unknowns corresponding to the interior equations) excluding the constraint unknowns (i.e. unknowns associated with the boundary conditions and interpolation equations). For homogeneous boundary conditions, this update is written in terms of the matrices 1 and 2 and the vector of active unknowns at time ,
For a given overset grid, the associated eigenvalues and eigenvectors can be determined and this shows whether discrete solutions grow in time or not. and ℎ = 0.5∕ . In our study of stability for a collection of overset grids, is held fixed and thus ℎ is also fixed. We then select a ratio of the grid spacings = ℎ ∕ℎ which implies = -1 + ℎ
. The value for , and the corresponding value for , is chosen to minimize the grid overlap while maintaining the assumption of an explicit interpolation as discussed in more detail below. The sampling of overset grids is performed for grids for a range of grid spacing ratios ∈ [ min , max ] as noted below.
In Stages 1 and 2 of the SPIE scheme in (30), the solutions for the left and right grids are advanced one time step (without dissipation) and the boundary/interpolation conditions are applied. This step can be written in matrix form as
where is a vector holding all of the unknowns on the two grids (including ghost points and interpolation points). The equations in (66) include the interior equations, boundary conditions and interpolation equations. The matrix 0 is equal to the identity matrix at active points using an explicit scheme, while the matrix has values corresponding to the implicit operators 2 in (50) and 4 in (57) at active points corresponding to the second and fourth-order accurate implicit schemes, respectively. The matrix 0 also includes the boundary conditions and interpolation equations; the corresponding rows in 1 and 2 are zero. For example, the boundary conditions on the left grid are (0)
where the odd symmetry conditions on the ghost points are determined from compatibility conditions. The boundary conditions on the right grid are similar. The values at interpolation points are found using Lagrange interpolation for a stencil of + 1 points. For example, an interpolation point on the left grid is found using a formula of the form
where denotes the left index of the interpolation stencil, chosen to make the interpolation as centered as possible, and , are interpolation weights. The values on the right hand side of (69) are known as donor points. The interpolation is taken to be explicit so that none of the donor points for one grid are interpolation points for the other grid. The interpolation stencil is of width + 1 as needed for a -order accurate scheme when the width of the overlap scales with the mesh spacing [41].
Explicit upwind dissipation is incorporated in Stage 3 of the SPIE scheme. Assuming u applications of the dissipation, the updates of the solution at this stage have the matrix form
where 0 , like 0 , includes the boundary/interpolation equations. Finally, the solution at the new time is
Combining the time step in Stages 1/2 and the corrections in Stage 3 leads to a three-level matrix equation of the form where 1 and 2 are coefficient matrices generated from the ones in (66) and (70). The constraint unknowns in (72) can be eliminated by row operations and this leads to the compressed form in (65). The correctness of the matrices 1 and 2 in (65) is checked by comparing, at each time-step, the solution computed using the SPIE scheme in (30) with the solution arising from the compressed form (65).
To investigate the growth of solutions to the discrete problem (65) we look for solutions of the form = 0 which leads to a quadratic eigenvalue problem for given by
This quadratic form can also be written as a regular eigenvalue problem of twice the dimension,
The eigenvalue problem in ( 74) is solved easily with standard software.
For the stability studies, we set = 10 for the right grid with fixed domain Ω = [0.5, 1]. This right grid represents the local boundary grid in a general overset grid. The grid spacing on the left grid is determined by ℎ = ℎ , where is the ratio of grid spacings. This ratio is varied from 1∕2 to 2 to represent typical overset grids where the grid spacings in the overlap are chosen to be nearly the same. For each value of , the parameters and for the left grid are determined based on the grid overlap as discussed above. Finally, the parameter for the scaled upwind dissipation coefficient
is varied between 0 and 1 to study how much dissipation is needed to stabilize the SPIE scheme for the different cases considered.
Fig. 5 shows some sample results for the SPIE2-UW-PC scheme. The left grid uses the EME2 explicit scheme with CFL = 0.9, while the right grid uses the IME2 implicit scheme. Explicit upwind dissipation is added in a corrector step with given by ( 29)
and u = 1 corrections. The safety factor for is chosen as = 0.9. Two grid cases are shown for = 0.3; one with grid-ratio = 0.5 and one with = 1.55. The grid plotted on the top of the figure is stable as shown in the middle left plot; all eigenvalues of satisfy | | ⩽ 1 + tol , where tol = 10 -8 . The grid on the bottom has two unstable modes, as illustrated on the middle right plot. The conclusion for this representative case is that there is insufficient dissipation for the SPIE2-UW-PC scheme with = 0.3 and u = 1.
Results are now presented using different combinations of explicit and implicit schemes; different orders of accuracy ( = 2 and = 4), and different numbers of upwind corrections. The following cases are considered 1. EMEp: explicit p th -order accurate ME schemes are used on the left and right grids. 2. IMEp: implicit p th -order accurate ME schemes are used on the left and right grids. 3. SPIEp: An EMEp scheme is used on the left grid and an IMEp scheme is used on the right grid.
Explicit upwind dissipation is used in all cases. Unless otherwise specified, the IME schemes use the parameters, 2 = 1∕4 and different overset grids). Fig. 6 shows results for the EMEp scheme. The time-step is chosen so that the CFL number is 0.9 on the side with the smallest grid spacing. The number of unstable grids is plotted versus , the scaling factor of the dissipation coefficient . For = 2 and = 0 (no dissipation), roughly 60% of the grids tested are unstable. This number drops to about 25% when = 0.1, and there are no unstable grids for ⩾ 0.3. For = 4 a value of about = 0.5 is sufficient to stabilize all the grids tested. Fig. 7 shows results for IMEp schemes. Note that these schemes use a single implicit solve over both grids (i.e. the implicit solves are coupled, not partitioned). For = 2 ( = 4), the time-step is chosen so that CFL number is 4.0 (5.0) on the side with the smallest grid spacing. The number of explicit upwind corrections is set to u = 4 for = 2 and u = 5 for = 5. This choice is made since the value of for the IME schemes scales with the inverse of the CFL. The results in Fig. 7 show that the schemes have no unstable grids for ⩾ 0.1.
Fig. 8 shows results for the implicit-explicit SPIE scheme. The left grid uses an explicit solver and the overall time-step is chosen to match a CFL number of 0.9 on this grid. The CFL number on the implicit grid varies between grids and reaches a maximum of 1.8. The This choice for 2 and 4 is made to reduce the magnitude of the coefficient in the leading term of the truncation error (see [43]). For a grid that uses an implicit method, the coefficient of upwind dissipation is chosen as
where is the safety-factor and where is the CFL number for grid , which on a Cartesian grid is given by
For grids using an explicit method we take
≈ 1 for such grids according to the CFL condition. Unless otherwise stated, all computations use a safety factor = 1 and = 1 upwind correction steps.
We begin with numerical results illustrating the accuracy and stability of the second and fourth-order accurate IME and SPIE schemes.
In this section, eigenmodes of the unit disk in two dimensions are computed. We look for time-periodic solutions to the wave equation. In polar coordinates ( , ), these solutions have the form Fig. 13. Left: Grid convergence for the IME-UW-PC scheme. Right: Grid convergence for the SPIE-UW-PC scheme.
We consider the scattering of a plane wave from a cylinder of radius in two dimensions. The incident field is taken to be
where is the wave number of the incident field in the reference direction given by . The exact solution is written in polar coordinates ( , ) with the usual assignment = cos . A homogeneous Dirichlet boundary condition on the cylinder is assumed so that the total field (incident plus scattered) is given by
where 0 = 1 and = 2 for > 0, and Fig. 12 shows the overset grid (2) scat and contours of the computed solution and errors at = 1 using the fourth-order accurate SPIE-UW-PC scheme. The grid (8) scat is used for this calculation with = 10. The errors are seen to be smooth. Fig. 13 shows grid convergence results at = 0.4 for an incident field with = 2 using the second and fourth-order accurate IME-UW-PC and SPIE-UW-PC schemes. Max-norm errors are plotted as a function of the grid spacing and the solutions are seen to be converging at close to the expected rates.
We now consider eigenmode solutions of a solid unit sphere assuming a homogeneous Dirichlet boundary condition on the surface of the sphere. Introduce spherical polar coordinates ( , , ), where is the radius, ∈ [0, 2 ] is the angle in the -plane and ∈ [0, ] the angle from the -axis. We assume time-periodic eigenmodes with frequency having the well known form , , ( , , , ) = -1∕2
where . The frequency of vibration is given by =
, . The initial conditions use the exact solution and its time derivative at = 0.
The composite grids for the solid sphere of radius one, denoted by ( ) , consist of four component grids, each with grid spacing approximately equal to Δ ( ) = 1∕ (10 ). The sphere is covered with three boundary-fitted patches near the surface as shown on the left in Fig. 14. There is one patch specified using spherical polar coordinates that covers much of the sphere except near the poles. To remove the polar singularities there are two patches that cover the north and south poles, defined by orthographic mappings. A background Cartesian grid (not shown) covers the interior of the sphere. The middle image in the figure shows the solution at = 0.5 for the eigenmode with ( , , ) = (2, 1, 1) and , ≈ 5.7634591968945. This solution is computed using the fourthorder accurate SPIE-UW-PC scheme and the grid (4) . The right image shows the max errors which are smooth (Fig. 14).
In this section we perform some very long-time simulations to confirm numerically that the solutions computed using the IME and SPIE schemes with upwinding remain stable and bounded. Initial conditions are chosen with random grid values on [0, 1] so that all eigenmodes, including any possible unstable ones, would be seeded with an order one amount of energy. The numerical schemes are integrated to very long times and the solutions are monitored for any growth. Due to the upwinding, the magnitude of the computed solutions for a stable scheme is expected to decay slowly to zero over time.
To assess the growth or decay of the solution we plot a discrete approximation to the energy given by
where ‖ ⋅ ‖ Ω denotes the 2 -norm over the domain Ω. We note that the energy defined in (84) remains constant in time for exact solutions of the wave equation on Ω assuming homogeneous Dirichlet or Neumann conditions specified on the boundary of Ω. For purposes of this study, a first-order accurate approximation to ( 84) is sufficient, denoted by ℎ . The backward time-difference - is used to approximate the time derivative in (84) and first-order accurate backward differences are used to approximate the spatial derivatives, for example on a Cartesian grid ≈ - . Note that the discrete energy ℎ would remain approximately constant if the scheme is stable, but with upwind dissipation included the discrete energy is expected to decay over time. Fig. 16 shows results from some long-time simulations for both the SPIE and IME schemes. In all cases the schemes remain stable.
The left plot shows the discrete energy ℎ ( ) over time for a computation on the disk grid (4) disk as described in Section 8.1.1. The final time is = 10 4 for the simulation and approximately 6 × 10 5 time-steps are used. Results are shown for the second and fourthorder accurate SPIE-UW-PC schemes. The discrete energy is seen to decay rapidly at first as the high-frequency components of the solution are damped by the high-order upwind dissipation. As time progresses the solution becomes smoother and the energy decays more slowly. The discrete energy for the fourth-order accurate scheme decays more slowly than the second-order accurate scheme since its dissipation scales as (ℎ 5 ) compared to (ℎ 3 ) for the second-order accurate scheme. The middle plot shows results for the three-dimensional solid sphere grid (4) as described in Section 8.1.3. In this case the final time is = 10 3 and the calculation requires approximately 10 5 time-steps. The results show that the discrete energy decays and schemes remain stable for the spherical case in qualitative agreement with the results for the disk case.
The right-most plot of Fig. 16 compares the energy decay for the (fully implicit) IME4-UW-PC scheme on the disk grid (4) disk for three different values of the CFL number, 1, 5 and 10. In each case the scheme remains stable and the discrete energy ℎ decays. The dissipation parameter is the same for each case. Note that the CFL=10 run takes 10 times fewer time-steps than the CFL=1 run, and thus the dissipation has fewer time-steps to act.
We now turn our attention to a set of examples that posses some geometric stiffness. For such problems it is demonstrated that the SPIE scheme can be much faster than the fully explicit schemes. Importantly, it is also shown that the accuracy of the computed solutions from the SPIE scheme are, in general, quite similar to the accuracy of the explicit ME solutions. Thus, at least for the cases shown here, taking a large CFL time-step using an implicit ME method in small parts of the domain where geometric stiffness occurs does not appear to have a significant effect on the overall accuracy.
This section studies the accuracy and performance of the SPIE scheme for the scattering of a plane wave from a small cylinder in two space dimensions. The incident wave, exact solution, and overset grid topology were described previously in Section 8.1.2. The results presented in this section show that the SPIE scheme can compute solutions much faster than the EME scheme but with similar errors. Fig. 17. Scattering from a small hole. Left: overset grid (2) scat and magnified views. Right: max-norm errors over time for the EME4 and SPIE4 schemes on grid (4) scat .
The SPIE4 scheme achieves similar errors to the EME4 scheme but at a factor 21 reduced CPU cost. A very small cylindrical hole of radius = 0.01 sits at the center of a square domain [-2, 2] 2 . The overset grid is shown in Fig. 17. An incident field with wave-number = 10 impinges on the hole where a homogeneous Dirichlet boundary condition is applied. The exact solution is imposed on the outer boundaries of the square. The solution is computed to a final time of = 10 using the EME4-UW-PC scheme and the SPIE4-UW-PC scheme with u = 2 upwind corrections. 5 For the SPIE scheme, implicit time-stepping is used for the boundary-fitted annular grid with radial stretching near the small hole, while explicit time-stepping is for the Cartesian background grid. Fig. 18 shows the computed solution at 1 for total field, , the scattered field, , and the error in . Note that there is a significant scattered field for this case even though the radius of the cylinder = 0.01 is fairly small compared to the wavelength, 2 ∕ ≈ 0.63, of the incident field. The error is seen to be smooth with the largest errors distributed throughout the domain; there are no particularly large errors in the vicinity of the hole.
The right graph in Fig. 17 compares the max-norm errors over time for the EME4 and SPIE4 schemes. The error in the EME scheme starts out smaller but then becomes similar in magnitude to the errors in the SPIE4 scheme. The time-step for the EME4 scheme is approximately 30 times smaller than that for the SPIE4 scheme, and the CPU time required to compute the solution at = 10 using the SPIE4 scheme is approximately 20 times smaller than that needed for the EME4 scheme.
These results thus illustrate that for longer times the problem can be considered geometrically stiff since an implicit algorithm can obtain the required accuracy faster than an explicit one. On the other hand, if high accuracy is needed over short times then a fully explicit scheme may be more efficient.
We consider the scattering of a modulated Gaussian plane wave from two different arrays of small holes. This example demonstrates an interesting scattering problem for a geometry with small geometric features for which the SPIE scheme gives a good speedup over the explicit scheme. The incident field consists of a modulated Gaussian plane wave traveling from left to right and given by the formula (85) where the Gaussian shape parameter is = 20, the modulation wave-number is 0 = 8, and the center of the pulse is at 0 = -2 initially. The initial conditions use (85) and its time derivative at = 0. Two configurations of holes are considered for the region [-3, 2] × [-2, 2]. The first, called the aligned-hole configuration, contains an array of × holes, each of radius = 0.01, with = 7 and = 26. The centers of the holes are located at , = -
where = 0.15 and = 0.15 denote the hole separations in the and directions, respectively. The second configuration, called the configuration, also = 7 columns of holes, but every second column is shifted vertically ∕2 and contains 27 holes instead of 26.
The overset grid for the aligned-hole configuration, denoted by ( ) h,a , consists of a background Cartesian grid for the region [-3, 2] × [-2, 2], together with small annular grids around the holes as shown in Fig. 19. The nominal grid spacing is Δ ( ) = 1∕ (10 ) with the grid lines on the annuli slightly smaller and clustered near the boundary as shown in the figure. The overset grid for the offsethole configuration, denoted by ( ) h,o , has a similar construction to that of the aligned-hole grid following its placement of holes. The boundary conditions are taken as Dirichlet on the holes, Dirichlet on the left and right ends of the outer rectangle and periodic in the y-direction of the outer rectangle.
Fig. 20 shows the solution at three times for the two grid configurations of aligned and offset holes. The solution is computed with the SPIE4-UW-PC scheme on grids (16) h,a and (16) h,o . Note that there are some edge effects in the solutions near the top and bottom periodic boundaries of the domain, due to the arrangement of the hole grids near these boundaries. The solution at = 1 shows the incident Gaussian plane wave just starting to impact the first column of holes. At = 2 the wave has traveled through most of the holes and a reflected wave is beginning to appear. By = 3.5 most of the incident wave has been reflected or transmitted, although some residual wave motion resides within the array of holes. Perhaps surprisingly, the transmitted wave is much stronger for the offset arrangement of holes. Returning to Fig. 19, the right plot compares contours of the solutions computed using the SPIE and EME schemes. The top half of the plot shows the SPIE4 solution, while the bottom half shows the EME4 solution. After accounting for the reflection symmetry about the horizontal centerline, the results are nearly indistinguishable. The speedup of the SPIE scheme over the EME scheme was about a factor of 2 for this case. The SPIE time-step is about 4 times that for the EME scheme. A better implementation of the implicit solvers should show an even bigger speedup of perhaps a factor of 3 or more (see the comments in Section 6.4).
In this example, a modulated Gaussian plane wave given by (85) travels from left to right and diffracts off a thin knife edge as shown in Figs. 1 and 21. This example demonstrates a problem that is geometrically stiff due to a sharp corner in the domain geometry, and one for which only a small portion of the overset grid is treated implicitly.
The overset grid for the geometry, denoted by ( ) ke , is shown in Fig. 1, and consists of four component grids. A background Cartesian grid covers the domain [-1.25, 1] × [0, 1]. Two other Cartesian grids lie adjacent to the lower sides of the knife edge which has a total height of 0.5. A curvilinear grid is used over the tip of the knife edge. The nominal grid spacing is Δ ( ) = 1∕ (10 ), although the tip grid uses a finer mesh with stretching to resolve the sharp tip of the knife edge. Fig. 21 shows contours of the solution for three times computed on grid (16) ke using the fourth-order accurate SPIE4-UW-PC scheme.
The Gaussian is centered at 0 = -0.75 initially and the Gaussian shape parameter is taken as = 80. Neumann boundary conditions are used on the outer boundaries of the domain and a Dirichlet condition is used on the knife edge. The initial conditions use (85) and its time derivative at = 0. The SPIE scheme is used with only the curvilinear tip grid treated implicitly. As a result, the scheme is able to use a time-step that is about 20 times larger than that required by the fully explicit EME-UW-PC scheme. The speedup factor Fig. 20. Scattering of a modulated Gaussian plane wave by small holes (the white dots are small holes with a grid around each as shown in Fig. 19). Top: Offset holes. Bottom: Aligned holes. over the fully explicit scheme is found to be about 11 for both the second and fourth-order accurate SPIE schemes. Note that the tip grid has just 1, 760 grid points out of a total of 928, 765, or 0.2% of the points. Obviously a more efficient implementation of the implicit solver should lead to speedups closer to a factor of 20, a task for future work.
We have described and analyzed a class of new implicit and implicit-explicit time-stepping methods for the numerical solution of the wave equation in second-order form. These single-step, three time-level, schemes are based on the modified equation (ME) approach. Second and fourth-order accurate schemes are developed, although the approach supports higher-order accurate schemes. The coefficient matrix implied by the implicit scheme is definite and well suited for solution by modern Krylov methods or multigrid. Conditions for accuracy and unconditional stability of the implicit ME (IME) schemes are derived. Several approaches for incorporating upwind dissipation into the IME schemes are discussed. A predictor-corrector approach that adds the upwinding in a separate explicit step appears to be quite useful. For problems on overset grids that are geometrically stiff due to locally small cells, we have developed
Many EME schemes take powers of an matrix operator (leading to a wider stencil) and the time-step restriction depends on the order of accuracy.
Note that the dispersion errors due to the large CFL time-step can be eliminated by an adjustment to the forcing frequency.
For interpretation of the references to colour please refer to the web version of this article.
= 1∕12. For each value of upwind scaling factor , the grid-ratio is varied form 0.25 to 2 using = 101 different values (i.e.
Using u = 1 is not sufficient for stability for this problem; as noted in Section 4.2, the dissipation coefficient on implicit grids decreases with increasing CFL number and for large CFL numbers it may be necessary to use a larger u .