# An Adaptive Time-Stepping Algorithm for the Allen–Cahn Equation

Authors: Chaeyoung Lee; Jintae Park; Soobin Kwak; Sangkwon Kim; Yongho Choi; Seokjun Ham; Junseok Kim
Journal: Journal of Function Spaces
 (2022)

---

## Abstract

In this paper, we present a simple and accurate adaptive time-stepping algorithm for the Allen–Cahn (AC) equation. The AC equation is a nonlinear partial differential equation, which was first proposed by Allen and Cahn for antiphase boundary motion and antiphase domain coarsening. The mathematical equation is a building block for modelling many interesting interfacial phenomena such as dendritic crystal growth, multiphase fluid flows, and motion by mean curvature. The proposed adaptive time-stepping algorithm is based on the Runge–Kutta–Fehlberg method, where the local truncation error is estimated by using fourth- and fifth-order numerical schemes. Computational experiments demonstrate that the proposed time-stepping technique is efficient in multiscale computations, i.e., both the fast and slow dynamics.


---

## Body

## 1. Introduction

We present a simple and accurate adaptive time-stepping algorithm for the Allen–Cahn (AC) equation [1–3]:
(1)∂ϕx,t∂t=−F′ϕx,tϵ2+Δϕx,t,x∈Ω,t>0,where Ω⊂ℝdd=1,2,3 is a bounded domain, ϕx,t is a compositional field, Fϕ=0.25ϕ2−12, and ϵ is a positive constant. The AC equation can be derived as the L2-gradient flow of the following Ginzburg–Landau energy functional:
(2)Eϕ=∫ΩFϕϵ2+12∇ϕ2dx.The AC equation was first proposed by Allen and Cahn, which was introduced as a phenomenological model for antiphase domain coarsening in a binary alloy. The AC equation and its various modified equations are used to deal with a wide range of problems, such as its phase transition [1], motion by mean curvature [4, 5], image study [6, 7], two-phase fluid flows [8], and crystal growth [9, 10]. To understand the dynamics of the AC equation and apply it to model scientific phenomena, it is essential to develop effective and accurate numerical methods for the AC equation. In general, if we use a small uniform time step, then the fast dynamics can be accurately captured. However, the computational cost is high. On the other hand, if we use a large time step, then we can save the computational cost with less accurate numerical results.Therefore, an adaptive time-stepping technique is needed in numerical simulations of the AC equation because it has multiple time scales. There are several studies on time stepping [11, 12]. In [13], the author presented the high-order time-adaptive method for solving the AC equation and showed the computational superiority that traditional methods do not have in solving the equation. Guillén-González and Tierra [14] developed an adaptive time-stepping technique based on a residual of the discrete energy law at each time step. Its fastness and efficiency were confirmed through various numerical experiments. Fu and Yang [15] presented an example of using adaptive time stepping to show unconditional energy stability. Fast and efficient computing results are slightly different from the result using one small time step. Karasözen et al. [16] studied the numerical solutions of the AC equation with constant or degenerate mobility and polynomial or logarithmic free energy function. To reduce computational time, an adaptive time step size method was applied, which is based on a local error estimator. They used the backward Euler method and the second-order average vector field (AVF) method to compute the local error estimator. Liao et al. [17] presented a second-order backward differentiation formula (BDF2) with a variable time step for the AC equation. In [18], the author obtained a solution for the time-dependent AC equation on surfaces using an explicit time splitting scheme. An adaptive finite element method (FEM) for the AC equation was developed in [19], which is based on a second-order accurate unconditionally stable FEM and a superconvergent cluster recovery- (SCR-) based error estimation. To show the effectiveness and robustness of their proposed method, various numerical tests were carried out. In [20, 21], the implicit integration factor (IIF) type methods were employed to numerically solve the space-fractional reaction-diffusion equations including the AC equation. Zhu [22] developed stable and effective exponential Runge–Kutta methods for parabolic equations such as the AC equation, which is easy to adopt the adaptive time step technique.The adaptive time-stepping strategy has also been applied to the modified AC equations. Fast explicit operator splitting spectral method was presented with a new adaptive time-stepping algorithm to solve the fractional nonlocal AC equation [23]. The fractional-in-space phase-field models were addressed including the AC equation in [24], applying the implicit-explicit (IMEX) methods as temporal discretization. With the adaptive time-stepping algorithm based on the difference between the first-order and second-order IMEX methods, they reduced the computational costs. Recently, the adaptive time-stepping strategy was also applied to the second-order maximum principle preserving scheme and was proposed to solve the time-fractional AC equations in [25]. In particular, the strategy is chosen to capture the early rapid changes in the performance of long-term simulations. In [26], thanks to these advantages, the adaptive time-stepping strategy was adopted considering two fast L1 time-stepping methods for the time-fractional AC equation. Moreover, the authors [27] proposed a variable time step BDF2 scheme for the fractional AC equation, which is energy stable and preserves maximum bound. Therefore, it is proved that the adaptive time-stepping algorithm is useful when we consider the local or nonlocal AC equations. In addition, there are many existing works on the adaptive time stepping for a variety of local and nonlocal phase-field models: the Cahn–Hilliard equation [28], molecular beam epitaxy (MBE) model [29, 30], phase-field crystal model [31], and so on.The main purpose of the paper is to propose a simple and accurate adaptive time-stepping algorithm for the AC equation. Our proposed algorithm is based on the Runge–Kutta–Fehlberg (RKF) method [32], where the fourth- and fifth-order numerical schemes are used for the local truncation error estimation.The contents of this paper are summarized as follows. In Section2, the proposed numerical solution algorithm is presented. In Section 3, numerical results using the proposed method are shown. Conclusions are made in Section 4.


## 2. Numerical Solution

For simplicity of notation, we rewrite Equation (1) as
(3)ϕt=ACϕ,where ACϕ=−F′ϕ/ϵ2+Δϕ. In this section, we describe the AC equation in a three-dimensional domain Ω=Lx,Rx×Ly,Ry×Lz,Rz. Let us discretize the domain Ω by a regular Cartesian grid Ωh with spatial step h=Rx−Lx/Nx=Ry−Ly/Ny=Rz−Lz/Nz, where Nx, Ny, and Nz are positive integers. That is, Ωh=xi,yj,zkxi=Lx+i−0.5h,yj=Ly+j−0.5h,zk=Lz+k−0.5h,i=1,⋯,Nx,j=1,⋯,Ny,k=1,⋯,Nz. Let ϕijkn be the approximation of ϕxi,yj,zk,nΔt, where Δt=T/Nt is a time step, T is the final time, Nt is the total number of time steps, and n is a nonnegative integer. We define the discrete maximum norm as
(4)ϕ∞=max1≤i≤Nx1≤j≤Ny1≤k≤Nzϕijk.To solve Equation (3), we adopt an explicit RKF method. Due to the explicit scheme, the linear term Δϕx,t in the AC Equation (1) restricts the time step for stability of the numerical solution as Δt<0.5h2/d, where d=1,2,3 is spatial dimension. However, the time step restriction is not severe in the case of the second-order partial differential equations, and sufficiently small time steps should be taken for numerical solutions with higher accuracy [33, 34]. Therefore, we take Δtmax=0.49h2/d as the maximum time step size. Let AChϕijkn=−F′ϕijkn/ϵ2+Δhϕijkn be the discrete AC operator, where Δhϕijkn=ϕi−1,jkn+ϕi+1,jkn+ϕi,j−1,kn+ϕi,j+1,kn+ϕij,k+1n+ϕij,k−1n−6ϕijkn/h2. Here, the homogeneous Neumann boundary condition [35] is used. Given final time T, tolerance tol, current time t, time step Δt, and numerical solution φijkn at t, the RKF method is as follows. The computation for the numerical solution is repeated for i=1,⋯,Nx,j=1,⋯,Ny, and k=1,⋯,Nz. Let us define the coefficient equations:
(5)k1=ΔtAChϕijkn,k2=ΔtAChϕijkn+14k1,k3=ΔtAChϕijkn+332k1+932k2,k4=ΔtAChϕijkn+19322197k1−72002197k2+72962197k3,k5=ΔtAChϕijkn+439216k1−8k2+3680513k3−8454104k4,k6=ΔtAChϕijkn−827k1+2k2−35442565k3+18594104k4−1140k5.Let(6)R=1Δt1360k1−1284275k3−219775240k4+150k5+255k6be an approximation for the local truncation error of the RKF method. The detailed explanation of the coefficient calculation, i.e., k1,⋯,k6, and the local truncation error given by Equation (6) can be found in [36].IfR∞≤tol, then set t=t+Δt and
(7)ϕijkn+1=ϕijkn+25216k1+14082565k3+21974104k4−15k5.Otherwise, i.e.,R∞>tol, then we do not update time t and repeat the computation using a new time step. To bound truncation error produced by applying the nth-order method with a new time step size δΔt by tol, we choose q so that δnR∞≤tol [36], that is,
(8)δ≤tolR∞1/n.Then, we can takeδ satisfying (8) with n=4 as follows:
(9)δ=0.84tolR∞1/4≤tol2R∞1/4.Next, we set a new time step as follows:(10)Δt=0.1Δtifδ<0.1,δΔtif0.1≤δ≤4,4Δtifδ>4.Finally, we setΔt=minΔt,Δtmax, and if t+Δt>T, then set Δt=t+Δt−T.


## 3. Numerical Results

In this section, numerical experiments in two- and three-dimensional (2D and 3D) spaces are performed to confirm the basic properties of the AC equation and to investigate the effect of the adaptive time-stepping technique.


### 3.1. Shrinking Circle and Sphere

Asϵ⟶0, the zero level set of ϕ moves according to motion by mean curvature [5]. Given a circle in 2D space and a sphere in 3D space, we observe that their radii decrease over time. First, we consider a circle with the initial radius R0. The radius Rt=R02−2t is the analytic solution at time t [33]. The initial condition is given as
(11)ϕx,y,0=tanhR0−x2+y22ϵ,on the computational domain Ω=−1,1×−1,1 with Nx×Ny mesh. The space step size h=2/Nx, tol=1.0e-5, and an initial radius R0=0.7 are used. We consider the effects of grid sizes and ϵ on the dynamics of the AC equation. Figure 1(a) shows the zero level contours with Nx=Ny=100 and ϵ=0.05. Here, the arrow indicates the direction of temporal evolution. Figure 1(b) shows the temporal evolution of the radius Rt with different grid sizes up to time t=0.245h2. Here, ϵ=0.05 is fixed. As the number of grid points increases, we can observe the convergence of the numerical solutions to the analytic solution Rt. Figure 1(c) shows the temporal evolution with different values of ϵ. Here, Nx=Ny=100 are used and the final time is 0.245. As shown in Figure 1(c), too small and too large values of ϵ decreases and increases the evolution of the interface, respectively. This is because the numerical spatial resolution is too low and the curvature is too large to capture the detailed structure compared to ϵ.Figure 1
Given a circle (8), temporal evolution of (a) the zero level contour and (b, c) its radiusRt with different grid sizes and ϵ values, respectively.
(a)(b)(c)We conduct the similar simulation with a sphere in 3D space. LetR0 be the initial radius of the sphere, then the radius Rt=R02−4t is the analytic solution at time t [33]. The initial condition is given as
(12)ϕx,y,z,0=tanhR0−x2+y2+z22ϵ,on the computational domain Ω=−1,1×−1,1×−1,1 with Nx×Ny×Nz mesh. The space step size h=2/Nx, tol=1.0e-5, and the initial radius R0=0.7 are used. Figures 2(a)–2(c) show the temporal evolution of sphere with Nx=Ny=Nz=75 and ϵ=0.05. Figure 2(d) represents the effect of grid size. Here, ϵ value is fixed to 0.05. It is observed that the larger grid size, the closer the numerical solution is to the analytic solution. Figure 2(e) illustrates the effect of ϵ value on the dynamics of the AC equation. Here, Nx=Ny=Nz=75 are used, and the numerical solution with ϵ=0.05 agrees well with the analytic solution.Figure 2
Given a sphere (9), (a)–(c) are the temporal evolution of the sphere and (d, e) are the changes in its radius according to different grid sizes andϵ values, respectively.
(a)
t=0(b)
t=0.05(c)
t=0.11(d)(e)


### 3.2. Comparison with Previous Results

In this section, we compare our numerical results with previous studies for the AC equation [37]. We consider a square shape as an initial condition,
(13)ϕx,y,0=+1,−2≤x≤2and−2≤y≤2,−1,otherwise,on 2D domain Ω=−5,5×−5,5 with ϵ=0.02 and homogeneous Neumann boundary conditions on the boundaries of domain ∂Ω. Figure 3 shows snapshots in 2D with square shape as initial condition (13) at time t=0, 0.2, 0.3, and 0.5 from left to right. The initial square shape transforms into a circular shape over time under the effect of the AC equation dynamics. Then, we compute the shrinking area of ϕ due to the motion by mean curvature and compare it with the theoretical area using the solution at a reference time tref=0.1. We can find theoretical area using analytic solution [33]. We obtain reference area Aref and reference radius Rref at t=tref as
(14)Aref=∑i=1Nx∑j=1Ny1+ϕijtref2,Rref=Arefπ.Figure 3
(a–d) Snapshots of the temporal evolution of an initially square shape withϵ=0.02 in square domain −5,5×−5,5 at time t=0, 0.2, 0.3, and 0.5, respectively.
(a)(b)(c)(d)UsingRt=Rref2−2t−tref, we obtain theoretical area as πR2t. Figure 4 compares the theoretical area with numerical one. From the computational results as shown in Figures 3 and 4, our test results are accurate and consistent with the results from the reference test [37].Figure 4
Temporal evolution of the theoretical and numerical areas of zero level set with reference timetref=0.1.


### 3.3. Traveling Wave Solutions

Let us consider traveling wave solutions of Equation (1) with an initial condition
(15)ϕx,y,0=121−tanhx22ϵ,on 2D domain Ω=−1,9×−0.01,0.01. Then, its closed-form solution on the infinite domain is
(16)ϕx,y,t=121−tanhx−st22ϵ,where s=3/2ϵ is the speed of the traveling wave [38]. Figure 5 shows the initial condition (dotted line); and the computational results with tol=1.0e-1, 1.0e-4, and 1.0e-7 with the exact solution at time t=0.2. Here, Nx=1000, Ny=2, and ϵ=0.06 are used. In this case of the traveling wave solution, the results are almost independent of the tolerance values because it has simple unique time scale.Figure 5
Initial condition (dotted line); and the numerical results withtol=1.0e-1, 1.0e-4, and 1.0e-7 with the exact solution (solid line) at time t=0.2.


### 3.4. Effect oftol

Let us consider the effect oftol on the size of the time step. The initial conditions are
(17)ϕx,y,0=1,ifx<0.8andy<0.8,−1,otherwise,on the domain Ω=−1,1×−1,1. We use Nx=Ny=Nz=100, h=2/Nx, and ϵ=0.05. Figure 6(a) shows the temporal evolution of the zero level contours of ϕx,y,t with tol=1.0e-7. Numerical tests are conducted to observe the changes in the time step Δt over time t as shown in Figure 6(b) according to the different tolerance values tol. When a relative large tolerance is used, the time step is Δt=Δtmax except the early times where there are sharp corners. However, when a small tolerance is used, the time step varies adaptively. It takes small time steps at not only the early times but also the later times with large curvature of interface. Figure 6(d) shows the temporal evolution of the discrete total energy Eht with different tolerance values. We can confirm the total energies decrease for all cases. Here, the discrete total energy is defined as
(18)Eht=∑i=1Nx∑j=1NyFϕijtϵ2+12ϕi+1,jt−ϕijth2+12ϕi,j+1t−ϕijth2h2,where ϕijt=ϕxi,yj,t and we have used the homogeneous Neumann boundary condition. In Figure 6(d), the maximum and minimum values of ϕ are given. It is observed that the numerical solutions are bounded in −1,1.Figure 6
Temporal evolution with the initial states (14): (a) zero level contour ofϕx,y,t with tol=1.0e-7 and (b–d) the changes in the time step size Δt, the discrete total energy Eht with different tol values, and the maximum and the minimum of ϕ.
(a)(b)(c)(d)Next, we consider another initial condition with the same parameters and domain as the above test:(19)ϕx,y,0=1ifx−0.3<0.5andy−0.3<0.5,orx+0.6<0.2andy+0.6<0.2,−1otherwise.Withtol=1.0e-7, the zero level contours of ϕx,y,t over time are shown in Figure 7(a). As mentioned above, according to the property of the AC equation called motion by mean curvature, the evolution is fast for large curvatures, i.e., small radii. On the other hand, for large radii, curvatures are small and their evolution is slow. In Figure 7(b), the changes in the time step Δt are shown with different tolerance values. Figures 7(c) and 7(d) show the discrete total energy Eht and the maximum and minimum values of ϕ according to different tolerance values, respectively. Therefore, the total energies decrease and the numerical solutions ϕ are bounded in −1,1. In addition, in Figure 7(d), we compare the numerical results obtained from the proposed method and an explicit Euler method with Δt=Δtmax. Both are explicit methods; however, when adopting the proposed method, we can use a slightly larger Δt than the explicit Euler method. Let us assume that the initial condition (19) is given with a single time step, not the adaptive time step. With a large time step, we might miss capturing fast dynamics, whereas using a small time step to capture the fast dynamics is inefficient in terms of overall computational cost. For this reason, therefore, our proposed method using multitime steps has an advantage.Figure 7
Temporal evolution with the initial states (16): (a) zero level contour ofϕx,y,t with tol=1.0e-7; (b–d) the changes in the time step size Δt, the discrete total energy Eht, and the maximum and the minimum of ϕ with different tol values.
(a)(b)(c)(d)


### 3.5. Shrinking Spirals

We simulate a shrinking test with a 2D spiral shape. Figure8(a) shows the initial condition which is a spiral shape expressed as a filled zero contour in 2D. The initial value of ϕ is ϕx,y,0=1 in the inside of spiral and otherwise ϕx,y,0=−1 on the computational domain Ω=−1,1×−1,1. Other parameters are Nx=Ny=100, h=2/Nx, and ϵ=0.03. In Figures 8(a)–8(d), the results show the phenomenon that the interface is evolved according to the motion by mean curvature. From left to right, the evolutionary times are t=0,0.01,0.07, and 0.16. Here, we set tol=1.0e-5. Figure 8(e) shows the temporal evolution of the adaptive time step sizes according to three different tolerance tol=1.0e-5, 1.0e-6, and 1.0e-7.Figure 8
(a–d) Snapshots of the 2D spiral shape at timet=0,0.01,0.07, and 0.16, respectively. (e) Changes of Δt over time t according to three different tol.
(a)
t=0(b)
t=0.01(c)
t=0.07(d)
t=0.16(e)Next, we make a 3D spiral shape on the computational domainΩ=−1,1×−1,1×−1,1 and simulate a similar shrinking test. The initial condition is shown in Figure 9(a). Initial value of ϕ is defined as ϕx,y,z,0=1 in the inside of the spiral and otherwise ϕx,y,z,0=−1. We use the parameters as Nx=Ny=Nz=100, h=2/Nx, and ϵ=0.03. Figures 9(b)–9(d) illustrate that the spiral shape is shrinking by motion by mean curvature. Here, tol=1.0e-5. Figure 9(e) shows the temporal evolution of the adaptive time step sizes according to three different tolerance tol=1.0e-5, 1.0e-6, and 1.0e-7.Figure 9
(a–d) Snapshots of the 3D spiral shape at timet=0,0.03,0.06, and 0.1, respectively. (e) Temporal evolution of the time step Δt according to three different tol.
(a)
t=0(b)
t=0.03(c)
t=0.06(d)
t=0.1(e)


## 3.1. Shrinking Circle and Sphere

Asϵ⟶0, the zero level set of ϕ moves according to motion by mean curvature [5]. Given a circle in 2D space and a sphere in 3D space, we observe that their radii decrease over time. First, we consider a circle with the initial radius R0. The radius Rt=R02−2t is the analytic solution at time t [33]. The initial condition is given as
(11)ϕx,y,0=tanhR0−x2+y22ϵ,on the computational domain Ω=−1,1×−1,1 with Nx×Ny mesh. The space step size h=2/Nx, tol=1.0e-5, and an initial radius R0=0.7 are used. We consider the effects of grid sizes and ϵ on the dynamics of the AC equation. Figure 1(a) shows the zero level contours with Nx=Ny=100 and ϵ=0.05. Here, the arrow indicates the direction of temporal evolution. Figure 1(b) shows the temporal evolution of the radius Rt with different grid sizes up to time t=0.245h2. Here, ϵ=0.05 is fixed. As the number of grid points increases, we can observe the convergence of the numerical solutions to the analytic solution Rt. Figure 1(c) shows the temporal evolution with different values of ϵ. Here, Nx=Ny=100 are used and the final time is 0.245. As shown in Figure 1(c), too small and too large values of ϵ decreases and increases the evolution of the interface, respectively. This is because the numerical spatial resolution is too low and the curvature is too large to capture the detailed structure compared to ϵ.Figure 1
Given a circle (8), temporal evolution of (a) the zero level contour and (b, c) its radiusRt with different grid sizes and ϵ values, respectively.
(a)(b)(c)We conduct the similar simulation with a sphere in 3D space. LetR0 be the initial radius of the sphere, then the radius Rt=R02−4t is the analytic solution at time t [33]. The initial condition is given as
(12)ϕx,y,z,0=tanhR0−x2+y2+z22ϵ,on the computational domain Ω=−1,1×−1,1×−1,1 with Nx×Ny×Nz mesh. The space step size h=2/Nx, tol=1.0e-5, and the initial radius R0=0.7 are used. Figures 2(a)–2(c) show the temporal evolution of sphere with Nx=Ny=Nz=75 and ϵ=0.05. Figure 2(d) represents the effect of grid size. Here, ϵ value is fixed to 0.05. It is observed that the larger grid size, the closer the numerical solution is to the analytic solution. Figure 2(e) illustrates the effect of ϵ value on the dynamics of the AC equation. Here, Nx=Ny=Nz=75 are used, and the numerical solution with ϵ=0.05 agrees well with the analytic solution.Figure 2
Given a sphere (9), (a)–(c) are the temporal evolution of the sphere and (d, e) are the changes in its radius according to different grid sizes andϵ values, respectively.
(a)
t=0(b)
t=0.05(c)
t=0.11(d)(e)


## 3.2. Comparison with Previous Results

In this section, we compare our numerical results with previous studies for the AC equation [37]. We consider a square shape as an initial condition,
(13)ϕx,y,0=+1,−2≤x≤2and−2≤y≤2,−1,otherwise,on 2D domain Ω=−5,5×−5,5 with ϵ=0.02 and homogeneous Neumann boundary conditions on the boundaries of domain ∂Ω. Figure 3 shows snapshots in 2D with square shape as initial condition (13) at time t=0, 0.2, 0.3, and 0.5 from left to right. The initial square shape transforms into a circular shape over time under the effect of the AC equation dynamics. Then, we compute the shrinking area of ϕ due to the motion by mean curvature and compare it with the theoretical area using the solution at a reference time tref=0.1. We can find theoretical area using analytic solution [33]. We obtain reference area Aref and reference radius Rref at t=tref as
(14)Aref=∑i=1Nx∑j=1Ny1+ϕijtref2,Rref=Arefπ.Figure 3
(a–d) Snapshots of the temporal evolution of an initially square shape withϵ=0.02 in square domain −5,5×−5,5 at time t=0, 0.2, 0.3, and 0.5, respectively.
(a)(b)(c)(d)UsingRt=Rref2−2t−tref, we obtain theoretical area as πR2t. Figure 4 compares the theoretical area with numerical one. From the computational results as shown in Figures 3 and 4, our test results are accurate and consistent with the results from the reference test [37].Figure 4
Temporal evolution of the theoretical and numerical areas of zero level set with reference timetref=0.1.


## 3.3. Traveling Wave Solutions

Let us consider traveling wave solutions of Equation (1) with an initial condition
(15)ϕx,y,0=121−tanhx22ϵ,on 2D domain Ω=−1,9×−0.01,0.01. Then, its closed-form solution on the infinite domain is
(16)ϕx,y,t=121−tanhx−st22ϵ,where s=3/2ϵ is the speed of the traveling wave [38]. Figure 5 shows the initial condition (dotted line); and the computational results with tol=1.0e-1, 1.0e-4, and 1.0e-7 with the exact solution at time t=0.2. Here, Nx=1000, Ny=2, and ϵ=0.06 are used. In this case of the traveling wave solution, the results are almost independent of the tolerance values because it has simple unique time scale.Figure 5
Initial condition (dotted line); and the numerical results withtol=1.0e-1, 1.0e-4, and 1.0e-7 with the exact solution (solid line) at time t=0.2.


## 3.4. Effect oftol

Let us consider the effect oftol on the size of the time step. The initial conditions are
(17)ϕx,y,0=1,ifx<0.8andy<0.8,−1,otherwise,on the domain Ω=−1,1×−1,1. We use Nx=Ny=Nz=100, h=2/Nx, and ϵ=0.05. Figure 6(a) shows the temporal evolution of the zero level contours of ϕx,y,t with tol=1.0e-7. Numerical tests are conducted to observe the changes in the time step Δt over time t as shown in Figure 6(b) according to the different tolerance values tol. When a relative large tolerance is used, the time step is Δt=Δtmax except the early times where there are sharp corners. However, when a small tolerance is used, the time step varies adaptively. It takes small time steps at not only the early times but also the later times with large curvature of interface. Figure 6(d) shows the temporal evolution of the discrete total energy Eht with different tolerance values. We can confirm the total energies decrease for all cases. Here, the discrete total energy is defined as
(18)Eht=∑i=1Nx∑j=1NyFϕijtϵ2+12ϕi+1,jt−ϕijth2+12ϕi,j+1t−ϕijth2h2,where ϕijt=ϕxi,yj,t and we have used the homogeneous Neumann boundary condition. In Figure 6(d), the maximum and minimum values of ϕ are given. It is observed that the numerical solutions are bounded in −1,1.Figure 6
Temporal evolution with the initial states (14): (a) zero level contour ofϕx,y,t with tol=1.0e-7 and (b–d) the changes in the time step size Δt, the discrete total energy Eht with different tol values, and the maximum and the minimum of ϕ.
(a)(b)(c)(d)Next, we consider another initial condition with the same parameters and domain as the above test:(19)ϕx,y,0=1ifx−0.3<0.5andy−0.3<0.5,orx+0.6<0.2andy+0.6<0.2,−1otherwise.Withtol=1.0e-7, the zero level contours of ϕx,y,t over time are shown in Figure 7(a). As mentioned above, according to the property of the AC equation called motion by mean curvature, the evolution is fast for large curvatures, i.e., small radii. On the other hand, for large radii, curvatures are small and their evolution is slow. In Figure 7(b), the changes in the time step Δt are shown with different tolerance values. Figures 7(c) and 7(d) show the discrete total energy Eht and the maximum and minimum values of ϕ according to different tolerance values, respectively. Therefore, the total energies decrease and the numerical solutions ϕ are bounded in −1,1. In addition, in Figure 7(d), we compare the numerical results obtained from the proposed method and an explicit Euler method with Δt=Δtmax. Both are explicit methods; however, when adopting the proposed method, we can use a slightly larger Δt than the explicit Euler method. Let us assume that the initial condition (19) is given with a single time step, not the adaptive time step. With a large time step, we might miss capturing fast dynamics, whereas using a small time step to capture the fast dynamics is inefficient in terms of overall computational cost. For this reason, therefore, our proposed method using multitime steps has an advantage.Figure 7
Temporal evolution with the initial states (16): (a) zero level contour ofϕx,y,t with tol=1.0e-7; (b–d) the changes in the time step size Δt, the discrete total energy Eht, and the maximum and the minimum of ϕ with different tol values.
(a)(b)(c)(d)


## 3.5. Shrinking Spirals

We simulate a shrinking test with a 2D spiral shape. Figure8(a) shows the initial condition which is a spiral shape expressed as a filled zero contour in 2D. The initial value of ϕ is ϕx,y,0=1 in the inside of spiral and otherwise ϕx,y,0=−1 on the computational domain Ω=−1,1×−1,1. Other parameters are Nx=Ny=100, h=2/Nx, and ϵ=0.03. In Figures 8(a)–8(d), the results show the phenomenon that the interface is evolved according to the motion by mean curvature. From left to right, the evolutionary times are t=0,0.01,0.07, and 0.16. Here, we set tol=1.0e-5. Figure 8(e) shows the temporal evolution of the adaptive time step sizes according to three different tolerance tol=1.0e-5, 1.0e-6, and 1.0e-7.Figure 8
(a–d) Snapshots of the 2D spiral shape at timet=0,0.01,0.07, and 0.16, respectively. (e) Changes of Δt over time t according to three different tol.
(a)
t=0(b)
t=0.01(c)
t=0.07(d)
t=0.16(e)Next, we make a 3D spiral shape on the computational domainΩ=−1,1×−1,1×−1,1 and simulate a similar shrinking test. The initial condition is shown in Figure 9(a). Initial value of ϕ is defined as ϕx,y,z,0=1 in the inside of the spiral and otherwise ϕx,y,z,0=−1. We use the parameters as Nx=Ny=Nz=100, h=2/Nx, and ϵ=0.03. Figures 9(b)–9(d) illustrate that the spiral shape is shrinking by motion by mean curvature. Here, tol=1.0e-5. Figure 9(e) shows the temporal evolution of the adaptive time step sizes according to three different tolerance tol=1.0e-5, 1.0e-6, and 1.0e-7.Figure 9
(a–d) Snapshots of the 3D spiral shape at timet=0,0.03,0.06, and 0.1, respectively. (e) Temporal evolution of the time step Δt according to three different tol.
(a)
t=0(b)
t=0.03(c)
t=0.06(d)
t=0.1(e)


## 4. Conclusion

In this study, we presented a simple and accurate adaptive time-stepping algorithm for the AC equation. The proposed adaptive time-stepping algorithm was based on the Runge–Kutta–Fehlberg method, where the local truncation error was estimated using fourth- and fifth-order numerical schemes. Computational experiments demonstrated that the proposed time-stepping technique was efficient in multiscale computations, i.e., both the fast and slow dynamics. In future work, we will apply finite element version such as [39].



---

*Source: 2731593-2022-07-16.xml*