The phase field crystal (PFC) model is an important approach to describe many physical processes and material properties, such as the formation of ordered structures, nucleation process, crystal growth, elastic and plastic deformations of the lattice, dislocations, etc chen2002phase ; provatas2010phase . More concretely, let the order parameter function be , the PFC model can be expressed by a free energy functional:
where are the physical parameters, is the interaction energy with polynomial type or log-type formulation and is the bulk energy that contains higher-order linear operators to form ordered structures brazovskii1975phase ; lifshitz1997theoretical ; swift1977hydrodynamic . A typical interaction potential function for a bounded domain is
which can be used to describe the pattern formation of periodic crystals, quasicrystals and multi-polynary crystals.
In order to understand the theory of PFC models as well as predict and guide experiments, it requires to find stationary states and construct phase diagrams of the energy functional (1). Mathematically, denote to be a feasible space, one should solve the minimization problem
with different physical parameters , which brings the tremendous computational burden. Therefore, within appropriate spatial discretization, the goal of this paper is to develop an efficient and robust numerical method for solving (3) with guaranteed convergence.
Most existing numerical methods for computing the stationary states of PFC model can be classified into two categories. The first class of numerical methods solves the steady nonlinear Euler-Lagrange equations of (3) through different spatial discretization approaches. The second class of numerical methods has been designed via the formulation of nonlinear gradient flow equations. In these numerical PDE approaches, the time-dependent nonlinear gradient flows are discretized in space via different numerical methods. In these time discretized approaches, great efforts have been made to keep the energy dissipation which is crucial for convergence. Typical energy stable schemes to the gradient flows include convex splitting and stabilized factor methods, and recently developed invariant energy quadrature, and scalar auxiliary variable approaches for a modified energy shen2018scalar . It is noted that the gradient flow approach is able to describe the quasi-equilibrium behavior of PFC systems. Numerically, the gradient flow is discretized in both space and time domain via different discretization techniques and the stationary state is obtained with a proper choice of initial data.
Under an appropriate spatial discretization scheme, the infinite dimensional problem (3) can be formulated as a minimization problem over a finite dimensional space. Thus, there may exist alternative numerical methods that can converge to the steady states quickly by using modern optimization techniques. Similar ideas have been shown success in computing steady states of the Bose-Einstein condensate wu2017regularized and the calculation of density functional theory ulbrich2015proximal ; liu2015analysis . In the PFC models, the discretized energy is nonlinear and non-convex which consists of two parts: bulk energy and interaction energy. Motivated by the semi-implicit scheme and the accelerated proximal gradient (APG) method beck2009fast ; tseng2008accelerated
which has been successfully applied in image processing and machine learning, we propose an efficient numerical method for calculating the steady states of (3). As the traditional APG method is proposed for convex problem and its oscillation phenomenon slows down the convergence o2015adaptive ; su2014differential , the restart scheme has been used for accelerating the convergence. Moreover, the numerical speed can be further accelerated by using the line search starting with Barzilai-Borwein steps bao2018coherence . The connection of classical explicit/implicit schemes in gradient flows and proximal gradient methods is also built by defining a generalized proximal operator. Extensive numerical experiments have demonstrated that our approach can quickly reach the vicinity of an optimal solution with moderately accuracy, even for very challenge cases. As a byproduct, our numerical result reaveals a deep physical intension of a simple PFC model, the Landau-Brazovskii (LB) model, by obtaining the sigma phase.
The rest of this paper is organized as follows. Different discretizations of the energy functional via the Fourier pseudospectral approach and the projection method are introduced in section 3. In section 4, we present the gradient type method and the adaptive APG method for solving the discretized minimization problem. The connection between our proposed approach and some existing time discretized schemes in these numerical methods for solving gradient-flow equations has been built in section 5. Numerical results are reported in section 6 to illustrate the efficiency and accuracy of our algorithms. Finally, some concluding remarks are given in section 7.
2 Physical models
Two classes of PFC models are considered in the paper. The first one is the Landau-Brazovskii (LB) model which describes periodic structures brazovskii1975phase
. The LB model was introduced to investigate the character of phases and phase transition of periodic crystals. It has been discovered in many different scientific fields, e.g., polymeric materialsshi1996theory . In particular, the energy functional of LB model is
where is a real-valued function which measures the order of system in terms of order parameter. is the system volume, is the bare correlation length, is the dimensionless reduced temperature, is phenomenological coefficient. Compared with double-well bulk energy, the cubic term in the LB functional helps us study the first-order phase transition.
The second one is the Lifshitz-Petrich (LP) model that can simulate quasiperiodic structures, such as the bi-frequency excited Faraday wave lifshitz1997theoretical , and the explanation of the stability of soft-matter quasicrystals lifshitz2007soft ; jiang2015stability . Before we present the LP model, an introduction of the average spacial integral, so-called almost periodic integral, is necessary. For a space-filling structure, such as the quasicrystal, the average spacial integral can be defined as
where is the ball centred at origin with radii . Using the above notation, the energy functional of LP model is given by
is the energy penalty, and are phenomenological coefficients.
Formally, the difference between the LB and LP energy functional is the number of length-scale governed by the differential term. The LB model has a one-length-scale which can be used to study the phase behavior of periodic structures brazovskii1975phase ; zhang2008efficient , while the LP model possesses a two-length-scale that can be used to studied the formation and stability of quasicrystals lifshitz1997theoretical ; lifshitz2007soft ; jiang2015stability ; dotera2014mosaic .
3 Discretization of the energy functional
In this section, we introduce different discretization schemes of the energy functionals (4) and (6), and reduce them to finite dimensional minimization problems. Two classes of stationary states are considered. The first class of stationary states is periodic in LB model which can be described in a bounded domain. Thus we can truncate the energy functional from the whole space to a bounded domain with periodic boundary condition. Then we employ Fourier pseudospectral method to discretize LB energy functional. The second class of stationary phases can be quasicrystals in LP model. For these structures, the discretization of the energy functional in a bounded domain results in a significant Diophantine approximation error. In this paper, we apply with the projection method jiang2014numerical , a high dimensional interpretation approach, to discretize the LP energy function (6), which can avoid the Diophantine approximation error.
3.1 Fourier pseudospectral discretization
Each of the -dimensional periodic system can be described by a Bravis lattice
where the vectorforms the primitive Bravis lattice . The smallest possible periodicity, or named the unit cell, of the system is
The associated reciprocal lattice is
The primitive reciprocal lattice vector satisfies the dual relationship
Then the periodic function on the Bravis lattice, i.e., , can be expanded as
where , is invertible. The coefficient, , satisfies
In numerical computations, we need to minimize the LB energy functional (4) in a finite dimensional subspace. More precisely, let , and
The number of elements in the set is . The order parameter can be projected into the finite dimensional space , i.e.,
Due to the orthonormal condition
the LB energy functional can be approximated as
where , Fourier coefficient , and
. The convolutions in the above expression can be calculated by Fourier pseudospectral method through the fast Fourier transform (FFT). Therefore, it reduces to a finite dimensional minimization problem:
where and are the discretized interaction and bulk energy, respectively. The gradient of is
where is a diagonal matrix with entries and is the discretized Fourier transform matrix.
3.2 Projection method discretization
For the -dimensional quasicrystals which are the space-filling structures, the spatial integral in the energy functional (1) shall be instead by the almost periodic spacial integral , as defined by Eq. (5). We immediately have the following orthonormal property:
For an almost periodic function, the average transformation is
and it is well defined katznelson2004anintroduction . In this paper, we carry out the above computation in a higher dimension using the projection method which is based on the fact that a -dimensional quasicrystal can be embedded into an -dimensional periodic structure () hiller1985crystallographic . Using the projection method, the order parameter is
where is invertible, related to the -dimensional primitive reciprocal lattice and the projection matrix depends on the property of quasicrystals, such as rotational symmetryhiller1985crystallographic . The Fourier coefficient satisfies
Again, in practice, we need to minimize the LP energy functional (4) in a finite dimensional subspace. More precisely, let , and
where , , , , . It is clear that the nonlinear (quadratic, cubic and cross) terms in Eq. (23) are -dimensional convolutions in the reciprocal space. A direct evaluation of these convolution terms is extremely expensive. Instead, these terms are simple multiplication in the -dimensional real space. Again, the efficient pseudospectral approach is applied to calculate these convolutions in Eq. (23) through the -dimensional FFT.
Therefore, it leads to the following finite dimensional minimization problem:
where and are the discretized interaction and bulk energies. The gradient of is
where is a diagonal matrix with entries . The is the discretized Fourier transform matrix and is the corresponding inverse discretized Fourier transform matrix. In the following, we will neglect the superscript of hat for simplicity.
4 The proposed numerical approach
In this section, we first review the classical semi-implicit method and accelerated proximal gradient (APG) method and then propose the adaptive APG method with proved convergence. Finally, the connection of generalized proximal operator with gradient flows approaches is present.
4.1 Semi-implicit scheme
The semi-implicit scheme is a simple but useful approach for finding the stationary state based on gradient flows. For example, the Allen-Cahn equation of the discretized energy functional is
with the periodic condition and is the spurious time variable. Given an initial value and the time step , the semi-implicit scheme is
where is the approximation of the solution at , i.e., . The semi-implicit scheme satisfies the following energy dissipation property.
Let . Assume that there exists a constant such that the bulk energy satisfies , and the time step length , then the solutions of (27) satisfy
From (27), it is easy to know
where the last inequality is from the Taylor expansion of and the boundedness constraint on . ∎
Therefore, to satisfy the energy dissipation law, the time step length depends on the Lipschitz constant . In a general PFC model, the universal Lipschitz constant may not exist or be very large in bounded domain which leads to a small time step and slows down the convergence speed. Despite its strict requirements on in theory, the semi-implicit scheme works well in practice which inspires us a further exploration of the semi-implicit scheme. In the following context, we will combine modern optimization approaches and the semi-implicit scheme to obtain a more efficient approach.
4.2 Accelerated proximal gradient (APG) method
where is the finite dimensional Hilbert space equipped with the inner product , and are both continuously convex and has a Lipschitz constant , i.e.
Given initializations and , the APG method consists of the following steps:
where and the mapping is defined as
It is noted that the proximal map in (31) is well defined as is convex. Moreover, the step size can be set adaptively as long as the following inequality holds:
The APG method has the attractive convergence property as follows.
4.3 Adaptive APG method
We omit the subscript for simplicity. However, there are two main obstacles when directly applying APG method for solving phase field models as is non-convex and has no universal Lipschitz constant. In this paper, we propose an efficient and convergent numerical algorithm for solving the discretized phase field model (16) by combining APG method with restart techniques.
The restart techniques for the APG method was proposed in o2015adaptive which has shown significant acceleration of the APG method by imposing the decreasing property of the objective value when solving convex problems. Furthermore, another restart strategy called speed restart is developed in su2014differential to ensure the linear convergence of the proposed restart APG method for strongly convex problems. In recent years, the restart techniques have been furtherer applied for solving non-convex composite problems in image processing bao2016image . We introduce the details of the proposed algorithm in the following context.
Let and be the current and previous states respectively and the extrapolation weight . We can obtain a candidate state by
It is noted that the proximal mapping in (34) is well defined as is convex. Different from the APG method, the restart technique is to determine whether we accept the result
as the new estimate. Inspired by the function value restart condition in o2015adaptive , we choose whenever the following condition holds:
for some . If the condition (36) does not hold, we restart the APG by setting . In this case, we have
the sub-level set is compact for any . Let be the closed ball that contains . From the smoothness of , is Lipschitz continuous in . Denote to be the Lipschitz constant of in the set , i.e.
Thus, we obtain the the next proposition that shows satisfying the sufficient decrease condition for all .
Given an initial point and the iterates with for . If , then
holds and thus the (37) is a safe-guard step which ensures energy dissipation. On the other hand, by Proposition 4.3, should be less than which might be very small. Thus, it only allows a small step size which may significantly slow down the convergence, and be always too conservative. By line search technique, we can adaptively estimate the step size which will be introduced as follows.
Estimation of step size . Define: , and . We initialize the search step by the Barzilai-Borwein (BB) method barzilai1988two , i.e.
and the proposed adaptive APG algorithm is present in Algorithm 2.
4.3.1 Convergence analysis
In this section, we show that our proposed method converges to a steady state of the original energy function. Firstly, we present a useful lemma for our analysis.
Lemma 4.4 (Uniformized Kurdyka-Lojasiewicz property bolte2014proximal .).
Let be a compact set and defined in (16) be bounded below. Assume that is constant on . Then, there exist , , and such that for all and all , one has,
where and .
The proof is based on the facts that and satisfy the so called Kurdyka-Lojasiewicz property on bolte2014proximal . ∎
and two sets and . It is noted that for any , we have .
Let , then there exists some such that for all as is increasing and is reset to at most every iteration. We show the following properties of the sequence generated by Algorithm 2.
Sufficient decrease property. If , we have
where . Since , there exists such that and as . This implies
Bounded the gradient. If , by the optimality condition of (34), we have
as where is a bounded set and is the Lipschitz constant of in . If , by the optimality condition of (43), we have
Then, , then
Subsequence convergence. Since which is compact, there exists a subsequence and such that
where the last two equalities are from the continuity of . Moreover, (46) implies
Then, we know from (49).
Finite length property. Let be the set of limiting points of the sequence starting from . By the boundedness of and the fact , it follows that is a non-empty and compact set. Moreover, from (45), we know is constant on , denoted by . If there exists some such that , then we have for all which is from (45). In the following proof, we assume that for all . Therefore, , there exists some such that for all , we have , i.e.
Applying lemma 4.4, for all we have
Form (49), it implies
By the convexity of , we have