A Second Order Energy Stable Scheme for the Cahn-Hilliard-Hele-Shaw Equations

We present a second-order-in-time finite difference scheme for the Cahn-Hilliard-Hele-Shaw equations. This numerical method is uniquely solvable and unconditionally energy stable. At each time step, this scheme leads to a system of nonlinear equations that can be efficiently solved by a nonlinear multigrid solver. Owing to the energy stability, we derive an $\ell^2 (0,T; H_h^3)$ stability of the numerical scheme. To overcome the difficulty associated with the convection term $\nabla \cdot (\phi \boldsymbol{u})$, we perform an $\ell^\infty (0,T; H_h^1)$ error estimate instead of the classical $\ell^\infty (0,T; \ell^2)$ one to obtain the optimal rate convergence analysis. In addition, various numerical simulations are carried out, which demonstrate the accuracy and efficiency of the proposed numerical scheme.


Introduction
The Cahn-Hilliard-Hele-Shaw (CHHS) diffuse interface model has attracted a lot of attention because it describes two phase flows in a simple way [17,18].It has been used to model spinodal decomposition of a binary fluid in a Hele-Shaw cell [14], tumor growth and cell sorting [11,24], and two phase flows in porous media [4].It describes the process of the phase separation of a viscous, binary fluid into domains.In this model, the Cahn-Hilliard (CH) energy of a binary fluid with a constant mass density is given by [2]: where Ω ⊂ R d (d = 2 or 3), φ : Ω → R is the concentration field, and ε is a constant.The phase equilibria are represented by the pure fluids φ = ±1.For simplicity, we assume that Ω = (0, L x ) × (0, L y ) × (0, L z ) and that ∂ n φ = 0 on ∂Ω, the latter condition representing local thermodynamic equilibrium on the boundary.The dynamic equations of CHHS model [17,18] are given by where γ > 0 is related to surface tension and the chemical potential is defined as u is the advective velocity; and p is the pressure.We assume no flux boundary condition, namely u • n = 0 and ∂ n µ = 0, with n the unit normal vector on ∂Ω: ∂φ ∂n = ∂µ ∂n = ∂p ∂n = 0 on ∂Ω T := ∂Ω × (0, T ], (1.6) The system (1.2)-(1.4) is mass conservative and energy dissipative, and the dissipation rate is readily found to be Another fundamental observation is that the energy (1.1) admits a splitting into purely convex and concave parts, i.e., E = E c − E e : where both E c and E e are convex.Based on this observation, a first order in time unconditionally energy stable finite difference for the CHHS equations was proposed in [23], and the detailed convergence analysis has become available in a more recent work [3].Meanwhile, Feng and Wise presented finite element analysis for the system (1.2)- (1.4), which arise as a diffuse interface model for the two phase Hele-Shaw flow in [10].Collins et al. proposed an unconditionally energy stable and uniquely solvable finite difference scheme for the Cahn-Hilliard-Brinkman (CHB) system, which is comprised of a CH-type diffusion equation and a generalized Brinkman equation modeling fluid flow.The detailed convergence analysis for the first order convex splitting scheme to the Cahn-Hilliard-Stokes (CHS) equation was provided in [5].In [13], Guo et al. presented an energy stable fully-discrete local discontinuous Galerkin (LDG) method for the CHHS equations.And also, Han proposed and analyzed a decoupled unconditionally stable numerical scheme for the CHHS equations with variable viscosity in [14].
Most of the existing schemes are of first oder accuracy in time.In this paper, we propose and analyze a second order convex splitting scheme for the system (1.2)- (1.4), which turns out to be uniquely solvable and unconditionally energy stable.A modified Crank-Nicholson approximation is applied to the nonlinear part of the chemical potential, an explicit Adams-Bashforth extrapolation is applied to the concave term, and an Adams-Moulton interpolation formula is applied to the highest order surface diffusion term.In more details, such an Adams-Moulton interpolation formula is applied at the time steps t n+1 and t n−1 (instead of the standard one at t n+1 and t n ), so that the diffusion coefficient at t n+1 dominates the others.This subtle fact will greatly facilitates the convergence analysis; see the related works for the pure CH flow: [12] with the finite difference spatial approximation, [7] with the finite element version.In addition, a semi-implicit approximation is applied to the nonlinear convection term, with the phase variable treated via extrapolation and the velocity field implicitly determined by the Darcy law at the numerical level.A careful analysis reveals a rewritten form of the numerical scheme as the gradient of strictly convex functional, so that both the unique solvability and unconditional energy stability could be theoretically justified.
Meanwhile, it is noted that an optimal rate convergence analysis for the second order scheme to the CHHS equation remains open.The key difficulty is associated with the high degree of nonlinearity of the convection term, ∇•(φu), with u the Helmholtz projection of −γφ∇µ.And also, the Darcy law in the fluid equation has also posed a serious challenge in the numerical analysis, in comparison with the CHS [5] or Cahn-Hilliard-Navier-Stokes (CHNS) model [6], in which a kinematic diffusion is available for the fluid.For the CHHS equation, even the highest order linear diffusion term could not directly control the error estimates for the nonlinear terms, due to the nonlinear convection.For the first order numerical scheme, the methodology to overcome such a difficulty was reported in a few recent works [3,19].However, these analysis techniques could not be directly applied to the second order scheme.In this article, we present a detailed analysis to establish the full order convergence of the proposed numerical scheme, with second order accuracy in both time and space.In more details, a nonlinear energy estimate by taking an inner product with µ k+1/2 , the numerical chemical potential at time instant t k+1/2 , gives an unconditional numerical stability.Moreover, a more careful analysis for the chemical potential gradient, in combination with the Sobolev inequalities at the discrete level, leads to an 2 (0, T ; H 3 h ) stability estimate of the numerical solution.On the other hand, a subtle observation indicates that the estimate for the nonlinear error associated with ∇ • (φu) cannot be carried out in a standard way, due to a broken structure for this nonlinear error function.As a result, an ∞ (0, T ; H 1 h ) error estimate has to be performed, instead of the classical ∞ (0, T ; 2 ) one, since the error term associated with the nonlinear convection has a non-positive inner product with the appropriate error test function.
In addition, the 2 (0, T ; H 3 h ) bound of the numerical solution plays a key role in the nonlinear error estimate, which enables us to apply the discrete Gronwall inequality to obtain the desired convergence result.
The remainder of this paper is organized as follows.In Section 2, we present the fully-discrete scheme for CHHS equations.The 2 (0, T ; H 3 h ) stability of the numerical scheme is further established in Section 3. In Section 4, we present the optimal rate convergence analysis with the help ∞ (0, T ; H 1 h ) error estimate.In Section 5, we provide some numerical results to validate our theoretical analysis and demonstrate the effectiveness of the proposed fully discrete finite difference method.To solve the nonlinear equations at each time step, the nonlinear multigrid solver is applied.Finally, we offer our concluding remarks in Section 6.

The fully discrete scheme and a-priori stabilities
In this section, we propose a second order in time fully discrete scheme for the system (1.2)- (1.4) with the discrete homogeneous Neumann boundary conditions (1.6).For simplicity, we consider M > 0 for some M ∈ N , be the time step size and t m = ms.We only consider the three-dimensional version of the fully discrete scheme for the CHHS system since an extension to the two-dimensional case is trivial.For convenience, some of the following notations are defined in Appendix A. For each integer m, 1 ≤ m ≤ M − 1, given on ∂Ω, where for any (ϕ, ψ) ∈ [C Ω ] 2 .Note that the three component variables of the velocity vector u m+1/2 ∈ E Ω are located at the staggered grid points.To facilitate the unique solvability analysis below, we could eliminate the velocity variable in the numerical scheme and rephrase it in terms of (φ m+1 , µ m+1/2 , p m+1/2 ) ∈ [C Ω ] 3 .In more details, we introduce M(φ) := 1 + γφ 2 and rewrite (2.1)- The symbol M(A h φ m+1/2 * )∇ h µ m+1/2 represents a discrete vector field.For instance, the y-component at a generic y-face grid point is given as The definitions of the discrete operators used above can be found in Appendix A.2 and are similar to those found in [23].
We now define a fully discrete energy that is consistent with the continuous space energy (1.1) as h → 0. In particular, the discrete energy We also define an alternate numerical energy via (2.9) We can not guarantee that the energy E h is non-increasing in time, but, we can guarantee the dissipation of auxiliary energy F h .
For our present and future use, we define the canonical grid projection operator and s → 0 for sufficiently regular u.The next theorem addresses the unique solvability and unconditional energy stability of the numerical solutions to the scheme (2.5) -(2.7): Theorem 2.1.Suppose that (φ e , µ e , u e ) is the sufficiently regular exact solution to the CHHS system (1.2)-(1.4).Take Φ i,j,k = P h φ e (•, t ) and suppose that the initial profile , there is a unique solution φ m+1 ∈ C Ω to the scheme (2.5) - (2.7).And also, the scheme (2.5) -(2.7), with starting values φ 0 and φ 1 , is unconditionally energy stable with respect to (2.9), i.e., for any s > 0 and h > 0, and any positive integer where C 0 > is a constant independent of s, h, and , and u m+1/2 ∈ E Ω is given by (2.3).
Proof.The unique solvability proof follows from the convexity analysis, as presented in [23] for the first order convex splitting scheme applied to the CHHS equation.We define a linear operator L as in which Following the arguments in [23], we are able to prove that L gives rise to a symmetric, coercive, and continuous bilinear form when the domain is restricted to CΩ := {µ ∈ C Ω : (µ, 1) = 0}; the details are skipped for the sake of brevity and left to interested readers.
Subsequently, an inner product on CΩ is introduced using L: let f µ and f ν ∈ CΩ and suppose µ, ν ∈ H1 h are the unique solutions to L(µ) = f µ and L(ν) = f ν .Then we define It is straightforward to verify that Next, we consider the following functional: with ) (2.17) The convexity of F c follows from the convexity of g c (φ) := 1 16 φ 4 + 1 12 φ m φ 3 + 1 8 (φ m ) 2 φ 2 (in terms of φ).And also, (•, •) L −1 is an inner product.Therefore, we conclude that G is convex.Moreover, G is coercive over the set of admissible functions Therefore it has a unique minimizer, and in particular the minimizer of G, which we denote as φ = φ m+1 , satisfied the discrete equation in which C is a constant and φ m+1 satisfies the homogeneous Neumann boundary condition.In other words, φ m+1 is a solution of which is equivalent to the numerical scheme (2.5) -(2.7).The proof of unique solvability is complete.
For the energy stability analysis, we look at the numerical scheme in the original formulation (2.1)- (2.3).Taking an inner product with µ m+1/2 (given by (2.6)) by (2.1) yields In more detail, the leading term has the expansion The estimate for I 1 , a convex term, is straightforward: For the second term I 2 of (2.22), a concave term, we see that where the Cauchy inequality was utilized in the last step.
The third term I 3 of (2.22), also a convex term, can be analyzed with the help of summation by parts: (2.25) The evaluation of I 3,1 is straightforward: (2.26) The estimate of the I 3,2 can be carried out in the following way: in which the Cauchy inequality was applied in the last step.Consequently, substituting (2.26) and (2.27) into (2.25)yields Finally, a combination of (2.22), (2.23), (2.24) and (2.28) results in (2.29) For the second term of (2.21), the boundary condition n•∇ h µ m+1/2 | ∂Ω = 0 leads to the following summation by parts: (2.30) The third term of (2.21) can be analyzed in a similar way.By a reformulated form of the second equation (2.3) we have in which the last step comes from ∇ h • u m+1/2 = 0 and u m+1/2 • n = 0 on ∂Ω.
As a result, a substitution of (2.29), (2.30) and (2.33) into (2.21)becomes By the definition of the alternate numerical energy, we arrive at This in turn shows that the modified energy is non-increasing in time.Summing over time for (2.35) yields Then we obtain the unconditional energy stability for the second order scheme (2.5)-(2.7).
Remark 2.2.It is observed that data for two initial time steps, either (φ 0 , φ 1 ) or (φ −1 , φ 0 ), are needed for (2.5) -(2.7), since ours is a two-step scheme.In this article, we take for simplicity of presentation.Other initialization choices, such as φ −1 = φ 0 , or computing φ 1 by a first order temporal scheme, could be taken.Moreover, the energy stability and the second order temporal convergence rate are also expected to be available for these initial data choices; see the related works for the Cahn-Hilliard model [7,12].
The ∞ (0, T ; H 1 h ) bound of the numerical solution could be derived based on the weak energy stability (2.36).The following quadratic inequality is observed: with the discrete H 1 h norm introduced in (A.26).Then we arrive at the following bound, for any φ ∈ C Ω : Consequently, its combination with (2.36) yields the following estimate: so that a uniform in time bound for φ in ∞ (0, T ; H 1 h ) is available: Theorem 3.1.Let φ m ∈ C Ω be the solution to the scheme (2.5) -(2.7), with sufficient regularity assumption for Φ 0 and Φ 1 , then for any 1 ≤ ≤ M − 1, we have where C 10 and C 11 , given by (3.14) and (3.15), respectively, only depend on L x , L y , L z , ε and several Sobolev embedding constants, and are independent of h, s and final time T .
Proof.We observe that in which a triangle inequality was applied in the last step.Furthermore, motivated by the quadratic For the last term in (3.7), by the definition of χ φ m+1 , φ m in (2.4), a detailed expansion and a careful application of discrete Hölder inequality shows that in which the uniform in time estimate (3.4) was used again in the last step.Moreover, the • ∞ bound of φ k can be obtained with an application of a discrete Gagliardo-Nirenberg type inequality: see [3] for a detailed proof.
This in turn gives the 2 (0, T ; H 3 h ) bound of the numerical solution.
Note that C 10 and C 11 , given by (3.14) and (3.15), respectively, only depend on L x , L y , L z , ε and several Sobolev embedding constants, and independent of final time T , h and s.

Optimal rate convergence analysis
The convergence analysis is carried out in three steps.Firstly, in section 4.1, we obtain error functions by using a standard consistence analysis.In the following, we provide an estimate for the nonlinear error term in section 4.2.Finally, we recover an a-priori error assumption and present the optimal rate error estimate in sections 4.3 and 4.4 , respectively.

Error equations and consistency analysis
We assume that the exact solution has regularity of class R: To facilitate our error analysis, we need to construct an approximate solution to the chemical potential via the exact solution φ e .In addition, we note that the exact velocity u e is not divergencefree at the discrete level (∇ h • u e = 0).To overcome this difficulty, we must also construct an approximate solution to the velocity vector (again through the exact solution), which satisfies the divergence-free conditions at the discrete level.Therefore, we define the cell-centered grid functions for 1 ≤ m ≤ M , where P h is the discrete Helmholtz projection defined in equations (3.2), (3.3) in [3].We need to enforce the discrete homogeneous Neumann boundary conditions for the chemical potential: n • ∇ h Γ m+1/2 = 0, for all 1 ≤ m ≤ M , so that, in particular, is well defined.
With the assumed regularities, the constructed approximations Γ m+1/2 and U m+1/2 obey the following estimates: for 0 ≤ m ≤ M , where the constant C 12 > 0 is independent of h > 0 and s > 0.
It follows that (Φ, Γ, U) satisfies the numerical scheme with an O(s 2 + h 2 ) truncation error: where the local truncation error satisfies with s • M = T , and C 13 independent of h and s.
The numerical error functions are denoted as where We also observe that φ0 = φ1 ≡ 0, due to our initial value choices φ 0 = Φ 0 , φ 1 = Φ 1 .This fact will facilitate the convergence analysis in later sections.
Lemma 4.1.Assume the exact solution is of regularity class R.Then, for any for some constant C 14 that is independent of s, h, and m.
Before we carry out the stability analysis for the numerical error functions, we assume that the exact solution Φ and the constructed solutions Γ, U have the following regularity: In addition, we also set the • ∞ and H 3 h norms (introduced by (A.23) and (A.27)) for the numerical solution φ m as Note that we have an ∞ (0, T ; H 1 h ) and 2 (0, T ; H 3 h ) bound for the numerical solution, as given by (3.3), (3.16), respectively.Meanwhile, its ∞ (0, T ; ∞ ) bound is not available at present.This where and the constants C 28 , C 29 , C 30 are given by (4.60)-(4.62),respectively.
Proof.Taking inner product of (4.11) with −2∆ h φm+1/2 gives The term associated with the local truncation error term I 4,1 in (4.22) can be bounded in a straightforward way: The regular diffusion term I 4,2 in (4.22) has the following decomposition: The concave term I 4,2,2 in (4.24) can be controlled by For the nonlinear error term I 4,2,1 in (4.22), we start from the following expansion An application of discrete Hölder's inequality to its gradient shows that with the ∞ (0, T ; H 1 h ) estimate (3.4) for the numerical solution, the regularity assumption (4.16) for the exact solution and the a-priori set up (4.17) used.Then we get with For the second part in (4.28) , we observe that the maximum norm of the numerical error can be analyzed by an application of Gagliardo-Nirenberg type inequality in 3-D, similar to (3.11): With an application of the Young inequality to the second part in (4.28), we arrive at for any α > 0. Furthermore, the last term appearing in (4.31) can also be handled by Young's inequality: . (4.32) We can always choose an α, such that 1 5 α ≤ 1 64 and 4 5 α ≤ 1 32 , so that the following bound is available:

.33)
A similar estimate for the third part in (4.28) can also be derived as Consequently, a combination of (4.28), (4.29), (4.33) and (4.34) yields Next we focus our attention on the terms associated with the convection term and the highest order nonlinear diffusion.The analysis of this part is highly non-trivial.
The I 4,3 in (4.22) can be bounded by For the term I 4,4 in (4.22), the expansion (4.13) for the velocity numerical error indicates that The first term I 4,2,1 in (4.38) can be estimated in a standard way: where ) in which we used the property ∀v ∈ L 2 , for the Helmholtz projection operator P h , in the second step [3].Note that (M m 0 ) 2 and (M m−1 0 ) 2 are involved in the growth coefficient.
The second term I 4,2,2 in (4.38) can be expanded as based on the identity (u, P h v) = (P h u, P h v) for any vector u, v ∈ L 2 .The above inequality is the key reason for an ∞ (0, T ; H 1 h ) error estimate instead of the standard ∞ (0, T ; 2 ) one.The analysis for second term I 4,4,2,2 in (4.40) is straightforward: For the first term I 4,4,2,1 of (4.40), we start from an application of Cauchy inequality and discrete Hölder's inequality: The rest estimates are very similar to those for the regular diffusion.The inequalities (4.27) shows that where The bound for the first term in (4.44) can be derived in the same manner as in (4.29) with The bound for the second and third terms appearing in (4.44) follows form the proof of (4.30)-(4.35).Hence, the following two estimates are available, and the details are skipped for simplicity of presentation: Going back to (4.43)-(4.44),we arrive at with Consequently, from (4.22), (4.23), (4.36), (4.37) and (4.50), we obtain On the other hand, a similar estimate as (3.8) could be carried out: Then we get For the sake of convenience, we now make the coefficient on the right side of (4.53) explicit to each time step.Define Then I 2 can be bounded as Recall the definition of C 25 in (4.45), the value of which can be controlled as As a result, the stability inequality (4.53) can be rewritten as (4.18), with the following constants: C 29 = 2C α,ε,γ 4 • max 8

The result of an a-priori error assumption
As discussed in [3], in which a first order numerical scheme was analyzed, the discrete Gronwall inequality could not be directly applied to derive an error estimate from the stability inequality as in the form of (4.18), since D m+1 We then use this a priori assumption to prove that s(C 29 D m+1 ) < 1, provided s is small enough.Then we conclude the induction argument by proving that the error estimate holds at the updated time step m + 1.
First, we need the following technical result, which is a direct result of Young's inequality.The proof is skipped for brevity.
where C 32 , C 33 > 0 may depend upon the final time T but are independent of s and h.Then Proof.As an application of (4.63), the non-leading terms appearing on the right hand side of (4.19) for the expansion of D m+1 1 can be bounded as follows: Then we get, for any 0 ≤ m ≤ M − 1, using the L 8 s (0, T ) bound for M m 0 := φ m ∞ in (4.64),where C 34 > 0 is a constant independent of h and s.Using the same skill, the non-leading order of D m+1 2 can be bounded as follows where C 18 > 0 is a constant that is independent of h and s.Now, the leading terms appearing on the right hand side of (4.19)-(4.20)cannot be bounded in this way.We divide them into two groups G 1 and G 2 as follows: We must, therefore, rely upon (4.65).This bound implies Consequently, we see that The first term on the right hand side can be handled in the same way as (4.67): where C 39 > 0 is independent of s and h.The details of the proof are skipped for the sake of brevity.
Therefore, a combination of (

The main result: an error estimate
The following theorem is the main theoretical result of this article.The basic idea is to extend the a-priori error estimate (4.65) by an induction argument.Then, provided s and h are sufficiently small, for all positive integers , such that s • ≤ T , we have where C > 0 is independent of s and h.
Proof.Suppose that m + 1 ≤ M .By summing (4.18) we obtain We proceed by induction.Namely, suppose that (4.65) holds.Then, if h and s are sufficiently small -as required in the proof of the last theorem -considering (4.66) and using φ−1 ≡ φ0 ≡ 0, we have Hence where C 40 > 0 is a constant that is independent of s and h.Using the discrete Gronwall inequality gives where C 41 > 0 is a constant that is independent of s and h.Consequently, the a priori assumption

Numerical Experiments
In this section, we perform some numerical tests in two-dimensional space to verify the accuracy and efficiency of the proposed numerical scheme (2.5)-(2.7).The coupled systems are solved by the Full Approximation Scheme (FAS) under the nonlinear multigrid framework in [8,23].Here we omit the details for brevity; more details in [4,23]

Convergence rate, energy dissipation and mass conservation test
To estimate the convergence rate, we perform the Cauchy-type convergence as in [1,4,9,16,21,22] on a square Ω = [0, ( The homogeneous Neumann boundary conditions are imposed for φ, µ and p.In this test, the Cauchy difference is defined as δ φ = φ h f − I f c φ hc , where h c = 2h f and I f c is a bilinear interpolation operator that maps the coarse grid approximation u hc onto the fine grid (we applied nearest matlab interpolation function).We take a liner refinement path, i.e. s = Ch.At the final time T = 0.8, we expect the global error to be O(s 2 ) + O(h 2 ) = O(h 2 ) under the 2 norm, as h, s → 0. The other parameters are given by L x = L y = 3.2, s = 0.05h, ε = 0.2 and γ = 2.The norms of Cauchy difference, the convergence rates, the average number of V-cycle and average CPU time for one time step can be found in Table 1, which confirms our second order convergence rate expectation and indicates the efficiency of the proposed numerical scheme.The evolutions of discrete energy and mass for the simulation, associated with Table 1 for the h = 3.2  512 , are presented in Figure 1.The energy dissipation property is clearly demonstrated in the evolutions of discrete energy in the figure.And also, the evolution of discrete mass indicates the mass conservative property, with Ω φ(x, y, 0)dx = −5.12.

Spinodal decomposition
In this test, we simulate the spinodal decomposition of a binary fluid in a Hele-Shaw cell and show the effect of γ on the phase decomposition.The simulation parameters are similar to those in [23], with the parameters given by L x = L y = 6.4,ε = 0.03, h = 6.4/512, and s = 0.01.The initial data for this simulation is taken as a random field values φ 0 i,j = φ + 0.05 • (2r i,j − 1) with an average  From Figure .2, we observe that the particles indeed have a smaller shape factor for γ = 4 than for γ = 0 at same time, which coincides with the real physical states.Since larger γ would improve the fluid flow and enhance the energy dissipation.The energy evolution plot in Figure 3 implies that the energy decay are almost the same in the early stages of decomposition.Meanwhile, it is not precisely clear from the energy inequality that the larger γ will result in a larger energy dissipation rate [14,17,25].

Conclusions
A second order accurate energy stable numerical scheme for the Cahn-Hilliard-Hele-Shaw equations is proposed and analyzed in this article.The unique solvability and unconditional energy stability are proved, based on a rewritten form of the scheme, following a convexity analysis.At each time step of this scheme, an efficient nonlinear multigrid solver could be applied to the nonlinear equations associated with the finite difference approximation.At the theoretical side, an 2 (0, T ; H 3 h ) stability of the numerical scheme is established, in addition to the leading order energy stability.As an outcome of this estimate, we perform an ∞ (0, T ; H 1 h ) error estimate for the numerical scheme, and an optimal rate convergence analysis is obtained.A few numerical simulation results are presented to demonstrate the accuracy and robustness of the proposed numerical scheme.

) with φ m+1/ 2 *
a known function, and homogeneous Neumann boundary conditions for both µ and p µ .
) in which the uniform in time estimate(3.4)  was utilized in the third step and another quadratic inequality (a + b) 2 ≤ 2(a 2 + b 2 ) was used in the last step.An application of the Cauchy inequality gives an estimate of the leading term in (3.7): bound will be justified by later analysis.The following theorem states the stability of the numerical error functions satisfying the error equations by (4.11) -(4.11).Theorem 4.2.Assume the exact solution is of regularity class R. Then the error function φm obeys the following discrete energy stability law: for any 1

1 and D m+1 2 do
not have a uniform bound.Instead, we have to use an induction argument to establish the convergence analysis.Specifically, we assume, as an induction hypothesis, that the desired error estimate holds at an arbitrary time step m (0 ≤ m ≤ M − 1).

Theorem 4 . 5 .
Suppose that h and s are sufficiently small and the following error estimate is valid up to the time step t m := m • s, for 2 ≤ m ≤ M − 1:

( 4 .
65) can be justified at time step t m+1 by taking C 32 = C 40 , C 33 = C 41 .This completes the induction argument, and the proof of Theorem 4.6 is finished.
are referred to the readers.In the following tests, all the numerical experiments were performed with Fortran90 on Thinkpad W541 running with Intel Core i7-4800MQ at 2.80Ghz with 7.4GB memory under the Ubuntu 14.04.The general parameters of FAS are finest grid 2 × 2, pre-and post-smooth steps ν 1 = ν 2 = 2 and stopping tolerance tol = 10 −10 .

Figure 1 :
Figure1: The evolutions of discrete energy and mass for the simulation depicted in Table1for the h = 3.2/512 case.

Table 1 :
Errors, convergence rates, average iteration numbers and average CPU time (in seconds) for each time step.Rate #V's T cpu (h f )