The fast scalar auxiliary variable approach with unconditional energy stability for nonlocal Cahn-Hilliard equation

TL;DR

This paper introduces a fast, energy-stable scalar auxiliary variable method for the nonlocal Cahn-Hilliard equation, improving computational efficiency and stability for modeling phase transitions with nonlocal effects.

Contribution

It develops an unconditionally energy-stable scheme for the nonlocal Cahn-Hilliard equation and proposes a fast computational method leveraging BTTB matrix structure.

Findings

Proved unconditional energy stability of the scheme.

Reduced computational work using BTTB matrix structure.

Numerical simulations confirm accuracy and efficiency.

Abstract

Comparing with the classical local gradient flow and phase field models, the nonlocal models such as nonlocal Cahn-Hilliard equations equipped with nonlocal diffusion operator can describe more practical phenomena for modeling phase transitions. In this paper, we construct an accurate and efficient scalar auxiliary variable approach for the nonlocal Cahn-Hilliard equation with general nonlinear potential. The first contribution is that we have proved the unconditional energy stability for nonlocal Cahn-Hilliard model and its semi-discrete schemes carefully and rigorously. Secondly, what we need to focus on is that the non-locality of the nonlocal diffusion term will lead the stiffness matrix to be almost full matrix which generates huge computational work and memory requirement. For spatial discretizaion by finite difference method, we find that the discretizaition for nonlocal operator will lead to a block-Toeplitz-Toeplitz-block (BTTB) matrix by applying four transformation operators. Based on this special structure, we present a fast procedure to reduce the computational work and memory requirement. Finally, several numerical simulations are demonstrated to verify the accuracy and efficiency of our proposed schemes.

Tables

Table 1. Table 1: SCHEM1: the L 2 subscript 𝐿 2 L_{2} errors, temporal convergence rates and CPU time for direct solver and fast solver with h r ​ e ​ f = 0.01 subscript ℎ 𝑟 𝑒 𝑓 0.01 h_{ref}=0.01 , Δ ​ t r ​ e ​ f = 2 − 14 ​ T Δ subscript 𝑡 𝑟 𝑒 𝑓 superscript 2 14 𝑇 \Delta t_{ref}=2^{-14}T for Example 1 with initial value ϕ 0 ​ ( x , y ) = 0.5 ​ sin ⁡ ( π ​ x ) ​ sin ⁡ ( π ​ y ) + 0.1 subscript italic-ϕ 0 𝑥 𝑦 0.5 𝜋 𝑥 𝜋 𝑦 0.1 \phi_{0}(x,y)=0.5\sin(\pi x)\sin(\pi y)+0.1 .

$\Delta t$	Direct Solver			Fast Solver
	$L_2$ error	Rate	CPU Time(s)	$L_2$ error	Rate	CPU Time(s)
$2^{-4}T$	2.5139e-3	-	538	2.5139e-3	-	10.88
$2^{-5}T$	1.5066e-3	0.7386	1035	1.5066e-3	0.7386	16.95
$2^{-6}T$	8.2386e-4	0.8708	2046	8.2386e-4	0.8708	27.04
$2^{-7}T$	4.2956e-4	0.9395	4045	4.2956e-4	0.9395	47.48
$2^{-8}T$	2.1803e-4	0.9783	7920	2.1803e-4	0.9783	82.00
$2^{-9}T$	1.0852e-4	1.0066	15758	1.0852e-4	1.0066	162

Table 2. Table 2: SCHEME2: the L 2 subscript 𝐿 2 L_{2} errors, temporal convergence rates and CPU time for direct solver and fast solver with h r ​ e ​ f = 0.01 subscript ℎ 𝑟 𝑒 𝑓 0.01 h_{ref}=0.01 , Δ ​ t r ​ e ​ f = 2 − 14 ​ T Δ subscript 𝑡 𝑟 𝑒 𝑓 superscript 2 14 𝑇 \Delta t_{ref}=2^{-14}T for Example 1 with initial value ϕ 0 ​ ( x , y ) = 0.5 ​ sin ⁡ ( π ​ x ) ​ sin ⁡ ( π ​ y ) + 0.1 subscript italic-ϕ 0 𝑥 𝑦 0.5 𝜋 𝑥 𝜋 𝑦 0.1 \phi_{0}(x,y)=0.5\sin(\pi x)\sin(\pi y)+0.1 .

$\Delta t$	Direct Solver			Fast Solver
	$L_2$ error	Rate	CPU Time(s)	$L_2$ error	Rate	CPU Time(s)
$2^{-4}T$	9.2979e-4	-	592	9.2979e-4	-	9.66
$2^{-5}T$	2.5391e-4	1.8726	1091	2.5391e-4	1.8726	15.65
$2^{-6}T$	6.6781e-5	1.9268	2103	6.6781e-5	1.9268	27.02
$2^{-7}T$	1.7155e-5	1.9608	4147	1.7155e-5	1.9608	47.98
$2^{-8}T$	4.3486e-6	1.9800	7965	4.3486e-6	1.9800	82.82
$2^{-9}T$	1.0932e-6	1.9920	15834	1.0932e-6	1.9920	140.6

Table 3. Table 3: SCHEM1: the L 2 subscript 𝐿 2 L_{2} errors, spatial convergence rates and CPU time for direct solver and fast solver with h r ​ e ​ f = 2 − 10 subscript ℎ 𝑟 𝑒 𝑓 superscript 2 10 h_{ref}=2^{-10} , Δ ​ t r ​ e ​ f = 5 ​ e − 5 Δ subscript 𝑡 𝑟 𝑒 𝑓 5 superscript 𝑒 5 \Delta t_{ref}=5e^{-5} for Example 1 with initial value ϕ 0 ​ ( x , y ) = 0.5 ​ sin ⁡ ( π ​ x ) ​ sin ⁡ ( π ​ y ) + 0.1 subscript italic-ϕ 0 𝑥 𝑦 0.5 𝜋 𝑥 𝜋 𝑦 0.1 \phi_{0}(x,y)=0.5\sin(\pi x)\sin(\pi y)+0.1 .

$h$	Direct Solver			Fast Solver
	$L_2$ error	Rate	CPU Time(s)	$L_2$ error	Rate	CPU Time(s)
$2^{-3}$	6.7255e-3	-	6.92	6.7255e-3	-	12.13
$2^{-4}$	1.4786e-3	2.1854	87.86	1.4786e-3	2.1854	19.42
$2^{-5}$	3.4718e-4	2.0905	2689	3.4718e-4	2.0905	72.54
$2^{-6}$	8.3250e-5	2.0602	12436	8.3250e-5	2.0602	280.8
$2^{-7}$	N/A	N/A	Out of Memory	1.9526e-5	2.0921	1527
$2^{-8}$	N/A	N/A	Out of Memory	3.8761e-6	2.3327	7979

Table 4. Table 4: SCHEM2: The L 2 subscript 𝐿 2 L_{2} errors, spatial convergence rates and CPU time for direct solver and fast solver with h r ​ e ​ f = 2 − 10 subscript ℎ 𝑟 𝑒 𝑓 superscript 2 10 h_{ref}=2^{-10} , Δ ​ t r ​ e ​ f = 5 ​ e − 5 Δ subscript 𝑡 𝑟 𝑒 𝑓 5 superscript 𝑒 5 \Delta t_{ref}=5e^{-5} for Example 1 with initial value ϕ 0 ​ ( x , y ) = 0.5 ​ sin ⁡ ( π ​ x ) ​ sin ⁡ ( π ​ y ) + 0.1 subscript italic-ϕ 0 𝑥 𝑦 0.5 𝜋 𝑥 𝜋 𝑦 0.1 \phi_{0}(x,y)=0.5\sin(\pi x)\sin(\pi y)+0.1 .

$h$	Direct Solver			Fast Solver
	$L_2$ error	Rate	CPU Time(s)	$L_2$ error	Rate	CPU Time(s)
$2^{-3}$	6.7282e-3	-	6.76	6.7282e-3	-	11.88
$2^{-4}$	1.4793e-3	2.1853	92.36	1.4793e-3	2.1853	16.28
$2^{-5}$	3.4733e-4	2.0905	2866	3.4733e-4	2.0905	68.86
$2^{-6}$	8.3287e-5	2.0601	12738	8.3287e-5	2.0601	275.4
$2^{-7}$	N/A	N/A	Out of Memory	1.9535e-5	2.0920	1531
$2^{-8}$	N/A	N/A	Out of Memory	3.8778e-6	2.3328	7233

Equations

$$\frac{\partial \varphi }{\partial t}+M(-\Delta )(-{ϵ}^2\Delta \varphi +F^{\prime }(\varphi ))=0,(\text{x},t)\in \Omega \times Q,$$

$$F(y)=\frac{1}{4}(y^2-1{)}^2,y\in (-\infty ,\infty ),$$

$$E(\varphi )={\int }_{\Omega }(\frac{1}{2}{ϵ}^2|\nabla \varphi {|}^2+F(\varphi ))d\text{x}.$$

$$E(\varphi )={\int }_{\Omega }\left(F(\varphi )+\frac{{ϵ}^2}{4}{\int }_{\Omega }J(\text{x}-\text{y}){\left[\varphi (\text{x})-\varphi (\text{y})\right]}^{2}d\text{y}\right)d\text{x},$$

$$\frac{\partial \varphi }{\partial t}+M\mathcal{L}(-{ϵ}^2\mathcal{L}\varphi +F^{\prime }(\varphi ))=0,(\text{x},t)\in \Omega \times J,$$

$$\mathcal{L}u(\text{x})={\int }_{\Omega }J(\text{x}-\text{y})[u(\text{x})-u(\text{y})]d\text{y},\forall \text{x}\in \Omega \subseteq {\mathbb{R}}^n.$$

$$(u,v)={\int }_{\Omega }uvd\text{x},\Vert u{\Vert }_2=(u,u{)}^{\frac{1}{2}},\forall u,v\in L^2(\Omega ).$$

$$\left(\mathcal{L}u(\text{x}),u(\text{x})\right)=2{\int }_{\Omega }{\int }_{\Omega }J(\text{x}-\text{y})[u(\text{x})-u(\text{y}){]}^{2}d\text{y}d\text{x}\ge 0.$$

$$\frac{\partial \varphi }{\partial t}=M\mathcal{G}\mu ,(\text{x},t)\in \Omega \times Q,$$

$$\frac{dE[\varphi (t)]}{dt}=\frac{\delta E(\varphi )}{\delta \varphi }\cdot \frac{\partial \varphi }{\partial t}=-(\mu ,M\mathcal{L}\mu )\le 0.$$

$$\frac{dE[\varphi (t)]}{dt}$$

$$={\int }_{\Omega }{F}^{\prime }(\varphi ){\varphi }_{t}dx+\frac{{ϵ}^2}{2}{\int }_{\Omega }{\int }_{\Omega }J(\text{x}-\text{y})\left[\varphi (\text{x})-\varphi (\text{y})\right]{\varphi }_t(\text{x})d\text{y}d\text{x}$$

$$-\frac{{ϵ}^2}{2}{\int }_{\Omega }{\int }_{\Omega }J(\text{x}-\text{y})\left[\varphi (\text{x})-\varphi (\text{y})\right]{\varphi }_t(\text{y})d\text{y}d\text{x}.$$

$${\int }_{\Omega }{\int }_{\Omega }J(\text{x}-\text{y})\left[\varphi (\text{x})-\varphi (\text{y})\right]{\varphi }_t(\text{y})d\text{y}d\text{x}$$

$$=-{\int }_{\Omega }{\int }_{\Omega }J(\text{x}-\text{y})\left[\varphi (\text{x})-\varphi (\text{y})\right]{\varphi }_t(\text{x})d\text{y}d\text{x}.$$

$$\frac{dE[\varphi (t)]}{dt}$$

$$\left\{\begin{array}[]{rlr}\displaystyle\frac{\partial\phi}{\partial t}&=-M\mathcal{L}\mu,&(\textbf{x},t)\in\Omega\times J,\\ \mu&=\displaystyle F^{\prime}(\phi)+\epsilon^{2}\mathcal{L}\phi(\textbf{x}),&(\textbf{x},t)\in\Omega\times J.\end{array}\right.$$

$$r(t)=\sqrt{E_1+C_0}.$$

$$\left\{\begin{array}[]{rl}\displaystyle\frac{\partial\phi}{\partial t}&=-M\mathcal{L}\mu,\\ \displaystyle\mu&\displaystyle=\epsilon^{2}\mathcal{L}\phi(\textbf{x})+\frac{r}{\sqrt{E_{1}(\phi)+C_{0}}}F^{\prime}(\phi),\\ \displaystyle\frac{dr}{dt}&\displaystyle=\frac{1}{2\sqrt{E_{1}(\phi)+C_{0}}}\int_{\Omega}F^{\prime}(\phi)\frac{\partial\phi}{\partial t}d\textbf{x}.\end{array}\right.$$

$$\frac{\partial }{\partial t}\left[{ϵ}^2(\mathcal{L}\varphi ,\varphi )+|r{|}^2\right]=-M(\mathcal{L}\mu ,\mu ),$$

$$\Delta t=T/N,t^n=n\Delta t,\text{for}n\le N.$$

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

$${\mu }^{n+1}={ϵ}^2\mathcal{L}{\varphi }^{n+1}+\frac{r^{n+1}}{\sqrt{E_1({\tilde{\varphi }}^{n+1})+C_0}}{F}^{\prime }({\tilde{\varphi }}^{n+1}),$$

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

$${\tilde{\varphi }}^{n+1}=2{\varphi }^n-{\varphi }^{n-1}.$$

$$\frac{{\tilde{\varphi }}^{n+1}-{\varphi }^n}{\Delta t}=-M\mathcal{L}({ϵ}^2\mathcal{L}{\tilde{\varphi }}^{n+1}+F^{\prime }({\varphi }^n)).$$

$$\frac{{ϵ}^2}{2}(\mathcal{L}{\varphi }^{n+1},{\varphi }^{n+1})+|r^{n+1}{|}^2\le \frac{{ϵ}^2}{2}(\mathcal{L}{\varphi }^n,{\varphi }^n)+|r^n{|}^2.$$

$${ϵ}^2(\mathcal{L}{\varphi }^{n+1},{\varphi }^{n+1}-{\varphi }^n)+2(r^{n+1},r^{n+1}-r^n)=-M\Delta t(\mathcal{L}{\mu }^{n+1},{\mu }^{n+1})\le 0.$$

$$\frac{{ϵ}^2}{2}(\mathcal{L}{\varphi }^{n+1},{\varphi }^{n+1})-\frac{{ϵ}^2}{2}(\mathcal{L}{\varphi }^n,{\varphi }^n)+\frac{{ϵ}^2}{2}(\mathcal{L}({\varphi }^{n+1}-{\varphi }^n),{\varphi }^{n+1}-{\varphi }^n)$$

$$+|r^{n+1}{|}^2-|r^n{|}^2+|r^{n+1}-r^n{|}^2=-M\Delta t(\mathcal{L}{\mu }^{n+1},{\mu }^{n+1})\le 0,$$

$$\frac{3{\varphi }^{n+1}-4{\varphi }^n+{\varphi }^{n-1}}{2\Delta t}=-M\mathcal{L}{\mu }^{n+1},$$

$${\mu }^{n+1}={ϵ}^2\mathcal{L}{\varphi }^{n+1}+\frac{r^{n+1}}{\sqrt{E_1({\tilde{\varphi }}^{n+\frac{1}{2}})+C_0}}{F}^{\prime }({\tilde{\varphi }}^{n+\frac{1}{2}}),$$

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

$$\frac{{ϵ}^2}{2}(\mathcal{L}{\varphi }^{n+1},{\varphi }^{n+1})+\frac{{ϵ}^2}{2}(\mathcal{L}(2{\varphi }^{n+1}-{\varphi }^n),2{\varphi }^{n+1}-{\varphi }^n)+|r^{n+1}{|}^2+|2r^{n+1}-r^n{|}^2$$

Full text

1 Introduction

Phase field models have been extensively applied to study the dynamics of different material phases via order parameters, such as phase transformations in binary alloys [19, 20]; epitaxial thin film growth [28, 29]; crystal faceting [6]; multi-phase fluid flow [1, 10]. They can describe the evolution of complex, morphology-changing surfaces and construct a general framework to take the consideration of more physical effects.

There exist many advantages in the phase field models from the mathematical point of view. Specially, since the phase field models are usually energy stable (thermodynamics-consistent) and well-posed, which is based on the energy variational approach, it is possible to perform effective numerical analysis and carry out reliable and accurate computer simulations. The significant goal is to preserve the energy stable property at the discrete level irrespectively of the coarseness of the discretization in time and space. Schemes with this property is extremely preferred for solving diffusive systems due to the fact that it is not only critical for the numerical scheme to capture the correct long time dynamics of the system, but also supplies sufficient flexibility for dealing with the stiffness issue. Moreover, the noncompliance of energy dissipation laws may lead to spurious numerical approximations if the mesh or time step sizes are not controlled carefully. However, due to the thin interface, it is a quite difficult issue to construct unconditionally energy stable schemes for phase field models such as the Cahn-Hilliard and Allen-Cahn equations with general nonlinear potential. Many efforts had been done in order to solve this problem (c.f. [9, 18, 22, 24, 26, 27, 31]).

Recently, nonlocal phase filed models such as nonlocal Cahn-Hilliard equation have attracted more and more attentions and been used in many fields involving physics, materials science, finance and image processing [7, 8, 13]. Many important phenomena are well described by nonlocal model with nonlocal diffusion term. Bates and Han [3, 4] investigated the well-posedness of the nonlocal Cahn-Hilliard equations equipped with Dirichlet or Neumann boundary condition. Stabilized linear semi-implicit schemes for the nonlocal Cahn-Hilliard equation were considered by Du et al. and the energy stabilities were established for two methods in the fully discrete sense [13]. Guan et al. [17] presented second-order accurate, unconditionally uniquely solvable and unconditionally energy stable schemes for the nonlocal Cahn-Hilliard and nonlocal Allen-Cahn equations for a large class of interaction kernels. Some other relative models such as nonlocal Cahn-Hilliard-Navier-Stokes [14, 15, 16] have also been considered and analyzed by many researchers.

The classical Cahn-Hilliard equation is a nonlinear, fourth order in space, parabolic partial differential equation which is often used as a diffuse interface model for the phase separation of a binary alloy:

[TABLE]

where $Q=(0,T]$, $M$ is the mobility constant, the interface width satisfies $O(\epsilon)$, which is small compared to the characteristic length of the laboratory scale. The phase-field $\phi$ represents the difference of local relative concentrations for the two components of the mixture such that $\phi=\pm 1$ correspond to the pure phases of the material while $\phi\in(-1,1)$ corresponds to the transition between the two phases. $\textbf{x}\in\Omega\subseteq\mathbb{R}^{d}$. $F(\phi)$ is the nonlinear bulk potential and the most commonly used form Ginzburg-Landau double-well type potential is defined as follows [21, 23, 30]:

[TABLE]

The free energy takes the form:

[TABLE]

By replacing the Laplacian in the above local energy by the nonlocal diffusion operator $\mathcal{L}$ in [11], one can obtain the nonlocal free energy functional as follows [2, 5, 13]:

[TABLE]

where the kernel $J$ satisfies the following conditions:

(i) $J(\textbf{x})\geq 0$, $\forall\textbf{x}\in\Omega$;

(ii) $J(\textbf{x}-\textbf{y})=J(\textbf{y}-\textbf{x})$;

(iii) $J$ is $\Omega$-periodic.

In addition, we can also replace the first Laplician in classical Cahn-Hilliard model by another nonlocal operator $\mathcal{L}$ to obtain a general nonlocal Cahn-Hilliard model:

[TABLE]

where the nonlocal operator $\mathcal{L}$ on the function $u(\textbf{x})$ has been introduced by many articles such as [11, 12]

[TABLE]

The main goal of this paper is to construct accurate and efficient linear algorithms for the general nonlocal Cahn-Hilliard equation with general nonlinear potential and prove the unconditional energy stability for its semi-discrete schemes carefully and rigorously. In addition, considering the huge computational work and memory requirement in solving the linear system, we analyze the structure of the stiffness matrix and seek some effective fast solution method to reduce the computational work and memory requirement. By applying four transformation operators $\mathcal{A}_{1}$, $\mathcal{A}_{2}$, $\mathcal{A}_{3}$, $\mathcal{A}_{4}$, we transform the stiffness matrix into a block-Toeplitz-Toeplitz-block (BTTB) matrix. Then, a fast solution technique which is based on a fast Fourier transform is presented to solve a new linear system with BTTB stiffness matrix. The overall computational cost of the fast conjugate gradient method is $O(N$log${}^{2}N)$, since the number of iterations is $O($log$N)$ where $N$ is the number of unknowns. What we need to focus is that if one uses the Gaussian elimination method straightforwardly to this linear system, then it requires $O(N^{3})$ complexity. In addition, since $N\times N$ BTTB matrix is determined by only $2N-1$ entries rather than $N^{2}$ entries, the fast solver will make memory requirement $O(N)$ instead of $O(N^{2})$. Finally, various 2D numerical simulations are demonstrated to verify the accuracy and efficiency of our proposed schemes.

The paper is organized as follows. In Sect.2, we provide the nonlocal Cahn-Hilliard model with general nonlinear potential by energy variational approach and present some notations. In Sect.3, the linear, first and second order numerical scalar auxiliary variable approaches to construct unconditionally energy stable schemes for the nonlocal Cahn-Hilliard model are considered. In Sect.4, we analyse the structure of the stiffness matrix and seek some effective fast solution methods to reduce the computational work and memory requirement. In Sect.5, various 2D numerical simulations are demonstrated to verify the accuracy and efficiency of our proposed schemes.

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

2 The nonlocal Cahn-Hilliard model and relative notations

In this section, we introduce the nonlocal Cahn-Hilliard model with general nonlinear potential by energy variational approach and present some notations which will be used in the later analysis.

First, the inner product and norm of $L^{2}(\Omega)$ are defined by

[TABLE]

It is easy to obtain that

[TABLE]

Next, we give a brief introduction about how the nonlocal Cahn-Hilliard model is resulted from the energetic variation of the energy functional (1.1). Denoting its variational derivative as $\mu=\frac{\delta E}{\delta\phi}$, the general form of the gradient flow model can be written as [25]

[TABLE]

where $Q=(0,T]$, $M$ is the mobility constant, $\mu$ is the chemical potential, and $f(\phi)=F^{\prime}(\phi)$. The initial condition is $\phi|_{t=0}=\phi_{0}$. The equation (2.3) will be Cahn-Hilliard type system if $\mathcal{G}=\Delta$. The equation is supplemented with the following boundary condition: periodic, or $\frac{\partial\phi}{\partial\textbf{n}}=\frac{\partial\mu}{\partial\textbf{n}}=0$, where n is the unit outward normal vector on the boundary $\partial\Omega$.

Noting that for Cahn-Hilliard model, the operator $\mathcal{G}=-\mathcal{L}$ is non-positive, the free energy $E(\phi)$ in equation (1.1) is non-increasing,

[TABLE]

By some simple calculations, we obtain

[TABLE]

Using the condition $(ii)$ of kernel $J$ and variable substitution, we obtain

[TABLE]

Combining equations $(\ref{section1_eneegy_compute})$ and equation $(\ref{section1_eneegy_compute2})$, we obtain

[TABLE]

Then, combining the equations (1.2), (2.3), (2.4), (2.7) with $\mu=\frac{\delta E}{\delta\phi}$ and $\mathcal{G}=-\mathcal{L}$, we obtain the nonlocal Cahn-Hilliard model with general nonlinear potential:

[TABLE]

3 The scalar auxiliary variable semi-implicit schemes for nonlocal Cahn-Hilliard model

In this section, we execute and analyze linear, first and second order (in time) numerical scalar auxiliary variable (SAV) approaches to construct unconditionally energy stable schemes for the nonlocal Cahn-Hilliard model with general nonlinear potential. The SAV approach is proposed for a large class of gradient flows that describes energy dissipative physical systems. It leads to numerical schemes that enjoy linear second-order unconditionally energy stability, which is very significant to solve the stiffness issue from the thin interface. For the SAV scheme, we only assume $E_{1}(\phi)=\int_{\Omega}F(\phi)d\textbf{x}$ to be bounded from below, i.e., $E_{1}(\phi)\geq-C_{0}$. Similar to [26], we introduce a scalar auxiliary variable (SAV):

[TABLE]

Then, the nonlocal Cahn-Hilliard model (2.8) can be transformed into the following formulation:

[TABLE]

Taking the inner products of the above equations with $\mu$, $\phi_{t}$ and $2r$, respectively, one can obtain the modified energy dissipation law:

[TABLE]

where $(\mathcal{L}\phi,\phi)=2\int_{\Omega}\int_{\Omega}J(\textbf{x}-\textbf{y})\left[\phi(\textbf{x})-\phi(\textbf{y})\right]^{2}d\textbf{y}d\textbf{x}\geq 0$.

Let $N>0$ be a positive integer and set

[TABLE]

A semi-implicit first order SAV scheme for (3.1) reads as

[TABLE]

where $\tilde{\phi}^{n+1}$ is any explicit $O(\Delta t)$ approximation for $\phi(t^{n+1})$, which can be flexible according to the problem. For instance, we may use an extrapolation as follows

[TABLE]

Besides, we can also use the following simple first-order scheme to obtain it:

[TABLE]

Theorem 1**.**

The scheme (3.2)-(3.4) for Cahn-Hilliard type system is unconditionally energy stable in the sense that

[TABLE]

Proof.

By taking the inner products with $\Delta t\mu^{n+1}$, $\phi^{n+1}-\phi^{n}$, $2\Delta tr^{n+1}$ for equations (3.2), (3.3), (3.4) respectively and some simple calculations, we obtain

[TABLE]

Using the identity $\textbf{x}\cdot(\textbf{x}-\textbf{y})=\frac{1}{2}|\textbf{x}|^{2}-\frac{1}{2}|\textbf{y}|^{2}+\frac{1}{2}|\textbf{x}-\textbf{y}|^{2}$, the equation (3.8) can be transformed as follows:

[TABLE]

which states that $\frac{\epsilon^{2}}{2}(\mathcal{L}\phi^{n+1},\phi^{n+1})+|r^{n+1}|^{2}\leq\frac{\epsilon^{2}}{2}(\mathcal{L}\phi^{n},\phi^{n})+|r^{n}|^{2}$. ∎

A semi-implicit second order SAV/BDF scheme for (3.1) reads as

[TABLE]

where $\tilde{\phi}^{n+\frac{1}{2}}$ is any explicit $O(\Delta t^{2})$ approximation for $\phi(t^{n+1})$, which can be flexible according to the problem.

Theorem 2**.**

The scheme (3.10)-(3.12) for Cahn-Hilliard type system is unconditionally energy stable in the sense that

[TABLE]

Proof.

By taking the inner products with $2\Delta t\mu^{n+1}$, $3\phi^{n+1}-4\phi^{n}+\phi^{n-1}$, $2\Delta tr^{n+1}$ for (3.10), (3.11), (3.12) respectively and some simple calculations, we can obtain

[TABLE]

Using the identity

[TABLE]

the equation (3.14) can be transformed as follows:

[TABLE]

which completes the proof. ∎

4 The SAV finite difference scheme for nonlocal Cahn-Hilliard model

In this section, we consider finite difference discretization for nonlocal Cahn-Hilliard model for some spatial operators in the two dimensional space with $\Omega=(-L,L)\times(-R,R)$. We give the linear second order (in time) SAV/BDF finite difference scheme for nonlocal Cahn-Hilliard model. The first order SAV fully discrete scheme can be obtained straightforwardly.

In order to discretize the system (3.10)-(3.12), we define $\Omega_{h}=\{(x_{i},y_{j})|x_{i}=-L+ih_{x},~{}y_{j}=-R+jh_{y},~{}h_{x}=2L/M_{x},~{}h_{y}=2R/M_{y},~{}0\leq i\leq M_{x},~{}0\leq j\leq M_{y}\}$ to be a uniform mesh of the domain $\Omega$.

From the analysis of the nonlocal operator $\mathcal{L}$ in [13], one can see that

[TABLE]

Then, for any $v$, $\mathcal{L}v$ can be discreted at $(t^{n},x_{i},y_{j}),$ $(0\leq n\leq N,~{}0\leq i\leq M_{x},~{}0\leq j\leq M_{y})$ as follows:

[TABLE]

where

[TABLE]

and

[TABLE]

Combining the semi-implicit scheme (3.10)-(3.10) with equation (4.1), we obtain the finite difference discretization for nonlocal Cahn-Hilliard model (3.1) as follows:

[TABLE]

Lemma 3**.**

[Solvability of the SAV finite difference scheme] For any $\Delta t$, $h_{x}$, $h_{y}>0$, the scheme (4.2)-(4.4) has a unique solution.

Proof.

Denote

[TABLE]

Then, the scheme (4.2)-(4.4) can be rewritten as the following formulation:

[TABLE]

Obviously, the stiffness matrix $A=\left(I_{h}+\frac{2}{3}M\Delta t\epsilon^{2}\mathcal{L}_{h}^{2}\right)$ is positive definite. From the above equation, we observe that $(\eta^{n+1},\Phi^{n+1})$ needs to be computed first. Multiplying the above equation with $A^{-1}$, and taking the inner product with $\eta^{n+1}$, we obtain

[TABLE]

where $\theta=(A^{-1}\mathcal{L}_{h}\eta^{n+1},\eta^{n+1})\geq 0$.

Noting that $(1+\frac{1}{6}M\Delta t\theta)\neq 0$ and the stiffness matrix $A$ is positive definite, one can conclude that the scheme (4.2)-(4.4) has a unique solution. ∎

5 The fast solution method

From the discrete formulation (4.1), one can see that the nonlocal diffusion term will lead the stiffness matrix to be an almost full matrix which requires huge computational work and large memory. In this case, the fast solution method for solving derived linear system will become very important and necessary. In this section, we will analyse the structure of the stiffness matrix and seek some effective fast solution method to reduce the computational work and memory requirement. This fast solution technique is based on a fast Fourier transform and depends on the special structure of coefficient matrices.

Without loss of generality, we assume $\Omega=[-L,L]\times[-L,L]$ and the partition is uniform and satisfies $h=h_{x}=h_{y}$. For any $Au=f$ in solving the discrete formulation (4.2)-(4.4), we all use conjugate gradient method to obtain the solution of linear system: Let $u_{0}$ be an initial guess. Then, compute $r_{0}=f-Au_{0}$, $d_{1}=r_{0}$ and

[TABLE]

From the algorithm of conjugate gradient method, one can see that for reducing the huge computational work and memory requirement, we only need to find fast and efficient procedure to accelerate the matrix-vector multiplication $Ad$ for any vector $d$ and store $A$ efficiently.

Based on the above analysis, we note that the stiffness matrix $A=\left(I_{h}+\frac{2}{3}M\Delta t\epsilon^{2}\mathcal{L}_{h}^{2}\right)$. So, we only need to analyze the structure of the matrix-vector multiplication $\mathcal{L}_{h}d$.

First, we can rewrite the discrete formulation $(J\ast v)_{i,j}$ as follows

[TABLE]

Define $\gamma_{i,j}=(J\ast v)_{i,j}$, ${\bf{\Lambda}}_{i}=(\gamma_{i,0},\gamma_{i,1},\cdots,\gamma_{i,M})$ and ${\bf{\Lambda}}=({\bf{\Lambda}}_{1},{\bf{\Lambda}}_{2},\cdots,{\bf{\Lambda}}_{M})^{T}$. Then, define four transformation operators $\mathcal{A}_{1}$, $\mathcal{A}_{2}$, $\mathcal{A}_{3}$, $\mathcal{A}_{4}$. For any $(M+1)^{2}$-vector $v=(v_{0},v_{1},\ldots,v_{M})^{T},$ $v_{j}=(v_{0,j},v_{1,j},\ldots,x_{M,j})$, the operators $\mathcal{A}_{i}$, $(i=1,2,3,4)$ satisfy:

[TABLE]

where for $0\leq j\leq M-1$ and $1\leq k\leq M$,

[TABLE]

Define the following four vectors $\alpha_{M}$, $\beta_{0}$, $\zeta_{M}$ and $\theta_{0}$ to be all zero $(1\times M)$-vectors. Then, the equation (5.1) can be rewritten as the following formulation

[TABLE]

By the above equation, we can compute the vector ${\bf{\Lambda}}$ by the following equation:

[TABLE]

Then, for any $(M+1)^{2}$-vector $v$, we obtain the matrix-vector multiplication $\mathcal{L}_{h}v$ by

[TABLE]

Combining equation (5.1) and equation (5.3), and noting the four vectors in (5.2), the matrix B can be written as follows:

[TABLE]

Define $x_{i,j}=J(ih,jh)$. By simple calculation, the block matrix $\textbf{S}_{j}$ can be expressed as the following Toeplitz formulation:

[TABLE]

Therefore, we note that the matrix B is a Block-Toeplitz-Toeplitz-Block (BTTB) matrix. Then, we can use fast Fourier transform (FFT) method to evaluate $\textbf{B}G$ for $G=f_{1}+f_{2}+f_{3}+f_{4}$ which can reduce the computational work and memory requirement effectively.

Firstly, the block Toeplitz matrix $\textbf{S}_{j}$ can be embedded into a $(2M+2)\times(2M+2)$ circulant matrix:

[TABLE]

Then, replacing the block matrix $\textbf{S}_{j}$ with the block circulant matrix $\textbf{C}_{j}$ in the BTTB matrix B, we obtain a Bolck-Toeplitz-Circulant-Block (BTCB) matrix $\widehat{\textbf{B}}$

[TABLE]

The BTCB matrix $\widehat{\textbf{B}}$ can be embedded into a $(2M+2)\times(2M+2)$ Bolck-Circulant-Circulant-Block (BCCB) matrix D as follows

[TABLE]

The BCCB matrix D has the following decomposition

[TABLE]

where $d$ is the first column vector of D and $F\otimes F$ is the two dimensional discrete Fourier transform matrix. Then, it is well known that the matrix-vector multiplication $Fw$ for $w\in\mathbb{R}^{(2M+2)^{2}}$ can be carried out in $O(M^{2}$log$M^{2})$ operations via the fast Fourier transform (FFT). Equation (5.4) shows that $\textbf{D}w$ can be evaluated in $O(M^{2}$log$M^{2})$ operations. Define $N=(M+1)^{2}$ and for $G=(g_{0},g_{1},\ldots,g_{M})^{T}$, define $\widehat{G}=(g_{0},0,g_{1},0,\ldots,g_{M},0)$. Then, we obtain that $\textbf{B}G$ can be evaluated in $O(N$log$N)$ operations for any $G\in\mathbb{R}^{N}$ by evaluating $\textbf{D}\widetilde{G}$ with FFT where $\widetilde{G}=(\widehat{G},0)$. The overall computational cost of the fast conjugate gradient method is $O(N$log${}^{2}N)$, since the number of iterations is $O($log$N)$. What we need to focus on is that if one uses the Gaussian elimination method straightforwardly to this linear system, then it requires $O(N^{3})$ complexity. In addition, since $N\times N$ BTTB matrix is determined by only $2N-1$ entries rather than $N^{2}$ entries, the fast solver will reduce memory requirement from $O(N^{2})$ to $O(N)$.

6 Numerical experiments

In this section, we present some numerical examples for the nonlocal Cahn-Hilliard equation in two dimension to test our theoretical analysis which contains energy stability and convergence rates of the proposed numerical schemes. We use the finite difference method for spatial discretization for all numerical examples. In all examples, we set the domain $\Omega=(-1,1)\times(-1,1)$. All the solvers are implemented using Matlab and all the numerical experiments are performed on a computer with 8-GB memory.

The Gaussian kernel $J_{\delta}$ will be given below [13]:

[TABLE]

From [13], one can see that for any $v\in C^{\infty}(\Omega)$, $\textbf{x}\in\Omega$, $\mathcal{L}_{\delta}v(\textbf{x})\rightarrow-\Delta v(\textbf{x})$ as $\delta\rightarrow 0$. It tells us that the nonlocal Cahn-Hilliard model equipped with above Gaussian kernel converges to the classical local Cahn-Hilliard model as $\delta\rightarrow 0$.

We first give an example to test convergence rates of the proposed schemes (first order: SCHEM1 and second order:SCHEM2 in time) for the nonlocal Cahn-Hilliard equation in two dimension and check the efficiency of our fast procedure.

Example 1: Consider the nonlocal Cahn-Hilliard equation with $\epsilon^{2}=0.1$, $M=1$, $\delta=\epsilon$, $T=0.05$ and the following initial condition:

[TABLE]

To observe the temporal convergence rate, we first calculate the reference solution $\phi_{ref}$ with $h_{ref}=0.01$ and $\Delta t_{ref}=0.05\times 2^{-14}$ since the exact solution is not known. Then, we use both direct solver (which means using the self-contained function $x=A\setminus b$ in Matlab R2015a) and fast solver (which means fast CG solver in Section 4) to obtain approximate solution with the time step sizes $\Delta t=0.05\times 2^{-4}$, $0.05\times 2^{-5}$, $0.05\times 2^{-6}$, $0.05\times 2^{-7}$, $0.05\times 2^{-8}$, $0.05\times 2^{-9}$ and $h=0.01$. Tables 1 and 2 show the discrete $L_{2}$ errors and the temporal convergence rates of numerical solutions and the CPU time costs for both direct and fast solvers. One can see that with the same mesh and time step sizes, both direct and fast solvers generate numerical solutions with almost the same errors and convergence rates. Furthermore, we again observe that fast solver saves more CPU time than direct solver to obtain same accuracy.

To observe spatial convergence rates, we calculate the reference solution $\phi_{ref}$ with $h_{ref}=2^{-10}$ and $\Delta t_{ref}=5e-5$. Then, we use $h=2^{-3},\ 2^{-4},\ 2^{-5},\ 2^{-6},\ 2^{-7},\ 2^{-8}$, and $\Delta t=5e-5$ to obtain the approximate solution at $t=T$ for direct solver and fast solver. Tables 3 and 4 show the discrete $L_{2}$ errors, the spatial convergence rates of numerical solutions and CPU time costs for both direct and fast solvers. It can be observed that the spatial error are almost $O(h^{2})$ for both SCHEM1 and SCHEM2. One can also obtain that the fast solver has almost the same error and convergence rate as direct solver. In addition, we again observe that fast solver saves more CPU time and memory requirement than direct solver to obtain similar accuracy. Thus, we can realize the numerical simulation under the finer meshes more quickly through the fast algorithm.

In the following examples, we study the phase separation behavior using the second order scheme SCHEM2 and obtain the numerical solution by fast solver.

Example 2: In the following, we take $\epsilon=0.02$, $M=1$, $\delta=0.02$. The initial condition is chosen as

[TABLE]

with the radius $R_{0}=0.36$, $(x_{1},y_{1})=(0.4,0)$ and $(x_{2},y_{2})=(-0.4,0)$. Initially, two bubbles, centered at $(0.4,0)$ and $(-0.4,0)$, respectively, are osculating or ”kissing”. In the simulation, we choose the mesh size $h=1/100$ and the time step $\Delta t=10^{-3}$. The process coalescence of two bubbles is demonstrated in Figure 1. Snapshots of the phase variable $\phi$ are taken at $t=0$, $0.05$, $0.1$, $1$, $5$, $10$ in Figure 1. In Figure 1, we show the evolutions of the phase field variable $\phi$ at various time by using the time step $\Delta t=1e-3$. We observe the coarsening effect that the small circle is absorbed into the big circle, and the total absorption happens at around $t=10$ which is consistent with classical Cahn-Hilliard model.

Example 3: We give the following example for nonlocal Cahn-Hilliard equation with $\epsilon=0.02$, $M=1$. The initial condition is

[TABLE]

where the $rand(x,y)$ is the random number in $[-1,1]$ with zero mean.

In Figure 2, we perform the simulations by using the time step $\Delta t=1e-3$. The figure shows the dynamical behaviors of the phase separation for the random initial value (6.3). In Figure 2, the snapshots of coarsening dynamics are taken at $t=0$, $0.1$, $0.5$, $1$, $2$, $10$ and the final steady shape forms several big drops. Next, we fix $\delta=0.05$ with $h=0.01$ and decrease $\epsilon$ from $0.1$ to $0.02$. The energies are plotted in Figure 3. One can find that the energy decay rates comply with the $-1/3$ power law quite well for all cases which is consistent to the classical local Cahn-Hilliard model.

7 Conclusion

In this paper, we develop accurate and efficient linear algorithms for the general nonlocal Cahn-Hilliard equation with general nonlinear potential and prove the unconditional energy stability for its semi-discrete schemes carefully and rigorously. we construct and analyze linear, first and second order (in time) numerical scalar auxiliary variable approaches to construct unconditionally energy stable schemes. In addition, considering the huge computational work and memory requirement in solving the linear system, we analyse the structure of the stiffness matrix and seek some effective fast solution method to reduce the computational work and memory requirement. This fast solution technique is based on a fast Fourier transform and depends on the special structure of coefficient matrices. In the future work, error estimates for the fully discrete schemes will be investigated.