Benchmark Problems for the Numerical Discretization of the Cahn–Hilliard Equation with a Source Term

In this paper, we present benchmark problems for the numerical discretization of the Cahn–Hilliard equation with a source term. If the source term includes an isotropic growth term, then initially circular and spherical shapes should grow with their original shapes. However, there is numerical anisotropic error and this error results in anisotropic evolutions. )erefore, it is essential to use isotropic space discretization in the simulation of growth phenomenon such as tumor growth. To test numerical discretization, we present two benchmark problems: one is the growth of a disk or a sphere and the other is the growth of a rotated ellipse or a rotated ellipsoid. )e computational results show that the standard discrete Laplace operator has severe grid orientation dependence. However, the isotropic discrete Laplace operator generates good results.


Introduction
We consider the Cahn-Hilliard (CH) equation with the reaction term: where ϕ(x, t) is the order parameter in the domain Ω ⊂ R d , d � 2, 3, F(ϕ) � 0.25(ϕ 2 − 1) 2 , and ϵ is a constant. e reaction term c(ϕ) is specified in Section 3. In this paper, we present two benchmark problems for the numerical discretization of the CH equation with a source term in 2D and 3D. e CH equation without a source term was derived for minimizing the interface in binary alloys [1,2]. is topic has been studied in various fields such as multiphase flows [3][4][5], image inpainting [6], phase-field model [7][8][9][10], and microstructures with elastic inhomogeneity [11]. Since the common behavior of the CH equation without source terms is minimizing the interface, the anisotropic space discretization error is not noticeable. However, in the case of a growth model, anisotropic discretization may have critical problems. In particular, tumor growth simulation with the CH equation is one of the cases. Several research studies have been conducted numerically by using the Galerkin finite element method [12][13][14][15]. Fakih [16,17] analyzed the asymptotic behavior of the CH equation with source in both Dirichlet and Neumann boundary conditions by using the P1 element method. Khain and Sander [18] studied the logistic growth model of the CH equation. In recent years, many studies for tumor growth simulation using the finite difference method have been actively conducted [19][20][21][22][23].
However, the conventional space discretization scheme has a high probability of having grid-related artifact because it has anisotropic properties. erefore, we need to use the isotropic discretization. For instance, Kumar introduced the isotropic finite difference method [24]. We use the method proposed by Kumar in this paper. e isotropic finite difference scheme was derived for symmetric dendritic solidification and reduces an observable computational bias of the numerical anisotropy. In [25], the authors introduced different stencils for the Laplacian, bi-Laplacian, and gradient Laplacian operators and emphasized the advantages of using the isotropic stencils. Assadi [26] applied the isotropic scheme to the Allen-Cahn-type equation.
e main purpose of this work is to present benchmark problems for the numerical discretization of the CH equation with a source term. e rest of this article is organized as follows. In Section 2, we describe two discretization schemes of the CH equation with a source term. We compare two schemes on several numerical experiments in Section 3. In Section 4, we provide a summary and present our conclusions.

Discretization Schemes
We consider the two discretization schemes of the CH equation with a source term. To solve the CH equation numerically, we first define a 2D discrete domain for where N x and N y are positive even integers, and h � (b − a)/N x � (d − c)/N y is the grid size. Let ϕ n ij and μ n ij be approximations of ϕ(x i , y j , t n ) and μ(x i , y j , t n ), respectively, where t n � nΔt and Δt is the time step. We use the periodic boundary condition for ϕ and μ as follows: For simplicity, we take Eyre's nonlinearly stabilized splitting scheme [27]. e first scheme is the following with the standard discrete Laplace operator: where the standard discrete Laplace operator Δ S is defined as e second scheme is the isotropic discretization of the Laplace operator [24]: where the isotropic discrete Laplace operator Δ I is defined as Similarly, we can define both discrete Laplace operators e standard discrete Laplace operator is defined as e isotropic discrete Laplace operator is defined as Discrete Dynamics in Nature and Society We use the nonlinear multigrid method [28][29][30] to solve the above discrete CH equations.

Convergence Test.
We start with numerical convergence tests of two discrete Laplace schemes. We consider the 2D heat equation with the periodic boundary condition (3) on a domain Ω � (0, 1) × (0, 1): For example, we use the following initial condition for equation (9): en, the exact solution is u(x, y, t) � e − 8 π 2 t cos (2πx)cos(2πy). Numerical convergence tests of two different schemes, standard and isotropic schemes, are per- be the numerical error in Ω h . We define the discrete l 2 -norm of the error as en, the rate of convergence is log 2 (‖e h ‖ 2 /‖e h/2 ‖ 2 ). Table 1 shows the discrete l 2 -norm of errors and the rates of convergence with the two different formulas. Numerical experimentation suggests that both schemes are second-order accurate in space.

Comparison of the Two Schemes for an Isotropic Initial
Shape. We consider the tumor growth simulation parameters in [12]: where λ g and λ d are the growth and death coefficients, respectively. We use λ g � 280 and λ d � 92. Figure 1 shows a plot of c(ϕ) We present the evolution of a disk by using the two types of schemes: standard and isotropic. e initial condition in the domain Ω � (−0.6, 0.6) × (−0.6, 0.6) is given by In this test, the uniform space step size h � 1.2/128, the time step size Δt � 10h 2 , and ϵ � 0.0125 are used. Figure 2 shows the circular evolution of the disk with standard and isotropic schemes, respectively, up to the final time t � 1200Δt at every 120Δt.
ere is a significant computational bias for diagonal directions in the case with the standard scheme. On the other hand, the computational bias in the case with the isotropic scheme seems not large. In order to compare both schemes in more detail, we compare the ratios such as where l implies the perimeter of zero-level contour and A implies the area enclosed by the contour. Note that r d ≈ 1 if an evolution maintains the form of the disk. Figure 3 shows the ratios r d of the standard scheme and the isotropic scheme. us, we can confirm that an initially circular shape actually evolves closer to the original shape with the isotropic scheme while the standard scheme does not.

3D
Case. We further present the evolution of sphere in Ω � (−0.6, 0.6) 3 by using both standard and isotropic schemes with the initial condition given by Figure 4 shows isosurfaces ϕ � 0 of the evolution of equation (6) with both schemes. Here, we use the uniform step size h � 1.2/128, the time step Δt � 10h 2 , ϵ � 0.0125, and the final time t � 1600Δt.
Obviously, as shown in Figure 4(b), it was tested using the isotropic scheme that grew closer to the shape of the initial condition. For the next step, we use the ratios r s such as each surface area over volume with appropriate scaling to compare evolutionary detail of both schemes, and hence r s is defined as where S is a surface area of sphere and V is a volume of sphere. Note that r s ≈ 1 if an evolution maintains the form of the sphere. Figure 5 shows the ratios r s of both schemes. us, we conclude that the isotropic scheme is relatively unbiased and has a correct evolution.

Comparison of the Two Schemes for an Anisotropic Initial
Shape. In Section 3.2, we have performed evolutions of a sphere by using standard and isotropic schemes. In order to Discrete Dynamics in Nature and Society get a more detailed look at the grid effect, we now perform a benchmark test on an anisotropic initial shape.

2D
Case. We investigate tumor growth of an elliptical shape in 2D with 0°and 45°rotations in clockwise direction in Ω � (−2.4, 2.4) × (−2.4, 2.4), where initial conditions are given, respectively, as follows: where X and Y are defined as     (14) with the standard scheme at t � 0 and t � 1600Δt. (b) Overlapped isosurfaces ϕ � 0 for evolution in equation (14) with the isotropic scheme at t � 0 and t � 1600Δt.
Here, θ � −π/4, h � 4.8/512, Δt � h 2 , and ϵ � 0.0125 are used. We show four stages of tumor growth: t � 0, 15000Δt, 35000Δt, and 50000Δt. We can see a difference of tumor growth in two kinds of space discretization. Figure 6 shows evolutions of ellipse with the standard scheme. Snapshots of zero-level contours for evolution over time in equation (16) are listed in the first row of Figure 6, and those in equation (17) are in the second row of Figure 6. e final row shows overlapped contours for each column at t � 0, 15000Δt, 35000Δt, and 50000Δt. e results under the same conditions by using the isotropic scheme are given in Figure 7.

3D
Case. We further investigate the evolution of ellipsoidal shape in Ω � (−1.0, 1.0) 3 , where initial conditions are given, respectively, as follows: where X, Y, and Z are defined as Note that R x , R y , and R z are rotation matrices as follows: 45°0°( l) Figure 6: (a-d) Snapshots of zero-level contours for evolution in (16) using the standard scheme at t � 0, 15000Δt, 35000Δt, and 50000Δt, respectively. (e-h) Snapshots of zero-level contours for evolution in (17) using the standard scheme at t � 0, 15000Δt, 35000Δt, and 50000Δt, respectively. (i-l) Overlapped zero-level contours for each column at t � 0, 15000Δt, 35000Δt, and 50000Δt, respectively. 6 Discrete Dynamics in Nature and Society where θ � −π/4 and rotations are based on x-, y-, and z-axes, respectively. We use the uniform step size h � 2.0/128, the time step size Δt � 10h 2 , ϵ � 0.0125, and the final time 1000Δt. Figure 8 shows evolutions of ellipsoid over time t � 0, 600Δt, and 1000Δt with the standard scheme. Snapshots of isosurface ϕ � 0 for evolution over time in equation (9) are listed in the first row of Figure 8, and those in equation (10) are listed in the second row of Figure 8. e final row shows overlapped isosurfaces ϕ � 0 for each column at t � 0, 600Δt, and 1000Δt, respectively. Note that we use rotation matrices to align the axes of θ � −π/4 case with those of θ � 0, so we could overlap isosurfaces on identical axes. e results under the same conditions by using the isotropic scheme are listed in Figure 9. Unlike the isotropic scheme, there is a computational bias in the standard scheme. is is the reason why we should use the isotropic scheme in tumor growth simulation. Figure 7: (a-d) Snapshots of zero-level contours for evolution in equation (16) using the isotropic scheme at t � 0, 15000Δt, 35000Δt, and 50000Δt, respectively. (e-h) Snapshots of zero-level contours for evolution in equation (17) using the isotropic scheme at t � 0, 15000Δt, 35000Δt, and 50000Δt, respectively. (i-l) Overlapped zero-level contours for each column at t � 0, 15000Δt, 35000Δt, and 50000Δt, respectively.

Conclusions
We presented two benchmark problems for the numerical discretization of the CH equation with a source term. e first benchmark problem is the growth of a disk or a sphere. If the source term is isotropic, then the growth should be isotropic. e second benchmark problem is the growth of a rotated ellipse or a rotated ellipsoid. Ideally, the growth should be independent of rotation. erefore, it is essential to satisfy these two benchmark problems if a numerical discretization is reliable and accurate.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.