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. 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.
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 , phase-field model [7–10], and microstructures with elastic inhomogeneity . 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  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 . 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 , the authors introduced different stencils for the Laplacian, bi-Laplacian, and gradient Laplacian operators and emphasized the advantages of using the isotropic stencils. Assadi  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 . 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 :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 :where and are the growth and death coefficients, respectively. We use and .
Figure 1 shows a plot of . Note that when and .
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 .
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.
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 .
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.
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.
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.
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.
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.
Junseok Kim was supported by the National Research Foundation (NRF), Korea, under project BK21 FOUR.
D. Anders and K. Weinberg, “Numerical simulation of diffusion induced phase separation and coarsening in binary alloys,” Computational Materials Science, vol. 50, no. 4, pp. 1359–1364, 2011.View at: Publisher Site | Google Scholar
J. W. Cahn and J. E. Hilliard, “Free energy of a nonuniform system. I. Interfacial free energy,” The Journal of Chemical Physics, vol. 28, no. 2, pp. 258–267, 1958.View at: Publisher Site | Google Scholar
F. Guillén-González and G. Tierra, “Unconditionally energy stable numerical schemes for phase-field vesicle membrane model,” Journal of Computational Physics, vol. 354, pp. 67–85, 2018.View at: Publisher Site | Google Scholar
B. S. Hosseini, S. Turek, M. Möller, and C. Palmes, “Isogeometric Analysis of the Navier-Stokes-Cahn-Hilliard equations with application to incompressible two-phase flows,” Journal of Computational Physics, vol. 348, pp. 171–194, 2017.View at: Publisher Site | Google Scholar
Y. Wang, C. Shu, H. B. Huang, and C. J. Teo, “Multiphase lattice Boltzmann flux solver for incompressible multiphase flows with large density ratio,” Journal of Computational Physics, vol. 280, pp. 404–423, 2015.View at: Publisher Site | Google Scholar
A. L. Bertozzi, S. Esedoglu, and A. Gillette, “Inpainting of binary images using the Cahn-Hilliard equation,” IEEE Transactions on Image Processing, vol. 16, no. 1, pp. 285–291, 2007.View at: Publisher Site | Google Scholar
X. Hu, Y. Li, and H. Ji, “A nodal finite element approximation of a phase field model for shape and topology optimization,” Applied Mathematics and Computation, vol. 339, pp. 675–684, 2018.View at: Publisher Site | Google Scholar
C. Liu, F. Frank, and B. M. Rivière, “Numerical error analysis for nonsymmetric interior penalty discontinuous Galerkin method of Cahn-Hilliard equation,” Numerical Methods for Partial Differential Equations, vol. 35, no. 4, pp. 1509–1537, 2019.View at: Publisher Site | Google Scholar
Z. Mu, Y. Gong, W. Cai, and Y. Wang, “Efficient local energy dissipation preserving algorithms for the Cahn-Hilliard equation,” Journal of Computational Physics, vol. 374, pp. 654–667, 2018.View at: Publisher Site | Google Scholar
J. M. Park, “Mathematical modeling and computational simulation of phase separation in ternary mixtures,” Applied Mathematics and Computation, vol. 330, pp. 11–22, 2018.View at: Publisher Site | Google Scholar
J. Zhu, L.-Q. Chen, and J. Shen, “Morphological evolution during phase separation and coarsening with strong inhomogeneous elasticity,” Modelling and Simulation in Materials Science and Engineering, vol. 9, no. 6, pp. 499–511, 2001.View at: Publisher Site | Google Scholar
A. C. Aristotelous, O. A. Karakashian, and S. M. Wise, “Adaptive, second–order in time, primitive–variable discontinuous Galerkin schemes for a Cahn–Hilliard equation with a mass source,” IMA Journal of Numerical Analysis, vol. 35, no. 3, pp. 167–1198, 2015.View at: Publisher Site | Google Scholar
A. Hawkins-Daarud, K. G. van der Zee, and J. Tinsley Oden, “Numerical simulation of a thermodynamically consistent four-species tumor growth model,” International Journal for Numerical Methods in Biomedical Engineering, vol. 28, no. 1, pp. 3–24, 2012.View at: Publisher Site | Google Scholar
V. Mohammadi and M. Dehghan, “Simulation of the phase field Cahn-Hilliard and tumor growth models via a numerical scheme: element-free Galerkin method,” Computer Methods in Applied Mechanics and Engineering, vol. 345, pp. 919–950, 2019.View at: Publisher Site | Google Scholar
J. Ren, D. Shi, and S. Vong, “High accuracy error estimates of a Galerkin finite element method for nonlinear time fractional diffusion equation,” Numerical Methods for Partial Differential Equations, vol. 36, no. 2, pp. 284–301, 2020.View at: Publisher Site | Google Scholar
H. Fakih, “A Cahn–Hilliard equation with a proliferation term for biological and chemical applications,” Asymptotic Analysis, vol. 94, no. 1-2, pp. 71–104, 2015.View at: Publisher Site | Google Scholar
H. Fakih, “Asymptotic behavior of a generalized Cahn-Hilliard equation with a mass source,” Applicable Analysis, vol. 96, no. 2, pp. 324–348, 2017.View at: Publisher Site | Google Scholar
E. Khain and L. M. Sander, “Generalized Cahn–Hilliard equation for biological applications,” Physical Review E., vol. 77, no. 5, Article ID 051129, 2008.View at: Publisher Site | Google Scholar
J. J. Benito, A. García, L. Gavete, M. Negreanu, F. Ureña, and A. M. Vargas, “On the convergence of the generalized finite difference method for solving a chemotaxis system with no chemical diffusion,” Computational Particle Mechanics, vol. 8, no. 3, pp. 625–636, 2021.View at: Publisher Site | Google Scholar
J. J. Benito, Á. García, M. L. Gavete, M. Negreanu, F. Ureña, and A. M. Vargas, “Convergence and numerical solution of a model for tumor growth,” Mathematics, vol. 9, no. 12, p. 1355, 2021.View at: Publisher Site | Google Scholar
M. I. P. Hidayat, “Meshless finite difference method with B-splines for numerical solution of coupled advection-diffusion-reaction problems,” International Journal of Thermal Sciences, vol. 165, Article ID 106933, 2021.View at: Publisher Site | Google Scholar
V. Mohammadi, M. Dehghan, and S. De Marchi, “Numerical simulation of a prostate tumor growth model by the RBF-FD scheme and a semi-implicit time discretization,” Journal of Computational and Applied Mathematics, vol. 388, Article ID 113314, 2021.View at: Publisher Site | Google Scholar
M. Smoka, M. Woniak, and R. Schaefer, “A three-level linearized time integration scheme for tumor simulations with Cahn-Hilliard equations,” in International Conference on Computational Science, pp. 173–185, Springer, Cham, Switzerland, 2021.View at: Google Scholar
A. Kumar, “Isotropic finite-differences,” Journal of Computational Physics, vol. 201, no. 1, pp. 109–118, 2004.View at: Publisher Site | Google Scholar
M. Patra and M. Karttunen, “Stencils with isotropic discretization error for differential operators,” Numerical Methods for Partial Differential Equations, vol. 22, no. 4, pp. 936–953, 2006.View at: Publisher Site | Google Scholar
H. Assadi, “A phase-field model for non-equilibrium solidification of intermetallics,” Acta Materialia, vol. 55, no. 15, pp. 5225–5235, 2007.View at: Publisher Site | Google Scholar
D. J. Eyre, “An unconditionally stable one-step scheme for gradient systems,” 1997.View at: Google Scholar
A. C. Aristotelous, O. Karakashian, O. Karakashian, and S. M. Wise, “A mixed discontinuous Galerkin, convex splitting scheme for a modified Cahn-Hilliard equation and an efficient nonlinear multigrid solver,” Discrete & Continuous Dynamical Systems-B, vol. 18, no. 9, pp. 2211–2238, 2013.View at: Publisher Site | Google Scholar
W. L. Briggs and S. F. McCormick, A Multigrid Tutorial, SIAM, Philadelphia, PA, USA, 2nd edition, 2000.
U. Trottenberg, A. Schüller, and C. Oosterlee, Multigrid, Academic Press, New York, NY, USA, 1st edition, 2000.