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:
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
where the local free energy is
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 
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
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
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
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
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
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
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
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
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
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
The first equation of (9) is derived as follows
The summation-by-parts formula in -direction can be obtained similarly.
It follows from the Cauchy–Schwarz inequality and (9) that
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
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
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
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
However, since is not a polynomial of , we cannot choose as a polynomial such that (13) exactly holds. Consider an alternative form of as
such that . Here and .
Using the trapezoidal rule at the half-time level, we obtain fully discretized scheme for the AC/CH system as
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
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
Here the truncation error as .
With the Taylor expansion (19), is approximated as
where and are given as follows
Here , , , and . Since , , , and , one can obtain as from Lemma 2.
With the approximation, scheme (17) for the AC/CH system is replaced with
where the approximate discrete variational derivative is defined as
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
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
Furthermore, there exists an integer such that