Fourier-Spectral Method for the Phase-Field Equations

Abstract

In this paper, we review the Fourier-spectral method for some phase-field models: Allen–Cahn (AC), Cahn–Hilliard (CH), Swift–Hohenberg (SH), phase-field crystal (PFC), and molecular beam epitaxy (MBE) growth. These equations are very important parabolic partial differential equations and are applicable to many interesting scientific problems. The AC equation is a reaction-diffusion equation modeling anti-phase domain coarsening dynamics. The CH equation models phase segregation of binary mixtures. The SH equation is a popular model for generating patterns in spatially extended dissipative systems. A classical PFC model is originally derived to investigate the dynamics of atomic-scale crystal growth. An isotropic symmetry MBE growth model is originally devised as a method for directly growing high purity epitaxial thin film of molecular beams evaporating on a heated substrate. The Fourier-spectral method is highly accurate and simple to implement. We present a detailed description of the method and explain its connection to MATLAB usage so that the interested readers can use the Fourier-spectral method for their research needs without difficulties. Several standard computational tests are done to demonstrate the performance of the method. Furthermore, we provide the MATLAB codes implementation in the Appendix A.

1. Introduction

In this paper, we review the unconditionally gradient stable Fourier-spectral method for several phase-field equations. The phase-field models are derived from the variational approach that minimizes the free energy of the system. Thus, the derived system follows a law of energy dissipation, which configures thermodynamic consistency, and hence leads to a mathematically well-posed model. The spectral methods, which belong to the collection of weighted residual methods, are originally derived to solve the spatial part of partial differential equations. For detailed, rigorous, and numerical aspect information on spectral methods, see [1,2] and the references therein. In spectral methods, one usually takes a finite set of the expansion functions (so called trial or basis functions) to represent the numerical solution as a linear combination of those functions. Choosing the expansion parts is important because they form a basis for entire spectral domain. Furthermore, we implicitly assume that those functions are smooth, therefore the most common choices are orthogonal polynomials and trigonometric functions. Since our focus is the Fourier-spectral method, the trigonometric basis functions are regarded as the trial functions. Moreover, instead of regarding entire continuous space, the representation is imposed only at discrete points; this is why this method is so called pseudo-spectral method. It is enable to one can save the computing resources in evaluating differentiation and employ the efficient algorithm such as the fast Fourier transform when the number of grid points is even. In addition, the error is decreasing exponentially when one employs the Fourier-spectral method while the finite difference methods lead to the algebraically decreasing order in fact. Figure 1 shows the order of spectral accuracy compared with the order of accuracy of finite difference. In this case, we use a simple function $u\left(x\right)=log(10+cosx)$; hence $u^{\prime }\left(x\right)=-sinx/(10+cosx)$. An error is defined as $\mathrm{Error}={max}_{1\le i\le N}|u^{\prime }\left(x_i\right)-D{u}_i|$ where $D{u}_i$ is an approximation of $u^{\prime }\left(x_i\right)$.

According to Figure 1, the error is decreasing exponentially via Fourier-spectral method indeed. Furthermore, we verify the efficiency of the Fourier-spectral method employing the following initial-boundary value problem in one-dimensional space with the $2\pi$-periodic boundary condition.

$$\begin{array}{c}\hfill u_t=u_{xx},u(x,0)=sin\left(x\right);(x,t)\in [0,2\pi ]\times (0,1].\end{array}$$

(1)

Then, the exact solution of Equation (1) is $u(x,t)=e^{-t}sin\left(x\right)$. Define $u_i^k$ as a numerical approximation to $u(x_i,k\Delta t)$ for $1\le k\le N_t$ where $\Delta t$ is a time step and $N_t$ is the number of iterations. Figure 2 depicts the convergence of second-order central finite difference scheme and Fourier-spectral method. Note that the backward difference method in time is employed to both cases and an error is defined as $\mathrm{Error}=\sqrt{{\sum }_{i=1}^N{(u(x_i,1)-u_i^{N_t})}^2/N}$ where N is the number of grid points, hence the number of modes, and it is easily verified that the Fourier-spectral method requires relatively few grid points within similar accuracy indeed.

Sometimes, an interpolant cannot apply to less flexible domains since the method employs implicitly periodic or homogeneous boundary conditions though there is an advantage in numerical accuracy. Therefore, a somewhat subtle additional process is required to employ complex or mixed boundary conditions. Moreover, it causes Gibbs phenomenon if there is a point of discontinuity, i.e., when it has a sharp gradient. Figure 3 shows the oscillations near discontinuities when one employs the Fourier interpolant.

Moreover, we set all the mobility terms in phase-field models as constant for convenience, more precisely set to 1 without loss of generality, since it requires further manipulation to handle a variable mobility in spectral methods [3].

The Allen–Cahn (AC) equation is discussed first. It describes the temporal evolution of a non-conserved order field during anti-phase domain coarsening [4]:

$$\begin{array}{c}\hfill \frac{\partial \varphi (\mathbf{x},t)}{\partial t}=-\frac{F^{\prime }\left(\varphi \right(\mathbf{x},t\left)\right)}{{ϵ}^2}+\Delta \varphi (\mathbf{x},t),\mathbf{x}\in \Omega ,t\ge 0,\end{array}$$

(2)

where $\varphi (\mathbf{x},t)$ is defined as the difference between the concentrations of the two components in a mixture, which varies on $-1$ to 1, $F\left(\varphi \right)=0.25{({\varphi }^2-1)}^2$ is the free energy potential, and $ϵ$ is a small constant related to the interfacial energy. Note that $F^{\prime }\left(\varphi \right)={\varphi }^3-\varphi$. The AC equation has a wide range of applications such as mean curvature flows, two-phase incompressible fluids, complex dynamics of dendritic growth, image inpainting, and image segmentation, see Figure 4 for some of these examples [5,6,7].

Lee and Lee [8] presented unconditionally energy stable first- and second-order semi-analytical Fourier-spectral methods for the AC equation to reduce the errors caused by large time step. Those methods were originated by decomposition of original equation into two parts, linear and nonlinear, that have closed-form solutions; hence operator splitting steps ensure high-order nonlinear solvers.

The next is a similar type to the previous one, the Cahn–Hilliard (CH) equation, which models the process of spinodal decomposition in conserved binary alloys [9]:

$$\begin{array}{c}\hfill \frac{\partial \varphi (\mathbf{x},t)}{\partial t}=\Delta \left(F^{\prime }\left(\varphi (\mathbf{x},t)\right)-{ϵ}^2\Delta \varphi \right),\mathbf{x}\in \Omega ,t\ge 0.\end{array}$$

(3)

The CH equation is widely used in applications such as phase separation, topology optimization, multiphase incompressible fluid flows, image inpainting, surface reconstruction, diblock copolymer, tumor growth simulation, and microstructures with elastic inhomogeneity, see Figure 5 for some of these examples [10].

Recently in [11], the discretization via nonlinear stabilized splitting scheme to the CH equation was reviewed and was solved by using a nonlinear multigrid method. Further in [12], Lee researched energy stability of the second-order strong-stability-preserving implicit-explicit Runge–Kutta methods for the CH equation. Christlieb et al. [13] presented the unconditionally gradient nonlinearly stabilized method for the CH equation, which is originally proposed by Eyre [14], within Fourier method and proposed an iterative scheme which is convergent for large time steps. There are a bunch of research related in AC and CH equations so far, some selected literatures are listed as follows. Montanelli and Bootland [15] proposed several exponential integration formula and compared their performance within stiff partial differential equations including AC and CH models. Such models are rewritten to sum of linear operator part with high-order terms and nonlinear operator part, and then Fourier-spectral method is applied in order to employ exponential integrator to this semilinear ordinary differential equations. Zhang and Liu [16] used several AC or CH type equations to represent the spatial patterns in ecological and biological system. Shen and Yang [17] presented numerical approximations of the AC and CH equations for semi-implicit or implicit schemes which are unconditionally energy stable, with stability analysis and error estimates based on spectral-Galerkin method. The results confirmed that spectral methods are suitable for diffusive interface models. Regarding on this respect especially, we introduce a nonlocal CH equation, which is appropriate to apply the Fourier-spectral method, that can model microphase separation in diblock copolymers, which consist of two different types of monomers [18] and an explicit form is listed as follows:

$$\begin{array}{c}\hfill \frac{\partial \varphi (\mathbf{x},t)}{\partial t}=\Delta \left(F^{\prime }\left(\varphi (\mathbf{x},t)\right)-{ϵ}^2\Delta \varphi \right)-\sigma \left(\varphi (\mathbf{x},t)-\overline{\varphi }\right),\mathbf{x}\in \Omega ,t\ge 0,\end{array}$$

(4)

where $\sigma$ is inversely proportional to the square of the total chain length of the copolymer and $\overline{\varphi }={\int }_{\Omega }\varphi (\mathbf{x},0)/|\Omega |$ is the average concentration over the domain $\Omega$. Block copolymer is a linear-chain molecule consists of at least two subchains connected to each other to make a polymer chain. A diblock copolymer exists if the subchain consists of two distinct monomer blocks. Related mathematical models have been developed in order to investigate the behaviors of phase separation of block copolymers and to find an available technique to manufacture nano-structured materials [19,20,21,22]. There is a direct applications of Equation (4) in [23] recently, where the authors employ the spectral method and see the references therein to check more details.

The Swift–Hohenberg (SH) equation was originally derived to model patterns from the influence of thermal fluctuations in hydrodynamics [24]:

$$\begin{array}{c}\hfill \frac{\partial \varphi (\mathbf{x},t)}{\partial t}=-\left({\varphi }^3+(1-ϵ)\varphi +2\Delta \varphi +{\Delta }^2\varphi \right),\mathbf{x}\in \Omega ,t\ge 0,\end{array}$$

(5)

where $\varphi$ is a density field and $ϵ$ is a temperature related positive constant. There are several applications to employ this model such as cellular materials, metallurgy, laser dynamics, electrohydrodynamics, crystallography, etc. See [25,26,27,28] and the references therein for more details. In addition, we introduce a classical phase field crystal (PFC) equation which takes account for atomic-crystallization growth [29]:

$$\begin{array}{c}\hfill \frac{\partial \varphi (\mathbf{x},t)}{\partial t}=\Delta \left[{\varphi }^3+(1-ϵ)\varphi +2\Delta \varphi +{\Delta }^2\varphi \right],\mathbf{x}\in \Omega ,t\ge 0.\end{array}$$

(6)

This model has a variety of applications such as crystallization in liquid-liquid interface with undercooled material, isotropic phase separation, etc. On the atomistic/molecular scale freezing, a theoretical approach to undercooled liquids crystallization has been studied in [30,31], and the efficient numerical methods based on the operator splitting method or spectral method are developed for the phase-field crystal model [32,33]. The operator splitting method with Fourier-spectral method can relax the time step restriction or shorten the computation time depending on which a solver is applied for each stage.

A molecular beam epitaxy (MBE) growth model describes a process in which a thin single crystal layer is deposited on a single crystal substrate using molecular beams [34,35]:

$$\begin{array}{c}\hfill \frac{\partial \varphi (\mathbf{x},t)}{\partial t}=-ϵ{\Delta }^2\varphi -\nabla ·\left[\left(\right|1-\nabla \varphi {|}^2)\nabla \varphi \right],\mathbf{x}\in \Omega ,t\ge 0.\end{array}$$

(7)

The MBE model is a substantially used approach for thin-film deposition of a surface or interface quality determined to single-monolayer precision. This procedure is widely applied in semiconductor heterostructures and persistently studied topic in material science. Consequently, considerable mathematical models have been evolved to study the epitaxy dynamics, covering from continuum models to molecular dynamical simulations. The readers are referred to the following references for more details [36,37,38,39,40,41,42,43].

The main purpose of this paper is to present brief reviews, to describe numerical solution algorithms, and to provide the MATLAB code implementations of the unconditionally gradient stable Fourier-spectral method for the several phase-field equations. In particular, we highlight the caution that needs to be taken when applying the MATLAB based fast Fourier transform to the Fourier-spectral method.

The outline of this paper is as follows. The numerical solutions in two- and three-dimensional cases of the above phase-field models are described in Section 2 and Section 4, respectively. In Section 3 and Section 5, we present the basic numerical simulations to the stated phase-field models in both two- and three-dimensional cases. We finalize the paper with the conclusion in Section 6. In the Appendix A, we provide the MATLAB codes for the numerical implementation of the presented equations.

2. Numerical Solutions in 2D

In this section, we present unconditionally stable Fourier-spectral methods for the phase-field models in two-dimensional space $\Omega =[l_x,r_x]\times [l_y,r_y]$. Let $N_x$, $N_y$ be positive even integers and $L_x=r_x-l_x$, $L_y=r_y-l_y$ be the length of each direction, respectively; hence define $h_x=L_x/N_x$ and $h_y=L_y/N_y$ as the spatial step size in each direction, respectively. We denote discretized points as $(x_m,y_n)=(l_x+m{h}_x,l_y+n{h}_y)$ where $0\le m\le N_x$ and $0\le n\le N_y$ are integers. Let ${\varphi }_{mn}^k$ be an approximation of $\varphi (x_m,y_n,t_k)$, where $t_k=k\Delta t$ and $\Delta t$ is the temporal step size. For the given data $\left\{{\varphi }_{mn}^k\right|m=1,\dots ,N_x$ and $n=1,\dots ,N_y\}$, the discrete Fourier transform is defined as

$$\begin{array}{c}\hfill {\widehat{\varphi }}_{pq}^k=\sum _{m=1}^{N_x}\sum _{n=1}^{N_y}{\varphi }_{mn}^k{e}^{-i({\xi }_p{x}_m+{\eta }_q{y}_n)},-\frac{N_x}{2}+1\le p\le \frac{N_x}{2},-\frac{N_y}{2}+1\le q\le \frac{N_y}{2},\end{array}$$

(8)

where ${\xi }_p=2\pi p/L_x$ and ${\eta }_q=2\pi q/L_y$. The inverse discrete Fourier transform is

$$\begin{array}{c}\hfill {\varphi }_{mn}^k=\frac{1}{N_x{N}_y}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}{\widehat{\varphi }}_{pq}^k{e}^{i({\xi }_p{x}_m+{\eta }_q{y}_n)}.\end{array}$$

(9)

Note that we can obtain spectral derivatives as if we perform analytic differentiations in the Fourier space. We assume that $\varphi (x,y,t)$ is sufficiently smooth and extended to continuous version of the numerical approximation ${\varphi }_{mn}^k$. The following shows step-by-step description of how the differentiation works in the Fourier transform with finite basis.

$$\begin{array}{c}\begin{array}{cc}\frac{\partial }{\partial x}\varphi (x,y,t)\hfill & =\frac{1}{N_x{N}_y}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\left(i{\xi }_p\right)\widehat{\varphi }({\xi }_p,{\eta }_q,t)e^{i({\xi }_{p}x+{\eta }_{q}y)},\hfill \\ \hfill \frac{\partial }{\partial y}\varphi (x,y,t)& =\frac{1}{N_x{N}_y}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\left(i{\eta }_q\right)\widehat{\varphi }({\xi }_p,{\eta }_q,t)e^{i({\xi }_{p}x+{\eta }_{q}y)}.\hfill \end{array}\end{array}$$

(10)

This process enables one can derive spectral derivatives in the Fourier space easily, not differentiate directly in the physical space. Therefore, we can represent the Laplacian to coefficients in the Fourier space as follows:

$$\begin{array}{c}\begin{array}{cc}\Delta \varphi (x,y,t)\hfill & =\frac{1}{N_x{N}_y}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\widehat{\Delta \varphi }({\xi }_p,{\eta }_q,t)e^{i({\xi }_{p}x+{\eta }_{q}y)}\hfill \\ & =\frac{1}{N_x{N}_y}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}-({\xi }_p^2+{\eta }_q^2)\widehat{\varphi }({\xi }_p,{\eta }_q,t)e^{i({\xi }_{p}x+{\eta }_{q}y)},\hfill \end{array}\hfill \end{array}$$

(11)

where the first line is the definition of the inverse Fourier transform and the second line is just applying Equation (10) twice to x- and y-direction to $\varphi (x,y,t)$ and its Fourier transform.

Now we present the numerical solutions of phase-field equations. First, we derive the numerical solution of the AC equation. We apply the linearly stabilized splitting scheme [14] to Equation (2).

$$\begin{array}{c}\hfill \frac{{\varphi }_{mn}^{k+1}-{\varphi }_{mn}^k}{\Delta t}=-\frac{2{\varphi }_{mn}^{k+1}+f\left({\varphi }_{mn}^k\right)}{{ϵ}^2}+{(\Delta {\varphi }^{k+1})}_{mn},\end{array}$$

(12)

where $f\left(\varphi \right)={\varphi }^3-3\varphi$. Thus, Equation (12) can be transformed into the discrete Fourier space as follows:

$$\begin{array}{c}\hfill \frac{{\widehat{\varphi }}_{pq}^{k+1}-{\widehat{\varphi }}_{pq}^k}{\Delta t}=-\frac{2{\widehat{\varphi }}_{pq}^{k+1}+{\widehat{f}}_{pq}^k}{{ϵ}^2}-\left({\xi }_p^2+{\eta }_q^2\right){\widehat{\varphi }}_{pq}^{k+1}.\end{array}$$

(13)

Therefore, we obtain the following discrete Fourier transform

$$\begin{array}{c}\hfill {\widehat{\varphi }}_{pq}^{k+1}=\frac{{ϵ}^2{\widehat{\varphi }}_{pq}^k-\Delta t{\widehat{f}}_{pq}^k}{{ϵ}^2+\Delta t[2+{ϵ}^2({\xi }_p^2+{\eta }_q^2)]}.\end{array}$$

(14)

Then, the updated numerical solution ${\varphi }_{mn}^{k+1}$ can be computed using Equation (9):

$$\begin{array}{c}\hfill {\varphi }_{mn}^{k+1}=\frac{1}{N_x{N}_y}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}{\widehat{\varphi }}_{pq}^{k+1}{e}^{i({\xi }_p{x}_m+{\eta }_q{y}_n)}.\end{array}$$

(15)

Next, we obtain the numerical solution of the CH equation. We employ the linearly stabilized splitting scheme [14] to Equation (3).

$$\begin{array}{c}\hfill \frac{{\varphi }_{mn}^{k+1}-{\varphi }_{mn}^k}{\Delta t}=\Delta {\left(2{\varphi }_{mn}^{k+1}-{ϵ}^2{(\Delta {\varphi }^{k+1})}_{mn}+f\left({\varphi }_{mn}^k\right)\right)}_{mn}.\end{array}$$

(16)

Thus, Equation (16) can be transformed into the discrete Fourier space as follows:

$$\begin{array}{c}\hfill \frac{{\widehat{\varphi }}_{pq}^{k+1}-{\widehat{\varphi }}_{pq}^k}{\Delta t}=-\left({\xi }_p^2+{\eta }_q^2\right)\left(2{\widehat{\varphi }}_{pq}^{k+1}+{ϵ}^2({\xi }_p^2+{\eta }_q^2){\widehat{\varphi }}_{pq}^{k+1}+{\widehat{f}}_{pq}^k\right).\end{array}$$

(17)

Therefore, we obtain the following discrete Fourier transform

$$\begin{array}{c}\hfill {\widehat{\varphi }}_{pq}^{k+1}=\frac{{\widehat{\varphi }}_{pq}^k-({\xi }_p^2+{\eta }_q^2)\Delta t{\widehat{f}}_{pq}^k}{1+\Delta t[2({\xi }_p^2+{\eta }_q^2)+{ϵ}^2{({\xi }_p^2+{\eta }_q^2)}^2]}.\end{array}$$

(18)

Then, the updated numerical solution ${\varphi }_{mn}^{k+1}$ can be computed using Equation (9):

$$\begin{array}{c}\hfill {\varphi }_{mn}^{k+1}=\frac{1}{N_x{N}_y}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}{\widehat{\varphi }}_{pq}^{k+1}{e}^{i({\xi }_p{x}_m+{\eta }_q{y}_n)}.\end{array}$$

(19)

Now, we present a numerical solution to the SH equation. In a similar manner, we discretize Equation (5) as follows:

$$\begin{array}{c}\hfill \frac{{\varphi }_{mn}^{k+1}-{\varphi }_{mn}^k}{\Delta t}=-(1-ϵ){\varphi }_{mn}^{k+1}-{\left({\Delta }^2{\varphi }^{k+1}\right)}_{mn}-g\left({\varphi }_{mn}^k\right)-2{(\Delta {\varphi }^k)}_{mn},\end{array}$$

(20)

where $g\left(\varphi \right)={\varphi }^3$. Then we transform Equation (20) as

$$\begin{array}{c}\hfill \frac{{\widehat{\varphi }}_{pq}^{k+1}-{\widehat{\varphi }}_{pq}^k}{\Delta t}=-(1-ϵ){\widehat{\varphi }}_{pq}^{k+1}-{({\xi }_p^2+{\eta }_q^2)}^2{\widehat{\varphi }}_{pq}^{k+1}-{\widehat{g}}_{pq}^k+2({\xi }_p^2+{\eta }_q^2){\widehat{\varphi }}_{pq}^k.\end{array}$$

(21)

Therefore, we have the following result

$$\begin{array}{c}\hfill {\widehat{\varphi }}_{pq}^{k+1}=\frac{{\widehat{\varphi }}_{pq}^k+\Delta t[2({\xi }_p^2+{\eta }_q^2){\widehat{\varphi }}_{pq}^k-{\widehat{g}}_{pq}^k]}{1+\Delta t[(1-ϵ)+{({\xi }_p^2+{\eta }_q^2)}^2]},\end{array}$$

(22)

and hence we have a numerical solution ${\varphi }_{mn}^{k+1}$ as follows:

$$\begin{array}{c}\hfill {\varphi }_{mn}^{k+1}=\frac{1}{N_x{N}_y}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}{\widehat{\varphi }}_{pq}^{k+1}{e}^{i({\xi }_p{x}_m+{\eta }_q{y}_n)}.\end{array}$$

(23)

Similarly, a numerical solution for the PFC model (6) is obtained by the same procedure

$$\begin{array}{c}\hfill \frac{{\varphi }_{mn}^{k+1}-{\varphi }_{mn}^k}{\Delta t}=\Delta {\left((1-ϵ){\varphi }_{mn}^{k+1}+{\left({\Delta }^2{\varphi }^{k+1}\right)}_{mn}+g\left({\varphi }_{mn}^k\right)+2{(\Delta {\varphi }^k)}_{mn}\right)}_{mn}.\end{array}$$

(24)

Then Equation (24) is transformed as

$$\begin{array}{c}\hfill \frac{{\widehat{\varphi }}_{pq}^{k+1}-{\widehat{\varphi }}_{pq}^k}{\Delta t}=-\left({\xi }_p^2+{\eta }_q^2\right)\left((1-ϵ){\widehat{\varphi }}_{pq}^{k+1}+{({\xi }_p^2+{\eta }_q^2)}^2{\widehat{\varphi }}_{pq}^{k+1}+{\widehat{g}}_{pq}^k-2({\xi }_p^2+{\eta }_q^2){\widehat{\varphi }}_{pq}^k\right).\end{array}$$

(25)

Therefore, we have the following result

$$\begin{array}{c}\hfill {\widehat{\varphi }}_{pq}^{k+1}=\frac{{\widehat{\varphi }}_{pq}^k+\Delta t[2{({\xi }_p^2+{\eta }_q^2)}^2{\widehat{\varphi }}_{pq}^k-({\xi }_p^2+{\eta }_q^2){\widehat{g}}_{pq}^k]}{1+\Delta t[(1-ϵ)({\xi }_p^2+{\eta }_q^2)+{({\xi }_p^2+{\eta }_q^2)}^3]}.\end{array}$$

(26)

Subsequently, we have a numerical solution ${\varphi }_{mn}^{k+1}$ as follows:

$$\begin{array}{c}\hfill {\varphi }_{mn}^{k+1}=\frac{1}{N_x{N}_y}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}{\widehat{\varphi }}_{pq}^{k+1}{e}^{i({\xi }_p{x}_m+{\eta }_q{y}_n)}.\end{array}$$

(27)

The remaining one is a numerical solution to the MBE growth model (7). Discretize Equation (7) with an expanded divergence term $\Delta \varphi$ and treat this in implicit way,

$$\begin{array}{c}\hfill \frac{{\varphi }_{mn}^{k+1}-{\varphi }_{mn}^k}{\Delta t}=-ϵ{\left({\Delta }^2{\varphi }^{k+1}\right)}_{mn}-{(\Delta {\varphi }^{k+1})}_{mn}+{\left(\nabla ·\left[{\left|{(\nabla {\varphi }^k)}_{mn}\right|}^2{(\nabla {\varphi }^k)}_{mn}\right]\right)}_{mn}.\end{array}$$

(28)

We define $b\left(\varphi \right)=\nabla ·[{|\nabla \varphi |}^2\nabla \varphi ]$. Then Equation (28) is transformed as

$$\begin{array}{c}\hfill \frac{{\widehat{\varphi }}_{pq}^{k+1}-{\widehat{\varphi }}_{pq}^k}{\Delta t}=-ϵ{({\xi }_p^2+{\eta }_q^2)}^2{\widehat{\varphi }}_{pq}^{k+1}+({\xi }_p^2+{\eta }_q^2){\widehat{\varphi }}_{pq}^{k+1}+{\widehat{b}}_{pq}^k,\end{array}$$

(29)

where ${\widehat{b}}_{pq}^k$ is

$$\begin{array}{ccc}{\widehat{b}}_{pq}^k\hfill & =& \left(\begin{array}{c}i{\xi }_p\\ i{\eta }_q\end{array}\right)·{\widehat{\mathbf{a}}}_{pq}^k,\hfill \end{array}$$

(30)

$$\begin{array}{ccc}{\widehat{\mathbf{a}}}_{pq}^k\hfill & =& {\left[\left({\left[\left(\begin{array}{c}i{\xi }_p{\widehat{\varphi }}_{pq}^k\\ i{\eta }_q{\widehat{\varphi }}_{pq}^k\end{array}\right)\right]}_{\mathrm{IFT}}·{\left[\left(\begin{array}{c}i{\xi }_p{\widehat{\varphi }}_{pq}^k\\ i{\eta }_q{\widehat{\varphi }}_{pq}^k\end{array}\right)\right]}_{\mathrm{IFT}}\right){\left[\left(\begin{array}{c}i{\xi }_p{\widehat{\varphi }}_{pq}^k\\ i{\eta }_q{\widehat{\varphi }}_{pq}^k\end{array}\right)\right]}_{\mathrm{IFT}}\right]}_{\mathrm{FT}}.\hfill \end{array}$$

(31)

Note that ${[·]}_{\mathrm{FT}}$ and ${[·]}_{\mathrm{IFT}}$ represent the discrete Fourier transform and the discrete inverse Fourier transform, respectively. Therefore, we have the following result

$$\begin{array}{c}\hfill {\widehat{\varphi }}_{pq}^{k+1}=\frac{{\widehat{\varphi }}_{pq}^k+\Delta t{\widehat{b}}_{pq}^k}{1+\Delta t[ϵ{({\xi }_p^2+{\eta }_q^2)}^2-({\xi }_p^2+{\eta }_q^2)]},\end{array}$$

(32)

and then we update a numerical solution ${\varphi }_{mn}^{k+1}$ as follows:

$$\begin{array}{c}\hfill {\varphi }_{mn}^{k+1}=\frac{1}{N_x{N}_y}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}{\widehat{\varphi }}_{pq}^{k+1}{e}^{i({\xi }_p{x}_m+{\eta }_q{y}_n)}.\end{array}$$

(33)

3. Numerical Experiments in 2D

We perform several numerical investigations in this section. Note that we employ the $2\pi$-periodic boundary condition for overall numerical simulations in two-dimensional case. Sometimes, we employ the initial conditions as random perturbations in domain since these are adequate to reaction-diffusion type models. For AC and CH equations, since there are equilibrium solutions in one-dimensional space as hyperbolic tangent profiles on infinite domain, we use those as initial conditions. Furthermore, we take the value of model parameter $ϵ$ as an approximate value close to ${ϵ}_m=hm/\left[2\sqrt{2}{tanh}^{-1}\left(0.9\right)\right]$, which is used in finite difference methods [44].

3.1. The AC Equation

Numerical simulations are conducted to verify the mean curvature flow of the AC equation at first. Initial conditions are given as follows:

$$\begin{array}{ccc}\hfill \varphi (x,y,0)& =& tanh\left(\frac{2-\sqrt{x^2+y^2}}{\sqrt{2}ϵ}\right),\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill \varphi (x,y,0)& =& tanh\left(\frac{1.7+1.2cos\left(6\theta \right)-\sqrt{x^2+y^2}}{\sqrt{2}ϵ}\right),\hfill \end{array}$$

(35)

where $\theta ={tan}^{-1}(y/x)$ for $(x,y)\in [-\pi ,\pi ]\times [-\pi ,\pi ]$. Figure 6 shows the numerical test results at $t=1$ to the AC equation using the Fourier-spectral method. Here, we use $N_x=N_y=128$, $h=2\pi /128$, $ϵ=0.05$, and $\Delta t=0.0001$. The final time is $T=1$. The results reflect the motion by mean curvature characteristic of the AC equation.

Figure 7 shows the zero-level contours over time with the initial conditions used in Figure 6.

Moreover, an energy dissipation of the AC model (2) is depicted in Figure 8, where a total free energy $\mathcal{E}\left(\varphi \right)$ is defined as

$$\begin{array}{c}\hfill \mathcal{E}\left(\varphi \right)={\int }_{\Omega }\left(\frac{F\left(\varphi \right)}{{ϵ}^2}+\frac{{|\nabla \varphi |}^2}{2}\right)d\mathbf{x}.\end{array}$$

(36)

Discretize Equation (36) with spectral derivatives, it yields

$$\begin{array}{c}\hfill {\mathcal{E}}^h\left({\varphi }^k\right)=\sum _{m=1}^{N_x}\sum _{n=1}^{N_y}\frac{F\left({\varphi }_{mn}^k\right)}{{ϵ}^2}{h}^2+\frac{h^2}{N_x{N}_y}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\frac{{\left|{\xi }_p{\widehat{\varphi }}_{pq}^k\right|}^2+{\left|{\eta }_q{\widehat{\varphi }}_{pq}^k\right|}^2}{2},\end{array}$$

(37)

since the following Parseval’s identity has been applied,

$$\begin{array}{c}\hfill \sum _{m=1}^{N_x}\sum _{n=1}^{N_y}{\left|{\varphi }_{mn}^k\right|}^2=\frac{1}{N_x{N}_y}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}{\left|\widehat{{\varphi }_{pq}^k}\right|}^2.\end{array}$$

(38)

3.2. The CH Equation

Next, we investigate the coarsening dynamics of the CH Equation (3) numerically with the following initial conditions,

$$\begin{array}{ccc}\hfill \varphi (x,y,0)& =& 0.05\mathrm{rand}(x,y),\hfill \end{array}$$

(39)

$$\begin{array}{ccc}\hfill \varphi (x,y,0)& =& -0.45+0.05\mathrm{rand}(x,y),\hfill \end{array}$$

(40)

where rand denotes a random number between $-1$ to 1 and $(x,y)\in [-0.5,0.5]\times [-0.5,0.5]$. Figure 9 shows the numerical results at $t=0.01$. We use $N_x=N_y=1000$, $h=1/1000$, $ϵ=0.0025$, and $\Delta t=0.00001$. The results represent well the coarsening dynamics of the CH equation.

Figure 10 shows the numerical simulation results at $t=2.5$ to the CH equation with the initial condition (40). Here, $N_x=N_y=1000$, $h=h_x=h_y$, $ϵ=0.0025$, and $\Delta t=0.001$ are used.

Furthermore, an energy dissipation of the CH equation is illustrated in Figure 11. We define a total free energy as follows:

$$\begin{array}{c}\hfill \mathcal{E}\left(\varphi \right)={\int }_{\Omega }\left(F\left(\varphi \right)+\frac{{ϵ}^2}{2}{|\nabla \varphi |}^2\right)d\mathbf{x}.\end{array}$$

(41)

Identical process represented in Equation (38) is adopted to Equation (41),

$$\begin{array}{c}\hfill {\mathcal{E}}^h\left({\varphi }^k\right)=\sum _{m=1}^{N_x}\sum _{n=1}^{N_y}F\left({\varphi }_{mn}^k\right)h^2+\frac{h^2}{N_x{N}_y}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\frac{{ϵ}^2}{2}({\left|{\xi }_p{\widehat{\varphi }}_{pq}^k\right|}^2+{\left|{\eta }_q{\widehat{\varphi }}_{pq}^k\right|}^2).\end{array}$$

(42)

We adopt the initial condition (39) and $\Delta t=0.0001$.

Further we present the phase separation behavior via nonlocal CH model for diblock copolymers. An initial condition is given by

$$\begin{array}{c}\hfill \varphi (x,y,0)=\overline{\varphi }+0.1\mathrm{rand}(x,y).\end{array}$$

(43)

A computational domain is given as $[0,2.4]\times [0,2.4]$. Here, we use $N_x=N_y=200$, $h=2.4/200$, $\overline{\varphi }=-0.3$, $\sigma =500$, $ϵ=h$, and $\Delta t=0.001$. The final time is $T=2$. Figure 12 shows the phase separation behaviors called the hexagonal patterns over time.

Another initial condition is given as follows:

$$\begin{array}{c}\hfill \varphi (x,y,0)=\overline{\varphi }+0.2\mathrm{rand}(x,y).\end{array}$$

(44)

A computational domain is given as $[0,3.2]\times [0,3.2]$. Parameter values are set to $N_x=N_y=300$, $h=3.2/300$, $\overline{\varphi }=0$, $\sigma =100$, $ϵ=3h$, and $\Delta t=0.1$. The final time is $T=1000$. Figure 13 describes the phase separation behaviors called the lamellar patterns over time.

Consequently, we illustrate an energy dissipation of the nonlocal CH model (4). A total free energy is given by

$$\begin{array}{c}\hfill \mathcal{E}\left(\varphi \right)={\int }_{\Omega }\left[\left(F\left(\varphi \right)+\frac{{ϵ}^2}{2}{|\nabla \varphi |}^2\right)+\frac{\sigma }{2}{|\nabla \psi |}^2\right]d\mathbf{x},\end{array}$$

(45)

where $\psi$ is obtained by solving $-\Delta \psi =\varphi -\overline{\varphi }$. Now Equation (45) is discretized as

$$\begin{array}{cc}{\mathcal{E}}^h\left({\varphi }^k\right)\hfill & ={\displaystyle \sum _{m=1}^{N_x}}{\displaystyle \sum _{n=1}^{N_y}}F\left({\varphi }_{mn}^k\right)h^2+\frac{h^2}{N_x{N}_y}{\displaystyle \sum _{p=-N_x/2+1}^{N_x/2}}{\displaystyle \sum _{q=-N_y/2+1}^{N_y/2}}\frac{{ϵ}^2}{2}({\left|{\xi }_p{\widehat{\varphi }}_{pq}^k\right|}^2+{\left|{\eta }_q{\widehat{\varphi }}_{pq}^k\right|}^2)\hfill \\ & +\frac{h^2}{N_x{N}_y}{\displaystyle \sum _{p=-N_x/2+1}^{N_x/2}}{\displaystyle \sum _{q=-N_y/2+1}^{N_y/2}}\frac{\sigma }{2}({\left|{\xi }_p{\widehat{\psi }}_{pq}^k\right|}^2+{\left|{\eta }_q{\widehat{\psi }}_{pq}^k\right|}^2),\hfill \end{array}$$

(46)

by using Equation (38). Figure 14 shows the energy dissipation of the nonlocal CH equation with the discrete energy (46). We employ the same condition in Figure 12 except for $\Delta t=0.0001$ and $T=0.1$.

3.3. The SH Equation

We present the pattern formation due to the SH Equation (5) in the rectangular domain $[0,120]\times [0,120]$. The initial condition is given as follows:

$$\begin{array}{c}\hfill \varphi (x,y,0)=\overline{\varphi }+\overline{\varphi }\mathrm{rand}(x,y),\end{array}$$

(47)

where $\overline{\varphi }=0.05$. Here, we choose parameter values as $N_x=N_y=140$, $h=120/140$, $ϵ=0.2$, and $\Delta t=0.01$. The final time is $T=1500$. Figure 15 depicts the pattern formation via instability of thermal fluctuations. Note that we fix the range of colors in $[\overline{\varphi }-0.6,\overline{\varphi }+0.6]$.

For the next step, we investigate an energy decay of the following functional

$$\begin{array}{c}\hfill \mathcal{E}\left(\varphi \right)={\int }_{\Omega }\left(\frac{1}{4}{\varphi }^4+\frac{(1-ϵ)}{2}{\varphi }^2-{|\nabla \varphi |}^2+\frac{1}{2}{(\Delta \varphi )}^2\right)d\mathbf{x}.\end{array}$$

(48)

Using Equation (38), a discrete version of Equation (48) is defined as

$$\begin{array}{cc}{\mathcal{E}}^h\left({\varphi }^k\right)=& \sum _{m=1}^{N_x}\sum _{n=1}^{N_y}\left(\frac{1}{4}{\left({\varphi }_{mn}^k\right)}^4+\frac{(1-ϵ)}{2}{\left({\varphi }_{mn}^k\right)}^2\right)\hfill \\ & +\frac{h^2}{N_x{N}_y}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\left(-({\left|{\xi }_p{\widehat{\varphi }}_{pq}^k\right|}^2+{\left|{\eta }_q{\widehat{\varphi }}_{pq}^k\right|}^2)+\frac{1}{2}{\left|({\xi }_p^2+{\eta }_q^2){\widehat{\varphi }}_{pq}^k\right|}^2\right).\end{array}$$

(49)

Figure 16 illustrates the energy dissipation of Equation (49) via SH equation.

3.4. Classical PFC Model

The phase separation behavior of the classical PFC model is described in this section. Numerical simulation is conducted in the rectangular domain $[0,120]\times [0,120]$. The initial condition is given as follows:

$$\begin{array}{c}\hfill \varphi (x,y,0)=\overline{\varphi }+\overline{\varphi }\mathrm{rand}(x,y),\end{array}$$

(50)

where $\overline{\varphi }=0.07$. Parameters used here are $N_x=N_y=120$, $h=120/120$, $ϵ=0.025$, and $\Delta t=0.001$. The final time is $T=1500$. Figure 17 illustrates the isotropic phase transition behaviors. Note that we restrict the range of colors in $[\overline{\varphi }-0.2,\overline{\varphi }+0.2]$.

Furthermore, we present an illustration of energy dissipation to this PFC model. Since both SH and PFC equations derived from the same total free energy, we employ Equation (49). Figure 18 depicts the energy dissipation over time with the corresponding time behaviors.

3.5. MBE Growth Model

The epitaxial growth of isotropic current MBE growth model is investigated in this section. Numerical simulation is conducted in the rectangular domain $[0,2\pi ]\times [0,2\pi ]$. The initial condition is given as follows:

$$\begin{array}{c}\hfill \varphi (x,y,0)=0.1(sin3xsin2y+sin5xsin5y).\end{array}$$

(51)

Parameters are used as $N_x=N_y=200$, $ϵ=0.1$, and $\Delta t=0.001$. The final time is $T=15$. Figure 19 illustrates the epitaxial growth in computational domain.

Further, we present an energy decay to the MBE growth model (7). An energy functional is given by

$$\begin{array}{c}\hfill \mathcal{E}\left(\varphi \right)={\int }_{\Omega }\left(\frac{ϵ}{2}{|\Delta \varphi |}^2+\frac{1}{4}{({|\nabla \varphi |}^2-1)}^2\right)d\mathbf{x}.\end{array}$$

(52)

Discretize Equation (52) using Equation (38),

$$\begin{array}{c}\hfill {\mathcal{E}}^h\left({\varphi }^k\right)=\frac{h^2}{N_x{N}_y}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\left(\frac{ϵ}{2}{|-({\xi }_p^2+{\eta }_q^2){\widehat{\varphi }}_{pq}^k|}^2+\frac{1}{4}{\left|{\widehat{g}}_{pq}^k\right|}^2\right),\end{array}$$

(53)

where ${\widehat{g}}_{pq}^k$ is a coefficient in Fourier space defined

$$\begin{array}{c}\hfill {\widehat{g}}_{pq}^k={\left[{\left({\left[\left(i{\xi }_p\right){\widehat{\varphi }}_{pq}^k\right]}_{IFT}\right)}^2+{\left({\left[\left(i{\eta }_q\right){\widehat{\varphi }}_{pq}^k\right]}_{IFT}\right)}^2-1\right]}_{FT}.\end{array}$$

(54)

Figure 20 depicts the discrete energy decay over time. Here, we employ the same condition used in Figure 19.

3.6. Logarithmic Free Energy

In this paper, we have used the polynomial form of double-well potential $F\left(\varphi \right)$. Note that it is not the only choice. Now, we consider the logarithmic Flory–Huggins free energy [45,46]

$$\begin{array}{c}\hfill F\left(\varphi \right)={\theta }_e\left[(1-\varphi )ln(1-\varphi )+(1+\varphi )ln(1+\varphi )\right]-{\theta }_c{\varphi }^2,\end{array}$$

(55)

for the AC and CH equations. Here, ${\theta }_e$ and ${\theta }_c$ is a positive constant, which are taken as 1 and 1.5, respectively, in this paper. Figure 21 illustrates the logarithmic free energy.

To solve the AC and CH equations with the logarithmic free energy (55), we also apply Eyre’s linearly stabilized splitting scheme [14]. First, we discretize the AC Equation (2) with a free energy (55):

$$\begin{array}{c}\hfill \frac{{\varphi }_{mn}^{k+1}-{\varphi }_{mn}^k}{\Delta t}=-\frac{2\alpha {\varphi }_{mn}^{k+1}+f\left({\varphi }_{mn}^k\right)}{{ϵ}^2}+{(\Delta {\varphi }^{k+1})}_{mn},\end{array}$$

(56)

where $f\left(\varphi \right)={\theta }_{e}ln\left(1-\varphi \right)-2(1+\alpha ){\theta }_c\varphi$, and $\alpha$ is an auxiliary constant for the splitting scheme [47]. Then, Equation (56) can be transformed into the discrete Fourier space as follows:

$$\begin{array}{c}\hfill \frac{{\widehat{\varphi }}_{pq}^{k+1}-{\widehat{\varphi }}_{pq}^k}{\Delta t}=-\frac{2\alpha {\widehat{\varphi }}_{pq}^{k+1}+{\widehat{f}}_{pq}^k}{{ϵ}^2}-\left({\xi }_p^2+{\eta }_q^2\right){\widehat{\varphi }}_{pq}^{k+1}.\end{array}$$

(57)

Therefore, we obtain the following discrete Fourier transform

$$\begin{array}{c}\hfill {\widehat{\varphi }}_{pq}^{k+1}=\frac{{ϵ}^2{\widehat{\varphi }}_{pq}^k-\Delta t{\widehat{f}}_{pq}^k}{{ϵ}^2+\Delta t[2\alpha {\theta }_e+{ϵ}^2({\xi }_p^2+{\eta }_q^2)]}.\end{array}$$

(58)

Then, the updated numerical solution ${\varphi }_{mn}^{k+1}$ can be computed using Equation (9):

$$\begin{array}{c}\hfill {\varphi }_{mn}^{k+1}=\frac{1}{N_x{N}_y}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}{\widehat{\varphi }}_{pq}^{k+1}{e}^{i({\xi }_p{x}_m+{\eta }_q{y}_n)}.\end{array}$$

(59)

We conduct a numerical test with the following initial condition is

$$\begin{array}{c}\hfill \varphi (x,y,0)=0.99tanh\left(\frac{1.7+1.2\cos \left(6\theta \right)-\sqrt{x^2{y}^2}}{\sqrt{2}ϵ}\right),\end{array}$$

(60)

where $\theta ={tan}^{-1}(y/x)$ for $(x,y)\in [\pi ,\pi ]\times [\pi ,\pi ]$. We set the parameters as $N_x=N_y=128,h=h_x=h_y,ϵ=0.05,\alpha =2,$ and $\Delta t=0.0001$. Figure 22 shows the zero-level contours over time with the initial condition (60) according to motion by mean curvature.

Next, we obtain the numerical solution of the CH Equation (3) with the logarithmic free energy (55). Applying the linearly stabilized splitting scheme [14] to the governing equation, we get

$$\begin{array}{c}\hfill \frac{{\varphi }_{mn}^{k+1}-{\varphi }_{mn}^k}{\Delta t}=\Delta \left(2\alpha {\varphi }_{mn}^{k+1}-{ϵ}^2{(\Delta {\varphi }^{k+1})}_{mn}+f\left({\varphi }_{mn}^k\right)\right).\end{array}$$

(61)

Thus, Equation (61) can be transformed into the discrete Fourier space as follows:

$$\begin{array}{c}\hfill \frac{{\widehat{\varphi }}_{pq}^{k+1}-{\widehat{\varphi }}_{pq}^k}{\Delta t}=-\left({\xi }_p^2+{\eta }_q^2\right)\left(2\alpha {\widehat{\varphi }}_{pq}^{k+1}+{ϵ}^2({\xi }_p^2+{\eta }_q^2){\widehat{\varphi }}_{pq}^{k+1}+{\widehat{f}}_{pq}^k\right).\end{array}$$

(62)

Therefore, we obtain the following discrete Fourier transform

$$\begin{array}{c}\hfill {\widehat{\varphi }}_{pq}^{k+1}=\frac{{\widehat{\varphi }}_{pq}^k-({\xi }_p^2+{\eta }_q^2)\Delta t{\widehat{f}}_{pq}^k}{1+\Delta t[2\alpha ({\xi }_p^2+{\eta }_q^2)+{ϵ}^2{({\xi }_p^2+{\eta }_q^2)}^2]}.\end{array}$$

(63)

Then, the updated numerical solution ${\varphi }_{mn}^{k+1}$ can be computed using Equation (9):

$$\begin{array}{c}\hfill {\varphi }_{mn}^{k+1}=\frac{1}{N_x{N}_y}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}{\widehat{\varphi }}_{pq}^{k+1}{e}^{i({\xi }_p{x}_m+{\eta }_q{y}_n)}.\end{array}$$

(64)

We perform the numerical test to the CH equation with the following initial condition (39) on the computational domain $\Omega =[-0.5,0.5]\times [-0.5,0.5]$. We use $N_x=N_y=1000$, $h=h_x=h_y$, $ϵ=0.0025$, and $\Delta t=0.00001$. Figure 23a–c shows the numerical results at $t=100\Delta t,300\Delta t,$ and $1000\Delta t$, respectively, and represent well the coarsening dynamics of the CH equation.

According to the numerical simulation results, the unconditionally stable Fourier-spectral method ensures the fast convergence that error is decaying exponentially since the collocation points in physical space guarantee the exact function values in the Fourier space and the method provides a sufficiently smooth basis function (an exponential basis function in this case), one can obtain high-order approximations of derivatives. Moreover, it has almost linearly increasing calculation speed when a large number of grid points are used; hence one can represent the spatial details of numerical results with sufficient grid points in a short time.

4. Numerical Solutions in 3D

We extend the Fourier-spectral method on two-dimensional space to three-dimensional space, $\Omega =[l_x,r_x]\times [l_y,r_y]\times [l_z,r_z]$ for the stated phase-field models. Let $N_x$, $N_y$, $N_z$ be positive even integers and $L_x=r_x-l_x$, $L_y=r_y-l_y$, $L_z=r_z-l_z$ be the length of each direction, respectively. We denote discretized points by $(x_m,y_n,z_o)=(l_x+m{h}_x,l_y+n{h}_y,l_z+o{h}_z)$ for $0\le m\le N_x$, $0\le n\le N_y$, $0\le o\le N_z$, where $h_x=L_x/N_x$, $h_y=L_y/N_y$, $h_z=L_z/N_z$ is the spatial step size in each direction, respectively. For $t_k=k\Delta t$, $\varphi (x_m,y_n,z_o,t_k)$ is denoted by ${\varphi }_{mno}^k$, where $\Delta t$ is the temporal step. For the given data $\left\{{\varphi }_{mno}^k\right|m=1,\dots ,N_x$, $n=1,\dots ,N_y$, $o=1,\dots ,N_z\}$, the discrete Fourier transform is defined as

$$\begin{array}{c}\hfill \begin{array}{c}{\widehat{\varphi }}_{pqr}^k=\sum _{m=1}^{N_x}\sum _{n=1}^{N_y}\sum _{o=1}^{N_z}{\varphi }_{mno}^k{e}^{-i({\xi }_p{x}_m+{\eta }_q{y}_n+{\omega }_r{z}_o)},\hfill \\ -\frac{N_x}{2}+1\le p\le \frac{N_x}{2},-\frac{N_y}{2}+1\le q\le \frac{N_y}{2},-\frac{N_z}{2}+1\le r\le \frac{N_z}{2},\hfill \end{array}\end{array}$$

(65)

where ${\xi }_p=2\pi p/L_x$, ${\eta }_q=2\pi q/L_y$, ${\omega }_r=2\pi r/L_z$. The inverse discrete Fourier transform is

$$\begin{array}{c}\hfill {\varphi }_{mno}^k=N_x{N}_y{N}_z\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\sum _{r=-N_z/2+1}^{N_z/2}{\widehat{\varphi }}_{pqr}^k{e}^{i({\xi }_p{x}_m+{\eta }_q{y}_n+{\omega }_r{z}_o)}.\end{array}$$

(66)

Similar to the two-dimensional case, we assume that $\varphi (x,y,z,t)$ is sufficiently smooth and extended to continuous version of the numerical approximation ${\varphi }_{mno}^k$. Therefore, we can obtain the following result.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\partial }{\partial x}\varphi (x,y,z,t)& =\frac{1}{N_x{N}_y{N}_z}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\sum _{r=-N_z/2+1}^{N_z/2}\left(i{\xi }_p\right)\widehat{\varphi }({\xi }_p,{\eta }_q,{\omega }_r,t)e^{i({\xi }_{p}x+{\eta }_{q}y+{\omega }_{r}z)},\hfill \\ \hfill \frac{\partial }{\partial y}\varphi (x,y,z,t)& =\frac{1}{N_x{N}_y{N}_z}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\sum _{r=-N_z/2+1}^{N_z/2}\left(i{\eta }_q\right)\widehat{\varphi }({\xi }_p,{\eta }_q,{\omega }_r,t)e^{i({\xi }_{p}x+{\eta }_{q}y+{\omega }_{r}z)},\hfill \\ \hfill \frac{\partial }{\partial z}\varphi (x,y,z,t)& =\frac{1}{N_x{N}_y{N}_z}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\sum _{r=-N_z/2+1}^{N_z/2}\left(i{\omega }_r\right)\widehat{\varphi }({\xi }_p,{\eta }_q,{\omega }_r,t)e^{i({\xi }_{p}x+{\eta }_{q}y+{\omega }_{r}z)}.\hfill \end{array}\end{array}$$

(67)

Consequently, we can write Laplacian as coefficients in the Fourier space as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \Delta \varphi (x,y,z,t)& =\frac{1}{N_x{N}_y{N}_z}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\sum _{r=-N_z/2+1}^{N_z/2}\widehat{\Delta \varphi }({\xi }_p,{\eta }_q,{\omega }_r,t)e^{i({\xi }_{p}x+{\eta }_{q}y+{\omega }_{r}z)}\hfill \\ & =\frac{1}{N_x{N}_y{N}_z}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\sum _{r=-N_z/2+1}^{N_z/2}-({\xi }_p^2+{\eta }_q^2+{\omega }_r^2)\widehat{\varphi }({\xi }_p,{\eta }_q,{\omega }_r,t)e^{i({\xi }_{p}x+{\eta }_{q}y+{\omega }_{r}z)},\hfill \end{array}\end{array}$$

(68)

where the first line is the definition of the inverse Fourier transform and the second line is just applying Equation (67) twice to x-, y-, and z-direction to $\varphi (x,y,z,t)$ and its Fourier transform. Now, we present the numerical solutions of three-dimensional phase-field models. Since most of the calculations are redundant, we just briefly list as follows by extending the solvers in two-dimensional space to those of three-dimensional space. Note that the functions f and g which are defined in Section 2 are simply extended to three-dimensional domain. The following is a numerical solution to the AC Equation (2).

$$\begin{array}{c}\hfill {\widehat{\varphi }}_{pqr}^{k+1}=\frac{{ϵ}^2{\widehat{\varphi }}_{pqr}^k-\Delta t{\widehat{f}}_{pqr}^k}{{ϵ}^2+\Delta t[2+{ϵ}^2({\xi }_p^2+{\eta }_q^2+{\omega }_r^2)]}.\end{array}$$

(69)

Then, the renewed numerical solution ${\varphi }_{mno}^{k+1}$ can be computed using Equation (66):

$$\begin{array}{c}\hfill {\varphi }_{mno}^{k+1}=\frac{1}{N_x{N}_y{N}_z}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\sum _{r=-N_z/2+1}^{N_z/2}{\widehat{\varphi }}_{pqr}^{k+1}{e}^{i({\xi }_p{x}_m+{\eta }_q{y}_n+{\omega }_r{z}_o)}.\end{array}$$

(70)

Next one is a numerical solution to the CH Equation (3).

$$\begin{array}{c}\hfill {\widehat{\varphi }}_{pqr}^{k+1}=\frac{{\widehat{\varphi }}_{pqr}^k-({\xi }_p^2+{\eta }_q^2+{\omega }_r^2)\Delta t{\widehat{f}}_{pqr}^k}{1+\Delta t[2({\xi }_p^2+{\eta }_q^2+{\omega }_r^2)+{ϵ}^2{({\xi }_p^2+{\eta }_q^2+{\omega }_r^2)}^2]}.\end{array}$$

(71)

Then, the renewed numerical solution ${\varphi }_{mno}^{k+1}$ can be calculated using Equation (66):

$$\begin{array}{c}\hfill {\varphi }_{mno}^{k+1}=\frac{1}{N_x{N}_y{N}_z}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\sum _{r=-N_z/2+1}^{N_z/2}{\widehat{\varphi }}_{pqr}^{k+1}{e}^{i({\xi }_p{x}_m+{\eta }_q{y}_n+{\omega }_r{z}_o)}.\end{array}$$

(72)

The following is a numerical solution to the SH Equation (5).

$$\begin{array}{c}\hfill {\widehat{\varphi }}_{pqr}^{k+1}=\frac{{\widehat{\varphi }}_{pqr}^k+\Delta t[2({\xi }_p^2+{\eta }_q^2+{\omega }_r^2){\widehat{\varphi }}_{pqr}^k-{\widehat{g}}_{pqr}^k]}{1+\Delta t[(1-ϵ)+{({\xi }_p^2+{\eta }_q^2+{\omega }_r^2)}^2]},\end{array}$$

(73)

and hence we have a numerical solution ${\varphi }_{mno}^{k+1}$ as follows:

$$\begin{array}{c}\hfill {\varphi }_{mno}^{k+1}=\frac{1}{N_x{N}_y{N}_z}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\sum _{r=-N_z/2+1}^{N_z/2}{\widehat{\varphi }}_{pqr}^{k+1}{e}^{i({\xi }_p{x}_m+{\eta }_q{y}_n+{\omega }_r{z}_o)}.\end{array}$$

(74)

Subsequently, we have a numerical solution to the classical PFC model (6) as follows.

$$\begin{array}{c}\hfill {\widehat{\varphi }}_{pqr}^{k+1}=\frac{{\widehat{\varphi }}_{pqr}^k+\Delta t[2{({\xi }_p^2+{\eta }_q^2+{\omega }_r^2)}^2{\widehat{\varphi }}_{pqr}^k-({\xi }_p^2+{\eta }_q^2+{\omega }_r^2){\widehat{g}}_{pqr}^k]}{1+\Delta t[(1-ϵ)({\xi }_p^2+{\eta }_q^2+{\omega }_r^2)+{({\xi }_p^2+{\eta }_q^2+{\omega }_r^2)}^3]}.\end{array}$$

(75)

Therefore, we have a numerical solution ${\varphi }_{mno}^{k+1}$ as follows:

$$\begin{array}{c}\hfill {\varphi }_{mno}^{k+1}=\frac{1}{N_x{N}_y{N}_z}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\sum _{r=-N_z/2+1}^{N_z/2}{\widehat{\varphi }}_{pqr}^{k+1}{e}^{i({\xi }_p{x}_m+{\eta }_q{y}_n+{\omega }_r{z}_o)}.\end{array}$$

(76)

5. Numerical Experiments in 3D

5.1. The AC Equation

We perform the numerical test to the AC Equation (2) with the following initial condition,

$$\begin{array}{ccc}\hfill \varphi (x,y,z,0)& =& tanh\left(\frac{1.3-\sqrt{x^2+y^2+z^2}}{\sqrt{2}ϵ}\right),\hfill \end{array}$$

(77)

Figure 24 shows the numerical solutions at $t=0,480\Delta t,800\Delta t$ to the 3D AC equation using the Fourier-spectral method for $(x,y,z)\in [-1.5,1.5]\times [-1.5,1.5]\times [-1.5,1.5]$. Here, we use $N_x=N_y=N_z=128$, $h=3/128$, $ϵ=2h$, and $\Delta t=0.001$. The final time is $T=0.8$. Similar to the two-dimensional case, the simulation results depict the motion by mean curvature property quite well.

We performed the numerical experiment for the complex shape to the 3D AC equation. We defined two initial conditions of complex shape as the following:

$$\begin{array}{ccc}\hfill \varphi (x,y,z,0)& =& tanh\left(\frac{1+0.2cos\left(6\theta \right)-\sqrt{x^2+2y^2+z^2}}{\sqrt{2}ϵ}\right),\hfill \end{array}$$

(78)

$$\begin{array}{ccc}\hfill \varphi (x,y,z,0)& =& tanh\left(\frac{{(x^2+9y^2/4+z^2-1)}^3-x^2{z}^3-9y^2{z}^3/80}{\sqrt{2}ϵ}\right),\hfill \end{array}$$

(79)

where $\theta ={tan}^{-1}(z/x)$ for $(x,y,z)\in [-1.5,1.5]\times [-1.5,1.5]\times [-1.5,1.5]$. Figure 25 shows the numerical test results at $t=0,500\Delta t,1000\Delta t$. Here, we use $N_x=N_y=N_z=128$, $h=3/128$, $ϵ=0.05$, and $\Delta t=0.0001$.

Furthermore, we present the energy dissipation to the AC equation in three-dimensional space. The discrete energy has a form as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{E}}^h\left({\varphi }^k\right)=& \sum _{m=1}^{N_x}\sum _{n=1}^{N_y}\sum _{o=1}^{N_z}\frac{F\left({\varphi }_{mno}^k\right)}{{ϵ}^2}{h}^3\hfill \\ & +\frac{h^3}{N_x{N}_y{N}_z}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\sum _{r=-N_z/2+1}^{N_z/2}\frac{{\left|{\xi }_p{\widehat{\varphi }}_{pqr}^k\right|}^2+{\left|{\eta }_q{\widehat{\varphi }}_{pqr}^k\right|}^2+{\left|{\omega }_r{\widehat{\varphi }}_{pqr}^k\right|}^2}{2}.\hfill \end{array}\end{array}$$

(80)

Figure 26 illustrates the non-increase discrete energy Equation (80) over time. Similar settings are given as above except for $T=0.6$.

5.2. The CH Equation

Next, we investigated the numerical test to the CH Equation (3) with the following initial conditions,

$$\begin{array}{ccc}\hfill \varphi (x,y,z,0)& =& 0.05\mathrm{rand}(x,y,z),\hfill \end{array}$$

(81)

$$\begin{array}{ccc}\hfill \varphi (x,y,z,0)& =& -0.45+0.05\mathrm{rand}(x,y,z),\hfill \end{array}$$

(82)

where rand implies a random number generating function ranged from $-1$ to 1 and $(x,y,z)\in [-0.5,0.5]\times [-0.5,0.5]\times [-0.5,0.5]$. Figure 27 shows the numerical investigation results at $t=100\Delta t,500\Delta t$, and $1000\Delta t$ to the 3D CH equation using the Fourier-spectral method with the initial condition (81). Here, we use $N_x=N_y=N_z=128$, $h=1/128$, $ϵ=h$, and $\Delta t=0.001$.

Figure 28 depicts the numerical simulation results to the 3D CH equation using the Fourier-spectral method with the initial condition (82). Here, we use $N_x=N_y=N_z=128$, $h=1/128$, $ϵ=h$, and $\Delta t=0.001$. The above results represent well the coarsening dynamics of the CH equation. Moreover, we present the energy dissipation to the CH equation in three-dimensional space. We extend Equation (42) to three-dimensional case as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{E}}^h\left({\varphi }^k\right)=& \sum _{m=1}^{N_x}\sum _{n=1}^{N_y}\sum _{o=1}^{N_z}F\left({\varphi }_{mno}^k\right)h^3\hfill \\ & +\frac{h^3}{N_x{N}_y{N}_z}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\sum _{r=-N_z/2+1}^{N_z/2}\frac{{ϵ}^2}{2}({\left|{\xi }_p{\widehat{\varphi }}_{pqr}^k\right|}^2+{\left|{\eta }_q{\widehat{\varphi }}_{pqr}^k\right|}^2+{\left|{\omega }_r{\widehat{\varphi }}_{pqr}^k\right|}^2).\hfill \end{array}\end{array}$$

(83)

Figure 29 describes the non-increase energy dissipation over time with the initial condition (81).

Further we present phase separation behaviors of the nonlocal CH equation in the domain $[0,1]\times [0,1]\times [0,1]$. An initial condition is given by

$$\begin{array}{c}\hfill \varphi (x,y,z,0)=\overline{\varphi }+0.15\mathrm{rand}(x,y,z).\end{array}$$

(84)

Figure 30 shows the snapshots of the spherical pattern formation via nonlocal CH equation. Here, we use Equation (84) with $N_x=N_y=N_z=80$, $h=1/80$, $ϵ=2h$, $\sigma =200$, $\overline{\varphi }=-0.3$, and $\Delta t=0.1$. The final time is $T=1000$.

Similarly, we conduct numerical simulations to generate the lamellar patterns. Figure 31 depicts the snapshots of the lamellar pattern formation via nonlocal CH equation with the initial condition (84). We take $N_x=N_y=N_z=80$, $h=1.2/80$, $ϵ=2.5h$, $\sigma =100$, $\overline{\varphi }=0$, and $\Delta t=0.1$. The final time is $T=1000$.

We present a result of long time energy decay in a numerical way. Figure 32 shows the energy decay of the nonlocal CH equation, where the discretization is simply extended from Equation (46) as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{E}}^h\left({\varphi }^k\right)=& \frac{h^3}{N_x{N}_y{N}_z}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\sum _{r=-N_z/2+1}^{N_z/2}[\frac{{ϵ}^2}{2}({\left|{\xi }_p{\widehat{\varphi }}_{pqr}^k\right|}^2+{\left|{\eta }_q{\widehat{\varphi }}_{pqr}^k\right|}^2+{\left|{\omega }_r{\widehat{\varphi }}_{pqr}^k\right|}^2)\hfill \\ & +\frac{\sigma }{2}({\left|{\xi }_p{\widehat{\psi }}_{pqr}^k\right|}^2+{\left|{\eta }_q{\widehat{\psi }}_{pqr}^k\right|}^2+{\left|{\omega }_r{\widehat{\psi }}_{pqr}^k\right|}^2)]+\sum _{m=1}^{N_x}\sum _{n=1}^{N_y}\sum _{o=1}^{N_z}F\left({\varphi }_{mno}^k\right)h^3.\hfill \end{array}\end{array}$$

(85)

5.3. The SH Equation

We present numerical simulation results to the SH equation in this section. A computational domain is $[0,100]\times [0,100]\times [0,100]$. An initial condition is given as

$$\begin{array}{c}\hfill \varphi (x,y,z,0)=\overline{\varphi }+\overline{\varphi }\mathrm{rand}(x,y,z),\end{array}$$

(86)

where $\overline{\varphi }=0.1$. Figure 33 shows the patterns due to the instability of thermal fluctuations. Here, we use $N_x=N_y=N_z=120$, $h=100/120$, $ϵ=0.25$, and $\Delta t=0.1$. The final time is $T=2400$.

As in the case of two-dimensional space, the reduction of the following discrete energy Equation (87) is represented in Figure 34.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{E}}^h\left({\varphi }^k\right)=& \frac{h^3}{N_x{N}_y{N}_z}\sum _{p=-N_x/2+1}^{N_x/2}\sum _{q=-N_y/2+1}^{N_y/2}\sum _{r=-N_z/2+1}^{N_z/2}[-({\left|{\xi }_p{\widehat{\varphi }}_{pqr}^k\right|}^2+{\left|{\eta }_q{\widehat{\varphi }}_{pqr}^k\right|}^2+{\left|{\omega }_r{\widehat{\varphi }}_{pqr}^k\right|}^2)\hfill \\ & +\frac{1}{2}{\left|({\xi }_p^2+{\eta }_q^2+{\omega }_r^2){\widehat{\varphi }}_{pqr}^k\right|}^2]+\sum _{m=1}^{N_x}\sum _{n=1}^{N_y}\sum _{o=1}^{N_z}\left(\frac{1}{4}{\left({\varphi }_{mno}^k\right)}^4+\frac{(1-ϵ)}{2}{\left({\varphi }_{mno}^k\right)}^2\right).\hfill \end{array}\end{array}$$

(87)

5.4. Classical PFC Model

We present numerical simulation results to the classical PFC model in the three-dimensional domain $[0,80]\times [0,80]\times [0,80]$. An initial condition is given as

$$\begin{array}{c}\hfill \varphi (x,y,z,0)=\overline{\varphi }+\overline{\varphi }\mathrm{rand}(x,y,z),\end{array}$$

(88)

where $\overline{\varphi }=0.07$. Figure 35 shows the phase transition behaviors of this PFC model. Here, we use $N_x=N_y=N_z=60$, $h=80/60$, $ϵ=0.025$, and $\Delta t=0.01$. The final time is $T=1500$.

We present the energy dissipation in three-dimensional case. Figure 36 shows the energy decay over time, Equation (87), via classical PFC model.

6. Conclusions

In this paper, we briefly investigate the Fourier-spectral method to the phase-field equations in two- and three-dimensional cases, especially focused on AC, CH, SH, PFC, MBE growth models which are very important parabolic partial differential equations and are applicable to a lot of interesting scientific problems. We present detailed descriptions and numerical solutions of the unconditionally gradient stable Fourier-spectral method and explain its connection to MATLAB usage so that the interested readers can simply implement the corresponding Fourier-spectral method for their research needs without difficulties. Numerical simulation results yield the corresponding method represents well the characteristics of each equation. Furthermore, we provide the MATLAB codes for these equations in the Appendix A.