# An accelerated method for computing stationary states of multicomponent phase field crystal model

Finding the stationary states of a free energy functional is an essential topic in multicomponent material systems. In this paper, we propose a class of efficient numerical algorithms to fast compute the stationary states of multicomponent phase field crystal model. Our approach formulates the problem as solving a constrained non-convex minimization problem. By using the block structure of multicomponent systems, we propose an adaptive block Bregman proximal gradient algorithm that updates each order parameter alternatively. The updating block can be chosen in a deterministic or random manner. The convergence property of the proposed algorithm is established without the requirement of global Lipschitz constant. The numerical results on computing stationary periodic crystals and quasicrystals in the multicomponent coupled-mode Swift-Hohenberg model have shown the significant acceleration over many existing methods.

## Authors

• 7 publications
• 16 publications
• 12 publications
• ### Efficient numerical methods for computing the stationary states of phase field crystal models

Finding the stationary states of a free energy functional is an importan...
02/23/2020 ∙ by Kai Jiang, et al. ∙ 0

• ### An efficient method for computing stationary states of phase field crystal models

Computing stationary states is an important topic for phase field crysta...
09/01/2019 ∙ by Kai Jiang, et al. ∙ 0

• ### Computing ground states of Bose-Einstein Condensates with higher order interaction via a regularized density function formulation

We propose and analyze a new numerical method for computing the ground s...
08/24/2019 ∙ by Weizhu Bao, et al. ∙ 0

• ### Leveraging Two Reference Functions in Block Bregman Proximal Gradient Descent for Non-convex and Non-Lipschitz Problems

In the applications of signal processing and data analytics, there is a ...
12/16/2019 ∙ by Tianxiang Gao, et al. ∙ 0

• ### Block Alternating Bregman Majorization Minimization with Extrapolation

In this paper, we consider a class of nonsmooth nonconvex optimization p...
07/09/2021 ∙ by Le Thi Khanh Hien, et al. ∙ 0

• ### On the Divergence of Decentralized Non-Convex Optimization

We study a generic class of decentralized algorithms in which N agents j...
06/20/2020 ∙ by Mingyi Hong, et al. ∙ 0

• ### Recursive multilevel trust region method with application to fully monolithic phase-field models of brittle fracture

The simulation of crack initiation and propagation in an elastic materia...
03/01/2019 ∙ by Alena Kopaničáková, et al. ∙ 0

##### 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

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

 E[{ϕi(\br)}si=1;Θ]=G[{ϕi(\br)}si=1;Θ]+F[{ϕi(\br)}si=1;Θ], (1)

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

 minE[{ϕi(\br)}si=1;Θ],~{} s.t.~{} ϕi(\br)∈Vi (i=1,2,⋯,s), (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 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

 \bbint=⎧⎨⎩1|Ω|∫Ω, for periodic % crystals,limR→∞1|BR|∫BR, for % quasicystals, (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

 E[{ϕj(\br)}sj=1]=\bbint{12s∑j=1[(∇2+q2j)ϕj(\br)]2+∑\calIs,nτi1,i2,⋯,iss∏j=1ϕijj(\br)}d\br, (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

 \calIs,n:={(i1,i1,⋯,is):ij∈N (j=1,2,⋯,s), 1≤s∑j=1ij≤n}.

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

 \bbintϕj(\br)d\br=0. (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

 Gj[ϕj]=\bbint12(∇2+q2j)ϕj(\br)]2d\br,  F[{ϕj(\br)}sj=1]=\bbint∑\calIs,nτi1,i2,⋯,iss∏j=1ϕijj(\br)d\br, (6)

we focus on solving the minimization:

 minE[{ϕj(\br)}sj=1]=s∑j=1Gj[ϕj(\br)]+F[{ϕj(\br)}sj=1]s.t.\bbintϕj(\br)d\br=0,j=1,2,⋯,s. (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

 ψ(\br)=∑\bh∈\bbZm\hpsi(\bh)ei[(\calP⋅\bB\bh)⊤⋅\br],\br∈\bbRd,

where is invertible, relevant to the -dimensional primitive reciprocal lattice. The associated Fourier coefficient can be obtained by

 \hpsi(\bh)=\bbintψ(\br)e−(\textcolorblack\calP⋅\bB\bh)⊤\brd\br. (8)

Similarly, the order parameter in multicomponent systems can be expanded as follows

 ϕj(\br)=∑\bh∈\bbZm\hphij(\bh)ei[(\calP⋅\bB\bh)⊤⋅\br],j=1,2,⋯,s. (9)

We truncate the Fourier coefficients by setting

 \hphij(\bh)=0,∀|hl|>Nl,j2,l=1,2,⋯,m.

Let with . Then the discretized energy functional can be stated as

 E({\bhphij}sj=1)=s∑j=1Gj(\bhphij)+F({\bhphij}sj=1), (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.

Let . Using the projection method discretization, the infinite dimensional problem (7) reduces to finite dimensional problem

 min\hPhi E(\hPhi)=s∑j=1Gj(\bhphij)+F(\hPhi)s.t.e⊤1\bhphij=0,j=1,2,⋯,s. (11)

In the next section, we aim at designing efficient algorithms for solving the non-convex and multi-block problem (11).

## 3 The proposed method

In this section, we consider the minimization problem in the form of

 minXE(X)=f(X)+s∑j=1gj(\bxj)s.t. \bxj∈\calSj,j=1,2,⋯,s. (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):

1. is proper and continuously differential on but non-convex.

2. is proper, lower semicontinuous and convex.

3. is a nonempty, closed and convex set.

4. is bounded below and level bounded.

5. 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.

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

 Dh(x,y)=h(x)−h(y)−⟨∇h(y),x−y⟩,∀ (x,y)∈\domh×\intdomh. (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

 f(x)−f(y)−⟨∇f(y),x−y⟩≤LfDh(x,y),∀ (x,y)∈\domh×\intdomh. (14)

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).

1. is -smooth relative to .

2. .

Given a feasible , the accelerated block Bregman proximal gradient (ABBPG) method can pick deterministically or randomly, then is updated as follows

 ⎧⎪⎨⎪⎩\bxk+1i=\bxki,if i≠bk\bxk+1i=\argmin\bz∈\calSi{gi(\bz)+⟨∇if(\byk,\bxk≠i),\bz−\byk⟩+1αkDhi(\bz,\byk)}if i=bk (15)

where is the step size and is the extrapolation

 \byk=(1+wk)\bxki−wk\bxprevi=(1+wk)\tbxnkii−wk\tbxnki−1i. (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.

Well-definedness of the proximal subproblem. Define the -th block Bregman proximal gradient mapping as

 Tjα(X):=\argmin\bz∈\calSj{gj(\bz)+⟨∇jf(X),\bz−\bxj⟩+1αDhj(\bz,\bxj)}. (17)

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

 Γj(\bz):=αgj(\bz)+α⟨∇jf(X),\bz−\bxj⟩+Dhj(\bz,\bxj).

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 well-defined 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

 \bzk=Γiαk(\byk,\bxk≠i)=\argmin\bz∈\calSi{gi(\bz)+⟨∇if(\byk,\bxk≠i),\bz−\byk⟩+1αkDhi(\bz,\byk)}. (18)

For some constant , we let , and set if the following inequality

 E(Xmk)−E(\bzk,\bxk≠i)≥σ∥\bxki−\bzk∥2 (19)

does not hold where is a small constant. Otherwise, we obtain via and update .

Step size estimation.

In each step, we add a non-monotone backtracking line search method for finding the appropriate step which is initialized by the similar idea of BB step barzilai1988two , i.e.

 αk=⎧⎪ ⎪⎨⎪ ⎪⎩α0,wk=0,⟨uki,uki⟩⟨uki,vki⟩ or ⟨vki,uki⟩⟨vki,vki⟩,wk≠0, (20)

where and . Let be a constant and is obtained form (18), we adopt the step size whenever the following inequality holds

 max(E(\byk,\bxk≠i),E(Xmk))−E(\bzk,\bxk≠i)≥η∥\byk−\bzk∥2. (21)

In summary, the ABBPG method is presented in Algorithm 2.

[!pbht] Algorithm 1 Estimation of at with respect to block [1] Inputs: , , and , . Initialize by (20); Calculate the smallest index such that (21) holds at defined in (18) with step size . Output: step size .

[!pbht] Algorithm 2 ABBPG method [1] Initialize , and , , . stopping criterion is not satisfied Pick in a deterministic or random manner; Update Obtain via Algorithm 1. Calculate via (18). (19) holds Set and choose . Restart by setting and . .

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

 E(\byk,\bxk≠i)=f(\byk,\bxk≠i)+gi(\byk)+∑j≠igj(\bxkj) =[f(\byk,\bxk≠i)+⟨∇if(\byk,\bxk≠i),\bz−\byk⟩+1αkDhi(\bz,\byk)+gi(\bz)]\bz=\byk+∑j≠igj(\bxkj) ≥f(\byk,\bxk≠i)+⟨∇if(\byk,\bxk≠i),\bzk−\byk⟩+1αkDhi(\bzk,\byk)+gi(\bzk)+∑j≠igj(\bxkj) ≥f(\bzk,\bxk≠i)−LjfD(\bzk,\byk)+1αkDhi(\bzk,\byk)+gi(\bzk)+∑j≠igj(\bxkj) =E(\bzk,\bxk≠i)+(1αk−Ljf)Dhi(\bzk,\byk)≥E(\bzk,\bxk≠i)+(1αk−Ljf)γi2∥\bzk−\byk∥2.

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

 E(Xmk)−E(Xk+1)≥0. (23)

If the non-restart condition (19) holds, we have

 E(Xmk)−E(Xk+1)≥σ∥Xk−Xk+1∥2≥0. (24)

Combing (23) with (24), we obtain

 E(Xmk+1)=max{E(Xj)|[k+1−M]∗≤j≤k+1)}≤max{E(Xj)|[k−M]∗≤j≤k)}=E(Xmk).

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

 ∣∣ ∣∣s∑j=1gj(\bxj)−s∑j=1gj(\byj)∣∣ ∣∣≤s∑j=1|gj(\bxj)−gj(\byj)|≤s∑j=1Ljg∥\bxj−\byj∥≤s(maxj=1,2,⋯sLjg)∥X−Y∥, ∀X,Y∈\calB(X0).

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

 limk→∞∥Xk+1−Xk∥=0,limk→∞E(Xk)=E∗. (25)

Since is non-increasing from Lemma 3.3 and is bounded below, there exists such that . Define and combine (23) with (24), we have

 E(Xmk)−E(Xk+1)≥σ∥dk∥. (26)

Assume , we prove by induction that the following relations hold for any finite :

 limk→∞∥dmk−j∥=0,limk→∞E(Xmk−j)=E∗. (27)

Substituting by in (26), we obtain

 σ∥dmk−1∥2 ≤E(Xmmk−1)−E(Xmk)≤E(Xmk−M−1)−E(Xmk),

where the last inequality holds since . Let , we get

 limk→∞∥dmk−1∥=0.

Since is -Lipschitz continuous on and , we have

 |E(Xmk−1)−E(Xmk)|≤LE∥Xmk−1−Xmk∥=LE∥dmk−1∥→0 (k→∞).

Thus, one has