
In this paper, we propose a novel, simple, efficient, and explicit numerical method for the Allen–Cahn (AC) equation on effective symmetric triangular meshes. First, we compute the net vector of all vectors starting from each node point to its one-ring neighbor vertices and virtually adjust the neighbor vertices so that the net vector is zero. Then, we define the values at the virtually adjusted nodes using linear and quadratic interpolations. Finally, we define a discrete Laplace operator on triangular meshes. We perform several computational experiments to demonstrate the performance of the proposed numerical method for the Laplace operator, the diffusion equation, and the AC equation on triangular meshes.
Citation: Youngjin Hwang, Seokjun Ham, Chaeyoung Lee, Gyeonggyu Lee, Seungyoon Kang, Junseok Kim. A simple and efficient numerical method for the Allen–Cahn equation on effective symmetric triangular meshes[J]. Electronic Research Archive, 2023, 31(8): 4557-4578. doi: 10.3934/era.2023233
[1] | Youngjin Hwang, Ildoo Kim, Soobin Kwak, Seokjun Ham, Sangkwon Kim, Junseok Kim . Unconditionally stable monte carlo simulation for solving the multi-dimensional Allen–Cahn equation. Electronic Research Archive, 2023, 31(8): 5104-5123. doi: 10.3934/era.2023261 |
[2] | Junseok Kim . Maximum principle preserving the unconditionally stable method for the Allen–Cahn equation with a high-order potential. Electronic Research Archive, 2025, 33(1): 433-446. doi: 10.3934/era.2025021 |
[3] | Yan Yang, Xiu Ye, Shangyou Zhang . A pressure-robust stabilizer-free WG finite element method for the Stokes equations on simplicial grids. Electronic Research Archive, 2024, 32(5): 3413-3432. doi: 10.3934/era.2024158 |
[4] | Kai Jiang, Wei Si . High-order energy stable schemes of incommensurate phase-field crystal model. Electronic Research Archive, 2020, 28(2): 1077-1093. doi: 10.3934/era.2020059 |
[5] | Yao Yu, Guanyu Xue . A nonlinear correction finite volume scheme preserving maximum principle for diffusion equations with anisotropic and discontinuous coefficient. Electronic Research Archive, 2025, 33(3): 1589-1609. doi: 10.3934/era.2025075 |
[6] | Shuhao Cao . A simple virtual element-based flux recovery on quadtree. Electronic Research Archive, 2021, 29(6): 3629-3647. doi: 10.3934/era.2021054 |
[7] | Hao Wu . A review on the Cahn–Hilliard equation: classical results and recent advances in dynamic boundary conditions. Electronic Research Archive, 2022, 30(8): 2788-2832. doi: 10.3934/era.2022143 |
[8] | Wenyan Tian, Yaoyao Chen, Zhaoxia Meng, Hongen Jia . An adaptive finite element method based on Superconvergent Cluster Recovery for the Cahn-Hilliard equation. Electronic Research Archive, 2023, 31(3): 1323-1343. doi: 10.3934/era.2023068 |
[9] | Jia-Min Luo, Hou-Biao Li, Wei-Bo Wei . Block splitting preconditioner for time-space fractional diffusion equations. Electronic Research Archive, 2022, 30(3): 780-797. doi: 10.3934/era.2022041 |
[10] | Derrick Jones, Xu Zhang . A conforming-nonconforming mixed immersed finite element method for unsteady Stokes equations with moving interfaces. Electronic Research Archive, 2021, 29(5): 3171-3191. doi: 10.3934/era.2021032 |
In this paper, we propose a novel, simple, efficient, and explicit numerical method for the Allen–Cahn (AC) equation on effective symmetric triangular meshes. First, we compute the net vector of all vectors starting from each node point to its one-ring neighbor vertices and virtually adjust the neighbor vertices so that the net vector is zero. Then, we define the values at the virtually adjusted nodes using linear and quadratic interpolations. Finally, we define a discrete Laplace operator on triangular meshes. We perform several computational experiments to demonstrate the performance of the proposed numerical method for the Laplace operator, the diffusion equation, and the AC equation on triangular meshes.
In this study, we present a novel explicit numerical method for solving the Allen–Cahn equation on triangular meshes. The method is characterized by its simplicity and efficiency. The AC equation, a partial differential equation (PDE), is employed in mathematical physics and materials science to represent the dynamics of phase transitions and pattern formation in binary systems [1]. The equation incorporates a free energy density function that characterizes the potential energy of the system and a linear diffusion term that drives the system toward a state of minimum energy [2]. In two-dimensional (2D) space, the AC equation with Dirichlet boundary condition on the domain Ω as follows
∂u(x,y,t)∂t=−F′(u(x,y,t))ϵ2+Δu(x,y,t),(x,y)∈Ω,t>0,u(x,y,t)=f(x,y) on ∂Ω, | (1.1) |
where u(x,y,t) is the phase-field function, ϵ is an interface transition layer related parameter, F(u)=0.25(u2−1)2 is a double-well potential energy functional, and Δ is a Laplace operator.
Choi et al. [1] developed an unconditionally gradient-stable numerical system for the AC equation on uniformly square grid and showed the same property for the corresponding discrete problem by using the eigenvalues of the Hessian matrix of the energy functional. In [3], the authors developed a simple and practical adaptive finite difference scheme for AC equations in 2D. The developed method uses a temporally adaptive narrow band domain in a uniform discrete rectangular domain. Xiao et al. [4] introduced the lumped mass surface finite element method on a triangular mesh to solve the surface AC equation. To maintain mass conservation, Rubinstein and Sternberg [5] introduced a conservative AC equation by incorporating a Lagrange multiplier into the AC equation. Construction of the conservative AC equation addresses the issue of the AC equation inability to conserve total mass. Sun and Zhang [6] developed a meshless radial basis function approach to solve the conservative AC equation on smooth compact surfaces embedded in R3. The developed method ensures mass conservation while employing an operator splitting scheme. The Laplace–Beltrami operator (LBO) is approximated iteratively by using radial basis functions. In addition, the time discretization of the diffusion equation is achieved using the Euler method. Xia et al. [7] developed a second-order accurate method in time and space for the AC and conservative AC equation using triangular discretization. The AC and conservative AC equation can also be solved using the explicit hybrid numerical method, as studied by Choi et al. [8]. Joshi and Jaiman [9] proposed an adaptive variational procedure for the AC equation with a two-phase field. This scheme is stable, robust, and general for the AC equation on a triangular mesh. Kwak et al. [10] proposed the conservative AC equation with a curvature-dependent Lagrange multiplier. Hong et al. [11] developed numerical algorithms with fully discrete structure-preserving of arbitrarily high order for the AC model with a nonlocal constraint subject to the Neumann boundary condition. Chai et al. [12] proposed a simple multiple-relaxation-time lattice-Boltzmann method for the local and nonlocal AC equation, presenting the property of mass conservation through numerical experiments. Inc et al. [13] performed Lie symmetry analysis, explicit solutions, and convergence analysis for the time-fractional Cahn–Allen equation. Kim et al. [14] proposed a temporally second-order unconditionally energy-stable computational method for the AC equation with high-order polynomial free energy potential. Bhatt et al. [15] devised an attractive and easy-to-implement alternative for integrating the multi-dimensional AC equation with no-flux boundary conditions by combining the Fourier spectral method with the strongly stable exponential time difference method. The method they developed has several key advantages, including the use of discrete fast Fourier transform for efficiency, the ability to extend to two and three spatial dimensions in a manner similar to 1D problems, and the capability to handle various boundary conditions. Sun and Gao [16] proposed a high order multiquadric trigonometric quasi-interpolation method for function approximation and derivative approximation based on periodic sampling data. Furthermore, the method was applied to solve different types of PDEs including the AC equation. Sun and Ling [17] developed a meshless conservative Galerkin method for solving Hamiltonian wave equations, which is a classical example of Hamiltionian PDEs. The equation was first discretized in space using radial basis functions in a Galerkin-type formulation. Their method used appropriate projection operators for the construction of the Galerkin equation and was shown to conserve global energies.
The AC equation without a nonlinear term is a diffusion equation that is a linear differential equation:
∂u(x,y,t)∂t=DΔu(x,y,t), | (1.2) |
where D is the diffusion coefficient parameter and we set D=1 in this paper for simplicity of exposition.
There exist meshes that do not satisfy the properties of the continuous Laplace operator for continuous functions using the discrete Laplace operator, and it is more difficult to satisfy the properties of the continuous Laplace operator for triangular meshes [18]. The method we propose requires defining a discrete Laplace operator to solve the diffusion term. Thus, a simple and efficient algorithm is included that numerically satisfies the simple properties of the Laplace operator.
Triangular meshes are widely used for representing complex domains or curved surfaces. Until recently, the discrete Laplace operator for triangular meshes has been a subject of numerous studies. Xu [19] proposed discretization schemes of the Laplace-Beltrami operator on triangulated surfaces with curvatures. Meanwhile, Caissard et al. [20] presented a new discretization of the Laplace-Beltrami operator over digital surfaces. Thampi et al. [21] demonstrated schemes for isotropic discrete Laplacian operators based on lattice hydrodynamics. McCartin [22] studied the discrete Laplacian on an equilateral triangle with both Dirichlet and Neumann boundary conditions using eigenvalues and eigenfunctions of the continuous Laplacian. A finite difference algorithm for solving the generalized Laplace equation on unstructured triangular meshes was developed by Ganzha et al. [23] using a support operator method combining discrete divergence and gradient operators. Yoon et al. [24] conducted numerical experiments on triangular meshes to investigate dendritic pattern formation in an isotropic crystal growth model on surfaces, and subsequently presented their findings. Paquet and Korikache [25] established an a priori error estimate for the fully discretized problem of the dual mixed method for the non-stationary heat equation in 2D polygonal domains. The Cahn–Hilliard (CH) equation is one of the famous phase-field models along with the AC equation [26]. Tian et al. [27] developed an adaptive finite element method based on superconvergent cluster recovery to solve the Laplace operator of the CH equation on a triangular mesh. Therefore, the construction of simple and efficient algorithms for solving Laplace operators on triangular meshes is important. We propose a numerical method involving the novel discrete Laplace operator and verify the proposed discrete Laplace operator through numerical experiments, including comparisons to general discrete Laplace operators on rectangular meshes or the LBO on triangular meshes.
The contents of this paper are as follows. In Section 2, the proposed numerical scheme is described. In Section 3, we present the numerical experiments. In Section 4, the conclusion is drawn.
In this section, we present a numerical solution algorithm. The proposed algorithm solves the AC equation by using an explicit Euler method. The numerical algorithm requires computing a linear term with the Laplace operator. We propose a novel Laplace operator in triangular meshes to numerically compute the AC equation. Firstly, we introduce the discrete LBO on triangular mesh and propose a novel discrete Laplace operator. Because the plane is a special case of the curved surfaces, let us consider the discretization of LBO [19] on curved surfaces. Let N be the number of node points on the triangular mesh, xk=(xk,yk), k=1,2,…,N be node points on the triangular mesh, Δt be temporal step size, and unk be an approximation of phase-field function u(xk,nΔt). In Figure 1, xk denotes k-th node point of a triangular mesh and {xk1,xk2,xk3,xk4,xk5} are one-ring node points neighboring xk. For simplicity, we assume that there are five node points neighboring xk. The discrete LBO at node point xk is given as
Δuk≈3A(xk)Nk∑m=1cotαkm+cotβkm2(ukm−uk), |
where A(xk) is the sum of areas of triangles Tm sharing the node point xk, Nk is number of node points neighboring xk, xk0=xkNk and xkNk+1=xk1. We also define angles αkm=∠xkxkm−1xkm and βkm=∠xkxkm+1xkm for m=1,⋯,Nk. For example, when m=5, αk5=∠xkxk4xk5 and βk5=∠xkxk6xk5 in triangles T5 and T6, respectively, as shown in Figure 1.
Now, we propose the novel discrete Laplace operator. Let N be the total number of node points in the triangular mesh and Nk be the number of one-ring neighboring node points of each node xk, for k=1,2,…,N in the triangular mesh. The neighboring node points of each node point xk are denoted as xkm, m=1,2,…,Nk. We define the net force Fk for each node point xk, k=1,2,…,N, which plays a crucial role in achieving symmetry within the proposed discrete Laplace operator.
Fk=Nk∑m=1(xkm−xk). |
The net force at each points in a symmetric object or system is zero because the forces acting on the object are symmetrically distributed. Additionally, the net force of the general five and nine stencil Laplace operators is 0. Therefore, we propose the following sub-algorithm using virtual nodes such that the net force Fk of each node xk, k=1,2,…,Nk becomes 0. We solve the AC equation on the irregular domain by using a sub-algorithm using virtual nodes to make each node point symmetric for the asymmetric triangular mesh generated on the irregular domain in 2D space.
The sub-algorithm of the proposed method is illustrated in Figure 2. As presented in Figure 2, for each node xk, one-ring neighboring node points are denoted by xk1, xk2, …, xk6. First, find all the one-ring neighboring node points for a given point xk. Figure 2(a) presents a given point xk and one-ring neighboring node points. Second, find the net force Fk and the triangular domain containing Fk. Here, we use the following method to find the triangle. In the case of Figure 2(b), let the vectors from xk of the found triangle to the other two node points xk3 and xk2 be u and v, respectively. Therefore, net force Fk is the linear combination of u and v with constants a and b as follows:
au+bv=F,a,b∈R, |
where u and v are 2D vectors. Let u=(u1,u2)T and v=(v1,v2)T. For the third step, we consider the following the matrix form to find values of a and b.
(u1v1u2v2)(ab)=F⇒(ab)=(u1v1u2v2)−1F. |
Then, subtract au and bv from u and v to reduce the net force Fk to 0 at xk. Let x∗k3 and x∗k2 be virtual node points corresponding to u−au and v−bv, respectively. In the context of the triangular mesh, it is important to note that the physical positions of the node points remain unchanged. Instead, a theoretical construct of imaginary node points is utilized to ensure that the net force acting on the mesh is balanced to zero. By introducing these virtual node points, it is possible to accurately calculate the net force and achieve a state of equilibrium, without the need to manipulate the physical node points. Fourth, we use interpolation for the points where the extension of the modified and original meshes meet on the opposite side of u and v, respectively. In Figure 2(e), the value of x∗∗k2 is obtained by the linear interpolation of the values of xk5 and xk6 as follows:
u(x∗∗k2)=u(xk6)−u(xk5)√(xk6−xk5)2+(yk6−yk5)2√(x∗∗k2−xk5)2+(y∗∗k2−yk5)2+u(xk5). |
Here, the virtual node point x∗∗k2 is the intersection point of the line passing through xk and xk2, and the line passing through xk5 and xk6. In addition, x∗∗k3 is obtained analogously. The net force at a point is zero in a symmetric object or system, as the forces acting on it are distributed symmetrically, ensuring that they cancel out each other. Therefore, when virtual node points are used, each node point of the triangular mesh becomes symmetrical with respect to one-ring neighboring points. Finally, the values of modified mesh points are obtained by quadratic interpolation. In Figure 2(f), the value of virtual node x∗k3 is obtained by the quadratic interpolation using the values of x∗∗k3, x3, and xk3. The value of virtual node x∗k2 is similarly obtained. Next, we define the discrete Laplace operator with the values of the neighboring node points set such that the net force Fk at a given point is zero. Let xkm for m=1,…,Nk be the one-ring neighboring node points at a point xk with net force Fk=0, then we propose the discrete Laplace operator as follows:
˜Δuk=∑Nkm=1(u(xkm)−u(xk))Sk, |
where Sk is a constant given for each xk. The general Laplace operator satisfies Δf(x,y)=1 if f(x,y)=0.25(x2+y2). Because the proposed Laplace operator should satisfy the characteristics of the general Laplace operator, we define Sk as
Sk=Nk∑m=1f(dkm). |
Here, f(x,y)=0.25(x2+y2) and dkm=xkm−xk. The reason for using this function is that it is a symmetric function with respect to the x- and y-axes, and unlike a linear function, a non-zero solution exists, and unlike a cubic function, the solution is uniquely determined. Therefore, the proposed discrete Laplace operator satisfies the characteristic of general Laplace operator, u(x,y)=0.25(x2+y2). Then we discretize Eq (1.1) in temporal and spatial as follows
un+1k−unkΔt=−F′(unk)ϵ2+˜Δunk=−F′(unk)ϵ2+∑Nkm=1(unkm−unk)Sk, | (2.1) |
where Sk=∑Nkm=1f(dkm). Therefore, we obtain the following equation by rearranging Eq (2.1) for un+1k.
un+1k=unk+Δt(−F′(unk)ϵ2+∑Nkm=1(unkm−unk)Sk). | (2.2) |
We consider the general five stencil, nine stencil and exact Laplace operators to compare the proposed Laplace operator with the conventional Laplace operators.
Δ1uij=ui−1,j+ui+1,j+ui,j−1+ui,j+1−4uijh2, | (3.1) |
Δ2uij=ui−1,j+1+ui+1,j+1+ui+1,j−1+ui−1,j−1−4uij2h2, | (3.2) |
Δ3uij=23Δ1uij+13Δ2uij=4(ui−1,j+ui+1,j+ui,j−1+ui,j+1)6h2+ui−1,j+1+ui+1,j+1+ui+1,j−1+ui−1,j−1−20uij6h2, | (3.3) |
which are the well-known standard discrete Laplace operators. Figures 3(a), (b), (c) represent the numerical stencils for Eqs (3.1), (3.2) and (3.3), respectively.
In this section, each discrete Laplace operator is compared with the C2 function to verify the performance of the proposed discrete Laplace operator. The proposed discrete Laplace operator ˜Δ at each node is made by neighboring nodes as given in Figure 3(c). This is one of the benefits of our proposed discrete Laplace operator on a triangular mesh. Therefore, we obtain the discrete Laplace equation as follows using the proposed method.
˜Δuij=(ui+1,j+ui−1,j+ui,j+1+ui,j−1+ui+1,j+1+ui+1,j−1+ui−1,j+1+ui−1,j−1−8uij)/3h2. |
We consider the following function on the domain Ω=(0,1)×(0,1):
u(x,y)=sin(πx)sin(πy),(x,y)∈Ω | (3.4) |
with zero Dirichlet boundary condition. Then, the exact solution of Δu is
Δu(x,y)=−2π2sin(πx)sin(πy). |
We use different parameter values Nx=9,17,33,65 to consider the convergence rates for different spatial mesh hx=hy=1/(Nx−1) and discrete domain xi=(i−1)hx for i=1,2,…,Nx and yj=(j−1)hy for j=1,2,…,Ny. Next, to calculate the convergence rate for Δuij, we define the discrete l2-error and discrete maximum-error in 2D space, respectively:
‖eNxNy‖l2=√1Nx1NyNx∑i=1Ny∑j=1(eij)2,‖eNxNy‖max=max1≤i≤Nx,1≤j≤Ny|eij|, |
where eij=ˆΔuij−Δuij for i=1,2,…,Nx and j=1,2,…,Ny. Figures 4(a), 4(b) present the l2-error and the maximum-error for different Laplace operators, respectively. We make a list of the convergence rates for the maximum-error and l2-error of the Laplace operator in Tables 1 and 2, respectively. We observed that the proposed discrete Laplace operator achieved second-order accuracy in spatial directions on a rectangular mesh. The proposed discrete Laplace operator can consider triangular meshes on irregular domains rather than rectangular meshes, while other discrete Laplace operators can only be used on rectangular meshes.
spatial step | h(=1/8) | rate | h/2 | rate | h/4 | rate | h/8 |
Δ1u | 2.5237×10−1 | 1.99 | 6.3336×10−2 | 2.00 | 1.5849×10−2 | 2.00 | 3.9633×10−3 |
Δ2u | 9.9404×10−1 | 1.98 | 2.5237×10−1 | 1.99 | 6.3336×10−2 | 2.00 | 1.5849×10−2 |
Δ3u | 4.9959×10−1 | 1.98 | 1.2635×10−1 | 2.00 | 3.1678×10−2 | 2.00 | 7.9253×10−3 |
˜Δu | 7.4682×10−1 | 1.98 | 1.8936×10−1 | 1.99 | 4.7507×10−2 | 2.00 | 1.1887×10−2 |
spatial step | h(=1/8) | rate | h/2 | rate | h/4 | rate | h/8 |
Δ1u | 1.1216×10−1 | 1.91 | 2.9805×10−2 | 1.96 | 7.6845×10−3 | 1.98 | 1.9511×10−3 |
Δ2u | 4.4180×10−1 | 1.90 | 1.1876×10−1 | 1.95 | 3.0708×10−2 | 1.98 | 7.8027×10−3 |
Δ3u | 2.2204×10−1 | 1.90 | 5.9457×10−2 | 1.95 | 1.5359×10−2 | 1.98 | 3.9017×10−3 |
˜Δu | 3.3192×10−1 | 1.90 | 8.9110×10−2 | 1.95 | 2.3034×10−2 | 1.98 | 5.8522×10−3 |
The proposed method is for 2D triangular mesh, therefore we need to define some conditions. Let N be the number of nodes in the triangular mesh and xk=(xk,yk) be k-th node in the triangular mesh for k=1,…,N. Set I be the set of indices of the interior nodes and Nk is the number of node points neighboring xk. The average of spatial mesh is
have=∑k∈I∑Nkm=1|xkm−xk|∑k∈INk, | (3.5) |
where xkm is neighboring node points of xk for m=1, …, Nk.
Let the initial function u(x,y)=ax+by+c for a,b,c∈R be a plane, then Δu=0. Therefore, we consider the initial condition on the domain Ω=(0,1)×(0,1):
u(x,y)=x+2y+3,(x,y)∈Ω. | (3.6) |
Here, we apply Dirichlet boundary conditions. To compute the discrete Laplace operator of u, we use a triangular mesh with have=0.1013. Figures 5(a), (b), (c) present the triangular mesh, u(x,y)=x+2y+3, and proposed discrete Laplace operator for u, respectively.
Next, we consider the initial function on the domain Ω=(0,1)×(0,1):
u(x,y)=sin(πx)sin(πy). | (3.7) |
Then we have
Δu(x,y)=−2π2sin(πx)sin(πy). | (3.8) |
We apply the zero Dirichlet boundary condition. To compute the discrete Laplace operator of u, we use the parameter have=0.1013. Figures 6(a), (b), (c) show the triangular mesh, u(x,y)=sin(πx)sin(πy), and the proposed discrete Laplace operator for u, i.e., ˜Δu(x,y), respectively.
In this section, we compare the numerical and analytical solutions to verify the proposed method. We consider a diffusion equation (1.2) with the following initial condition on the domain Ω=(0,1)×(0,1):
u(x,y,0)=sin(πx)sin(πy) | (3.9) |
with zero Dirichlet boundary condition. Then, the exact solution is
uexact(x,y,t)=sin(πx)sin(πy)e−2π2t. |
Let unk be an approximation of u(xk,yk,nΔt), where Δt is a time step. Let N be the number of nodes of the triangular mesh, xk=(xk,yk) be k-th node in triangular mesh for k=1,…,N. For this test, we use different triangular meshes with have=0.1254,0.1013,0.0762. Using the explicit Euler method and the proposed discrete Laplace operator, Eq (1.2) is as follows:
un+1k−unkΔt=˜Δunk. | (3.10) |
Let us define the l2-error and maximum-error as follows:
‖eNtN‖l2=√1NN∑k=1(eNtk)2,‖eNtN‖max=max1≤k≤N|eNtk|, |
where eNtk=uNtk−uexact(xk,yk,T), k=1,2,…,N. Figure 7 presents the triangular mesh and numerical solution using the proposed method for each have. We list the maximum-error and l2-error of the diffusion equation with T=0.02, Nt=2048, Δt=T/Nt in Table 3.
Spatial step | have=0.1254 | rate | have=0.1013 | rate | have=0.0762 |
Maximum-error | 1.8338×10−2 | 1.19 | 1.4246×10−2 | 1.92 | 8.2534×10−3 |
l2-error | 4.2098×10−3 | 1.64 | 2.9665×10−3 | 1.77 | 1.7904×10−3 |
Next, we perform numerical experiments for a larger time evolution with a final time T=0.1 using the same triangular meshes as in the previous test. Figure 8 shows the numerical solutions for different spatial steps at the final time T=0.1. Table 4 lists the maximum-error and l2-error for the proposed discrete Laplace operator. We observed that the l2-error and maximum-error of the proposed discrete Laplace operator decreased as the spatial step decreased for a larger final time T=0.1.
Spatial step | have=0.1254 | rate | have=0.1013 | rate | have=0.0762 |
Maximum-error | 4.6855×10−3 | 2.07 | 3.0121×10−3 | 1.98 | 1.7142×10−3 |
l2-error | 2.0184×10−3 | 2.89 | 1.0894×10−3 | 1.09 | 7.9989×10−4 |
We perform a convergence test on the discrete LBO to compare the discrete LBO and the proposed discrete Laplace operator. The initial conditions given are the same as in the previous convergence test with the final time T=0.1. Figure 9 shows the numerical solutions obtained using the discrete LBO for different spatial steps. Table 5 lists the maximum-error and l2-error for the discrete LBO.
Spatial step | have=0.1254 | rate | have=0.1013 | rate | have=0.0762 |
Maximum-error | 5.3441×10−3 | 2.11 | 3.4101×10−3 | –7.49 | 2.8734×10−2 |
l2-error | 2.3112×10−3 | 2.01 | 1.5061×10−3 | –7.16 | 1.1569×10−2 |
We observed that the error of the discrete LBO is larger than that of the proposed discrete Laplace operator and does not converge as the spatial step becomes smaller.
In this section, we consider a diffusion equation on a circular domain using the proposed discrete Laplace operator. We use the circular domain of radius 1 and with center at the origin, i.e., Ω={(x,y)|x2+y2<1}. We use the zero Dirichlet boundary condition and u(x,y,0)=1−x2−y2 as the initial condition.
∂u(x,y,t)∂t=Δu(x,y,t),(x,y)∈Ω,t>0,u(x,y,t)=0,(x,y)∈∂Ω,u(x,y,0)=1−x2−y2,(x,y)∈Ω. | (3.11) |
Let us consider a benchmark problem using a manufactured solution [28], i.e., we assume that the solution of the benchmark problem is
u(x,y,t)=(1−x2−y2)e−t. | (3.12) |
Then, we consider the following benchmark problem:
∂u(x,y,t)∂t=Δu(x,y,t)+g(x,y,t),(x,y)∈Ω,t>0,u(x,y,t)=0,(x,y)∈∂Ω,u(x,y,0)=1−x2−y2,(x,y)∈Ω, | (3.13) |
where g(x,y,t)=∂u(x,y,t)/∂t−Δu(x,y,t)=(3+x2+y2)e−t is the source term. To construct the numerical solution, we use explicit Euler method as follows:
un+1k=unk+Δt(˜Δuk+gn+12k), | (3.14) |
where gn+12k=(3+x2k+y2k)e−(n+12)Δt.
Let us consider three different spatial meshes with have=0.5315,0.2133 and 0.1020. The total time is T=0.1 and time step Δt=0.001h2ave with total iteration Nt=T/Δt. Using the proposed discrete Laplace operator, we present the numerical solution at final time T and compare it with the exact solution, Eq (3.12). Figure 10 presents the triangular mesh, exact solution and numerical results, and the absolute error between the results.
In this section, we consider the AC Eq (1.1) on plane with triangular mesh. Let Ω be the domain in 2D space and N be the number of nodes in domain. For k=1,…,N, let xk=(xk,yk) be the k-th node in domain. On the computational triangular mesh, let unk be an approximation of u(xk,yk,nΔt), where Δt is a time step.
In Figure 12, we verify the mean curvature flow of the AC equation using following initial condition on Ω=(−1,1)×(−1,1):
u(x,y,0)=tanh(R0−√x2+y2√2ϵ),(x,y)∈Ω, | (3.15) |
where R0=0.4 is the initial radius with Dirichlet boundary condition u(x,y,t)=−1 on ∂Ω. The analytic radius is given as R(t)=√R02−2t for time t [29]. The numerical radius is computed using the average of distances between zero level points of the numerical solution and the center point (0,0). As shown in Figure 12, the numerical radii shrink along with the analytic radius.
Next, we consider three triangular meshes generated by different spatial steps, have=0.1524, 0.0997, 0.0499. The initial condition is Eq (3.15) and interface parameter is ϵ=1.3have. Figures 11(a), (b), (c) show the effect of motion by mean curvature for the AC equation with the different spatial steps. Figure 12 shows temporal evolutions of analytic and numerical radii for different spatial steps. The time steps are used as Δt=0.2h2ave for each have. We can observe that the numerical results agree with the analytic exact solution as the spatial step decreases to zero.
In this section, we observe the dynamics of the AC equation on various domains. We consider rectangular domain with circular hole Ωrectangle=((−1.3,1.3)×(−1.3,1.3))∩B2(0.5)c and circular domain with circular hole Ωcircle=B1(1.3)∩B2(0.5)c to verify the performance of the proposed method according to the shape of the domain, where B1(r)={(x,y)|x2+y2<r} and B2(r)={(x,y)|x2+y2≤r}. The initial condition is that all internal node points in the domain are 1 with Dirichlet boundary condition u(x,y,t)=1 in ∂Ωinner and u(x,y,t)=−1 in ∂Ωouter. Here, ∂Ωinner and ∂Ωouter are the inner and outer boundaries of the domain, respectively. We use the parameters have≈0.2,0.15,0.1, ϵ=0.15, Δt=0.001, and Nt=T/Δt to compute Eq (2.2). Here, we use T=5.6, which is a sufficient time to reach an numerical equilibrium state for all parameters. The value of have is below each figure. Figures 13 and 14 show the dynamics of the AC equation on Ωrectangle and Ωcircle, respectively, for the given initial conditions.
In this study, we proposed a novel, effective, simple, and explicit numerical method for solving the diffusion and AC equations on a triangular mesh. In the proposed algorithm, we defined a novel discrete Laplace operator. The proposed discrete Laplace operator has a property that the adjusted net vector of each point using virtual nodes is zero vector on the triangular mesh. This means that the forces acting at each point of the triangular mesh are symmetrically distributed, thus demonstrating the symmetry of the proposed method. To investigate the efficiency and stability of the proposed numerical method, we presented the numerical experiments for solving the diffusion and AC equations in various domains on triangular mesh. We verified the proposed Laplace operator through numerical experiments on general discrete Laplace operators and discrete LBOs. Additionally, we verified the proposed operator division method through numerical experiments involving motion by mean curvature using the AC equation and solving the AC equation in irregular domains. Results indicated that the proposed method can be served in various applications in 2D domain with triangular mesh. In future work, we will solve the AC equation on complex surfaces in 3D space and conduct theoretical analysis on the convergence of the proposed method.
The authors have not used Artificial Intelligence (AI) tools in the creation of this article.
The corresponding author (J. S. Kim) expresses thanks for the support from the BK21 FOUR program. The authors express their sincere gratitude to the reviewers for providing valuable feedback on this revised version. Their constructive comments have significantly contributed to the improvement of the manuscript.
The authors declare there is no conflicts of interest.
[1] |
J. W. Choi, H. G. Lee, D. Jeong, J. Kim, An unconditionally gradient stable numerical method for solving the Allen–Cahn equation, Phys. A, 388 (2009), 1791–1803. https://doi.org/10.1016/j.physa.2009.01.026 doi: 10.1016/j.physa.2009.01.026
![]() |
[2] |
S. M. Allen, J. W. Cahn, A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening, Acta Metall., 27 (1979), 1085–1095. https://doi.org/10.1016/0001-6160(79)90196-2 doi: 10.1016/0001-6160(79)90196-2
![]() |
[3] |
D. Jeong, Y. Li, Y. Choi, C. Lee, J. Yang, J. Kim, A practical adaptive grid method for the Allen–Cahn equation, Phys. A, 573 (2021), 125975. https://doi.org/10.1016/j.physa.2021.125975 doi: 10.1016/j.physa.2021.125975
![]() |
[4] |
X. Xiao, R. He, X. Feng, Unconditionally maximum principle preserving finite element schemes for the surface Allen–Cahn type equations, Numer. Meth. Part Differ. Equations, 36 (2020), 418–438. https://doi.org/10.1002/num.22435 doi: 10.1002/num.22435
![]() |
[5] |
J. Rubinstein, P. Sternberg, Nonlocal reaction—diffusion equations and nucleation, IMA J. Appl. Math., 48 (1992), 249–264. https://doi.org/10.1093/imamat/48.3.249 doi: 10.1093/imamat/48.3.249
![]() |
[6] |
Z. Sun, S. Zhang, A radial basis function approximation method for conservative Allen–Cahn equations on surfaces, Appl. Math. Lett., 143 (2023), 108634. https://doi.org/10.1016/j.aml.2023.108634 doi: 10.1016/j.aml.2023.108634
![]() |
[7] |
B. Xia, Y. Li, Z. Li, Second-order unconditionally stable direct methods for Allen–Cahn and conservative Allen–Cahn equations on Surfaces, Mathematics, 8 (2020), 1486. https://doi.org/10.3390/math8091486 doi: 10.3390/math8091486
![]() |
[8] |
Y. Choi, Y. Li, C. Lee, H. Kim, J. Kim, Explicit hybrid numerical method for the Allen–Cahn type equations on curved surfaces, Numer. Math. Theory Methods Appl., 14 (2021), 797–810. https://doi.org/10.4208/nmtma.OA-2020-0155 doi: 10.4208/nmtma.OA-2020-0155
![]() |
[9] |
V. Joshi, R. K. Jaiman, An adaptive variational procedure for the conservative and positivity preserving Allen–Cahn phase-field model, J. Comput. Phys., 366 (2018), 478–504. https://doi.org/10.1016/j.jcp.2018.04.022 doi: 10.1016/j.jcp.2018.04.022
![]() |
[10] |
S. Kwak, J. Yang, J. Kim, A conservatice Allen–Cahn equation with a curvature-dependent Lagrange multiplier, Appl. Math. Lett., 126 (2022), 107838. https://doi.org/10.1016/j.aml.2021.107838 doi: 10.1016/j.aml.2021.107838
![]() |
[11] |
Q. Hong, Y. Gong, J. Zhao, Q. Wang, Arbitrarily high order structure-preserving algorithms for the Allen–Cahn model with a nonlocal constraint, Appl. Numer. Math., 170 (2021), 321–339. https://doi.org/10.1016/j.apnum.2021.08.002 doi: 10.1016/j.apnum.2021.08.002
![]() |
[12] |
Z. Chai, D. Sun, H. Wang, B. Shi, A comparative study of local and nonlocal Allen–Cahn equations with mass conservation, Int. J. Heat Mass Transf., 122 (2018), 631–642. https://doi.org/10.1016/j.ijheatmasstransfer.2018.02.013 doi: 10.1016/j.ijheatmasstransfer.2018.02.013
![]() |
[13] |
M. Inc, A. Yusuf, A. I. Aliyu, D. Baleanu, Time-fractional Cahn–Allen and time-fractional Klein–Gordon equations: Lie symmetry analysis, explicit solutions and convergence analysis, Phys. A, 493 (2018), 94–106. https://doi.org/10.1016/j.physa.2017.10.010 doi: 10.1016/j.physa.2017.10.010
![]() |
[14] |
J. Kim, H. Lee, Unconditionally energy stable second-order numerical scheme for the Allen–Cahn equation with a high-order polynomial free energy, Adv. Differ. Equations, 2021 (2021), 1–13. https://doi.org/10.1186/s13662-021-03571-x doi: 10.1186/s13662-021-03571-x
![]() |
[15] |
H. Bhatt, J. Joshi, I. Argyros, Fourier spectral high-order time-stepping method for numerical simulation of the multi-dimensional Allen–Cahn equations, Symmetry, 13 (2021), 245. https://doi.org/10.3390/sym13020245 doi: 10.3390/sym13020245
![]() |
[16] |
Z. Sun, Y. Gao, High order multiquadric trigonometric quasi-interpolation method for solving time-dependent partial differential equations, Numer. Algorithms, (2022), 1–21. https://doi.org/10.1007/s11075-022-01486-6 doi: 10.1007/s11075-022-01486-6
![]() |
[17] |
Z. Sun, L. Ling, A kernel-based meshless conservative Galerkin method for solving Hamiltonian wave equations, SIAM J. Sci. Comput., 44 (2022), A2789–A2807. https://doi.org/10.1137/21M1436981 doi: 10.1137/21M1436981
![]() |
[18] | M. Wardetzky, S. Mathur, F. Kälberer, E. Grinspun, Discrete Laplace operators: No free lunch, Eurographics Symp. Geom. Process., (2007), 33–37. |
[19] |
G. Xu, Discrete Laplace–Beltrami operators and their convergence, Comput. Aided Geom. Des., 21 (2004), 767–784. https://doi.org/10.1016/j.cagd.2004.07.007 doi: 10.1016/j.cagd.2004.07.007
![]() |
[20] |
T. Caissard, D. Coeurjolly, J. O. Lachaud, T. Roussillon, Laplace-beltrami operator on digital surfaces, J. Math. Imaging Vis., 61 (2019), 359–379. https://doi.org/10.1007/s10851-018-0839-4 doi: 10.1007/s10851-018-0839-4
![]() |
[21] |
S. P. Thampi, S. Ansumali, R. Adhikari, S. Succi, Isotropic discrete Laplacian operators from lattice hydrodynamics, J. Comput. Phys., 234 (2013), 1–7. https://doi.org/10.1016/j.jcp.2012.07.037 doi: 10.1016/j.jcp.2012.07.037
![]() |
[22] | B. J. McCartin, Eigenstructure of the discrete Laplacian on the equilateral triangle: the Dirichlet & Neumann problems, Appl. Math. Sci., 4 (2010), 2633–2646. |
[23] | V. Ganzha, R. Liska, M. Shashkov, C. Zenger, Support operator method for Laplace equation on unstructured triangular grid, Selcuk J. Appl. Math., 3 (2002), 21–48. |
[24] |
S. Yoon, J. Park, J. Wang, C. Lee, J. Kim, Numerical simulation of dendritic pattern formation in an isotropic crystal growth model on curved surfaces, Symmetry, 12 (2020), 1155. https://doi.org/10.3390/sym12071155 doi: 10.3390/sym12071155
![]() |
[25] |
L. Paquet, R. Korikache, The complete discretization of the dual mixed method for the heat diffusion equation in a polygonal domain, Math. Comput. Simul., 186 (2021), 145–160. https://doi.org/10.1016/j.matcom.2020.09.023 doi: 10.1016/j.matcom.2020.09.023
![]() |
[26] |
J. W. Cahn, J. E. Hilliard, Free energy of a nonuniform system. Ⅰ. Interfacial free energy, J. Chem. Phys., 28 (1958), 258–267. https://doi.org/10.1063/1.1744102 doi: 10.1063/1.1744102
![]() |
[27] |
W. Tian, Y. Chen, Z. Meng, H. Jia, An adaptive finite element method based on superconvergent cluster recovery for the Cahn–Hilliard equation, Electron. Res. Arch., 31 (2023), 1323–1343. https://doi.org/10.3934/era.2023068 doi: 10.3934/era.2023068
![]() |
[28] |
Y. Hwang, C. Lee, S. Kwak, Y. Choi, S. Ham, S. Kang, et al., Benchmark problems for the numerical schemes of the phase-field equations, Discrete Dyn. Nat. Soc., 2022 (2022), 1–10. https://doi.org/10.1155/2022/2751592 doi: 10.1155/2022/2751592
![]() |
[29] |
D. Jeong, J. Kim, An explicit hybrid finite difference scheme for the Allen–Cahn equation, J. Comput. Appl. Math., 340 (2018), 247–255. https://doi.org/10.1016/j.cam.2018.02.026 doi: 10.1016/j.cam.2018.02.026
![]() |
1. | Junxiang Yang, Jian Wang, Soobin Kwak, Seokjun Ham, Junseok Kim, A modified Allen–Cahn equation with a mesh size-dependent interfacial parameter on a triangular mesh, 2024, 304, 00104655, 109301, 10.1016/j.cpc.2024.109301 | |
2. | Chaeyoung Lee, Yunjae Nam, Minjoon Bang, Seokjun Ham, Junseok Kim, Numerical investigation of the dynamics for a normalized time-fractional diffusion equation, 2024, 9, 2473-6988, 26671, 10.3934/math.20241297 | |
3. | Junseok Kim, Modified Wave-Front Propagation and Dynamics Coming from Higher-Order Double-Well Potentials in the Allen–Cahn Equations, 2024, 12, 2227-7390, 3796, 10.3390/math12233796 | |
4. | Soobin Kwak, Yongho Choi, Jian Wang, Yunjae Nam, Junseok Kim, Phase-field modeling for curvature-dependent tissue growth on surfaces, 2025, 171, 09557997, 106090, 10.1016/j.enganabound.2024.106090 | |
5. | Seokjun Ham, Jaeyong Choi, Yunjae Nam, Junseok Kim, Stability analysis of a numerical method for the 3D high-order Allen–Cahn equation, 2025, 15, 2158-3226, 10.1063/5.0248165 | |
6. | Yongho Choi, Youngjin Hwang, Soobin Kwak, Seokjun Ham, Hyundong Kim, Junseok Kim, A cell structure implementation of the multigrid method for the two-dimensional diffusion equation, 2025, 15, 2158-3226, 10.1063/5.0247042 | |
7. | Soobin Kwak, Seokjun Ham, Jian Wang, Hyundong Kim, Junseok Kim, Effective perpendicular boundary conditions in phase-field models using Dirichlet boundary conditions, 2025, 0177-0667, 10.1007/s00366-025-02109-z | |
8. | Junxiang Yang, Seungyoon Kang, Sangkwon Kim, Youngjin Hwang, Soobin Kwak, Seokjun Ham, Junseok Kim, An efficient computational method for simulating incompressible fluid flows on a virtual cubic surface, 2025, 144, 10075704, 108676, 10.1016/j.cnsns.2025.108676 | |
9. | Mengyu Luo, Chaeyoung Lee, Jian Wang, Soobin Kwak, Hyundong Kim, Junseok Kim, Mengxin Chen, Designing Team Projects for Envy‐Free Group Collaboration to Overcome Free‐Rider Problem, 2025, 2025, 1026-0226, 10.1155/ddns/3370833 | |
10. | Mengxin Chen, Canrong Tian, Seokjun Ham, Hyundong Kim, Junseok Kim, Impact of prey-taxis on a harvested intraguild predation predator–prey model, 2025, 0956-7925, 1, 10.1017/S0956792525000087 | |
11. | Seokjun Ham, Junxiang Yang, Youngjin Hwang, Junseok Kim, A novel multi-component Allen–Cahn system for reducing the vacuum phenomenon at the triple junction, 2025, 15, 2158-3226, 10.1063/5.0261749 |
spatial step | h(=1/8) | rate | h/2 | rate | h/4 | rate | h/8 |
Δ1u | 2.5237×10−1 | 1.99 | 6.3336×10−2 | 2.00 | 1.5849×10−2 | 2.00 | 3.9633×10−3 |
Δ2u | 9.9404×10−1 | 1.98 | 2.5237×10−1 | 1.99 | 6.3336×10−2 | 2.00 | 1.5849×10−2 |
Δ3u | 4.9959×10−1 | 1.98 | 1.2635×10−1 | 2.00 | 3.1678×10−2 | 2.00 | 7.9253×10−3 |
˜Δu | 7.4682×10−1 | 1.98 | 1.8936×10−1 | 1.99 | 4.7507×10−2 | 2.00 | 1.1887×10−2 |
spatial step | h(=1/8) | rate | h/2 | rate | h/4 | rate | h/8 |
Δ1u | 1.1216×10−1 | 1.91 | 2.9805×10−2 | 1.96 | 7.6845×10−3 | 1.98 | 1.9511×10−3 |
Δ2u | 4.4180×10−1 | 1.90 | 1.1876×10−1 | 1.95 | 3.0708×10−2 | 1.98 | 7.8027×10−3 |
Δ3u | 2.2204×10−1 | 1.90 | 5.9457×10−2 | 1.95 | 1.5359×10−2 | 1.98 | 3.9017×10−3 |
˜Δu | 3.3192×10−1 | 1.90 | 8.9110×10−2 | 1.95 | 2.3034×10−2 | 1.98 | 5.8522×10−3 |
Spatial step | have=0.1254 | rate | have=0.1013 | rate | have=0.0762 |
Maximum-error | 1.8338×10−2 | 1.19 | 1.4246×10−2 | 1.92 | 8.2534×10−3 |
l2-error | 4.2098×10−3 | 1.64 | 2.9665×10−3 | 1.77 | 1.7904×10−3 |
Spatial step | have=0.1254 | rate | have=0.1013 | rate | have=0.0762 |
Maximum-error | 4.6855×10−3 | 2.07 | 3.0121×10−3 | 1.98 | 1.7142×10−3 |
l2-error | 2.0184×10−3 | 2.89 | 1.0894×10−3 | 1.09 | 7.9989×10−4 |
Spatial step | have=0.1254 | rate | have=0.1013 | rate | have=0.0762 |
Maximum-error | 5.3441×10−3 | 2.11 | 3.4101×10−3 | –7.49 | 2.8734×10−2 |
l2-error | 2.3112×10−3 | 2.01 | 1.5061×10−3 | –7.16 | 1.1569×10−2 |
spatial step | h(=1/8) | rate | h/2 | rate | h/4 | rate | h/8 |
Δ1u | 2.5237×10−1 | 1.99 | 6.3336×10−2 | 2.00 | 1.5849×10−2 | 2.00 | 3.9633×10−3 |
Δ2u | 9.9404×10−1 | 1.98 | 2.5237×10−1 | 1.99 | 6.3336×10−2 | 2.00 | 1.5849×10−2 |
Δ3u | 4.9959×10−1 | 1.98 | 1.2635×10−1 | 2.00 | 3.1678×10−2 | 2.00 | 7.9253×10−3 |
˜Δu | 7.4682×10−1 | 1.98 | 1.8936×10−1 | 1.99 | 4.7507×10−2 | 2.00 | 1.1887×10−2 |
spatial step | h(=1/8) | rate | h/2 | rate | h/4 | rate | h/8 |
Δ1u | 1.1216×10−1 | 1.91 | 2.9805×10−2 | 1.96 | 7.6845×10−3 | 1.98 | 1.9511×10−3 |
Δ2u | 4.4180×10−1 | 1.90 | 1.1876×10−1 | 1.95 | 3.0708×10−2 | 1.98 | 7.8027×10−3 |
Δ3u | 2.2204×10−1 | 1.90 | 5.9457×10−2 | 1.95 | 1.5359×10−2 | 1.98 | 3.9017×10−3 |
˜Δu | 3.3192×10−1 | 1.90 | 8.9110×10−2 | 1.95 | 2.3034×10−2 | 1.98 | 5.8522×10−3 |
Spatial step | have=0.1254 | rate | have=0.1013 | rate | have=0.0762 |
Maximum-error | 1.8338×10−2 | 1.19 | 1.4246×10−2 | 1.92 | 8.2534×10−3 |
l2-error | 4.2098×10−3 | 1.64 | 2.9665×10−3 | 1.77 | 1.7904×10−3 |
Spatial step | have=0.1254 | rate | have=0.1013 | rate | have=0.0762 |
Maximum-error | 4.6855×10−3 | 2.07 | 3.0121×10−3 | 1.98 | 1.7142×10−3 |
l2-error | 2.0184×10−3 | 2.89 | 1.0894×10−3 | 1.09 | 7.9989×10−4 |
Spatial step | have=0.1254 | rate | have=0.1013 | rate | have=0.0762 |
Maximum-error | 5.3441×10−3 | 2.11 | 3.4101×10−3 | –7.49 | 2.8734×10−2 |
l2-error | 2.3112×10−3 | 2.01 | 1.5061×10−3 | –7.16 | 1.1569×10−2 |