Arbitrarily High-order Unconditionally Energy Stable Schemes for Thermodynamically Consistent Gradient Flow Models

TL;DR

This paper develops high-order, unconditionally energy stable numerical schemes for gradient flow models using the energy quadratization method, enabling accurate and stable simulations of thermodynamically consistent PDEs.

Contribution

It introduces two novel classes of energy stable schemes for gradient flows, employing prediction-correction and Gaussian collocation methods, with rigorous stability proofs and broad applicability.

Findings

Schemes are unconditionally energy stable.

Numerical experiments confirm high accuracy.

Methods are generalizable to other PDE models.

Abstract

We present a systematical approach to developing arbitrarily high order, unconditionally energy stable numerical schemes for thermodynamically consistent gradient flow models that satisfy energy dissipation laws. Utilizing the energy quadratization (EQ) method, We formulate the gradient flow model into an equivalent form with a corresponding quadratic free energy functional. Based on the equivalent form with a quadratic energy, we propose two classes of energy stable numerical approximations. In the first approach, we use a prediction-correction strategy to improve the accuracy of linear numerical schemes. In the second approach, we adopt the Gaussian collocation method to discretize the equivalent form with a quadratic energy, arriving at an arbitrarily high-order scheme for gradient flow models. Schemes derived using both approaches are proved rigorously to be unconditionally energy stable. The proposed schemes are then implemented in four gradient flow models numerically to demonstrate their accuracy and effectiveness. Detailed numerical comparisons among these schemes are carried out as well. These numerical strategies are rather general so that they can be readily generalized to solve any thermodynamically consistent PDE models.

Equations

$$\frac{\partial }{\partial t}\Phi =\mathcal{G}\frac{\delta F}{\delta \Phi },$$

$$\frac{dF}{dt}=\left(\frac{\delta F}{\delta \Phi },\frac{\partial \Phi }{\partial t}\right)=\left(\frac{\delta F}{\delta \Phi },\mathcal{G}\frac{\delta F}{\delta \Phi }\right)\le 0,$$

$$F=\frac{1}{2}(\Phi,\mathcal{L}\Phi)+\big{(}f(\Phi),1\big{)}.$$

$$\mathcal{F}=\frac{1}{2}(\Phi ,\mathcal{L}\Phi )+\frac{1}{2}\Vert q{\Vert }^2-A,$$

$$F=\frac{{\epsilon }^2}{2}\Vert \nabla \varphi {\Vert }^2+\frac{1}{4}\Vert {\varphi }^2-1{\Vert }^2,$$

$$F=\frac{1}{2}\Big{(}-\varepsilon^{2}\Delta\phi+\gamma_{0}\phi,\phi\Big{)}+\frac{1}{4}\|\phi^{2}-1-\gamma_{0}\|^{2}+\frac{\varepsilon^{2}}{2}\int_{\partial\Omega}\phi\nabla\phi\cdot\mathbf{n}dS-\left(\frac{\gamma_{0}}{2}+\frac{\gamma_{0}^{2}}{4}\right)|\Omega|,$$

$$\left\{\begin{array}[]{l}\frac{\partial}{\partial t}\Phi=\mathcal{G}\Big{(}\mathcal{L}\Phi+q\frac{\partial g}{\partial\Phi}\Big{)},\\ \frac{\partial}{\partial t}q=\frac{\partial g}{\partial\Phi}\cdot\frac{\partial\Phi}{\partial t},\end{array}\right.$$

$$\frac{\partial }{\partial t}\Psi =\mathcal{N}(\Psi )\mathcal{B}\Psi ,$$

$$\mathcal{A}={\left({\mathbf{I}}_s\frac{\partial g}{\partial \Phi }\right)}_{s\times (s+1)},\mathcal{B}=\text{diag}(\mathcal{L},1{)}_{(s+1)\times (s+1)}.$$

$$\Vert \Psi {\Vert }_{\mathcal{B}}=\sqrt{(\Psi ,\mathcal{B}\Psi )},$$

$$\frac{d\mathcal{F}}{dt}=\left(\frac{\delta \mathcal{F}}{\delta \Psi },\frac{\partial \Psi }{\partial t}\right)=\left(\mathcal{B}\Psi ,\mathcal{N}(\Psi )\mathcal{B}\Psi \right)\le 0.$$

$$\left\{\begin{array}[]{l}\frac{\partial}{\partial t}\Phi=\mathcal{G}\Big{(}\mathcal{L}\Phi+q\frac{\partial g}{\partial\Phi}-\nabla\cdot(q\frac{\partial g}{\partial\nabla\Phi})\Big{)},\\ \frac{\partial}{\partial t}q=\frac{\partial g}{\partial\Phi}\cdot\frac{\partial\Phi}{\partial t}+\frac{\partial g}{\partial\nabla\Phi}\cdot\nabla\frac{\partial\Phi}{\partial t}.\end{array}\right.$$

$$\mathcal{A}=\left(\mathbf{I}_{s}\quad\frac{\partial g}{\partial\Phi}-\nabla\cdot\frac{\partial g}{\partial\nabla\Phi}\right),\quad\mathcal{A}^{*}=\left(\begin{array}[]{l}\mathbf{I}_{s}\\ \Big{(}\frac{\partial g}{\partial\Phi}+\frac{\partial g}{\partial\nabla\Phi}\cdot\nabla\Big{)}^{T}\end{array}\right).$$

$$\frac{\Psi^{n+1}-\Psi^{n}}{\Delta t}=\mathcal{N}\Big{(}\frac{\Psi^{n+1}_{*}+\Psi^{n}}{2}\Big{)}\mathcal{B}\frac{\Psi^{n+1}+\Psi^{n}}{2};$$

$$\frac{3{\Psi }^{n+1}-4{\Psi }^n+{\Psi }^{n-1}}{2\Delta t}=\mathcal{N}({\Psi }_{\ast }^{n+1})\mathcal{B}{\Psi }^{n+1}.$$

$${\Psi }_{\ast }^{n+1}=2{\Psi }^n-{\Psi }^{n-1},$$

$$\frac{\Psi^{n+1}_{i+1}-\Psi^{n}}{\Delta t}=\mathcal{N}\Big{(}\frac{\Psi^{n+1}_{i}+\Psi^{n}}{2}\Big{)}\mathcal{B}\frac{\Psi^{n+1}_{i+1}+\Psi^{n}}{2};$$

$$\frac{3{\Psi }_{i+1}^{n+1}-4{\Psi }^n+{\Psi }^{n-1}}{2\Delta t}=\mathcal{N}({\Psi }_i^{n+1})\mathcal{B}{\Psi }_{i+1}^{n+1}.$$

$$\frac{\Psi^{n+1}_{i+1}-\Psi^{n}}{\Delta t}=\mathcal{N}_{1}\mathcal{B}\frac{\Psi^{n+1}_{i+1}+\Psi^{n}}{2}+\mathcal{N}_{2}\Big{(}\frac{\Psi^{n+1}_{i}+\Psi^{n}}{2}\Big{)}\mathcal{B}\frac{\Psi^{n+1}_{i}+\Psi^{n}}{2};$$

$$\frac{3{\Psi }_{i+1}^{n+1}-4{\Psi }^n+{\Psi }^{n-1}}{2\Delta t}={\mathcal{N}}_1\mathcal{B}{\Psi }_{i+1}^{n+1}+{\mathcal{N}}_2({\Psi }_i^{n+1})\mathcal{B}{\Psi }_i^{n+1}.$$

$$\frac{\Psi^{1}-\Psi^{0}}{\Delta t}=\mathcal{N}\Big{(}\frac{\Psi^{1}_{*}+\Psi^{0}}{2}\Big{)}\mathcal{B}\frac{\Psi^{1}+\Psi^{0}}{2},$$

$$\frac{{\Psi }_{\ast }^1-{\Psi }^0}{\Delta t}={\mathcal{N}}_1\mathcal{B}\frac{{\Psi }_{\ast }^1+{\Psi }^0}{2}+{\mathcal{N}}_2({\Psi }^0)\mathcal{B}{\Psi }^0.$$

$$\frac{F^{n+1}-F^{n}}{\Delta t}=\left(\mathcal{B}\Psi^{n+\frac{1}{2}},\mathcal{N}\Big{(}\frac{\Psi^{n+1}_{*}+\Psi^{n}}{2}\Big{)}\mathcal{B}\Psi^{n+\frac{1}{2}}\right)\leq 0,~{}\forall n\geq 1,$$

$$\frac{F^{n+\frac{3}{2}}-F^{n+\frac{1}{2}}+\widetilde{F}^{n+1}}{\Delta t}=\Big{(}\mathcal{B}\Psi^{n+1},\mathcal{N}(\Psi_{*}^{n+1})\mathcal{B}\Psi^{n+1}\Big{)}\leq 0,~{}\forall n\geq 1,$$

$${\Lambda }_k{\Psi }^{n+1}=\sum _{i=0}^k{\lambda }_i^k{\Psi }^{n+1-i},$$

$$\frac{1}{\Delta t}{\Lambda }_k{\Psi }^{n+1}=\mathcal{N}({\overline{\Psi }}^{n+1})\mathcal{B}{\Psi }^{n+1},$$

$$\frac{\partial }{\partial t}\Psi =\mathcal{N}(\Psi )\mathcal{B}\Psi .$$

$$\begin{array}[]{l}k_{i}=\mathcal{N}\Big{(}\Psi^{n}+\Delta t\sum\limits_{j=1}^{s}a_{ij}k_{j}\Big{)}\mathcal{B}\Big{(}\Psi^{n}+\Delta t\sum\limits_{j=1}^{s}a_{ij}k_{j}\Big{)},\quad i=1,\cdots,s,\\ \Psi^{n+1}=\Psi^{n}+\Delta t\sum\limits_{i=1}^{s}b_{i}k_{i}.\end{array}$$

$$\begin{array}[]{c|c}\mathbf{c}&\mathbf{A}\\ \hline\cr&\mathbf{b}^{T}\\ \end{array}=\begin{array}[]{c|ccc}c_{1}&a_{11}&\cdots&a_{1s}\\ \vdots&\vdots&&\vdots\\ c_{s}&a_{s1}&\cdots&a_{ss}\\ \hline\cr&b_{1}&\cdots&b_{s}\\ \end{array},$$

$$\begin{array}[]{l}u(t_{n})=\Psi^{n},\\ \\ \dot{u}(t_{n}+c_{i}\Delta t)=\mathcal{N}\big{(}u(t_{n}+c_{i}\Delta t)\big{)}\mathcal{B}\big{(}u(t_{n}+c_{i}\Delta t)\big{)},\quad i=1,\cdots,s,\end{array}$$

Full text

1 Introduction

In nonequilibrium thermodynamics, transient dynamics is stipulated primarily following the ”linear response theory”, collectively known as the Onsager principle, or equivalently the second law of thermodynamics [38, 36, 54]. The second law of thermodynamics has been used to develop non-equilibrium models for a wide range of material systems, ranging from life science, materials science, industrial processing, and engineering. The gradient flow model describing relaxation dynamics is one of the well-known examples, and other examples include most differential constitutive laws for complex fluid flows [2]. The models derived from the second law of thermodynamics are referred to as the thermodynamically consistent models.

For a given thermodynamically consistent, gradient flow model in the form of partial differential equations (PDEs), a computationally efficient, energy-dissipation-rate preserving, and high order numerical approximation is always desired. In the literature, the numerical scheme that preserves the energy dissipation property is known as the energy stable scheme [13]. If the energy-dissipation-rate preserving property doesn’t depend on the time-step size, the schemes are called unconditionally energy stable.

In the past, two widely used, distinct strategies for developing (unconditionally) energy stable schemes have been proposed, which are the convex splitting method [13, 49, 39, 48, 31, 51, 6] and the stabilizing method [58, 35, 43, 45, 44, 59, 17]. The convex-splitting strategy relies on the existence of a pair of convex components that give rise to the free energy functional as a difference of the two convex functionals. If such a splitting exists in the free energy functional, a nonlinear scheme can be devised to render an unconditionally energy stable scheme. The convex splitting schemes have been widely applied to the gradient flow models [28, 13, 49, 39, 48, 31, 52, 14, 8, 16, 7]. The stabilizing approach augments discretized equations by high order terms to turn the scheme into an energy stable scheme. Usually, this is accomplished by adding additional dissipation to the numerical scheme. Both strategies can yield a dissipative scheme, but do not guarantee to preserve the energy dissipation rate. Some other related work include [31, 19, 20, 15, 29, 32, 27, 26] .

Recently, Badia, Guillen-Gonzales, Gutierres-Santacreu, and Tierra pioneered a new idea of transforming the free energy into a quadratic functional to derive energy stable schemes [1, 27]. This is amplified and systematically applied to many specific thermodynamic and hydrodynamic models by Yang, Zhao, Shen and Wang, etc. [53, 55, 57, 59, 60, 23, 25]. Yang, Zhao, and Wang coined the name Invariant Energy Quadratization (IEQ) method for this class of methods. Later, we abbreviated the name to simply Energy Quadratization (EQ) method, which is more appropriate. This strategy bypasses the complications in the other methods to derive linear, second-order energy stable schemes in time readily. It is so general that the EQ approach has little restriction on the specific expression of the free energy functional. Recently, Shen et al. [42, 41] implemented the idea of EQ using a scalar auxiliary variable and called it the SAV method. However, using either EQ or SAV strategy, one has only designed and proved rigorously unconditionally energy stable numerical approximations up to second-order in time so far. This may not be sufficient for some gradient flow problems with long time simulations or sharp transition dynamics. Recently in [46], the authors propose a high order Runge-Kutta (RK) method for gradient flow problems based on convex splitting methods, but it requires too many stages even to reach 3rd order accuracy, and it does not work for gradient flows with variable mobilities.

In this paper, we develop systematically two classes of numerical approximations exploiting the EQ approach for thermodynamically consistent gradient flow models. In the first strategy, we introduce a prediction step to correct the extrapolation values (in the explicit terms of the linear schemes resulted from the EQ approach). Specifically, instead of using extrapolations to obtain explicit terms at the desired time level [55, 57, 59, 60, 23, 25], we use a fixed-point iteration to predict and correct the terms using values from previous iteration steps. In numerical experiments, we show that this strategy can improve the accuracy of numerical schemes significantly. In the second strategy, we extend the second-order EQ method [55, 57, 59, 60, 23, 25] into an arbitrarily high-order energy quadratization (HEQ) method. The HEQ strategy consists of two parts. Firstly, by introducing auxiliary variables, we transform the original gradient flow model into a gradient flow model with an energy functional consisting of only quadratic terms of the unknown variables. The resulting model is referred to as the quadratized gradient flow model. Secondly, we apply the Gaussian collocation method to the quadratized gradient flow model to produce unconditionally energy stable numerical schemes. It turns out that the classical Crank-Nicolson and backward differential formula schemes are special cases of the newly proposed general scheme. Our new numerical strategy provides an elegant solution for developing arbitrarily high order and unconditionally energy stable numerical schemes for gradient flow models. Moreover, the schemes preserve the energy dissipation rate in the transformed variables accurately. Due to their high-order accuracy and unconditional energy stability, these schemes can allow large time steps, making the numerical approximations especially appealing for long time computations.

We organize the rest of the paper as follows. In §2, we reformulate gradient flow models by using the EQ method. In §3, we present two classes of unconditionally energy stable, numerical approximations for the reformulated EQ model. In §4, the Fourier pseudospectral method is employed to give rise to the spatial discretization. Then four numerical examples of gradient flow models are shown to validate the efficiency and accuracy of our proposed schemes in §5. Finally, we give a concluding remark in the last section.

2 Gradient Flow Models and Their EQ Reformulation

We present the general thermodynamically consistent, gradient flow model firstly. Then, we reformulate the general gradient flow system into an equivalent form with a quadratic energy functional using the energy quadratization technique, called EQ reformulation. The EQ reformulation for this class of gradient flow models provides an elegant platform for developing arbitrarily high-order unconditionally energy stable schemes. In this paper, we adopt periodic boundary conditions for simplicity. The results can be readily applied to gradient flow problems with physical boundary conditions so long as the spatial discretization respects the integration-by-parts or summation-by-parts formula.

2.1 Gradient flow models

Mathematically, the general form of the governing system of equations of a gradient flow model is given by [61, 41]

[TABLE]

where $\Phi=(\phi_{1},\cdots,\phi_{s})^{T}$ are the state variables, $\mathcal{G}$ is an $s$-by-$s$ mobility matrix operator which is negative semi-definite and may depend on $\Phi$. Here $F$ is the effective free energy of the material system, and $\frac{\delta F}{\delta\Phi}$ is the variational derivative of the free energy functional with respect to the state variables, known as the chemical potential. Then, the triple $(\Phi,\mathcal{G},F)$ uniquely defines a gradient flow model. One intrinsic property of (2.1) is the energy dissipation law

[TABLE]

where the inner product is defined by $({\bf f},{\bf g})=\sum\limits_{i=1}^{s}\int_{\Omega}f_{i}g_{i}d\mathbf{x}$, $\forall\mathbf{f},\mathbf{g}\in(L^{2}(\Omega))^{s}$, and $\Omega$ is the material domain. Note that the energy dissipation law (2.2) holds only for suitable boundary conditions. These boundary conditions include the periodic boundary conditions and the boundary conditions that make the boundary integrals resulted during the integration by parts vanish.

2.2 Model reformulation using the EQ approach

We reformulate the general gradient flow model (2.1) by first transforming the free energy into a quadratic form. We illustrate the idea using a simple case, where the free energy is given by

[TABLE]

Here $\mathcal{L}$ is a linear, self-adjoint, positive definite operator and $f$ is the bulk part of the free energy density, which is bounded from below for physically accessible state of $\Phi$. Then the free energy $F$ can be rewritten into

[TABLE]

where q=\sqrt{2\Big{(}f(\Phi)+\frac{A}{|\Omega|}\Big{)}}, and $A$ is a constant large enough to make $q$ well-defined. Here $\|\bullet\|$ represents the $L^{2}$ norm, i.e. $\|f\|=\sqrt{\int_{\Omega}f^{2}d\mathbf{x}}$, $\forall f\in L^{2}(\Omega)$.

For instance, given a Ginzburg-Landau free energy

[TABLE]

if we conduct an integration by part once, it converts to

[TABLE]

we identify $\mathcal{L}=-\varepsilon^{2}\Delta+\gamma_{0}$ and $q=\frac{1}{\sqrt{2}}(\phi^{2}-1-\gamma_{0})$, assuming the boundary integral term vanishes.

Denote $g(\Phi)=\sqrt{2\left(f(\Phi)+\frac{A}{|\Omega|}\right)}$. Then we reformulate model (2.1) to an equivalent form

[TABLE]

where $\mathbf{a}\cdot\mathbf{b}=\sum\limits_{i=1}^{s}a_{i}b_{i}.$ Letting \Psi=\left(\begin{array}[]{l}\Phi\\ q\end{array}\right), system (2.7) can be written in the following compact form

[TABLE]

where $\mathcal{N}(\Psi)=\mathcal{A}^{*}\mathcal{G}\mathcal{A}$ is an $(s+1)\times(s+1)$ matrix operator depending on $\Psi$, $\mathcal{A}^{*}$ is the adjoint operator of $\mathcal{A},$ and

[TABLE]

In this case, we have $\mathcal{A}^{*}=\mathcal{A}^{T}.$ Since $\mathcal{G}$ is negative semi-definite, it can be shown easily that $\mathcal{N}(\Psi)$ is also negative semi-definite for any $\Psi$. In addition, $\mathcal{B}$ is a linear, self-adjoint, positive definite operator thanks to the property of $\mathcal{L}$. Define the $\mathcal{B}$-norm as

[TABLE]

then the free energy (2.4) is rewritten as $\mathcal{F}=\frac{1}{2}\|\Psi\|_{\mathcal{B}}^{2}-A.$ System (2.8) preserves the following energy dissipation law

[TABLE]

The original gradient flow system defined by $(\Phi,\mathcal{G},F)$ is transformed into a new system $(\Psi,\mathcal{N},\mathcal{F}),$ where the free energy $\mathcal{F}$ in the new system is regarded as a quadratic functional of $\Psi$. Therefore, this process is called energy quadratization reformulation.

Remark 2.1**.**

More generally, if $q=g(\Phi,\nabla\Phi)=\sqrt{2\left(f(\Phi,\nabla\Phi)+\frac{A}{|\Omega|}\right)}$, then the equivalent EQ form is given by

[TABLE]

The system can again be written into the form (2.8) with the operators $\mathcal{A}$ and $\mathcal{A}^{*}$ given by

[TABLE]

Remark 2.2**.**

*In general, the EQ reformulation approach can be applied to any gradient flow models with free energies of high-order spatial derivatives so long as they are thermodynamically consistent and energy dissipative. *

We next discuss how to design accurate and energy stable schemes for gradient flow models. We will demonstrate that the equivalent form in (2.8) can be handled more easily than its original system given in (2.1).

3 Temporal Discretization

In this section, a class of linear second-order prediction-correction schemes and a class of arbitrarily high-order Gaussian collocation schemes are proposed respectively, where all the schemes are shown to be unconditionally energy stable, i.e., the energy dissipation property is conserved for any time step sizes at the semi-discrete level.

3.1 Linear energy stable schemes

As we know, the linear-implicit Crank-Nicolson (LCN) scheme and the linear-implicit second-order backward differentiation/extrapolation (LBDF2) method can be applied directly for discretizing the reformulated EQ system (2.8) in time to obtain linear unconditionally energy stable schemes [23, 56, 27, 57]. Even though both schemes are second-order accurate, the rate of convergence usually can be reached when the time step is small. Here, we propose a new class of prediction-correction schemes motivated by the works in [40, 18, 37]. Employing the prediction-correction strategy for the CN scheme or the BDF2 scheme, we obtain the following prediction-correction schemes:

Scheme 3.1** (Linear Prediction-Correction Scheme).**

Given $\Psi^{n-1}$ and $\Psi^{n}$, $\forall n\geq 1$, we obtain $\Psi^{n+1}$ through the following two steps:

Prediction: predict $\Psi^{n+1}_{*}$ via some efficient and at least second order accurate numerical schemes. 2. 2.

CN correction:

[TABLE]

or BDF2 correction:

[TABLE]

For each time step, several prediction strategies can be devised so long as the predicted value $\Psi^{n+1}_{*}$ is at least second-order consistent in time. As an illustration, we list some choices for the prediction step below.

$\bullet$ Case 1: if we set

[TABLE]

then scheme (3.1) and (3.2) reduce to LCN and LBDF2 schemes studied in [23], respectively.

$\bullet$ Case 2: we set $\Psi^{n+1}_{0}=2\Psi^{n}-\Psi^{n-1}.$ For $i=0$ to $N-1$, we compute $\Psi^{n+1}_{i+1}$ using

[TABLE]

or

[TABLE]

If $\|\Psi^{n+1}_{i+1}-\Psi^{n+1}_{i}\|_{\infty}<\varepsilon_{0}$ and $i+1<N$, we stop the iteration and set $\Psi^{n+1}_{*}=\Psi^{n+1}_{i+1}$; otherwise, we set $\Psi^{n+1}_{*}=\Psi^{n+1}_{N}.$

$\bullet$ Case 3: we rewrite $\mathcal{N}(\Psi)$ into the sum of $\mathcal{N}_{1}$, a linear operator of constant coefficients and a nonlinear operator $\mathcal{N}_{2}(\Psi)$, $\mathcal{N}(\Psi)=\mathcal{N}_{1}+\mathcal{N}_{2}(\Psi)$. Set $\Psi^{n+1}_{0}=2\Psi^{n}-\Psi^{n-1}.$ For $i=0$ to $N-1$, we compute $\Psi^{n+1}_{i+1}$ by

[TABLE]

or

[TABLE]

If $\|\Psi^{n+1}_{i+1}-\Psi^{n+1}_{i}\|_{\infty}<\varepsilon_{0}$ and $i+1<N$, we stop the iteration and set $\Psi^{n+1}_{*}=\Psi^{n+1}_{i+1}$; otherwise, we set $\Psi^{n+1}_{*}=\Psi^{n+1}_{N}.$

Remark 3.1**.**

*In some practical implementations, Case 3 is more efficient than Case 2 because $\mathcal{N}_{1}$ is a linear operator of constant coefficients so that the Fast Fourier transform (FFT) can be readily applied to the linear part of the scheme. *

Remark 3.2**.**

*If $N$ is large enough, the prediction-correction scheme (3.1) with (3.4) (or (3.6)) approximates to the fully implicit CN scheme while the prediction-correction scheme (3.2) with (3.5) (or (3.7)) approximates to the traditional BDF2 scheme. There is no theoretical result on the choice of iteration step $N$. From our numerical tests, several iteration steps $N\leq 5$ would improve the accuracy noticeably. *

Remark 3.3**.**

In Scheme 3.1, the initial second level datum $\Psi^{1}$ is computed by

[TABLE]

where $\Psi^{1}_{*}$ is given by

[TABLE]

Similar to the linear unconditionally energy stable schemes in [23], we have the following theorems.

Theorem 3.1**.**

Scheme (3.1) satisfies the following energy dissipation law

[TABLE]

*with $F^{n}=\frac{1}{2}\|\Psi^{n}\|_{\mathcal{B}}^{2}-A$ and $\Psi^{n+\frac{1}{2}}=(\Psi^{n+1}+\Psi^{n})/2.$ *

Theorem 3.2**.**

Scheme (3.2) satisfies the following energy identity

[TABLE]

*with F^{n+\frac{1}{2}}=\frac{1}{4}\Big{(}\|\Psi^{n}\|_{\mathcal{B}}^{2}+\|2\Psi^{n}-\Psi^{n-1}\|_{\mathcal{B}}^{2}\Big{)}-A and $\widetilde{F}^{n+1}=\frac{1}{4}\|\Psi^{n+1}-2\Psi^{n}+\Psi^{n-1}\|_{\mathcal{B}}^{2}.$ *

Remark 3.4**.**

*Eqs. (3.10) and (3.11) imply that the two schemes given by (3.1) and (3.2) are unconditionally energy stable, i.e. they possess the energy decay property, respectively, $F^{n+1}\leq F^{n},~{}\forall n\geq 1$ and $F^{n+\frac{3}{2}}\leq F^{n+\frac{1}{2}},~{}\forall n\geq 1.$ Therefore, they are unconditionally energy stable. In addition, Scheme 3.1 preserves the energy dissipation rate. *

3.2 Arbitrarily high-order, unconditionally energy stable schemes

The linear prediction-correction semi-discrete schemes discussed above are at most second-order accurate in time. To derive arbitrarily high-order linear schemes, one can combine high-order backward differentiation formula with the matched extrapolation. Define the backward difference operator for the $k$th order derivative as follows:

[TABLE]

where $\Lambda_{0}\Psi^{n+1}:=\Psi^{n+1}$. The values of $\{\lambda_{i}^{k}\}$ could be found in [30]. Then, we propose the linear-implicit BDF-$k$ scheme (where $k$ is the order of the scheme)

Scheme 3.2** (BDF-$k$ Scheme).**

Given $\Psi^{n+1-k},\Psi^{n+2-k},\cdots,\Psi^{n}$, we obtain $\Psi^{n+1}$ using

[TABLE]

*where $\overline{(\bullet)}^{n+1}$ is a matched extrapolation with values from previous time steps. *

Unfortunately, we are not able to prove energy stability for Scheme 3.2 presently, although in practice, Scheme 3.2 is usually shown to deliver an energy decay numerical result with reasonable time steps.

Given that we can’t prove energy stability for the high order schemes obtained using the BDF method in time, we turn to another time discretization strategy for developing energy-stable schemes in time. Starting from the reformulated EQ system (2.8), we apply the Gaussian collocation methods to construct arbitrarily high-order schemes in time. Then, we can prove the obtained schemes are unconditionally energy stable rigorously.

Recall the energy-quadratized system (2.8), reformulated from the general gradient flow model (2.1), as follows

[TABLE]

Firstly, we briefly recall the Runge-Kutta (RK) and collocation method (see Chapter II of [30] for detailed discussions). Applying an $s$-stage RK method to solve equation (3.14), we obtain the following high-order energy quadratization (HEQ) scheme.

Scheme 3.3** ($s$-stage HEQ-RK Method).**

Let $b_{i}$, $a_{ij}$ ($i,j=1,\cdots,s$) be real numbers and let $c_{i}=\sum\limits_{j=1}^{s}a_{ij}$. Given $\Psi^{n}$, $\Psi^{n+1}$ is calculated by

[TABLE]

The coefficients are given by a Butcher table

[TABLE]

where $\mathbf{A}\in\mathbb{R}^{s,s}$, $\mathbf{b}\in\mathbb{R}^{s}$, and $\mathbf{c}=\mathbf{A}\mathbf{l}$, with $\mathbf{l}=(1,1,\cdots,1)^{T}\in\mathbb{R}^{s}$. Applying an $s$-stage collocation method to (3.14), we obtain the following scheme.

Scheme 3.4** ($s$-stage HEQ Collocation Method).**

Let $c_{1},\cdots,c_{s}$ be distinct real numbers ($0\leq c_{i}\leq 1$). Given $\Psi^{n}$, the collocation polynomial $u(t)$ is a polynomial of degree $s$ satisfying

[TABLE]

*and the numerical solution is defined by $\Psi^{n+1}=u(t_{n}+\Delta t)$. *

Theorem 1.4 on page 31 of [30] indicates that the collocation method yields a special RK method. If the collocation points $c_{1},\cdots,c_{s}$ are chosen as the Gaussian quadrature nodes, i.e., the zeros of the $s$-th shifted Legendre polynomial \frac{d^{s}}{dx^{s}}\Big{(}x^{s}(x-1)^{s}\Big{)}, Scheme 3.4 is called the Gaussian collocation method. Based on the Gaussian quadrature nodes, the interpolating quadrature formula has order $2s$, and the Gaussian collocation method shares the same order $2s$. Collocation points for Gaussian collocation methods of order 4 and 6 are given explicitly in [30].

For conservative systems with quadratic invariants, the Gaussian collocation methods have been proven to conserve the corresponding discrete quadratic invariants. For more details, please refer to the book [30]. Applying the theory to our reformulated EQ system (2.8), we have the following theorem.

Theorem 3.3**.**

The $s$-stage HEQ Gaussian collocation Scheme 3.4 is unconditionally energy stable, i.e., it satisfies the following energy dissipation law

[TABLE]

*where $F^{n}=\frac{1}{2}\|\Psi^{n}\|_{\mathcal{B}}^{2}-A,$ $c_{i}$ ($i=1,\cdots,s$) is the Gaussian quadrature nodes, $b_{i}\geq 0$ ($i=1,\cdots,s$) are the Gauss-Legendre quadrature weights, $u(t)$ be the collocation polynomial of the Gaussian collocation methods. *

Proof.

Denoting $\Psi^{n}=u(t_{n})$ and $\Psi^{n+1}=u(t_{n+1}),$ we have

[TABLE]

where the self-adjoint property of $\mathcal{B}$ was used. The integrand \Big{(}\dot{u}(t),\mathcal{B}u(t)\Big{)} is a polynomial of degree $2s-1$, which is integrated without error by the $s$-stage Gaussian quadrature formula. It therefore follows from the collocation condition that

[TABLE]

which leads to (3.17). This completes the proof.

For general RK methods, we have the following theorem.

Theorem 3.4**.**

If the coefficients of a HEQ-RK method satisfy

[TABLE]

then it is unconditionally energy stable, i.e., it satisfies the following energy dissipation law

[TABLE]

*where $F^{n}=\frac{1}{2}\|\Psi^{n}\|_{\mathcal{B}}^{2}-A,$ $\Psi_{i}=\Psi^{n}+\Delta t\sum\limits_{j=1}^{s}a_{ij}k_{j}.$ *

Proof.

Denoting $\Psi^{n+1}=\Psi^{n}+\Delta t\sum\limits_{i=1}^{s}b_{i}k_{i}$ and noticing that operator $\mathcal{B}$ is linear and self-adjoint, we have

[TABLE]

which implies

[TABLE]

Applying $\Psi^{n}=\Psi_{i}-\Delta t\sum\limits_{j=1}^{s}a_{ij}k_{j}$ to (3.23), we obtain

[TABLE]

Combining (3.20), (3.24) and $k_{i}=\mathcal{N}(\Psi_{i})\mathcal{B}\Psi_{i}$, we obtain (3.21). This completes the proof.

Remark 3.5**.**

Due to the collocation method reduces to a special RK method, we have to solve the nonlinear system (3.15), which will be implemented by using the following simple fixed-point iteration method. Denote $\mathcal{N}(\Psi)=\mathcal{N}_{1}+\mathcal{N}_{2}(\Psi)$, where $\mathcal{N}_{1}$ and $\mathcal{N}_{2}(\Psi)$ are the linear part and nonlinear part of $\mathcal{N}(\Psi)$, respectively. At time step $n$, the nonlinear system for $k_{i}$ in (3.15) is first computed by

[TABLE]

where we take the initial iteration $k_{i}^{0}=0$ for simplicity. We iterate the solution until the following criteria is satisfied

[TABLE]

*We note that the FFT algorithm is applied for solving the linear equation system (3.25). Then we obtain $\Psi^{n+1}$ by (3.15). *

Remark 3.6**.**

*Though the energy (2.3) and (2.4) are equivalent in the continuum form, the proposed schemes only satisfy a discrete energy dissipation law in term of the reformulated energy (2.4), instead of the original energy (2.3). However, we point out the discrete version of (2.4) is a high-order approximation of (2.4), i.e. (2.3). In addition, the stabilization technique is introduced to inherit the stability of the original energy as much as possible. For instance, in the Ginzburg-Landau free energy (2.5), the modified energy is $F=\frac{\varepsilon^{2}}{2}\|\nabla\phi\|^{2}+\frac{\gamma_{0}}{2}\|\phi\|^{2}+\frac{1}{2}\|q\|^{2}-A$, which implies $H^{1}$ boundedness of phase variable $\phi$. More details, please refer to [3, 41]. *

Remark 3.7**.**

*We note that the EQ approach is employed for developing second-order linear or arbitrarily high-order schemes in this paper, for respecting the energy-dissipation law. For the gradient flow models equipped with singular energy potentials, such as the Cahn-Hilliard equation with the Flory-Huggins free energy, we can’t prove that the proposed methods preserve the positivity. For some seminal work of the positive-preserving algorithms for gradient flow models, please refer to [5]. *

Remark 3.8**.**

*For gradient flow models, a class of high-order convex splitting Runge-Kutta schemes have been developed recently [46]. However, these schemes only work when one can find a convex splitting of the free energy, which may not always be true. Besides, they require too many multi-stages to reach even the 3rd order accuracy. The existence of higher order convex splitting RK coefficients is not guaranteed. Moreover, they don’t work for the case where the mobility is variable. The HEQ schemes introduced in this paper do not have these constraints so that they can be applied to a much broader class of problems. *

4 Spatial discretization

Next, we present structure-preserving spacial discretizations for the quadratized gradient flow models. For the spatial discretization, one idea inspired by [9] is to preserve the negative semi-definite property of operator $\mathcal{N}(\Psi)$ and the self-adjoint, positive definite property of operator $\mathcal{B}$. Another idea is to develop spatial discretization methods that preserve the discrete integration-by-parts formulae (please see [21, 23, 24, 25, 22, 60] for details). Based on these ideas, we apply the Fourier pseudospectral method in space to (2.8), which leads to an ODE system that preserves the spatial semi-discrete energy dissipation law. Then we apply the methods discussed in the previous section to the ODE system to obtain fully discrete energy stable schemes.

Let $N_{x},N_{y}$ be two positive even integers. The spatial domain $\Omega=[0,L_{x}]\times[0,L_{y}]$ is uniformly partitioned with mesh size $h_{x}=L_{x}/N_{x},h_{y}=L_{y}/N_{y}$ and

[TABLE]

We define $S_{N}=\textrm{span}\{X_{j}(x)Y_{k}(y),j=0,1,\ldots,N_{x}-1;k=0,1,\ldots,N_{y}-1\}$ as the interpolation space, where $X_{j}(x)$ and $Y_{k}(y)$ are trigonometric polynomials of degree $N_{x}/2$ and $N_{y}/2$, given respectively by

[TABLE]

where $a_{m}=\begin{cases}1,|m|<N_{x}/2,\\ 2,|m|=N_{x}/2,\end{cases}\mu_{x}=2\pi/L_{x},~{}b_{m}=\begin{cases}1,|m|<N_{y}/2,\\ 2,|m|=N_{y}/2,\end{cases}\mu_{y}=2\pi/L_{y}.$ We define the interpolation operator $I_{N}:C(\Omega)\rightarrow S_{N}$ as follows:

[TABLE]

where $u_{j,k}=u(x_{j},y_{k})$. The key of spatial Fourier pseudospectral discretization is to obtain derivative $\partial_{x}^{s_{1}}\partial_{y}^{s_{2}}I_{N}u(x,y)$ at collocation points. Then, we differentiate (4.2) and evaluate the resulting expressions at point $(x_{j},y_{k})$ as follows

[TABLE]

where $\mathbf{D}_{s_{1}}^{x}$ and $\mathbf{D}_{s_{2}}^{y}$ are $N_{x}\times N_{x}$ and $N_{y}\times N_{y}$ matrices, respectively, with elements given by

[TABLE]

Here we note that the Fourier pseudospectral method preserves discrete integration-by-parts formulae. For more details, please refer to our previous work [23].

Applying the Fourier pseudospectral method to (2.8), we obtain

[TABLE]

where $\mathcal{N}_{d}(\Psi)$ is a discrete negative semi-definite operate that approximates to $\mathcal{N}(\Psi)$, $\mathcal{B}_{d}$ is a constant self-adjoint, positive definite operate that approximates to $\mathcal{B}.$ It is readily to show that the system (4.3) possesses the discrete energy dissipation law

[TABLE]

where $F_{h}=\frac{1}{2}\left(\Psi,\mathcal{B}_{d}\Psi\right)_{h}-A,$ $(\cdot,\cdot)_{h}$ is the corresponding discrete inner product. Then the linear prediction-correction schemes and Gaussian collocation methods proposed in Section 3 can be applied directly for (4.3) to obtain fully discrete energy stable schemes.

Next, we specifically apply our spatial discretization method to two examples. Firstly, we consider the following Cahn-Hilliard equation

[TABLE]

where the free energy functional $F$ is given by

[TABLE]

By introducing $q=\frac{1}{\sqrt{2}}(\phi^{2}-a-\gamma_{0}),$ the free energy can be rewritten as

[TABLE]

where $\mathcal{L}=\gamma_{0}-b\Delta+c(1+\Delta)^{2}$ and $A=\frac{1}{4}(a+\gamma_{0})^{2}|\Omega|.$ Assume the boundary conditions annihilate the boundary terms. Then Eq. (4.5) reduces to the following system

[TABLE]

which can be written into the compact form (2.8) with $\Psi=(\phi,q)^{T},$ $\mathcal{N}(\Psi)=(1,\sqrt{2}\phi)^{T}\lambda\Delta(1,\sqrt{2}\phi)$ and $\mathcal{B}=\textrm{diag}(\mathcal{L},1).$ Following the notations of Ref. [23] and applying the Fourier pseudospectral method to (4.8), we obtain

[TABLE]

where $\Delta_{h}=\mathbf{D}_{2}^{x}ⓧ+\mathbf{D}_{2}^{y}ⓨ$ and $\mathcal{L}_{d}=\gamma_{0}-b\Delta_{h}+c(1+\Delta_{h})^{2}.$ System (4.9) can be written into the form of (4.3) with discrete operators $\mathcal{N}_{d}(\Psi)=(1,\sqrt{2}\phi\odot)^{T}\lambda\Delta_{h}(1,\sqrt{2}\phi\odot)$ and $\mathcal{B}_{d}=\textrm{diag}(\mathcal{L}_{d},1).$

Next we consider the molecular beam epitaxy model

[TABLE]

where the free energy functional $F$ is given by

[TABLE]

Letting q=\frac{1}{\sqrt{2}}\big{(}|\nabla\phi|^{2}-1-\gamma_{0}\big{)}, we obtain the corresponding free energy

[TABLE]

where $\mathcal{L}=\varepsilon^{2}\Delta^{2}-\gamma_{0}\Delta$ and $A=\frac{1}{4}(2\gamma_{0}+\gamma_{0}^{2})|\Omega|.$ Assume the boundary conditions annihilates the boundary term. Then Eq. (4.10) is written equivalently as

[TABLE]

which can be written into the compact form of (2.8) with $\Psi=(\phi,q)^{T},$ \mathcal{N}(\Psi)=-\lambda\left(\begin{array}[]{c}1\\ \sqrt{2}\nabla\phi\cdot\nabla\end{array}\right)(1,-\nabla\cdot\sqrt{2}\nabla\phi) and $\mathcal{B}=\textrm{diag}(\mathcal{L},1).$ Applying the Fourier pseudospectral method for (4.13), we obtain

[TABLE]

where $\mathcal{L}_{d}=\varepsilon^{2}\Delta_{h}^{2}-\gamma_{0}\Delta_{h}.$ System (4.14) can be written into the form of (4.3) with discrete operators \mathcal{N}_{d}(\Psi)=-\lambda\left(\begin{array}[]{c}1\\ \sqrt{2}(\mathbf{D}_{1}^{x}ⓧ\phi)\odot\mathbf{D}_{1}^{x}ⓧ+\sqrt{2}(\mathbf{D}_{1}^{y}ⓨ\phi)\odot\mathbf{D}_{1}^{y}ⓨ\end{array}\right)\Big{(}1,-\sqrt{2}\mathbf{D}_{1}^{x}ⓧ\big{(}\mathbf{D}_{1}^{x}ⓧ\phi\big{)}\odot-\sqrt{2}\mathbf{D}_{1}^{y}ⓨ\big{(}\mathbf{D}_{1}^{y}ⓨ\phi\big{)}\odot\Big{)} and $\mathcal{B}_{d}=\textrm{diag}(\mathcal{L}_{d},1).$

5 Numerical Results

In this section, we apply the proposed numerical schemes to several gradient flow problems. For simplicity, we assume periodic boundary conditions in the numerical experiments. Note that our schemes can be easily applied to gradient flow models with physical boundary conditions as long as they annihilate the boundary terms in the energy dissipation rate formula. Systematical comparisons among the prediction-correction scheme 3.1, the BDF-k scheme 3.2, and HEQ schemes 3.3 and 3.4 are presented below.

We fix the number of iteration $N=5$ and $\varepsilon_{0}=10^{-12}$ in the rest of this paper. In addition, for easily referring the schemes, we name (3.1) with predictor (3.3) as the LCN (linear CN) scheme; (3.1) with predictor (3.6) as the ICN (implicit CN) scheme; (3.2) with predictor (3.3) as the LBDF2 scheme; (3.2) with predictor (3.7) as the IBDF2 scheme.

Example 1: Cahn-Hilliard equation. In the first example, we consider the Cahn-Hilliard equation with the free energy functional containing a double-well bulk term:

[TABLE]

where $\varepsilon$ is a small parameter. The Cahn-Hilliard equation is given as follows

[TABLE]

where $\lambda$ is the mobility parameter. If we introduce the auxiliary variable $q=\frac{1}{\sqrt{2}}(\phi^{2}-1-\gamma_{0})$, where $\gamma_{0}>0$ is a constant, we have

[TABLE]

Then we denote $\Psi=(\phi,q)^{T}$ and rewrite the Cahn-Hilliard equation into a prototypical form

[TABLE]

In the numerical experiment, we choose a square domain with $\Omega=[0,1]^{2}$, and the parameters are chosen as $\gamma_{0}=1$, $\varepsilon=0.01$ and $\lambda=10^{-3}$. Using initial condition $\phi(x,y)=\sin(2\pi x)\sin(2\pi y)$, we carry out time-step refinement tests. By choosing the numerical solution calculated by the $6$th order Gaussian collocation method with time step $\Delta t=10^{-5}$ as the “exact” solution, the errors are calculated as the differences between the numerical solutions at different time steps and the “exact” solution, respectively.

The $\log$-$\log$ plots of $L^{2}$ errors against time step $\Delta t$ are depicted in Figure 5.1. From Figure 5.1(a), we observe that the HEQ-RK scheme 3.3 with the 4th and 6th order Gaussian collocation points reaches the 4th and 6th order of accuracy, respectively. From Figure 5.1(b), we observe that both the LCN scheme and the ICN scheme reach the second order accuracy, while the ICN scheme reaches it at a larger time step. The ICN scheme has a significantly smaller error than the LCN scheme by comparing the errors of the two schemes at the same time step size, with similar CPU time cost, since the ICN/IBDF2 schemes only require up to 5 iterations, where each iteration is solved using an FFT solver. As shown in Figure 5.1(c), the BDF schemes have similar results as the CN schemes do. Besides, the BDF-k schemes also reach the expected order of convergence when the time step is small enough, as shown in Figure 5.1(d). Comparing between Figure 5.1(a) and 5.1(d), we observe the HEQ schemes are more accurate than the BDF-k schemes when using the same time steps.

Next, we study the coarsening dynamics, and compare the accuracy of the proposed schemes by comparing the numerically calculated energies using different schemes with various time-step sizes. To have a more detailed comparison, we consider coarsening dynamics of a binary mixture (which demonstrates a dramatic energy change when the system coarsens). We use $256\times 256$ equal distanced meshes in space to discretize the domain $\Omega=[0\,\,\,4\pi]\times[0\,\,\,4\pi]$. The parameters are chosen as $\lambda=0.1$, $\epsilon=0.05$, $\gamma_{0}=1$, and initial condition is \phi(x,y)=0.001\big{(}2~{}\textrm{rand}(x,y)-1\big{)}.

In Figure 5.2(a-b), the energy computed using the CN schemes in time $[0\,\,\,2]$ are plotted. We observe that for the LCN scheme, when $\Delta t<0.0125$, the calculated energy converges to the accurate energy. For ICN scheme, it has a dramatic improvement in accuracy, where the numerical energy calculated with $\Delta t=0.1$ is already very accurate, where the CPU time cost is negligible as only up to 5 times of an FFT solver is applied to (3.6) in the computation in each step. As shown in Figure 5.2(c-d), the BDF4 and BDF6 schemes provide accurate results given the time step is small enough. However, when the time step is large, it violates the energy dissipation law. The results using the HEQ Gaussian collocation method of order four and order six are summarized in Figure 5.2(e-f). We observe that the HEQ Gaussian collocation scheme provides much better accuracy. Both 4th order and the 6th order schemes with a large time step size, such as $\Delta t=0.2$, give very accurate energy predictions, which is much larger than the proper time steps for the BDF4 and BDF6 schemes respectively.

One typical simulation using the 6th order Gaussian collocation method with a time step size $\Delta t=0.01$ is depicted in Figure 5.3, where we observe fairly accurate prediction of the coarsening dynamics in various times.

Example 2: Allen-Cahn Equation. Here we test the proposed numerical schemes on the Allen-Cahn equation

[TABLE]

Introducing the same auxiliary variable, $q=\frac{1}{\sqrt{2}}(\phi^{2}-1-\gamma_{0}),$ we obtain the reformulated EQ system

[TABLE]

which can be written into

[TABLE]

Here we test the numerical schemes via a benchmark example used in [4]. We choose the parameter values as $\lambda=1,\varepsilon=1$, $L_{x}=L_{y}=256$. We use $256\times 256$ mesh points to reduce the spatial error. The initial condition is chosen as a disk,

[TABLE]

It has been shown that the radius $R$ of the disk at time $t$ is given as $R=\sqrt{R_{0}^{2}-2t}$ [4], where $R_{0}$ is the initial radius. In other words, the volume of the disk is given as $V=\pi R_{0}^{2}-2\pi t$. To test it, we implement and compare the proposed schemes with different time steps. The calculated results using the second-order CN schemes are summarized in Figure 5.4(a-b), where we observe that when $\Delta t=0.0125$, the LCN scheme can predict the correct $-2\pi$ slope for the volume decreasing rate. For the ICN scheme, even with $\Delta t=0.5$, it predicts the correct volume decreasing rate. For the BDF4 and BDF6 schemes, they require approximately $\Delta t=0.025$ to predict the volume decreasing rate accurately, as shown in Figure 5.4(c-d). The calculated volumes using the HEQ schemes with different time step sizes are summarized in Figure 5.4(e-f). It shows the Gaussian collocation method is more accurate such that even with time step size $\Delta t=5$, it predicts a very accurate volume decaying rate $-2\pi$.

In particular, the temporal evolution of the disk using ICN scheme with time step size $\Delta t=0.1$ is plotted in Figure 5.5.

Example 3: Phase Field Crystal Growth Model. Next we solve the phase field crystal growth model by the proposed schemes. The phase field crystal growth model was introduced in [10, 11, 47] in the form of a Cahn-Hilliard equation (4.5) with the free energy

[TABLE]

We set the parameter values $a=0.325,c=1.$ We use the initial conditions and parameter values given in [47], i.e. $\phi_{0}(\mathbf{x})=\overline{\phi}+\omega(\mathbf{x})(A\phi_{s}(\mathbf{x}))$, $\phi_{s}(\mathbf{x})=\cos(\frac{k}{\sqrt{3}}y)\cos(kx)-\frac{1}{2}\cos(\frac{2k}{\sqrt{3}}y)$, where $k$ represents a wavelength related to the lattice constant, $A$ represents an amplitude of the fluctuations in density, and the scaling function $\omega(\mathbf{x})$ is defined as

[TABLE]

$d_{0}$ is a prescribed parameter. In the simulation, we choose $\Omega=[0,\,\,150]^{2}$, $\varepsilon=0.325$, $\overline{\phi}=\frac{\sqrt{\varepsilon}}{2}$, A=\frac{4}{5}\Big{(}\overline{\phi}+\frac{\sqrt{15\varepsilon-36\overline{\phi}^{2}}}{3}\Big{)}, $d_{0}=25$, $k=\frac{\sqrt{3}}{2}$ and $512\times 512$ meshes. A numerical simulation with the 4th order Gaussian collocation method is summarized in Figure 5.6, where we observe qualitatively similar predictions as reported in [50, 12, 53].

Example 4: Molecular Beam Epitaxy Model. Next, we study the molecular beam epitaxy (MBE) model (4.10). With $\phi$ representing the scaled height function of the thin film, the continuum MBE model reads as

[TABLE]

By introducing q=\frac{1}{\sqrt{2}}\Big{(}|\nabla\phi|^{2}-1-\gamma_{0}\Big{)}, the schemes proposed in previous sections could be readily applied to this model. Here we omit the details for simplicity.

To compare the schemes, we carry out a standard benchmark problem used in [34, 33]. Given the domain $[0,2\pi]^{2}$ and parameter values $\lambda=1$, $\epsilon^{2}=0.1$, and initial profile $\phi(x,y)=0.1(\sin(2x)\sin(2y)+\sin(5x)\sin(5y))$. We use $256\times 256$ meshes. The comparisons are summarized in Figure 5.7. It shows that even though the LCN scheme is unconditionally energy stable, it predicts erratic coarsening dynamics with large time steps. Also, the BDF-6 scheme requires relatively small time step to predict the accurate energy curve. On the contrary, the ICN scheme and the HEQ Gaussian collocation scheme can predict robust dynamics with larger time steps. As we have alluded to earlier that the ICN scheme only increases the CPU time slightly for up to 5 times of an FFT solve in each predictor-corrector step (3.6).

6 Conclusion

In this paper, we have demonstrated two general strategies to derive unconditionally energy-stable numerical approximations for thermodynamically consistent gradient flow models. In the first strategy, we present a prediction-correction approach to derive effective, linear schemes to solve the nonlinear thermodynamically consistent gradient flow models and in the meantime improve the accuracy of the schemes. The new schemes retain second-order accurate in time, and their errors are much smaller than those of the linear schemes. Moreover, the 2nd order convergence rate is attained even with a large time step. In the second strategy, we propose a novel idea utilizing the quadratic-invariant preserving Runge-Kutta multi-stage discretization, resulting in arbitrary order, unconditionally energy stable numerical approximations for a general class of gradient flow models, which we named HEQ schemes. Unconditional energy stability is established rigorously.

Several numerical experiments for gradient flow problems are presented to illustrate the accuracy and efficiency of the numerical schemes. The gradient flow models include the Cahn-Hilliard equation with the Ginzburg-Landau double well free energy, the Allen-Cahn equation with the same energy functional, the crystal growth model, and the molecular beam epitaxy grwoth model. These models are solved with the proposed schemes, and numerical comparisons are presented. Through numerical experiments, it is clear that the prediction-correction CN scheme yields a significantly smaller error than the linear CN scheme, while the prediction-correction BDF2 scheme demonstrates the similar improvement over the lower order method at large time steps. This newly proposed prediction-correction schemes and the high-order schemes based on Runge-Kutta methods are rather general that they can be readily applied to a large class of thermodynamically consistent models. Also, instead of using the EQ approach, the proposed numerical strategies here could easily be applied to the reformulated models using the SAV approach to derive high-order, efficient energy stable schemes.

Acknowledgment

Yuezheng Gong’s work is partially supported by the Natural Science Foundation of Jiangsu Province (Grant No. BK20180413) and the National Natural Science Foundation of China (Grant No. 11801269). Qi Wang’s research is partially supported by NSF-DMS-1517347, DMS-1815921, and OIA-1655740 award, NSFC awards #11571032, #91630207 and NSAF-U1530401. Jia Zhao’s work is partially supported by National Science Foundation under grant number NSF DMS-1816783.