# Fast and Accurate Numerical Solution of Allen–Cahn Equation

Authors: Yongho Kim; Gilnam Ryu; Yongho Choi
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Simulation speed depends on code structures. Hence, it is crucial how to build a fast algorithm. We solve the Allen–Cahn equation by an explicit finite difference method, so it requires grid calculations implemented by many for-loops in the simulation code. In terms of programming, many for-loops make the simulation speed slow. We propose a model architecture containing a pad and a convolution operation on the Allen–Cahn equation for fast computation while maintaining accuracy. Also, the GPU operation is used to boost up the speed more. In this way, the simulation of other differential equations can be improved. In this paper, various numerical simulations are conducted to confirm that the Allen–Cahn equation follows motion by mean curvature and phase separation in two-dimensional and three-dimensional spaces. Finally, we demonstrate that our algorithm is much faster than an unoptimized code and the CPU operation.


---

## Body

## 1. Introduction

The Allen–Cahn (AC) equation is a reaction-diffusion equation composed of the reaction term−F′ϕx,t/ϵ2 and the diffusion term Δϕx,t:(1)ϕtx,t=−F′ϕx,tϵ2+Δϕx,t,x∈Ω,t>0,where ϕx,t is the order parameter which is defined as the difference in concentration of the two components in a mixture; Fϕ is double well potential energy function with minimum values at −1 and 1, and its form is Fϕ=0.25ϕ2−12. ϵ is the thickness of the transition layer which is a small positive constant value. The AC equation is first introduced in a research on the phase separation of binary iron alloys [1]. The AC equation was studied and applied to various fields such as image inpainting [2–4], image segmentation [5, 6], and crystal growth [7–9]. Also, there are various methods for numerically solving the reaction-diffusion type equations. In [10–12], the authors use the fourth-order exponential time differencing Runge–Kutta method to solve reaction-diffusion type equation as Burgers–Huxley and Gray–Scott model, respectively. There is a study of solving equations using an adaptive method [13–15]. The authors of [16, 17] considered the Fourier spectral method to solve system equations.With the development of computer hardware such as GPU (graphics processing unit) and memory cards, neural networks are applied in a wide range of research areas such as computer vision, natural language processing, and numerical analysis. As GPU operations outperform CPU (central processing unit) performance in multitasks and high-dimensional problems, open source machine learning libraries such as Pytorch provide a variety of neural networks using GPU and are useful to build the architecture combined with neural networks and numerical methods. Thus, many research studies use machine learning libraries. Raissi et al. [18] proposed physics-informed neural networks combined by multilayer perceptrons and numerical methods to solve nonlinear partial differential equations. Karumuri et al. [19] introduced a solver-free approach for stochastic partial differential equations, and Yang et al. [20] proposed a Bayesian physics-informed neural network. In this paper, we propose a structure using padding and convolution operation for the GPU calculation of the AC equation and demonstrate the validity of the proposed structure by verifying the result with the Python code that has the same mathematical meaning.This paper is organized as follows. In Section2, we present an explicit finite difference method to solve the AC equation which is implemented by CPU and GPU algorithms. In Section 3, the numerical simulations including a motion by mean curvature, phase separation, and temporal evolutions of various initial shapes are introduced as well as the runtime results between CPU and GPU operations are compared. Finally, conclusions are drawn in Section 4.


## 2. Numerical Solutions

In this section, we present an explicit finite difference method to solve AC equation (1). Also, we give an explanation of each algorithm for CPU and GPU computing. For simplicity of expression, we describe a numerical scheme for the AC equation in two dimensions (2D), and the definition in the three-dimensional (3D) space can be easily extended and considered. A computational domain is defined using a uniform grid of size h=1/Nx and Ωh=xi,yj=a+i−0.5h,c+j−0.5h for 1≤i≤Nx, 1≤j≤Ny is the set of cell centers. Here, NxandNy are mesh sizes on computational domain a,b×c,d. For the definition of the boundary condition, we define the extended computational domain as follows:(2)Ωh=xi,yj=a+i−1.5h,c+j−1.5h,for 1≤i≤Nx+2 and 1≤j≤Ny+2. Let ϕijn be approximations of ϕxi,yj,nΔt, where Δt=0.1h2 is the temporal step size, T is a final time, and Nt is the total number of time steps. The boundary condition is zero Neumann boundary condition:(3)ϕi,1n=ϕi,2n,ϕi,Ny+2n=ϕi,Ny+1n,1≤i≤Nx,ϕ1,jn=ϕ2,jn,ϕNx+2,jn=ϕNx+1,jn,1≤j≤Ny.We define the thickness of the transition layerϵ in equation (1) as ϵm [21]:(4)ϵm=hm22tanh−10.9,where m is the number of grids representing the thickness.


### 2.1. Numerical Solutions on CPU (Baseline)

The AC equation (1) is discretized using the explicit finite difference method as(5)ϕijn+1−ϕijnΔt=ϕijn−ϕijn3ϵ2+Δhϕijn,ϕijn+1=ϕijn+Δtϕijn−ϕijn3ϵ2+Δhϕijn,ϕijn+1=1+αϕijn−αϕijn3+ΔtΔhϕijn,where α=Δt/ϵ2 and Δhϕijn=ϕi−1,jn+ϕi+1,jn+ϕi,j−1n+ϕi,j+1n−4ϕijn/h2. The baseline algorithm is implemented in equation (5) using Numpy (CPU array).


### 2.2. Numerical Solutions on GPU (Pytorch)

For GPU computing, the AC equation (1) can be expressed using Pytorch, and the algorithm can be represented as Figure 1.where f is a Pytorch model. The main Algorithm 1 is a set of steps in equation (6) using Pytorch.Figure 1
Schematic of the process to obtainf.(6)ϕn+1−ϕnΔt=ϕn−ϕn3ϵ2+Δhϕn,ϕn+1=ϕn+Δtϕn−ϕn3ϵ2+Δhϕn,ϕn+1=1+αϕn−αϕn3+convϕnα=Δtϵ2,convϕn=ΔtΔhϕn,ϕn+1=fϕn,Algorithm 1: Our CNN model for solving AC equation (6).
class Net (nn.Module):def __init__(self, h2, dt, eps, device):super (Net, self).__init__()# 2nd order differencing filterself.lap = torch.Tensor ([[[[0., 1., 0.], [1., −4., 1], [0., 1., 0.]]]]).to (device)self.pad = nn.ReplicationPad2d (1) #Replication pad for boundary condition.self.alpha = dt/eps∗∗ 2self.beta = dt/h2def forward (self,x):u_pad = self.pad (x) #boundary conditionreaction = F.conv2d (u_pad, self.lap) #reaction termx = (1 + self.alpha) ∗x-self.alpha ∗x∗∗ 3 + self.beta ∗ reactionreturnxIn this Algorithm1, nn.ReplicationPad2d (or nn.ReplicationPad3d) is applied to satisfy the boundary condition by padding the ϕn input using replication of the input boundary. Also, the convolutional operator F.conv2d (or F.conv3d) with the 2nd order differencing filter is used to calculate the diffusion term Δϕn. The model can choose an operation mode between GPU and CPU by to (device) in Algorithm 1. If device = cuda:0, the model is implemented on GPU, otherwise on CPU. All the codes are available from the first author’s GitHub web page (https://github.com/kimy-de/gpuallencahn) and the corresponding author’s web page (https://sites.google.com/view/yh-choi/code).


## 2.1. Numerical Solutions on CPU (Baseline)

The AC equation (1) is discretized using the explicit finite difference method as(5)ϕijn+1−ϕijnΔt=ϕijn−ϕijn3ϵ2+Δhϕijn,ϕijn+1=ϕijn+Δtϕijn−ϕijn3ϵ2+Δhϕijn,ϕijn+1=1+αϕijn−αϕijn3+ΔtΔhϕijn,where α=Δt/ϵ2 and Δhϕijn=ϕi−1,jn+ϕi+1,jn+ϕi,j−1n+ϕi,j+1n−4ϕijn/h2. The baseline algorithm is implemented in equation (5) using Numpy (CPU array).


## 2.2. Numerical Solutions on GPU (Pytorch)

For GPU computing, the AC equation (1) can be expressed using Pytorch, and the algorithm can be represented as Figure 1.where f is a Pytorch model. The main Algorithm 1 is a set of steps in equation (6) using Pytorch.Figure 1
Schematic of the process to obtainf.(6)ϕn+1−ϕnΔt=ϕn−ϕn3ϵ2+Δhϕn,ϕn+1=ϕn+Δtϕn−ϕn3ϵ2+Δhϕn,ϕn+1=1+αϕn−αϕn3+convϕnα=Δtϵ2,convϕn=ΔtΔhϕn,ϕn+1=fϕn,Algorithm 1: Our CNN model for solving AC equation (6).
class Net (nn.Module):def __init__(self, h2, dt, eps, device):super (Net, self).__init__()# 2nd order differencing filterself.lap = torch.Tensor ([[[[0., 1., 0.], [1., −4., 1], [0., 1., 0.]]]]).to (device)self.pad = nn.ReplicationPad2d (1) #Replication pad for boundary condition.self.alpha = dt/eps∗∗ 2self.beta = dt/h2def forward (self,x):u_pad = self.pad (x) #boundary conditionreaction = F.conv2d (u_pad, self.lap) #reaction termx = (1 + self.alpha) ∗x-self.alpha ∗x∗∗ 3 + self.beta ∗ reactionreturnxIn this Algorithm1, nn.ReplicationPad2d (or nn.ReplicationPad3d) is applied to satisfy the boundary condition by padding the ϕn input using replication of the input boundary. Also, the convolutional operator F.conv2d (or F.conv3d) with the 2nd order differencing filter is used to calculate the diffusion term Δϕn. The model can choose an operation mode between GPU and CPU by to (device) in Algorithm 1. If device = cuda:0, the model is implemented on GPU, otherwise on CPU. All the codes are available from the first author’s GitHub web page (https://github.com/kimy-de/gpuallencahn) and the corresponding author’s web page (https://sites.google.com/view/yh-choi/code).


## 3. Numerical Experiments

In this section, we perform the following numerical tests in 2D and 3D:ϵm effect, the motion by mean curvature effect with various initial shapes (circle (sphere in 3D), dumbbell, star, torus, maze), and phase separation. In the ϵm effect test, we compare the numerical solutions to analytic values to find a proper ϵm value. We simulate a phenomenon that follows motion by mean curvature in various initial shapes and phase separation with random initial condition. We perform the simulation on the following specifications: Intel (R) Core (TM) i9-9900K CPU @3.60 GHz, 32 GB RAM/NVIDIA GeForce RTX 2080 Super.


### 3.1. Convergence Tests

We perform convergence tests for time and space. We define the discreteL2 norm as(7)eL2d=h2∑i=1Nx∑j=1Nyeijn2,where eijn=ϕijn−ϕrefn and ϕrefn is reference solution. We use the following initial condition:(8)ϕx,y,0=0.1cos2πxcos2πy,with zero Neumann boundary condition on the computational domain 0,1×0,1.First, to show the convergence rate with respect to time, the space step sizeh=0.0025 (i.e., Nx=Ny=400) is fixed, and we use the reference solution with Nx=Ny=400 and Δt=0.001h2. Then, the test proceeds to time t=0.0001 for time steps Δt=0.2h2,0.1h2, and 0.05h2. The temporal convergence rate is defined as log2eΔtL2d/eΔt/2L2d. Table 1 shows that our proposed algorithm has first-order accuracy for time.Table 1
TemporalL2 errors and convergence rates at t=0.0001.
Δt0.2h20.1h20.05h2L2 error9.282e − 34.759e − 32.487e − 3Rate0.9630.936Next, we show the spatial accuracy of our proposed algorithm. To show the spatial convergence rate, the time step was fixed toΔt=0.1h2 for each mesh size h=0.01,0.005, and 0.0025. For the reference solution, the time step was set to 0.001h2 for each mesh size. The spatial convergence rate is defined as log2ehL2d/eh/2L2d. Table 2 shows that our proposed algorithm has between second- and third-order for the spatial accuracy.Table 2
SpatialL2 errors and convergence rates at t=0.0001.
h0.010.0050.0025L2 error1.011e − 46.722e − 044.760e − 3Rate2.7332.824


### 3.2. Energy Dissipation and Maximum Principle

AC equation (1) follows the energy dissipation law, and the equation is derived from(9)Eϕ=∫ΩFϕϵ2+12∇ϕ2dx.Eϕ is decreasing with time(10)∂∂tEϕ=−∫Ω∂ϕ∂t2dx≤0.To check a discrete energy, we can rewrite equation (9) as Edϕ. Figure 2 shows mesh plots and a discrete energy graph when time evolution is performed with random initial condition. The tested initial condition is used as equation (17a) and other parameters are used Nx = Ny = 200, space step size h = 1/Nx, time step size Δt = 0.1h2, and thickness of transition layer ϵ10 on the computational domain Ω = (0, 1)2.Figure 2
(a–d) Mesh plots at the marked points from left to right of the energy graph (e). (e) Time-dependent normalized discrete total energyEdϕt/Edϕ0. (f) Maximum and minimum values of ϕ over time evolution.
(a)(b)(c)(d)(e)(f)Figure 3
(a–d) The changes in the circle over time are shown, and each time is described in the figures. (e) Change of radius for variousϵm.
(a)(b)(c)(d)(e)Also, we can obtain that the AC equation does not conserve the initial mass.(11)ddt∫Ωϕdx=∫Ωϕtdx=∫Ω−F′ϕϵ2+Δϕdx=−∫ΩF′ϕϵ2dx≤0,with zero Neumann boundary condition.


### 3.3. Initial Conditions

We consider the 2D and 3D initial conditions introduced in this section: to check theϵm effect, we measure the circle’s (sphere in 3D) radius that changes with time and take the initial conditions:(12a)ϕx,y,0=tanhR0−x−0.52+y−0.522ϵ,in2D,(12b)ϕx,y,z,0=tanhR0−x−0.52+y−0.52+z−0.522ϵ,in3D,where R0 is the initial radius of a circle (sphere in 3D).In various tests to check the motion by the mean curvature, the initial conditions for the circle and sphere refer to equations (12a) and (12b), respectively. A dumbbell is constructed by following the initial conditions (o/w: otherwise)(13a)ϕx,y,0=1.0,if0.4<x<1.6and0.4<y<0.6,1+tanhR0−x−0.32+Y2ϵ+tanhR0−x−1.72+Y2ϵ,otherwise,in2D,(13b)ϕx,y,z,0=1.0,if0.4<x<1.6and0.4<y,z<0.6,1+tanhR0−x−0.32+YZ2ϵ+tanhR0−x−1.72+YZ2ϵ,otherwise,in3D.where R0 is the initial radius of both sides of dumbbell’s circle (sphere in 3D) and for simplicity of expression, Y=y−0.52 and YZ=y−0.52+z−0.52.The initial conditions of a star shape are defined as(14a)ϕx,y,0=tanh0.25+0.1cos6θ−x−0.52+y−0.522ϵ,in2D,(14b)ϕx,y,z,0=tanh0.7+0.2cos6θ−x2+y2+z22ϵ,in3D,where(15)θ=tan−1y−0.5x−0.5,ifx>0.5,π+tan−1y−0.5x−0.5,otherwise,in2D,θ=tan−1zx,ifx>0.5,π+tan−1zx,otherwise.in3D.and we use different domain sizes in 2D and 3D, so the center of the star depends on the dimensions.And, a torus shape is given by(16a)ϕx,y,0=−1+tanhR1−XY2ϵ−tanhR2−XY2ϵ,in2D,(16b)ϕx,y,z,0=z2+x2+y2−R12−R2,in3D,where R1 and R2 are the radius of major (outside) and minor (inside) circles, respectively. And, for simplicity of expression, XY=x−0.52+y−0.52.The initial conditions of a maze shape are complicated to describe its equation, so refer to the code (https://github.com/kimy-de/gpuallencahn, https://sites.google.com/view/yh-choi/code).And, the last initial condition is random for confirming the phase separation of the AC equation:(17a)ϕx,y,0=0.1randx,y,in2D,(17b)ϕx,y,z,0=0.1randx,y,z,in3D.Here, the function randx,y has a random value between −1 and 1.


### 3.4. Simulations in the 2-Dimensional Space

Unless otherwise stated, we use the following parameters: mesh sizeNx=Ny=200, space step size h=1/Nx, time step size Δt=0.1h2, and computational domain Ω=0,1×0,1. We should find a proper thickness of the transition layer as defined in equation (4). First, we find an appropriate ϵm by comparing the numerical solution with the exact solution for the radius that decreases due to the motion by mean curvature in the circle initial condition in equation (12a).Figures3(a)–3(d) show a circle shrinking with motion by mean curvature based on ϵ10. The exact solution of the radius r decreasing with time evolution can be calculated [22]:(18)rt=r02+21−dt,where r0 is the initial radius of the circle, d is the dimension, and t is time. We simulate various ϵm until the final time T=0.03. As shown in Figure 3, when ϵ10 is used, it is found that the exact solution and the numerical solution are most similar in the experiment. Therefore, the other tests are performed using ϵ10 except for the case of changing the grid in 2D.Figure4 shows the temporal evolution of the dumbbell shape, and the initial condition is shown in equation (13a). For resolution, we use Nx=400 and Ny=200 on the computational domain Ω=0,2×0,1. The final time T of the simulation is 0.0094 and R0=0.2. In Figure 4(e), the changing direction by the motion by mean curvature is indicated by arrows.Figure 4
(a–d) Time evolution of dumbbell shape, and each time is described in the figures. (e) Contour lines over time are shown by overlapping. Changing by the mean curvature flow can be seen.
(a)(b)(c)(d)(e)Figure5 shows the evolution of the star shape created by equation (14a). The parameters are used as mentioned at the beginning of this section and T=0.0325. As shown in Figure 5, the tips of the star move inward and the gaps between the tips move outward. When it changes to the shape of a circle, the change of the radius can be predicted as shown in Figure 3.Figure 5
Time evolution of a star shape. As in the shape of a dumbbell, it can be seen that it changes with the mean curvature flow. (a)t = 0, (b) t = 0.0025, (c) t = 0.0188, (d) t = 0.0325.
(a)(b)(c)(d)Figure6 shows the evolution of the torus shape of equation (16a) with T=0.0575, R1=0.4, and R2=0.3. Because the inner circle has a larger curvature, it shrinks faster than the outer circle, and after the inner circle disappears, the change of radius over time can be measured as shown in Figure 3.Figure 6
Time evolution of a torus shape, and each time is described in the figures. The inner circle has larger curvature than the outer circle, so the inner circle shrinks faster than the outer circle. (a)t = 0, (b) t = 0.0188, (c) t = 0.0375, (d) t = 0.0575.
(a)(b)(c)(d)Figure7 shows the evolution of a maze shape. We use Nx=Ny=100, ϵ5, and T=0.04. As shown in Figure 7, we obtain the results of shrinking while maintaining its initial shape.Figure 7
Time evolution of a maze shape on100×100 mesh size with ϵ5, and T=0.04. (a) t = 0, (b) t = 0.01, (c) t = 0.025, (d) t = 0.04.
(a)(b)(c)(d)The last simulation in 2D is phase separation with a random initial condition equation (17a). In Figure 8, starting with random values with 0.1 amplitude, but over time, phase separation occurs with values −1 to 1.Figure 8
Time evolution of phase separation with a random initial condition. (a)t = 0, (b) t = 0.0025, (c) t = 0.0038, (d) t = 0.0125.
(a)(b)(c)(d)According to the results in Table3, the speed gap between CPU and GPU is significant. In 2D, it is up to 251.6 times the difference between the Python:CPU and the Pytorch:GPU codes. Also, Pytorch:GPU tensors make the model up to 4.73 times faster than Pytorch:CPU tensors in the same code.Table 3
Runtime Result in 2D(sec).
IterationsInitial valueCircleDumbbellStarTorusMazeRandom12001150011300123001400112001CPU:Python680.09 (185.31)1736.05 (251.6)729.74 (128.47)1326.82 (165.65)60.15 (52.3)674.63 (119.62)CPU:Pytorch14.73 (4.01)30.67 (4.44)17.60 (3.1)37.85 (4.73)2.98 (2.59)17.79 (3.15)GPU:Pytorch3.676.905.688.011.155.64The values in parentheses describe how many times the difference is based on GPU:Pytorch time for each test.


### 3.5. Simulations in the 3-Dimensional Space

Three-dimensional simulations can be considered as an extension of two-dimensional tests. Unless otherwise stated, we use the following parameters: mesh sizeNx=Ny=Nz=100, space step size h=1/Nx, time step size Δt=0.1h2, and computational domain Ω=0,1×0,1×0,1. As in the 2D simulation, we find an appropriate ϵm by comparing the numerical solution with the exact solution (using in equation (18)) for the radius of the sphere initial condition in equation (12b).Figures9(a)–9(d) show the time evolution results for ϵ12, and (e) presents the comparison of the numerical and exact solutions of radius for various epsilons. As shown in Figure 9, when ϵ12 is used, it is found that the exact solution and the numerical solution are most similar. Therefore, the other tests are performed using ϵ12.Figure 9
(a–d) The changes in the circle over time are shown, and each time is described in the figures. (e) Change of radius for variousϵm.
(a)(b)(c)(d)(e)Figure10 shows the temporal evolution of the dumbbell shape, and the initial condition is shown in equation (13b). For resolution, we use Nx=200 and Ny=Nz=100 on the computational domain Ω=0,2×0,1×0,1. The final time T of the simulation is T=0.0025 and R0=0.25. The tendency of the 3D dumbbell motion is different from the result of 2D dumbbell. The reason is that, in the initial shape in Figure 10(a), the radius of the handle is much smaller than the radius of the spheres at both ends. Therefore, the handle part shrinks faster than the end of spheres and breaks as in Figure 10.Figure 10
(a–d) Time evolution of the dumbbell shape. The handle shrinks quickly and breaks occur because the curvature of the handle is larger than that of both spheres. (a)t = 0, (b) t = 0.0006, (c) t = 0.0013, (d) t = 0.0025.
(a)(b)(c)(d)Figure11 shows the evolution of the star shape shown in equation (14b) on the computational domain Ω=−1,1×−1,1×−1,1. The parameters are used as mentioned at the beginning of this section and T=0.02. As shown in Figure 11, similar to the 2D result, the tips of the star move inward, and the gaps between the tips move outward.Figure 11
Time evolution of a star shape. The shape changes with the mean curvature flow. (a)t = 0, (b) t = 0.0025, (c) t = 0.0075, (d) t = 0.02.
(a)(b)(c)(d)Figure12 shows the evolution of the torus shape with initial condition in equation (16b) using R1=0.3 and R2=0.3 and T=0.01 on the computational domain Ω=−1,1×−1,1×−1,1. Contrary to the 2D result, the inner circle becomes larger because the radius of the minor circle exists. Therefore, the radius of the major circle and the minor circle each has different curvature. Since the radius of the minor circle is smaller than major circle, the mean curvature drives the motion into the inside of the torus as shown in Figure 12.Figure 12
Time evolution of a torus shape. Contrary to the 2D results, the inner circle increases due to the effect of the mean curvature. (a)t = 0, (b) t = 0.005, (c) t = 0.008, (d) t = 0.01.
(a)(b)(c)(d)Figure13 shows the evolution of a maze shape. The initial condition is described in code (https://github.com/kimy-de/gpuallencahn, https://sites.google.com/view/yh-choi/code). In this simulation, we use T=0.0175 and computational domain Ω=−1,1×−1,1×−1,1. As shown in Figure 13, we obtain the results of shrinking while maintaining its initial shape in Figure 13(b), and then merges and shrinks in Figures 13(c) and 13(d). If we use different ϵm, it can shrink while preserving its initial shape.Figure 13
Time evolution of a maze shape onΩ=−1,1×−1,1×−1,1, and T=0.0175. (a) t = 0, (b) t = 0.005, (c) t = 0.0125, (d) t = 0.0175.
(a)(b)(c)(d)The last simulation in 3D is phase separation with random initial condition in equation (17b). In Figure 14, starting with random values with 0.1 amplitude, but over time, phase separation occurs with values −1 to 1.Figure 14
Time evolution of phase separation with a random initial condition. (a)t = 0, (b) t = 0.0025, (c) t = 0.005, (d) t = 0.01.
(a)(b)(c)(d)To estimate the error of the CPU and GPU codes, we use the following definedErr:(19)Err=1N∑t=1Navgϕgput−ϕcput2,where N is the total number of iterations, avgX is the average of elements in an array X, and t means evolution time (i.e., t=nΔt, n=1,2,…,N). We show the error results obtained by equation (19) to check the difference between the numerical results of CPU:Python (baseline) and GPU:Pytorch. In Table 4, all the errors for any cases are less than 1.0e−6.Table 4
Errors of various numerical simulations with the baseline and ours.
Initial valueDimensionRandomDumbbellCircle/sphereMazeStarTorus2D1.75×10−67.03×10−75.51×10−72.22×10−76.55×10−77.52×10−73D3.01×10−61.20×10−61.11×10−63.36×10−61.56×10−61.91×10−6According to the results in Table5, the speed gap between CPU and GPU operations is significant. Of course, it will be faster by performing GPU calculations, but we proposed a structure using padding and convolution operation in performing GPU calculations on AC equations. And, the results are demonstrated by verifying the results (Table 4) with Python code which has the same mathematical meaning. In 3D, the GPU performance is much more overwhelming than the CPUs. For instance, the GPU code is 4766 times faster than the Python code in the torus problem. Also, GPU tensors make the model up to 76 times faster than CPU tensors in the same code. The values in parentheses describe how many times the difference is based on GPU:Pytorch time for each test.Table 5
Runtime result in 3D (sec).
IterationsInitial valueSphereDumbbellStarTorusMazeRandom200120012001120124012001CPU:Python4022.38 (3944)8050.45 (4087)4014.83 (4015)2478.22 (4766)4844.24 (3814)4163.17 (4042)CPU:Pytorch65.95 (65)79.23 (40)65.22 (65)39.56 (76)78.72 (62)65.80 (64)GPU:Pytorch1.021.971.000.521.271.03The values in parentheses describe how many times the difference is based on GPU:Pytorch time for each test.


## 3.1. Convergence Tests

We perform convergence tests for time and space. We define the discreteL2 norm as(7)eL2d=h2∑i=1Nx∑j=1Nyeijn2,where eijn=ϕijn−ϕrefn and ϕrefn is reference solution. We use the following initial condition:(8)ϕx,y,0=0.1cos2πxcos2πy,with zero Neumann boundary condition on the computational domain 0,1×0,1.First, to show the convergence rate with respect to time, the space step sizeh=0.0025 (i.e., Nx=Ny=400) is fixed, and we use the reference solution with Nx=Ny=400 and Δt=0.001h2. Then, the test proceeds to time t=0.0001 for time steps Δt=0.2h2,0.1h2, and 0.05h2. The temporal convergence rate is defined as log2eΔtL2d/eΔt/2L2d. Table 1 shows that our proposed algorithm has first-order accuracy for time.Table 1
TemporalL2 errors and convergence rates at t=0.0001.
Δt0.2h20.1h20.05h2L2 error9.282e − 34.759e − 32.487e − 3Rate0.9630.936Next, we show the spatial accuracy of our proposed algorithm. To show the spatial convergence rate, the time step was fixed toΔt=0.1h2 for each mesh size h=0.01,0.005, and 0.0025. For the reference solution, the time step was set to 0.001h2 for each mesh size. The spatial convergence rate is defined as log2ehL2d/eh/2L2d. Table 2 shows that our proposed algorithm has between second- and third-order for the spatial accuracy.Table 2
SpatialL2 errors and convergence rates at t=0.0001.
h0.010.0050.0025L2 error1.011e − 46.722e − 044.760e − 3Rate2.7332.824


## 3.2. Energy Dissipation and Maximum Principle

AC equation (1) follows the energy dissipation law, and the equation is derived from(9)Eϕ=∫ΩFϕϵ2+12∇ϕ2dx.Eϕ is decreasing with time(10)∂∂tEϕ=−∫Ω∂ϕ∂t2dx≤0.To check a discrete energy, we can rewrite equation (9) as Edϕ. Figure 2 shows mesh plots and a discrete energy graph when time evolution is performed with random initial condition. The tested initial condition is used as equation (17a) and other parameters are used Nx = Ny = 200, space step size h = 1/Nx, time step size Δt = 0.1h2, and thickness of transition layer ϵ10 on the computational domain Ω = (0, 1)2.Figure 2
(a–d) Mesh plots at the marked points from left to right of the energy graph (e). (e) Time-dependent normalized discrete total energyEdϕt/Edϕ0. (f) Maximum and minimum values of ϕ over time evolution.
(a)(b)(c)(d)(e)(f)Figure 3
(a–d) The changes in the circle over time are shown, and each time is described in the figures. (e) Change of radius for variousϵm.
(a)(b)(c)(d)(e)Also, we can obtain that the AC equation does not conserve the initial mass.(11)ddt∫Ωϕdx=∫Ωϕtdx=∫Ω−F′ϕϵ2+Δϕdx=−∫ΩF′ϕϵ2dx≤0,with zero Neumann boundary condition.


## 3.3. Initial Conditions

We consider the 2D and 3D initial conditions introduced in this section: to check theϵm effect, we measure the circle’s (sphere in 3D) radius that changes with time and take the initial conditions:(12a)ϕx,y,0=tanhR0−x−0.52+y−0.522ϵ,in2D,(12b)ϕx,y,z,0=tanhR0−x−0.52+y−0.52+z−0.522ϵ,in3D,where R0 is the initial radius of a circle (sphere in 3D).In various tests to check the motion by the mean curvature, the initial conditions for the circle and sphere refer to equations (12a) and (12b), respectively. A dumbbell is constructed by following the initial conditions (o/w: otherwise)(13a)ϕx,y,0=1.0,if0.4<x<1.6and0.4<y<0.6,1+tanhR0−x−0.32+Y2ϵ+tanhR0−x−1.72+Y2ϵ,otherwise,in2D,(13b)ϕx,y,z,0=1.0,if0.4<x<1.6and0.4<y,z<0.6,1+tanhR0−x−0.32+YZ2ϵ+tanhR0−x−1.72+YZ2ϵ,otherwise,in3D.where R0 is the initial radius of both sides of dumbbell’s circle (sphere in 3D) and for simplicity of expression, Y=y−0.52 and YZ=y−0.52+z−0.52.The initial conditions of a star shape are defined as(14a)ϕx,y,0=tanh0.25+0.1cos6θ−x−0.52+y−0.522ϵ,in2D,(14b)ϕx,y,z,0=tanh0.7+0.2cos6θ−x2+y2+z22ϵ,in3D,where(15)θ=tan−1y−0.5x−0.5,ifx>0.5,π+tan−1y−0.5x−0.5,otherwise,in2D,θ=tan−1zx,ifx>0.5,π+tan−1zx,otherwise.in3D.and we use different domain sizes in 2D and 3D, so the center of the star depends on the dimensions.And, a torus shape is given by(16a)ϕx,y,0=−1+tanhR1−XY2ϵ−tanhR2−XY2ϵ,in2D,(16b)ϕx,y,z,0=z2+x2+y2−R12−R2,in3D,where R1 and R2 are the radius of major (outside) and minor (inside) circles, respectively. And, for simplicity of expression, XY=x−0.52+y−0.52.The initial conditions of a maze shape are complicated to describe its equation, so refer to the code (https://github.com/kimy-de/gpuallencahn, https://sites.google.com/view/yh-choi/code).And, the last initial condition is random for confirming the phase separation of the AC equation:(17a)ϕx,y,0=0.1randx,y,in2D,(17b)ϕx,y,z,0=0.1randx,y,z,in3D.Here, the function randx,y has a random value between −1 and 1.


## 3.4. Simulations in the 2-Dimensional Space

Unless otherwise stated, we use the following parameters: mesh sizeNx=Ny=200, space step size h=1/Nx, time step size Δt=0.1h2, and computational domain Ω=0,1×0,1. We should find a proper thickness of the transition layer as defined in equation (4). First, we find an appropriate ϵm by comparing the numerical solution with the exact solution for the radius that decreases due to the motion by mean curvature in the circle initial condition in equation (12a).Figures3(a)–3(d) show a circle shrinking with motion by mean curvature based on ϵ10. The exact solution of the radius r decreasing with time evolution can be calculated [22]:(18)rt=r02+21−dt,where r0 is the initial radius of the circle, d is the dimension, and t is time. We simulate various ϵm until the final time T=0.03. As shown in Figure 3, when ϵ10 is used, it is found that the exact solution and the numerical solution are most similar in the experiment. Therefore, the other tests are performed using ϵ10 except for the case of changing the grid in 2D.Figure4 shows the temporal evolution of the dumbbell shape, and the initial condition is shown in equation (13a). For resolution, we use Nx=400 and Ny=200 on the computational domain Ω=0,2×0,1. The final time T of the simulation is 0.0094 and R0=0.2. In Figure 4(e), the changing direction by the motion by mean curvature is indicated by arrows.Figure 4
(a–d) Time evolution of dumbbell shape, and each time is described in the figures. (e) Contour lines over time are shown by overlapping. Changing by the mean curvature flow can be seen.
(a)(b)(c)(d)(e)Figure5 shows the evolution of the star shape created by equation (14a). The parameters are used as mentioned at the beginning of this section and T=0.0325. As shown in Figure 5, the tips of the star move inward and the gaps between the tips move outward. When it changes to the shape of a circle, the change of the radius can be predicted as shown in Figure 3.Figure 5
Time evolution of a star shape. As in the shape of a dumbbell, it can be seen that it changes with the mean curvature flow. (a)t = 0, (b) t = 0.0025, (c) t = 0.0188, (d) t = 0.0325.
(a)(b)(c)(d)Figure6 shows the evolution of the torus shape of equation (16a) with T=0.0575, R1=0.4, and R2=0.3. Because the inner circle has a larger curvature, it shrinks faster than the outer circle, and after the inner circle disappears, the change of radius over time can be measured as shown in Figure 3.Figure 6
Time evolution of a torus shape, and each time is described in the figures. The inner circle has larger curvature than the outer circle, so the inner circle shrinks faster than the outer circle. (a)t = 0, (b) t = 0.0188, (c) t = 0.0375, (d) t = 0.0575.
(a)(b)(c)(d)Figure7 shows the evolution of a maze shape. We use Nx=Ny=100, ϵ5, and T=0.04. As shown in Figure 7, we obtain the results of shrinking while maintaining its initial shape.Figure 7
Time evolution of a maze shape on100×100 mesh size with ϵ5, and T=0.04. (a) t = 0, (b) t = 0.01, (c) t = 0.025, (d) t = 0.04.
(a)(b)(c)(d)The last simulation in 2D is phase separation with a random initial condition equation (17a). In Figure 8, starting with random values with 0.1 amplitude, but over time, phase separation occurs with values −1 to 1.Figure 8
Time evolution of phase separation with a random initial condition. (a)t = 0, (b) t = 0.0025, (c) t = 0.0038, (d) t = 0.0125.
(a)(b)(c)(d)According to the results in Table3, the speed gap between CPU and GPU is significant. In 2D, it is up to 251.6 times the difference between the Python:CPU and the Pytorch:GPU codes. Also, Pytorch:GPU tensors make the model up to 4.73 times faster than Pytorch:CPU tensors in the same code.Table 3
Runtime Result in 2D(sec).
IterationsInitial valueCircleDumbbellStarTorusMazeRandom12001150011300123001400112001CPU:Python680.09 (185.31)1736.05 (251.6)729.74 (128.47)1326.82 (165.65)60.15 (52.3)674.63 (119.62)CPU:Pytorch14.73 (4.01)30.67 (4.44)17.60 (3.1)37.85 (4.73)2.98 (2.59)17.79 (3.15)GPU:Pytorch3.676.905.688.011.155.64The values in parentheses describe how many times the difference is based on GPU:Pytorch time for each test.


## 3.5. Simulations in the 3-Dimensional Space

Three-dimensional simulations can be considered as an extension of two-dimensional tests. Unless otherwise stated, we use the following parameters: mesh sizeNx=Ny=Nz=100, space step size h=1/Nx, time step size Δt=0.1h2, and computational domain Ω=0,1×0,1×0,1. As in the 2D simulation, we find an appropriate ϵm by comparing the numerical solution with the exact solution (using in equation (18)) for the radius of the sphere initial condition in equation (12b).Figures9(a)–9(d) show the time evolution results for ϵ12, and (e) presents the comparison of the numerical and exact solutions of radius for various epsilons. As shown in Figure 9, when ϵ12 is used, it is found that the exact solution and the numerical solution are most similar. Therefore, the other tests are performed using ϵ12.Figure 9
(a–d) The changes in the circle over time are shown, and each time is described in the figures. (e) Change of radius for variousϵm.
(a)(b)(c)(d)(e)Figure10 shows the temporal evolution of the dumbbell shape, and the initial condition is shown in equation (13b). For resolution, we use Nx=200 and Ny=Nz=100 on the computational domain Ω=0,2×0,1×0,1. The final time T of the simulation is T=0.0025 and R0=0.25. The tendency of the 3D dumbbell motion is different from the result of 2D dumbbell. The reason is that, in the initial shape in Figure 10(a), the radius of the handle is much smaller than the radius of the spheres at both ends. Therefore, the handle part shrinks faster than the end of spheres and breaks as in Figure 10.Figure 10
(a–d) Time evolution of the dumbbell shape. The handle shrinks quickly and breaks occur because the curvature of the handle is larger than that of both spheres. (a)t = 0, (b) t = 0.0006, (c) t = 0.0013, (d) t = 0.0025.
(a)(b)(c)(d)Figure11 shows the evolution of the star shape shown in equation (14b) on the computational domain Ω=−1,1×−1,1×−1,1. The parameters are used as mentioned at the beginning of this section and T=0.02. As shown in Figure 11, similar to the 2D result, the tips of the star move inward, and the gaps between the tips move outward.Figure 11
Time evolution of a star shape. The shape changes with the mean curvature flow. (a)t = 0, (b) t = 0.0025, (c) t = 0.0075, (d) t = 0.02.
(a)(b)(c)(d)Figure12 shows the evolution of the torus shape with initial condition in equation (16b) using R1=0.3 and R2=0.3 and T=0.01 on the computational domain Ω=−1,1×−1,1×−1,1. Contrary to the 2D result, the inner circle becomes larger because the radius of the minor circle exists. Therefore, the radius of the major circle and the minor circle each has different curvature. Since the radius of the minor circle is smaller than major circle, the mean curvature drives the motion into the inside of the torus as shown in Figure 12.Figure 12
Time evolution of a torus shape. Contrary to the 2D results, the inner circle increases due to the effect of the mean curvature. (a)t = 0, (b) t = 0.005, (c) t = 0.008, (d) t = 0.01.
(a)(b)(c)(d)Figure13 shows the evolution of a maze shape. The initial condition is described in code (https://github.com/kimy-de/gpuallencahn, https://sites.google.com/view/yh-choi/code). In this simulation, we use T=0.0175 and computational domain Ω=−1,1×−1,1×−1,1. As shown in Figure 13, we obtain the results of shrinking while maintaining its initial shape in Figure 13(b), and then merges and shrinks in Figures 13(c) and 13(d). If we use different ϵm, it can shrink while preserving its initial shape.Figure 13
Time evolution of a maze shape onΩ=−1,1×−1,1×−1,1, and T=0.0175. (a) t = 0, (b) t = 0.005, (c) t = 0.0125, (d) t = 0.0175.
(a)(b)(c)(d)The last simulation in 3D is phase separation with random initial condition in equation (17b). In Figure 14, starting with random values with 0.1 amplitude, but over time, phase separation occurs with values −1 to 1.Figure 14
Time evolution of phase separation with a random initial condition. (a)t = 0, (b) t = 0.0025, (c) t = 0.005, (d) t = 0.01.
(a)(b)(c)(d)To estimate the error of the CPU and GPU codes, we use the following definedErr:(19)Err=1N∑t=1Navgϕgput−ϕcput2,where N is the total number of iterations, avgX is the average of elements in an array X, and t means evolution time (i.e., t=nΔt, n=1,2,…,N). We show the error results obtained by equation (19) to check the difference between the numerical results of CPU:Python (baseline) and GPU:Pytorch. In Table 4, all the errors for any cases are less than 1.0e−6.Table 4
Errors of various numerical simulations with the baseline and ours.
Initial valueDimensionRandomDumbbellCircle/sphereMazeStarTorus2D1.75×10−67.03×10−75.51×10−72.22×10−76.55×10−77.52×10−73D3.01×10−61.20×10−61.11×10−63.36×10−61.56×10−61.91×10−6According to the results in Table5, the speed gap between CPU and GPU operations is significant. Of course, it will be faster by performing GPU calculations, but we proposed a structure using padding and convolution operation in performing GPU calculations on AC equations. And, the results are demonstrated by verifying the results (Table 4) with Python code which has the same mathematical meaning. In 3D, the GPU performance is much more overwhelming than the CPUs. For instance, the GPU code is 4766 times faster than the Python code in the torus problem. Also, GPU tensors make the model up to 76 times faster than CPU tensors in the same code. The values in parentheses describe how many times the difference is based on GPU:Pytorch time for each test.Table 5
Runtime result in 3D (sec).
IterationsInitial valueSphereDumbbellStarTorusMazeRandom200120012001120124012001CPU:Python4022.38 (3944)8050.45 (4087)4014.83 (4015)2478.22 (4766)4844.24 (3814)4163.17 (4042)CPU:Pytorch65.95 (65)79.23 (40)65.22 (65)39.56 (76)78.72 (62)65.80 (64)GPU:Pytorch1.021.971.000.521.271.03The values in parentheses describe how many times the difference is based on GPU:Pytorch time for each test.


## 4. Conclusions

In this paper, we proposed a structure using padding and convolution operation for the GPU calculation of the Allen–Cahn equation. We increased the simulation speed and demonstrated the validity of the proposed structure by verifying the result with the Python code that has the same mathematical meaning. We solved the Allen–Cahn equation by the explicit finite difference method and compared the runtime results between CPU and GPU algorithms. The errors of CPU:Python and GPU:Pytorch are less than1.0e−6 for the given initial conditions. Also, we showed that our GPU code is up to 251.6 and 4765.81 times faster than the CPU codes in 2D and 3D, respectively. By showing these results, accuracy and efficiency have been demonstrated. Various numerical simulations were presented to confirm that the Allen–Cahn equation follows the motion by mean curvature and phase separation in the 2D and 3D space. In this way, we can build a fast algorithm for any differential equations using the finite difference method by efficient programmatic code structures.



---

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