1 Introduction↩︎

Many physical problems such as interfacial dynamics [1], [2], crystal growth [3], [4], image inpainting [5], [6], thin films [7], [8], and polymers [9] can be modeled by gradient flow models which comply with the second law of thermodynamics. As a special kind of nonlinear partial differential equations, the gradient flow model can be derived from the mobility and the variational derivative of free energy when the total free energy is known. Generally, consider the spatial-temporal domain \(\Omega_t:=\Omega \times \left(0,T\right]\), \(\Omega \subset \mathbb{R}^d, d=1,2,3\) and the state variable \(\phi =\phi \left(\boldsymbol{x} , t\right) :\Omega _t\mapsto \mathbb{R}\), the total free energy \(E(\phi)\) consists of a quadratic term \(E_{\mathcal{L}}(\phi)\) and a strongly nonlinear term \(E_{\mathcal{N}}(\phi)\), namely \[\begin{align} \label{E0} E\left( \phi \right) =E_{\mathcal{L}}\left( \phi \right) +E_{\mathcal{N}}\left( \phi \right) =\frac{1}{2}\int_{\Omega}{\phi \mathcal{L}\phi}\, \text{d}\boldsymbol{x}+\int_{\Omega}{F\left( \phi \right)}\,\text{d}\boldsymbol{x}, \end{align}\tag{1}\] where \(F(\phi)\) is a nonlinear potential function, and \(\mathcal{L}\) is a linear self-adjoint non-negative definite operator. Then the corresponding general gradient flow system of 1 reads \[\label{PF} \begin{align} {}&\frac{\partial \phi}{\partial t}=-M\mathcal{G}\mu ,\\ {}&\mu =\mathcal{L}\phi +f\left( \phi \right) , \end{align}\tag{2}\] with the initial condition \(\phi(\boldsymbol{x},0)\) and thermodynamically consistent boundary conditions (such as periodic and homogeneous Neumann boundaries), where \(f(\phi)=F'(\phi)\), \(\mu =\frac{\delta E}{\delta \phi}\) is the chemical potential, \(M\) is a mobility constant and \(\mathcal{G}\) is a non-negative definite operator reflecting the dissipation path. The gradient flow system \(\eqref{PF}\) preserves the energy dissipation law \[\begin{align} \label{dE} \frac{dE\left( \phi \right)}{dt}=\left( \frac{\delta E}{\delta \phi},\frac{\partial \phi}{\partial t} \right) =\left( \mathcal{L}\phi +f\left( \phi \right) ,-M\mathcal{G}\mu \right) =-\left( \mu ,M\mathcal{G}\mu \right) \le 0, \end{align}\tag{3}\] where \(\left( \cdot ,\cdot \right)\) denotes the \(L^2\) product, i.e. \(\left( f,g \right) =\int_{\Omega}{fg}\,\text{d}\boldsymbol{x}, \forall f,g\in L^2\left( \Omega \right)\).

In order to eliminate non-physical numerical solutions, the energy dissipation law must be preserved at the discrete level when designing numerical schemes for nonlinear gradient flow models. In recent years, there has been tremendous interest in designing efficient and accurate energy stable schemes for the nonlinear dissipation systems. To this end, many efficient schemes have been developed, including convex splitting method [10], [11], the exponential time differencing (ETD) [12], [13] method, stabilized semi-implicit (SSI) approach [14], [15], and the new Lagrange multiplier method [16]. Recently, a numerical approach named invariant energy quadratization (IEQ) or energy quadratization (EQ) method is proposed to solve thermodynamically consistent PDE problems, and by utilizing auxiliary variables, the original PDE problems are reformulated into equivalent PDE problems. Based on the baseline IEQ method, linear, second-order even any high-order, easy-to-implement, and unconditionally energy stable schemes are constructed for complex gradient flow models [17]–[25]. Subsequently, the scalar auxiliary variable (SAV)-type approaches [26]–[29] are proposed and can also be applied to solve a class of gradient flow models.

However, for the baseline IEQ method, the truncation errors are introduced in the long-time numerical simulation process, so that the numerical solutions of auxiliary variables are no longer equivalent to their original continuous definitions. We know that the baseline IEQ method preserves the modified energy dissipation law with respect to the auxiliary variables instead of the original variables, which means the original energy corresponding to the numerical solutions is not necessarily monotonically decreasing. This inconsistency also leads to a loss of accuracy when the time step is not small enough. In order to overcome the adverse effects of this inconsistency, the pioneering relaxation techniques [29]–[32] are used to correct the numerical solutions of the auxiliary variables, so that the modified energy is closer to the original energy. In this paper, inspired by a novel and essential energy-optimal technique based on the SAV method, i.e. EOP-SAV method [33], we propose an energy-optimized approach based on the baseline IEQ method, and construct first- and second-order (in time) unconditionally energy stable schemes which can solve a large class of gradient flow problems. We effectively penalize the inconsistency between the numerical solution of auxiliary variables and their continuous definitions by using the proposed energy-optimized IEQ (EOP-IEQ) method and the main advantages of our approach are listed as follows:

• The EOP-IEQ method inherits all the advantages of the baseline IEQ method [17] and relaxed IEQ (REQ) method [30];

• Our energy-optimized technique makes the calculation process of updating auxiliary variables more simpler and efficient;

• Compared with the IEQ and REQ methods, the modified energy obtained by the energy-optimized technique is closer to the original energy.

All numerical schemes based on the EOP-IEQ method are linear unconditionally energy stable and easy-to-implement. In numerical examples of various gradient flow models, the modified energy corresponding of the numerical solutions based on the auxiliary variables almost satisfy the original energy dissipation law.

The rest of this paper is organized as follows. In Section 2, we briefly review the baseline IEQ and REQ methods for gradient flow models. In Section 3, we propose the EOP-IEQ method for general dissipative systems, and prove that the new first- and second-order (in time) numerical schemes are unconditionally energy stable. In Section 4, we verify its optimality by comparing the EOP-IEQ method with the baseline IEQ and REQ methods, and provide a large number of numerical experiments to verify the accuracy, efficiency and energy stability of the proposed schemes.

2 A brief review of the IEQ-type methods↩︎

In this section, we review the baseline IEQ method which considered by Yang in [17] and the relaxed IEQ (REQ) approach proposed in [30] by Zhao.

2.1 The baseline IEQ method↩︎

Firstly, we need to assume that the bulk free energy \(F(\phi)\) is bounded from below which means \(F(\phi) > - C_*\) for a positive constant \(C_*\). Introduce an auxiliary variable \[\begin{align} \label{eq95q} q\left( \boldsymbol{x},t \right) := Q\left( \phi\left( \boldsymbol{x},t \right) \right) =\sqrt{F\left( \phi\left( \boldsymbol{x},t \right) \right) +C_0},\quad C_0>C_*, \end{align}\tag{4}\] where the positive constant \(C_0\) makes sure \(q(\boldsymbol{x},t)\) a well defined real function. Then we rewrite the gradient flow system 2 as the following equivalent system \[\label{IEQ} \begin{align} {}&\frac{\partial \phi}{\partial t}=-M\mathcal{G}\mu ,\\ {}&\mu =\mathcal{L}\phi +\frac{q}{Q(\phi)}f\left( \phi \right) ,\\ {}&\frac{\partial q}{\partial t}=\frac{f\left( \phi \right)}{2Q(\phi)}\frac{\partial \phi}{\partial t}. \end{align}\tag{5}\] Using the IEQ method above, we convert the free energy 1 into the following quadratic form \[\begin{align} \label{E_{1}^{n+1}} \bar{E}( \phi ,q ) =\frac{1}{2}( \mathcal{L}\phi ,\phi ) +(q,q)-C_0|\Omega|, \end{align}\tag{6}\] It is not difficult to obtain the following modified energy dissipation law \[\begin{align} \label{dE_{1}^{n+1}} \frac{d\bar{E}\left( \phi,q \right)}{dt}=-M\left( \mathcal{G}\mu ,\mu \right) \le 0. \end{align}\tag{7}\]

Now we introduce some notations that will be used throughout the paper. Let \(L^p(\Omega)\) denote the usual Lebesgue space on \(\Omega\) with the norm \(\|\cdot\|_{L^p}\). \(W^{k,p}(\Omega)\) stands for the standard Sobolev spaces equipped with the standard Sobolev norms \(\|\cdot\|_{k,p}\). For \(p = 2\), we write \(H^k(\Omega)\) for \(W^{k,2}(\Omega)\) and the corresponding norm is \(\|\cdot\|_k\), where \(\|\cdot\|_0\) is simply represented by \(\|\cdot\|\).

We discretize time domain \([0,T]\) into uniform meshes \(0=t_0<t_1<\cdots<t_N=T\) with \(t_n=n\tau\), \(\tau=T/N\), and use \((\bullet )^{n+1}\) to denote the numerical approximation of \((\bullet )\) at \(t_{n+1}\). With above notations, we obtain the following semi-discrete \(k\)th-order IEQ/BDF\(k\) schemes based on the backward Euler formulas.

2.2 The REQ method↩︎

To overcome the inconsistency issue of the modified energy and original energy at the discrete level, Zhao [30] proposes a relaxation technique based on the baseline IEQ method, which called REQ method to overcome the numerical difference between \(q^{n+1}\) and \(Q(\phi^{n+1})\). The semi-implicit relaxed IEQ Crank-Nicolson (REQ/CN) scheme is as follows.

::: {#th-E_{1}^{n+1} .lemma} Lemma 2 ([30]). The REQ/CN scheme 2 is unconditionally energy stable in the sense that \[\frac{1}{2}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) +\lVert q^{n+1} \rVert ^2-\frac{1}{2}\left( \mathcal{L}\phi ^n,\phi ^n \right) -\lVert q^n \rVert ^2=-\tau M \left( \mathcal{G}\mu ^{n+\frac{1}{2}},\mu ^{n+\frac{1}{2}} \right) \le 0.\] :::

3 The EOP-IEQ method↩︎

The major issue is that the quxiliary variable \(q(t_{n+1})\) is no longer equal to \(Q(\phi(t_{n+1}))\) numerically. Thus the modified energy law [dE_{1}^{n+1}] does not necessarily satisfy the original energy law 3 at the discrete level. To better restrain the inconsistency issue between \(q(t_{n+1})\) and \(Q(\phi(t_{n+1}))\) (that are supposed to be equal as introduced in 4 ), Liu et.al. propose a novel modified SAV method [33], named energy-optimized SAV (EOP-SAV) approach which is unconditionally energy stable with respect to a modified energy that is closer to the original energy than the baseline SAV[26] and relaxed SAV [31] approaches. Inspired by this method, the EOP-IEQ method is considered in this paper, which has the optimal correction on modified energy compared to the IEQ and REQ methods. The main idea of the energy-optimized (EOP) technique is to correct the auxiliary variable \(q^n\) from the perspective of original nonlinear energy \(E_{\mathcal{N}}\left( \phi^n \right)\), so that the modified energy is closest to the original energy at the discrete level, and the numerical schemes preserve the energy dissipation law.

3.1 The EOP-IEQ/CN scheme.↩︎

Now we consider the following second-order Crank-Nicolson scheme based on the novel EOP-IEQ approach.

::: {#th-E_{2}^{n+1} .theorem} Theorem 1. The energy-optimized technique of \(q^{n+1}\) ?? in \(\boldsymbol{Step~2}\) is an optimal choice to correct the modified energy dissipative law of the baseline IEQ method [17]. :::

Now we establish the consistency estimate for the EOP-IEQ/CN scheme 3. The consistency errors \(\varepsilon _{1}^{n+1}\), \(\varepsilon _{2}^{n+1}\), \(\varepsilon _{3}^{n+1}\) for the EOP-IEQ/CN method are determined by the following: \[\label{trun} \begin{align} {}&\frac{\phi _{\star}^{n+1}-\phi _{\star}^{n}}{\tau}=-M\mathcal{G}\left( \mathcal{L}\phi _{\star}^{n+\frac{1}{2}}+H\left( \phi _{\star}^{n+\frac{1}{2}} \right) \hat{q}_{\star}^{n+\frac{1}{2}} \right) +\varepsilon _{1}^{n+1},\\ {}&\frac{\hat{q}_{\star}^{n+1}-q_{\star}^{n}}{\tau}=\frac{H\left( \phi _{\star}^{n+\frac{1}{2}} \right)}{2}\left( \frac{\phi _{\star}^{n+1}-\phi _{\star}^{n}}{\tau} \right) +\varepsilon _{2}^{n+1},\\ {}&q_{\star}^{n+1}=\lambda ^{n+1}Q\left( \phi _{\star}^{n+1} \right)+\varepsilon _{3}^{n+1} ,\quad \lambda ^{n+1}=\min \left\{ 1,\sqrt{\frac{E^{n+1}_2}{E^{n+1}_1}} \right\}, \end{align}\tag{13}\] where \(\phi _{\star}^{n+1}=\phi \left( t_{n+1} \right)\), \(q_{\star}^{n+1}=q\left( t_{n+1} \right)\), \(Q_{\star}^{n+1}=Q\left( \phi(t_{n+1}) \right)\).

Suppose that the analytical solutions (\(\phi,q\)) of the equivalent gradient flow model 5 possess the following regularity condition \[\label{zheng} \left\{ \begin{array}{l} \phi \in L^{\infty}\left( 0,T;H^2\left( \Omega \right) \right) ,\quad q\in L^{\infty}\left( 0,T;L^{\infty}\left( \Omega \right) \right) ,\\ \phi _t\in L^{\infty}\left( 0,T;L^{\infty}\left( \Omega \right) \right) ,\quad q_{tt},\phi _{tt}\in L^2\left( 0,T;L^2\left( \Omega \right) \right). \end{array} \right.\tag{14}\] We first give the following lemma to establish the quantitative relation between the \(L^2\) norm of \(H( \phi _{\star}^{n+1}) -H( \phi _{\star}^{n} )\) and \(\phi _{\star}^{n+1}-\phi _{\star}^{n}\) under the regularity condition 14 .

3.2 The EOP-IEQ/BDF\(k\) schemes↩︎

If we use the backward Euler formulas for the time discretization, we will have the semi-implicit EOP-IEQ/BDF\(k\) \((k=1,2)\) schemes as below.

3.3 energy-optimized techniques for the variations of the baseline IEQ method↩︎

Due to the wide applicability of the baseline IEQ method, several variants have been proposed to be suitable for different situations [34], [35]. Since our proposed energy-optimized technique does not restrict the form of auxiliary variables, our EOP-IEQ approach can be easily extended to gradient flow models that contain multiple auxiliary variables. Below we describe the EOP-IEQ method applied to a gradient flow with one phase field variable \(\phi\). We consider a more general form of the free energy \[\begin{align} \label{ME} \mathcal{E}\left( \phi \right) =\frac{1}{2}\left( \mathcal{L}\phi ,\phi \right) +\sum_{i=1}^k{\left( F_i\left( \phi \right) ,1 \right)}, \end{align}\tag{18}\] where the operator \(\mathcal{L}\) is a linear self-adjoint non-negative definite operator and \(F_i\left( \phi \right)\), \(i=1\), \(2\), \(\cdots\), \(k\) are the bulk potentials. The general gradient flow system of the free energy 18 as follows \[\label{MPF} \begin{align} {}&\frac{\partial \phi}{\partial t}=-M\mathcal{G}\mu ,\\ {}&\mu =\mathcal{L}\phi +\sum_{i=1}^k{F_{i}^{'}\left( \phi \right)}, \end{align}\tag{19}\] which has the following energy dissipation law \[\frac{dE\left( \phi \right)}{dt}=-\mathcal{G}\left( \mu ,\mu \right)\le 0 ,\] where \(M\) is a mobility constant and the operator \(\mathcal{G}\) is a non-negative definite operator. For the system 19 , assume that each nonlinear bulk potentials have a lower bound \(F_i(\phi)\ge C^*_i\), \(i=1\), \(2\), \(\cdots\), \(k\). Denote \(\boldsymbol{q}=\left( q_1,q_2,\dots ,q_k \right)\), we introduce multiple auxiliary variables \[q_i( \boldsymbol{x},t ) =:Q_i(\phi( \boldsymbol{x},t )) =\sqrt{F_i( \phi( \boldsymbol{x},t ) ) +C^0_i},\quad C^0_i>C^*_i,i=1,2,\dots,k.\] Then the equivalent form of the gradient flow system 19 can be written as \[\label{MIEQ} \begin{align} {}&\frac{\partial \phi}{\partial t}=-M\mathcal{G}\mu ,\\ {}&\mu =\mathcal{L}\phi +\sum_{i=1}^k{q_iH_i( \phi )},,\\ {}&\frac{\partial q_j}{\partial t}=\frac{H_j( \phi )}{2}\frac{\partial \phi}{\partial t},\quad j=1,2,\dots ,k. \end{align}\tag{20}\] where \(H_i( \phi ) =\frac{f_i( \phi )}{Q_i( \phi )},f_i=F_{i}^{'}\left( \phi \right),i=1,2,\dots ,k.\) Then we can get the modified free energy as \[\begin{align} \bar{\mathcal{E}}\left( \phi,\boldsymbol{q} \right) =\frac{1}{2}\left( \mathcal{L}\phi ,\phi \right) +\sum_{i=1}^k{\lVert q_i \rVert ^2}-\sum_{i=1}^k{C^0_i|\Omega|}. \end{align}\] For the reformulated system 20 , it has the following modified energy dissipation law \[\frac{d\bar{E}( \phi ,\boldsymbol{q} )}{dt}=-M\mathcal{G}( \mu ,\mu ) \le 0.\] We introduce the energy-optimized technique into the system 20 to improve the inconsistency between the modified energy and the original energy after discretization. Then a second-order EOP-IEQ scheme based on Crank-Nicolson formula can be written as follows.

The proof of this theorem is not particularly difficult but will not be reproduced here.

4 Numerical implementation and results↩︎

In this section, we adopt some widely used gradient flow models including Allen-Cahn (AC) equation, Cahn-Hilliard (CH) equation, molecular beam epitaxy (MBE) model, and phase field crystal (PFC) model to verify the generality of EOP-IEQ method. Furthermore, the corresponding numerical algorithms are implemented to verify the above theoretical results. For simplicity, we consider periodic boundary condition and use the Fourier spectrum method in space in all examples, where \(N_x\) and \(N_y\) represent the number of Fourier modes along the \(x\) and \(y\) axis respectively. It is worth noting that the proposed schemes can also be applied to other thermodynamically consistent boundary conditions which satisfy energy dissipation law.

4.1 Allen-Cahn equation↩︎

To begin with, we consider the free energy [36] \[\begin{align} \label{E-ACH} E\left( \phi \right) =\int_{\Omega}{\frac{\alpha_0}{2}\left| \nabla \phi \right|^2}+\frac{\alpha_0}{\varepsilon ^2}F\left( \phi \right) \text{d}\boldsymbol{x,}\quad F\left( \phi \right) =\frac{\left( \phi ^2-1 \right) ^2}{4}, \end{align}\tag{21}\] where the interface parameters \(0<\) \(\alpha_0\), \(\varepsilon\) \(\ll 1\).

The Allen-Cahn (AC) equation [37] can be expressed in the energy variational form 2 of the free energy 21 with \(\mathcal{G} = MI\) as follows \[\begin{align} \label{AC} \frac{\partial \phi}{\partial t}=-M\left( -\alpha_0 \Delta \phi +\frac{\alpha_0}{\varepsilon ^2}f\left( \phi \right) \right). \end{align}\tag{22}\] where \(f\left( \phi \right) =F'(\phi)=\phi ^3-\phi\). Now we introduce the auxiliary variable \(q:=Q\left( \phi \right) =\sqrt{F\left( \phi \right) +C_0}\), then the general model 5 is specified as \[\label{AC-IEQ} \begin{align} {}&\frac{\partial \phi}{\partial t}=-M\left( -\alpha_0 \Delta \phi +\frac{\alpha_0 }{\varepsilon ^2}\frac{q}{Q\left( \phi \right)} f\left( \phi \right)\right) ,\\ {}&\frac{\partial q}{\partial t}=\frac{f\left( \phi \right)}{2Q\left( \phi \right)}\frac{\partial \phi}{\partial t}. \end{align}\tag{23}\]

Case A. We check the convergence rate in time of the proposed schemes for the AC equation under the \(2D\) computational region \(\Omega = [-0.5,0.5]^2\). The initial data is chosen to be smooth \[\begin{align} \label{AC-init} \phi \left( x,y,0 \right) =\tanh\frac{a+b\cos \left( 6\theta \right) -c \sqrt{x^2+y^2}}{\sqrt{2}\varepsilon}, \quad \theta =\arctan \frac{y}{x}, \end{align}\tag{24}\] where \(a,b,c\) are positive constants. We set \((N_x, N_y) =(64,64)\), so that the spatial discrete error is negligible compared to the temporal discrete error. The parameters are \((a,b,c) = (1.5,1.2,2\pi)\), \(\alpha_0=0.01\), \(\varepsilon =0.01\), \(M=0.6\), \(C_0=100\). In addition, since the exact solution is unknown, we choose the EOP-IEQ/CN scheme with time step \(\tau = 1e-6\) as the reference solution for the calculation error. In Fig. 1, we obtain the numerical errors of \(L^2\) norm and convergence rates at \(T = 0.1\) using various first- and second-order schemes. We observe that (i) the temporal convergence rates obtained by baseline IEQ and EOP-IEQ methods are consistent with theoretical expectations in all cases; (ii) the \(L^2\) errors of the numerical solutions \(\phi\) and the auxiliary variable \(q\) in the EOP-IEQ schemes are significantly smaller than those in the baseline IEQ schemes.

Case B. We also choose the initial condition 24 to compare the accuracy of the baseline IEQ, REQ and the EOP-IEQ methods for solving the AC equation. We adjust some parameters \((N_x, N_y) =(256,256)\), \(\eta=0.5\), \(M=0.51\), \(C_0=1\), \(\tau=1e-2\) and \(T=10.5\) and use the results of the EOP-IEQ/BDF2 scheme of \(\tau=1e-6\) as the reference solution. The profiles of \(\phi\) at various times are summarized in Fig. 2. A comparison of the modified energy (left), the energy error (middle), and the update parameters (right) for the IEQ/CN, REQ/CN, and EOP-IEQ/CN schemes shown in Fig. [F14], where the update parameters include relaxation parameters \(\xi^n\) and energy-optimized parameters \(\lambda^n\), where \(\lambda_1^n=|\lambda^n-1|\). We observe that the energy-optimized parameters \(\lambda_1^n=0, \forall n\ge 1\) in the EOP-IEQ/CN scheme, which means \(\lambda^n=1\) and \(q^{n}=Q(\phi^n), \forall n\ge 1\). It is essential to preserve the consistency of the modified energy with the original energy. Meanwhile, Fig. [F14] (b) shows that the errors between the corrected energy and the reference energy using the IEQ/CN and REQ/CN schemes are significantly larger. It is verified that our proposed EOP-IEQ approach ensure the modified energy closer to original energy than the REQ approach.

4.2 Cahn-Hilliard equation↩︎

Now we consider the Cahn-Hilliard (CH) equation given as \[\label{CH} \begin{align} {}&\frac{\partial \phi}{\partial t}=-M\Delta \mu ,\\ {}&\mu =-\alpha _0\Delta \phi +\frac{\alpha _0}{\varepsilon ^2}f\left( \phi \right) . \end{align}\tag{25}\] It can be expressed in the energy variational form 2 of the free energy 21 with \(\mathcal{G} = -M\Delta\). We introduce the auxiliary variable \(q:=Q\left( \phi \right) =\sqrt{F\left( \phi \right) -\kappa\phi ^2/2+C_0}\), then the general model in 5 is specified as \[\begin{align} {}&\frac{\partial \phi}{\partial t}=-M\Delta \left( -\alpha _0\Delta \phi +\frac{\alpha _0}{\varepsilon ^2}\kappa \phi +\frac{q}{Q\left( \phi \right)}\frac{\alpha _0}{\varepsilon ^2}f_{\kappa}\left( \phi \right) \right) ,\\ {}&\frac{\partial q}{\partial t}=\frac{f_{\kappa}\left( \phi \right)}{2Q\left( \phi \right)}\frac{\partial \phi}{\partial t}, \end{align}\] where \(f_{\kappa}\left( \phi \right) =\phi^3-\phi-\kappa\phi\). In addition, the CH equation satisfies the law of conservation of mass \[\begin{align} \label{MM} \frac{dm\left( t \right)}{dt}=0,\quad m\left( t \right) =\int_{\Omega}{\phi \left( \boldsymbol{x,}t \right)}\,\text{d}\boldsymbol{x}. \end{align}\tag{26}\]

Case A. The initial conditions [38] is considered as the following discrete form \[\begin{align} \label{CH95init} \phi \left( x,y,0 \right) =0.25+0.4rand\left( x,y \right) , \end{align}\tag{27}\] where \(rand(x,y)\) represents a pseudo-random value derived from a uniform distribution within \((-1,1)^2\). We choose the numerical parameters \(\Omega = [-0.5,0.5]^2, (N_x,N_y)=(64,64)\), \(\alpha_0=0.01^2\), \(\varepsilon = 0.01\), \(\kappa=4\), \(M=0.001\), \(C_0=100\), \(T=0.1\), and the reference solution is computed by the EOP-IEQ/BDF2 scheme for \(\tau=1e-6\). The numerical errors and convergence rates of EOP-IEQ method are shown in Fig 4. We observed that the numerical results of the numerical solution \(\phi\) and the auxiliary variable \(q\) in \(L^2\) errors obtained by using the proposed schemes are in agreement with expectations.

Case B. Considering that coarsening effect is an important phenomenon in phase field simulation, this example is used to verify the validity of using energy-optimized technique to simulates long-time dynamics. The initial data is set to be the same as 27 and we adjust some parameters to \((N_x,N_y)=(128,128)\), \(\tau=1e-3\), \(M=0.1\) and \(T=100\). The results of the EOP-IEQ/CN scheme of \(\tau=1e-5\) are calculated as the reference solution. Fig. 5 shows several profiles of coarsening process under the EOP-IEQ/CN scheme at \(t=1, 10, 90\). We use the IEQ/CN and EOP-IEQ/CN schemes in Fig. 6 to investigate the modified energy, the ratio of nonlinear free energy and mass profiles during the evolution process. We observe that the modified energy obtained by the EOP-IEQ method is closer to the original energy than that obtained by the baseline IEQ method. It indicates the energy-optimized technique improves the accuracy significantly.

Case C. In the following examples, we will use the second-order EOP-IEQ/CN scheme [EOP47CN] for numerical simulation if not specified. We study the evolution of isolated ellipses and circles [32]. The initial condition is \[\phi \left( x,y,0 \right) =\tanh\left( \frac{dist\left( \left( x,y \right) ,\varGamma \right)}{\sqrt{2}\varepsilon} \right) ,\] where \(dist\left( \left( x,y \right) ,\varGamma \right)\) represents the signed distance from point \(\left( x,y \right)\) to the ellipse and circular interfaces \(\varGamma\). The center of the ellipse is located at \((-0.1, -0.1 )\), its major and minor axes are \(\sqrt{2}/5\), \(\sqrt{2}/10\), the central position of the circle is \((0.25,0.25)\), and the radius is set to \(0.1\). Other parameters are \((N_x,N_y)=(128,128)\), \(\alpha_0=0.01^2\), \(\tau=1e-4\), \(\varepsilon = 0.01\), \(\kappa=4\), \(M=1\), \(C_0=100\) and \(T=1\). We use the results of the EOP-IEQ/CN scheme of \(\tau=1e-6\) as the reference solution. Fig. 7 shows the morphological evolution under the IEQ/CN scheme with or without energy-optimized technique, and it can also be seen that the small circle is gradually absorbed by the ellipse until it completely disappears. We observe that the numerical results of the EOP-IEQ/CN scheme are closer to the reference solution from the enlarged figures. At the same time, it can be observed that the modified energy using EOP-IEQ/CN scheme is closer to the reference energy from the energy curve in Fig. 8 (a), (b), while the baseline IEQ method has significant deviation. In addition, we observe that the baseline IEQ method can maintain the conservation of mass regardless of the energy-optimized technique in Fig. 8 (c).

4.3 Phase field crystal model↩︎

Next, we investigate the phase field crystal (PFC) model [39], [40], considering the free energy \[\begin{align} \label{E-PFC} E\left( \phi \right) =\int_{\Omega}{\left( \frac{1}{4}\phi ^4+\frac{\beta-\varepsilon }{2}\phi ^2-\beta\left| \nabla \phi \right|^2+\frac{1}{2}\left( \Delta \phi \right) ^2 \right) \,\text{d}\boldsymbol{x}}, \end{align}\tag{28}\] where \(\beta\) the positive parameter and \(\beta>\varepsilon\). Then the PFC equation can be expressed in the energy variational form 2 of the free energy 28 with \(\mathcal{G} = -M\Delta\) as follows \[\begin{align} {}&\phi _t=-M\Delta \mu, \\ {}&\mu =\phi ^3+\left( \beta-\varepsilon \right) \phi +2\beta\Delta \phi +\Delta ^2\phi . \end{align}\] For the baseline IEQ method, we introduce the auxiliary variable \(q:=Q\left( \phi \right) =\phi^2\), then the general system 5 is specified as \[\label{PFC-IEQ} \begin{align} {}&\frac{\partial \phi}{\partial t}=M\Delta \mu ,\\ {}&\mu =q\phi +\left( \beta-\varepsilon \right) \phi +2\beta\Delta \phi +\Delta ^2\phi ,\\ {}&\frac{\partial q}{\partial t}=2\phi \frac{\partial \phi}{\partial t}. \end{align}\tag{29}\] In addition, the PFC equation also satisfies the law of conservation of mass, i.e. 26 is valid.

We simulate the complex dynamics of the growth of a polycrystal in a two-dimensional supercooled liquid by using the EOP-IEQ method to further highlight its power for simulating long-time dynamics, involving the motion of the liquid crystal interface and the grain boundaries separating the crystals. To define the initial configuration, we set up three microcrystals with different orientations as shown in Fig. 9 (a) in the computing domain \([0,800]^2\). Other parameters are \((N_x,N_y)=(512,512)\), \(\varepsilon =0.25\), \(\beta=1\), \(M=1\), \(\tau=1e-2\) and \(T=1500\). We use the results of the EOP-IEQ/CN scheme of \(\tau=1e-4\) as the reference solution. Fig. 9 shows microstructure evolution dynamics for triangular crystal phases driven by the PFC model using the EOP-IEQ/CN scheme. The results also show that the different arrangements of crystals lead to defects and dislocations, similar to those reported in the literature [27]. The comparison of the numerical results of discrete modified energies is shown in Fig. 10 (a). The EOP-IEQ/CN scheme is more accurate than the baseline IEQ/CN scheme in Fig. 10 (b). Most importantly, the ratio of nonlinear free energy \(\lVert Q\left( \phi ^n \right) \rVert ^2/\lVert q^n \rVert ^2=1\) by using the energy-optimized technique in Fig. 10 (b), which shows that the EOP-IEQ method preserves the consistency of \(Q( \phi ^n)\) and \(q^n\). Finally, Fig. 10 (c) indicates that the IEQ/CN scheme can ensure the conservation of mass regardless of the energy-optimized technique.

4.4 Molecular beam epitaxy model↩︎

Finally, we focus on the molecular beam epitaxy (MBE) model [41]–[43] without slope section. Considering the Ehrlich-Schwoebel energy \[\begin{align} \label{E-MBE} E\left( \phi \right) =\int_{\Omega}{\left( \frac{\varepsilon ^2}{2}\left( \Delta \phi \right) ^2-\frac{1}{2}\ln \left( 1+\left| \nabla \phi \right|^2 \right) \right) \,\text{d}\boldsymbol{x}}. \end{align}\tag{30}\] Then the MBE equation can be expressed in the energy variational form 2 of the free energy 30 with \(\mathcal{G} = MI\), which reads as \[\frac{\partial \phi}{\partial t}=-M\left( \varepsilon ^2\Delta ^2\phi +\nabla \cdot \left( \frac{1}{1+\left| \nabla \phi \right|^2} \nabla \phi \right) \right) .\] By introducing the auxiliary variable \(q :=Q\left( \phi \right) =\sqrt{\ln \left( 1+\left| \nabla \phi \right|^2 \right) +C_0}\), the MBE model can be equivalently written in quadratic form as follows \[\label{MBE-IEQ} \begin{align} {}&\frac{\partial \phi}{\partial t}=-M\left( \varepsilon ^2\Delta ^2\phi +\nabla \cdot \left( q\boldsymbol{H}(\phi) \right) \right) ,\\ {}&\frac{\partial q}{\partial t}=\boldsymbol{H}(\phi)\cdot \nabla \frac{\partial \phi}{\partial t}, \end{align}\tag{31}\] where \[\boldsymbol{H}(\phi)=\frac{\nabla \phi}{\left( 1+\left| \nabla \phi \right|^2 \right) Q(\phi) }.\] Then the modified energy is \[\bar{E}\left( \phi ,q \right) =\frac{\varepsilon ^2}{2}\lVert \Delta \phi \rVert ^2-\frac{1}{2}\lVert q \rVert ^2+\frac{1}{2}C_0\left| \Omega \right|.\]

The energy-optimized technique for the MBE model without slope section 31 is somewhat different from that of the general scheme [EOP47CN], so we reconsider the EOP-IEQ method for the equivalent system 31 based on the Crank-Nicolson formula.

The proof of the Theorem 7 is not particularly difficult but will be omitted. We define the roughness measurement function \(W(t)\) [43] as follows \[W\left( t \right) =\frac{1}{\left| \Omega \right|^{\frac{1}{2}}}\lVert \phi \left( \boldsymbol{x},t \right) -\bar{\phi}\left( t \right) \rVert , \quad \bar{\phi}\left( t \right) =\frac{1}{|\Omega|}\int_{\Omega}{\phi \left( \boldsymbol{x},t \right)}\,\text{d}\boldsymbol{x}.\]

In the last example, we also use the EOP-IEQ method and the baseline IEQ method to solve the MBE equation. Numerical simulations are carried out in the rectangular region \([-\pi, \pi] ^2\). We start from the classical benchmark problem \[\phi(x,y,0) = 0.1(\sin{3x}\sin{3y}+\cos{5x}\cos{5y}).\] The parameters are \((Nx,Ny) = (128,128)\), \(M=0.5\), \(\varepsilon = 0.1\), \(C_0=1\), \(\tau = 1e-2\) and \(T=50\). We use the results of the EOP-IEQ/CN scheme of \(\tau=1e-4\) as the reference solution. Several profiles of simulation process under the EOP-IEQ/CN scheme at \(t=1, 10, 45\) in Fig. 11. We summarize the normalized numerical results of the modified energy evolution in Fig. 12 (a). It can be observed that the EOP-IEQ/CN scheme provides significantly more accurate results than the baseline IEQ/CN scheme in Fig. 12 (b). Furthermore, Fig. 12 shows that the modified energy and roughness decay very quickly at the beginning, and after a long period of time, the roughness begins to increase. In general, the results are consistent with those in [41].

5 Conclusion↩︎

In this paper, we propose an EOP-IEQ method from the perspective of original energy, and construct a new class of first- and second-order (in time) unconditionally energy stable schemes, which inherit the advantages of the baseline IEQ and REQ schemes. The novel energy-optimized technique avoids the nonlinear optimization problems caused by the correction of auxiliary variables in REQ method, so our process of updating auxiliary variables is simpler and more efficient. Ample numerical examples of gradient flow models, including AC equation, CH equation, PFC equation and MBE equation, verify the accuracy, efficiency and energy stability of the proposed schemes. Furthermore, our proposed numerical schemes are not only linearly solvable, but also have no significant deviation between the numerical solutions and the reference solutions. On the one hand, we observe that the modified energy using energy-optimized technique is significantly more closer to the original energy than the relaxation technique and baseline calculated step. On the other hand, in most cases of long-term numerical simulation, the ratio of nonlinear free energy \(q^n=Q( \phi ^n ), n\ge 0\) by using the EOP-IEQ method, indicating that the numerical solutions obtained by the EOP-IEQ method preserve the original energy dissipation law. Overall, the main conclusion of this paper is that the energy-optimized technique is applicable to all available IEQ schemes in the literature with the natures of ease to use and improved accuracy.

CRediT authorship contribution statement↩︎

Xiaoqing Meng: Methodology, Software, Validation, Formal analysis, Investigation, Writing-Original Draft, Visualization. Aijie Cheng: Validation, Formal analysis, Resources, Writing-Review & Editing, Supervision. Zhengguang Liu: Conceptualization, Methodology, Validation, Formal analysis, Writing-Review & Editing, Supervision.

Declaration of competing interest↩︎

The authors declare that there are no conflicts of interest.

Data availability↩︎

Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

Acknowledgment↩︎

This work was supported by Natural Science Outstanding Youth Fund of Shandong Province (Grant No: ZR2023YQ007).

References↩︎