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

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. -e oneand two-dimensional computational experiments confirm the accuracy and efficiency of the proposed scheme as the benchmark problems.


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. e first equation is the Allen-Cahn (AC) equation [1,2]: , for x ∈ Ω, t > 0. (1) with the Neumann boundary condition n · ∇ϕ � 0 on zΩ, where n is the exterior normal vector to the domain boundary zΩ. e phase-field ϕ(x, t) is the difference between the concentrations of the two mixtures components, the free energy F(ϕ) � 0.25(ϕ 2 − 1) 2 is double-well potential, and ϵ is a positive constant related to the thickness of the interfacial transition layer. e second equation is the Cahn-Hilliard (CH) equation [3]: ϕ t (x, t) � Δ F ′ (ϕ(x, t)) − ϵ 2 Δϕ(x, t) , for x ∈ Ω, t > 0. (2) with the Neumann boundary condition n · ∇ϕ � n · ∇ · Δϕ � 0 on zΩ. 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. erefore, to validate the accuracy of the newly proposed numerical methods, it is essential to test the developed methods with benchmark problems. e 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.
e 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]. e improved SAV approach has all the advantages of classical SAV approach, further simplifies the algorithm, and easily constructs the temporally firstorder 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]. e 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. ey 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. ey 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.
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. e layout of this paper is as follows. In Section 2, the procedure of finding benchmark problems is described. In Section 3, the numerical experiments are performed. In Section 4, the conclusions are given.

Proposed Benchmark Problems
In this section, the proposed benchmark problems are presented for the numerical methods of the one-and twodimensional 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
Discrete Dynamics in Nature and Society

Benchmark Problems for the AC Equation.
First, let us consider the benchmark problems for the one-and twodimensional AC equation.

One-Dimensional AC Equation.
e one-dimensional AC equation is as follows: with the Neumann boundary condition ϕ x (0, t) � ϕ x (2π, 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 ′ (ϕ) ≈ − ϕ. erefore, the linearized AC equation is as follows: For a positive integer k, we take Φ(x, t) � α(t)cos(kx) as a solution for (10), where α(t) is an amplitude. Substituting above Φ(x, t) into (10), we have Dividing cos(kx) on the both sides of (11), we obtain en, the solution of the ordinary differential equation (ODE) (12) is as follows: Let Be a benchmark problem solution for a one-dimensional modified AC equation, where α(t) is given by (13). Here, Φ(x, 0) � α(0)cos(kx) 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: where 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).

Two-Dimensional AC Equation.
e two-dimensional AC equation is as follows: with the Neumann boundary condition n · ∇ϕ � 0 on zΩ, where n is the exterior normal vector to the domain boundary zΩ. 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. e nonlinear term F′(ϕ) � ϕ 3 − ϕ is linearized by applying the Taylor expansion, and then the linearized term F ′ (ϕ) ≈ − ϕ is obtained. erefore, the linearized AC equation is as follows: For positive integers k x and k y , let us consider Dividing cos(k x x)cos(k y y) on the both sides of (19), we obtain en, the solution of the ODE (19) is as follows: be a benchmark problem solution for the modified AC equation, where α(t) is given by (21). e initial condition is given by Φ(x, y, 0) � α(0)cos(k x x)cos(k y y). Finally, when we perform convergence tests for numerical schemes of the AC equation, the following modified AC equation with a source term is taken: where 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 k x and k y satisfying k 2 x + k 2 y < 1/ε 2 , then the numerical solution grows as time evolves as seen from (21).

Benchmark Problems for the CH Equation.
Next, let us consider the benchmark problems for the one-and twodimensional CH equation.

One-Dimensional CH Equation.
e one-dimensional CH equation is as follows: 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. en, the linearized CH equation is as follows: (26) For the positive integer k, we consider Dividing cos(kx) on the both sides of (27), we obtain en, the solution of the ODE (28) is as follows: be a benchmark problem solution for the CH equation, where β(t) is given by (29). e initial condition is given by Φ(x, 0) � β(0)cos(kx). Finally, to conduct convergence tests for numerical schemes of the CH equation, the following modified CH equation with a source term is considered: where the triple angle identity for cosine is applied, that is, cos 3 (kx) � [cos(3kx) + 3 cos(kx)]/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).

(33)
With the Neumann boundary condition n · ∇ϕ � 0 on zΩ, where n is the exterior normal vector to the domain boundary zΩ. A linear stability analysis is conducted around a spatially constant critical composition solution ϕ � 0. en, the linearized CH equation is as follows: For the positive integers k x and k y , Φ(x, y, t) � β(t)cos(k x x)cos(k y y) is given as a solution for (34), where β(t) is an amplitude. Substituting above Φ(x, y, t) into (34), we have Dividing cos(k x x)cos(k y y) on the both sides of (35), we obtain en, the solution of the ODE (36) is as follows: Let Φ(x, y, t) � β(t)cos k x x cos k y y , be a benchmark problem solution for the CH equation, where β(t) is given by (37). e initial condition is given by Φ(x, y, 0) � β(0)cos(k x x)cos(k y y). Finally, when conducting convergence tests for the numerical schemes of the Discrete Dynamics in Nature and Society 5 CH equation, we consider the following modified CH equation with a source term: 16 3 cos k x x cos 3k y y k 2 x + 9k 2 y + 3 cos 3k x x cos k y y 9k 2 x + k 2 y +9 cos k x x cos k y y k 2 x + k 2 y + 9 cos 3k x x cos 3k y y k 2 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 k x and k y are given satisfying k 2 x + k 2 y < 1/ε 2 , then the numerical solution grows as time evolves as seen from (37).

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][20][21] is used to solve the discrete equations. Let the l 2 -norm errors (l 2 -error) be defined as ‖e . . , N x , j � 1, . . . , N y in two-dimensional space.

Convergence Test for One-Dimensional AC Equation.
A convergence test of the proposed method is conducted for the one-dimensional AC equation. Let ϕ n i be approximation (1) is discretized by applying the θ -method [22] as follows: where Δ h ϕ n i � (ϕ n i− 1 − 2ϕ n i + ϕ n i+1 )/h 2 . Here, if θ � 0.5, then (41) is the Crank-Nicolson (CN) method [23]; else if θ � 1, then it is the fully implicit Euler method.
e 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. Table 1 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/N t where N t � 8, 16, 32, and 64. Other parameters are fixed as follows: N x � 2048, h � 2π/N x , and ε � h. Table 2 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π/N x where N x � 8, 16, 32, and 64. Other parameters are fixed as follows: N t � 2048, Δt � T/N t , and ε � π/32. Table 3 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/N t where N t � 8, 16, 32, and 64. Other parameters are fixed as follows: N x � 2048, h � 2π/N x , and ε � h. Table 4 shows the convergence test results of the CN method for space step, with T � 1.0e -5 and various spatial step sizes h � 2π/N x where N x � 8, 16, 32, and 64. Other parameters are fixed as follows: N t � 2048, Δt � T/N t , and ε � π/32.

Convergence Test for Two-Dimensional AC Equation.
Let ϕ n i be approximation of ϕ(x i , y j , t n ), where Discrete Dynamics in Nature and Society two-dimensional AC equation is discretized by applying the θ -method as follows: We set the ghost points as ϕ n 0,j � ϕ n 1,j , ϕ n N x +1,j � ϕ n N x ,j , j � 1, 2, . . . , N y and ϕ n i,0 � ϕ n i,1 , ϕ n i,N y +1 � ϕ n i,N y , i � 1, 2, . . . , N x for all n � 0, 1, . . ., because the Neumann boundary condition is adopted. e discretized governing equation is solved by using the multigrid method. e initial condition is ϕ(x, y, 0) � 0.2cos(2x)cos(2y) on Ω � (0, 2π) × (0, 2π), thus, Φ(x, y, T) � 0.2e (1/ε 2 − 8)T cos(2x)cos(2y) on Ω � (0, 2π) × (0, 2π) is the benchmark problem solution for the two-dimensional AC equation. e 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. Table 5 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/N t where N t � 8, 16, 32, and 64. Other parameters are fixed as follows: N x � N y � 2048, h � 2π/N x , and ε � h. Table 6 shows the convergence test results of the fully implicit method for space step, with T � 1.0e -5 and various spatial 16, 32, and 64. Other parameters are fixed as follows: N t � 2048, Δt � T/N t , and ε � π/32. It is demonstrated that the calculated numerical results are appropriate benchmark problems for verifying the correctness of the numerical methods for the AC equation.
e 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. Table 9 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/N t where N t � 8, 16, 32, and 64. Other parameters are fixed as follows: N x � 2048, h � 2π/N x , and ε � 8h. Table 10 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π/N x where

Convergence Test for Two-Dimensional CH Equation Based the SAV Approach.
e two-dimensional CH equation is discretized by applying the SAV approach [5]. e 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 erefore, we can rewrite (44) as where i � 1, 2, . . . , N x , j � 1, 2, . . . , N y and S is a positive constant that plays a role of stabilization. e ghost points are set as ϕ . . , N x for all n � 1, 2, . . ., and ghost points for ϕ n , μ n are set similarly to ϕ n , ϕ n , respectively. e multigrid algorithm is used to solve CH equation with SAV approach. We use ϕ(x, y, 0) � 0.2cos(2x)cos(2y) on Ω � (0, 2π) × (0, 2π), therefore Φ(x, y, T) � 0.2e (8− 64ϵ 2 )T cos(2x)cos(2y) on Ω � (0, 2π) × (0, 2π) is the benchmark problem solution for the CH equation and stabilization constant is S � 2. e 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. Table 13 shows the first-order convergence test results of the SAV approach for time step, with T � 1.0e -2 and various      Table 14 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π/N x � 2π/N y where N x � N y � 8, 16, 32, and 64. Other parameters are fixed as follows: N t � 2048, Δt � T/N t , and ε � π/32. It 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.

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. e computational results confirmed the accuracy and efficiency of the proposed method. e 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.

Data Availability
e data used to support the findings of this study are included in the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.    Discrete Dynamics in Nature and Society 9