Unformatted Document Excerpt
>> New Jersey
Course Hero has millions of student submitted documents similar to the one
below including study guides, practice problems, reference materials, practice exams, textbook help and tutor support.
Course Hero has millions of student submitted documents similar to the one below including study guides, practice problems, reference materials, practice exams, textbook help and tutor support.
J. SIAM NUMER. ANAL. Vol. 38, No. 6, pp. 17631783
c 2001 Society for Industrial and Applied Mathematics
PARTIALLY IMPLICIT BDF2 BLENDS FOR CONVECTION DOMINATED FLOWS
WILLEM HUNDSDORFER Abstract. In this paper we consider various blends of implicit and explicit time integration schemes, based on the well-known BDF2 method, applied to convection-diffusion problems with dominating convection. A fully implicit treatment of convection terms is often not very efficient. We shall deal with second order schemes that are implicit in the convection terms only locally in space, without introducing the internal inconsistencies that are common with many time-splitting methods. Along with implementation aspects of the implicit relations, we shall discuss accuracy of the schemes, positivity and monotonicity properties. Key words. numerical analysis, initial-boundary value problems, BDF methods, implicitexplicit methods, splitting methods AMS subject classifications. 65M06, 65M12, 65M20 PII. S0036142999364741
1. Introduction. When adopting the method of lines approach, space discretization of multidimensional, time-dependent partial differential equations (PDEs) results in large systems of ordinary differential equations (ODEs) which are to be integrated in time by an appropriate time stepping scheme. Frequently in such applications one is confronted with problems having both stiff and nonstiff parts. Diffusion, for example, leads to stiff terms that need implicit treatment. Convection terms can usually be taken explicitly, but if we have locally large convective velocities an explicit treatment is unfavorable due to the CFL restrictions on stability, whereas a fully implicit approach leads to systems of algebraic equations that are rather difficult to solve numerically. Here we shall deal with partial implicit treatment of convective terms in such a way that the resulting scheme is fully implicit only in those spatial regions where the solution is smooth and the convective velocities are large. The focus in this paper is on convection dominated equations. First, consider the convection equation without any diffusion, (1.1) ut + (q(x, t)f (u)) = 0, x , t 0,
on a spatial domain Rd with appropriate initial and boundary conditions. Here q(x, t) Rd is a given velocity and f is a scalar flux function. Discretization of the spatial derivatives leads to a large system of ODEs, the so-called semidiscrete system, (1.2) w (t) = F (t, w(t)), t 0,
where F contains the discretized convective terms, and an initial value w0 = w(0) is given. We consider numerical time integration schemes with step size > 0, yielding approximations wn w(tn ) at the time levels tn = n . For spatial discretization we shall deal with limited second order finite volume or finite difference formulas. The dimension of the semidiscrete system is proportional to the number of grid points, and components wi (tn ) of w(tn ) refer to approximations at the grid point xi or to an
Received by the editors November 29, 1999; accepted for publication (in revised form) August 21, 2000; published electronically January 5, 2001. http://www.siam.org/journals/sinum/38-6/36474.html CWI, P.O. Box 94079, 1090 GB, Amsterdam, The Netherlands (Willem.Hundsdorfer@cwi.nl).
average value on a cell i around xi . With multidimensional problems i will denote a multi-index. One of the most popular implicit methods for solving (1.2) is the second order BDF2 method (1.3)
3 2 wn
- 2wn-1 + 1 wn-2 = F (tn , wn ) 2
with n 2; see . Along with w0 , this two-step method needs w1 as starting value. It can be computed by a one-step method, for instance, implicit Euler. The popularity of this BDF2 method is due to its stability and damping properties; see , for instance. These are crucial properties for efficient solution of diffusion equations. Convection equations, on the other hand, are often treated more efficiently by an explicit method. Here we consider the related second order scheme (1.4)
3 2 wn 1 - 2wn-1 + 2 wn-2 = F (tn , wn ),
wn = 2wn-1 - wn-2 ,
to which we shall refer as the explicit BDF2 method. Note that wn = 2wn-1 - wn-2 is just an explicit prediction by linear extrapolation. As with any standard explicit method, we now have a CFL condition for stability. Therefore, if we deal with large velocities or fine spatial grids, very small time steps have to be taken. As we shall see, the fully implicit method also gives us difficulties when applied to large Courant numbers. This is due to slow convergence of the Newton iterations for the implicit relations but also due to loss of monotonicity. In this paper we therefore consider a partially implicit convection treatment, where only those parts in the domain with little spatial variation in the solution are treated implicitly. The resulting formula is (1.5)
3 2 wn 1 - 2wn-1 + 2 wn-2 = F (tn , wn + (I - )wn ),
where is a diagonal matrix with entries i = 0 if the convection term is taken explicitly at the grid point xi , and i (0, 1] otherwise. The actual choice for the i is discussed in section 4. With convection-diffusion problems, (1.6) ut + (q(x, t)f (u)) = (D(x, t, u) u),
the resulting semidiscrete system will be of the form (1.7) w (t) = F (t, w(t)) + G(t, w(t)), t 0,
where F contains the convective terms and G denotes discretized diffusion. The above formula (1.5) for the convection part can be well combined with implicit treatment of the diffusion term by considering (1.8)
3 2 wn
- 2wn-1 + 1 wn-2 = F (tn , wn + (I - )wn ) + G(tn , wn ), 2
so that we obtain a formula that benefits from the damping properties of the fully implicit BDF2 scheme for the diffusion part. If = O this is an implicit-explicit method of the type that was introduced by Crouzeix  and Varah . Stability results can be found in [1, 5, 8, 21], for example, and a practical application in the field of air pollution was discussed in . In general, the stability of this method is completely determined by the CFL restriction for the explicit convection part.
PARTIALLY IMPLICIT BDF2 BLENDS
Note that all of the above methods are different from the usual time-splitting techniques, where different subproblems, such as v (t) = F (t, v(t)) and v (t) = G(t, v(t)), are solved subsequently on small time intervals. This leads to intermediate results which have little physical meaning, since they are not consistent with the total equation. Boundary conditions or interface conditions are usually lacking for these intermediate results. With the above BDF2-type methods we use only fully consistent approximations wn and no intermediate results. Further we note that if = I the above formula (1.5) is a two-step extension of the more familiar -method (1.9) wn - wn-1 = F (tn+ , wn + (1 - )wn-1 )
with the explicit Euler method, = 0, and the implicit Euler method, = 1, as boundary cases. We shall not consider these methods here since both the implicit and explicit Euler methods are not well suited for convection problems. The implicit Euler method is much too diffusive, whereas the explicit Euler method is unstable for the spatial convection discretizations considered in this paper. Actually, in a method of lines setting, the explicit Euler method is unstable for all well-known spatial convection discretizations except for the (diffusive) first order upwind discretization. A related method has been formulated by Blunt and Rubin  for one-dimensional (1D) problems, where the implicit Euler scheme was combined with an explicit, direct space-time scheme (LaxWendroff-type) with limiting. However, for multidimensional problems this combined scheme needs dimensional splitting since the formulation of such a direct space-time scheme for multidimensional problems is different than with the implicit Euler scheme; see also . Moreover, due to the use of implicit Euler, the order is 1 at most. In this paper we shall consider the second order BDF2 blends (1.5) mainly for purely convective problems. If diffusion is added as in (1.8), the method becomes implicit over the whole spatial domain, but in those regions where the entries i are zero the implicit relations have a nice symmetric, diagonally dominant structure, so that standard linear solvers, such as conjugate gradients, will be very efficient. Spatial discretization of the convective terms will be done by limiting in order to avoid oscillations and negative solution values. In section 2 we discuss by means of 1D examples implementation issues and qualitative behavior. As we shall see, the standard implicit BDF2 method (1.3) becomes rather expensive, and, more importantly, the results are also rather disappointing with respect to qualitative behavior and accuracy. This is due to the poor monotonicity properties of the standard implicit BDF2 method. In section 3 we consider formula (1.5) with = I, with the aim of selecting values of with better monotonicity properties than = 1. To obtain theoretical results we shall concentrate on positivity for linear systems. The results in this section can be regarded as an extension of the positivity theory of Bolley and Crouzeix . In section 4 we consider implementations of (1.5) with variable entries i . The actual choices will be motivated by the preceding results. We shall discuss the accuracy of the schemes with variable entries in some detail in section 5, since the standard local truncation error no longer gives proper information about the accuracy of these schemes. This is similar to the situation for stiff ODEs as considered in Hundsdorfer and Steininger . Numerical results will be presented in section 6 for a test example from reservoir simulation, where we have locally large convective velocities q near injection and production wells and moderate or small velocities elsewhere in
the spatial region. It will be seen that the locally implicit schemes can be much more efficient than the fully implicit counterparts such as (1.3), whereas this locally implicit approach allows step sizes much larger than with explicit schemes such as (1.4). 2. One-dimensional examples. In this paper we shall deal with convectiondiffusion discretizations for 1D or two-dimensional (2D) problems. For ease of presentation we first consider the 1D convection problem (2.1) ut + (q(x, t)f (u))x = 0,
on = [0, 1], with monotonically increasing flux function f . Further, it is assumed that an initial profile u(x, 0) and appropriate boundary conditions are given. In this section we shall discuss the advantages and disadvantages of the implicit BDF2 method (1.3) compared to its explicit counterpart (1.4). 2.1. The spatial discretizations. For the spatial derivative in (2.1) we consider discretizations in flux form on a uniform mesh, (2.2) wi = 1 1 q 1 f (wi- 1 ) - qi+ 1 f (wi+ 2 ) , 2 2 h i- 2
1 with grid points xi = ih and qi1/2 = q(xi 2 h, t). Here wi = wi (t) stands for a semidiscrete approximation to the average value of u(x, t) over the cell i = [xi - 1 1 2 h, xi + 2 h]. The choice for the cell boundary values wi1/2 determines the actual discretization. It is well known that the first order upwind approximation wi+1/2 = wi , for q > 0, gives very inaccurate and diffusive results. On the other hand, higher order linear discretizations, such as second order central wi+1/2 = 1 (wi + wi+1 ) or second 2 order upwind wi+1/2 = 1 (-wi-1 + 3wi ), give results that are very oscillatory. For 2 that reason, discretizations with limiters have become increasingly popular. In the following, let
wi - wi-1 . wi+1 - wi
In (2.2) we shall deal with limited approximations for the cell boundary values of the form (2.3) w
i+ 1 2
1 wi + 2 (i )(wi+1 - wi ) 1 2 (1/i+1 )(wi
- wi+1 )
if qi+ 1 0, 2
1 if qi+ 2 < 0,
where is the limiter function. For this limiter function two choices are considered, (2.4) () = + || , 1 + || 2 1 + , 2 3 3 .
() = max 0, min 2,
The first limiter is due to van Leer , the second to Koren . The limiters provide a suitable balance between the monotone first order upwind flux and higher order fluxes. Formal statements on accuracy are difficult, due to the built-in switches,
PARTIALLY IMPLICIT BDF2 BLENDS
but simple numerical tests for smooth solutions show that the spatial discretizations are approximately second order in the L2 -norm. With both limiters we have w(t) 0 whenever w(0) 0, together with monotonicity properties such as the total variation diminishing (TVD) property; see, for instance, [15, 17] for more details. For points adjacent to the boundaries, some of the wj values that are needed in (2.3) might be missing, and for those, constant extrapolation is used, which means that we switch locally to first order upwind. The above discretizations extend easily to more dimensions on Cartesian meshes. We observed that the explicit BDF2 method (1.4) is stable with these spatial discretizations up to Courant number 1/2, approximately. This is an experimental bound; precise results can be obtained for the corresponding linear nonlimited discretizations; see [8, 22]. 2.2. Implementation. For test purposes we consider the linear 1D convection problem, (2.1), with (2.6) f (u) = u, q 1.
Note that even for this linear problem the resulting semidiscrete system will be nonlinear, due to the limiter. Therefore, with implicit time integration, some form of Newton iteration is required, which in turn needs an approximation to the Jacobian 3 matrix A w F (t, w). Within the Newton iteration for (1.3) the matrix 2 I - A is used. The first choice to be considered is the first order upwind approximation A = A1 q h 1 -1 0
in stencil notation. The resulting iteration scheme is related to the defect correction approach used in [6, 18], for instance. Other choices for the Jacobian approximation can be obtained by realizing that the above flux formulas are nonlinear counterparts of formulas obtained by linearizing around = 1 (replacement of () in (2.3) by (1) + (1)( - 1)). For the van Leer limiter (2.4) this leads to A = A2 q 4h -1 4 -1 -2 0
corresponding to the linear Fromm scheme. For the Koren limiter (2.5) we get A = A3 q 6h -1 6 -3 -2 0 ,
which corresponds to the well-known linear third order upwind-biased scheme. Finally, we also consider the choice A = 0, which gives standard functional iteration. In Table 2.1 the average number of Newton iterations per step are listed for the implicit BDF2 method (1.3) with these various choices and several Courant numbers = /h. As starting procedure to calculate w1 , the implicit Euler method was taken. The solutions were calculated on the spatial interval [0,1] with periodicity. The results are given here for an initial block-profile u(x, 0) = and for a smooth initial profile u(x, 0) = sin2 (x). 0 1
1 for 0 < x < 2 , otherwise,
In this test, the mesh width has been chosen as h = 1/100 and output time is T = 1 4 . The convergence criterion for the Newton iteration is that the max-norm of the residual should be less than 10-6 . This is rather strict, but accurate solution of the implicit relations is necessary to maintain the monotonicity of the limiting procedure. The maximum number of Newton iterations per step is set to 100. If convergence is still not reached, then the calculations are aborted and is used for the corresponding entry in Table 2.1. Actually, with A = 0, = 1 this means genuine divergence, with the other cases in the table extremely slow convergence.
Table 2.1 Linear convection test (2.1), (2.6) with implicit BDF 2 method. The entries are the average number of Newton iterations per step with block-profile and sin2 -profile, respectively.
Limiter (2.4) (2.4) (2.4) (2.5) (2.5) (2.5)
A A1 A2 0 A1 A3 0
=1 10.8 - 8.0 13.9 - 11.0 - 14.7 - 11.0 - -
= 1/2 8.6 - 6.6 9.3 - 6.5 23.7 - 13.5 - 7.5 24.6 - 12.2 -
= 1/4 6.9 - 4.5 6.8 - 4.3 7.9 - 5.7 8.4 - 4.9 9.5 - 5.3 9.1 - 6.8
The first observation from Table 2.1 is that the choices A = A2 and A = A3 do not perform well. Especially with (2.5) and A = A3 we get a convergence behavior that is hardly better than with functional iteration. The only choice that does perform reasonably here is A = A1 . Moreover, we see that the algebraic relations with limiter (2.4) are easier to solve than with (2.5). It should be noted that the latter gives slightly better results with respect to accuracy, with somewhat less numerical diffusion, but the differences are small. Even with explicit methods the limiter (2.5) is more expensive than (2.4), due to the max-min calculations. Therefore we consider in the following only the limiter (2.4) with first order upwind approximation for the Jacobian. This implementation seems quite robust. For example, in the above test, if only one time step is performed, = T , = 25, the Newton process still converges (with 16 iterations for both profiles). Moreover, with first order upwind approximations for the Jacobian the resulting linear system is diagonally dominant, which is of importance in more space dimensions in connection with iterative linear solvers. However, even with this choice a rather large number of Newton iterations is needed per step. Note that in the above test, the explicit version of the BDF2 method could be used up to Courant number = 1/2, and with this explicit method the CPU time per step is much smaller than with the implicit scheme. Therefore, we can conclude that accurately solving the implicit relations with limiting is expensive in terms of CPU time. Some gain could be achieved by setting the tolerance in the convergence criterion to less strict values, but it was observed that even with small Courant numbers negative values arise that are of the same order of magnitude as this tolerance. Numerical tests in 1D with Burgers and BuckleyLeverett equations gave results comparable to those in Table 2.1. With multidimensional problems we shall adopt the same implementation as above. The Jacobian required in the Newton iteration is approximated by the Jacobian that corresponds to first order upwind spatial discretization.
PARTIALLY IMPLICIT BDF2 BLENDS
2.3. Qualitative behavior. The advantage of an implicit time stepping method is the possibility of taking large step sizes without introducing instabilities. However, in several numerical tests we observed that the quality of the implicit solutions are rather poor with large, or even moderately large, Courant numbers if the solution has steep gradients. As an example, consider the 1D BuckleyLeverett equation given by (2.1) with (2.7) and initial block-profile u(x, 0) =
1 0 for 0 < x < 2 , 1 otherwise.
f (u) =
3u2 , 3u2 + (1 - u)2
At the inflow the boundary condition is u(0, t) = 1. For the mesh width we take h = 1/100 and the endpoint in time is T = 1 . In the following figures the numerical 4 solutions are plotted with solid lines. Dashed lines are used to indicate the reference solution that uses the same mesh width h but computed with a very small time step; this corresponds to the exact solution of the semidiscrete system. In Figure 2.1 the implicit (1.3) and explicit (1.4) numerical solutions are plotted as function of x with 100 time steps, = 1/400. There is little difference between the two solutions and they are close to the reference solution.
1 Fig. 2.1. Numerical solutions at T = 4 with BuckleyLeverett equation, h = Left picture: explicit method (1.4), right picture: implicit BDF 2 method (1.3).
1 , 100
1 . 400
If the number of time steps is decreased to 50, = 1/200, we see from Figure 2.2 that now the explicit solution becomes unstable, but at the same time the implicit solution becomes very inaccurate. Both the shock speed and the shock height are no longer correct. With linear convection f (u) = u, the same phenomenon was observed: if the solution has steep gradients, then the implicit method gives poor results whenever the step sizes are significantly larger than those that can be taken with the explicit method. As we shall see in the following section, this disappointing qualitative behavior of the implicit BDF2 method is due to loss of monotonicity for large step sizes. Although this can be somewhat improved with variants of the implicit BDF2 method (see next section), tests with other implicit schemes of RungeKutta or linear multistep-type consistently showed a similar behavior. This means that implicit methods can be
1 Fig. 2.2. Numerical solutions at T = 4 with BuckleyLeverett equation, h = Left picture: explicit method (1.4), right picture: implicit BDF 2 method (1.3).
1 , 100
1 . 200
used well only with large Courant numbers if the solution has little temporal or spatial variation. In case this is valid, an implicit treatment will be more efficient than an explicit one. In the following sections we shall consider combinations of the implicit and explicit BDF2 methods with the aim of combining the favorable aspects of these two methods. 3. The -BDF2 methods. As a first step to combine the implicit and explicit methods we consider the following class of methods, with parameter [0, 1]: (3.1)
3 2 wn 1 - 2wn-1 + 2 wn-2 = F (tn , wn + (1 - )wn ),
where as before wn = 2wn-1 - wn-2 . Clearly, for = 0 and = 1 we reobtain the methods (1.3), (1.4), respectively. As we shall see later on, the above methods 3 have order 2 for any choice of . Moreover, the methods are A-stable for 4 and consequently we then have unconditional stability for convection-diffusion problems. 3 In fact, if = 4 the stability region consists precisely of the left half complex plane. With this value of the method has no inherent damping. For diffusion problems the fully implicit BDF2 method with = 1 is therefore preferred. For convection, on the other hand, damping is not necessarily a favorable property and we shall see 3 that = 4 has better monotonicity properties, and consequently it gives a better qualitative behavior for convection problems. 3.1. Positivity properties. We shall consider monotonicity and positivity properties of the -BDF2 method (3.1) for linear equations (3.2) w (t) = Aw(t) + g(t).
In the following we shall write v 0 for a vector v if all its components are nonnegative. It will be assumed that the matrix A = (aij ) Rmm has no real positive eigenvalues and (3.3) aij 0 (for i = j), aii - (for all i),
with > 0. The class of matrices satisfying this condition is denoted by M . By a continuity argument (on > 0) it can be shown that for any A M (3.4) (I - A)-1 0 for all > 0.
PARTIALLY IMPLICIT BDF2 BLENDS
Further we consider g(t) 0 for all t 0 in (3.2). Under these assumptions it holds that (3.5) w(t) 0 whenever t 0 and w(0) 0,
irrespective of the value of R; see . We note that for linear systems w (t) = Aw(t) with the property that Ae = 0 for e = (1, 1, . . . , 1)T , it easily follows that the solution will also satisfy a maximum principle min wi (0) wj (t) max wi (0).
A rational function is said to be absolutely monotonic on the interval [-, 0] if and all its derivatives are nonnegative on this interval. It was shown by Bolley and Crouzeix  that ( A) 0 for all A M iff is absolutely monotonic on [- , 0].
This result gives necessary and sufficient conditions for having wn 0, n = 1, 2, . . . whenever w0 0
with one-step time discretizations, such as RungeKutta methods. The condition of absolute monotonicity is already necessary for A = h-1 (E - I) Rmm with backward shift operator E Rmm , = m = h-1 , provided that the dimension m is sufficiently large. Note that this is simply the semidiscrete system obtained from ut + ux = 0 with first order upwind discretization in space and homogeneous Dirichlet condition at the inflow boundary. In particular, for the one-step -method (1.9), we get the condition on the step size 1 . 1-
Therefore, with the implicit Euler method there is no step size restriction for positivity. With all other well-known methods we do get a restriction on the allowable step sizes, since unconditional positivity implies that the order of the method is at most 1; see . Application of method (3.1) to the linear system (3.2) gives the recursion (3.6) wn = 1 ( A)wn-1 + 2 ( A)wn-2 + ( A) g(tn )
with rational functions (3.7) 1 (z) = 4(1 + (1 - )z) -(1 + 2(1 - )z) 2 , 2 (z) = , (z) = . 3 - 2z 3 - 2z 3 - 2z
Positivity results with arbitrary nonnegative starting values w0 , w1 were derived by Bolley and Crouzeix  for a class of linear multistep methods (see also Spijker  and Shu  for related results). These results, however, require that 1 ( A), 2 ( A), ( A) 0, and therefore they are not applicable to the BDF schemes. Due to the fact that 2 (0) = - 1 one never has w2 0 for arbitrary starting values w0 , w1 0. 3 We shall derive positivity results for the -BDF2 methods (3.1) under the assumption that w1 is obtained by a suitable starting procedure from w0 , for instance, by Euler's method. The derivation of these results is partly based on discussions with M. van Loon (1996, private communications). Results of this type for general multistep methods seem unknown.
3.2. The threshold function. The positivity results will be obtained by considering the above recursion (3.6) with suitable linear combinations wn - wn-1 . In this subsection some technical results will be derived. The final result is given in Theorem 3.1. In the following we denote C(z) = Then V C(z)V -1 = with 1 (z) = 1 (z) - , 2 (z) = 1 (z) + 2 (z) -
1 (z) 1
2 (z) 0
1 (z) 1
We shall determine > 0 such that the entries of V C(z)V -1 are absolutely monotonic on the interval [-, 0] with as large as possible. Since the j are fractional linear (i.e., rational with linear denominator and numerator), it follows that this is equivalent to j (0) 0 and j (z) 0 for z [-, 0], j = 1, 2. It is straightforward to verify that j (0) 0 and j (0) 0 for j = 1, 2 iff (3.8)
1 3 - 2 , 3 6 - 2
Further we want j (z) 0. As we consider z 0, this is seen to be equivalent with (3.9) where r( ) = 2 + 4(1 - ) p( ) = (1 - )(3 - 1), q( ) = 2
|z| r( ),
q( )|z| p( ), 4-3 ,
+ 4(1 - ) - 2(1 - ).
The optimal choice for will depend on the location of the largest zero 2 of q( ). We have 1 - 1 1,2 = - q( ) = 2( - 1 )( - 2 ), 1 - . Note that r( ) is monotonically decreasing in , and to satisfy |z| r( ) for z [-, 0] with as large as possible, we should take [ 0 , 1] as small as possible, but of course within the second constraint of (3.9). 1 First, assume that 2 3 , that is, 3 . Then q( ) 0 for [ 1 , 2 ], and thus 4 3 the second constraint in (3.9) will be automatically satisfied for these . Therefore we can choose = 0 , yielding the restriction r( 0 ). Thus the optimal is given by (3.10) () = 15 - 2 , 24 - 26 + 42
3 For the second case 2 < 1 , that is, 4 , we get the condition 3
1max min r( ),
p( ) q( )
PARTIALLY IMPLICIT BDF2 BLENDS
By some tedious calculations it can be shown that the second constraint is now the dominating one and that the above condition is restrictive least with = [(3 - 2) + (4 - 3)]/(6 - 2). This leads to the optimal given by 3 + 2 - 3 4 - 3 () = (3.11) , > 3. 4 2(6 - 5) + 2 4 - 3 The threshold function () from (3.10), (3.11) is plotted in Figure 3.1. In the next subsections the relevance of this function is discussed.
Fig. 3.1. Positivity threshold function () versus [0, 1] according to (3.10), (3.11).
3.3. Results for linear systems. From the calculations in section 3.2 it is easy to obtain positivity results for linear systems. In the following, () refers to the threshold function given by (3.10), (3.11) and = () stands for the optimal value such that V C(z)V -1 0 for all z [-(), 0]. Theorem 3.1. Consider the linear semidiscrete system (3.2) with A M and g(t) 0. Then wn 0 whenever (), w0 0, and w1 - w0 0. Proof. Denote Wn = wn wn-1 , C( A) = 1 ( A) I 2 ( A) O , Gn = ( A)g(tn ) . 0
Recursion (3.6) can be written as Wn = C( A)Wn-1 + Gn . We consider Un = V Wn Then Un = V C( A)V -1 Un-1 + Gn . with V = I O - I I .
From the results in section 3.2 it follows that the entries of the block matrix V C( A)V -1 are nonnegative provided that (). Further we have Gn 0 and U0 0. Therefore Un 0 for all n, and consequently the same holds for the Wn . Whether the condition w1 - w0 0 is satisfied will of course depend on the starting procedure used to calculate w1 . It will hold if w1 is calculated from one implicit Euler step. However, if = 0 it is more natural to use an explicit Euler step. Since = 1 if = 0, we then get 2
1 w1 - w0 = w1 - 2 w0 = 1 w0 + Aw0 + g(0), 2 1 and this is guaranteed to be nonnegative only if 2 . This condition is slightly 5 more restrictive than with the threshold value (0) = 8 for the explicit BDF2 method itself. This extra time step restriction due to the explicit Euler start can be easily avoided by calculating w1 by another starting procedure; for example, w1 = w0 + F (t0 , w0 ), w2 = w1 + F (t1 , w1 ), 1 w1 = 2 (w0 + w2 ),
in which case it is seen that w1 - 1 w0 = 1 w2 0 whenever 1. 2 2
3.4. Test with the van Leer limiter. The above theoretical results give sufficient conditions for nonnegative solutions with linear problems. To test the relevance with the nonlinear semidiscrete systems obtained with limited spatial discretization 1 (2.2)(2.4), we consider once more the 1D test equation ut + ux = 0, 0 t 4 with a 1 block-function as initial profile and h = 100 . In Table 3.1 we have listed the minimal number of steps N (r) needed to obtain numerical solutions with minimum larger than -10-r with r = 3, 4. As before, the convergence criterion in the Newton iteration was that the max-norm of the residual should be less than 10-6 (same results with smaller tolerances), and the starting value w1 was computed with the implicit Euler method.
Table 3.1 Linear convection test (2.1), (2.6) with -BDF 2 methods. Number of steps required for (almost) nonnegative solutions. The Courant numbers are /h = 25/N .
N (4) N (3)
0 40 39
.7 21 21
.74 21 21
.75 24 23
.76 31 26
.8 46 38
1 75 63
For the larger values of the number of steps needed to achieve minimal values larger than -10-4 and -10-3 are relatively far apart; we do not have an explanation for this. We see from Table 3.1 that the theoretical results obtained for the linear class of problems do have a relevance for the van Leer limiter. In particular, if is close to 0.75, we can take significantly larger steps than with equal to 0 or 1. On the other hand, in this test the largest step sizes could be taken with values of slightly less than 0.75, in contrast to Figure 3.1. Also, the allowable step size with = 0 seems somewhat larger than one would expect on the basis of Figure 3.1 in comparison with equal to 0.75 or 1.
PARTIALLY IMPLICIT BDF2 BLENDS
It should be noted that the semidiscrete system obtained here with limiting can be written in the quasi-linear form wi = 1 q ai (w)(wi-1 - wi ) with h 0 ai (w) 2;
see . The results for the linear systems therefore suggest positivity if the Courant numbers = q /h are not larger than 1 (). In the above experiment this condition 2 indeed seems sufficient, but it also seems a bit too strict, probably due to the fact that the limiter switches locally to first order upwind discretization for which the condition () is sufficient (and necessary). Similar to  for explicit Runge Kutta methods, we can conclude that the linear theory does give reasonable qualitative predictions for more difficult, nonlinear situations, but these predictions should not be taken too literally. As noted before, the -BDF2 methods are unconditionally stable for convectiondiffusion problems iff 3 . Based on the linear theory and practical experience, 4 we do prefer the implicit method with = 3 over the standard fully implicit BDF2 4 method with = 1 for convection. For instance, with the 1D BuckleyLeverett 3 test problem (2.7) the choice = 4 still gives accurate results with = 1/200, h = 1/100 for which the standard BDF2 method produces qualitatively poor results; see Figure 2.2. Note, however, that basically we still have the same problems as with = 1, namely, the high cost of solving the implicit relations and the fact that large Courant numbers lead to loss of monotonicity. Therefore we would like to apply this method with = 3 only if the temporal or spatial variation in the solution is not too 4 large. This will be achieved by considering different values for in different parts of the spatial domain. 4. The -BDF2 scheme. To combine implicit and explicit formulas we shall allow to vary over the spatial grid. Let in the following = diag(i ), where i will correspond with grid point xi . We consider once more (1.5) but now with specification of i , 1 3 2 wn - 2wn-1 + 2 wn-2 = F tn , wn + (I - )(2wn-1 - wn-2 ) , (4.1) 0 if i , i = otherwise, with i denoting the local Courant number at grid point xi . We choose = 3 since 4 this appeared the best choice to aim for with respect to stability and positivity, and = 1 since the explicit scheme appears to be stable and positive for i 1 . 2 2 Note that for 1D problems (2.1) the local Courant number is given by i = |q(xi )f (wi )|/hi , where hi is the length of the cell i around xi . For multidimensional problems on Cartesian grids, i is taken as the sum of the 1D contributions. When implemented with variable time steps the matrix will also become variable in time even for linear convection with constant velocities. In section 6 we shall consider a simple variable step size selection procedure that essentially limits the max-norm of the displacement wn - wn-1 . As a consequence, the scheme will be implicit only in those spatial regions where the velocities are large, but the solution is smooth. With the above choice for we apply the explicit scheme as much as possible 3 within the stability constraint, and we switch to = 4 elsewhere. With this choice there are abrupt changes in the values of the i over the grid. The effect of this on the accuracy is discussed next.
First we take a look at the truncation error of (4.1). Let w(tn ) = 2w(tn-1 ) - w(tn-2 ). Insertion of the exact solution of (1.2) into the scheme gives (4.2)
3 2 w(tn ) 1 - 2w(tn-1 ) + 2 w(tn-2 ) = F tn , w(tn ) + (I - )w(tn ) + rn 1 - 2w(tn-1 ) + 1 w(tn-2 ) = w (tn ) - 3 3 w (tn ) + O( 4 ), 2
with truncation error rn . By a Taylor expansion we obtain
3 2 w(tn )
w(tn ) + (I - )w(tn ) = w(tn ) - (I - ) 2 w (tn ) + O( 3 ), and hence (4.3)
with Jacobian matrix An = w F (tn , w(tn )). If a diffusion term is added as in (1.8) this formula for the truncation error is still valid. The truncation error is often a good measure of the accuracy. Indeed, if we are dealing with a fixed ODE system, then the truncation error is O( 2 ), reflecting the second order accuracy of the formula. However, in our situation where the ODE system is a semidiscrete PDE, the function F and its derivatives will contain negative powers of the mesh width h. In particular, the term 2 An (I - )w (tn ) in (4.3) will be only a genuine O( 2 ) term if is sufficiently smooth in space. With the choice (4.1) this does not hold. Yet, as we shall see, the accuracy is not affected by this. Instead of looking only at the truncation error, a more refined error analysis is needed. This will be presented in the next section for linear systems. We note that in (4.1) the linear combination with is taken "within" the function F to ensure mass conservation. The related method 3 1 2 wn - 2wn-1 + 2 wn-2 = F (tn , wn ) (4.4) + 2(I - )F (tn-1 , wn-1 ) - (I - )F (tn-2 , wn-2 ) 1 rn = - 3 2 w (tn ) + 2 An (I - )w (tn ) + O( 3 )
has smaller truncation errors in general. By Taylor expansion it is easily seen that the truncation error of (4.4) is equal to
1 3 2 w(tn ) 1 - 2w(tn-1 ) + 2 w(tn-2 ) - w (tn ) - 2(I - )w (tn-1 ) 2 3I
+ (I - )w (tn-2 ) = 2
- w (tn ) + O( 3 ).
Therefore, as far as local accuracy is concerned, the form (4.4) is better than (4.1) in general. This is similar to genuine multistep formulas versus the so-called one-leg formulations; see . However, the form (4.4) is not mass conserving. Suppose that the discrete mass is given by T w(t) = i wi (t) with components i denoting the length of grid cell i , or area or volume in more dimensions; then mass conservation of the semidiscrete system (1.2) means that T w(t) should remain constant in time for all starting values w(0). This is equivalent to the condition T F (t, w) = 0 for all t, w. Now, suppose that T w0 = T w1 . Then with (4.1) it easily follows by induction that we will have T wn = T w0 for all n. With formula (4.4), however, this will hold only if = I, that is, constant over the space. Therefore, even though (4.4) has smaller truncation errors in general, we shall continue with the form (4.1).
PARTIALLY IMPLICIT BDF2 BLENDS
5. Global accuracy results. In this section an error analysis for the -BDF2 scheme (4.1) will be presented for linear systems (5.1) w (t) = Aw(t) + g(t),
where the matrix A is assumed to be a finite difference approximation to a convective operator. Stability results with a that varies over the space according to (4.1) are not available. The variation in over space has as a consequence that the standard von Neumann analysis, based on Fourier decompositions, is no longer applicable. In the numerical tests the scheme (4.1) never encountered stability problems. In the following it will therefore simply be assumed that the scheme is stable in a given norm for the above linear system, and we will consider global accuracy of the scheme under this assumption. Let n = w(tn ) - wn be the global discretization error. From (1.5) and (4.2) we obtain the error recursion (5.2)
1 2 4 n - 3 n-1 + 3 n-2 = 2 Z(n + (I - )(2n-1 - n-2 )) + 3 rn , 3
where Z = A and rn is the local truncation error. This can be written in the more transparent form (5.3) with matrices 1 =
4 3 2 I - 3 Z -1
n = 1 n-1 + 2 n-2 + n
I + Z(I - ) ,
2 = - 1 I - 2 Z 3 3
I + 2Z(I - )
determining the propagation of previous errors, and with n the local discretization error introduced in the step from tn-1 to tn ,
2 n = I - 3 Z -1 2 3 rn .
For the linear system (5.1) this local discretization error equals (5.4)
2 2 n = (I - 2 Z)-1 (- 9 3 w (tn ) + 3 2 Z(I - )w (tn )) + O( 3 ). 3
Here the last term contains only genuine O( 3 ) terms; there are no hidden negative powers of h in the constant. Our tacit stability assumption can now be specified: we assume that from the error recursion (5.3) it can be concluded that
with C > 0 a moderate stability constant, independent of the mesh width h. In particular, this assumption implies that 1 and 2 are bounded, from which it 2 easily follows that terms like (I - 2 Z)-1 and (I - 3 Z)-1 Z are also bounded 3 (by moderate constants, independent of h). It thus follows from (5.4) that n = O( 2 ). Note that this deviates from the estimate that would be obtained in the standard ODE case with a fixed, bounded matrix A. In that case Z = O( ) and consequently n = O( 3 ).
n C 0 + 1 +
Since we are dealing with semidiscrete systems arising from PDEs, where A will contain negative powers of h, the local error n is merely O( 2 ) in general. Thus one might expect the global errors to be first order only. However, similar to  and [10, sect.V.7] for stiff ODEs, it will be shown here that due to cancellation and damping effects we still have global convergence with order 2. To demonstrate this second order convergence, define (5.6) = n + 2 (I - )w (tn ), n
which will turn out to behave more regular than n . By observing that I - 1 - 2 = - 2 I - 2 Z 3 3
it follows that these transformed errors satisfy the recursion n
= 1 + 2 + n n n-1 n-2
with transformed local error
n = n - 2 (I - )w (tn ) + 1 2 (I - )w (tn-1 ) 2 + 2 2 (I - )w (tn-2 ) = - I - 3 Z -1 2 3 9 w
(tn ) + O( 3 ).
This transformed local error is genuinely of order 3, independent of the mesh width h. The stability argument applied to the recursion of the transformed errors now yields in a standard way order 2 convergence for the . Hence it follows that we also have n for our original errors n = O( 2 ), uniformly for tn T , independent of the mesh width h. Although this is not a complete convergence proof, since we had to assume that the scheme is stable, it does show that the choice for in (4.1), with abrupt changes in i over the grid, will not lead to an order reduction. Remark. The above analysis carries over to linear systems w (t) = Aw(t) + Bw(t) + g(t), where B is a diffusion term that is treated fully implicitly as in (1.8). The transformed errors should then be defined as = n + 2 X(I - )w (tn ) n with X = (A + B)-1 A. We then obtain second order convergence provided that X = O(1) uniformly in h. 6. Numerical results. In this section numerical results are presented for a 2D test convection problem arising from the quarter of five spots problem in reservoir simulations; see [7, 18], for example. On a square region = [0, 1]2 we have a source 1 term at the point x = (0, 0), with volumetric rate = 4 , and a sink term - at x = (1, 1), corresponding to an injection and production well, respectively. It is assumed here that the permeability K and viscosity in the actual reservoir problem are constant, say, K/ = 1. The velocity q and pressure p are then given by (6.1) q = - p, p + s = 0,
PARTIALLY IMPLICIT BDF2 BLENDS
with s = s(x) describing the sources and sinks and with homogeneous Neumann boundary conditions for p. This determines p up to an additive constant. The resulting convection problem is (6.2) ut + .(qf (u)) = s+ + s- u,
where s+ = max(s, 0) and s- = min(s, 0). The initial condition is u 0. For the flux function f we shall consider both the linear flux function (2.6) and the Buckley Leverett flux function (2.7). These are simplified model situations for miscible and immiscible reservoir flows. Illustrations for the behavior of the solutions on = [0, 1]2 are given in Figure 6.1.
T = 1/2 , miscible 1 0.8 0.6 0.4 0.2 0.8 0.6 0.4 0.2 0 0.8 0.6 0.4 0.2 T = 1 , miscible 1 0.8 0.6 0.4 0.2 0
T = 1/2 , immiscible
T = 1 , immiscible
0.8 0.6 0.4 0.2
0.8 0.6 0.4 0.2 0
0.8 0.6 0.4 0.2
0.8 0.6 0.4 0.2 0
Fig. 6.1. Numerical solutions at T = 1 and T = 1 on 50 50 grids for the miscible model 2 with linear convection (top pictures) and the immiscible model with Buckley-Leverett fluxes (bottom pictures).
In the numerical tests, the pressure equation was solved using standard second order finite differences on a uniform m m grid, mesh width h = 1/m, resulting in a first order approximation of the velocities at the cell edges. The injection well was modelled as a source term /h2 in the lower left grid block. Likewise, for the production well we get a sink term -wm,m /h2 at the upper right grid block. For real reservoir simulations the pressure equations are usually solved in a more sophisticated manner; see, for instance, the contribution of Russell and Wheeler in . With the above test problem the pressure could even be calculated analytically, but numerical solution directly leads to approximations for the velocities that are divergence-free in a discrete fashion. The convection terms in (6.2) are discretized on the same uniform grid with the van Leer limiter as described in section 2; see also Molenaar .
The velocities are large only at the corners where the wells are located, approx1 imately 2 log r with distance r to the well near (0, 0) and (1, 1), respectively. Due to the injection at x = (0, 0) a front has formed at t = 0, which is roughly halfway to the production well at time t = 1 ; see Figure 6.1. Therefore, in the vicinity of the 2 sharp front we could then use the explicit BDF2 method. Near the wells the solution is smooth, so that there an implicit method could easily be applied. A combination of this is provided by the blended scheme (4.1). The time integrations in the numerical tests were started with a small initial time 1 step 0 = 100 h2 and subsequently a simple variable step size selection was used, (6.3) tn+1 = tn + n , n = n-1 , = min(2, tol un
un - un-1
The variable step size form of the -BDF2 method was taken as (6.4) (1 + 2)wn+1 - (1 + )2 wn + 2 wn-1 = (1 + )n F n wn+1 + (I - n )((1 + )wn - wn-1 ) , where the coefficients are similar to the standard implicit BDF2 method; see , for example. The initial step is taken with the Euler method, implicit if > 0 and explicit if = 0. We note that the step size selection used here is the same as in . Results with a more sophisticated selection, based on an estimate of higher derivatives, gave comparable results. Since the focus here is on the methods and not on step size selections, only the results for the above implementation are presented. The implicit relations were solved with a modified Newton iteration, using first order upwind discretizations for Jacobian approximations, as described in section 2. In the Newton iteration the initial guess for wn+1 in (6.4) was taken as n wn+1 + (I - n ) (1 + )2 2 1+ wn - wn-1 + n F (tn+1 , wn+1 ) . 1 + 2 1 + 2 1 + 2
To solve the arising linear systems the Bi-CGSTAB method  was used without preconditioning. Note that due to the first order upwind approximation in the Newton iteration the linear system is diagonally dominant. This choice for the linear solver was guided by experiments in , where several linear solvers were compared for more general porous media equations. Both the Newton iteration and the Bi-CGSTAB iteration were stopped as soon as the norm of the residue was below 10-6 . The norm used here is the maximum norm, as in the step size selection, instead of the more common weighted L2 -norm as in , since we also want to resolve the steep solutions gradients accurately. 1 In Tables 6.1 and 6.2 the statistics are presented for output time T = 2 with the implicit, explicit, and blended scheme (4.1). Along with a CPU timing in seconds on a SUN sparc4 workstation, also given are the average number of Newton iterations per step (N-it) and the average number of Bi-CGSTAB iterations per Newton iteration (L-it). In the step size selection we used tol = 0.1 for the implicit and partially implicit scheme, and tol = 0.01 for the explicit scheme. With the explicit scheme this smaller value of tol was needed to avoid oscillations (mild instabilities) near the inflow well. With this choice, the accuracy of the various schemes was very similar; the spatial discretization errors are the dominating ones. Since the errors of the three methods were similar in the experiments, the CPU time is a measure of efficiency here. Obviously this is most favorable with the blended
PARTIALLY IMPLICIT BDF2 BLENDS Table 6.1 Statistics for 2D linear convection at T = 1 on 50 50 and 100 100 grid. 2 Implicit Blended Explicit Implicit Blended Explicit 1 .75 0 1 .75 0 0 .5 0 0 .5 0 tol .1 .1 .01 .1 .1 .01 Grid 50 50 50 50 50 50 100 100 100 100 100 100 Steps 218 226 2142 340 364 4016 cpu (s) 217 44 131 2205 413 963 N-it 3.34 0.25 3.92 0.51 L-it 2.52 1.14 4.19 2.37 -
Table 6.2 Statistics for 2D BuckleyLeverett at T = 1 on 50 50 and 100 100 grid. 2 Implicit Blended Explicit Implicit Blended Explicit 1 .75 0 1 .75 0 0 .5 0 0 .5 0 tol .1 .1 .01 .1 .1 .01 Grid 50 50 50 50 50 50 100 100 100 100 100 100 Steps 292 280 2985 531 498 5515 cpu (s) 288 65 227 2318 445 1603 N-it 3.57 0.21 3.90 0.24 L-it 1.55 1.00 1.60 0.99 -
method. It should be noted, however, that the explicit scheme also performs quite well. With the step size selection described above, the maximal Courant numbers are much larger than unity without introducing instabilities. There are still some small oscillations with the explicit method near the inflow corner, but on the scale of Figure 6.1 these are not visible. Apparently, relatively large Courant numbers can be taken here with the explicit scheme since the velocities are large only near the wells and possible instabilities are transported to the interior domain where they are damped. However, the step sizes that can be taken with the implicit and blended scheme are much larger, but the fully implicit scheme is not efficient due to the amount of work that has to be performed in solving the algebraic relations. The blended scheme is initially fully explicit, since the step sizes selected according to (6.3) are small if the sharp front is in a region with large velocities. After awhile this scheme becomes implicit near the wells, but then the implicit relations are easy to solve since the solution does not vary much anymore near the wells. It should be noted that the performance of the explicit scheme will decrease if a local grid refinement is used near the wells. This is often done in practice to capture small-scale geological features. In such a situation a more pronounced advantage of the blended scheme can be expected. This has not been tested, since for the present model problem such a grid refinement would be very artificial. Numerical tests with small diffusion terms added to the convection equation, implemented as in (1.8), did give very similar results. Finally it should be noted that our implementation of the blended scheme in the above experiments was not very sophisticated. For example, the whole function F was calculated in each Newton iteration step, whereas this is not necessary inside the region where = 0 (more precisely, at those grid points where i = 0 for the grid point itself, its neighbors and
adjacent points). For ease of programming it was decided to use the same subroutines as for the fully implicit scheme. In view of these experiments, we conclude that the blended scheme works very well for problems of the above type, where there are locally large velocities. If the size of the velocities is more or less uniform, and the solution is not very smooth, an explicit treatment of the convective terms will be more efficient in general. Fully implicit methods seem to be efficient only if the solution is sufficiently smooth in space, but with convection dominated flows steep gradients in the solution are the generic case.
 U. M. Ascher, S. J. Ruuth, and B. T. R. Wetton, Implicit-explicit methods for timedependent partial differential equations, SIAM J. Numer. Anal., 32 (1995), pp. 797823.  C. Bolley and M. Crouzeix, Conservation de la positivit lors de la discrtisation des e e probl`mes d'volution paraboliques, RAIRO Anal. Numer., 12 (1978), pp. 237245. e e  J. G. Blom, J. G. Verwer, and R. A. Trompert, A comparison between direct and iterative methods to solve the linear systems arising from a time-dependent 2D groundwater flow model, Comp. Fluid Dyn., 1 (1993), pp. 95113.  M. Blunt and B. Rubin, Implicit flux limiting schemes for petroleum reservoir simulation, J. Comput. Phys., 102 (1992), pp. 194210.  M. Crouzeix, Une mthode multipas implicite-explicite pour l' approximation des quations e e d' volution paraboliques, Numer. Math., 35 (1980), pp. 257276. e  J.-A. Desideri and P. W. Hemker, Convergence analysis of the defect-correction iteration for hyperbolic problems, SIAM J. Sci. Comput., 16 (1995), pp. 88118.  R. E. Ewing, ed., The Mathematics of Reservoir Simulation, Frontiers Appl. Math. 1, SIAM, Philadelphia, 1984.  J. Frank, W. Hundsdorfer, and J. G. Verwer, On the stability of implicit-explicit linear multistep methods, Appl. Numer. Math., 25 (1997), pp. 193205.  E. Hairer, S. P. Nrsett, and G. Wanner, Solving Ordinary Differential Equations I-- Nonstiff Problems, Springer-Verlag, Berlin, 1987.  E. Hairer and G. Wanner, Solving Ordinary Differential Equations II--Stiff and DifferentialAlgebraic Problems, Springer-Verlag, Berlin, 1991.  W. Hundsdorfer, B. Koren, M. van Loon, and J. G. Verwer, A positive finite-difference advection scheme, J. Comput. Phys., 117 (1995), pp. 3546.  W. Hundsdorfer and B. I. Steininger, Convergence of linear multistep and one-leg methods for stiff nonlinear initial value problems, BIT, 31 (1991), pp. 124143.  W. Hundsdorfer and R. Trompert, Method of lines and direct discretization: A comparison for linear advection, Appl. Numer. Math., 13 (1994), pp. 469490.  B. Koren, A robust upwind discretization for advection, diffusion and source terms, in Numerical Methods for Advection-Diffusion Problems, C. B. Vreugdenhil and B. Koren, eds., Notes Numer. Fluid Mech. 45, Vieweg, Braunschweig, 1993.  D. Kroner, Numerical Schemes for Conservation Laws, Wiley and Teubner, Chichester, UK, Stuttgart, 1997.  B. van Leer, Towards the ultimate conservative difference scheme II. Monotonicity and conservation combined in a second order scheme, J. Comput. Phys., 14 (1974), pp. 361370.  R. J. LeVeque, Numerical methods for conservation laws, Lectures Math. ETH Zrich, u Birkhuser-Verlag, Basel, 1992. a  J. Molenaar, Multigrid methods for high-order accurate fully implicit simulations of flow in porous media, in Proceedings of the Second ECCOMAS Conference on Numerical Methods in Engineering, J.-A. Dsidri, P. Le Tallec, E. Onate, J. Priaux, and E. Stein, eds., John e e e Wiley, 1996.  M. N. Spijker, Contractivity in the numerical solution of initial boundary value problems, Numer. Math., 42 (1983), pp. 271290.  C.-W. Shu, Total-variation-diminishing time discretizations, SIAM J. Sci. Statist. Comput., 9 (1988), pp. 10731084.
PARTIALLY IMPLICIT BDF2 BLENDS
 J. M. Varah, Stability restrictions on second order, three level finite difference schemes for parabolic equations, SIAM J. Numer. Anal., 17 (1980), pp. 300309.  J. G. Verwer, J. G. Blom, and W. Hundsdorfer, An implicit-explicit approach for atmospheric transport-chemistry problems, Appl. Numer. Math., 20 (1996), pp. 191209.  H. A. van der Vorst, Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems, SIAM J. Sci. Statist. Comput., 13 (1992), pp. 631644.