An Explicit Hybrid Method for the Nonlocal Allen–Cahn Equation

: We extend the explicit hybrid numerical method for solving the Allen–Cahn (AC) equation to the scheme for the nonlocal AC equation with isotropically symmetric interfacial energy. The proposed method combines the previous explicit hybrid method with a space-time dependent Lagrange multiplier which enforces conservation of mass. We perform numerical tests for the area-preserving mean curvature ﬂow, which is the basic property of the nonlocal AC equation. The numerical results show good agreement with the theoretical solutions. Furthermore, to demonstrate the usefulness of the proposed method, we perform a cell growth simulation in a complex domain. Because the proposed numerical scheme is explicit, it is remarkably simple to implement the numerical solution algorithm on complex discrete domains.


Introduction
The phase-field model is one of the representatives of the interface capturing approach and has been widely investigated in order to interpret the interfacial dynamics [1]. The Allen-Cahn (AC) equation, which is a phenomenological model for antiphase domain coarsening in a binary mixture [2], is one of the most popular phase-field models and is given as where Ω is a bounded domain in R d (d = 1, 2, 3), t is time, φ(x, t) is an order parameter, F(φ) = 0.25(φ 2 − 1) 2 , and is a positive constant related to the interfacial energy. Here, we use the homogenous Neumann boundary condition, i.e., n · ∇φ(x, t) = 0 for x ∈ ∂Ω, where n is the unit outer normal vector on the domain boundary ∂Ω. The AC equation has been applied to phase transition, image processing, motion by mean curvature, multiphase flows, and dendritic growth (see e.g., [3][4][5][6][7][8] and the references therein). However, there are some mathematical problems that cannot be solved using the original form of the classical AC equation. For instance, the long range interactions in the interfacial dynamics is difficult to investigate using the original AC equation, therefore the AC equation with nonlocal diffusion (fractional) operator was analyzed in [9]. Meanwhile, the AC equation is not conservative, and Brassel and Bretin [10] proposed the nonlocal AC equation with the time-dependent Lagrange multiplier which preserves the shape of interface in local coordinates unlike the multiplier presented earlier in [11]. This equation has both nonlocal and local effects, even though the mass conservation property can be achieved in the local AC equation; hence we adopt the nonlocal 2 of 14 equation in [10] and propose an explicit hybrid method in this paper for the following nonlocal AC equation with isotropically symmetric interfacial energy: ∂φ(x, t) ∂t = − F (φ(x, t)) 2 + ∆φ(x, t) + β(t) F(φ(x, t)), x ∈ Ω, t > 0, (2) where β(t) = Ω F (φ(x, t)) dx/[ 2 Ω F(φ(x, t)) dx] to possess the total mass conservation property [12], i.e., Please note that if β(t) = 0, then Equation (2) becomes the classical AC Equation (1). The nonlocal AC equation with mass conservation has been investigated numerically and analytically [13][14][15][16][17][18][19]. In particular, the AC equation is of second order; hence the nonlocal AC model can partially replace the Cahn-Hilliard type model, which is of fourth order and conserves mass but is relatively costly to obtain numerical solutions. Therefore, the nonlocal AC equation can be used to apply various mathematical models. In this regard, there is an advantage of presence of the mass conservation property to implement mass conservative models such as the lattice Boltzmann (LB) method [13]. Proposing two simple LB methods, the authors in [13] conducted comparative studies of both the local and nonlocal AC equations with mass conservation. Moreover, Zhang et al. [14] proposed an unconditionally energy stable second-order non-iterative method to the anisotropic AC equation based on scalar auxiliary variable approach. The authors treated the nonlinear terms as a source, this linear stabilization provides stability for a large time step and hence computational efficiency. Guan et al. [15] proposed the concrete convergence analysis about the discrete second-order convex splitting scheme for the nonlocal AC equation. One of the widely used numerical methods to solve the AC equation with mass conservation is the Fourier spectral method. The authors in [16] employed the Fourier spectral method for nonlocal diffusion operators. Using the Fourier spectral approximation, a hybrid algorithm was presented and applied to the nonlocal AC equation. Furthermore, a mass conserving high-order method for solving the nonlocal AC equation was developed in [17]. The author employed the spectral method in space and three stage multiple-order semi-implicit Runge-Kutta method in time, and emphasized that the existing methods did not actually achieve mass conservation property by numerical simulations. The authors in [18] presented a second-order fast explicit operator splitting (FEOS) method based on the Fourier spectral method. The authors further provided the convergence analysis, investigated the discrete maximum principle of the proposed scheme, and presented various numerical experiments. Zhai et al. [19] solved the fractional nonlocal AC equation numerically using the FEOS method with adaptive time-stepping algorithm. It was taken the advantage of being able to adjust time step to multiscale temporal characteristic of the phase separation and coarsening process. The nonlocal AC model can be applied to multiphase models coupled with incompressible Navier-Stokes (NS) equations [20][21][22]. Specifically, under the framework of the LB method, Ren et al. [20] incorporated the incompressible hydrodynamic equations with the conservative AC equation and provided an effective solution for binary flow modeling, which had been difficult due to the interface limitations and numerical dispersion. Aihara et al. [21] developed a new multiphase method using the conservative AC equation and showed the accurate evaluation in the movement of bubbles interacting with the liquid-liquid interface. Furthermore, Joshi and Jaiman [22] proposed a nonlinear adaptive variational method to solve the coupled AC and NS equations for fluid-fluid phase flows, which is formulated by the finite element approach. The proposed algorithm was proved to be energy stable and nonlocal mass conserving through the spinodal decomposition.
While much of the literature dealt with implicit schemes, the main purpose of this article is to propose an explicit hybrid numerical method for Equation (2). Basically, it is easy to implement and can be applied directly to various types of domain, thus it consumes relatively less computational cost than implicit methods unless the long term simulation is conducted. Furthermore, the advantages to employ the proposed explicit hybrid method can be maximized since the AC model is of second order as mentioned above.
The article is organized as follows. In Section 2, we describe a numerical algorithm using an operator splitting method. Several numerical results that demonstrate the accuracy and usefulness of the proposed numerical method are presented in Section 3. Conclusions are made in Section 4.

Numerical Algorithm
We present the explicit hybrid numerical solution algorithm for solving the nonlocal AC equation, which is an extension of the algorithm for the AC equation [4]. Let us consider a computational domain where N x and N y are the numbers of cells in xand ydirections, respectively. We denote by φ n ij the numerical approximations of φ( . . , N x and j = 1, . . . , N y and ∆t is the time step. Let us rewrite Equation (2) before applying an operator splitting method as follows: The explicit hybrid method consists of three steps: Step 1.

Numerical Results
In this section, we perform numerical experiments such as the temporal evolution of three disjoint disks and a growing cell inside a complex domain. For all tests, we , where qh is the width of the transition layer [24].

Evolution of Disks
The evolution equations for the radii of spheres in d-dimensional geometric flows are given in [25]. Let R i (t) be the radius of the ith-sphere at time t for i = 1, 2, . . . , m. Then, the evolution equations are given by where κ i is the sum of principal curvatures of interface Γ i and |Γ i | is the perimeter. In this test, we consider three disjoint disks in two-dimensional space, i.e., d = 2 and m = 3. From Equation (10), we have Let us assume 0 . Applying the explicit Euler method to Equations (11)-(13) with a temporal step size δt, we get For simplicity of notation, let R k i = R i (kδt), (i = 1, 2, 3). To investigate the relationship among R 1 , R 2 , and R 3 , we first consider By mathematical induction, if R n 1 > 0, then it can be shown that Let us show Equation (18) is valid for n = 1.
where we used the arithmetic-geometric-harmonic mean (AGHM). Next, assume that Equation (18) is true for n = k. Using the similar process in Equation (19), we get where 2R k 2 > R k 1 + R k 3 is used. Furthermore, we get the other inequality as follows: with the similar process by using the AGHM. Therefore, Equation (18) holds for n = k + 1. That is, (17) is satisfied. As δt approaches 0, Equation (18) also holds for the analytical solution of Equations (11)-(13) as follows: To show Equation (18) numerically, we perform a test. In this test, we compare the result between the numerical solution of Equation (2) and the reference solution of Equations (11)- (13). We use a temporal step size ∆t = 0.1h 2 on Ω = (0, 1.2) × (0, 1.2) with a mesh grid 256 × 256 and = 6 and set the initial condition for three disjoint disks as follows: with R 0 1 = 0.175, R 0 2 = 0.205, and R 0 3 = 0.225. For the reference solutions of R 1 (t), R 2 (t), and R 3 (t), we solve Equations (11)-(13) numerically by adopting the explicit Euler method with a small temporal step size δt = 0.01h 2 . We note that a high-order method such as fourth-order Runge-Kutta method [26] can be used as a solver for ordinary differential systems. Figure 1a shows evolutions of the radii of the three distinct disks at t = 0, 17,500∆t, and 22,100∆t. The disk with radius R 1 (t) becomes smaller as time goes on, while the others get larger. Figure 1b shows that the reference solution and the numerical solutions of three distinct radii up to t = 62,500∆t. As shown in Figure 1b, R 1 (t) disappears at about t = 22,500∆t and R 2 (t) increases until t = 22,500∆t. Subsequently, R 2 (t) is strictly decreasing. Another initial condition is Using the same process of Equations (19)- (22), we obtain We take the numerical test for Equation (25) with the initial condition (24) where R 0 1 = 0.175, R 0 2 = 0.2, and R 0 3 = 0.225 under the same conditions above. Figure 2a shows evolutions of the radii of three distinct disks at t = 0, 19,000∆t, and 22,500∆t. The evolution of disks is similar to Figure 1a. Figure 2b shows that the reference solution and the numerical solutions of three distinct radii up to time t = 55,000∆t. As shown in Figure 2b, R 1 (t) disappears at about t = 25,000∆t and R 2 (t) increases until t = 25,000∆t. On the other hand, the initial condition results in the opposite direction. That is, for all n ≥ 1. We conduct the numerical test for Equation (27) with the initial condition (24) where R 0 1 = 0.175, R 0 2 = 0.18, and R 0 3 = 0.225 under the same conditions above. Figure 3a shows evolutions of the radii of three distinct disks at t = 0, 13,000∆t, and 26,000∆t. The disks with radius R 1 (t) and R 2 (t) become smaller as time goes on, while the other gets larger. In Figure 3b, the reference solution and the numerical solution of three distinct radii are shown up to time t = 33,000∆t. The radius R 3 (t) grows monotonically with our numerical scheme while R 1 (t) and R 2 (t) decrease as time goes on.

Comparison Test with a Conventional Method
We present a comparison test with an implicit hybrid scheme [27]: ).
Here, β n+1,2 is the same as Equation (9). The authors in [27] used the multigrid method to solve Equation (29) implicitly. We use a tolerance, tol = 10 −10 , for the V-cycle convergence in the multigrid method. Generally, the implicit methods can use large time steps. However, when the time step is large, then it is not accurate. We consider the temporal evolution of two circles with different radii on Ω = (0, 2) × (0, 2) with a mesh grid 128 × 128. The reference time step size ∆t ref = 0.1h 2 and = 6 are used. The initial condition is given by Figure 4a shows the evolution of two disks using the proposed method. We use the explicit Euler method with a small temporal step size δt = 0.01h 2 to find the exact evolution of radii. The explicit hybrid scheme and implicit hybrid scheme using time step ∆t = ∆t ref show highly accurate results compared to the exact radii. Next, we measure the CPU times to compare two schemes, the implicit hybrid scheme and the proposed scheme. The simulations are performed on Intel Core i5-6400 CPU @ 2.70 GHz processor and 4 GM RAM. The CPU time by using the proposed scheme with ∆t = ∆t ref and 10, 000 time steps is approximately 43.946 s. Table 1  However, the results are less accurate than the explicit scheme as shown in Table 1. Therefore, the results suggest that the proposed method is more accurate than the implicit scheme under similar computational cost as shown in Figure 4b in the case of ∆t = 4∆t ref .

Cell Growth in a Complex Domain
In this section, we consider a growing cell inside a non-rectangular domain. To simulate cell growth, we set the time-dependent parameter β(t) in Equation (2) as constant β. That is, where Ω in is a complex domain. Figure 5 illustrates  We present a numerical scheme for Equation (30) in a complex domain with the Dirichlet boundary condition, i.e., φ n ij = −1 for (x i , y j ) ∈ Ω h out . We solve Equation (30) only in Ω h in by using the operator splitting method (6)- (8). Here, we replace Equation (8) by Figure 6a-c shows the temporal evolutions of the initial disk in a complex domain which is star-shaped at t = 0, 800∆t, and 1100∆t. The shape of the domain is created with the image in Figure 6d. The initial condition on the embedding domain Ω = (0, 1) × (0, 1) is defined by Here, we use N x = N y = 125, h = 1/124, ∆t = 0.1h 2 , = 12 , and β = 2h to make the cell grow. In Figure 7a-c, we illustrate the temporal evolutions in a drop-shaped domain which we created with Figure 7d at t = 0, 550∆t, and 750∆t. We use N x = 92, N y = 124, h = 1/(N x − 1), ∆t = 0.1h 2 , = 12 , and β = 2h to make the cell grow. The initial condition on the embedding domain Ω = (0, (N x − 1)h) × (0, (N y − 1)h) is defined by We can obtain a similar computational result with a much simpler proposed scheme compared to the previous complex method [29]. Moreover, while there is a restriction on the ratio of domain sizes when we use a multigrid method, we do not have that restriction in the proposed method.

Conclusions
We proposed the explicit hybrid numerical method for the nonlocal AC equation with isotropically symmetric interfacial energy in this paper. In general, an explicit time-stepping scheme is not suitable for solving large systems because it requires strong time step restrictions and this further causes a stability problem. Because the AC model is of the second order, however, this gives a rationale for solving such drawbacks. For example, the nonlocal AC equation can be applied to many fields by replacing the Cahn-Hillard equation which is a fourth-order equation [30]. The AC model does not take much time to evaluate numerical solutions; hence numerical solutions can be obtained at a relatively adequate time as depicted in Table 1 even though one uses the explicit hybrid scheme. Accordingly, robustness and accuracy can be also achieved simultaneously when a sufficiently small size of time step is applied. Numerical tests are presented to demonstrate the basic property of the nonlocal AC equation and usefulness of the proposed method. Moreover, we performed numerical experiments for the area-preserving mean curvature flow. There have been good agreements between the theoretical solutions and the computational results. Compared to the fully implicit scheme, the proposed method is more accurate with similar elapsed time though the implicit scheme is still suitable for long time simulation. In the last section, numerical simulations confirm that the proposed method is not affected by the shapes of domain. Considering that it is common to employ square grids when an implicit numerical scheme is used for complex domains, the proposed algorithm is much simpler and faster. This is the reason the proposed method is preferred over the existing methods, and there are some direct applications of it such as smoothing the surface with volume objects [31], direct computation in narrow band corresponding to surface [32], etc. Therefore, the proposed method can be used efficiently with little input cost and get accurate numerical solutions.