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

by   Liupeng Wang, et al.

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.


page 21

page 25


Mixed finite element method for a second order Dirichlet boundary control problem

The main aim of this article is to analyze mixed finite element method f...

An adaptive edge-based smoothed finite element method (ES-FEM) for phase-field modeling of fractures at large deformations

This work presents the Griffith-type phase-field formation at large defo...

A vertex scheme for two-phase flow in heterogeneous media

This paper presents the numerical solution of immiscible two-phase flows...

Phase field predictions of microscopic fracture and R-curve behaviour of fibre-reinforced composites

We present a computational framework to explore the effect of microstruc...

Scalar auxiliary variable finite element scheme for the parabolic-parabolic Keller-Segel model

We describe and analyze a finite element numerical scheme for the parabo...

Thermodynamically Consistent Diffuse Interface Model for Cell Adhesion and Aggregation

A thermodynamically consistent phase-field model is introduced for simul...

Microstructure-informed reduced modes synthesized with Wang tiles and the Generalized Finite Element Method

A recently introduced representation by a set of Wang tiles – a generali...

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[4]. 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

[14]. Compared with the typical Swift-Hohenberg (SH) model with double-well bulk energy [8], 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[29]. 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[41]. The analysis of energy stability and convergence of fully discretized SAV block-centered finite difference method has been established for gradient flows[40]. 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[13]


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


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




therefore, the LB model becomes


The free energy take a derivative with respect to time is

here we introduce two Neumann boundary conditions




Define a function space as


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

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


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


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


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


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

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[29]. 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


where .

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


Remark 3.1 The higher-order schemes based on SAV technique can be easily constructed, see Ref.[29] 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


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


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

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


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

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

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


According to Eqn. (17b), we have


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


Setting to Eqn. (20) and using the identity

we have


Multiplying Eqn. (17c) with , we obtain


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


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


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,


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 ,


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


Integrating Eqn. (13b) with , we obtain


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


Thanks to Eqn. (26), then


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


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


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


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

and .

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


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


For , such that



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

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


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




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




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


Using Hölder’s inequality, we get


Due to