# Benchmark Problems for the Numerical Schemes of the Phase-Field Equations

Authors: Youngjin Hwang; Chaeyoung Lee; Soobin Kwak; Yongho Choi; Seokjun Ham; Seungyoon Kang; Junxiang Yang; Junseok Kim
Journal: Discrete Dynamics in Nature and Society
 (2022)

---

## Abstract

In this study, we present benchmark problems for the numerical methods of the phase-field equations. To find appropriate benchmark problems, we first perform a linear stability analysis and then take a growth mode solution as the benchmark problem, which is closely related to the dynamics of the original governing equations. As concrete examples, we perform convergence tests of the numerical methods of the Allen–Cahn (AC) and Cahn–Hilliard (CH) equations using the proposed benchmark problems. The one- and two-dimensional computational experiments confirm the accuracy and efficiency of the proposed scheme as the benchmark problems.


---

## Body

## 1. Introduction

In this study, we present benchmark problems for the numerical schemes of the phase-field equations. Two famous phase-field models are chosen as examples. The first equation is the Allen–Cahn (AC) equation [1, 2]:(1)ϕtx,t=−F′ϕx,tε2+Δϕx,t,forx∈Ω,t>0.with the Neumann boundary conditionn⋅∇ϕ=0 on ∂Ω, where n is the exterior normal vector to the domain boundary ∂Ω. The phase-field ϕx,t is the difference between the concentrations of the two mixtures components, the free energy Fϕ=0.25ϕ2−12 is double-well potential, and ϵ is a positive constant related to the thickness of the interfacial transition layer. The second equation is the Cahn–Hilliard (CH) equation [3]:(2)ϕtx,t=ΔF′ϕx,t−ϵ2Δϕx,t,forx∈Ω,t>0.with the Neumann boundary conditionn⋅∇ϕ=n⋅∇⋅Δϕ=0 on ∂Ω. Because there are no closed-form analytic solutions for the general initial and boundary conditions of the AC and CH equations, it is required to approximate the solutions of the equations using numerical methods. Therefore, to validate the accuracy of the newly proposed numerical methods, it is essential to test the developed methods with benchmark problems. The logarithmic Flory–Huggins potential is one of the useful potentials among phase-field models. In [4], the authors developed an appropriate temporal discretization method for the nonlinear term of the fourth-order CH equation with concentration dependent mobility and logarithmic Flory–Huggins potentials. The developed Invariant Energy Quadratization (IEQ) method is linear and unconditionally stable. Shen et al. [5] proposed a scalar auxiliary variable (SAV) approach based on the IEQ method to construct accurate and efficient time discretization methods for a large class of gradient flows. Yang et al. [6] developed an improved SAV approach inspired by a step-by-step solving scheme based on the SAV approach for gradient flow problems which includes the AC and CH equations [7]. The improved SAV approach has all the advantages of classical SAV approach, further simplifies the algorithm, and easily constructs the temporally first-order and second-order methods. To solve the AC equation with Flory–Huggins potential, a novel energy factorization approach is used based on a stabilization technique called the stabilized energy factorization approach [8]. The discrete maximum principle and unconditional energy stability of the novel energy factorization approach have been rigorously proved using the discrete variational principle. A novel linear, energy stable, and maximum principle preserving scheme was developed to semi-implicitly treat double-well potentials in the AC equation using the energy factorization approach [9]. A variety of benchmark problems have been proposed for the phase-field models and used to validate the numerical schemes. Jokisaari et al. [10] proposed two benchmark problems: solute diffusion and second-phase growth and coarsening. Jeong et al. [11] presented two benchmark problems, which are shrinking annulus and spherical shell. They considered the CH equation in radially and spherically symmetric forms to obtain simple benchmark solutions. For the AC and CH equations, Church et al. [12] provided four benchmark problems. Zhang et al. [13] tested a class of linear numerical schemes for the functionalized CH equation using stabilized scalar auxiliary variable (SAV) method. Wu et al. [14] presented two benchmark problems on homogeneous and heterogeneous nucleation. Li et al. [15] suggested the numerical benchmark solution of the CH equation. They adopted the fourth-order Runge–Kutta method and finite difference method for the integration in time and the spatial differential operator, respectively, with a cosine initial condition. Hug et al. [16] proposed a benchmark problem for brittle fracture.For example, Jeong et al. [11] used the following initial state for a benchmark problem:(3)ϕr,0=tanh0.1−∣r−0.75∣2ε.In [12], the authors used the following four benchmark problems:(4)ϕx,0=cos2x+0.01ecosx+0.1,(5)ϕx,y,0=0.05cos3xcos4y+cos24xcos23y+cosx−5ycos2x−y,(6)ϕx,y,0=tanhx−π2+y−π2−22ε,(7)ϕx,y,0=−1+∑i=17fx−xi2+y−yi2−ri,where fs=2e−ϵ2/s2 if s<0 and fs=0 otherwise. In [13], the authors used the following initial condition for the benchmark problem:(8)ϕx,y,0=2esinx+siny−2+2.2e−sinx−siny−2−1.Figure1 shows some of the initial conditions for the benchmark problems.Figure 1
(a–f) The different initial conditions for benchmark problems from equations (3) to (8), respectively.
(a)(b)(c)(d)(e)(f)However, most of the previous benchmark problems for the phase-field equations are chosen without any concrete theoretical basis. In [17], the authors developed a mathematical model including heat equation and optimization of the experimental parameters and temperature profiles, and the mathematical model was validated according to the numerical results without suggesting a benchmark problem. Benchmark problems with concrete theoretical basis are an important basis for validating the accuracy of mathematical models and numerical methods.In this paper, we propose appropriate benchmark problems to verify the accuracy of the numerical methods based on the theoretical basis through linear stability analysis using simple initial conditions. First, we perform a linear stability analysis and then take a growth mode solution [18] as the benchmark problem, which is closely related to the dynamics of the original governing equations.The layout of this paper is as follows. In Section2, the procedure of finding benchmark problems is described. In Section 3, the numerical experiments are performed. In Section 4, the conclusions are given.


## 2. Proposed Benchmark Problems

In this section, the proposed benchmark problems are presented for the numerical methods of the one- and two-dimensional AC and CH phase-field models. We first conduct a linear stability analysis and then take a growth mode solution as the benchmark problem, which is closely related to the dynamics of the AC equation or the CH equation.


### 2.1. Benchmark Problems for the AC Equation

First, let us consider the benchmark problems for the one- and two-dimensional AC equation.


#### 2.1.1. One-Dimensional AC Equation

The one-dimensional AC equation is as follows:(9)ϕtx,t=−ϕ3x,t−ϕx,tε2+ϕxxx,t,forx∈Ω=0,2π,t>0.with the Neumann boundary conditionϕx0,t=ϕx2π,t=0. To find appropriate benchmark problems for the AC equation, a linear stability analysis is first conducted around a spatially constant critical composition solution ϕ=0. We linearize the nonlinear term F′ϕ=ϕ3−ϕ by applying the Taylor expansion and then get linearized term F′ϕ≈−ϕ. Therefore, the linearized AC equation is as follows:(10)ϕtx,t=ϕx,tε2+ϕxxx,t.For a positive integerk, we take Φx,t=αtcoskx as a solution for (10), where αt is an amplitude. Substituting above Φx,t into (10), we have(11)α′tcoskx=αtcoskxε2−k2αtcoskx.Dividingcoskx on the both sides of (11), we obtain(12)α′t=1ε2−k2αt.Then, the solution of the ordinary differential equation (ODE) (12) is as follows:(13)αt=α0exp1ε2−k2t.Let(14)Φx,t=αtcoskx.Be a benchmark problem solution for a one-dimensional modified AC equation, whereαt is given by (13). Here, Φx,0=α0coskx is the initial condition. Next, when performing the convergence tests for the numerical schemes of the AC equation, we consider the following modified AC equation with a source term:(15)ϕtx,t=−F′ϕx,tϵ2+ϕxxx,t+sx,tfort>0,where(16)sx,t=Φtx,t+F′Φx,tϵ2−Φxxx,t=α′t+1ϵ2α3tcos2kx−αt+k2αtcoskx.Here, the value ofαt is given in (13) and α′t can be found by taking the first derivative of αt with respect to t. Note that if k is chosen to satisfy k<1/ε, then the numerical solution grows as time evolves as seen from (13).


#### 2.1.2. Two-Dimensional AC Equation

The two-dimensional AC equation is as follows:(17)ϕtx,y,t=−ϕ3x,y,t−ϕx,y,tε2+Δϕx,y,t,forx,y∈Ω=0,2π×0,2π,t>0.with the Neumann boundary conditionn⋅∇ϕ=0 on ∂Ω, where n is the exterior normal vector to the domain boundary ∂Ω. To find appropriate benchmark problems for the two-dimensional AC equation, we first conduct a linear stability analysis around a spatially constant critical composition solution ϕ=0. The nonlinear term F′ϕ=ϕ3−ϕ is linearized by applying the Taylor expansion, and then the linearized term F′ϕ≈−ϕ is obtained. Therefore, the linearized AC equation is as follows:(18)ϕtx,y,t=ϕx,y,tϵ2+ϕxxx,y,t+ϕyyx,y,t.For positive integerskx and ky, let us consider Φx,y,t=αtcoskxxcoskyy as a two-dimensional benchmark solution for (18), where αt is an amplitude. Substituting above Φx,y,t into (18), then(19)α′tcoskxxcoskyy=αtcoskxxcoskyyϵ2−kx2+ky2αtcoskxxcoskyy.Dividingcoskxxcoskyy on the both sides of (19), we obtain(20)α′t=1ϵ2−kx2−ky2αt.Then, the solution of the ODE (19) is as follows:(21)αt=α0exp1ε2−kx2−ky2t.Let(22)Φx,y,t=αtcoskxxcoskyy.be a benchmark problem solution for the modified AC equation, whereαt is given by (21). The initial condition is given by Φx,y,0=α0coskxxcoskyy. Finally, when we perform convergence tests for numerical schemes of the AC equation, the following modified AC equation with a source term is taken:(23)ϕtx,y,t=−F′ϕx,y,tε2+Δϕx,y,t+sx,y,t,forx,y∈Ω,t>0,where(24)sx,y,t=Φtx,y,t+F′Φx,y,tε2−ΔΦx,y,t=α′t+1ε2α3tcos2kxxcos2kyy−αt+kx2+ky2αtcoskxxcoskyy.Here, the value ofαt is given in (21) and α′t can be found by taking the first derivative of αt with respect to t. Note that if we set kx and ky satisfying kx2+ky2<1/ε2, then the numerical solution grows as time evolves as seen from (21).


### 2.2. Benchmark Problems for the CH Equation

Next, let us consider the benchmark problems for the one- and two-dimensional CH equation.


#### 2.2.1. One-Dimensional CH Equation

The one-dimensional CH equation is as follows:(25)ϕtx,t=ϕ3x,t−ϕx,t−ϵ2ϕxxx,txx,forx∈Ω=0,2π,t>0.with the Neumann boundary conditionϕx0,t=ϕx2π,t=ϕxxx0,t=ϕxxx2π,t=0. Similar to the benchmark problem solution of the AC equation, a linear stability analysis is first conducted around a spatially constant critical composition solution ϕ=0. Then, the linearized CH equation is as follows:(26)ϕtx,t=−ϕxxx,t−ϵ2ϕxxxxx,t.For the positive integerk, we consider Φx,t=βtcoskx as a solution for (26), where βt is an amplitude. Substituting above Φx,t into (26), we have(27)β′tcoskx=k2βtcoskx−k4ϵ2βtcoskx.Dividingcoskx on the both sides of (27), we obtain(28)β′t=k21−kϵ2βt.Then, the solution of the ODE (28) is as follows:(29)βt=β0expk21−kϵ2t.Let(30)Φx,t=βtcoskx.be a benchmark problem solution for the CH equation, whereβt is given by (29). The initial condition is given by Φx,0=β0coskx. Finally, to conduct convergence tests for numerical schemes of the CH equation, the following modified CH equation with a source term is considered:(31)ϕtx,t=F′ϕx,t−ϵ2ϕxxx,txx+sx,tfort>0,where(32)sx,t=Φtx,t−F′Φx,t−ϵ2Φxxx,txx=β′t+3k2β3t4−k2βt+k4ϵ2βtcoskx+9k2β3t4cos3kx,where the triple angle identity for cosine is applied, that is, cos3kx=cos3kx+3coskx/4. Here, the value of βt is given in (29) and β′t can be found by taking the first derivative of βt with respect to t. Note that if k is chosen to satisfy k<1/ε, then the numerical solution grows as time evolves as seen from (29).


#### 2.2.2. Two-Dimensional CH Equation

The two-dimensional CH equation is as follows:(33)ϕtx,y,t=Δϕ3x,y,t−ϕx,y,t−ε2Δϕx,y,t,forx,y∈Ω=0,2π×0,2π,t>0.With the Neumann boundary conditionn⋅∇ϕ=0 on ∂Ω, where n is the exterior normal vector to the domain boundary ∂Ω. A linear stability analysis is conducted around a spatially constant critical composition solution ϕ=0. Then, the linearized CH equation is as follows:(34)ϕtx,y,t=Δ−ϕx,y,t−ϵ2Δϕx,y,t.For the positive integerskx and ky, Φx,y,t=βtcoskxxcoskyy is given as a solution for (34), where βt is an amplitude. Substituting above Φx,y,t into (34), we have(35)β′tcoskxxcoskyy=kx2+ky21−kx2+ky2ϵ2βtcoskxxcoskyy.Dividingcoskxxcoskyy on the both sides of (35), we obtain(36)β′t=kx2+ky21−kx2+ky2ϵ2βt.Then, the solution of the ODE (36) is as follows:(37)βt=β0expkx2+ky21−kx2+ky2ϵ2t.Let(38)Φx,y,t=βtcoskxxcoskyy,be a benchmark problem solution for the CH equation, where βt is given by (37). The initial condition is given by Φx,y,0=β0coskxxcoskyy. Finally, when conducting convergence tests for the numerical schemes of the CH equation, we consider the following modified CH equation with a source term:(39)ϕtx,y,t=ΔF′ϕx,y,t−ϵ2Δϕx,y,t+sx,y,t,where(40)sx,y,t=Φtx,y,t−ΔF′Φx,y,t−ε2ΔΦx,y,t=β3t163coskxxcos3kyykx2+9ky2+3cos3kxxcoskyy9kx2+ky2+9coskxxcoskyykx2+ky2+9cos3kxxcos3kyykx2+ky2+β′t−βt1−ε2kx2+ky2kx2+ky2coskxxcoskyy.Here, the value ofβt is given in (37) and β′t can be found by taking the first derivative of βt with respect to t. Note that if kx and ky are given satisfying kx2+ky2<1/ε2, then the numerical solution grows as time evolves as seen from (37).


## 2.1. Benchmark Problems for the AC Equation

First, let us consider the benchmark problems for the one- and two-dimensional AC equation.


### 2.1.1. One-Dimensional AC Equation

The one-dimensional AC equation is as follows:(9)ϕtx,t=−ϕ3x,t−ϕx,tε2+ϕxxx,t,forx∈Ω=0,2π,t>0.with the Neumann boundary conditionϕx0,t=ϕx2π,t=0. To find appropriate benchmark problems for the AC equation, a linear stability analysis is first conducted around a spatially constant critical composition solution ϕ=0. We linearize the nonlinear term F′ϕ=ϕ3−ϕ by applying the Taylor expansion and then get linearized term F′ϕ≈−ϕ. Therefore, the linearized AC equation is as follows:(10)ϕtx,t=ϕx,tε2+ϕxxx,t.For a positive integerk, we take Φx,t=αtcoskx as a solution for (10), where αt is an amplitude. Substituting above Φx,t into (10), we have(11)α′tcoskx=αtcoskxε2−k2αtcoskx.Dividingcoskx on the both sides of (11), we obtain(12)α′t=1ε2−k2αt.Then, the solution of the ordinary differential equation (ODE) (12) is as follows:(13)αt=α0exp1ε2−k2t.Let(14)Φx,t=αtcoskx.Be a benchmark problem solution for a one-dimensional modified AC equation, whereαt is given by (13). Here, Φx,0=α0coskx is the initial condition. Next, when performing the convergence tests for the numerical schemes of the AC equation, we consider the following modified AC equation with a source term:(15)ϕtx,t=−F′ϕx,tϵ2+ϕxxx,t+sx,tfort>0,where(16)sx,t=Φtx,t+F′Φx,tϵ2−Φxxx,t=α′t+1ϵ2α3tcos2kx−αt+k2αtcoskx.Here, the value ofαt is given in (13) and α′t can be found by taking the first derivative of αt with respect to t. Note that if k is chosen to satisfy k<1/ε, then the numerical solution grows as time evolves as seen from (13).


### 2.1.2. Two-Dimensional AC Equation

The two-dimensional AC equation is as follows:(17)ϕtx,y,t=−ϕ3x,y,t−ϕx,y,tε2+Δϕx,y,t,forx,y∈Ω=0,2π×0,2π,t>0.with the Neumann boundary conditionn⋅∇ϕ=0 on ∂Ω, where n is the exterior normal vector to the domain boundary ∂Ω. To find appropriate benchmark problems for the two-dimensional AC equation, we first conduct a linear stability analysis around a spatially constant critical composition solution ϕ=0. The nonlinear term F′ϕ=ϕ3−ϕ is linearized by applying the Taylor expansion, and then the linearized term F′ϕ≈−ϕ is obtained. Therefore, the linearized AC equation is as follows:(18)ϕtx,y,t=ϕx,y,tϵ2+ϕxxx,y,t+ϕyyx,y,t.For positive integerskx and ky, let us consider Φx,y,t=αtcoskxxcoskyy as a two-dimensional benchmark solution for (18), where αt is an amplitude. Substituting above Φx,y,t into (18), then(19)α′tcoskxxcoskyy=αtcoskxxcoskyyϵ2−kx2+ky2αtcoskxxcoskyy.Dividingcoskxxcoskyy on the both sides of (19), we obtain(20)α′t=1ϵ2−kx2−ky2αt.Then, the solution of the ODE (19) is as follows:(21)αt=α0exp1ε2−kx2−ky2t.Let(22)Φx,y,t=αtcoskxxcoskyy.be a benchmark problem solution for the modified AC equation, whereαt is given by (21). The initial condition is given by Φx,y,0=α0coskxxcoskyy. Finally, when we perform convergence tests for numerical schemes of the AC equation, the following modified AC equation with a source term is taken:(23)ϕtx,y,t=−F′ϕx,y,tε2+Δϕx,y,t+sx,y,t,forx,y∈Ω,t>0,where(24)sx,y,t=Φtx,y,t+F′Φx,y,tε2−ΔΦx,y,t=α′t+1ε2α3tcos2kxxcos2kyy−αt+kx2+ky2αtcoskxxcoskyy.Here, the value ofαt is given in (21) and α′t can be found by taking the first derivative of αt with respect to t. Note that if we set kx and ky satisfying kx2+ky2<1/ε2, then the numerical solution grows as time evolves as seen from (21).


## 2.1.1. One-Dimensional AC Equation

The one-dimensional AC equation is as follows:(9)ϕtx,t=−ϕ3x,t−ϕx,tε2+ϕxxx,t,forx∈Ω=0,2π,t>0.with the Neumann boundary conditionϕx0,t=ϕx2π,t=0. To find appropriate benchmark problems for the AC equation, a linear stability analysis is first conducted around a spatially constant critical composition solution ϕ=0. We linearize the nonlinear term F′ϕ=ϕ3−ϕ by applying the Taylor expansion and then get linearized term F′ϕ≈−ϕ. Therefore, the linearized AC equation is as follows:(10)ϕtx,t=ϕx,tε2+ϕxxx,t.For a positive integerk, we take Φx,t=αtcoskx as a solution for (10), where αt is an amplitude. Substituting above Φx,t into (10), we have(11)α′tcoskx=αtcoskxε2−k2αtcoskx.Dividingcoskx on the both sides of (11), we obtain(12)α′t=1ε2−k2αt.Then, the solution of the ordinary differential equation (ODE) (12) is as follows:(13)αt=α0exp1ε2−k2t.Let(14)Φx,t=αtcoskx.Be a benchmark problem solution for a one-dimensional modified AC equation, whereαt is given by (13). Here, Φx,0=α0coskx is the initial condition. Next, when performing the convergence tests for the numerical schemes of the AC equation, we consider the following modified AC equation with a source term:(15)ϕtx,t=−F′ϕx,tϵ2+ϕxxx,t+sx,tfort>0,where(16)sx,t=Φtx,t+F′Φx,tϵ2−Φxxx,t=α′t+1ϵ2α3tcos2kx−αt+k2αtcoskx.Here, the value ofαt is given in (13) and α′t can be found by taking the first derivative of αt with respect to t. Note that if k is chosen to satisfy k<1/ε, then the numerical solution grows as time evolves as seen from (13).


## 2.1.2. Two-Dimensional AC Equation

The two-dimensional AC equation is as follows:(17)ϕtx,y,t=−ϕ3x,y,t−ϕx,y,tε2+Δϕx,y,t,forx,y∈Ω=0,2π×0,2π,t>0.with the Neumann boundary conditionn⋅∇ϕ=0 on ∂Ω, where n is the exterior normal vector to the domain boundary ∂Ω. To find appropriate benchmark problems for the two-dimensional AC equation, we first conduct a linear stability analysis around a spatially constant critical composition solution ϕ=0. The nonlinear term F′ϕ=ϕ3−ϕ is linearized by applying the Taylor expansion, and then the linearized term F′ϕ≈−ϕ is obtained. Therefore, the linearized AC equation is as follows:(18)ϕtx,y,t=ϕx,y,tϵ2+ϕxxx,y,t+ϕyyx,y,t.For positive integerskx and ky, let us consider Φx,y,t=αtcoskxxcoskyy as a two-dimensional benchmark solution for (18), where αt is an amplitude. Substituting above Φx,y,t into (18), then(19)α′tcoskxxcoskyy=αtcoskxxcoskyyϵ2−kx2+ky2αtcoskxxcoskyy.Dividingcoskxxcoskyy on the both sides of (19), we obtain(20)α′t=1ϵ2−kx2−ky2αt.Then, the solution of the ODE (19) is as follows:(21)αt=α0exp1ε2−kx2−ky2t.Let(22)Φx,y,t=αtcoskxxcoskyy.be a benchmark problem solution for the modified AC equation, whereαt is given by (21). The initial condition is given by Φx,y,0=α0coskxxcoskyy. Finally, when we perform convergence tests for numerical schemes of the AC equation, the following modified AC equation with a source term is taken:(23)ϕtx,y,t=−F′ϕx,y,tε2+Δϕx,y,t+sx,y,t,forx,y∈Ω,t>0,where(24)sx,y,t=Φtx,y,t+F′Φx,y,tε2−ΔΦx,y,t=α′t+1ε2α3tcos2kxxcos2kyy−αt+kx2+ky2αtcoskxxcoskyy.Here, the value ofαt is given in (21) and α′t can be found by taking the first derivative of αt with respect to t. Note that if we set kx and ky satisfying kx2+ky2<1/ε2, then the numerical solution grows as time evolves as seen from (21).


## 2.2. Benchmark Problems for the CH Equation

Next, let us consider the benchmark problems for the one- and two-dimensional CH equation.


### 2.2.1. One-Dimensional CH Equation

The one-dimensional CH equation is as follows:(25)ϕtx,t=ϕ3x,t−ϕx,t−ϵ2ϕxxx,txx,forx∈Ω=0,2π,t>0.with the Neumann boundary conditionϕx0,t=ϕx2π,t=ϕxxx0,t=ϕxxx2π,t=0. Similar to the benchmark problem solution of the AC equation, a linear stability analysis is first conducted around a spatially constant critical composition solution ϕ=0. Then, the linearized CH equation is as follows:(26)ϕtx,t=−ϕxxx,t−ϵ2ϕxxxxx,t.For the positive integerk, we consider Φx,t=βtcoskx as a solution for (26), where βt is an amplitude. Substituting above Φx,t into (26), we have(27)β′tcoskx=k2βtcoskx−k4ϵ2βtcoskx.Dividingcoskx on the both sides of (27), we obtain(28)β′t=k21−kϵ2βt.Then, the solution of the ODE (28) is as follows:(29)βt=β0expk21−kϵ2t.Let(30)Φx,t=βtcoskx.be a benchmark problem solution for the CH equation, whereβt is given by (29). The initial condition is given by Φx,0=β0coskx. Finally, to conduct convergence tests for numerical schemes of the CH equation, the following modified CH equation with a source term is considered:(31)ϕtx,t=F′ϕx,t−ϵ2ϕxxx,txx+sx,tfort>0,where(32)sx,t=Φtx,t−F′Φx,t−ϵ2Φxxx,txx=β′t+3k2β3t4−k2βt+k4ϵ2βtcoskx+9k2β3t4cos3kx,where the triple angle identity for cosine is applied, that is, cos3kx=cos3kx+3coskx/4. Here, the value of βt is given in (29) and β′t can be found by taking the first derivative of βt with respect to t. Note that if k is chosen to satisfy k<1/ε, then the numerical solution grows as time evolves as seen from (29).


### 2.2.2. Two-Dimensional CH Equation

The two-dimensional CH equation is as follows:(33)ϕtx,y,t=Δϕ3x,y,t−ϕx,y,t−ε2Δϕx,y,t,forx,y∈Ω=0,2π×0,2π,t>0.With the Neumann boundary conditionn⋅∇ϕ=0 on ∂Ω, where n is the exterior normal vector to the domain boundary ∂Ω. A linear stability analysis is conducted around a spatially constant critical composition solution ϕ=0. Then, the linearized CH equation is as follows:(34)ϕtx,y,t=Δ−ϕx,y,t−ϵ2Δϕx,y,t.For the positive integerskx and ky, Φx,y,t=βtcoskxxcoskyy is given as a solution for (34), where βt is an amplitude. Substituting above Φx,y,t into (34), we have(35)β′tcoskxxcoskyy=kx2+ky21−kx2+ky2ϵ2βtcoskxxcoskyy.Dividingcoskxxcoskyy on the both sides of (35), we obtain(36)β′t=kx2+ky21−kx2+ky2ϵ2βt.Then, the solution of the ODE (36) is as follows:(37)βt=β0expkx2+ky21−kx2+ky2ϵ2t.Let(38)Φx,y,t=βtcoskxxcoskyy,be a benchmark problem solution for the CH equation, where βt is given by (37). The initial condition is given by Φx,y,0=β0coskxxcoskyy. Finally, when conducting convergence tests for the numerical schemes of the CH equation, we consider the following modified CH equation with a source term:(39)ϕtx,y,t=ΔF′ϕx,y,t−ϵ2Δϕx,y,t+sx,y,t,where(40)sx,y,t=Φtx,y,t−ΔF′Φx,y,t−ε2ΔΦx,y,t=β3t163coskxxcos3kyykx2+9ky2+3cos3kxxcoskyy9kx2+ky2+9coskxxcoskyykx2+ky2+9cos3kxxcos3kyykx2+ky2+β′t−βt1−ε2kx2+ky2kx2+ky2coskxxcoskyy.Here, the value ofβt is given in (37) and β′t can be found by taking the first derivative of βt with respect to t. Note that if kx and ky are given satisfying kx2+ky2<1/ε2, then the numerical solution grows as time evolves as seen from (37).


## 2.2.1. One-Dimensional CH Equation

The one-dimensional CH equation is as follows:(25)ϕtx,t=ϕ3x,t−ϕx,t−ϵ2ϕxxx,txx,forx∈Ω=0,2π,t>0.with the Neumann boundary conditionϕx0,t=ϕx2π,t=ϕxxx0,t=ϕxxx2π,t=0. Similar to the benchmark problem solution of the AC equation, a linear stability analysis is first conducted around a spatially constant critical composition solution ϕ=0. Then, the linearized CH equation is as follows:(26)ϕtx,t=−ϕxxx,t−ϵ2ϕxxxxx,t.For the positive integerk, we consider Φx,t=βtcoskx as a solution for (26), where βt is an amplitude. Substituting above Φx,t into (26), we have(27)β′tcoskx=k2βtcoskx−k4ϵ2βtcoskx.Dividingcoskx on the both sides of (27), we obtain(28)β′t=k21−kϵ2βt.Then, the solution of the ODE (28) is as follows:(29)βt=β0expk21−kϵ2t.Let(30)Φx,t=βtcoskx.be a benchmark problem solution for the CH equation, whereβt is given by (29). The initial condition is given by Φx,0=β0coskx. Finally, to conduct convergence tests for numerical schemes of the CH equation, the following modified CH equation with a source term is considered:(31)ϕtx,t=F′ϕx,t−ϵ2ϕxxx,txx+sx,tfort>0,where(32)sx,t=Φtx,t−F′Φx,t−ϵ2Φxxx,txx=β′t+3k2β3t4−k2βt+k4ϵ2βtcoskx+9k2β3t4cos3kx,where the triple angle identity for cosine is applied, that is, cos3kx=cos3kx+3coskx/4. Here, the value of βt is given in (29) and β′t can be found by taking the first derivative of βt with respect to t. Note that if k is chosen to satisfy k<1/ε, then the numerical solution grows as time evolves as seen from (29).


## 2.2.2. Two-Dimensional CH Equation

The two-dimensional CH equation is as follows:(33)ϕtx,y,t=Δϕ3x,y,t−ϕx,y,t−ε2Δϕx,y,t,forx,y∈Ω=0,2π×0,2π,t>0.With the Neumann boundary conditionn⋅∇ϕ=0 on ∂Ω, where n is the exterior normal vector to the domain boundary ∂Ω. A linear stability analysis is conducted around a spatially constant critical composition solution ϕ=0. Then, the linearized CH equation is as follows:(34)ϕtx,y,t=Δ−ϕx,y,t−ϵ2Δϕx,y,t.For the positive integerskx and ky, Φx,y,t=βtcoskxxcoskyy is given as a solution for (34), where βt is an amplitude. Substituting above Φx,y,t into (34), we have(35)β′tcoskxxcoskyy=kx2+ky21−kx2+ky2ϵ2βtcoskxxcoskyy.Dividingcoskxxcoskyy on the both sides of (35), we obtain(36)β′t=kx2+ky21−kx2+ky2ϵ2βt.Then, the solution of the ODE (36) is as follows:(37)βt=β0expkx2+ky21−kx2+ky2ϵ2t.Let(38)Φx,y,t=βtcoskxxcoskyy,be a benchmark problem solution for the CH equation, where βt is given by (37). The initial condition is given by Φx,y,0=β0coskxxcoskyy. Finally, when conducting convergence tests for the numerical schemes of the CH equation, we consider the following modified CH equation with a source term:(39)ϕtx,y,t=ΔF′ϕx,y,t−ϵ2Δϕx,y,t+sx,y,t,where(40)sx,y,t=Φtx,y,t−ΔF′Φx,y,t−ε2ΔΦx,y,t=β3t163coskxxcos3kyykx2+9ky2+3cos3kxxcoskyy9kx2+ky2+9coskxxcoskyykx2+ky2+9cos3kxxcos3kyykx2+ky2+β′t−βt1−ε2kx2+ky2kx2+ky2coskxxcoskyy.Here, the value ofβt is given in (37) and β′t can be found by taking the first derivative of βt with respect to t. Note that if kx and ky are given satisfying kx2+ky2<1/ε2, then the numerical solution grows as time evolves as seen from (37).


## 3. Numerical Experiments

In this section, we compute the numerical solutions of the two phase-field equations and verify the accuracy of the proposed method. A multigrid algorithm [19–21] is used to solve the discrete equations. Let the l2 -norm errors (l2-error) be defined as eNxNt2=1/Nx∑i=1NxeiNt2, where eiNt=ϕiNt−Φxi,T for i=1,…,Nx in one-dimensional space, and eNx,NyNt2=1/NxNy∑i=1Nx∑j=1NyeijNt2, where eijNt=ϕijNt−Φxi,yj,T for i=1,…,Nx,j=1,…,Ny in two-dimensional space.


### 3.1. Convergence Test for One-Dimensional AC Equation

A convergence test of the proposed method is conducted for the one-dimensional AC equation. Letϕin be approximation of ϕxi,tn, where xi=i−0.5h, i=1,2,…,Nx and h=2π/Nx. Equation (1) is discretized by applying the θ -method [22] as follows:(41)ϕin+1−ϕinΔt=1−θ−ϕin3−ϕinε2+Δhϕin+θ−ϕin+13−ϕin+1ε2+Δhϕin+1+sin+1/2,where Δhϕin=ϕi−1n−2ϕin+ϕi+1n/h2. Here, if θ=0.5, then (41) is the Crank–Nicolson (CN) method [23]; else if θ=1, then it is the fully implicit Euler method. The Neumann boundary condition is applied, thus, ϕ0n=ϕ1n and ϕNx+1n=ϕNxn for all n=0,1,…. The discretized (41) is solved by using the multigrid algorithm. The nonlinear part ϕin+13 of F′ϕin+1 is linearized at ϕim, therefore, ϕin+13≈ϕim3+3ϕim2ϕin+1−ϕim. Here, m is the Gauss–Seidel relaxation step in the multigrid algorithm.The initial condition isϕx,0=0.2cos2x on Ω=0,2π, thus, Φx,T=0.2e1/ε2−4Tcos2x on Ω=0,2π is a benchmark problem solution for the AC equation. The multigrid algorithm parameters are set as follows: the number of Gauss–Seidel relaxation iteration =3, the tolerance =1.0e -8, and the maximum number of iteration =300.


#### 3.1.1. Fully Implicit Method

Table1 shows the convergence test results of the fully implicit method θ=1 for time step, with T=1.0e −5 and various temporal step sizes Δt=T/Nt where Nt=8,16,32, and 64. Other parameters are fixed as follows: Nx=2048, h=2π/Nx, and ε=h.Table 1
Convergence test results of the fully implicit method for time step.
Case(h, Δt)Rate(h, Δt/2)Rate(h, Δt/4)Rate(h, Δt/8)l2-error1.6730e-21.028.2512e-31.014.0976e-31.002.0419e-3Table2 shows the convergence test results of the fully implicit method for space step, with T=1.0e−5 and various spatial step sizes h=2π/Nx where Nx=8,16,32, and 64. Other parameters are fixed as follows: Nt=2048, Δt=T/Nt, and ε=π/32.Table 2
Convergence test results of the fully implicit method for space step.
Case(h, Δt)Rate(h/2, Δt)Rate(h/4, Δt)Rate(h/8, Δt)l2-error1.0727e-61.912.8517e-71.987.2424e-81.991.8199e-8


#### 3.1.2. Crank–Nicolson Method

Table3 shows the convergence test results of the CN method θ=0.5 for time step, with T=1.0e -5 and various temporal step sizes Δt=T/Nt where Nt=8,16,32, and 64. Other parameters are fixed as follows: Nx=2048, h=2π/Nx, and ε=h.Table 3
Convergence test results of the CN method for time step.
Case(h, Δt)Rate(h, Δt/2)Rate(h, Δt/4)Rate(h, Δt/8)l2-error3.5295e-41.998.8705e-52.002.2206e-52.005.5533e-6Table4 shows the convergence test results of the CN method for space step, with T=1.0e -5 and various spatial step sizes h=2π/Nx where Nx=8,16,32, and 64. Other parameters are fixed as follows: Nt=2048, Δt=T/Nt, and ε=π/32.Table 4
Convergence test results of the CN method for space step.
Case(h, Δt)Rate(h/2, Δt)Rate(h/4, Δt)Rate(h/8, Δt)l2-error1.0726e-61.912.8514e-71.987.2393e-81.991.8168e-8


### 3.2. Convergence Test for Two-Dimensional AC Equation

Letϕin be approximation of ϕxi,yj,tn, where xi=i−0.5h, i=1,2,…,Nx, yj=j−0.5h, j=1,2,…,Ny and h=2π/Nx=2π/Ny. The two-dimensional AC equation is discretized by applying the θ -method as follows:(42)ϕijn+1−ϕijnΔt=1−θ−F′ϕijnε2+Δhϕijn+θ−F′ϕijn+1ε2+Δhϕijn+1+sijn+1/2,where 0≤θ≤1, Δhϕijn=ϕi+1n+ϕi−1n+ϕi,j+1n+ϕi,j−1n−4ϕijn/h2. We set the ghost points as ϕ0,jn=ϕ1,jn, ϕNx+1,jn=ϕNx,jn, j=1,2,…,Ny and ϕi,0n=ϕi,1n, ϕi,Ny+1n=ϕi,Nyn, i=1,2,…,Nx for all n=0,1,…, because the Neumann boundary condition is adopted. The discretized governing equation is solved by using the multigrid method.The initial condition isϕx,y,0=0.2cos2xcos2y on Ω=0,2π×0,2π, thus, Φx,y,T=0.2e1/ε2−8Tcos2xcos2y on Ω=0,2π×0,2π is the benchmark problem solution for the two-dimensional AC equation. The multigrid algorithm parameters are taken as follows: the number of Gauss–Seidel relaxation iteration =3, the tolerance =1.0e -8, and the maximum number of iteration =300.


#### 3.2.1. Fully Implicit Method

Table5 shows the convergence test results of the fully implicit method for time step, with T=1.0e -5 and various temporal step sizes Δt=T/Nt where Nt=8,16,32, and 64. Other parameters are fixed as follows: Nx=Ny=2048, h=2π/Nx, and ε=h.Table 5
Convergence test results of the fully implicit method for time step.
Case(h, Δt)Rate(h, Δt/2)Rate(h, Δt/4)Rate(h, Δt/8)l2-error1.4487e-21.047.0306e-31.023.4646e-31.011.7199e-3Table6 shows the convergence test results of the fully implicit method for space step, with T=1.0e -5 and various spatial step sizes h=2π/Nx=2π/Ny where Nx=Ny=8,16,32, and 64. Other parameters are fixed as follows: Nt=2048, Δt=T/Nt, and ε=π/32.Table 6
Convergence test results of the fully implicit method for space step.
Case(h, Δt)Rate(h/2, Δt)Rate(h/4, Δt)Rate(h/8, Δt)l2-error1.5169e-61.914.0326e-71.981.0240e-71.992.5714e-8


#### 3.2.2. Crank–Nicolson Method

Table7 shows the convergence test results of the CN method for time step, with T=1.0e -5 and various temporal step sizes Δt=T/Nt where Nt=8,16,32, and 64. Other parameters are fixed as follows: Nx=Ny=2048, h=2π/Nx, and ϵ=h.Table 7
Convergence test results of the CN method for time step.
Case(h, Δt)Rate(h, Δt/2)Rate(h, Δt/4)Rate(h, Δt/8)l2-error2.1410e-42.005.3658e-52.001.3423e-52.003.3562e-6Table8 shows the convergence test results of the CN method for space step, with T=1.0e -5 and various spatial step sizes h=2π/Nx where Nx=8,16,32, and 64. Other parameters are fixed as follows: Nt=2048, Δt=T/Nt, and ε=π/32.Table 8
Convergence test results of the CN method for space step.
Case(h, Δt)Rate(h/2, Δt)Rate(h/4, Δt)Rate(h/8, Δt)l2-error1.5169e-61.914.0324e-71.981.0238e-71.992.5693e-8It is demonstrated that the calculated numerical results are appropriate benchmark problems for verifying the correctness of the numerical methods for the AC equation.


### 3.3. Convergence Test for One-Dimensional CH Equation

We consider the convergence tests of numerical schemes for CH equation. Equation (2) is discretized by applying the nonlinearly convex splitting method [24] as follows:(43)ϕin+1−ϕinΔt=Δhμin+1,μin+1=ϕin+13−ϕin−ϵ2Δhϕin+1+sin+12,where i=1,2,…,Nx. Here, the ghost points are set as ϕ0n=ϕ1n, ϕNx+1n=ϕNxn, μ0n=μ1n and μNx+1n=μNxn for all n=0,1,…, because the Neumann boundary condition is applied. The multigrid algorithm is used to solve (43)The initial condition isϕx,0=0.2cos2x on Ω=0,2π, thus, Φxi,T=0.2e4−16ϵ2Tcos2x on Ω=0,2π is the benchmark problem solution for the CH equation. The multigrid algorithm parameters are set as follows: the number of Gauss–Seidel relaxation iteration =3, the tolerance =1.0e -8, and the maximum number of iteration =100.Table9 shows the convergence test results of the unconditionally stable method [25] for time step, with T=1.0e -2 and various temporal step sizes Δt=T/Nt where Nt=8,16,32, and 64. Other parameters are fixed as follows: Nx=2048, h=2π/Nx, and ε=8h.Table 9
Convergence test results of the unconditionally stable method for time step.
Case(h, Δt)Rate(h, Δt/2)Rate(h, Δt/4)Rate(h, Δt/8)l2-error1.6655e-51.008.3544e-61.004.1886e-60.992.1017e-6Table10 shows the convergence test results of the unconditionally stable method for space step, with T=1.0e -2 and various spatial step sizes h=2π/Nx where Nx=8,16,32, and 64. Other parameters are fixed as follows: Nt=2048, Δt=T/Nt, and ε=π/32.Table 10
Convergence test results of the unconditionally stable method for space step.
Case(h, Δt)Rate(h/2, Δt)Rate(h/4, Δt)Rate(h/8, Δt)l2-error1.4948e-32.103.4801e-41.889.4222e-51.972.4093e-05


### 3.4. Convergence Test for Two-Dimensional CH Equation

Next, the two-dimensional CH equation is also discretized by applying the nonlinearly convex splitting method as follows:(44)ϕijn+1−ϕijnΔt=Δhμijn+1,μijn+1=ϕijn+13−ϕijn−ε2Δhϕijn+1+sijn+1/2,where i=1,2,…,Nx, j=1,2,…,Ny. The ghost points are set as ϕ0,jn=ϕ1,jn, ϕNx+1,jn=ϕNx,jn, μ0,jn=μ1,jn, μNx+1,jn=μNx,jn, j=1,2,…,Ny and ϕi,0n=ϕi,1n, ϕi,Nx+1n=ϕi,Nxn, μi,0n=μi,1n, μi,Nx+1n=μi,Nxn, i=1,2,…,Nx, for all n=1,2,…, because of the Neumann boundary condition. The multigrid algorithm is used to solve (44)The initial condition isϕx,y,t=0.2cos2xcos2y on Ω=0,2π×0,2π, thus Φx,y,T=0.2e8−64ε2Tcos2xcos2y on Ω=0,2π×0,2π is the benchmark problem solution for the CH equation. The multigrid algorithm parameters are taken as follows: the number of Gauss–Seidel relaxation iteration =3, the tolerance =1.0e -8, and the maximum number of iteration =100.Table11 shows the convergence test results of the unconditionally stable method for time step, with T=1.0e -2 and various temporal step sizes Δt=T/Nt where Nt=8,16,32, and 64. Other parameters are fixed as follows: Nx=Ny=2048, h=2π/Nx, and ε=8h.Table 11
Convergence test results of the unconditionally stable method for time step.
Case(h, Δt)Rate(h, Δt/2)Rate(h, Δt/4)Rate(h, Δt/8)l2-error9.5984e-20.994.8199e-21.002.4169e-21.001.2119e-2Table12 shows the convergence test results of the unconditionally stable method for space step, with T=1.0e -2 and various spatial step sizes h=2π/Nx=2π/Ny where Nx=Ny=8,16,32, and 64. Other parameters are fixed as follows: Nt=2048, Δt=T/Nt, and ε=π/32.Table 12
Convergence test results of the unconditionally stable method for space step.
Case(h, Δt)Rate(h/2, Δt)Rate(h/4, Δt)Rate(h/8, Δt)l2-error1.7444e-32.084.1116e-41.931.0805e-41.982.7467e-5It is demonstrated that the calculated numerical results are appropriate benchmark problems for verifying the correctness of the nonlinear convex splitting method for the CH equation.


### 3.5. Convergence Test for Two-Dimensional CH Equation Based the SAV Approach

The two-dimensional CH equation is discretized by applying the SAV approach [5]. The SAV approach is a step-by-step solving approach, and it can achieve both temporally first- and second-order accuracy. In this section, we consider the first-order SAV approach and perform the convergence test. For a detailed description of the SAV approach, see [5]. First, we define a time-dependent variable U as Ut=∫ΩFϕdx,Fϕ=0.25ϕ2−12.Next, letϕn+1=ϕ¯n+1+Un+1ϕ^n+1, μn+1=μ¯n+1+Un+1μ^n+1. Therefore, we can rewrite (44) as(45)ϕ¯n+1−ϕnΔt=Δμ¯n+1+sijn+1/2,μ¯n+1=−ε2Δϕ¯n+1+Sϕ¯n+1−ϕn,(46)ϕ^n+1Δt=Δμ^n+1,μ^n+1=−ϵ2Δϕ^n+1+F′ϕn∫ΩFϕndx+Sϕ^n+1,where i=1,2,…,Nx, j=1,2,…,Ny and S is a positive constant that plays a role of stabilization. The ghost points are set as ϕ¯0,jn=ϕ¯1,jn, ϕ¯Nx+1,jn=ϕ¯Nx,jn, μ¯0,jn=μ¯1,jn, μ¯Nx+1,jn=μ¯Nx,jn, j=1,2,…,Ny and ϕ¯i,0n=ϕ¯i,1n, ϕ¯i,Nx+1n=ϕ¯i,Nxn, μ¯i,0n=μ¯i,1n, μ¯i,Nx+1n=μ¯i,Nxn, i=1,2,…,Nx for all n=1,2,…, and ghost points for ϕ^n, μ^n are set similarly to ϕ¯n, ϕ¯n, respectively. The multigrid algorithm is used to solve CH equation with SAV approach. We use ϕx,y,0=0.2cos2xcos2y on Ω=0,2π×0,2π, therefore Φx,y,T=0.2e8−64ϵ2Tcos2xcos2y on Ω=0,2π×0,2π is the benchmark problem solution for the CH equation and stabilization constant is S=2. The multigrid algorithm parameters are taken as follows: the number of Gauss–Seidel relaxation iteration =3, the tolerance =1.0e -8, and the maximum number of iteration =100.Table13 shows the first-order convergence test results of the SAV approach for time step, with T=1.0e -2 and various temporal step sizes Δt=T/Nt, where Nt=8,16,32, and 64. Other parameters are fixed as follows: Nx=Ny=2048, h=2π/Nx, and ϵ=8h.Table 13
Convergence test results of the SAV approach for time step.
Case(h, Δt)Rate(h, Δt/2)Rate(h, Δt/4)Rate(h, Δt/8)l2-error1.3104e-030.986.6381e-040.993.3418e-040.991.6772e-04Table14 shows the second-order convergence test results of the SAV approach for space step, with T=1.0e -2 and various spatial step sizes h=2π/Nx=2π/Ny where Nx=Ny=8,16,32, and 64. Other parameters are fixed as follows: Nt=2048, Δt=T/Nt, and ε=π/32.Table 14
Convergence test results of the SAV approach for space step.
Case(h, Δt)Rate(h/2, Δt)Rate(h/4, Δt)Rate(h/8, Δt)l2-error1.0963e-022.082.5863e-031.926.8208e-041.961.7592e-04It is demonstrated that the calculated numerical experiment results are appropriate benchmark problems for verifying the convergence rates of the SAV approach for the CH equation.


## 3.1. Convergence Test for One-Dimensional AC Equation

A convergence test of the proposed method is conducted for the one-dimensional AC equation. Letϕin be approximation of ϕxi,tn, where xi=i−0.5h, i=1,2,…,Nx and h=2π/Nx. Equation (1) is discretized by applying the θ -method [22] as follows:(41)ϕin+1−ϕinΔt=1−θ−ϕin3−ϕinε2+Δhϕin+θ−ϕin+13−ϕin+1ε2+Δhϕin+1+sin+1/2,where Δhϕin=ϕi−1n−2ϕin+ϕi+1n/h2. Here, if θ=0.5, then (41) is the Crank–Nicolson (CN) method [23]; else if θ=1, then it is the fully implicit Euler method. The Neumann boundary condition is applied, thus, ϕ0n=ϕ1n and ϕNx+1n=ϕNxn for all n=0,1,…. The discretized (41) is solved by using the multigrid algorithm. The nonlinear part ϕin+13 of F′ϕin+1 is linearized at ϕim, therefore, ϕin+13≈ϕim3+3ϕim2ϕin+1−ϕim. Here, m is the Gauss–Seidel relaxation step in the multigrid algorithm.The initial condition isϕx,0=0.2cos2x on Ω=0,2π, thus, Φx,T=0.2e1/ε2−4Tcos2x on Ω=0,2π is a benchmark problem solution for the AC equation. The multigrid algorithm parameters are set as follows: the number of Gauss–Seidel relaxation iteration =3, the tolerance =1.0e -8, and the maximum number of iteration =300.


### 3.1.1. Fully Implicit Method

Table1 shows the convergence test results of the fully implicit method θ=1 for time step, with T=1.0e −5 and various temporal step sizes Δt=T/Nt where Nt=8,16,32, and 64. Other parameters are fixed as follows: Nx=2048, h=2π/Nx, and ε=h.Table 1
Convergence test results of the fully implicit method for time step.
Case(h, Δt)Rate(h, Δt/2)Rate(h, Δt/4)Rate(h, Δt/8)l2-error1.6730e-21.028.2512e-31.014.0976e-31.002.0419e-3Table2 shows the convergence test results of the fully implicit method for space step, with T=1.0e−5 and various spatial step sizes h=2π/Nx where Nx=8,16,32, and 64. Other parameters are fixed as follows: Nt=2048, Δt=T/Nt, and ε=π/32.Table 2
Convergence test results of the fully implicit method for space step.
Case(h, Δt)Rate(h/2, Δt)Rate(h/4, Δt)Rate(h/8, Δt)l2-error1.0727e-61.912.8517e-71.987.2424e-81.991.8199e-8


### 3.1.2. Crank–Nicolson Method

Table3 shows the convergence test results of the CN method θ=0.5 for time step, with T=1.0e -5 and various temporal step sizes Δt=T/Nt where Nt=8,16,32, and 64. Other parameters are fixed as follows: Nx=2048, h=2π/Nx, and ε=h.Table 3
Convergence test results of the CN method for time step.
Case(h, Δt)Rate(h, Δt/2)Rate(h, Δt/4)Rate(h, Δt/8)l2-error3.5295e-41.998.8705e-52.002.2206e-52.005.5533e-6Table4 shows the convergence test results of the CN method for space step, with T=1.0e -5 and various spatial step sizes h=2π/Nx where Nx=8,16,32, and 64. Other parameters are fixed as follows: Nt=2048, Δt=T/Nt, and ε=π/32.Table 4
Convergence test results of the CN method for space step.
Case(h, Δt)Rate(h/2, Δt)Rate(h/4, Δt)Rate(h/8, Δt)l2-error1.0726e-61.912.8514e-71.987.2393e-81.991.8168e-8


## 3.1.1. Fully Implicit Method

Table1 shows the convergence test results of the fully implicit method θ=1 for time step, with T=1.0e −5 and various temporal step sizes Δt=T/Nt where Nt=8,16,32, and 64. Other parameters are fixed as follows: Nx=2048, h=2π/Nx, and ε=h.Table 1
Convergence test results of the fully implicit method for time step.
Case(h, Δt)Rate(h, Δt/2)Rate(h, Δt/4)Rate(h, Δt/8)l2-error1.6730e-21.028.2512e-31.014.0976e-31.002.0419e-3Table2 shows the convergence test results of the fully implicit method for space step, with T=1.0e−5 and various spatial step sizes h=2π/Nx where Nx=8,16,32, and 64. Other parameters are fixed as follows: Nt=2048, Δt=T/Nt, and ε=π/32.Table 2
Convergence test results of the fully implicit method for space step.
Case(h, Δt)Rate(h/2, Δt)Rate(h/4, Δt)Rate(h/8, Δt)l2-error1.0727e-61.912.8517e-71.987.2424e-81.991.8199e-8


## 3.1.2. Crank–Nicolson Method

Table3 shows the convergence test results of the CN method θ=0.5 for time step, with T=1.0e -5 and various temporal step sizes Δt=T/Nt where Nt=8,16,32, and 64. Other parameters are fixed as follows: Nx=2048, h=2π/Nx, and ε=h.Table 3
Convergence test results of the CN method for time step.
Case(h, Δt)Rate(h, Δt/2)Rate(h, Δt/4)Rate(h, Δt/8)l2-error3.5295e-41.998.8705e-52.002.2206e-52.005.5533e-6Table4 shows the convergence test results of the CN method for space step, with T=1.0e -5 and various spatial step sizes h=2π/Nx where Nx=8,16,32, and 64. Other parameters are fixed as follows: Nt=2048, Δt=T/Nt, and ε=π/32.Table 4
Convergence test results of the CN method for space step.
Case(h, Δt)Rate(h/2, Δt)Rate(h/4, Δt)Rate(h/8, Δt)l2-error1.0726e-61.912.8514e-71.987.2393e-81.991.8168e-8


## 3.2. Convergence Test for Two-Dimensional AC Equation

Letϕin be approximation of ϕxi,yj,tn, where xi=i−0.5h, i=1,2,…,Nx, yj=j−0.5h, j=1,2,…,Ny and h=2π/Nx=2π/Ny. The two-dimensional AC equation is discretized by applying the θ -method as follows:(42)ϕijn+1−ϕijnΔt=1−θ−F′ϕijnε2+Δhϕijn+θ−F′ϕijn+1ε2+Δhϕijn+1+sijn+1/2,where 0≤θ≤1, Δhϕijn=ϕi+1n+ϕi−1n+ϕi,j+1n+ϕi,j−1n−4ϕijn/h2. We set the ghost points as ϕ0,jn=ϕ1,jn, ϕNx+1,jn=ϕNx,jn, j=1,2,…,Ny and ϕi,0n=ϕi,1n, ϕi,Ny+1n=ϕi,Nyn, i=1,2,…,Nx for all n=0,1,…, because the Neumann boundary condition is adopted. The discretized governing equation is solved by using the multigrid method.The initial condition isϕx,y,0=0.2cos2xcos2y on Ω=0,2π×0,2π, thus, Φx,y,T=0.2e1/ε2−8Tcos2xcos2y on Ω=0,2π×0,2π is the benchmark problem solution for the two-dimensional AC equation. The multigrid algorithm parameters are taken as follows: the number of Gauss–Seidel relaxation iteration =3, the tolerance =1.0e -8, and the maximum number of iteration =300.


### 3.2.1. Fully Implicit Method

Table5 shows the convergence test results of the fully implicit method for time step, with T=1.0e -5 and various temporal step sizes Δt=T/Nt where Nt=8,16,32, and 64. Other parameters are fixed as follows: Nx=Ny=2048, h=2π/Nx, and ε=h.Table 5
Convergence test results of the fully implicit method for time step.
Case(h, Δt)Rate(h, Δt/2)Rate(h, Δt/4)Rate(h, Δt/8)l2-error1.4487e-21.047.0306e-31.023.4646e-31.011.7199e-3Table6 shows the convergence test results of the fully implicit method for space step, with T=1.0e -5 and various spatial step sizes h=2π/Nx=2π/Ny where Nx=Ny=8,16,32, and 64. Other parameters are fixed as follows: Nt=2048, Δt=T/Nt, and ε=π/32.Table 6
Convergence test results of the fully implicit method for space step.
Case(h, Δt)Rate(h/2, Δt)Rate(h/4, Δt)Rate(h/8, Δt)l2-error1.5169e-61.914.0326e-71.981.0240e-71.992.5714e-8


### 3.2.2. Crank–Nicolson Method

Table7 shows the convergence test results of the CN method for time step, with T=1.0e -5 and various temporal step sizes Δt=T/Nt where Nt=8,16,32, and 64. Other parameters are fixed as follows: Nx=Ny=2048, h=2π/Nx, and ϵ=h.Table 7
Convergence test results of the CN method for time step.
Case(h, Δt)Rate(h, Δt/2)Rate(h, Δt/4)Rate(h, Δt/8)l2-error2.1410e-42.005.3658e-52.001.3423e-52.003.3562e-6Table8 shows the convergence test results of the CN method for space step, with T=1.0e -5 and various spatial step sizes h=2π/Nx where Nx=8,16,32, and 64. Other parameters are fixed as follows: Nt=2048, Δt=T/Nt, and ε=π/32.Table 8
Convergence test results of the CN method for space step.
Case(h, Δt)Rate(h/2, Δt)Rate(h/4, Δt)Rate(h/8, Δt)l2-error1.5169e-61.914.0324e-71.981.0238e-71.992.5693e-8It is demonstrated that the calculated numerical results are appropriate benchmark problems for verifying the correctness of the numerical methods for the AC equation.


## 3.2.1. Fully Implicit Method

Table5 shows the convergence test results of the fully implicit method for time step, with T=1.0e -5 and various temporal step sizes Δt=T/Nt where Nt=8,16,32, and 64. Other parameters are fixed as follows: Nx=Ny=2048, h=2π/Nx, and ε=h.Table 5
Convergence test results of the fully implicit method for time step.
Case(h, Δt)Rate(h, Δt/2)Rate(h, Δt/4)Rate(h, Δt/8)l2-error1.4487e-21.047.0306e-31.023.4646e-31.011.7199e-3Table6 shows the convergence test results of the fully implicit method for space step, with T=1.0e -5 and various spatial step sizes h=2π/Nx=2π/Ny where Nx=Ny=8,16,32, and 64. Other parameters are fixed as follows: Nt=2048, Δt=T/Nt, and ε=π/32.Table 6
Convergence test results of the fully implicit method for space step.
Case(h, Δt)Rate(h/2, Δt)Rate(h/4, Δt)Rate(h/8, Δt)l2-error1.5169e-61.914.0326e-71.981.0240e-71.992.5714e-8


## 3.2.2. Crank–Nicolson Method

Table7 shows the convergence test results of the CN method for time step, with T=1.0e -5 and various temporal step sizes Δt=T/Nt where Nt=8,16,32, and 64. Other parameters are fixed as follows: Nx=Ny=2048, h=2π/Nx, and ϵ=h.Table 7
Convergence test results of the CN method for time step.
Case(h, Δt)Rate(h, Δt/2)Rate(h, Δt/4)Rate(h, Δt/8)l2-error2.1410e-42.005.3658e-52.001.3423e-52.003.3562e-6Table8 shows the convergence test results of the CN method for space step, with T=1.0e -5 and various spatial step sizes h=2π/Nx where Nx=8,16,32, and 64. Other parameters are fixed as follows: Nt=2048, Δt=T/Nt, and ε=π/32.Table 8
Convergence test results of the CN method for space step.
Case(h, Δt)Rate(h/2, Δt)Rate(h/4, Δt)Rate(h/8, Δt)l2-error1.5169e-61.914.0324e-71.981.0238e-71.992.5693e-8It is demonstrated that the calculated numerical results are appropriate benchmark problems for verifying the correctness of the numerical methods for the AC equation.


## 3.3. Convergence Test for One-Dimensional CH Equation

We consider the convergence tests of numerical schemes for CH equation. Equation (2) is discretized by applying the nonlinearly convex splitting method [24] as follows:(43)ϕin+1−ϕinΔt=Δhμin+1,μin+1=ϕin+13−ϕin−ϵ2Δhϕin+1+sin+12,where i=1,2,…,Nx. Here, the ghost points are set as ϕ0n=ϕ1n, ϕNx+1n=ϕNxn, μ0n=μ1n and μNx+1n=μNxn for all n=0,1,…, because the Neumann boundary condition is applied. The multigrid algorithm is used to solve (43)The initial condition isϕx,0=0.2cos2x on Ω=0,2π, thus, Φxi,T=0.2e4−16ϵ2Tcos2x on Ω=0,2π is the benchmark problem solution for the CH equation. The multigrid algorithm parameters are set as follows: the number of Gauss–Seidel relaxation iteration =3, the tolerance =1.0e -8, and the maximum number of iteration =100.Table9 shows the convergence test results of the unconditionally stable method [25] for time step, with T=1.0e -2 and various temporal step sizes Δt=T/Nt where Nt=8,16,32, and 64. Other parameters are fixed as follows: Nx=2048, h=2π/Nx, and ε=8h.Table 9
Convergence test results of the unconditionally stable method for time step.
Case(h, Δt)Rate(h, Δt/2)Rate(h, Δt/4)Rate(h, Δt/8)l2-error1.6655e-51.008.3544e-61.004.1886e-60.992.1017e-6Table10 shows the convergence test results of the unconditionally stable method for space step, with T=1.0e -2 and various spatial step sizes h=2π/Nx where Nx=8,16,32, and 64. Other parameters are fixed as follows: Nt=2048, Δt=T/Nt, and ε=π/32.Table 10
Convergence test results of the unconditionally stable method for space step.
Case(h, Δt)Rate(h/2, Δt)Rate(h/4, Δt)Rate(h/8, Δt)l2-error1.4948e-32.103.4801e-41.889.4222e-51.972.4093e-05


## 3.4. Convergence Test for Two-Dimensional CH Equation

Next, the two-dimensional CH equation is also discretized by applying the nonlinearly convex splitting method as follows:(44)ϕijn+1−ϕijnΔt=Δhμijn+1,μijn+1=ϕijn+13−ϕijn−ε2Δhϕijn+1+sijn+1/2,where i=1,2,…,Nx, j=1,2,…,Ny. The ghost points are set as ϕ0,jn=ϕ1,jn, ϕNx+1,jn=ϕNx,jn, μ0,jn=μ1,jn, μNx+1,jn=μNx,jn, j=1,2,…,Ny and ϕi,0n=ϕi,1n, ϕi,Nx+1n=ϕi,Nxn, μi,0n=μi,1n, μi,Nx+1n=μi,Nxn, i=1,2,…,Nx, for all n=1,2,…, because of the Neumann boundary condition. The multigrid algorithm is used to solve (44)The initial condition isϕx,y,t=0.2cos2xcos2y on Ω=0,2π×0,2π, thus Φx,y,T=0.2e8−64ε2Tcos2xcos2y on Ω=0,2π×0,2π is the benchmark problem solution for the CH equation. The multigrid algorithm parameters are taken as follows: the number of Gauss–Seidel relaxation iteration =3, the tolerance =1.0e -8, and the maximum number of iteration =100.Table11 shows the convergence test results of the unconditionally stable method for time step, with T=1.0e -2 and various temporal step sizes Δt=T/Nt where Nt=8,16,32, and 64. Other parameters are fixed as follows: Nx=Ny=2048, h=2π/Nx, and ε=8h.Table 11
Convergence test results of the unconditionally stable method for time step.
Case(h, Δt)Rate(h, Δt/2)Rate(h, Δt/4)Rate(h, Δt/8)l2-error9.5984e-20.994.8199e-21.002.4169e-21.001.2119e-2Table12 shows the convergence test results of the unconditionally stable method for space step, with T=1.0e -2 and various spatial step sizes h=2π/Nx=2π/Ny where Nx=Ny=8,16,32, and 64. Other parameters are fixed as follows: Nt=2048, Δt=T/Nt, and ε=π/32.Table 12
Convergence test results of the unconditionally stable method for space step.
Case(h, Δt)Rate(h/2, Δt)Rate(h/4, Δt)Rate(h/8, Δt)l2-error1.7444e-32.084.1116e-41.931.0805e-41.982.7467e-5It is demonstrated that the calculated numerical results are appropriate benchmark problems for verifying the correctness of the nonlinear convex splitting method for the CH equation.


## 3.5. Convergence Test for Two-Dimensional CH Equation Based the SAV Approach

The two-dimensional CH equation is discretized by applying the SAV approach [5]. The SAV approach is a step-by-step solving approach, and it can achieve both temporally first- and second-order accuracy. In this section, we consider the first-order SAV approach and perform the convergence test. For a detailed description of the SAV approach, see [5]. First, we define a time-dependent variable U as Ut=∫ΩFϕdx,Fϕ=0.25ϕ2−12.Next, letϕn+1=ϕ¯n+1+Un+1ϕ^n+1, μn+1=μ¯n+1+Un+1μ^n+1. Therefore, we can rewrite (44) as(45)ϕ¯n+1−ϕnΔt=Δμ¯n+1+sijn+1/2,μ¯n+1=−ε2Δϕ¯n+1+Sϕ¯n+1−ϕn,(46)ϕ^n+1Δt=Δμ^n+1,μ^n+1=−ϵ2Δϕ^n+1+F′ϕn∫ΩFϕndx+Sϕ^n+1,where i=1,2,…,Nx, j=1,2,…,Ny and S is a positive constant that plays a role of stabilization. The ghost points are set as ϕ¯0,jn=ϕ¯1,jn, ϕ¯Nx+1,jn=ϕ¯Nx,jn, μ¯0,jn=μ¯1,jn, μ¯Nx+1,jn=μ¯Nx,jn, j=1,2,…,Ny and ϕ¯i,0n=ϕ¯i,1n, ϕ¯i,Nx+1n=ϕ¯i,Nxn, μ¯i,0n=μ¯i,1n, μ¯i,Nx+1n=μ¯i,Nxn, i=1,2,…,Nx for all n=1,2,…, and ghost points for ϕ^n, μ^n are set similarly to ϕ¯n, ϕ¯n, respectively. The multigrid algorithm is used to solve CH equation with SAV approach. We use ϕx,y,0=0.2cos2xcos2y on Ω=0,2π×0,2π, therefore Φx,y,T=0.2e8−64ϵ2Tcos2xcos2y on Ω=0,2π×0,2π is the benchmark problem solution for the CH equation and stabilization constant is S=2. The multigrid algorithm parameters are taken as follows: the number of Gauss–Seidel relaxation iteration =3, the tolerance =1.0e -8, and the maximum number of iteration =100.Table13 shows the first-order convergence test results of the SAV approach for time step, with T=1.0e -2 and various temporal step sizes Δt=T/Nt, where Nt=8,16,32, and 64. Other parameters are fixed as follows: Nx=Ny=2048, h=2π/Nx, and ϵ=8h.Table 13
Convergence test results of the SAV approach for time step.
Case(h, Δt)Rate(h, Δt/2)Rate(h, Δt/4)Rate(h, Δt/8)l2-error1.3104e-030.986.6381e-040.993.3418e-040.991.6772e-04Table14 shows the second-order convergence test results of the SAV approach for space step, with T=1.0e -2 and various spatial step sizes h=2π/Nx=2π/Ny where Nx=Ny=8,16,32, and 64. Other parameters are fixed as follows: Nt=2048, Δt=T/Nt, and ε=π/32.Table 14
Convergence test results of the SAV approach for space step.
Case(h, Δt)Rate(h/2, Δt)Rate(h/4, Δt)Rate(h/8, Δt)l2-error1.0963e-022.082.5863e-031.926.8208e-041.961.7592e-04It is demonstrated that the calculated numerical experiment results are appropriate benchmark problems for verifying the convergence rates of the SAV approach for the CH equation.


## 4. Conclusion

In this study, a benchmark problem is presented for the numerical technique of phase-field equations. To show the effectiveness of the proposed method, two famous phase-field equations were considered. To design an appropriate benchmark problem for the AC equation and the CH equation on the one- and two-dimensional, we first performed a linear stability analysis and then took a growth mode solution as the benchmark problem, which is closely related to the dynamics of the original governing equation. For each phase-field equation, convergence tests were conducted of the numerical schemes. The computational results confirmed the accuracy and efficiency of the proposed method. The proposed method is general; therefore, it is possible to design benchmark problems for various phase-field equations or directly extend to higher dimensions such as three-dimensional space.



---

*Source: 2751592-2022-01-21.xml*