# Error analysis of SAV finite element method to phase field crystal model

In this paper, we construct and analyze an energy stable scheme by combining the latest developed scalar auxiliary variable (SAV) approach and linear finite element method (FEM) for phase field crystal (PFC) model, and show rigorously that the scheme is first-order in time and second-order in space for the L 2 and H -1 gradient flow equations. To reduce efficiently computational cost and capture accurately the phase interface, we give a simple adaptive strategy, equipped with a posteriori gradient estimator, i.e. L 2 norm of the recovered gradient. Extensive numerical experiments are presented to verify our theoretical results and to demonstrate the effectiveness and accuracy of our proposed method.

## Authors

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

The phase field crystal (PFC) model[1, 2] was proposed as an approach to simulate crystals at the atomic scale but on a coarse-grained diffusive time scale. Many physical processes, such as the formation of ordered structures, phase separation of polynary systems, can be described using this model. The PFC model can also explain elastic and plastic deformations of the lattice, dislocations, grain boundaries, multiple crystal orientations and many other observable phenomena[4, 6].

There are several kinds of PFC models. In general, they can be classified into two classes according to characteristic length scale: one-length-scale and multi-length-scale. One-length-scale PFC models can be used to describe the phase behavior of periodic structures

[7, 8, 9]. Accordingly, multi-length-scale PFC models can be employed to explain the formation of quasicrystals[10, 11]. In this work, we focus on the development of numerical methods of one-length-scale PFC model. In particular, the classic Landau-Brazovskii (LB) model[7, 12, 13]

will be used to demonstrate our proposed method. The LB model was built to investigate the character of phase transition. It has been discovered in many scientific fields. For example, the LB model can be derived from more complicated self-consistent field theory of diblock copolymers

. Compared with the typical Swift-Hohenberg (SH) model with double-well bulk energy , LB energy functional includes a cubic term which can be used to study the first-order phase transition.

The (Allen-Cahn) or

(Cahn-Hilliard) gradient flow equation is usually adopted to describe the dynamic behavior of the phase-field or PFC model. These dynamic equations are time-dependent nonlinear partial differential equations (PDEs). It is hard to find non-trivial analytical solutions. Therefore, numerically solving these nonlinear PDEs is an efficient approach. To guarantee convergence, numerical schemes of these equations are required to satisfy the energy dissipation property. Meanwhile, an accurate and efficient approach should be designed to deal with nonlinear terms. In terms of time discretization, there have been several effective methods which can preserve energy dissipation law, including the convex splitting method

[4, 15, 16, 17, 18, 19, 20, 21], stabilized approach[22, 23, 25, 24], invariant energy quadratization (IEQ) method[26, 19] and recently developed scalar auxiliary variable (SAV) approach. By introduing a scale auxiliary variable to the nonlinear part of energy functional, the SAV approach has a modified energy dissipation property for a large class of gradient flows. The convergent and error analysis of semi-discrete SAV scheme has been given by Shen and Xu. The analysis of energy stability and convergence of fully discretized SAV block-centered finite difference method has been established for gradient flows. More studies about the PFC problem can be found in recent literature [16, 43, 17, 44, 45].

In the study to the PFC model, finite difference methods[2, 3, 4, 5, 46] or spectral methods[27, 28, 29] are limited to regular regions, such as two-dimensional square region or three-dimensional cube region. For complex geometries, finite element method (FEM)[30, 55, 42] is a better choice. Furthermore, the FEM can be further combined with adaptive technologies, which are well suitable for the phase behavior of PFC models, such as the formation of ordered structures, phase transition processes, and coarse-grained processes. The adaptive method can effectively decrease the cost of computing and accurately capture the phase interface.

In this work, we will combine SAV time discretization and FEM spatial discretization to solve the gradient flow equation of LB model. Based on the energy dissipation and the SAV scheme, the derivation process of bounds of the solution is shown in detail. For our fully discrete scheme, we demonstrate its energy stability, and carry out error estimate. Applying our method, we can effectively simulate the mesoscale self-assembly of the diblock copolymer system in two-dimensional convex geometries. In addition, we will consider an adaptive FEM for the PFC model. There are many adaptive finite element methods for phase field equation[49, 50, 51, 52, 53]. To reduce computational cost, we first apply an adaptive method which is effective against phase field equation to the PFC model. Numerical results demonstrate that directly using the gradient as the indicator is more efficient than the posterior error estimator does in solving this problem. Since the gradient obtained from the numerical solution may be discontinuous, therefore, a smooth recovered gradient is employed as the adaptive indicator in our adaptive FEM.

The rest of this paper is organized as follows. In Sec. 2, we introduce the LB free energy functional and take the gradient flow as an example to derive its Allen-Cahn dynamic equation. Sec. 3 details our numerical method, which consists of discretization schemes, energy dissipation and error estimate for and dynamical equations. In Sec. 4, numerical experiments are given to illustrate the accuracy and effectiveness of our scheme. Several standard ordered structures on two-dimensional convex regions are also obtained in this section. Sec. 5 gives a simple but efficient adaptive FEM to PFC model. In Sec. 6, we give conclusions and outlooks.

## 2 Physical model

The dimensionless free energy functional of LB model is

 E(ϕ(r))=∫Ω{ξ22[(∇2+1)ϕ(r)]2+α2[ϕ(r)]2−γ3![ϕ(r)]3+14![ϕ(r)]4}dr. (1)

where is a two dimensional bounded domain with Lipschitz boundary , is the density deviation of a kind of monomer from the disordered phase, and are the parameters of the model, is the Laplace operator.

The gradient flow of LB model is

 ϕt=LδEδϕ, (2)

where is a negative operator. For gradient flow , while for gradient flow, . For brevity, in the following, we take the gradient flow as an example to derive the Allen-Cahn dynamical equation. The gradient flow equation can be derived similarly.

Now we will give the boundary conditions to guarantee the energy dissipation property of the Allen-Cahn dynamical equation. Denote that

 N(ϕ):=α2ϕ2−γ3!ϕ3+14!ϕ4, (3)

then

 N′(ϕ)=αϕ−γ2ϕ2+13!ϕ3,   N′′(ϕ)=α−γϕ+12ϕ2, (4)

therefore, the LB model becomes

 E(ϕ(r))=∫Ω{ξ22[(∇2+1)ϕ]2+N(ϕ)}dr. (5)

The free energy take a derivative with respect to time is

 ∂E(ϕ)∂t = ∫Ω{ξ2[(∇2+1)ϕ(∇2+1)ϕt]+N′(ϕ)ϕt}dr = ∫Ω[ξ2(∇2+1)ϕ+N′(ϕ)]ϕtdr+ξ2∫Ω(∇2+1)ϕ∇2ϕtdr = ∫Ω[ξ2(∇2+1)ϕ+N′(ϕ)]ϕtdr+ξ2∫∂Ω(∇2+1)ϕ∇ϕt⋅ndS −ξ2∫Ω∇(∇2+1)ϕ⋅∇ϕtdr = ∫Ω[ξ2(∇2+1)ϕ+N′(ϕ)]ϕtdr+ξ2∫∂Ω(∇2+1)ϕ∇ϕt⋅ndS

here we introduce two Neumann boundary conditions

 ∇ϕ⋅n|∂Ω=0,∇(∇2+1)ϕ⋅n|∂Ω=0, (6)

then

 ∂E(ϕ)∂t=∫Ω[ξ2(∇2+1)2ϕ+N′(ϕ)]ϕtdr. (7)

Define a function space as

 W(Ω):={w∈H2(Ω):∇w⋅n|∂Ω=0,∇(∇2+1)w⋅n|∂Ω=0}. (8)

In the sense of the Gateaux differential for all , we have the following equation

 (δEδϕ,v) = ddθE(ϕ+θv)|θ=0 = limθ→01θ[E(ϕ+θv)−E(ϕ)] = (ξ2(∇2+1)2ϕ+N′(ϕ),v),

where , denote the inner product. According to the variational principle, we have

 δEδϕ=ξ2(∇2+1)2ϕ+N′(ϕ). (9)

From Eqns. (2), (7) and (9), it is easy to verify

 ∂E(ϕ)∂t=−∫Ω[ξ2(∇2+1)2ϕ+N′(ϕ)]2dr≤0. (10)

Therefore, Allen-Cahn equation (2) satisfies energy dissipation with the Neumann boundary conditons (6).

Combining Eqns. (2), (3), (6) and (9), the governing equation can be written as

 ϕt=−ξ2(∇2+1)2ϕ−N′(ϕ), (11a) ∇ϕ⋅n|∂Ω=0,∇(∇2+1)ϕ⋅n|∂Ω=0. (11b)

By introducing a new function , we can write Eqns. (11) as the coupled system

 ϕt=−ξ2(∇2+1)ψ−N′(ϕ), (12a) ψ=(∇2+1)ϕ, (12b) ∇ϕ⋅n|∂Ω=0,∇ψ⋅n|∂Ω=0. (12c)

Remark 2.1 The splitting technique used in Eqns. (12) is valid for convex regions.

## 3 Numerical Methods

The main aim of this section is to propose numerical methods to solve the gradient flows of LB model. For brevity, for the Allen-Cahn equation, we first present the discretization scheme, prove the energy stability, and give the error estimate in detail. Then the corresponding results about the Cahn-Hilliard equation are also given.

### 3.1 SAV discretization

Recently, Shen et al. proposed an efficient time discretization scheme, i.e., the SAV scheme, to a class of gradient flows. Here, we shall apply the idea of the SAV approach to discretize Eqns. (12) in time direction.

Let . Then we introduce the scalar auxiliary variable , where is a constant to ensure , and write Eqns. (12) as

 ϕt=−{ξ2(∇2+1)ψ+u(ϕ)s}, (13a) ψ=(∇2+1)ϕ, (13b) st=12∫Ωu(ϕ)ϕtdr, (13c)

where .

Assume is a fixed time step and is the approximation of at time , we can construct the first-order SAV scheme:

 ϕn+1−ϕnΔt =−{ξ2(∇2+1)ψn+1+u(ϕn)sn+1}, (14a) ψn+1 =(∇2+1)ϕn+1, (14b) sn+1−sn =12(u(ϕn),ϕn+1−ϕn). (14c)

Remark 3.1 The higher-order schemes based on SAV technique can be easily constructed, see Ref. for more details.

### 3.2 FEM discretization

We discretize Eqns. (14) in space using the FEM. Let denote both the trial and test function spaces

 V(Ω):={w∈H1(Ω):∇w⋅n|∂Ω=0}. (15)

The corresponding Galerkin form of Eqns. (14) can be stated as follows: for , find such that:

 (ϕn+1−ϕnΔt,v) =−{ξ2(ψn+1,v)−ξ2(∇ψn+1,∇v)+(u(ϕn),v)sn+1}, (16a) (ψn+1,v) (16b) sn+1−sn =12(u(ϕn),ϕn+1−ϕn). (16c)

Let be a comforming mesh of with , be line segment in 1D or triangle in 2D, and be the linear finite element space over defined by

 Vh:={w∈H1(Ω):w|τ∈P1(τ),∀k∈Th}.

Thus, Eqns. (16) is transformed as follow: find , such that for :

 (ϕn+1h−ϕnhΔt,vh) =−{ξ2(ψn+1h,vh)−ξ2(∇ψn+1h,∇vh)+(u(ϕnh),vh)sn+1h}, (17a) (ψn+1h,vh) =(ϕn+1h,vh)−(∇ϕn+1h,∇vh), (17b) sn+1h−snh =12(u(ϕnh),ϕn+1h−ϕnh), (17c)

### 3.3 Energy stability

We shall prove the energy stability of Eqns. (17). The norm of is denoted by .

Theorem 3.1 If we denote the modified energy

 ~E(ψ,s):=ξ22∥ψ∥2+s2,

then Eqns. (17) is unconditionally energy stable with the modified energy.

Proof: We take in Eqn. (17a) , and find

 (18) +(u(ϕnh),ϕn+1h−ϕnh)sn+1h=0.

According to Eqn. (17b), we have

 (ψnh,vh)=(ϕnh,vh)−(∇ϕnh,∇vh). (19)

Substracting Eqn. (17b) by Eqn. (19), we have

 (20)

Setting to Eqn. (20) and using the identity

 (a−b,2a)=|a|2−|b|2+|a−b|2,

we have

 (ϕn+1h−ϕnh,ψn+1h)−(∇(ϕn+1h−ϕnh),∇ψn+1h) (21) =(ψn+1h−ψnh,ψn+1h) =12(∥ψn+1h∥2−∥ψnh∥2+∥ψn+1h−ψnh∥2).

Multiplying Eqn. (17c) with , we obtain

 (u(ϕnh),ϕn+1h−ϕnh)sn+1h (22) =2(sn+1h−snh)sn+1h =(sn+1h)2−(snh)2+(sn+1h−snh)2.

Substituting Eqn. (21) and Eqn. (22) into Eqn. (18), we have

 1Δt∥ϕn+1h−ϕnh∥2+ξ22∥ψn+1h∥2−ξ22∥ψnh∥2 (23) +ξ22∥ψn+1h−ψnh∥2+(sn+1h)2−(snh)2+(sn+1h−snh)2=0.

Then the discretized energy dissipative property is satisfied, i.e.

 ~E(ψn+1h,sn+1h)−~E(ψnh,snh) (24) =−{1Δt∥ϕn+1h−ϕnh∥2+ξ22∥ψn+1h−ψnh∥2+(sn+1h−snh)2} ≤0.

The energy stability of Eqns. (17) is derived for the modified free energy , not for the original one , owing to the introduction of a function and a scaler variable in the constructing processes of Eqns. (17). The modified free energy plays an important role in the error analysis of the discrete scheme.

### 3.4 Error estimate

We first give some lemmas below.

Lemma 3.1 Assume . Let be solutions of Eqns. (13). There exists a constant depending only on and such that,

 ∥ϕ(t)∥H2≤C. (25)

Proof: It is easy to know that Eqns. (13) are unconditionally energy stability with modified energy (24). There exists a constant depending only on and such that ,

 ξ22∥ψ(tn)∥2+|s(tn)|2≤C. (26)

We take the inner product of Eqn. (13a) with , then for all ,

 (ϕt(tn),ϕ(tn))+ξ2((∇2+1)ψ(tn),ϕ(tn))+(u(ϕ(tn−1))s(tn),ϕ(tn))=0. (27)

Integrating Eqn. (13b) with , we obtain

 (ψ(tn),(∇2+1)ϕ(tn))=∥ψ(tn)∥2. (28)

Combining Eqn. (27) and Eqn. (28), and using Hölder’s inequality and Young’s inequality, we have

 12d∥ϕ(tn)∥2dt+∥ψ(tn)∥2 =−(u(ϕ(tn−1))s(tn),ϕ(tn)) (29) ≤|s(tn)|∥u(ϕ(tn−1))∥∥ϕ(tn)∥.

Thanks to Eqn. (26), then

 12d∥ϕ(tn)∥dt ≤C∥u(ϕ(tn−1))∥. (30)

By Minkowski’s inequality and Sobolve embedding, and choosing the appropriate such that we have

 ∥ϕ(tn)∥ ≤CΔt∥u(ϕ(tn−1))∥+∥ϕ(tn−1)∥ (31) =CΔt∥N′(ϕ(tn−1))√E1(ϕ(tn−1))+D0∥+∥ϕ(tn−1)∥ ≤CΔt∥N′(ϕ(tn−1))∥+∥ϕ(tn−1)∥ =CΔt∥αϕ(tn−1)−γ2(ϕ(tn−1))2+16(ϕ(tn−1))3∥+∥ϕ(tn−1)∥ ≤CΔt(α∥ϕ(tn−1)∥+γ2∥(ϕ(tn−1))2∥+16∥(ϕ(tn−1))3∥)+∥ϕ(tn−1)∥ ≤CΔt(α∥ϕ(tn−1)∥+c1γ2∥ϕ(tn−1)∥2L6+16∥ϕ(tn−1)∥3L6)+∥ϕ(tn−1)∥ ≤CΔt(α∥ϕ(tn−1)∥+c1γ2∥ϕ(tn−1)∥2H1+16∥ϕ(tn−1)∥3H1)+∥ϕ(tn−1)∥.

Next we shall prove using mathematical induction.
When , .
If , , by Minkowski’s inequality and Hölder’s inequality, we have the results as follow:

 ∥∇2ϕ(tk)∥=∥(∇2+1)ϕ(tk)−ϕ(tk)∥≤∥(∇2+1)ϕ(tk)∥+∥ϕ(tk)∥≤C, (32) ∥∇ϕ(tk)∥2=−(∇2ϕ(tk),ϕ(tk))≤∥∇2ϕ(tk)∥∥ϕ(tk)∥≤C.

Note that and , using Eqns. (31), we deduce that .
Thus, . In addition to Eqns. (32), we can deduce Eqn. (33).

For Eqns. (17), we can obtain Lemma 3.2, as its proofs essentially the same as for Lemma 3.1.

Lemma 3.2 Assume . Let be solutions of Eqns. (17). Then for all and , we have

 ∥ϕnh∥H2≤C. (33)

where and is the discrete Laplace operator.

Here, we shall derive error estimates of Eqns. (17). Denote , , , where .

Let be the standard elliptic (Ritz) projection operator, satisfying

 (∇(Rhw(t)−w(t)),∇vh)=0,∀vh∈Vh,

and .

Lemma 3.3[38, 39] If are sufficiently smooth, there exist a positive constant independent of , such that

 ∥ρ(t)∥+h∥∇ρ(t)∥≤Ch2∥w∥H2, (34) ∥ρt(t)∥+h∥∇ρt(t)∥≤Ch2∥wt∥H2. (35)

Theorem 3.2 Let and be solutions of Eqns. (13) and Eqns. (17), respectively. Assume . In addition, we assume that

 ϕ∈C(0,T;H2),ϕt∈L∞(0,T;L2)∩L2(0,T;L2)∩C(0,T;H2),ϕtt∈L2(0,T;L2),

and

 ψ∈C(0,T;H2),ψt∈C(0,T;H2)

For , such that

 ∥enϕ,h∥2+ξ22∥enψ,h∥2+(ens,h)2≤C(K1Δt2+K2h4). (36)

where

 K1 =∥ϕt∥4L2(0,T;L2)+∥ϕtt∥2L2(0,T;L2), K2 =∥ϕ∥2C(0,T;H2)+∥ϕt∥2C(0,T;H2)+∥ψ∥2C(0,T;H2)+∥ψt∥2C(0,T;H2).

Proof: By Lemma 3.1 and Lemma 3.2, we know that

 ∥ϕ(t)∥H2,∥ϕnh∥H2≤C.

Note that . Therefore, we can find a constant such that

 |N(ϕ(t))|,|N′(ϕ(t))|,|N′′(ϕ(t))|≤C, |N(ϕnh)|,|N′(ϕnh)|,|N′′(ϕnh)|≤C. (37)

Subtracting the variational formulation of Eqn. (13a) from Eqn. (17a), we have

 (Dten+1ϕ,h+ξ2en+1ψ,h,vh)−ξ2(∇en+1ψ,h,∇vh) =−(u(ϕnh)sn+1h−u(ϕ(tn))s(tn+1),vh)−(Dtϕ(tn+1)−ϕt(tn+1),vh). (38)

Then

 (Dtθn+1ϕ,h+ξ2θn+1ψ,h,vh)−ξ2(∇θn+1ψ,h,∇vh) =−(u(ϕnh)sn+1h−u(ϕ(tn))s(tn+1),vh)−(Dtϕ(tn+1)−ϕt(tn+1),vh) −(Dtρn+1ϕ,h+ξ2ρn+1ψ,h,vh)−ξ2(∇ρn+1ψ,h,∇vh) =−(u(ϕnh)sn+1h−u(ϕ(tn))s(tn+1),vh)−(Dtϕ(tn+1)−ϕt(tn+1),vh) −(Dtρn+1ϕ,h+ξ2ρn+1ψ,h,vh). (39)

Subtracting the variational formulation of Eqn. (13b) from Eqn. (17b), we get

 (en+1ψ,h,vh)=(en+1ϕ,h,vh)−(∇en+1ϕ,h,∇vh). (40)

Then

 (θn+1ϕ,h,vh)−(∇θn+1ϕ,h,∇vh) =(θn+1ψ,h,vh)+(ρn+1ψ,h,vh)−(ρn+1ϕ,h,vh)+(∇ρn+1ϕ,h,∇vh) =(θn+1ψ,h,vh)+(ρn+1ψ,h−ρn+1ϕ,h,vh). (41)

Taking in Eqn. (3.4) and in Eqn. (3.4), we obtain

 12Dt∥θn+1ϕ,h∥2+∥θn+1ϕ,h−θnϕ,h∥22Δt+ξ2∥θn+1ψ,h∥2 =−(u(ϕnh)sn+1h−u(ϕ(tn))s(tn+1),θn+1ϕ,h) −(Dtϕ(tn+1)−ϕt(tn+1),θn+1ϕ,h) −(Dtρn+1ϕ,h+ξ2ρn+1ψ,h,θn+1ϕ,h) −ξ2(ρn+1ψ,h−ρn+1ϕ,h,θn+1ψ,h) =(u(ϕnh)(sn+1h−s(tn+1)),θn+1ϕ,h) +s(tn+1)(u(ϕnh)−u(ϕ(tn)),θn+1ϕ,h) −(Dtϕ(tn+1)−ϕt(tn+1),θn+1ϕ,h) −(Dtρn+1ϕ,h+ξ2ρn+1ψ,h,θn+1ϕ,h) −ξ2(ρn+1ψ,h−ρn+1ϕ,h,θn+1ψ,h) :=I1+I2+⋯+I5. (42)

Using Hölder’s inequality, we get

 I1 =(u(ϕnh)en+1s,h,θn+1ϕ,h)≤C(en+1s,h)2+∥θn+1ϕ,h∥2. (43)

Due to

 u(ϕnh)−u(ϕ(tn)) =N′(ϕnh)√E1(ϕnh)+D0−N′(ϕ(tn))√E1(ϕ(tn
{"ft_lang_label":"__label__en","ft_lang_prob":0.88128704,"math_prob":0.99839294,"size":11134,"snap":"2021-31-2021-39","text_gpt3_token_len":2776,"char_repetition_ratio":0.123809524,"word_repetition_ratio":0.02111111,"special_character_ratio":0.25139213,"punctuation_ratio":0.1471541,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.9985471,"pos_list":[0,1,2],"im_url_duplicate_count":[null,null,null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2021-09-25T11:55:51Z\",\"WARC-Record-ID\":\"<urn:uuid:1bcfff3a-2b8f-45be-998b-7be1a0435a8d>\",\"Content-Length\":\"1048897\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:88ad9113-2471-4e1f-8d63-4b187bf8f4f8>\",\"WARC-Concurrent-To\":\"<urn:uuid:8e250b41-79bd-4c90-9758-8341e545c633>\",\"WARC-IP-Address\":\"44.238.89.211\",\"WARC-Target-URI\":\"https://deepai.org/publication/error-analysis-of-sav-finite-element-method-to-phase-field-crystal-model\",\"WARC-Payload-Digest\":\"sha1:UEDEEGTU2TMIOPPLKKV3ELC6MJIGA6F7\",\"WARC-Block-Digest\":\"sha1:FHMSZK57HF6QMMN3SPT6BQUK5TQL6MXF\",\"WARC-Truncated\":\"length\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2021/CC-MAIN-2021-39/CC-MAIN-2021-39_segments_1631780057622.15_warc_CC-MAIN-20210925112158-20210925142158-00543.warc.gz\"}"}