Multicomponent systems, such as alloys, soft matters, are an important class of materials, in particular for technical applications and processes. The microstructure formation of a material plays a central role for a broad range of industrial application, such as the mechanical property of the quality and the durability, optical device, high-capacity data storage devices garcke1999multiphase ; nestler2005multicomponent ; nestler2011phase ; ofori2013multicomponent ; taha2019phase . Advances in modelling and computation have significantly improved the understanding of the fundamental nature of microstructure and phase selection processes. Notable contributions have been made through using the phase field methodology, which has been successful at examining mesoscale microstructure evolution over diffusive time scales. Recently, phase field crystal (PFC) models have been proposed to efficiently simulate eutectic solidification and elastic anisotropy, solute drag, quasicrystal formation, solute clustering and precipitation mechanisms chen2002phase ; nestler2005multicomponent ; steinbach2013phase .
The PFC model for a general class of multicomponent systems is formulated consisting components in dimensional space. The concentrations of the components are described by vector-valued functions . The variable , so-called the order parameter, denotes the local fraction of phase . The free energy functional of PFC model of a -component system can be described by two contributions, a bulk free energy and an interaction potential , which drive the density fields to become ordered by creating minimal in the free energy for these states. Formally, we can write the free energy functional of the multicomponent system as
where are relevant physical parameters. can be polynomial type or log-type formulation brazovskii1975phase ; swift1977hydrodynamic ; cross1993pattern and is the interaction potential, such as high-order differential terms or convolution terms brazovskii1975phase ; savitz2018multiple . Usually, some constraint conditions shall be imposed on the PFC model, such as the mass conservation or incompressibility which means the order parameter belong to a feasible space.
To understand the fundamental nature of multicomponent systems, it often involves finding stationary states corresponding to ordered structures and comparing free energy values of different candidate structures to construct phase diagrams. Denote to be a feasible space of -th order parameter, the above problem is transformed into solving the minimization problem
with different physical parameters
, which brings a tremendous computational burden. Most existing numerical methods for computing the stationary states of PFC models can be classified into two categories. One is to solve the steady nonlinear Euler-Lagrange system of (2) through different spatial discretization approaches. The another class aims at solving the nonlinear gradient flow equations which is able to describe the quasi-equilibrium phase behavior. In these approaches, great efforts have been made for keeping the energy dissipation which is crucial for convergence. Typical energy stable schemes to gradient flows include convex splitting and stabilized factor methods, and recently developed invariant energy quadrature, and scalar auxiliary variable approaches for a modified energy xu2006stability ; wise2009energy ; shen2010numerical ; yang2016linear ; shen2018scalar . Numerically, a gradient flow equation is discretized in both space and time domains via different discretization techniques and the stationary state is obtained with a proper choice of initialization. To the best of our knowledge, most works focus to develop these methods for binary or ternary multicomponent phase field models badalassi2003computation ; boyer2011numerical ; yang2017numerical . Recently, the scalar auxiliary variable method is introduced for general multiple energy functionals shen2018scalar and successfully applied for two-phase flow yang2019unconditionally ; qiao2019new . Besides, in the time evolution process, a suitable time step has a great impact on the performance of the methods. A large time step may cause the divergence, while a small one usually leads to slow convergence. In the most existing methods, the time step is usually fixed, or obtained by some empirical formulas qiao2011adaptive . In this work, we focus on the multicomponent PFC model and propose a class of efficient, fast and robust numerical methods for solving (2) with desired dissipation and conservation properties. Meanwhile, we propose an adaptive time step approach by using the line search technique.
Under an appropriate spatial discretization scheme, the infinite dimensional problem (2) can be formulated as a minimization problem in a finite dimensional space. Thus, there may exist alternative numerical methods that can converge to the steady states quickly by using modern optimization techniques. For example, similar ideas have shown success in computing steady states of the Bose-Einstein condensate wu2017regularized , one component PFC models jiang2020efficient , and the calculation of density functional theory ulbrich2015proximal ; liu2015analysis . Due to the block structure in the multicomponent PFC model, we propose an accelerated block Bregman proximal gradient method for computing the stationary states of (2). Our method updates each order parameter function, corresponding to the block, with an adaptive time step size by line search method initialized with Barzilai-Borwein step barzilai1988two . The sequence of blocks can be updated either deterministically cyclic or randomly shuffled for each iteration. The convergence analysis only needs that each block is updated at lease once in every fixed number of iterations. Meanwhile, the convergence of the proposed algorithm is provided without the requirement of the global Lipschitz restriction on the gradient of bulk energy.
The rest of this paper is organized as follows. Section 2 presents a concrete multicomponent PFC model, i.e. the coupled-mode Swift-Hohenberg (CMSH) model, and a spacial discretization formulation based on the projection method. In section 3, we present the adaptive block Bregman proximal gradient (ABBPG) method, including non-monotone and monotone versions, for solving the constrained non-convex multi-block problems with proved convergence. In section 4, the proposed approaches are applied to the CMSH model with two choices of Bregman distance. Numerical results are reported in section 5 to illustrate the efficiency and accuracy of our algorithms.
2 Problem formulation
2.1 Coupled-mode Swift-Hohenberg model
In this work, we consider the coupled-mode Swift-Hohenberg (CMSH) model of multicomponent systems which extends classical Swift-Hohenberg model from one length scale to multiple length scales swift1977hydrodynamic ; jiang2020stability ; jiang2016stability . The CMSH model allows the study of the formation and relative stability of periodic crystals and quasicrystals. Define the integral average to be
where is a bounded domain and is a ball centered at origin with radii . The free energy of the CMSH model for component system is
where is the -the order paramter, is the -th characteristic length scale, is interaction intensity related to the physical conditions, and is the index set defined as
Moreover, to conserve the number of particles in the system, we are working in a canonical ensemble. That is, the average of each order parameter satisfies
Theoretically, the ordered patterns including periodic and quasiperiodic structures correspond to local minimizers of the free energy functional (4) with respect to order parameters . Thus, denote
we focus on solving the minimization:
Next, under an appropriate spatial discretization, we formulate this infinite dimensional problem as a minimization problem over a finite dimensional space.
2.2 Projection method discretization
In this subsection, we apply the projection method jiang2014numerical , which is a general numerical framework to study the quasicrystals and periodic crystals, to discretize the CMSH free energy functional. By the projection method, any -dimensional quasicrystal can be embedded into an -dimensional periodic structure (). Using the projection matrix , any -dimensional quasiperiodic structure over the whole space can be expanded as
where is invertible, relevant to the -dimensional primitive reciprocal lattice. The associated Fourier coefficient can be obtained by
Similarly, the order parameter in multicomponent systems can be expanded as follows
We truncate the Fourier coefficients by setting
Let with . Then the discretized energy functional can be stated as
where and is a diagonal matrix with nonnegative entries . are -dimensional convolutions in the reciprocal space. A direct evaluation of the nonlinear term is extremely expensive. Thanks that is a simple multiplication in the -dimensional real space. The pseudospectral method takes the advantage of this observation by evaluating in the Fourier space and
in the real space via the Fast Fourier Transformation algorithm. As a result, it provides an efficient technique to reduce the computation cost.
3 The proposed method
In this section, we consider the minimization problem in the form of
where , is the feasible space of variable , with . It is easy to know that the problem (11) can be reduced to (12) by setting , and . We make the following assumptions on the objective function (12) (the notations can be found in subsection 3.1 for details):
is proper and continuously differential on but non-convex.
is proper, lower semicontinuous and convex.
is a nonempty, closed and convex set.
is bounded below and level bounded.
For all , where is the closed ball that contains the sub-level set .
Before proceeding to present the numerical algorithm and analysis for solving (12), we first introduce some notations and useful definitions used in the following analysis. Then, we give an abstract framework of the first order method with proved convergence. Two kinds of concrete numerical algorithms to the CMSH model based on the abstract formulation for solving (11) will be presented in the section 4.
3.1 Notations and definitions
We denote and the projection operator is defined as . For a subset , . Let be the -th continuously differential functions on . The domain of a function is defined as and the relative interior of is defined as , where is the smallest affine set that contains and . is proper if and . For , is the -(sub)level set of . We say that is level bounded if is bounded for all . is lower semicontinuous if all level set of is closed. The subgradient of at is defined as . For , we denote . Table 1 summarizes the notations used in this section.
|the total number of blocks|
|the update block selected at the -th iteration|
|the number of updates to within the first iterations|
|the value of after the -th iteration|
|the value of after -th update|
|the value of extrapolation point at the th iteration|
|the extrapolation weight used at the th iteration|
|the value of|
|the partial gradient of with respect to|
Throughout this paper, we assume that is a strongly convex function. Then we have the following definitions. [Bregman divergence li2019provable ] The Bregman divergence with respect to is defined as
It is noted that and if and only if due to the strongly convexity of . Moreover, is also strongly convex for any fixed .
[Relative smoothness bauschke2016descent ] A function is called -smooth relative to if there exists such that is convex for all . The next Lemma 1 establishes the relationship between the Bregman divergence and the relative smoothness condition. [bauschke2016descent ] If is -smooth relative to , then
To deal with the problem of form (12), we generalize the definition of relative smoothness to block-wise function as the definition 2.4 in ahookhosh2019multi . [Block-wise relative smoothness] For a block-wise function , we call is -smooth relative to if for each -th block and fixed ,
is -strongly convex and .
is -smooth relative to with respect to .
3.2 ABBPG algorithm
For the next discussion, we make an additional assumption of problem (12).
is -smooth relative to .
Given a feasible , the accelerated block Bregman proximal gradient (ABBPG) method can pick deterministically or randomly, then is updated as follows
where is the step size and is the extrapolation
The extrapolation weight for some and is the value of before it was updated to . The definitions of and can be found in Table 1.
Well-definedness of the proximal subproblem. Define the -th block Bregman proximal gradient mapping as
The next lemma shows the well-possedness of . Under Assumption 12 and Assumption 3.2, the map defined in (17) is nonempty and single-valued from to . By assumptions, and , then and is well-defined. Let
By the convexity of and , we know that is strongly convex. Thus, is coercive (bauschke2011convex , Corollary 11.17). According to the Corollary 3.23 in brezis2010functional , achieves its minimum on and the strongly convexity implies the uniqueness of minimum.
Restart technique. The extrapolation in (15) may generates the oscillating phenomenon of the objective value that slows down the convergence. We use an adaptive restart technique when the sharp oscillation is detected. More precisely, given and , define
For some constant , we let , and set if the following inequality
does not hold where is a small constant. Otherwise, we obtain via and update .
where and . Let be a constant and is obtained form (18), we adopt the step size whenever the following inequality holds
In summary, the ABBPG method is presented in Algorithm 2.
The Algorithm 2 becomes a monotone decreasing scheme if as the (19) becomes for some . As Algorithm 1 contains a line search scheme, we will show that it can be terminated in finite iterations by using the following descent property. [Sufficient decrease] Under Assumption 12 and Assumption 3.2. Let , then generated by (18) satisfies
where is the strong convexity coefficient of . By the convexity of and Lemma 3.2, it follows that . Then
The first inequality is obtained by the definition of Bregman proximal subproblem, the second inequality holds since the block-wise relative smoothness of and the last inequality holds from the -convexity of . Let and . Lemma 2 shows that the line search condition (21) is satisfied once we set . Moreover, a direct consequence of Lemma 2 is that if , then . In fact, we have and , then the line search inequality (21) reduces to a stronger condition than the non-restart condition (19). The above arguments guarantee that Algorithm 2 is well defined.
3.3 Convergence analysis
To ensure the convergence, we make a mild assumption, for the update manner of Algorithm 2. A similar assumption has been used in xu2017globally , With any consecutive iterations, each block should be updated at least once, i.e. for any , it has . Let be the sequence generated by Algorithm 2. Then, we have and is non-increasing. By the definition of , it is easy to know . If the non-restart condition (19) does not hold, we know which implies
If the non-restart condition (19) holds, we have
Since , the sequence generated by Algorithm 2 is contained in . From the Assumption 12, we know is compact. Together with Lemma 3.2, we know where is the closed ball that contains . The next lemma establishes that is Lipschitz continuous on .
Suppose Assumption 12 and Assumption 3.2 hold. Then, there exists such that is -Lipschitz continuous on for all . Since is a compact subset of , then is closed. It’s easy to know that is bounded by the fact that . Thus, is -Lipschitz continuous on (Corollary 8.41 in bauschke2011convex ). As a result, we get
Together with the fact that , we conclude that there exists suth that is -Lipschitz continuous on the compact set .
Assume , we prove by induction that the following relations hold for any finite :
Substituting by in (26), we obtain
where the last inequality holds since . Let , we get
Since is -Lipschitz continuous on and , we have