An explicit numerical method for the conservative Allen–Cahn equation on a cubic surface

1.   Introduction

The classic Allen–Cahn (AC) equation ^[1] is a prototypical example of gradient flow associated with a functional:

where $\Omega \subset {\mathbb R}^d$, $u({\bf x}, t)$ is an order parameter at spatial point ${\bf x}$ and time $t$, $\epsilon$ represents the interface width parameter, and $F(u) = 0.25(u^2-1)^2$. Thus, the AC equation can be derived through a variational approach as $\partial u / \partial t = - \delta \mathcal{E} / \delta u$, which can be expressed as follows.

The AC equation is a fundamental phase-field model in a wide range of applications such as materials science ^[2,3,4], physics ^[5], classification ^[6,7], and image processing ^[8,9]. Its versatility lies in its ability to describe complex phenomena such as the evolution of phase transitions, where distinct phases within a material change over time, and the dynamics of interfaces that separate these phases. These characteristics make the AC equation an essential tool in understanding interface motion and phase separation in various scientific and engineering contexts ^{[10,11,12,13]}. The equation captures critical aspects of interface dynamics, particularly the behavior of systems as they evolve toward a minimal energy state, driven by forces such as curvature. This makes it a key model for studying processes where the morphology of interfaces plays a central role. The conventional AC equation is not conservative but preserves the maximum principle. Recently, several studies have focused on the maximum-principle-preserving characteristic of the AC equation. Zhang et al. ^[14] developed and analyzed the maximum principle preserving method for solving the AC equation. Sun et al. ^[15] presented a class of up to temporally eighth-order maximum principle preserving method for the AC equation. Choi et al. ^[16] researched stability analysis and maximum principle property of the fully explicit finite difference scheme for the high-order AC equation. Kim et al. ^[17] presented a hybrid numerical method on nonuniform grids for solving the AC equation.

However, in the situations where mass conservation is critical, a conservative variant of the AC equation is used, which ensures that the total quantity of the evolving order parameter remains constant over time ^[18]. To make the AC equation conservative in nature, several different Lagrange multipliers have been proposed to satisfy the mass conservation law for the AC equation. For example, Rubinstein and Sternberg ^[19] introduced a nonlocal space-independent Lagrange multiplier into the AC equation. Brassel and Bretin ^[20] proposed a local-nonlocal space-time dependent Lagrange multiplier for mass conservation of the AC equation. The two conservative forms of the AC equation have received significant attention, and the conservative AC (CAC) equation satisfies mass conservation while effectively modeling the dynamics of interfaces. Furthermore, compared to the fourth-order Cahn–Hilliard equation, the CAC equation is a second-order partial differential equation (PDE), which makes it simpler to solve ^[21]. These advantages have motivated a wide range of theoretical and numerical studies on the CAC equation ^[22,23,24]. Xia et al. ^[25] proposed temporally second-order unconditionally stable direct schemes for the AC and CAC equations on surfaces. They found that the proposed method is unconditionally energy-stable. Sun and Zhang ^[26] developed a meshless radial basis function approximation method for the CAC equation on surfaces based on an operator splitting scheme. The developed numerical method satisfies the mass conservation law of the CAC equation. The authors demonstrated the accuracy and mass conservation law of the developed method through various numerical experiments on different smooth and compact surfaces. Liu et al. ^[27] presented multi-physical structure-preserving method for the CAC equation. Yang and Kim ^[28] developed a numerical method for the conservative Allen–Cahn–Navier–Stokes model on arbitrarily curved surfaces in three-dimensional (3D) space. Previously, the authors proposed a numerical solution for the Cahn–Hilliard–Navier–Stokes system on arbitrarily curved surfaces in 3D space ^[29]. However, simulating mass-conserving binary flow on such surfaces can be solved more efficiently and straightforwardly using the CAC equation instead of the Cahn–Hilliard equation ^[30]. This is because the Cahn–Hilliard equation is a fourth-order PDE, whereas the conservative Allen–Cahn equation is second-order, which results in reduced computational cost and complexity. Figure 1 shows the solution obtained using the numerical method developed by Yang and Kim ^[28] for the conservative Allen–Cahn–Navier–Stokes model. The our proposed method focuses on efficiently and rapidly computing the CAC equation in two-dimensionality (2D) on the surface of a cuboid that approximates a spherical surface in 3D.

The analysis of the numerical methods based on the finite difference method for the CAC equation is highly challenging. Consequently, various numerical methods have either replaced the analysis with numerical verification of properties such as the maximum principle, energy stability, and convergence or modified the CAC equation to facilitate analysis. Choi and Kim ^[31] presented a novel CAC equation with a new Lagrange multiplier that ensures strict preservation of the maximum principle and an unconditionally stable numerical method, which was validated through computational experiments. Moreover, solving this equation on complex geometries, such as curved surfaces, presents challenges due to the need to accurately capture surface geometry while ensuring numerical stability and conservation properties. Kim et al. ^[32] developed the unconditionally stable hybrid scheme for the CAC equation with a space-time dependent Lagrange multiplier. The authors validated the numerical methods through various numerical experiments, such as practical unconditional stability tests.

In this paper, we focus on a simple and effective finite difference method (FDM) for solving the CAC equation on cubic surfaces. In other words, the key aspect is to consider the unique characteristics of cubic surfaces. The proposed method is simple to implement and can serve as a foundation for various methods such as implicit and scalar auxiliary variable approaches. The cubic surfaces, which are characterized by their sharp edges and non-smooth features, pose additional difficulties for numerical approximation methods ^[33]. Traditional schemes may struggle with maintaining mass conservation and interface accuracy in such settings, making the development of specialized algorithms critical. Therefore, the development of finite difference schemes designed for the CAC equation on surfaces such as cubic geometries is an area of growing interest.

The structure of this paper is as follows. Section 2 presents the governing CAC equation. Section 3 presents the proposed computational scheme. Section 4 presents numerical tests. In Section 6, conclusions are presented.

2.   Governing equation

We consider the following CAC equation ^[20]:

where $\Omega$ is a bounded domain in $\mathbb{R}^d\; (d = 1, 2, 3)$, $u({\bf x}, t)$ is conserved order parameter, and $\beta(t)$ is given by

We define the unit cubic surface domain $\hat \Omega$ in 3D space as follows to solve the problem on the unit cubic surface.

Then, the cubic surface domain $\hat \Omega$ is decomposed into sub-domains $\hat \Omega_k$, $k = 1, 2, \ldots, 6$.

where $\overline {\hat \Omega}_k$ denotes the closure of the set ${\hat \Omega}_k$ for $k = 1, \ldots, 6$. To define the unfolded cubic surface domain, we define the sub-domains in 2D space as follows.

Thus, the unfolded cubic surface domain is defined by $\Omega = \cup_{k = 1}^6{\overline \Omega_k}$. Figure 2 displays the schematic illustration of the folded and unfolded cubic surface domains. We can observe that the proposed numerical scheme can be easily extended to surfaces of arbitrarily sized cuboids. Unless stated otherwise, we consider the unit cubic surface.

We consider the folded cubic surface in 3D space and apply proper boundary conditions to each sub-domain $\Omega_k$ for $k = 1, 2, \ldots, 6$. For example, we focus a sub-domain $\Omega_4$, and the boundary conditions are shown Figure 3(a).

Similarly, we apply the boundary condition for the other sub-domain $\Omega_k$, $k = 1, 2, 3, 5, 6$ considering the folded cubic surface domain. To prove the mass conservation law of the CAC equation on a cubic surface, we define the total mass, denoted by $\mathcal{M}(u)$, as

which can be written as

Therefore, it is equivalent to state that $\int_{\Omega} u({\bf x}, t) d{\bf x}$ is constant, which implies that $\mathcal{M}(u)$ is constant. Taking the time derivative of Eq (2.3), we have

Considering that the unfolded domain $\Omega$ consists of a combination of six sub-domains $\Omega_k$, $k = 1, 2, \ldots, 6$, we obtain the following equation:

By using the equation above and Green's theorem, we can rewrite Eq (2.5) as follows.

where $\partial \Omega_k$, $k = 1, 2, \ldots, 6$ are curves oriented in the positive direction. We calculate Eq (2.6) by breaking up $\partial \Omega_k$ for $k = 1, 2, \ldots, 6$ as the union of the four curves $\partial \Omega_{k, 1}$, $\partial \Omega_{k, 2}$, $\partial \Omega_{k, 3}$, and $\partial \Omega_{k, 4}$, as shown for $\Omega_4$ in Figure 3(b).

For each $\partial \Omega_{k, c}$, $k = 1, \ldots, 6$, $c = 1, 2, \ldots, 4$, there exist line integrals with opposite directions that cancel each other out. Therefore, we have

which shows that the mass conservation law of the CAC equation on the unfolded cubic surface is satisfied.

3.   Numerical method

Now, we describe the numerical method of the CAC equation on the cubic surface. Let the number of grid $N$ for a one sub-domain be a positive integer. Then, we define the set of indexes for the sub-domain $\Omega_k$, $k = 1, 2, \ldots, 6$ and $\Omega$.

The discretization domains for the sub-domains $\Omega_k$, $k = 1, \cdots, 6$ in 2D space are defined by $\Omega^d_k = \{(x_i = (i-0.5)h, y_j = (j-0.5)h)\; |\; (i, j) \in I_k\}$, where $h = 1/N$ is a space step size. We discretize unfolded cubic surface domain $\Omega$ as $\Omega^d = \cup_{k = 1}^6 \Omega_k^d$, see Figure 4(a).

Next, we consider the boundary conditions for the domain $\Omega^d$. For clarity, we focus only on one sub-domain $\Omega_4^d$, within the unfolded cubic surface domain $\Omega^d$. The boundary conditions at $x = 2$ and $x = 3$ for $\Omega_4^d$ are appropriately defined by the sub-domains $\Omega_3^d$ and $\Omega_5^d$ of the discrete unfolded cubic surface domain $\Omega^d$. However, the boundary conditions at $y = 1$ and $y = 2$ for $\Omega_4^d$ are defined by rotating $\Omega_1^d$ and $\Omega_6^d$ counterclockwise and clockwise by $180$ degrees, respectively, and then placing them below and above $\Omega_4^d$ respectively. These boundary conditions are schematically shown in Figure 4(b). Thus, considering the folded cubic surface, we apply boundary conditions for each subdomain $\Omega_k^d$, $k = 1, 2, \ldots, 6$ of the discrete unfolded cubic surface domain. For $m = 1, \dots, N$,

The discrete total mass is defined by

The proposed algorithm consists of two steps. First, we compute

Thus, we have

Next, we compute the conservation term using the solution of Eq (3.1) to solve the CAC equation.

where

Theorem 1. The proposed method satisfies the discrete mass conservation law of the CAC equation.

Proof. For any nonnegative integer $n$, we can rewrite Eq (3.2) by the definition of $\beta^*$ as

Thus, we obtain the following equation:

Hence, the proposed scheme satisfies the mass conservation law of the CAC equation by the above equation.

□

4.   Numerical experiments

Now, we perform the numerical simulation for the CAC equation using the proposed method. We define the three-dimensional domain ${\hat \Omega}$ as the cubic surface with a side length $L = 5$ and then $\Omega^d$ is defined as the corresponding discrete unfolded cubic surface domain. The normalized discrete total mass is defined by $\mathcal{M}^h(u^n)/L^2$. The interface layer parameter with the spatial step size is defined as

which means there are $m$ grid points in the interface layer. First, we consider the following initial condition on $\Omega^d$.

where rand$(x_i, y_j)$ is random value from $-1$ to $1$ at the point $(x_i, y_j)$. The parameters used are $N = 128$, $\epsilon = \epsilon_8$, $\Delta t = 0.5\epsilon^2 h^2/(2\epsilon^2+h^2)$, and $T = 100.788$. We define the maximum and minimum values of the computational solutions at $n\Delta t$ as

respectively. Figure 5 displays the temporal evolution of the computational solution for the CAC equation. The computational results from the numerical simulation show a distinct process of domain coarsening, where larger domains grow progressively larger while smaller domains shrink and eventually disappear. This phenomenon reflects the tendency of the system to minimize its interfacial energy, which leads to the dominance of larger regions and the elimination of smaller ones over time.

We observed that the computational solution calculated using the proposed algorithm satisfies the discrete mass conservation. Furthermore, we maintain the properties of the CAC equation when the interface of numerical solutions of the CAC equation solved by the proposed scheme crosses the boundary between different $\Omega_k^d$.

Next, the following initial condition is considered to observe the effect of the unit folded cubic surface domain.

The parameters used are $N = 128$, $\epsilon = \epsilon_8$, $\Delta t = 0.5\epsilon^2 h^2/(2\epsilon^2+h^2)$, and $T = 7.727$. Figure 6 shows the computational solutions for the CAC equation at time $t = 0$, 1.176, 4.368, and 7.727.

We can observe that the interface layer of computational solutions at $\Omega_k$ for $k = 2, 3, 4, 5$ tends to approach the boundary with $\Omega_6$ by the CAC equation, while the characteristics of folded cubic surfaces result in the interface layer crossing the boundary and forming a circle according to the geometric properties of the CAC equation. We use the following initial condition, similar to the above, except with some different phase positions.

The parameters used are $N = 128$, $\epsilon = \epsilon_8$, $\Delta t = 0.5\epsilon^2 h^2/(2\epsilon^2+h^2)$, and $T = 1.680$. Figure 7 shows the numerical solutions for the CAC equation on the folded cubic surface using the proposed scheme.

We observe that the computational solution quickly becomes equilibrium in $\Omega_k$ for $k = 2, 3, 4, 5$ since the interface layer does not cross the boundary of $\Omega_6$, unlike the results in Figure 6.

Next, the two initial conditions on $\Omega^d$ are given by

where $w = 0.4$ or $w = 0.5$. We used parameters $N = 128$, $\epsilon = \epsilon_8$, $\Delta t = 0.5\epsilon^2h^2/(2\epsilon^2+h^2)$, $T = 1.680$. The top and bottom rows in Figure 8 show the temporal evolution of the numerical solutions of the CAC equation with $w = 0.4$ and $w = 0.5$, respectively, at times $t = 0$, $0.504$, $0.840$, and $1.680$.

We investigate the property of motion by mass conserving mean curvature flow on a cubic surface. The initial condition is defined as follows:

where $x \bmod 5$ represents the remainder when $x$ is divided by $5$. The parameters used are $N = 100$, $\epsilon = \epsilon_8$, $\Delta t = 0.5\epsilon^2h^2/(2\epsilon^2+h^2)$, and $T = 1.927$. Figure 9 shows the temporal evolution of $u$. In Fig. 9(a), the phase of $u$ presents like a cross-shaped structure with sharp interfaces. As time progresses, the interface gradually smooths out due to mass conserving mean curvature flow. As shown in Figure 9(d), the interface eventually evolves into a stable, circular shape near the top of the domain.

5.   Discussion

In this section, we discuss the application of an implicit scheme for the CAC equation on cubic surfaces, along with its potential benefits and limitations of the implicit scheme, and compare the implicit scheme with the proposed explicit numerical method through numerical experiments. Generally, the implicit scheme provides higher stability than the explicit scheme. This means that, compared to the time step size limitation of the explicit scheme, the implicit scheme can use relatively larger time step sizes without causing the numerical solution to blow up, thus maintaining stability ^[34,35]. However, since the implicit scheme requires solving a nonlinear system at each time step, the computations can become more complex and time-consuming depending on the given discretization grid ^[36]. To compare the proposed method with the implicit scheme, a nonlinear convex splitting method is applied to the CAC equation ^[37].

Let $u_{ij}^{*, s}$ and $u_{ij}^{*, s+1}$ be the approximations of $u_{ij}^*$ before and after a Gauss–Seidel iteration, respectively. We linearize the nonlinear term $(u_{ij}^*)^3$ in the Gauss–Seidel iteration method as

Thus, Eq. (5.1) can be rewritten in Gauss–Seidel form as follows.

where

The discrete $l_2$-norm is defined as $\Vert u^n \Vert_2 = \sqrt{\sum_{(i, j)\in I^d}(u_{ij}^n)^2/(6N^2)}$. We calculate Eq. (5.3) repeatedly until the $l_2$-norm of the consecutive error $\Vert u_{ij}^{*, s+1}-u_{ij}^{*, s}\Vert_2$ is less than a given tolerance $tol = 10^{-6}$. For numerical simulation, the initial condition on $\Omega^d$ is given by

The parameters used are $N = 128$, $\epsilon = \epsilon_8$, and the final time $T = 12.5$. The implicit convex splitting scheme allows for the use of larger time step sizes. Therefore, we use $\Delta t = 0.00025$ for the proposed method, while the implicit method uses a time step size that is $8$ times larger, $\Delta t = 0.002$. Figure 10 shows the zero-contour of the computational solutions of the CAC equation using the proposed and implicit convex splitting methods. The solid line shows the computational solution from the proposed scheme, while the dashed line displays that from the implicit scheme.

We observed that the numerical solution using the implicit scheme is less affected by the motion by mass conserving mean curvature compared to the numerical solution using the proposed method. Table 1 lists the central processing unit (CPU) time for obtaining a numerical solution for time $t = 12.5$ of the CAC equation using the explicit and implicit convex splitting methods. We observed that although the time step size used in the implicit method was $8$ times larger than that used in the proposed method, the CPU time for the implicit method was more than $2$ times as long.

6.   Conclusions

In this work, we introduced a fully numerical method for solving the CAC equation on a cubic surface. The method is structured in two steps: first, solving the AC equation using an efficient explicit finite difference approach, followed by incorporating a conservation term to ensure consistency with the conservation principles of the CAC equation. Numerical experiments confirm that the proposed scheme effectively conserves discrete mass, a critical requirement for the accurate solution of the CAC equation. Furthermore, the solution showed constrained motion by mass conserving mean curvature, another essential property of the CAC equation, which validates the robustness and accuracy of the proposed method. These findings demonstrate that the scheme is both reliable and computationally efficient for solving the CAC equation on complex geometries such as cubic surfaces. In future work, we will consider unconditionally stable and high-order accurate numerical methods for the CAC equation on cubic surfaces and perform an analysis of the stability, convergence, and consistency, etc., of the numerical methods.

Author contributions

Youngjin Hwang: Conceptualization, Methodology, Software, Visualization, Formal analysis, Validation, Investigation, Writing-original draft; Jyoti: Validation, Investigation, Funding acquisition, Writing-original draft; Soobin Kwak: Visualization, Validation, Investigation, Writing-original draft; Hyundong Kim: Validation, Investigation, Funding acquisition, Writing-original draft; Junseok Kim: Conceptualization, Methodology, Project administration, Supervision, Funding acquisition, Writing-original draft;

Acknowledgments

The corresponding author (J.S. Kim) was supported by the National Research Foundation of Korea(NRF) grant funded by the Korea government(MSIT) (No. 2022R1A2C1003844). Jyoti was supported by Brain Pool program funded by the Ministry of Science and ICT through the National Research Foundation of Korea (2022H1D3A2A02081237). Hyundong Kim was supported by Basic Science Research Program through the National Research Foundation of Korea(NRF) funded by the Ministry of Education (2021R1A6A1A03044326). The authors extend their thanks to the reviewers for the valuable and constructive input they provided during the revision of the article.

Conflict of interest

Professor Hyundong Kim is the Guest Editor for AIMS Mathematics and was not involved in the editorial review or the decision to publish this article. All authors declare that there are no competing interests.

Appendix

The following code is the main program with a random initial condition, which is also available from the corresponding author's webpage:

https://mathematicians.korea.ac.kr/cfdkim/open-source-codes/