Nonlinear Multigrid Implementation for the Two-Dimensional Cahn–Hilliard Equation

We present a nonlinear multigrid implementation for the two-dimensional Cahn–Hilliard (CH) equation and conduct detailed numerical tests to explore the performance of the multigrid method for the CH equation. The CH equation was originally developed by Cahn and Hilliard to model phase separation phenomena. The CH equation has been used to model many interface-related problems, such as the spinodal decomposition of a binary alloy mixture, inpainting of binary images, microphase separation of diblock copolymers, microstructures with elastic inhomogeneity, two-phase binary fluids, in silico tumor growth simulation and structural topology optimization. The CH equation is discretized by using Eyre’s unconditionally gradient stable scheme. The system of discrete equations is solved using an iterative method such as a nonlinear multigrid approach, which is one of the most efficient iterative methods for solving partial differential equations. Characteristic numerical experiments are conducted to demonstrate the efficiency and accuracy of the multigrid method for the CH equation. In the Appendix, we provide C code for implementing the nonlinear multigrid method for the two-dimensional CH equation.

Here, n is the unit normal vector on the domain boundary ∂Ω. The first boundary condition (1) implies that the interface contacts the domain boundary at a 90 • angle. The second boundary condition (2) implies that the total mass is conserved.
We can derive the CH equation from the following total free energy functional Taking the variational derivative of E with respect to φ, we define the chemical potential: Conservation of mass implies the following CH equation where the flux is given by F = −M∇µ. If we differentiate E (φ) and Ω φ dx with respect to time t, then we have and which imply that the total energy is decreasing and that the total mass is conserved in time, respectively. The CH equation was originally developed by Cahn and Hilliard to model spinodal decomposition in a binary alloy. The CH equation has been used to address many major problems such as the spinodal decomposition of a binary alloy mixture [3,4], inpainting of binary images [5,6], microphase separation of diblock copolymers [7,8], microstructures with elastic inhomogeneity [9,10], two-phase binary fluids [11,12], tumor growth models [13][14][15] and structural topology optimization [14,16]. Further details regarding the basic principles and practical applications of the CH equation are available in a recent review [14]. Thus, knowing how to implement a discrete scheme for the CH equation in detail is very useful because this equation is a building-block equation for many applications. The CH equation is discretized by using Eyre's unconditionally gradient stable scheme [17] and is solved by using a nonlinear multigrid technique [1], which is one of the most efficient iterative methods for solving partial differential equations. Several studies have used the nonlinear multigrid method for the CH-type equations [18][19][20][21][22][23]. However, details regarding the implementation, multigrid performance, and source codes have not been provided. Therefore, the main purpose of this paper is to describe a detailed multigrid implementation of the two-dimensional CH equation, evaluate its performance and provide its C programming language source code.
The remainder of this paper is organized as follows. In Section 2, we describe the numerical solution in detail. In Section 3, we describe the characteristic numerical experiments that are conducted to demonstrate the accuracy and efficiency of the multigrid method for the CH equation. In Section 4, we provide a conclusion. In the Appendix A, we provide the C code for implementing the nonlinear multigrid technique for the two-dimensional CH equation.

Numerical Solution
We consider a finite difference approximation for the CH equation. An unconditionally gradient energy stable scheme, which was proposed by Eyre, is applied to the time discretization. A nonlinear multigrid technique [1] is applied to solve the resulting system at an implicit time level.

Discretization
We discretize the CH equation in the two-dimensional space Ω = (a, b) × (c, d). Let N x = 2 p and N y = 2 q be the numbers of mesh points with integers p and q. Let ∆x = (b − a)/N x and ∆y = (d − c)/N y be the mesh size. Let ij and µ n ij be approximations of φ(x i , y j , t n ) and µ(x i , y j , t n ), respectively. Here, t n = n∆t and ∆t represent the temporal step. We assume a uniform mesh grid h = ∆x = ∆y and a constant mobility M = 1. Using the nonlinear stabilized splitting scheme of Eyre's unconditionally gradient stable scheme, the CH equation is discretized as where the discrete Laplace operator is defined by The homogeneous Neumann boundary conditions (1) and (2) are discretized as We define the discrete residual as For each element a ij of size N x × N y in matrix A, we define the Frobenius norm with a scaling and infinite norm as respectively. The discretizations (5) and (6) are conservative, that is, To show this conservation property, we take the summation of Equation (5) Here, we used the homogenous Neumann boundary conditions (7) and (8). Therefore, Equation (11) holds. We define the discrete energy functional as where we used the homogenous Neumann boundary conditions (7) and (8). We also define the discrete total mass as Then, the unconditionally gradient stable scheme satisfies the reduction in the discrete total energy [24]: which implies the pointwise boundedness of the numerical solution: The proof of Equation (16) can be found in Reference [25]. We provide the proof herein for the sake of completeness. We show that a constant K exists for all n values that satisfy the following inequality: Let us assume that there is an integer n K that is dependent on K such that φ n K ∞ > K for any K. Then, φ n K ij exists such that |φ n K ij | > K. Let K be the largest solution of E h (φ 0 ) = h 2 F(K), that is, where we utilize the fact that the total energy is decreasing and F(φ) is a strictly increasing function on (K, ∞). Equation (18) leads to a contradiction. Therefore, Equation (17) should be satisfied.

Multigrid V-Cycle Algorithm
We use the nonlinear full approximation storage (FAS) multigrid method to solve the nonlinear discrete systems (5) and (6). For simplicity, we define the discrete domains, Ω 2 , Ω 1 , and Ω 0 , which represent a hierarchy of meshes (Ω 2 , Ω 1 , and Ω 0 ) created by successively coarsening the original mesh Ω 2 , as shown in Figure 2. We summarize here the nonlinear multigrid method for solving the discrete CH system as follows: First, let us rewrite Equations (5) and (6) as where the linear operator NSO is defined as and the source term is denoted by Next, we describe the multigrid method, which includes the pre-smoothing, coarse grid correction and post-smoothing steps. We denote a mesh grid Ω k as the discrete domain for each multigrid level k. Note that a mesh grid Ω k contains 2 k × 2 k grid points. Let k min be the coarsest multigrid level. We now introduce the SMOOTH and V-cycle functions. Given the number ν 1 of pre-smoothing and ν 2 of post-smoothing relaxation sweeps, the V-cycle is used as an iteration step in the multigrid method.
(1) Pre-smoothing which represents ν 1 smoothing steps with the initial approximations φ m k , µ m k , source terms ξ n k , ψ n k and the SMOOTH relaxation operator to obtain the approximationsφ m k ,μ m k . One SMOOTH relaxation operator step consists of solving the systems (22) and (23), given as follows by 2 × 2 matrix inversion for each i and j. Here, we derive the smoothing operator in two dimensions. Rewriting Equation (5), we obtain: After substituting of this into (6), we obtain Next, we replace φ n+1 α,β and µ n+1 α,β in Equations (20) and (21) withφ m α,β andμ m α,β for α ≤ i and β ≤ j, otherwise with φ m α,β and µ m α,β , that is, (2) Compute the defect The restriction operator I k−1 k maps k-level functions to (k − 1)-level functions.
(4) Compute the right-hand side If k = 1, we explicitly invert the 2 × 2 matrix to obtain the solution. If k > 1, we solve Equation (24) by performing a FAS k-grid cycle using {φ m k−1 ,μ m k−1 } as the initial approximation: Here, the coarse values are simply transferred to the four nearby fine grid points, that is,

Further Numerical Schemes for the CH Equation
Previous studies have described the numerical solution of the CH equation with a variable mobility [19], the adaptive mesh refinement technique [26,27], the Neumann boundary condition in complex domains [20], the Dirichlet boundary conditions in complex domains [28], contact angle boundary [29], parallel multigrid method [30] and fourth-order compact scheme [31].
All computational simulations described in this section are performed on an Intel Core i5-6400 CPU @ 2.70 GHz with 4 GB of RAM.

Phase Separation
For the first numerical test, we consider spinodal decomposition in binary alloys. This decomposition is a process by which a mixture of binary materials separates into distinct regions with different material concentrations [2]. Figure 5a-c show snapshots of the phase-field φ at t = 100∆t, 200∆t and 1000∆t, respectively. The initial condition is φ(x, y, 0) = 0.1(1 − 2rand(x, y)) on Ω = (0, 1) × (0, 1), where rand(x, y) is a random value between 0 and 1. The parameters = 4 , h = 1/64, ∆t = 0.1h 2 and a tolerance of tol = 1.0 × 10 −10 are used.  Figure 6 shows the time evolution of the normalized discrete total energy

Non-Increase in Discrete Energy and Mass Conservation
where rand(x, y) is a random value between 0 and 1. We use the simulation parameters, = 4 , h = 1/64, ∆t = 0.1h 2 and tol = 1.0 × 10 −10 . The energy is non-increasing and the average concentration is conserved. These numerical results agree well with the total energy dissipation property (3) and the conservation property (4). The inscribed small figures are the concentration fields at the indicated times.
of the numerical solutions with the initial state (25).

Convergence Test
We consider the convergence of the Frobenius norm with a scaling of the residual error with respect to the grid size. The initial condition on the domain Ω = (0, 1) × (0, 1) is given as We fix = 0.06, ∆t = 0.01 and tol = 1.0 × 10 −15 . Here, we use the V(2, 2) scheme with a Gauss-Seidel relaxation, where (2, 2) indicates 2 pre-and 2 post-correction relaxation sweeps. We define the residual after m V-cycles as Table 1 shows the residual norm r m F after each V-cycle. Because no closed-form analytical solution exists for this problem, we define the Frobenius norm with a scaling of the residual error The grid sizes are set as 32 × 32, 64 × 64 and 128 × 128. The error norms and ratios of residual between successive V-cycle are shown in Table 1. As we have expected, the residual error decreases successively along with the V-cycle. The sharp increase in the residual norm ratio during the last few cycles reflects the fact that the numerical approximation is already accurate to near machine precision.

Effects of the Smooth Relaxation Numbers ν 1 and ν 2
We investigate the effects of the SMOOTH relaxation numbers ν 1 (pre-relaxation) and ν 2 (post-relaxation) on the CPU time. In this test, we perform a numerical simulation with the initial condition φ(x, y, 0) = 0.1 cos(πx) cos(πy) on Ω = (0, 1) × (0, 1), h = 1/128, 4 , ∆t = 0.1h 2 and tol = 1.0 × 10 −10 . Table 2 lists the average CPU times and average numbers of V-cycles for various pre-and post-relaxation numbers after 100 time steps. The relaxation numbers are rounded off to the nearest integer. Figure 8 shows the average CPU times for different pre-and post-relaxation numbers. We observe that the average CPU time is the lowest when the numbers of pre-and post-relaxation iterations are ν 1 = 2 and ν 2 = 4, respectively. Table 2. Average CPU times and average numbers of V-cycles (given in parentheses) for various preand post-relaxation numbers after 100 time steps. Next, we investigate the effect of SMOOTH relaxation numbers on the finest multigrid level. We perform a numerical simulation with = 0.0038 and ∆t = 1.0 × 10 −7 . The SMOOTH relaxation numbers ν 1 0 and ν 2 0 (on the finest multigrid level) are taken to be 1, 2, 3, 4, 5, 6, 7, 8, 9 and 10. In addition, other multigrid levels ν 1 k = ν 2 k are 1, 2, and 3 in the 128 × 128 and 512 × 512 mesh sizes. The other parameters are the same as those used previously. Table 3 shows the variations in average CPU time for different relaxation numbers with a 128 × 128 mesh size. Figure 9 illustrates the results in Table 3. Table 3. Average CPU times for various relaxation numbers on finest multigrid level (ν 1 0 , ν 2 0 ) after 10 time steps. In other grids, ν 1 k = ν 2 k , (1 ≤ k) relaxation number is fixed at 1, 2 and 3.  Figure 9. Average CPU times for various relaxation numbers on finest multigrid level (ν 1 0 , ν 2 0 ) and fixed relaxation number 1, 2 and 3 on other grids with a 128 × 128 mesh size. Table 4 lists the variations in average CPU time for different relaxation numbers with a 512 × 512 mesh size. Figure 10 illustrates the results in Table 4. Table 4. Average CPU times for various relaxation numbers on finest multigrid level (ν 1 0 , ν 2 0 ) after 10 time steps. In other grids, the ν 1 k = ν 2 k , (1 ≤ k) relaxation number is fixed at 1, 2 and 3.  Figure 10. Average CPU times for various relaxation numbers on finest multigrid level (ν 1 0 , ν 2 0 ) and fixed relaxation numbers 1, 2 and 3 on other grids with a 512 × 512 mesh size.

Effect Of V-Cycle
Next, we investigate the effect of V-cycle by changing the multigrid levels. In this test, we use the parameters ∆t = 0.01, = 0.06, SMOOTH relaxation = 2, tol = 1.0 × 10 −10 and h = 1/128 on Ω = (0, 1) × (0, 1) with a initial condition φ(x, y, 0) = 0.1 cos(πx) cos(πy). The highest number of V-cycles is taken to be 10000. We use level 2 , level 3 , level 4 , level 5 , level 6 and level 7 in a single time step as examples to illustrate the effect of the V-cycle. We calculate the CPU time for each level after 100∆t, as listed in Table 5. The number of the multigrid level and CPU time shown in Figure 11 indicate that a greater number of the multigrid level leads to a obvious decrease in CPU time. It is important to select an appropriate multigrid level for a specific mesh size.

Comparison between Gauss-Seidel and Multigrid Algorithms
We compare the average CPU times required to perform the Gauss-Seidel algorithm and multigrid algorithm. In this test, the initial condition is taken to be φ(x, y, 0) = cos(πx) cos(πy). The following parameters are used-∆t = 0.01, T = 10∆t, the SMOOTH relaxation = 2 and = 0.06. The highest number of V-cycle is taken to be 10000. The mesh sizes are 32 × 32, 64 × 64 and 128 × 128. The tolerances are 1.0 × 10 −3 , 1.0 × 10 −4 and 1.0 × 10 −5 . Table 6 shows the average CPU times for these two methods after 10 time steps. We observe that the multigrid method require less CPU time than the Gauss-Seidel method does.

Effects of tol and ∆t on the V-Cycle
In this test, we study the effects of tol and ∆t on the V-cycle with the initial condition being φ(x, y, 0) = 0.1 cos(πx) cos(πy), h = 1/128, = 0.06. The highest number of the V-cycle is taken to be 10,000. Table 7 shows the number of V-cycle for various tol and ∆t after a single time step. We can find that lower values of tol lead to an increase in the V-cycle. For different values of tol, it is essential to choose an appropriate ∆t to reduce the number of V-cycle.

Comparison of the Jacobi, Red-Black and Gauss-Seidel
We compare the performance of three relaxation methods: Jacobi, Red-Black and Gauss-Seidel. The initial condition is φ(x, y, 0) = 0.1 cos(πx) cos(πy) on Ω = (0, 1) × (0, 1). The parameters are h = 1/128, = 0.06, ∆t = 10 −7 , T = 100∆t and tol = 1.0 × 10 −10 . The SMOOTH relaxation numbers on the finest multigrid level (i.e., ν 1 = ν 2 ) are taken to be from 1 to 5 and those on the other multigrid levels are selected as 2 with the 128 × 128 mesh size. The relaxation numbers are rounded off to the nearest integer. Table 8 shows the average number of V-cycles for different relaxation numbers with the three methods. The relationship between the average numbers of V-cycles and ν 1 = ν 2 with the Jacobi, Red-Black and Gauss-Seidel method is plotted in Figure 12. The Gauss-Seidel method is observed to be the fastest. In the parallel multigrid method, the relaxation options are either Jacobi or Red-Black [33]. The Jacobi method requires approximately twice as many V-cycles as the Red-Black method does. Table 8. Average numbers of V-cycles for various relaxation numbers. The SMOOTH relaxation numbers on the finest multigrid level (i.e., ν 1 = ν 2 ) are taken to be from 1 to 5.

Effect of
Next, we investigate the effect of = m , which is related to the interface width. In this test, we perform a numerical simulation with the initial condition φ(x, y, 0) = 1 if 0.15 ≤ x ≤ 0.85 and 0.15 ≤ y ≤ 0.85, −1 otherwise, on Ω = (0, 1) × (0, 1). We use h = 1/128, ∆t = h, SMOOTH relaxation = 2, tol = 1.0 × 10 −10 and T = 1000∆t. Figure 13 presents the evolution of the CH equation with the three values 4 , 8 and 16 . As we have expected, the lower value of leads to a narrower interface width.

Effect of mesh size, N x × N y
In this test, we compare the CPU times with different mesh sizes N x × N y . The initial condition is φ(x, y, 0) = 0.1 cos(πx) cos(πy) on Ω = (0, N x /32) × (0, N y /32). The parameters are h = 1/32, ∆t = h, T = 100∆t, = 0.06, SMOOTH relaxation = 2 and tol = 1.0 × 10 −10 . Table 9 shows the CPU times and their ratios (that is, the ratio of the CPU time with the mesh size 2N x × 2N y to the CPU time with N x × N y ). We observe that the values converge to 4.

Conclusions
In this paper, we presented a nonlinear multigrid implementation for the CH equation in a two-dimensional space. Eyre's unconditionally gradient stable scheme was used to discretize the governing equation. The resulting discretizing equations were solved using the nonlinear multigrid method. We described the implementation of our numerical scheme in detail. We numerically showed the decrease in discrete total energy and the convergence of discrete total mass. We took a convergence test by studying the reductions in residual error on various mesh sizes in a single time step. The results of various numerical experiments were presented to demonstrate the effects of tolerance, SMOOTH relaxation, V-cycle and . The provided multigrid source code will be useful to beginners who needs the numerical implementation of the nonlinear multigrid method for the CH equation.
Author Contributions: All authors, C.L., D.J., J.Y., and J.K., contributed equally to this work and critically reviewed the manuscript. All authors have read and agreed to the published version of the manuscript. Acknowledgments: The authors thank the editor and the reviewers for their constructive and helpful comments on the revision of this article.

Conflicts of Interest:
The authors declare no conflicts of interest.