1 Introduction
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, highcapacity 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 vectorvalued functions . The variable , socalled 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
(1) 
where are relevant physical parameters. can be polynomial type or logtype formulation brazovskii1975phase ; swift1977hydrodynamic ; cross1993pattern and is the interaction potential, such as highorder 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
(2) 
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 EulerLagrange 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 quasiequilibrium 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 twophase 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 BoseEinstein 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 BarzilaiBorwein 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 coupledmode SwiftHohenberg (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 nonmonotone and monotone versions, for solving the constrained nonconvex multiblock 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 Coupledmode SwiftHohenberg model
In this work, we consider the coupledmode SwiftHohenberg (CMSH) model of multicomponent systems which extends classical SwiftHohenberg 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
(3) 
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
(4) 
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
(5) 
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
(6) 
we focus on solving the minimization:
(7) 
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
(8) 
Similarly, the order parameter in multicomponent systems can be expanded as follows
(9) 
We truncate the Fourier coefficients by setting
Let with . Then the discretized energy functional can be stated as
(10) 
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
(12) 
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 nonconvex.

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 sublevel 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.
Notion  Definition 

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
(13) 
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
(14) 
To deal with the problem of form (12), we generalize the definition of relative smoothness to blockwise function as the definition 2.4 in ahookhosh2019multi . [Blockwise relative smoothness] For a blockwise 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
(15) 
where is the step size and is the extrapolation
(16) 
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.
Welldefinedness of the proximal subproblem. Define the th block Bregman proximal gradient mapping as
(17) 
The next lemma shows the wellpossedness of . Under Assumption 12 and Assumption 3.2, the map defined in (17) is nonempty and singlevalued from to . By assumptions, and , then and is welldefined. 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.
In Assumption 12, is convex for all . Thus, Lemma 3.2 implies that the iteration in (15) is welldefined as long as is set in the initialization.
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
(18) 
For some constant , we let , and set if the following inequality
(19) 
does not hold where is a small constant. Otherwise, we obtain via and update .
Step size estimation.
In each step, we add a nonmonotone backtracking line search method for finding the appropriate step which is initialized by the similar idea of BB step barzilai1988two , i.e.(20) 
where and . Let be a constant and is obtained form (18), we adopt the step size whenever the following inequality holds
(21) 
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
(22) 
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 blockwise 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 nonrestart 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 nonincreasing. By the definition of , it is easy to know . If the nonrestart condition (19) does not hold, we know which implies
(23) 
If the nonrestart condition (19) holds, we have
(24) 
Combing (23) with (24), we obtain
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 .
Suppose Assumption 12 and Assumption 3.2 hold. Let obtained by Algorithm 2, there exists such that
(25) 
Since is nonincreasing from Lemma 3.3 and is bounded below, there exists such that . Define and combine (23) with (24), we have
(26) 
Assume , we prove by induction that the following relations hold for any finite :
(27) 
Substituting by in (26), we obtain
where the last inequality holds since . Let , we get
Since is Lipschitz continuous on and , we have
Thus, one has , which implies that (27) holds for . Suppose now that (27) holds for some , we show that it also holds for . Substituting by in (26), it gives