Abstract
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. Therefore, 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. The computational results show that the standard discrete Laplace operator has severe grid orientation dependence. However, the isotropic discrete Laplace operator generates good results.
1. Introduction
We consider the Cahn–Hilliard (CH) equation with the reaction term:where is the order parameter in the domain , , , and is a constant. The reaction term 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 and .
The CH equation without a source term was derived for minimizing the interface in binary alloys [1, 2]. This topic has been studied in various fields such as multiphase flows [3–5], image inpainting [6], phase-field model [7–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–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–23].
However, the conventional space discretization scheme has a high probability of having grid-related artifact because it has anisotropic properties. Therefore, 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. The 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.
The main purpose of this work is to present benchmark problems for the numerical discretization of the CH equation with a source term.
The 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.
2. 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 discrete domain for , i.e., , where and are positive even integers, and is the grid size. Let and be approximations of and , respectively, where and 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]. The first scheme is the following with the standard discrete Laplace operator:where the standard discrete Laplace operator is defined as
The second scheme is the isotropic discretization of the Laplace operator [24]:where the isotropic discrete Laplace operator is defined as
Similarly, we can define both discrete Laplace operators in a discrete domain for , where , and are positive even integers, and is the grid size. The standard discrete Laplace operator is defined as
The isotropic discrete Laplace operator is defined as
We use the nonlinear multigrid method [28–30] to solve the above discrete CH equations.
3. Numerical Results and Discussion
3.1. Convergence Test
We start with numerical convergence tests of two discrete Laplace schemes. We consider the heat equation with the periodic boundary condition (3) on a domain :
For example, we use the following initial condition for equation (9):
Then, the exact solution is . Numerical convergence tests of two different schemes, standard and isotropic schemes, are performed with and , . The final time is for , and 6. Here, is used. Let be the numerical error in . We define the discrete -norm of the error as . Then, the rate of convergence is . Table 1 shows the discrete -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.
3.2. Comparison of the Two Schemes for an Isotropic Initial Shape
We consider the tumor growth simulation parameters in [12]:where and are the growth and death coefficients, respectively. We use and .
Figure 1 shows a plot of . Note that when and .
3.2.1. Case
We present the evolution of a disk by using the two types of schemes: standard and isotropic. The initial condition in the domain is given by
In this test, the uniform space step size , the time step size , and are used. Figure 2 shows the circular evolution of the disk with standard and isotropic schemes, respectively, up to the final time at every .
(a)
(b)
There 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 aswhere implies the perimeter of zero-level contour and implies the area enclosed by the contour. Note that if an evolution maintains the form of the disk. Figure 3 shows the ratios of the standard scheme and the isotropic scheme.
Thus, 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.
3.2.2. Case
We further present the evolution of sphere in by using both standard and isotropic schemes with the initial condition given by
Figure 4 shows isosurfaces of the evolution of equation (6) with both schemes. Here, we use the uniform step size , the time step , , and the final time .
(a)
(b)
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 such as each surface area over volume with appropriate scaling to compare evolutionary detail of both schemes, and hence is defined aswhere is a surface area of sphere and is a volume of sphere. Note that if an evolution maintains the form of the sphere. Figure 5 shows the ratios of both schemes. Thus, we conclude that the isotropic scheme is relatively unbiased and has a correct evolution.
3.3. 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 get a more detailed look at the grid effect, we now perform a benchmark test on an anisotropic initial shape.
3.3.1. Case
We investigate tumor growth of an elliptical shape in with and rotations in clockwise direction in , where initial conditions are given, respectively, as follows:where and are defined as
Here, , , , and are used. We show four stages of tumor growth: , , , and . 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. The final row shows overlapped contours for each column at , , , and . The results under the same conditions by using the isotropic scheme are given in Figure 7.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
3.3.2. Case
We further investigate the evolution of ellipsoidal shape in , where initial conditions are given, respectively, as follows:where , , and are defined as
Note that , , and are rotation matrices as follows:where and rotations are based on -, -, and -axes, respectively. We use the uniform step size , the time step size , , and the final time . Figure 8 shows evolutions of ellipsoid over time , and with the standard scheme. Snapshots of isosurface 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. The final row shows overlapped isosurfaces for each column at , , and , respectively. Note that we use rotation matrices to align the axes of case with those of , so we could overlap isosurfaces on identical axes. The 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. This is the reason why we should use the isotropic scheme in tumor growth simulation.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
4. Conclusions
We presented two benchmark problems for the numerical discretization of the CH equation with a source term. The first benchmark problem is the growth of a disk or a sphere. If the source term is isotropic, then the growth should be isotropic. The second benchmark problem is the growth of a rotated ellipse or a rotated ellipsoid. Ideally, the growth should be independent of rotation. Therefore, it is essential to satisfy these two benchmark problems if a numerical discretization is reliable and accurate.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
Junseok Kim was supported by the National Research Foundation (NRF), Korea, under project BK21 FOUR.