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

Authors: Sungha Yoon; Hyun Geun Lee; Yibao Li; Chaeyoung Lee; Jintae Park; Sangkwon Kim; Hyundong Kim; Junseok Kim
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## 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.


---

## Body

## 1. Introduction

We consider the Cahn–Hilliard (CH) equation with the reaction term:(1)∂ϕx,t∂t=Δμx,t+γϕx,t,x∈Ω,t>0,μx,t=F′ϕx,t−ϵ2Δϕx,t,where ϕx,t is the order parameter in the domain Ω⊂ℝd, d=2,3, Fϕ=0.25ϕ2−12, 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 2D and 3D.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 Section2, 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 a2D discrete domain for Ω=a,b×c,d, i.e., Ωh=xi,yj:xi=a+i−0.5h,yj=c+j−0.5h,1≤i≤Nx,1≤j≤Ny, where Nx and Ny are positive even integers, and h=b−a/Nx=d−c/Ny is the grid size. Let ϕijn and μijn be approximations of ϕxi,yj,tn and μxi,yj,tn, respectively, where tn=nΔt and Δt is the time step. We use the periodic boundary condition for ϕ and μ as follows:(2)ϕi,0=ϕi,Ny,ϕi,Ny+1=ϕi,1,ϕ0,j=ϕNx,j,ϕNx+1,j=ϕ1,j,μi,0=μi,Ny,μi,Ny+1=μi,1,μ0,j=μNx,j,μNx+1,j=μ1,j.For simplicity, we take Eyre’s nonlinearly stabilized splitting scheme [27]. The first scheme is the following with the standard discrete Laplace operator:(3)ϕijn+1−ϕijnΔt=ΔSμijn+1+γϕijn,μijn+1=ϕijn+13−ϕijn−ϵ2ΔSϕijn+1,where the standard discrete Laplace operator ΔS is defined as(4)ΔSϕij=ϕi−1,j+ϕi+1,j+ϕi,j−1+ϕi,j+1−4ϕijh2.The second scheme is the isotropic discretization of the Laplace operator [24]:(5)ϕijn+1−ϕijnΔt=ΔIμijn+1+γϕijn,μijn+1=ϕijn+13−ϕijn−ϵ2ΔIϕijn+1,where the isotropic discrete Laplace operator ΔI is defined as(6)ΔIϕij=4ϕi+1,j+ϕi−1,j+ϕi,j+1+ϕi,j−1+ϕi+1,j+1+ϕi−1,j+1+ϕi−1,j−1+ϕi+1,j−1−20ϕij6h2.Similarly, we can define both discrete Laplace operators in a3D discrete domain Ωh=xi,yj,zk:xi=a+i−0.5h,yj=c+j−0.5h,zk=e+k−0.5h,1≤i≤Nx,1≤j≤Ny,1≤k≤Nz for Ω=a,b×c,d×e,f, where Nx,Ny, and Nz are positive even integers, and h=b−a/Nx=d−c/Ny=f−e/Nz is the grid size. The standard discrete Laplace operator is defined as(7)ΔSϕijk=ϕi−1,j,k+ϕi+1,j,k+ϕi,j−1,k+ϕi,j+1,k+ϕi,j,k−1+ϕi,j,k+1−6ϕijkh2.The isotropic discrete Laplace operator is defined as(8)ΔIϕijk=20ϕi−1,j,k+ϕi+1,j,k+ϕi,j−1,k+ϕi,j+1,k+ϕi,j,k−1+ϕi,j,k+1+6ϕi−1,j−1,k+ϕi−1,j+1,k+ϕi+1,j−1,k+ϕi+1,j+1,k+ϕi,j−1,k−1+ϕi,j−1,k+1+ϕi,j+1,k−1+ϕi,j+1,k+1+ϕi−1,j,k−1+ϕi−1,j,k+1+ϕi+1,j,k−1+ϕi+1,j,k+1+ϕi−1,j−1,k−1+ϕi−1,j−1,k+1+ϕi−1,j+1,k−1+ϕi−1,j+1,k+1+ϕi+1,j−1,k−1+ϕi+1,j−1,k+1+ϕi+1,j+1,k−1+ϕi+1,j+1,k+1−200ϕijk/48h2−1.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 the2D heat equation with the periodic boundary condition (3) on a domain Ω=0,1×0,1:(9)utx,y,t=Δux,y,t,forx,y∈Ω,t>0.For example, we use the following initial condition for equation (9):(10)ux,y,0=cos2πxcos2πy.Then, the exact solution isux,y,t=e−8π2tcos2πxcos2πy. Numerical convergence tests of two different schemes, standard and isotropic schemes, are performed with Nx=Ny=2n and h=1/Nx, Δt=10−9. The final time is T=NtΔt for n=3,4,5, and 6. Here, Nt=20 is used. Let eijh:=uijNt−uxi,yj,T be the numerical error in Ωh. We define the discrete l2-norm of the error as e2=∑i=1Nx∑j=1Nyeij2/NxNy. Then, the rate of convergence is log2eh2/eh/22. Table 1 shows the discrete l2-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.Table 1
l2-norm of the errors and convergence rates of standard and isotropic schemes.
StandardIsotropichl2 errorOrderl2 errorOrder1/83.976×10−87.636×10−81/161.009×10−81.9781.998×10−81.9341/322.533×10−91.9945.054×10−91.9831/646.340×10−101.9991.267×10−91.996


### 3.2. Comparison of the Two Schemes for an Isotropic Initial Shape

We consider the tumor growth simulation parameters in [12]:(11)γϕ=λg1+ϕ21−ϕ2−λd21+ϕ,where λg and λd are the growth and death coefficients, respectively. We use λg=280 and λd=92.Figure1 shows a plot of γϕ=2801+ϕ21−ϕ2−461+ϕ. Note that γϕ=0 when ϕ=−1 and ϕ≈−0.957109,0.688032,1.26908.Figure 1
Plot ofγϕ=2801+ϕ21−ϕ2−461+ϕ.


#### 3.2.1.2D Case

We present the evolution of a disk by using the two types of schemes: standard and isotropic. The initial condition in the domainΩ=−0.6,0.6×−0.6,0.6 is given by(12)ϕx,y,0=tanh0.1−x2+y22ϵ.In this test, the uniform space step sizeh=1.2/128, the time step size Δt=10h2, 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.Figure 2
(a, b) Zero-level contours for the circular evolution in equation (12) with standard and isotropic schemes up to t=1200Δt at every 120Δt, respectively.
(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 as(13)rd=l24πA,where l implies the perimeter of zero-level contour and A implies the area enclosed by the contour. Note that rd≈1 if an evolution maintains the form of the disk. Figure 3 shows the ratios rd of the standard scheme and the isotropic scheme.Figure 3
Plot ofrd among standard and isotropic schemes.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.3D Case

We further present the evolution of sphere inΩ=−0.6,0.63 by using both standard and isotropic schemes with the initial condition given by(14)ϕx,y,z,0=tanh0.1−x2+y2+z22ϵ.Figure4 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=10h2, ϵ=0.0125, and the final time t=1600Δt.Figure 4
(a) Overlapped isosurfacesϕ=0 for evolution in equation (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.
(a)(b)Obviously, as shown in Figure4(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 rs such as each surface area over volume with appropriate scaling to compare evolutionary detail of both schemes, and hence rs is defined as(15)rs=S336πV2,where S is a surface area of sphere and V is a volume of sphere. Note that rs≈1 if an evolution maintains the form of the sphere. Figure 5 shows the ratios rs of both schemes. Thus, we conclude that the isotropic scheme is relatively unbiased and has a correct evolution.Figure 5
Plot ofrs from standard and isotropic schemes.


### 3.3. Comparison of the Two Schemes for an Anisotropic Initial Shape

In Section3.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.2D Case

We investigate tumor growth of an elliptical shape in2D 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:(16)ϕx,y,0=tanh0.1−2x2+0.25y22ϵ,(17)ϕx,y,0=tanh0.1−2X2+0.25Y22ϵ,where X and Y are defined as(18)XY=cosθsinθ−sinθcosθxy.Here,θ=−π/4, h=4.8/512, Δt=h2, 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. The final row shows overlapped contours for each column at t=0, 15000Δt, 35000Δt, and 50000Δt. The results under the same conditions by using the isotropic scheme are given in Figure 7.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.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)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.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)


#### 3.3.2.3D Case

We further investigate the evolution of ellipsoidal shape inΩ=−1.0,1.03, where initial conditions are given, respectively, as follows:(19)ϕx,y,z,0=tanh0.2−0.5x2+2y2+2z22ϵ,(20)ϕx,y,z,0=tanh0.2−0.5X2+2Y2+2Z22ϵ,where X, Y, and Z are defined as(21)XYZ=RxRyRzxyz.Note thatRx, Ry, and Rz are rotation matrices as follows:(22)Rx=1000cosθsinθ0−sinθcosθ,Ry=cosθ0−sinθ010sinθ0cosθ,Rz=cosθsinθ0−sinθcosθ0001,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=10h2, ϵ=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. The 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. 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.Figure 8
(a–c) Snapshots of isosurfaceϕ=0 for evolution in equation (9) using the standard scheme at t=0, 600Δt, and 1000Δt, respectively. (d–f) Snapshots of isosurface ϕ=0 for evolution in equation (10) using the standard scheme at t=0, 600Δt, and 1000Δt, respectively. (g–i) Overlapped isosurfaces ϕ=0 for each column at t=0, 600Δt, and 1000Δt, respectively. Note that the axes are limited to show only −0.6,0.63.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 9
(a–c) Snapshots of isosurfaceϕ=0 for evolution in equation (9) using the isotropic scheme at t=0, 600Δt, and 1000Δt, respectively. (d–f) Snapshots of isosurface ϕ=0 for evolution in equation (10) using the isotropic scheme at t=0, 600Δt, and 1000Δt, respectively. (g–i) Overlapped isosurfaces ϕ=0 for each column at t=0, 600Δt, and 1000Δt, respectively. Note that the axes are limited to show only −0.6,0.63.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 3.1. Convergence Test

We start with numerical convergence tests of two discrete Laplace schemes. We consider the2D heat equation with the periodic boundary condition (3) on a domain Ω=0,1×0,1:(9)utx,y,t=Δux,y,t,forx,y∈Ω,t>0.For example, we use the following initial condition for equation (9):(10)ux,y,0=cos2πxcos2πy.Then, the exact solution isux,y,t=e−8π2tcos2πxcos2πy. Numerical convergence tests of two different schemes, standard and isotropic schemes, are performed with Nx=Ny=2n and h=1/Nx, Δt=10−9. The final time is T=NtΔt for n=3,4,5, and 6. Here, Nt=20 is used. Let eijh:=uijNt−uxi,yj,T be the numerical error in Ωh. We define the discrete l2-norm of the error as e2=∑i=1Nx∑j=1Nyeij2/NxNy. Then, the rate of convergence is log2eh2/eh/22. Table 1 shows the discrete l2-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.Table 1
l2-norm of the errors and convergence rates of standard and isotropic schemes.
StandardIsotropichl2 errorOrderl2 errorOrder1/83.976×10−87.636×10−81/161.009×10−81.9781.998×10−81.9341/322.533×10−91.9945.054×10−91.9831/646.340×10−101.9991.267×10−91.996


## 3.2. Comparison of the Two Schemes for an Isotropic Initial Shape

We consider the tumor growth simulation parameters in [12]:(11)γϕ=λg1+ϕ21−ϕ2−λd21+ϕ,where λg and λd are the growth and death coefficients, respectively. We use λg=280 and λd=92.Figure1 shows a plot of γϕ=2801+ϕ21−ϕ2−461+ϕ. Note that γϕ=0 when ϕ=−1 and ϕ≈−0.957109,0.688032,1.26908.Figure 1
Plot ofγϕ=2801+ϕ21−ϕ2−461+ϕ.


### 3.2.1.2D Case

We present the evolution of a disk by using the two types of schemes: standard and isotropic. The initial condition in the domainΩ=−0.6,0.6×−0.6,0.6 is given by(12)ϕx,y,0=tanh0.1−x2+y22ϵ.In this test, the uniform space step sizeh=1.2/128, the time step size Δt=10h2, 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.Figure 2
(a, b) Zero-level contours for the circular evolution in equation (12) with standard and isotropic schemes up to t=1200Δt at every 120Δt, respectively.
(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 as(13)rd=l24πA,where l implies the perimeter of zero-level contour and A implies the area enclosed by the contour. Note that rd≈1 if an evolution maintains the form of the disk. Figure 3 shows the ratios rd of the standard scheme and the isotropic scheme.Figure 3
Plot ofrd among standard and isotropic schemes.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.3D Case

We further present the evolution of sphere inΩ=−0.6,0.63 by using both standard and isotropic schemes with the initial condition given by(14)ϕx,y,z,0=tanh0.1−x2+y2+z22ϵ.Figure4 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=10h2, ϵ=0.0125, and the final time t=1600Δt.Figure 4
(a) Overlapped isosurfacesϕ=0 for evolution in equation (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.
(a)(b)Obviously, as shown in Figure4(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 rs such as each surface area over volume with appropriate scaling to compare evolutionary detail of both schemes, and hence rs is defined as(15)rs=S336πV2,where S is a surface area of sphere and V is a volume of sphere. Note that rs≈1 if an evolution maintains the form of the sphere. Figure 5 shows the ratios rs of both schemes. Thus, we conclude that the isotropic scheme is relatively unbiased and has a correct evolution.Figure 5
Plot ofrs from standard and isotropic schemes.


## 3.2.1.2D Case

We present the evolution of a disk by using the two types of schemes: standard and isotropic. The initial condition in the domainΩ=−0.6,0.6×−0.6,0.6 is given by(12)ϕx,y,0=tanh0.1−x2+y22ϵ.In this test, the uniform space step sizeh=1.2/128, the time step size Δt=10h2, 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.Figure 2
(a, b) Zero-level contours for the circular evolution in equation (12) with standard and isotropic schemes up to t=1200Δt at every 120Δt, respectively.
(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 as(13)rd=l24πA,where l implies the perimeter of zero-level contour and A implies the area enclosed by the contour. Note that rd≈1 if an evolution maintains the form of the disk. Figure 3 shows the ratios rd of the standard scheme and the isotropic scheme.Figure 3
Plot ofrd among standard and isotropic schemes.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.3D Case

We further present the evolution of sphere inΩ=−0.6,0.63 by using both standard and isotropic schemes with the initial condition given by(14)ϕx,y,z,0=tanh0.1−x2+y2+z22ϵ.Figure4 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=10h2, ϵ=0.0125, and the final time t=1600Δt.Figure 4
(a) Overlapped isosurfacesϕ=0 for evolution in equation (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.
(a)(b)Obviously, as shown in Figure4(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 rs such as each surface area over volume with appropriate scaling to compare evolutionary detail of both schemes, and hence rs is defined as(15)rs=S336πV2,where S is a surface area of sphere and V is a volume of sphere. Note that rs≈1 if an evolution maintains the form of the sphere. Figure 5 shows the ratios rs of both schemes. Thus, we conclude that the isotropic scheme is relatively unbiased and has a correct evolution.Figure 5
Plot ofrs from standard and isotropic schemes.


## 3.3. Comparison of the Two Schemes for an Anisotropic Initial Shape

In Section3.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.2D Case

We investigate tumor growth of an elliptical shape in2D 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:(16)ϕx,y,0=tanh0.1−2x2+0.25y22ϵ,(17)ϕx,y,0=tanh0.1−2X2+0.25Y22ϵ,where X and Y are defined as(18)XY=cosθsinθ−sinθcosθxy.Here,θ=−π/4, h=4.8/512, Δt=h2, 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. The final row shows overlapped contours for each column at t=0, 15000Δt, 35000Δt, and 50000Δt. The results under the same conditions by using the isotropic scheme are given in Figure 7.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.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)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.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)


### 3.3.2.3D Case

We further investigate the evolution of ellipsoidal shape inΩ=−1.0,1.03, where initial conditions are given, respectively, as follows:(19)ϕx,y,z,0=tanh0.2−0.5x2+2y2+2z22ϵ,(20)ϕx,y,z,0=tanh0.2−0.5X2+2Y2+2Z22ϵ,where X, Y, and Z are defined as(21)XYZ=RxRyRzxyz.Note thatRx, Ry, and Rz are rotation matrices as follows:(22)Rx=1000cosθsinθ0−sinθcosθ,Ry=cosθ0−sinθ010sinθ0cosθ,Rz=cosθsinθ0−sinθcosθ0001,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=10h2, ϵ=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. The 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. 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.Figure 8
(a–c) Snapshots of isosurfaceϕ=0 for evolution in equation (9) using the standard scheme at t=0, 600Δt, and 1000Δt, respectively. (d–f) Snapshots of isosurface ϕ=0 for evolution in equation (10) using the standard scheme at t=0, 600Δt, and 1000Δt, respectively. (g–i) Overlapped isosurfaces ϕ=0 for each column at t=0, 600Δt, and 1000Δt, respectively. Note that the axes are limited to show only −0.6,0.63.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 9
(a–c) Snapshots of isosurfaceϕ=0 for evolution in equation (9) using the isotropic scheme at t=0, 600Δt, and 1000Δt, respectively. (d–f) Snapshots of isosurface ϕ=0 for evolution in equation (10) using the isotropic scheme at t=0, 600Δt, and 1000Δt, respectively. (g–i) Overlapped isosurfaces ϕ=0 for each column at t=0, 600Δt, and 1000Δt, respectively. Note that the axes are limited to show only −0.6,0.63.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 3.3.1.2D Case

We investigate tumor growth of an elliptical shape in2D 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:(16)ϕx,y,0=tanh0.1−2x2+0.25y22ϵ,(17)ϕx,y,0=tanh0.1−2X2+0.25Y22ϵ,where X and Y are defined as(18)XY=cosθsinθ−sinθcosθxy.Here,θ=−π/4, h=4.8/512, Δt=h2, 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. The final row shows overlapped contours for each column at t=0, 15000Δt, 35000Δt, and 50000Δt. The results under the same conditions by using the isotropic scheme are given in Figure 7.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.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)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.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)


## 3.3.2.3D Case

We further investigate the evolution of ellipsoidal shape inΩ=−1.0,1.03, where initial conditions are given, respectively, as follows:(19)ϕx,y,z,0=tanh0.2−0.5x2+2y2+2z22ϵ,(20)ϕx,y,z,0=tanh0.2−0.5X2+2Y2+2Z22ϵ,where X, Y, and Z are defined as(21)XYZ=RxRyRzxyz.Note thatRx, Ry, and Rz are rotation matrices as follows:(22)Rx=1000cosθsinθ0−sinθcosθ,Ry=cosθ0−sinθ010sinθ0cosθ,Rz=cosθsinθ0−sinθcosθ0001,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=10h2, ϵ=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. The 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. 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.Figure 8
(a–c) Snapshots of isosurfaceϕ=0 for evolution in equation (9) using the standard scheme at t=0, 600Δt, and 1000Δt, respectively. (d–f) Snapshots of isosurface ϕ=0 for evolution in equation (10) using the standard scheme at t=0, 600Δt, and 1000Δt, respectively. (g–i) Overlapped isosurfaces ϕ=0 for each column at t=0, 600Δt, and 1000Δt, respectively. Note that the axes are limited to show only −0.6,0.63.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 9
(a–c) Snapshots of isosurfaceϕ=0 for evolution in equation (9) using the isotropic scheme at t=0, 600Δt, and 1000Δt, respectively. (d–f) Snapshots of isosurface ϕ=0 for evolution in equation (10) using the isotropic scheme at t=0, 600Δt, and 1000Δt, respectively. (g–i) Overlapped isosurfaces ϕ=0 for each column at t=0, 600Δt, and 1000Δt, respectively. Note that the axes are limited to show only −0.6,0.63.
(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.



---

*Source: 1290895-2021-12-06.xml*