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.
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]:
with the Neumann boundary condition on , where is the exterior normal vector to the domain boundary . The phase-field is the difference between the concentrations of the two mixtures components, the free energy 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]:
with the Neumann boundary condition 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:
In [12], the authors used the following four benchmark problems:where if and otherwise. In [13], the authors used the following initial condition for the benchmark problem:
Figure 1 shows some of the initial conditions for the benchmark problems.

(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 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.
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:
with the Neumann boundary condition . To find appropriate benchmark problems for the AC equation, a linear stability analysis is first conducted around a spatially constant critical composition solution . We linearize the nonlinear term by applying the Taylor expansion and then get linearized term . Therefore, the linearized AC equation is as follows:
For a positive integer , we take as a solution for (10), where is an amplitude. Substituting above into (10), we have
Dividing on the both sides of (11), we obtain
Then, 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 is given by (13). Here, 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 is given in (13) and can be found by taking the first derivative of with respect to . Note that if is chosen to satisfy , 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:
with the Neumann boundary condition on , where 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 . The nonlinear term is linearized by applying the Taylor expansion, and then the linearized term is obtained. Therefore, the linearized AC equation is as follows:
For positive integers and , let us consider as a two-dimensional benchmark solution for (18), where is an amplitude. Substituting above into (18), then
Dividing on the both sides of (19), we obtain
Then, the solution of the ODE (19) is as follows:
Let
be a benchmark problem solution for the modified AC equation, where is given by (21). The initial condition is given by . 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 is given in (21) and can be found by taking the first derivative of with respect to . Note that if we set and satisfying , 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:
with the Neumann boundary condition . Similar to the benchmark problem solution of the AC equation, a linear stability analysis is first conducted around a spatially constant critical composition solution . Then, the linearized CH equation is as follows:
For the positive integer , we consider as a solution for (26), where is an amplitude. Substituting above into (26), we have
Dividing on the both sides of (27), we obtain
Then, the solution of the ODE (28) is as follows:
Let
be a benchmark problem solution for the CH equation, where is given by (29). The initial condition is given by . Finally, to conduct convergence tests for numerical schemes of the CH equation, the following modified CH equation with a source term is considered:wherewhere the triple angle identity for cosine is applied, that is, . Here, the value of is given in (29) and can be found by taking the first derivative of with respect to . Note that if is chosen to satisfy , 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:
With the Neumann boundary condition on , where is the exterior normal vector to the domain boundary . A linear stability analysis is conducted around a spatially constant critical composition solution . Then, the linearized CH equation is as follows:
For the positive integers and , is given as a solution for (34), where is an amplitude. Substituting above into (34), we have
Dividing on the both sides of (35), we obtain
Then, the solution of the ODE (36) is as follows:
Letbe a benchmark problem solution for the CH equation, where is given by (37). The initial condition is given by . Finally, when conducting convergence tests for the numerical schemes of the CH equation, we consider the following modified CH equation with a source term:where
Here, the value of is given in (37) and can be found by taking the first derivative of with respect to . Note that if and are given satisfying , 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 -norm errors (-error) be defined as , where for in one-dimensional space, and , where for 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 be approximation of , where , and . Equation (1) is discretized by applying the -method [22] as follows:where . Here, if , then (41) is the Crank–Nicolson (CN) method [23]; else if , then it is the fully implicit Euler method. The Neumann boundary condition is applied, thus, and for all . The discretized (41) is solved by using the multigrid algorithm. The nonlinear part of is linearized at , therefore, . Here, is the Gauss–Seidel relaxation step in the multigrid algorithm.
The initial condition is on , thus, on is a benchmark problem solution for the AC equation. The multigrid algorithm parameters are set as follows: the number of Gauss–Seidel relaxation iteration , the tolerance -8, and the maximum number of iteration .
3.1.1. Fully Implicit Method
Table 1 shows the convergence test results of the fully implicit method for time step, with −5 and various temporal step sizes where , and 64. Other parameters are fixed as follows: , , and .
Table 2 shows the convergence test results of the fully implicit method for space step, with −5 and various spatial step sizes where , and 64. Other parameters are fixed as follows: , , and .
3.1.2. Crank–Nicolson Method
Table 3 shows the convergence test results of the CN method for time step, with -5 and various temporal step sizes where , and 64. Other parameters are fixed as follows: , , and .
Table 4 shows the convergence test results of the CN method for space step, with -5 and various spatial step sizes where , and 64. Other parameters are fixed as follows: , , and .
3.2. Convergence Test for Two-Dimensional AC Equation
Let be approximation of , where , , , and . The two-dimensional AC equation is discretized by applying the -method as follows:where , . We set the ghost points as , , and , , for all , because the Neumann boundary condition is adopted. The discretized governing equation is solved by using the multigrid method.
The initial condition is on , thus, on 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 , the tolerance -8, and the maximum number of iteration .
3.2.1. Fully Implicit Method
Table 5 shows the convergence test results of the fully implicit method for time step, with -5 and various temporal step sizes where , and 64. Other parameters are fixed as follows: , , and .
Table 6 shows the convergence test results of the fully implicit method for space step, with -5 and various spatial step sizes where , and 64. Other parameters are fixed as follows: , , and .
3.2.2. Crank–Nicolson Method
Table 7 shows the convergence test results of the CN method for time step, with -5 and various temporal step sizes where , and 64. Other parameters are fixed as follows: , , and .
Table 8 shows the convergence test results of the CN method for space step, with -5 and various spatial step sizes where , and 64. Other parameters are fixed as follows: , , and .
It 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:where . Here, the ghost points are set as , , and for all , because the Neumann boundary condition is applied. The multigrid algorithm is used to solve (43)
The initial condition is on , thus, on is the benchmark problem solution for the CH equation. The multigrid algorithm parameters are set as follows: the number of Gauss–Seidel relaxation iteration , the tolerance -8, and the maximum number of iteration .
Table 9 shows the convergence test results of the unconditionally stable method [25] for time step, with -2 and various temporal step sizes where , and 64. Other parameters are fixed as follows: , , and .
Table 10 shows the convergence test results of the unconditionally stable method for space step, with -2 and various spatial step sizes where , and 64. Other parameters are fixed as follows: , , and .
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:where , . The ghost points are set as , , , , and , , , , , for all , because of the Neumann boundary condition. The multigrid algorithm is used to solve (44)
The initial condition is on , thus on is the benchmark problem solution for the CH equation. The multigrid algorithm parameters are taken as follows: the number of Gauss–Seidel relaxation iteration , the tolerance -8, and the maximum number of iteration .
Table 11 shows the convergence test results of the unconditionally stable method for time step, with -2 and various temporal step sizes where , and 64. Other parameters are fixed as follows: , , and .
Table 12 shows the convergence test results of the unconditionally stable method for space step, with -2 and various spatial step sizes where , and 64. Other parameters are fixed as follows: , , and .
It 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 as
Next, let , . Therefore, we can rewrite (44) as
where , and is a positive constant that plays a role of stabilization. The ghost points are set as , , , , and , , , , for all , and ghost points for , are set similarly to , , respectively. The multigrid algorithm is used to solve CH equation with SAV approach. We use on , therefore on is the benchmark problem solution for the CH equation and stabilization constant is . The multigrid algorithm parameters are taken as follows: the number of Gauss–Seidel relaxation iteration , the tolerance -8, and the maximum number of iteration .
Table 13 shows the first-order convergence test results of the SAV approach for time step, with -2 and various temporal step sizes , where , and 64. Other parameters are fixed as follows: , , and .
Table 14 shows the second-order convergence test results of the SAV approach for space step, with -2 and various spatial step sizes where , and 64. Other parameters are fixed as follows: , , and .
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.
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.
Data Availability
The data used to support the findings of this study are included in the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The first author (Y. Hwang) was supported by the National Research Foundation (NRF), Korea, under project BK21 FOUR. The corresponding author (J. S. Kim) was supported by Korea University Grant.