Energy Stability and Convergence of SAV Block-centered Finite Difference Method for Gradient Flows

TL;DR

This paper develops and analyzes a second-order accurate block-centered finite difference method for gradient flows, demonstrating its stability, convergence, and efficiency through theoretical proofs and numerical experiments.

Contribution

It introduces a novel SAV/CN-BCFD scheme with rigorous stability and convergence analysis for gradient flows, enhancing accuracy and efficiency.

Findings

Second-order accuracy in time and space confirmed

Scheme is stable and convergent under various norms

Numerical tests verify robustness and efficiency

Abstract

We present in this paper construction and analysis of a block-centered finite difference method for the spatial discretization of the scalar auxiliary variable Crank-Nicolson scheme (SAV/CN-BCFD) for gradient flows, and show rigorously that scheme is second-order in both time and space in various discrete norms. When equipped with an adaptive time strategy, the SAV/CN-BCFD scheme is accurate and extremely efficient. Numerical experiments on typical Allen-Cahn and Cahn-Hilliard equations are presented to verify our theoretical results and to show the robustness and accuracy of the SAV/CN-BCFD scheme.

Figures

Tables

Table 1. Table 1: Errors and convergence rates of Example 1.

$h$	${\Vert e_Z\Vert }_{\infty ,2}$	Rate	${\Vert e_{\text{d}Z}\Vert }_{\infty ,2}$	Rate	${\Vert e_W\Vert }_{\infty }$	Rate
$1/10$	6.36E-3	—	5.96E-2	—	5.93E-3	—
$1/20$	1.59E-3	2.00	1.57E-2	1.93	1.47E-3	2.01
$1/40$	3.98E-4	2.00	3.98E-3	1.98	3.69E-4	2.00
$1/80$	9.96E-5	2.00	9.98E-4	1.99	9.23E-5	2.00

Table 2. Table 2: Errors and convergence rates of example 2.

$h$	${\Vert e_Z\Vert }_{\infty ,2}$	Rate	${\Vert e_{\text{d}Z}\Vert }_{\infty ,2}$	Rate	${\Vert e_R\Vert }_{\infty }$	Rate
$1/10$	5.49E-3	—	2.78E-2	—	4.88E-3	—
$1/20$	1.36E-3	2.01	6.91E-3	2.01	1.20E-3	2.02
$1/40$	3.41E-4	2.00	1.73E-3	2.00	3.00E-4	2.00
$1/80$	8.51E-5	2.00	4.31E-4	2.00	7.49E-5	2.00

Table 3. Table 3: Errors and convergence rates of example 2.

$h$	${\Vert e_W\Vert }_{2,2}$	Rate	${\Vert e_{\text{d}W}\Vert }_{2,2}$	Rate
$1/10$	2.50E-2	—	2.18E-1	—
$1/20$	6.11E-3	2.03	5.46E-2	2.00
$1/40$	1.52E-3	2.01	1.37E-2	2.00
$1/80$	3.79E-4	2.00	3.42E-3	2.00

Equations

$$E(\varphi )={\int }_{\Omega }(\frac{\lambda }{2}{\varphi }^2+\frac{1}{2}|\nabla \varphi {|}^2)d\text{x}+E_1(\varphi ),$$

$$\left\{\begin{array}[]{l}\displaystyle\frac{\partial\phi}{\partial t}=M\mathcal{G}\mu,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}in~{}\Omega\times J,\\ \displaystyle\mu=-\Delta\phi+\lambda\phi+F^{\prime}(\phi),~{}~{}~{}~{}in~{}\Omega\times J,\end{array}\right.$$

$$\left\{\begin{array}[]{l}\partial_{\textbf{n}}\phi|_{\partial\Omega}=0,~{}~{}\partial_{\textbf{n}}\mu|_{\partial\Omega}=0,\\ \phi|_{t=0}=\phi_{0},\end{array}\right.$$

$$\frac{dE}{dt}={\int }_{\Omega }\frac{\partial \varphi }{\partial t}\mu d\text{x}=M{\int }_{\Omega }\mu \mathcal{G}\mu d\text{x}\le 0.$$

$$\frac{\partial \varphi }{\partial t}=M\mathcal{G}\mu ,$$

$$\mu =-\Delta \varphi +\lambda \varphi +\frac{r}{\sqrt{E_1(\varphi )}}{F}^{\prime }(\varphi ),$$

$$r_t=\frac{1}{2\sqrt{E_1(\varphi )}}{\int }_{\Omega }{F}^{\prime }(\varphi ){\varphi }_{t}d\text{x},$$

$$\frac{{\varphi }^{n+1}-{\varphi }^n}{\Delta t}=M\mathcal{G}{\mu }^{n+1/2},$$

$${\mu }^{n+1/2}=-\Delta {\varphi }^{n+1/2}+\lambda {\varphi }^{n+1/2}+\frac{r^{n+1/2}}{\sqrt{E_1({\tilde{\varphi }}^{n+1/2})}}{F}^{\prime }({\tilde{\varphi }}^{n+1/2}),$$

$$\frac{r^{n+1}-r^n}{\Delta t}=\frac{1}{2\sqrt{E_1({\tilde{\varphi }}^{n+1/2})}}{\int }_{\Omega }{F}^{\prime }({\tilde{\varphi }}^{n+1/2})\frac{{\varphi }^{n+1}-{\varphi }^n}{\Delta t}d\text{x},$$

$${\tilde{\varphi }}^{n+1/2}=\frac{1}{2}(3{\varphi }^n-{\varphi }^{n-1}).$$

$$\Vert v{\Vert }_{L^m(\Omega )}={\left({\int }_{\Omega }|v{|}^{m}d\Omega \right)}^{1/m}.$$

$$(f,g)=(f,g{)}_{L^2(\Omega )}={\int }_{\Omega }fgd\Omega$$

$$W^{k,p}(\Omega )=\{g:\Vert g{\Vert }_{W_p^k(\Omega )}<\infty \},$$

$$\Vert g{\Vert }_{W^{k,p}(\Omega )}={\left(\sum _{|\alpha |\le k}\Vert D^{\alpha }g{\Vert }_{L^p(\Omega )}^p\right)}^{1/p}.$$

$$(x_{i+1/2},y_{j+1/2}),i=0,...,N_x,j=0,...,N_y,$$

$$x_i=(x_{i-\frac{1}{2}}+x_{i+\frac{1}{2}})/2,i=1,...,N_x,$$

$$h_x=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}},i=1,...,N_x,$$

$$y_j=(y_{j-\frac{1}{2}}+y_{j+\frac{1}{2}})/2,j=1,...,N_y,$$

$$h_y=y_{j+\frac{1}{2}}-y_{j-\frac{1}{2}},j=1,...,N_y,$$

$$(f,g{)}_m=\sum _{i=1}^{N_x}\sum _{j=1}^{N_y}{h}_x{h}_y{f}_{i,j}{g}_{i,j},$$

$$(f,g{)}_x=\sum _{i=1}^{N_x-1}\sum _{j=1}^{N_y}{h}_x{h}_y{f}_{i+\frac{1}{2},j}{g}_{i+\frac{1}{2},j},$$

$$(f,g{)}_y=\sum _{i=1}^{N_x}\sum _{j=1}^{N_y-1}{h}_x{h}_y{f}_{i,j+\frac{1}{2}}{g}_{i,j+\frac{1}{2}},$$

$$(\text{v},\text{r}{)}_{TM}=(v_1,r_1{)}_x+(v_2,r_2{)}_y.$$

$$[d_{x}g{]}_{i+\frac{1}{2},j}=(g_{i+1,j}-g_{i,j})/h_x,$$

$$[d_{y}g{]}_{i,j+\frac{1}{2}}=(g_{i,j+1}-g_{i,j})/h_y,$$

$$[D_{x}g{]}_{i,j}=(g_{i+\frac{1}{2},j}-g_{i-\frac{1}{2},j})/h_x,$$

$$[D_{y}g{]}_{i,j}=(g_{i,j+\frac{1}{2}}-g_{i,j-\frac{1}{2}})/h_y,$$

$$[d_{t}g{]}_{i,j}^n=(g_{i,j}^n-g_{i,j}^{n-1})/\Delta t.$$

$$(q,D_x{w}_1{)}_m=-(d_{x}q,w_1{)}_x,$$

$$(q,D_y{w}_2{)}_m=-(d_{y}q,w_2{)}_y.$$

$$[d_{t}Z{]}_{i,j}^{n+1}=M[D_x{d}_{x}W+D_y{d}_{y}W{]}_{i,j}^{n+1/2},$$

$$W_{i,j}^{n+1/2}=-[D_x{d}_{x}Z+D_y{d}_{y}Z{]}_{i,j}^{n+1/2}+\lambda Z_{i,j}^{n+1/2}$$

$$+\frac{R^{n+1/2}}{\sqrt{E_1^h({\tilde{Z}}^{n+1/2})}}{F}^{\prime }({\tilde{Z}}_{i,j}^{n+1/2}),$$

$$d_t{R}^{n+1}=\frac{1}{2\sqrt{E_1^h({\tilde{Z}}^{n+1/2})}}(F^{\prime }({\tilde{Z}}^{n+1/2}),d_t{Z}^{n+1}{)}_m,$$

$$E_1^h({\tilde{Z}}^{n+1/2})=\sum _{i=1}^{N_x}\sum _{j=1}^{N_y}{h}_x{h}_{y}F({\tilde{Z}}_{i,j}^{n+1/2}).$$

$$\left\{\begin{array}[]{l}\displaystyle[d_{x}Z]_{1/2,j}^{n}=[d_{x}Z]_{N_{x}+1/2,j}^{n}=0,\quad 1\leq j\leq N_{y},\\ \displaystyle[d_{y}Z]_{i,1/2}^{n}=[d_{y}Z]_{i,N_{y}+1/2}^{n}=0,\quad 1\leq i\leq N_{x},\\ \displaystyle[d_{x}W]_{1/2,j}^{n}=[d_{x}W]_{N_{x}+1/2,j}^{n}=0,\quad 1\leq j\leq N_{y},\\ \displaystyle[d_{y}W]_{i,1/2}^{n}=[d_{y}W]_{i,N_{y}+1/2}^{n}=0,\quad 1\leq i\leq N_{x},\\ \displaystyle Z_{i,j}^{0}=\phi_{0,i,j},1\leq i\leq N_{x},1\leq j\leq N_{y}.\end{array}\right.$$

$$[d_{t}Z{]}_{i,j}^{n+1}=-M{W}_{i,j}^{n+1/2},$$

$$W_{i,j}^{n+1/2}=-[D_x{d}_{x}Z+D_y{d}_{y}Z{]}_{i,j}^{n+1/2}+\lambda Z_{i,j}^{n+1/2}$$

Full text

{AMS}

65M06, 65M12, 65M15, 35K20, 35K35, 65Z05

1 Introduction

Gradient flows are widely used in mathematical models for problems in many fields of science and engineering, particularly in materials science and fluid dynamics, cf. [1, 2, 31, 21] and the references therein. Therefore it is important to develop efficient and accurate numerical schemes for their simulation. There exists an extensive literature on the numerical analysis of gradient flows, see for instance [3, 14, 7, 10, 23, 9, 15] and the references therein.

In the algorithm design of gradient flows, an important goal is to guarantee the energy stability at the discrete level, in order to capture the correct long-time dynamics of the system and provide enough flexibility for dealing with the stiffness problem induced by the thin interface. Many schemes for gradient flows are based on the traditional fully-implicit or explicit discretization for the nonlinear term, which may suffer from harsh time step constraint due to the thin interfacial width [11, 22]. In order to deal with this problem, the convex splitting approach [18, 24, 16] and linear stabilization approach [17, 22, 28, 33] have been widely used to construct unconditionally energy stable schemes. However, the convex splitting approach usually leads to nonlinear schemes and linear stabilization approach is usually limited to first-order accuracy.

Recently, a novel numerical method, the so called invariant energy quadratization (IEQ), was proposed in [29, 32, 30]. This method is a generalization of the method of Lagrange multipliers or of auxiliary variable. The IEQ approach is remarkable as it permits us to construct linear, unconditionally stable, and second-order unconditionally energy stable schemes for a large class of gradient flows. However, it leads to coupled systems with variable coefficients that may be difficult or expensive to solve. The scalar auxiliary variable (SAV) approach [21, 20] was inspired by the IEQ approach, which inherits its main advantages but overcomes many of its shortcomings. In particular, in a recent paper [19], the authors established the first-order convergence and error estimates for the semi-discrete SAV scheme.

In this paper, we construct a SAV/CN scheme with block-centered finite differences for gradient flows, carried out a rigorous stability and error analysis, and implemented an adaptive time stepping strategy so that the time step is only dictated by accuracy rather than by stability. The block-centered finite difference method can be thought as the lowest order Raviart-Thomas mixed element method with a suitable quadrature. Its main advantage over using a regular finite difference method is that it can approximate both the phase function and chemical potential with Neumann boundary conditions in the mixed formulation to second-order accuracy, and it guarantees local mass conservation. Our approach for error estimates here is very different from that in [19] which is based on deriving $H^{2}$ bounds for the numerical solution. However, this approach can not be used in the fully discrete case with finite-differences in space. The essential tools used in the proof are the summation-by-parts formulae both in space and time to derive energy stability, and an induction process to show that the discrete $L^{\infty}$ norm of the numerical solution is uniformly bounded, without assuming a uniform Lipschitz condition on the nonlinear potential. To the best of the authors’ knowledge, this is the first paper with rigorous proof of second-order convergence both in time and space for a linear scheme to a class of gradient flows without assuming a uniform Lipschitz condition for the nonlinear potential.

The paper is organized as follows. In Section 2, we describe our numerical scheme, including the temporal discretization and spacial discretization. In Section 3, we demonstrate the energy stability for our SAV/CN-BCFD scheme. In Section 4, we carry out error estimates for the SAV/CN-BCFD schemes. In Section 5, we present some numerical experiments to verify the energy stability and accuracy of the proposed schemes.

Throughout the paper we use $C$, with or without subscript, to denote a positive constant, which could have different values at different places.

2 The SAV/CN-BCFD scheme

Given a typical energy functional [19]:

[TABLE]

where $\Omega$ is a rectangular domain in $\mathbb{R}^{2}$, $\lambda\geq 0$ and $E_{1}(\phi)=\int_{\Omega}F(\phi)d\textbf{x}\geq-c_{0}$ for some $c_{0}>0$, i.e., it is bounded from below. We consider the following gradient flow:

[TABLE]

$J=(0,T]$, and $T$ denotes the final time. $M$ is the mobility constant which is positive. The chemical potential $\mu=\frac{\delta E}{\delta\phi}$. $\mathcal{G}=-1$ for the $L^{2}$ gradient flow and $\mathcal{G}=\Delta$ for the $H^{-1}$ gradient flow. $F(\phi)$ is the nonlinear free energy density and we focus on as an example, when $E_{1}(\phi)=\int_{\Omega}\alpha(1-\phi^{2})^{2}d\textbf{x}$, the $L^{2}$ and $H^{-1}$ gradient flows are the well-known Allen-Cahn and Cahn-Hilliard equations, respectively.

The boundary and initial conditions are as follows.

[TABLE]

where n is the unit outward normal vector of the domain $\Omega$. The equation satisfies the following energy dissipation law:

[TABLE]

2.1 The semi discrete SAV/CN scheme

We recall the SAV/CN scheme introduced in [21] first.

Let $C_{0}>c_{0}$ so that $E_{1}(\phi)+C_{0}>0$. Without loss of generality, we substitute $E_{1}$ with $E_{1}+C_{0}$ without changing the gradient flow. Then $E_{1}$ has a positive lower bound $\hat{C}_{0}=C_{0}-c_{0}$, which we still denote as $C_{0}$ for simplicity.

In the SAV approach, a scalar variable $r(t)=\sqrt{E_{1}(\phi)}$ is introduced, and the system (2.2) can be transformed into:

[TABLE]

Then, the SAV/CN scheme is given as follows:

[TABLE]

where $\phi^{n+1/2}=\frac{1}{2}(\phi^{n}+\phi^{n+1})$, $r^{n+1/2}=\frac{1}{2}(r^{n}+r^{n+1})$, $\tilde{\phi}^{n+1/2}$ can be any explicit approximation of $\phi(t^{n+1/2})$ with an error of $O(\Delta t^{2})$. For instance, we may let $\tilde{\phi}^{n+1/2}$ be the extrapolation by

[TABLE]

2.2 Spacial discretization

we apply the BCFD method on the staggered grids for the spacial discretization.

First we give some preliminaries. Let $L^{m}(\Omega)$ be the standard Banach space with norm

[TABLE]

For simplicity, let

[TABLE]

denote the $L^{2}(\Omega)$ inner product, $\|v\|_{\infty}=\|v\|_{L^{\infty}(\Omega)}.$ And $W^{k,p}(\Omega)$ be the standard Sobolev space

[TABLE]

where

[TABLE]

The grid points are denoted by

[TABLE]

and the notations similar to those in [26] are used.

[TABLE]

where $h_{x}$ and $h_{y}$ are grid spacings in $x$ and $y$ directions, and $N_{x}$ and $N_{y}$ are the number of grids along the $x$ and $y$ coordinates, respectively.

Let $g_{i,j},~{}g_{i+\frac{1}{2},j},~{}g_{i,j+\frac{1}{2}}$ denote $g(x_{i},y_{j}),~{}g(x_{i+\frac{1}{2}},y_{j}),~{}g(x_{i},y_{j+\frac{1}{2}}).$ Define the discrete inner products and norms as follows,

[TABLE]

For simplicity, from now on we always omit the superscript $n$ (the time level) if the omission does not cause conflicts. Define

[TABLE]

The following discrete-integration-by-part lemma [26] plays an important role in the analysis.

Lemma 2.1**.**

Let $q_{i,j},w_{1,i+1/2,j}~{}and~{}w_{2,i,j+1/2}$ be any values such that $w_{1,1/2,j}=w_{1,N_{x}+1/2,j}=w_{2,i,1/2}=w_{2,i,N_{y}+1/2}=0$, then

[TABLE]

[TABLE]

2.2.1 SAV/CV-BCFD scheme for $H^{-1}$ gradient flow

Let us denote by $\{Z^{n},W^{n},R^{n}\}_{n=0}^{N}$ the BCFD approximations to $\{\phi^{n},\mu^{n},r^{n}\}_{n=0}^{N}$. The scheme for $H^{-1}$ gradient flow is as follows: for $1\leq i\leq N_{x},~{}1\leq j\leq N_{y}$,

[TABLE]

where $\tilde{Z}^{n+1/2}$ is an approximation of $\tilde{\phi}^{n+1/2}$, and

[TABLE]

The boundary and initial approximations as follows.

[TABLE]

Remark. The solution procedure of the above scheme is described in detail in [21, 20], and hence is omitted here.

2.2.2 SAV/CV-BCFD scheme for $L^{2}$ gradient flow

Let us denote by $\{Z^{n},W^{n},R^{n}\}_{n=0}^{N}$ the BCFD approximations to $\{\phi^{n},\mu^{n},r^{n}\}_{n=0}^{N}$. The scheme for $L^{2}$ gradient flow is as follows: for $1\leq i\leq N_{x},~{}1\leq j\leq N_{y}$,

[TABLE]

where $\tilde{Z}^{n+1/2}$ is an approximation of $\tilde{\phi}^{n+1/2}$. The boundary and initial conditions are given in (2.16).

3 Unconditional energy stability

We demonstrate below that the full discrete SAV/CN-BCFD schemes are unconditionally energy stable with the discrete energy functional

[TABLE]

where $\textbf{d}Z=(d_{x}Z,d_{y}Z)$.

3.1 $H^{-1}$ gradient flow

Theorem 3.1**.**

The scheme (2.13)-(2.15) is unconditionally stable and the following discrete energy law holds for any $\Delta t$:

[TABLE]

Proof 3.2**.**

Multiplying equation (2.13) by $W_{i,j}^{n+1/2}h_{x}h_{y}$, and making summation on $i,j$ for $1\leq i\leq N_{x},~{}1\leq j\leq N_{y}$, we have

[TABLE]

Using Lemma 2.1, equation (3.3) can be transformed into the following:

[TABLE]

Multiplying equation (2.14) by $d_{t}Z^{n+1}_{i,j}h_{x}h_{y}$, and making summation on $i,j$ for $1\leq i\leq N_{x},~{}1\leq j\leq N_{y}$, we have

[TABLE]

Using Lemma 2.1 again, the first term on the right hand side of equation (3.5) can be written as:

[TABLE]

Multiplying equation (2.15) by $R^{n+1}+R^{n}$ leads to

[TABLE]

Combining equation (3.7) with equations (3.4) - (3.6) gives that

[TABLE]

which implies the desired results (3.2).

3.2 $L^{2}$ gradient flow

For $L^{2}$ gradient flow, we shall only state the result, as its proof is essentially the same as for the $H^{-1}$ gradient flow.

Theorem 3.3**.**

The scheme (2.17)-(2.19) is unconditionally stable and the following discrete energy law holds for any $\Delta t$:

[TABLE]

4 Error estimates

In this section, we derive our main results of this paper, i.e., error estimates for the fully discrete SAV/CN-BCFD schemes.

For simplicity, we set

[TABLE]

4.1 $H^{-1}$ gradient flow

We shall first derive error estimates for the case of $H^{-1}$ gradient flow.

Theorem 4.1**.**

We assume that $F(\phi)\in C^{3}(\mathbb{R})$ and $\phi\in W^{1,\infty}(J;W^{4,\infty}(\Omega))\cap W^{3,\infty}(J;W^{1,\infty}(\Omega)),~{}\mu\in L^{\infty}(J;W^{4,\infty}(\Omega))$. Let $\Delta t\leq C(h_{x}+h_{y})$, then for the discrete scheme (2.13)-(2.15), there exists a positive constant $C$ independent of $h_{x}$, $h_{y}$ and $\Delta t$ such that

[TABLE]

We shall split the proof of the above results into three lemmas below.

Lemma 4.2**.**

Under the condition of Theorem 4.1, there exists a positive constant $C$ independent of $h_{x}$, $h_{y}$ and $\Delta t$ such that

[TABLE]

Proof 4.3**.**

Denote

[TABLE]

Subtracting equation (2.5) from equation (2.13), we obtain

[TABLE]

where

[TABLE]

[TABLE]

Subtracting equation (2.6) from equation (2.14) leads to

[TABLE]

where

[TABLE]

Subtracting equation (2.7) from equation (2.15) gives that

[TABLE]

where

[TABLE]

Multiplying equation (4.3) by $e_{\mu,i,j}^{n+1/2}h_{x}h_{y}$, and making summation on $i,j$ for $1\leq i\leq N_{x},~{}1\leq j\leq N_{y}$, we have

[TABLE]

Using Lemma 2.1, we can write the first term on the right hand side of equation (4.10) as:

[TABLE]

Thanks to Cauchy-Schwarz inequality, the last two terms on the right hand side of equation (4.11) can be transformed into:

[TABLE]

Multiplying equation (4.6) by $d_{t}e_{\phi,i,j}^{n+1}h_{x}h_{y}$, and making summation on $i,j$ for $1\leq i\leq N_{x},~{}1\leq j\leq N_{y}$, we have

[TABLE]

Similar to the estimate of equation (3.6), the first term on the right hand side of equation (4.13) can be transformed into the following:

[TABLE]

The second term on the right hand side of equation (4.13) can be rewritten as follows:

[TABLE]

Recalling equation (4.3), the first term on the right hand side of equation (4.15) can be transformed into the following:

[TABLE]

Next, we shall first make the hypothesis that there exists a positive constant $C_{*}$ such that

[TABLE]

This hypothesis will be verified in Lemma 4.6 using a bootstrap argument.

Since $F(\phi)\in C^{3}(\mathbb{R})$, we have

[TABLE]

Using above and the Cauchy-Schwartz inequality, we can deduce that

[TABLE]

Similarly we can obtain

[TABLE]

Then equation (4.16) can be estimated by:

[TABLE]

Similar to (4.16), the second term on the right hand side of equation (4.15) can be controlled by:

[TABLE]

The third term on the right hand side of equation (4.13) can be estimated by:

[TABLE]

Multiplying equation (4.8) by $e_{r}^{n+1}+e_{r}^{n}$ leads to

[TABLE]

The first and second terms on the right hand side of equation (4.24) can be transformed into:

[TABLE]

Since $F(\phi)\in C^{3}(\mathbb{R})$, we have that

[TABLE]

*Recalling the midpoint approximation property of the rectangle quadrature formula, we can obtain that *

[TABLE]

*Combining equation (4.24) with equations (4.10)-(4.27) and using Cauchy-Schwarz inequality result in *

[TABLE]

From the discrete-integration-by-parts,

[TABLE]

we find

[TABLE]

Similarly we have

[TABLE]

Multiplying equation (4.29) by $\Delta t$, summing over $n,~{}n=0,1,\ldots,k$ and combining with equations (4.31) and (4.32), we can obtain (4.2).

Lemma 4.4**.**

Under the condition of Theorem 4.1, there exists a positive constant $C$ independent of $h_{x}$, $h_{y}$ and $\Delta t$ such that

[TABLE]

Proof 4.5**.**

Multiplying equation (4.3) by $e_{\phi,i,j}^{n+1/2}h_{x}h_{y}$, and making summation on $i,j$ for $1\leq i\leq N_{x},~{}1\leq j\leq N_{y}$, we have

[TABLE]

Using Lemma 2.1, the first term on the right hand side of equation (4.34) can be transformed into the following:

[TABLE]

The first term on the right hand side of equation (4.35) can be estimated as:

[TABLE]

In the $y$ direction, we have the similar estimates. Then the left hand side in (4.35) can be bounded by:

[TABLE]

Thanks to (4.6) and (4.15), the first term on the right hand side of (4.37) can be estimated as follows:

[TABLE]

Combining equation (4.34) with equations (4.37) and (4.38) and multiplying equation (4.29) by $2\Delta t$, summing over $n,~{}n=0,1,\ldots,k$ lead to (4.33).

Lemma 4.6**.**

Under the condition of Theorem 4.1, there exists a positive constant $C_{*}$ independent of $h_{x}$, $h_{y}$ and $\Delta t$ such that

[TABLE]

Proof 4.7**.**

We proceed in two steps.

Step 1 (Definition of $C_{*}$): Using the scheme (2.13)-(2.15) for $n=0$ and applying the inverse assumption, we can get the approximation $Z^{1}$ with the following property:

[TABLE]

where $h=\max\{h_{x},h_{y}\}$ and $\Pi_{h}$ is an bilinear interpolant operator with the following estimate [6]:

[TABLE]

Thus we can choose the positive constant $C_{*}$ independent of $h$ and $\Delta t$ such that

[TABLE]

Step 2 (Induction): By the definition of $C_{*}$, it is trivial that hypothesis (4.17) holds true for $l=1$. Supposing that $\|Z^{l-1}\|_{\infty}\leq C_{*}$ holds true for an integer $l=1,\cdots,k+1$, with the aid of the estimate (4.42), we have that

[TABLE]

Next we prove that $\|Z^{l}\|_{\infty}\leq C_{*}$ holds true. Since

[TABLE]

Let $\Delta t\leq C_{2}h$ and a positive constant $h_{1}$ be small enough to satisfy

[TABLE]

Then for $h\in(0,h_{1}],$ we derive from (4.40) that

[TABLE]

This completes the induction.

We are now in position to prove our main results.

Proof 4.8** (Proof of Theorem 4).**

Thanks to the above three lemmas, we can obtain

[TABLE]

Finally applying the discrete Gronwall’s inequality, we arrive at the desired result:

[TABLE]

Thus, the proof of Theorem 4 is complete.

4.2 $L^{2}$ gradient flow

For the $L^{2}$ gradient flow, we shall only state the error estimates below, as their proofs are essentially the same as for the $H^{-1}$ gradient flow.

Theorem 4.9**.**

We assume that $F(\phi)\in C^{3}(\mathbb{R})$ and $\phi\in W^{1,\infty}(J;W^{4,\infty}(\Omega))\cap W^{3,\infty}(J;W^{1,\infty}(\Omega))$ and $\Delta t\leq C(h_{x}+h_{y})$. Then for the discrete scheme (2.17)-(2.19), there exists a positive constant $C$ independent of $h_{x}$, $h_{y}$ and $\Delta t$ such that

[TABLE]

5 Numerical simulations

We present in this section various numerical experiments to verify the energy stability and accuracy of the proposed numerical schemes.

5.1 Accuracy test for Allen-Cahn and Cahn-Hilliard equations

We consider the free energy

[TABLE]

and for better accuracy, rewrite it as

[TABLE]

where $\beta$ is a positive number to be chosen. To apply our schemes (2.13)-(2.15) or (2.17)-(2.19) to the system (2.2), we drop the constant in the free energy and specify the operator $\mathcal{G}$, the energy $E_{1}(\phi)$ and $\lambda$ as follows:

[TABLE]

The system (2.2) becomes the standard Allen-Cahn equation with $s=0$, and the standard Cahn-Hilliard equation with $s=1$.

We denote

[TABLE]

where $q=\frac{1}{2},~{}1$ and $X=m,~{}TM$.

In the following simulations, we choose $\Omega=(0,1)\times(0,1)$ and $C_{0}=0$.

5.1.1 Convergence rates of the SAV/CN-BCFD scheme for Allen-Cahn equation

Example 1. We take $T=0.5,~{}~{}\mathcal{G}=-1,~{}~{}\beta=0,~{}~{}M=0.01,~{}\epsilon=0.08,\ \Delta t=5E-4$, and the initial solution $\phi_{0}=\cos(\pi x)\cos(\pi y)$. To get around the fact that we do not have possession of exact solution, we measure Cauchy error, which is similar to [5, 27, 7]. Specifically, the error between two different grid spacings $h$ and $\frac{h}{2}$ is calculated by $\|e_{\zeta}\|=\|\zeta_{h}-\zeta_{h/2}\|$.

The numerical results are listed in Table 1. we observe the second-order convergence predicted by the error estimates in Theorem 4.9.

5.1.2 Convergence rates of SAV/CN-BCFD scheme for Cahn-Hilliard equation

Example 2. We take $T=0.5,~{}~{}\mathcal{G}=\Delta,~{}~{}\beta=0,~{}~{}M=0.01,~{}\epsilon=0.2,\ \Delta t=5E-4$, with the same initial solution as in Example 1. The numerical results are listed in Tables 2 and 3. Again, we observe the expected second-order convergence rate in various discrete norms.

5.2 Coarsening dynamics and adaptive time stepping

In this example, we simulate the coarsening dynamics of the Cahn-Hilliard equation.

Since the scheme (2.13)-(2.15) is unconditionally energy stable, we can choose time steps according to accuracy only with an adaptive time stepping. Actually in many situations, the energy and solution of gradient flows can vary drastically in certain time intervals, but only slightly elsewhere. In order to maintain the desired accuracy, we adjust the time sizes based on an adaptive time-stepping strategy below (Ref. [13, 20]).

We update the time step size by using the formula

[TABLE]

where $\rho$ is a default safety coefficient, $tol$ is a reference tolerance, and $e$ is the relative error at each time level. In this simulation, we take

[TABLE]

with a random initial condition with values in $[-0.05,0.05]$, and the initial time step is taken as $\Delta t_{min}$.

To demonstrate the effectivity of the SAC/CN-BCFD scheme with adaptive time stepping, we compute the reference solutions with a small uniform time step $\Delta t=10^{-5}$ and a large uniform time step $\Delta t=10^{-3}$ respectively. Characteristic evolutions of the phase functions are presented in Fig. 1(c). We also present in Fig. 2 the energy evolutions and the roughness of interface, where the roughness measure function $R(t)$ is defined as follows:

[TABLE]

with $\bar{\phi}=\frac{1}{|\Omega|}\int_{\Omega}\phi d\Omega$. One observes that the solution obtained with adaptive time steps is consistent with the reference solution obtained with a small time step, while the solution with large time step deviates from the reference solution. This is also verified by both the energy evolutions and roughness measure function $R(t)$. We present in Fig. 3 the adaptive time steps for different $\epsilon=0.02,\ 0.01,\ 0.005$. We observe that there are about two-orders of magnitude variation in the time steps with the adaptive time stepping, which indicates that the adaptive time stepping for the SAV/CN-BCFD scheme is very efficient.

Acknowledgments

X. Li thanks for the financial support from China Scholarship Council. The work of J. Shen is supported in part by NSF grants DMS-1620262, DMS-1720442 and AFOSR grant FA9550-16-1-0102. The work of H. Rui is supported by the National Natural Science Foundation of China grant 11671233.