Comparison of operator splitting schemes for the numerical solution of the Allen-Cahn equation

In this paper, we have analyzed the first- and second-order operator splitting schemes for the numerical solution of the Allen-Cahn equation. Different possibilities are considered for the derivation of the schemes. The numerical errors are computed in order to verify the effectiveness of each scheme. Validation of the schemes is provided by numerically solving some well-known examples.


I. INTRODUCTION
The problem of interfacial dynamics in multiphase flows is multidisciplinary with several applications such as the impact of a droplet on a solid surface, drop coalescence, realistic interfaces in computer graphics and animation movies, tumor/cancer growth, and wound healing and dendritic crystal growth. The Allen-Cahn equation was first introduced by Allen and Cahn to describe the phase transition and coarsening phenomena in material sciences. 1 Mathematically, the Allen-Cahn equation is a second-order stiff nonlinear partial differential equation, which is used as a reactiondiffusion equation in material sciences and as a convection-diffusion equation in computational fluid dynamics. [2][3][4] This equation is extensively studied and applied to many kinds of moving interface problems. Since the exact solution of the Allen-Cahn equation is not possible always, due to its stiff nature, numerical methods play an essential role in solving this equation.
In this paper, we consider the following Allen-Cahn equation: ∂ϕ(x, y, t) ∂t = Δϕ(x, y, t) − 1 ϵ 2 f (ϕ), ϕ : [0, T] × Ω → R, ϕ(x, y, 0) = ϕ 0 (x, y) (1) subject to the appropriate boundary conditions. Here, ϕ is defined as the difference between the concentrations of the two components in a mixture (e.g., ϕ = m α −m β m α +m β , where mα and m β are the mole fractions of phases α and β). The phase variable ϕ is also known as the order parameter, which represents the local state of the entire system. For example, ϕ = 1 is in one phase, while ϕ = −1 is in the other phase, which is separated by a diffuse interface of thickness ϵ. The small positive constant ϵ is the gradient energy coefficient related to the interfacial free energy. The interface between two phases is defined by Γt = {x| ϕ(x, t) = 0} at any time t. In Eq. (1), Δ is the Laplace operator, Ω is a two-dimensional bounded domain, and f (ϕ) = F ′ (ϕ) = ϕ 3 − ϕ. The Allen-Cahn equation is the L 2 -gradient flow of the following total free energy functional: 5 The function F(ϕ) = 1 4 (ϕ 2 − 1) 2 is the Helmholtz free energy density for ϕ, which has a double-well form with two minima at ϕ = ±1. Differentiating the energy functional E(ϕ, ∇ϕ) with respect to t gives where the integration by parts and the homogeneous Neumann boundary conditions are used. Therefore, the total energy is ARTICLE scitation.org/journal/adv nonincreasing in time. The Allen-Cahn equation and its various modified forms have been applied in addressing a range of problems, such as phase transitions, 1 image analysis, 6,7 motion by mean curvature, [8][9][10] two-phase flows, 11 and crystal growth. [12][13][14] Therefore, several numerical methods have been developed for efficient and accurate solution of this equation. Numerous studies show that the explicit schemes have a severe time step restriction for stability reasons due to the nonlinear term F ′ (ϕ) = ϕ 3 − ϕ, while implicit schemes suffer from a solvability problem and error may increase with large time step. One of the semi-implicit methods proposed by Eyre 15 is an unconditionally gradient stable method, which is firstorder accurate in time. Moreover, first-and second-order stabilized semi-implicit methods have also been proposed. 16,17 In recent years, operator splitting schemes are used for solving the nonlinear partial differential equations 16,18,19 and evolution equations. [20][21][22][23][24] The main idea behind the method is that the original nonlinear problem is decomposed into simpler subproblems. These subproblems are then solved in a given sequential order to get their exact or approximate solutions. The solutions obtained are then used to compute the approximate solution of the original problem. The main attractive features of the operator splitting schemes are their easy implementation and computational efficiency. Due to the unconditional stability of each substep, the first-and the second-order operator splitting schemes are unconditionally stable. 18,19,23 In this paper, we provide the numerical stability and accuracy of the first-and secondorder operator splitting schemes for solving the two-dimensional Allen-Cahn equation. At each time level, the Allen-Cahn equation is first decomposed into a linear Laplace equation and the nonlinear partial differential equation. The linear Laplace equation is then solved using the Fourier spectral method while the nonlinear equation is solved analytically. Based on the computed solution, the stability and error estimate of the numerical schemes are discussed in L 2 -norm. The rest of this paper is organized as follows: In Sec. II, a brief review of the operator splitting schemes is given. Derivation and stability of the first-and second-order operator splitting schemes are provided in Secs. III and IV, respectively. Numerical tests to validate the proposed schemes are given in Sec. V. Section VI concludes this paper.

II. REVIEW OF THE OPERATOR SPLITTING SCHEMES
In this section, we briefly review the first-and second-order operator splitting schemes for a general time evolution equation. It is easy to construct a first-order solution ϕ(t n+1 ) of the following time evolution equation: by computing ϕ(t n+1 ) ≈ (S Δt 1 S Δt 2 )ϕ(t n ), where S Δt 1 and S Δt 2 are the solution operators for ∂ϕ/∂t = f 1 (ϕ) and ∂ϕ/∂t = f 2 (ϕ), respectively. Then, the operators S 1 and S 2 satisfy the semigroup properties. A second-order scheme can be derived simply by symmetrizing the first-order scheme, The choices of S 1 and S 2 [or f 1 (ϕ) and f 2 (ϕ)] are arbitrary and is the basis of this work. In Secs. III-V, we provide effects on the accuracy, stability, and convergence of a numerical method by choosing S 1 and S 2 in a different order.

III. THE FIRST-ORDER SCHEMES
In this section, two unconditionally stable first-order operator splitting schemes are given for solving the Allen-Cahn equation. For the space discretization, the Fourier spectral method is used. A spatial grid Ω h of size h = b−a N is considered, where N is a positive even integer, and is the set of cell centers. Let ϕ m ij = ϕ(xi, yj, mΔt), where Δt = T M is the time step and M is the total number of time steps. First, we give a definition, which is necessary for the derivation and analysis of the numerical schemes. 25 Definition. Suppose the 2D Laplacian operator (−Δ) has a complete set of orthonormal eigenfunctions upq corresponding to the eigenvalues λpq on a bounded region [a, b] 2 , i.e., (−Δ)upq = λpqupq. Let Then, for any ϕ ∈ Φ, it is easy to have where λpq and upq are given by Here, Cp and Cq are defined as Then, according to the theory of spectral methods, 26 and where C * p and C * q are defined as

ARTICLE scitation.org/journal/adv
Now, we use the first-order operator splitting scheme to decompose the original equation (1) into linear and nonlinear subequations as follows: and First, we find the pqth Fourier mode of Eq. (7) from Eq. (4) as Second, Eq. (9) is solved analytically by the method of separation of variables 27 using the initial conditionφ m pq aŝ Hence, we have where ϕ * is an intermediate solution, C represents the discrete cosine transform, and its inverse is denoted by C −1 . The analytical solution 18 of Eq. (8) is given as follows: So to update the solution ϕ m+1 from ϕ m , we use two first-order operator splitting schemes as follows: and G2 : Next, we show the unconditional stability of the scheme as follows: In formulation G1, the solution of the heat equation with initial condition ϕ o (x, y) is a function defined as a cosine series, thus satisfying the maximum principle 19 Therefore, maxp,q|ϕ * p,q | = maxx,y|ϕ * (x, y)| ≤ maxx,y|ϕ m (x, y)| = maxp,q|ϕ m p,q |. Now, assuming that |ϕ * p,q | ≤ 1, for each p, q in the second step, we, therefore, have Thus, scheme G1 is unconditionally stable since the solution remains bounded by 1. The unconditional stability of scheme G2 can also be described similarly.

IV. THE SECOND-ORDER SCHEMES
In this section, we describe two unconditionally stable secondorder operator splitting schemes for solving the Allen-Cahn equation. Both schemes are based on the Fourier spectral method for space discretization (as in the case of first-order schemes) and Strang's splitting scheme for time discretization. The spatial and temporal grids are selected in the same way as in the case of firstorder schemes described in Eqs. (12) and (13). Again, the core idea behind the given schemes is to decompose the Allen-Cahn equation into linear and nonlinear subequations given as Eqs. (7) and (8). Let S 1 and S 2 be the exact solution operators associated with Eq. (7) and Eq. (8), respectively. The solution of the Allen-Cahn equation can, therefore, be evolved in time in three substeps by the following two different ways: or This three-step splitting algorithm is second-order accurate based on the Strang's operator splitting scheme. 20 Now, the exact solution operators S 1 and S 2 are replaced by their numerical approximations S 1h and S 2h . As discussed in Sec. III, the solutions of Eqs. (14) and (15) give .
In order to study the effect of the choice of operators S 1h and S 2h , we consider the two second-order splitting methods each consisting of three steps as follows: and H2 : The stability and accuracy of the schemes H1 and H2 are studied for getting a numerical solution accurate at least within some required tolerance of the true solution. For this, we provide the following stability statement with its proof: 29 Proof. For convenience, we only give the proof for the method H1, and the case of H2 can be discussed similarly. Using Parseval's identity and the fact that e −λ pq Δt 2 < 1 ∀(p, q), it follows from the first formula of method H1 that Now, the second formula of method H1 is It follows that For the last formula of H1, an argument similar to the one used in Eq. (20) shows that Therefore, from Eqs. (20)- (22), we obtain the expected result, i.e., where T denotes the final time. This completes the proof.

V. RESULTS AND DISCUSSIONS
In this section, we compare the stability and accuracy between formulations G1 and G2 based on the operator splitting scheme for the classical benchmark problem in Ref. 28. Similarly, we also compare schemes H1 and H2 using the cosine function as an initial condition. Moreover, we have applied schemes H1 and H2 on randomly perturbed data and analyzed them.

A. Accuracy of the first-order schemes
To show that the proposed methods have the first-order accuracy in time, we solve a mean curvature problem. 28 Initially, a radially symmetric circular interface centered at (0.5, 0.5) is considered on the domain Ω = [0, 1] × [0, 1], i.e., The initial radius is Ro = 0.25, ϵ = 0.01, and the computational domain is of size 128 × 128 spatial grids. As time evolves, the radius of the circle shrinks at the rate of the curvature of the circle and eventually disappears. The moving interface velocity V(t) is given by The analytical solution is where d is the number of dimensions. The circle will shrink with time to the center of the circle and disappears for times larger than interfacial thickness, i.e., ϵ = 0.015 and ϵ = 0.01 using methods G1 and G2 and the results are recorded in Tables I and II,   where O(A) = order of scheme A, ϕe = exact solution, ϕ f = solution at finer grids, and ϕc = solution at the coarse grid. The computed results are given in Table III. It is observed that both methods yield temporal approximation errors close to one. Stability test. In this section, we take the final time T f = 0.024 with different time steps as Δt = 1 ⋅ 10 −6 , Δt = 1 ⋅ 10 −4 , Δt = 2 ⋅ 10 −4 , and Δt = 4 ⋅ 10 −4 while keeping the spatial step size uniform, i.e., h = 1/128. The radii of a circle as a function of time obtained from G1 and G2 at different time steps are shown in Figs. 3(a) and 3(b), respectively. Both the schemes have stable results even at large time steps showing that they are unconditionally stable.
These two equations are solved in sequential order corresponding to scheme G1. First, find the pqth Fourier mode of Eq. (30) from Eq. (4) as Then, Eq. (32) is solved using the initial conditionφ m pq to give Now, we solve the nonlinear equation using the solution obtained from Eq. (33) as the initial condition. The nonlinear term in Eq. (31) is solved pseudospectrally: as the solution in the physical space changes, the nonlinear term is calculated at each time step by using the discrete values in the first step to form a power ϕ 3 in the physical space and then applying the cosine transform and the Laplacian operators in Fourier space. Finally, convert back to the physical space. Then, the spatial discretization of the nonlinear partial differential equation ref splitCH2 can be written as For the time integration of this equation, we use the forward Euler method, which gives We run the simulation for a square simulation cell having N = 128 and h = 0.1. All of the parameters used in the simulation are dimensionless. 32 The initial composition was modulated by introducing a random noise term of 0.02, defined as ϕ(x, y, 0) = 0.02 (0.5 − rand(x, y)). temporal evolution of the initial composition due to phase separation is summarized in Figs. 4(a)-4(f). As can be seen at time t = 0, the microstructure is relatively fine and contains a large amount of precipitate. As time passes, it is easy to infer from the graph that the second phase is roughened by migration, dissolution, coalescence, and cracking of the phase boundary. As shown in Fig. 5, this phase separation process is the result of a reduction in total free energy.  Table IV. The spatial mesh is fixed as h = 1/64, and homogeneous Neumann boundary conditions are used. Since the analytical solution is unknown, therefore, we take a reference solution ϕe at the finest mesh as the exact one and compute the relative error as Err = ∥ϕ−ϕ e ∥ ∥ϕ e ∥ and convergence order is given by Eq. (27). Table IV shows the error and the convergence order of both schemes H1 and H2. It is observed that both methods yield temporal approximation orders close to 2. We can also see from the results reported in Table IV that H2 is more accurate than H1.

C. Spinodal decomposition
In order to further verify the effectiveness of the proposed schemes for solving the AC equation, we consider an example of spinodal decomposition simulation using the scheme H2 (other schemes can be used as well). A 256 × 256 mesh is used on the domain Ω = [−1, 1] 2 with ϵ = 0.01. The nondimensional time increment per time step was set as Δt = 0.01, and simulation was carried out up to 20 000 time steps. The initial condition is taken to be a randomly perturbed distribution, i.e., ϕ(x, y, 0) = 0.01 (rand(x, y) − 0.5).
Here, rand(x, y) is a random number generator between −1 and 1. The time evolution of the morphological patterns is summarized in Fig. 7. The time values quoted in Fig. 7 are in nondimensional form.

VI. CONCLUSIONS
In this work, we have proposed some unconditionally stable operator splitting schemes for solving stiff nonlinear differential equations. The error estimates and order of accuracy of the proposed schemes are given. For validation, some benchmark problems are simulated numerically and the results are illustrated graphically. It is observed that schemes G1 and H2 have competitive advantages in accuracy and can be used for solving the nonlinear and stiff partial differential equations.