    # Parallel energy-stable solver for a coupled Allen-Cahn and Cahn-Hilliard system

In this paper, we study numerical methods for solving the coupled Allen-Cahn/Cahn-Hilliard system associated with a free energy functional of logarithmic type. To tackle the challenge posed by the special free energy functional, we propose a method to approximate the discrete variational derivatives in polynomial forms, such that the corresponding finite difference scheme is unconditionally energy stable and the energy dissipation law is maintained. To further improve the performance of the algorithm, a modified adaptive time stepping strategy is adopted such that the time step size can be flexibly controlled based on the dynamical evolution of the problem. To achieve high performance on parallel computers, we introduce a domain decomposition based, parallel Newton-Krylov-Schwarz method to solve the nonlinear algebraic system constructed from the discretization at each time step. Numerical experiments show that the proposed algorithm is second-order accurate in both space and time, energy stable with large time steps, and highly scalable to over ten thousands processor cores on the Sunway TaihuLight supercomputer.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

To describe the evolution of coexistent phases in binary alloy systems that exhibit simultaneous phase separation and one or more order-disorder transitions, it is often of great interest to study the solution of an Allen–Cahn/Cahn–Hilliard (AC/CH) system, originally introduced by Cahn and Novick–Cohen . The AC/CH system can be obtained by taking the quasi-continuum limits for the free energy defined on the lattice, so that the dynamics of a binary alloy is reduced to a gradient flow as follows:

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩∂u∂t=∇⋅c(u,v)∇δEδu,∂v∂t=−c(u,v)ρδEδv, (1)

where , are functions on , and , are the variational derivatives (in the inner product) of the total free energy functional . The first equation in (1) is the Cahn–Hilliard equation , in which represents a conserved concentration field for the phase separation. The second equation in (1) is the Allen–Cahn equation , in which denotes a non-conserved order parameter for the anti-phase coarsening. In the AC/CH system, is the mobility, which is degenerate at pure phases, and the density  is a positive constant. The total free energy functional can be formulated as

 E(u,v)=∫Ωe(u,v)dx, (2)

where the local free energy is

 e(u,v)=α2u(1−u)−β2v2+θ(Φ(u+v)+Φ(u−v))+γ2(|∇u|2+|∇v|2)

and . Here positive constants , , , are corresponding to the coefficients of the gradient energy, the entropy, the nearest neighbors pairwise energetic interactions, and next-nearest neighbors pairwise energetic interactions, respectively. We consider periodic boundary conditions or the following homogeneous Neumann boundary conditions 

 n⋅∇u|∂Ω=n⋅∇v|∂Ω=n⋅c∇δEδu∣∣∂Ω=0,

where is the outward normal of .

To solve the AC/CH system (1), many early works [3, 4, 18, 19, 30] made use of explicit time stepping schemes, in which very small time steps were usually used due to the severe stability restriction, thus making the explicit methods impractical, especially for long time simulations. To relax the restriction of time step size, an implicit finite element framework together with Jacobian-free Newton Krylov method were introduced to solve the coupled AC/CH system in [25, 27], in which the free energy defined in (2) was replaced by simplified ones, such as the double well free energy functional. In , a fully implicit method based on Newton–Krylov–Schwarz (NKS) algorithm was employed to solve the coupled AC/CH system with the free energy defined in (2). Those implicit methods can use relatively larger time steps, but no energy stability analysis was provided. Therefore, the choice of time step size could easily violate the free-energy dissipation law of the AC/CH systems. It is worth pointing out that only one and two-dimensional numerical experiments were conducted in the aforementioned works, except for . Therefore, it is of great interest to study how to efficiently solve the three-dimensional AC/CH system, with a scheme that is energy stable and obeys the free-energy dissipation law.

As is well known, the AC/CH system can be viewed as a gradient flow driven by a free energy. For gradient flows, there exist several popular approaches, such as convex splitting [5, 13, 14, 21], stabilization [35, 23], exponential time differencing , invariant energy quadratization [32, 34], scalar auxiliary variable , and so on. However, due to the existence of the logarithmic function in the free energy of the AC/CH system, it is quite difficult to extend the above approaches to the AC/CH system. Recently, a new approach for gradient flow problems, namely the discrete variational derivative (DVD) method, has been proposed and successfully applied to Cahn–Hilliard, Allen–Cahn, and phase field crystal equations [16, 28]. The most important advantage of DVD method is that the numerical scheme, by careful design, is able to keep several important properties of the original system, for instance, free-energy dissipation, free-energy conservation, mass conservation, etc.

However, with the existence of the logarithmic term in the free energy functional, the discrete variational derivatives obtained by a DVD method will introduce terms in rational forms, of which the denominators are usually very close to zero. Due to this difficulty, numerical calculation of the discrete variational derivatives is often unstable and inaccurate. As a result, the DVD method is not directly applicable to the AC/CH system. In this work, we propose a method to approximate the discrete variational derivatives of the AC/CH system so that the numerical instability of the DVD discretization is avoided. With the approximation, an implicit scheme is constructed, which is unconditionally energy stable and therefore obeys the energy dissipative law. Here, unconditionally energy stability means that the free energy is nonincreasing in time, regardless of the time step size [26, 29]. To solve the large sparse nonlinear algebraic system arising at every implicit time step, we present a parallel, highly scalable, NKS algorithm . Due to the multiple time scales exhibited by the AC/CH system, using fixed time step size is no longer practical. To that end, we propose an adaptive time stepping strategy modified from [28, 33], and verify the efficiency of the proposed method by a series of experiments.

The remainder of this paper is organized as follows. In Sec. 2, an implicit unconditionally energy stable scheme for the AC/CH is constructed. In Sec. 3, we introduce the NKS algorithm together with the adaptive time stepping strategy to solve the nonlinear system. Experiment results on several two or three dimensional test cases are reported in Sec. 4 and some concluding remarks are given in Sec. 5.

## 2 Discretization for the Allen–Cahn/Cahn–Hilliard system

First, we rewrite the AC/CH system as

 ∂U∂t=−AG, (3)

in which , , and represents the variational derivative of the free energy functional with respect to . The local free energy is decomposed into three parts, which are

 e1(U)=α2u(1−u)−β2v2, (4) e2(U)=θ(Φ(u+v)+Φ(u−v)), e3(∇UT)=γ2(|∇u|2+|∇v|2).

In (4), and represent the bulk energy and the interfacial energy, respectively. The existence of the logarithmic term in is due to an ideal mixing of entropy for binary alloy . Here and should satisfy the following restrictions: , , and . In the AC/CH system, is defined as

 δEδu=∂e∂u−∇⋅∂e∂(∇u)=−α(u−1/2)+θΦ′(u+v)+θΦ′(u−v)−γΔu, (5) δEδv=∂e∂v−∇⋅∂e∂(∇v)=−βv+θΦ′(u+v)−θΦ′(u−v)−γΔv.

It is easy to check that is a semi-positive operator for , .

By denoting , the following integral relationship between the variational derivative and the free energy holds

 E(U+δU,∇UT+δ∇UT)−E(U,∇UT)=∫Ωeδ−e(U,∇UT)dx (6) ≈∫Ω(∂e1∂U⋅δU+∂e2∂U⋅δU+(∂e3∂∇u⋅δ∇u+∂e3∂∇v⋅δ∇v))dx =∫Ω(∂e1∂U+∂e2∂U−γ(Δu,Δv)T)⋅δUdx+B =∫ΩG⋅δUdx+B.

The first equality is obtained by using Taylor expansion, and the second equality comes from the integration-by-parts formula. Here represents the boundary terms from the integration-by-parts formula and equals zero due to the given boundary conditions. Equation (6) shows the connection between the variational derivative and the free energy, and plays an important role in the construction of the energy-stable numerical scheme for the AC/CH system.

The free energy of the AC/CH system satisfies the following equation

 ddtE(U,∇UT) =ddt∫Ωe(U,∇UT)dx (7) =∫Ω(∂e∂U⋅∂U∂t+∂e∂(∇u)⋅∂(∇u)∂t+∂e∂(∇v)⋅∂(∇v)∂t)dx =∫ΩG⋅∂U∂tdx+B1=−∫ΩG⋅AGdx+B1,

where is the boundary terms due to integration-by-parts, which vanishes with the given boundary conditions. Since is a semi-positive operator, it follows that . Thus we have

 ddtE(U,∇UT)≤0, (8)

which demonstrates the energy dissipation of the AC/CH system.

Without loss of generality, consider discretizing the AC/CH system on , which is covered by a uniform mesh with mesh sizes and . The temporal interval is split by a set of nonuniform points with time steps . Denote as the approximate solution of AC/CH system at the -th time step, where , . Throughout the paper, notations with tilde are the corresponding approximate solutions or functions at the discrete level. Let us introduce some useful notations as following

 [D+x~ϕ]ni,j=~ϕni+1,j−~ϕni,jΔx,  [D−x~ϕ]ni,j=~ϕni,j−~ϕni−1,jΔx,
 [Dx~ϕ]ni,j=~ϕni+12,j−~ϕni−12,jΔx,  [(D±x~ϕ)2]ni,j=([D+x~ϕ]ni,j)2+([D−x~ϕ]ni,j)22.

are defined similarly. In what follows, let us denote and as the approximations of scalar value functions and , respectively. We denote the discrete gradient operator as . The Laplacian operator is discretized by . The operator in is discretized by where the approximated value of the mobility is calculated as

 ~cni+12,j=c(~uni+1,j,~vni+1,j)+c(~uni,j,~vni,j)2.

The discrete operator for the AC/CH system is then defined as . With the periodic boundary conditions or the homogeneous Neumann boundary conditions, we present two vital formulas

 −Nx∑i=1~φni,j[(Dx~cDx)~ϕ]ni,j=12Nx∑i=1~cni+12,j[D+x~φ]ni,j[D+x~ϕ]ni,j+~cni−12,j[D−x~φ]ni,j[D−x~ϕ]ni,j, (9) −Ny∑i=1~φni,j[(Dy~cDy)~ϕ]ni,j=12Ny∑i=1~cni,j+12[D+y~φ]ni,j[D+y~ϕ]ni,j+~cni,j−12[D−y~φ]ni,j[D−y~ϕ]ni,j.

Next, we prove (9) in the case of homogeneous Neumann boundary conditions. The proof with periodic boundary conditions can be obtained similarly. Let us denote and be ghost points outside the computational domain. Assume , . The homogeneous Neumann boundary conditions in the -direction are discretized as

 ~ϕn1,j−~ϕn0,j=0,~ϕnNx+1,j−~ϕnNx,j=0.

The first equation of (9) is derived as follows

 −Nx∑i=1~ϕni,j[(Dx~cDx)~ϕ]ni,j=−Nx∑i=1~ϕni,j~cni+12,j(~ϕni+1,j−~ϕni,j)−~cni−12,j(~ϕni,j−~ϕni−1,j)(Δx)2 =~ϕn1,j~cn12,j(~ϕn1,j−~ϕn0,j)(Δx)2−~ϕnNx,j~cnNx+12,j(~ϕnNx+1,j−~ϕnNx,j)(Δx)2 +12Nx−1∑i=1~cni+12,j(~ϕni+1,j−~ϕni,j)2(Δx)2+12Nx∑i=2~cni−12,j(~ϕni,j−~ϕni−1,j)2(Δx)2 =12Nx∑i=1~cni+12,j[D+x~ϕ]ni,j[D+x~ϕ]ni,j+~cni−12,j[D−x~ϕ]ni,j[D−x~ϕ]ni,j.

The summation-by-parts formula in -direction can be obtained similarly.

It follows from the Cauchy–Schwarz inequality and (9) that

 (10)

where . Thus, we conclude that the discrete operator for the AC/CH system is semi-positive.

With the aforementioned notations, the discretizations of the local free energy and the free energy at time are respectively defined as

 [~e]ni,j= e1(~uni,j,~vni,j)+e2(~uni,j,~vni,j)+γ2([(D±x~u)2]ni,j+[(D±y~u)2]ni,j) (11) +γ2([(D±x~v)2]ni,j+[(D±y~v)2]ni,j).

and

 ~En:=⟨[~e]ni,j⟩=Nx∑i=1Ny∑j=1[~e]ni,jΔxΔy, (12)

where is defined as .

To obtain a full discretization scheme for the AC/CH system, a discrete form of the variational derivative is needed. By taking and , we obtain . We then choose the discrete variational derivative such that the following summation formula exactly holds

 (13)

where expresses the discrete variational derivative. Equation (13) can be viewed as a discrete form of (6). According to (4) and (13), we can finally reduce the discrete variational derivative to the following form

 ~Gi,j:=G1(~Un+1i,j,~Uni,j)+G2(~Un+1i,j,~Uni,j)+G3([∇~UT]n+1i,j,[∇~UT]ni,j),

where are derived from and in local free energy, respectively. Since both and are polynomials, we can obtain that and are also in polynomial forms as follows

 G1(~Un+1i,j,~Uni,j):=⎛⎜ ⎜⎝−α2(~un+1i,j+~uni,j−1)−β2(~vn+1i,j+~vni,j)⎞⎟ ⎟⎠, (14)
 G3([∇~UT]n+1i,j,[∇~UT]ni,j):=⎛⎜ ⎜⎝−γ2([Δ~u]n+1i,j+[Δ~u]ni,j)−γ2([Δ~v]n+1i,j+[Δ~v]ni,j)⎞⎟ ⎟⎠. (15)

However, since is not a polynomial of , we cannot choose as a polynomial such that (13) exactly holds. Consider an alternative form of as

 G2(~Un+1i,j,~Uni,j)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝θΦ(pn+1ij)−Φ(pnij)pn+1ij−pnij+θΦ(qn+1ij)−Φ(qnij)qn+1ij−qnijθΦ(pn+1ij)−Φ(pnij)pn+1ij−pnij−θΦ(qn+1ij)−Φ(qnij)qn+1ij−qnij⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠ (16)

such that . Here and .

Using the trapezoidal rule at the half-time level, we obtain fully discretized scheme for the AC/CH system as

 ~Un+1i,j−~Uni,jΔtn=−~A~Gi,j. (17)

The stability of the proposed scheme (17) is given by the following theorem. For any given time step , the numerical scheme (17) is unconditionally energy stable and the solution of scheme (17) satisfies the energy dissipative law Similar to the proof of Theorem 2 in , according to (10), (13), and (17), the discrete free energy induced by scheme (17) satisfies

 ~En+1−~EnΔtn=⟨~Gi,j⋅~Un+1i,j−~Uni,jΔtn⟩=−⟨~Gi,j⋅~A~Gi,j⟩≤0. (18)

Since the total free energy functional is non-convex, the existence and uniqueness of the solution for system (17) (especially for large values of ) are not immediate. Nonetheless, the proof of Theorem 2.1 did not require a unique solution to (17). In fact, even if one could prove that scheme (17) is unconditionally energy stable, it is numerically unstable when applied to solve the AC/CH system. The reason is that the numerical computation of fraction is unstable and inaccurate when is close to zero. Unfortunately, in the AC/CH system, or is often close to zero, which indicates that the numerical calculation of by (16) is numerically unstable and inaccurate in this situation. As a result, the scheme (17) may lead to an inaccurate or non-physical solution with the existence of the complicated function . To overcome this difficulty, we propose an approach to calculate , which is numerically stable and highly accurate. In this approach, and are approximately calculated by Taylor expansion at the point for very small . Analogous methods can be applied to terms related to and . The accuracy of the approximation is guaranteed by the following lemma, which can be verified directly from the Taylor expansion. Let us assume . If , we have the following expansion

 f(x,σ) =ln(x)−16(σx)2−⋯−1(2S+1)2S(σx)2S+O(1S2(σx)2S+2) (19) :=fS(x,σ)+RS.

Here the truncation error as .

With the Taylor expansion (19), is approximated as

 GS2(~Un+1i,j,~Uni,j)=(θϕSi,j+θψSi,jθϕSi,j−θψSi,j):=G2(~Un+1i,j,~Uni,j)−RSi,j, (20)

where and are given as follows

 ϕSi,j:=fS(ζp,σp)+fS(1−ζp,−σp),ψSi,j:=fS(ζq,σq)+fS(1−ζq,−σq). (21)

Here , , , and . Since , , , and , one can obtain as from Lemma 2.

With the approximation, scheme (17) for the AC/CH system is replaced with

 ~Un+1i,j−~Uni,jΔtn=−~A~GSi,j, (22)

where the approximate discrete variational derivative is defined as

 ~GSi,j:=G1(~Un+1i,j,~Uni,j)+GS2(~Un+1i,j,~Uni,j)+G3([∇~UT]n+1i,j,[∇~UT]ni,j). (23)

It should be noted that the solution obtained by using scheme (22) is usually different from the one by scheme (17) due to the different approximation made by the two schemes. For simplicity, we still denote the solution of scheme (22) as . According to (13) and (20), we have

 ~En+1−~En=⟨~Gi,j⋅(~Un+1i,j−~Uni,j)⟩=⟨(~GSi,j+RSi,j)⋅(~Un+1i,j−~Uni,j)⟩. (24)

The stability of the proposed scheme (22) is given by the following theorem.

Given , scheme (22) is unconditionally energy stable and the solution satisfies the following energy dissipation relationship with a cut off error

 ^En+1≤^En+Δtn4⟨RSi,j⋅~ARSi,j⟩. (25)

Furthermore, there exists an integer such that

 ^En+1≤^En, (26)

holds for any , which implies that the solution of scheme (22) obeys the energy dissipative law. According to (22), and (24), the discrete free energy decided by scheme (22) satisfies

 ~En+1−~EnΔtn=⟨~Gi,j⋅~Un+1i,j−~Uni,jΔtn⟩=−⟨~Gi,j⋅~A~GSi,j⟩ (27) =−⟨(~Gi,j−12RSi,j)⋅~A(~Gi,j−12RSi,j)−14RSi,j⋅~ARSi,j⟩ ≤⟨14RSi,j⋅~ARSi,j⟩,

which completes the proof of (25). Fu