Stability and Error estimates of the SAV Fourier-spectral method for the Phase Field Crystal Equation
SSTABILITY AND ERROR ESTIMATES OF THE SAV FOURIER-SPECTRAL METHODFOR THE PHASE FIELD CRYSTAL EQUATION ∗ XIAOLI LI † AND
JIE SHEN ‡ Abstract.
We consider fully discrete schemes based on the scalar auxiliary variable (SAV) approach and stabilized SAVapproach in time and the Fourier-spectral method in space for the phase field crystal (PFC) equation. Unconditionally energystability is established for both first- and second-order fully discrete schemes. In addition to the stability, we also provide arigorous error estimate which shows that our second-order in time with Fourier-spectral method in space converges with order O (∆ t + N − m ), where ∆ t , N and m are time step size, number of Fourier modes in each direction, and regularity index inspace, respectively. We also present numerical experiments to verify our theoretical results and demonstrate the robustness andaccuracy of the schemes. Key words.
Phase field crystal, Fourier-spectral method, scalar auxiliary variable (SAV), energy stability, error estimates
AMS subject classifications.
1. Introduction.
The phase field crystal equation, which was developed in [2, 3], has been frequentlyused in the study of the microstructural evolution of crystal growth on atomic length and diffusive timescales. It is well known that the crystal growth is the major stage in crystallization, which is an importantstep in the purification of solid compounds. The PFC equation is a sixth-order nonlinear parabolic equationand the phase field variable is introduced to describe the phase transition from the liquid phase to the crystalphase. The PFC equation has been employed to simulate a number of physical phenomena, including crystalgrowth in a supercooled liquid, dendritic and eutectic solidification, epitaxial growth, material hardness andreconstructive phase transitions.It is challenging to develop efficient and accurate numerical schemes for the PFC equation due to thesix-order spacial derivative and its nonlinearity. Gomez and Nogueira [4] proposed a numerical algorithmfor the phase field crystal equation which is second-order time-accurate and unconditionally stable. Localdiscontinuous Galerkin method has been developed by Guo and Xu [5] for the the PFC equation, which isbased on the first order and second order convex splitting principle. Li and Kim [8] studied an efficient andstable compact fourth-order finite difference scheme for the phase field crystal equation. It is worth notingthat all these schemes are nonlinear, so their implementations are relatively complex and costly comparedto linear schemes. To obtain linear schemes for this model, the main difficulty is how to discretize thequartic potential. Yang and Han [17] constructed linearly unconditionally energy stable schemes for thePFC equation by adopting the ”Invariant Energy Quadratization” (IEQ) approach. They established theunconditionally energy stability. We should point out that although there are many works on the numericalsimulation of the PFC model, very few are with convergence analysis and error estimates. Note that Wise,Wang and Lowengrub [15] proved the error estimates for the nonlinear first order finite difference methodbased on the convex splitting method.The main goals of this paper are to construct linear and unconditionally energy stable schemes basedon the recently proposed scalar auxiliary variable (SAV) approach [12], and provide rigorous error analysisfor them. The work presented in this paper for the PFC model is unique in the following aspects. First, weconstruct two linear, unconditional energy stable schemes for the PFC model based on the stabilized scalarauxiliary variable (S-SAV) approach in time and Fourier-spectral method in space, where extra stabilizedterms are added, compared with [17, 15], while keeping the required accuracy. Second, we carry out rigorouserror analysis, which is made possible by the uniform bound of the discrete solutions that we derive thanksto the unconditional energy stability. We believe that this is the first such result for any fully discrete linearschemes for the PFC model. ∗ The work of X. Li is supported by the Postdoctoral Science Foundation of China under grant numbers BX20190187 and2019M650152. The work of J. Shen is supported in part by NSF grants DMS-1620262, DMS-1720442 and AFOSR grantFA9550-16-1-0102. † School of Mathematical Sciences and Fujian Provincial Key Laboratory on Mathematical Modeling and High PerformanceScientific Computing, Xiamen University, Xiamen, Fujian, 361005, China. Email: [email protected]. ‡ Corresponding Author. Department of Mathematics, Purdue University, West Lafayette, IN 47907, USA. Email:[email protected]. 1
This manuscript is for review purposes only. a r X i v : . [ m a t h . A P ] J u l XIAOLI LI AND JIE SHEN
The paper is organized as follows. In Section 2 we introduce the governing system and some preliminaries.In Section 3 we present the fully discrete schemes using the S-SAV approach in time and Fourier-spectralmethod in space, and we prove their unconditional energy stability In Section 4. In Section 5 we providerigorous error estimate for our fully discrete schemes. In Section 6 we present some numerical experimentsto verify the accuracy of the proposed numerical schemes. Some concluding remarks are given in the lastsection.
2. Governing system and some preliminaries.
Consider the free energy of Swift-Hohenberg type(cf. [14, 17, 6]) E ( φ ) = (cid:90) Ω (cid:18) φ (∆ + β ) φ + 14 φ − (cid:15) φ (cid:19) d x , (1)where the phase field variable φ is the atomic density field and β and (cid:15) are two positive constants such that (cid:15) < β and (cid:15) (cid:28)
1. We assume Ω = (0 , L x ) × (0 , L y ) and φ is Ω-periodic. Then the PFC equation, whichdescribes the phenomena of crystal growth on the atomic length and diffusive time scales, can be modeledby: ∂φ∂t = M ∆ µ, x ∈ Ω , t > ,µ = (∆ + β ) φ + φ − (cid:15)φ, x ∈ Ω , t > ,φ ( x ,
0) = φ ( x ) , (2)where M is the mobility function, µ = δE ( φ ) δφ is the chemical potential. We impose the periodic boundaryconditions, then it can be easily obtained that the PFC equation (2) is mass-conservative in the sense that ddt (cid:82) Ω φd x = 0. Besides, we can derive that the system satisfies the following energy law by using integrationby parts: ddt E ( φ ) = −(cid:107)√ M ∇ µ (cid:107) ≤ . (3)Now we give some notations that will be used later. For each s ≥
0, Let ( · , · ) H s and (cid:107) · (cid:107) H s be the H s (Ω) inner product and norm, respectively. Note that H (Ω) = L (Ω). In particular, we use ( · , · ) and (cid:107) · (cid:107) to denote the L inner product. We define Sobolev spaces L per (Ω) = { v ∈ L (Ω) | v is periodic on Ω } and H sper (Ω) = { v ∈ H s (Ω) | v is periodic on Ω } . Besides, let N T > J = (0 , T ] in thispaper. Set ∆ t = T /N T , t n = n ∆ t, for n ≤ N T , where T is the final time.Throughout the paper we use C , with or without subscript, to denote a positive constant, which couldhave different values at different appearances.
3. Fully discrete schemes by the stabilized SAV Fourier-spectral method.
In this section, wefirst construct semi-discrete stabilized SAV schemes with the first- and second-order accuracy for the PFCequation, followed by the fully discretization with Fourier-spectral method in space.
To construct SAV schemes with a linear stabilization, we recastthe second equation in the PFC model (2) by µ = (∆ + β ) φ + λφ + F (cid:48) ( φ ) , (4)where λ is a positive constant and F ( φ ) = φ − (cid:15) + λ φ . In the SAV approach, a scalar variable r ( t ) = (cid:112) E ( φ )is introduced, where E ( φ ) = (cid:82) Ω F ( φ ) d x + C and C ≥ E ( φ ) >
0. Then thePFC model can be transformed into the following system: ∂φ∂t = M ∆ µ, (5a) µ = (∆ + β ) φ + λφ + r ( t ) (cid:112) E ( φ ) F (cid:48) ( φ ) , (5b) r t = 12 (cid:112) E ( φ ) (cid:90) Ω F (cid:48) ( φ ) φ t d x . (5c) This manuscript is for review purposes only. he SAV and Fourier-spectral method for the PFC Equation Scheme I (first-order accuracy):
Assuming φ n and R n are known, we update φ n +1 and R n +1 bysolving φ n +1 − φ n = M ∆ t ∆ µ n +1 , (6a) µ n +1 = (∆ + β ) φ n +1 + λφ n +1 − S ∆( φ n +1 − φ n ) + R n +1 (cid:112) E ( φ n ) F (cid:48) ( φ n ) , (6b) R n +1 − R n = 12 (cid:112) E ( φ n ) (cid:90) Ω F (cid:48) ( φ n )( φ n +1 − φ n ) d x , (6c)where S >
Scheme II (second-order accuracy):
Assuming φ n , R n and φ n − , R n − are known, then we update φ n +1 and R n +1 by solving φ n +1 − φ n = M ∆ t ∆ µ n +1 / , (7a) µ n +1 / = (∆ + β ) φ n +1 / + λφ n +1 / − S ∆( φ n +1 − φ n + φ n − )+ R n +1 / (cid:113) E ( ˜ φ n +1 / ) F (cid:48) ( ˜ φ n +1 / ) , (7b) R n +1 − R n = 12 (cid:113) E ( ˜ φ n +1 / )) (cid:90) Ω F (cid:48) ( ˜ φ n +1 / ))( φ n +1 − φ n ) d x , (7c)where ˜ φ n +1 / = (3 φ n − φ n − ) / S > n = 0, we can computer˜ φ / by the first order scheme. We first describe theFourier-spectral framework. We partition the domain Ω = (0 , L x ) × (0 , L y ) uniformly with size h x = L x /N x , h y = L y /N y where N x and N y are positive integers. The Fourier approximation space is S N = span { e iξ k x e iη l y : − N x ≤ k ≤ N x − , − N y ≤ l ≤ N y − } , where i = √− ξ k = 2 πk/L x and η l = 2 πl/L y . Then any function u ( x, y ) ∈ L (Ω) can be approximated by u ( x, y ) ≈ u N ( x, y ) = Nx − (cid:88) k = − Nx Ny − (cid:88) l = − Ny ˆ u k,l e iξ k x e iη l y , (8)where the Fourier coefficients are denoted asˆ u k,l = < u, e iξ k x e iη l y > = 1 | Ω | (cid:90) Ω ue − i ( ξ k x + η l y ) d x . In what follows, we take N = N x = N y for simplicity. Then the fully discrete schemes with Fourierspectral method in space based on the mixed formulation can be constructed as follows: Scheme I (first-order accuracy):
Assuming φ nN and R n are known, then we update φ n +1 N and R n +1 by solving ( φ n +1 N − φ nN , q ) + M ∆ t ( ∇ µ n +1 N , ∇ q ) = 0 , ∀ q ∈ S N , (9a)( µ n +1 N , Ψ) = (cid:0) (∆ + β ) φ n +1 N , (∆ + β )Ψ (cid:1) + λ ( φ n +1 N , Ψ)+ S (cid:0) ∇ ( φ n +1 N − φ nN ) , ∇ Ψ (cid:1) + R n +1 (cid:112) E ( φ nN ) ( F (cid:48) ( φ nN ) , Ψ) , ∀ Ψ ∈ S N , (9b) R n +1 − R n = 12 (cid:112) E ( φ nN ) ( F (cid:48) ( φ nN ) , φ n +1 N − φ nN ) , (9c) This manuscript is for review purposes only.
XIAOLI LI AND JIE SHEN where φ N is the L -orthogonal projection of φ , which will be defined late. Scheme II (second-order accuracy):
Assuming φ nN , R n and φ n − N , R n − are known, then we update φ n +1 N and R n +1 by solving( φ n +1 N − φ nN , q ) + M ∆ t ( ∇ µ n +1 / N , ∇ q ) = 0 , ∀ q ∈ S N , (10a)( µ n +1 / N , Ψ) = (cid:16) (∆ + β ) φ n +1 / N , (∆ + β )Ψ (cid:17) + λ ( φ n +1 / N , Ψ) (10b)+ S (cid:0) ∇ ( φ n +1 N − φ nN + φ n − N ) , ∇ Ψ (cid:1) + R n +1 / (cid:113) E ( ˜ φ n +1 / N ) ( F (cid:48) ( ˜ φ n +1 / N ) , Ψ) , ∀ Ψ ∈ S N ,R n +1 − R n = 12 (cid:113) E ( ˜ φ n +1 / N ) ( F (cid:48) ( ˜ φ n +1 / N ) , φ n +1 N − φ nN ) . (10c)
4. Unconditional energy stability.
In this section, we prove the unconditional energy stability forthe first- and second-order fully discrete schemes. Same results can be established for their semi-discreteversions using a similar approach so we omit the detail here.We define the discrete energy as E ( φ nN , R n ) = 12 (cid:107) (∆ + β ) φ nN (cid:107) + λ (cid:107) φ nN (cid:107) + R − C . (11) We have the following results for the first-order scheme.
Theorem Let S ≥ . The first-order fully discrete scheme (9) is unconditionally energy stable in thesense that the following discrete energy law holds for any ∆ t : E ( φ n +1 N , R n +1 ) − E ( φ nN , R n ) + S (cid:107)∇ ( φ n +1 N − φ nN ) (cid:107) ≤ − M ∆ t (cid:107)∇ µ n +1 N (cid:107) . (12) Proof.
Taking q = µ n +1 N and Ψ = φ n +1 N − φ nN in (9) and multiplying (9c) with 2 R n +1 lead to( φ n +1 N − φ nN , µ n +1 N ) + M ∆ t (cid:107)∇ µ n +1 N (cid:107) = 0 . (13)( µ n +1 N , φ n +1 N − φ nN ) = (cid:0) (∆ + β ) φ n +1 N , (∆ + β )( φ n +1 N − φ nN ) (cid:1) + λ ( φ n +1 N , φ n +1 N − φ nN ) + S (cid:107)∇ ( φ n +1 N − φ nN ) (cid:107) + R n +1 (cid:112) E ( φ nN ) ( F (cid:48) ( φ nN ) , ( φ n +1 N − φ nN )) . (14)( R n +1 − R n , R n +1 ) = R n +1 (cid:112) E ( φ nN ) ( F (cid:48) ( φ nN ) , φ n +1 N − φ nN ) . (15)Noting the identity 2( a − b, a ) = a − b + ( a − b ) and combing the above equations, we can obtain12 (cid:107) (∆ + β ) φ n +1 N (cid:107) + λ (cid:107) φ n +1 N (cid:107) + ( R n +1 ) + S (cid:107)∇ ( φ n +1 N − φ nN ) (cid:107) − (cid:18) (cid:107) (∆ + β ) φ nN (cid:107) + λ (cid:107) φ nN (cid:107) + ( R n ) (cid:19) ≤ − M ∆ t (cid:107)∇ µ n +1 N (cid:107) , (16)which implies the result (12) after we drop some positive terms. This manuscript is for review purposes only. he SAV and Fourier-spectral method for the PFC Equation We consider now the second-order scheme.
Theorem Let S ≥ . The second-order fully discrete scheme (10) is unconditionally energy stable inthe sense that the following discrete energy law holds for any ∆ t : ˜ E ( φ n +1 N , R n +1 ) − ˜ E ( φ nN , R n ) = − M ∆ t (cid:107)∇ µ n +1 / N (cid:107) , (17) where ˜ E ( φ n +1 N , R n +1 ) = E ( φ n +1 N , R n +1 ) + S (cid:107)∇ ( φ n +1 N − φ nN ) (cid:107) Proof.
Taking q = µ n +1 / N and Ψ = φ n +1 N − φ nN in (10) and multiplying (10c) with 2 R n +1 / lead to( φ n +1 N − φ nN , µ n +1 / N ) + M ∆ t (cid:107)∇ µ n +1 / N (cid:107) = 0 . (18)( µ n +1 / N , φ n +1 N − φ nN ) = (cid:16) (∆ + β ) φ n +1 / N , (∆ + β )( φ n +1 N − φ nN ) (cid:17) + λ ( φ n +1 / N , φ n +1 N − φ nN ) + S (cid:0) ∇ ( φ n +1 N − φ nN + φ n − N ) , ∇ ( φ n +1 N − φ nN ) (cid:1) + R n +1 / (cid:113) E ( ˜ φ n +1 / N ) ( F (cid:48) ( ˜ φ n +1 / N ) , ( φ n +1 N − φ nN )) . (19)( R n +1 − R n , R n +1 / ) = R n +1 / (cid:113) E ( ˜ φ n +1 / N ) ( F (cid:48) ( ˜ φ n +1 / N ) , φ n +1 N − φ nN ) . (20)Then we can obtain the following equation by combing the above equations:12 (cid:107) (∆ + β ) φ n +1 N (cid:107) + λ (cid:107) φ n +1 N (cid:107) + ( R n +1 ) + S (cid:107)∇ ( φ n +1 N − φ nN ) (cid:107) − (cid:18) (cid:107) (∆ + β ) φ nN (cid:107) + λ (cid:107) φ nN (cid:107) + ( R n ) + S (cid:107)∇ ( φ nN − φ n − N ) (cid:107) (cid:19) = − M ∆ t (cid:107)∇ µ n +1 N (cid:107) , (21)which leads to the desired result. Remark We observe that both schemes are unconditionally energy stable for all S ≥ , which impliesthat even without the stabilization term (i.e., S = 0 ), both schemes are also unconditionally energy stable.However, as we shall demonstrate through numerical results later, the stabilization terms are essential toobtain accurate results without using exceedingly small time steps.
5. Error estimates.
In this section, we provide rigorous error estimates for the second-order fullydiscrete scheme (10). Since the proofs for the first-order scheme (9) are essentially the same as for thesecond order scheme, we skip it for brevity.Define the L -orthogonal projection operator Π N : L (Ω) → S N by(Π N u − u, Ψ) = 0 , ∀ Ψ ∈ S N , u ∈ L (Ω) , (22)The following results hold (cf. [10, 9, 1]): For any 0 ≤ µ ≤ m , there exists a constant C such that (cid:107) Π N u − u (cid:107) µ ≤ C (cid:107) u (cid:107) m N µ − m , ∀ u ∈ H mper (Ω) , (23)where for m ≥ k = 0 , · · · , m − H mper (Ω) = { u ∈ H m (Ω) : u ( k ) (0 , · ) = u ( k ) ( L x , · ) , u ( k ) ( · ,
0) = u ( k ) ( · , L y ) } . (24)Moreover, the operator Π N commutes with the derivation on H per (Ω), i.e.Π N ∇ u = ∇ Π N u, ∀ u ∈ H per (Ω) . (25)We first demonstrate that energy stability leads to the H boundedness of the discrete solutions. This manuscript is for review purposes only.
XIAOLI LI AND JIE SHEN
Lemma Let φ nN be the solution of (10) . We have (cid:107) φ nN (cid:107) H ≤ C for all n and N .Proof. Recalling the energy stability (17), we have (cid:107) (∆ + β ) φ nN (cid:107) ≤ C, (cid:107) φ nN (cid:107) ≤ C. (26)Then we have (cid:107) ∆ φ nN (cid:107) ≤ C . Using integration by parts and Cauchy inequality yields (cid:107)∇ φ nN (cid:107) = − (∆ φ nN , φ nN ) ≤ ζ (cid:107) ∆ φ nN (cid:107) + ζ (cid:107) φ nN (cid:107) ≤ C, (27)which implies the desired result.For simplicity, we set e nφ = φ nN − Π N φ n + Π N φ n − φ n = ¯ e nφ + ˇ e nφ ,e nµ = µ nN − Π N µ n + Π N µ n − µ n = ¯ e nµ + ˇ e nµ ,e nr = R n − r n . Theorem Let S ≥ . Suppose that φ ∈ L ∞ (0 , T ; H per (Ω)) (cid:84) L ∞ (0 , T ; H m +1 per (Ω)) , r ∈ H (0 , T ) , ∂ φ∂t ∈ L (0 , T ; H per (Ω)) and ∂ φ∂t ∈ L (0 , T ; H − (Ω)) . If S > , we assume ∂ φ∂t ∈ L (0 , T ; H per (Ω)) additionally. Then for the second-order fully discrete scheme (10) , we have (cid:107) φ kN − φ k (cid:107) + ( R k − r k ) ≤ C ∆ t (cid:90) t k ( (cid:107) ( − ∆) − / ∂ φ∂t (cid:107) + (cid:107) ∂ φ∂t (cid:107) H + (cid:107) ∂ φ∂t (cid:107) H + | d rdt | ) ds + C (cid:107) φ (cid:107) L ∞ ( J ; H m +1 per (Ω)) N − m , (28) where k ≤ N T and the constant C is independent on N and ∆ t .Proof. Subtracting (5) from (10) at t n +1 / , we find(¯ e n +1 φ − ¯ e nφ , q ) + M ∆ t ( ∇ ¯ e n +1 / µ , ∇ q ) = ( Q n +1 / , q ) , (29)(¯ e n +1 / µ , Ψ) = (cid:16) (∆ + β )¯ e n +1 / φ , (∆ + β )Ψ (cid:17) + λ (¯ e n +1 / φ , Ψ)+ R n +1 / (cid:113) E ( ˜ φ n +1 / N ) ( F (cid:48) ( ˜ φ n +1 / N ) , Ψ) − r n +1 / (cid:112) E ( φ n +1 / ) ( F (cid:48) ( φ n +1 / ) , Ψ)+ S (cid:16) ∇ (¯ e n +1 φ − e nφ + ¯ e n − φ ) , ∇ Ψ (cid:17) + S (cid:0) ∇ ( φ n +1 − φ n + φ n − ) , ∇ Ψ (cid:1) , (30) e n +1 r − e nr = 12 (cid:113) E ( ˜ φ n +1 / N ) ( F (cid:48) ( ˜ φ n +1 / N ) , φ n +1 N − φ nN ) − (cid:112) E ( φ n +1 / ) (cid:90) Ω F (cid:48) ( φ n +1 / ) φ n +1 / t ∆ td x + Q n +1 / , (31)where the truncation errors are give by the Taylor expansion: Q n +1 / =∆ t ∂φ n +1 / ∂t − ( φ n +1 − φ n )= 12 (cid:90) t n +1 / t n ( t n − s ) ∂ φ∂t ( s ) ds + 12 (cid:90) t n +1 / t n +1 ( t n +1 − s ) ∂ φ∂t ( s ) ds. (32) Q n +1 / =∆ tr n +1 / t − ( r n +1 − r n )= 12 (cid:90) t n +1 / t n ( t n − s ) d rdt ( s ) ds + 12 (cid:90) t n +1 / t n +1 ( t n +1 − s ) d rdt ( s ) ds. (33) This manuscript is for review purposes only. he SAV and Fourier-spectral method for the PFC Equation q = ¯ e n +1 / µ and Ψ = ¯ e n +1 φ − ¯ e nφ in (29) and (30) respectively, and multiplying (31) with 2 e n +1 / r yield(¯ e n +1 φ − ¯ e nφ , ¯ e n +1 / µ ) + M ∆ t (cid:107)∇ ¯ e n +1 / µ (cid:107) = ( Q n +1 / , ¯ e n +1 / µ ) . (34)(¯ e n +1 / µ , ¯ e n +1 φ − ¯ e nφ ) = 12 ( (cid:107) (∆ + β )¯ e n +1 φ (cid:107) − (cid:107) (∆ + β )¯ e nφ (cid:107) ) + λ (cid:107) ¯ e n +1 φ (cid:107) − (cid:107) ¯ e nφ (cid:107) )+ R n +1 / (cid:113) E ( ˜ φ n +1 / N ) ( F (cid:48) ( ˜ φ n +1 / N ) , ¯ e n +1 φ − ¯ e nφ ) − r n +1 / (cid:112) E ( φ n +1 / ) ( F (cid:48) ( φ n +1 / ) , ¯ e n +1 φ − ¯ e nφ )+ S (cid:16) ∇ (¯ e n +1 φ − e nφ + ¯ e n − φ ) , ∇ (¯ e n +1 φ − ¯ e nφ ) (cid:17) + S (cid:16) ∇ ( φ n +1 − φ n + φ n − ) , ∇ (¯ e n +1 φ − ¯ e nφ ) (cid:17) . (35)( e n +1 r ) − ( e nr ) = e n +1 / r (cid:113) E ( ˜ φ n +1 / N ) ( F (cid:48) ( ˜ φ n +1 / N ) , φ n +1 N − φ nN ) − e n +1 / r (cid:112) E ( φ n +1 / ) (cid:90) Ω F (cid:48) ( φ n +1 / ) φ n +1 / t ∆ td x + 2( Q n +1 / , e n +1 / r ) . (36)The term on the right-hand side of (34) can be estimated by( Q n +1 / , ¯ e n +1 / µ ) ≤ M ∆ t (cid:107)∇ ¯ e n +1 / µ (cid:107) + C ∆ t (cid:107) ( − ∆) − / Q n +1 / (cid:107) ≤ M ∆ t (cid:107)∇ ¯ e n +1 / µ (cid:107) + C ∆ t (cid:90) t n +1 t n (cid:107) ( − ∆) − / ∂ φ∂t ( s ) (cid:107) ds, (37)where taking notice of ( Q n +1 / ,
1) = 0, the operator ( − ∆) − / , which is the power of − ∆, can be welldefined by the spectral theory of self-adjoint operators.The third and fourth terms on the right-hand side of (35) can be transformed into R n +1 / (cid:113) E ( ˜ φ n +1 / N ) ( F (cid:48) ( ˜ φ n +1 / N ) , ¯ e n +1 φ − ¯ e nφ ) − r n +1 / (cid:112) E ( φ n +1 / ) ( F (cid:48) ( φ n +1 / ) , ¯ e n +1 φ − ¯ e nφ )= e n +1 / r (cid:113) E ( ˜ φ n +1 / N ) ( F (cid:48) ( ˜ φ n +1 / N ) , ¯ e n +1 φ − ¯ e nφ ) + r n +1 / ( F (cid:48) ( ˜ φ n +1 / N ) (cid:113) E ( ˜ φ n +1 / N ) , ¯ e n +1 φ − ¯ e nφ )) − r n +1 / ( F (cid:48) ( φ n +1 / ) (cid:112) E ( φ n +1 / ) , ¯ e n +1 φ − ¯ e nφ ) . (38)Assuming F ( φ ) ∈ C ( R ) and noting (23) and (25), the last two terms on the right-hand side of (38) can This manuscript is for review purposes only.
XIAOLI LI AND JIE SHEN be controlled, similar to the estimates in [11], by r n +1 / ( F (cid:48) ( ˜ φ n +1 / N ) (cid:113) E ( ˜ φ n +1 / N ) , ¯ e n +1 φ − ¯ e nφ )) − r n +1 / ( F (cid:48) ( φ n +1 / ) (cid:112) E ( φ n +1 / ) , ¯ e n +1 φ − ¯ e nφ )= r n +1 / M ∆ t ( F (cid:48) ( ˜ φ n +1 / N ) (cid:113) E ( ˜ φ n +1 / N ) − F (cid:48) ( φ n +1 / ) (cid:112) E ( φ n +1 / ) , ∆¯ e n +1 / µ )+ r n +1 / ( F (cid:48) ( ˜ φ n +1 / N ) (cid:113) E ( ˜ φ n +1 / N ) − F (cid:48) ( φ n +1 / ) (cid:112) E ( φ n +1 / ) , Q n +1 / ) ≤ M ∆ t (cid:107)∇ ¯ e n +1 / µ (cid:107) + C ∆ t (cid:107) ∇ F (cid:48) ( ˜ φ n +1 / N ) (cid:113) E ( ˜ φ n +1 / N ) − ∇ F (cid:48) ( φ n +1 / ) (cid:112) E ( φ n +1 / ) (cid:107) + C ∆ t (cid:107) ( − ∆) − / Q n +1 / (cid:107) ≤ M ∆ t (cid:107)∇ ¯ e n +1 / µ (cid:107) + C ∆ t ( (cid:107) ¯ e nφ (cid:107) + (cid:107) ¯ e n − φ (cid:107) + (cid:107)∇ ¯ e nφ (cid:107) + (cid:107)∇ ¯ e n − φ (cid:107) )+ C ∆ t (cid:90) t n +1 t n (cid:107) ( − ∆) − / ∂ φ∂t ( s ) (cid:107) ds + C ∆ t (cid:90) t n +1 t n (cid:107) ∂ φ∂t ( s ) (cid:107) H ds + C (cid:107) ˜ φ n +1 / (cid:107) m +1 N − m ∆ t, (39)where the last inequality holds by the fact that one can find a constant C such that | F (cid:48) ( φ nN ) | ≤ C , | F (cid:48)(cid:48) ( φ nN ) | ≤ C by using Lemma 4 and the Sobolev embedding theorem H ⊆ L ∞ .We also have the following using the integration by parts and Cauchy-Schwartz inequality: (cid:107)∇ ¯ e nφ (cid:107) = − (¯ e nφ , ∆¯ e nφ ) ≤ C ( (cid:107) ¯ e nφ (cid:107) + (cid:107) ∆¯ e nφ (cid:107) ) ≤ C ( (cid:107) ¯ e nφ (cid:107) + (cid:107) (∆ + β )¯ e nφ (cid:107) ) . (40)Similar to the estimate in (39), the last term on the right-hand side of (35) can be directly controlledby using integration by parts and Cauchy-Schwarz inequality: S (cid:16) ∇ ( φ n +1 − φ n + φ n − ) , ∇ (¯ e n +1 φ − ¯ e nφ ) (cid:17) = − S (cid:16) ∆( φ n +1 − φ n + φ n − ) , M ∆ t ∆¯ e n +1 / µ + Q n +1 / (cid:17) ≤ M ∆ t (cid:107)∇ ¯ e n +1 / µ (cid:107) + C ∆ t (cid:90) t n +1 t n (cid:107) ∂ φ∂t ( s ) (cid:107) H ds + C ∆ t (cid:90) t n +1 t n (cid:107) ( − ∆) − / ∂ φ∂t ( s ) (cid:107) ds. (41)The first two terms on the right-hand side of (36) can be recast as e n +1 / r (cid:113) E ( ˜ φ n +1 / N ) ( F (cid:48) ( ˜ φ n +1 / N ) , φ n +1 N − φ nN ) − e n +1 / r (cid:112) E ( φ n +1 / ) (cid:90) Ω F (cid:48) ( φ n +1 / ) φ n +1 / t ∆ td x = e n +1 / r (cid:113) E ( ˜ φ n +1 / N ) ( F (cid:48) ( ˜ φ n +1 / N ) , ¯ e n +1 φ − ¯ e nφ ) + e n +1 / r (cid:113) E ( ˜ φ n +1 / N ) ( F (cid:48) ( ˜ φ n +1 / N ) , ˘ e n +1 φ − ˘ e nφ )+ e n +1 / r ( F (cid:48) ( ˜ φ n +1 / N ) (cid:113) E ( ˜ φ n +1 / N ) − F (cid:48) ( φ n +1 / ) (cid:112) E ( φ n +1 / ) , φ n +1 − φ n ) − e n +1 / r ( F (cid:48) ( φ n +1 / ) (cid:112) E ( φ n +1 / ) , Q n +1 / ) . (42) This manuscript is for review purposes only. he SAV and Fourier-spectral method for the PFC Equation e n +1 / r ( F (cid:48) ( ˜ φ n +1 / N ) (cid:113) E ( ˜ φ n +1 / N ) − F (cid:48) ( φ n +1 / ) (cid:112) E ( φ n +1 / ) , φ n +1 − φ n ) ≤ C ∆ t (cid:107) φ t (cid:107) L ∞ ( J ; H − ) (cid:0) ( e n +1 r ) + ( e nr ) + (cid:107) ¯ e nφ (cid:107) + (cid:107) ¯ e n − φ (cid:107) + (cid:107)∇ ¯ e nφ (cid:107) + (cid:107)∇ ¯ e n − φ (cid:107) )+ C ∆ t (cid:90) t n +1 t n (cid:107) ∂ φ∂t ( s ) (cid:107) H ds + C (cid:107) ˜ φ n +1 / (cid:107) H m +1 N − m ∆ t. (43) − e n +1 / r ( F (cid:48) ( φ n +1 / ) (cid:112) E ( φ n +1 / ) , Q n +1 / ) ≤ C ∆ t ( e n +1 / r ) + C ∆ t (cid:90) t n +1 t n (cid:107) ( − ∆) − / ∂ φ∂t ( s ) (cid:107) ds. (44)The last term on the right-hand side of (36) can be controlled by2( Q n +1 / , e n +1 / r ) ≤ C ∆ t ( e n +1 / r ) + C ∆ t (cid:90) t n +1 t n | d rdt ( s ) | ds. (45)By combining the above equations, we can obtain M ∆ t (cid:107)∇ ¯ e n +1 / µ (cid:107) + 12 ( (cid:107) (∆ + β )¯ e n +1 φ (cid:107) − (cid:107) (∆ + β )¯ e nφ (cid:107) ) + λ (cid:107) ¯ e n +1 φ (cid:107) − (cid:107) ¯ e nφ (cid:107) )+ S ( (cid:107)∇ (¯ e n +1 φ − ¯ e nφ ) (cid:107) − (cid:107)∇ (¯ e nφ − ¯ e n − φ ) (cid:107) ) + ( e n +1 r ) − ( e nr ) ≤ M ∆ t (cid:107)∇ ¯ e n +1 / µ (cid:107) + C ∆ t ( (cid:107) ¯ e nφ (cid:107) + (cid:107) ¯ e n − φ (cid:107) + (cid:107) (∆ + β )¯ e nφ (cid:107) + (cid:107) (∆ + β )¯ e n − φ (cid:107) )+ C ∆ t (cid:107) φ t (cid:107) L ∞ ( J ; H − ) (cid:0) ( e n +1 r ) + ( e nr ) (cid:1) + C (cid:107) ˜ φ n +1 / (cid:107) m +1 N − m ∆ t + C ∆ t (cid:90) t n +1 t n (cid:107) ( − ∆) − / ∂ φ∂t ( s ) (cid:107) ds + C ∆ t (cid:90) t n +1 t n (cid:107) ∂ φ∂t ( s ) (cid:107) H ds + C ∆ t (cid:90) t n +1 t n (cid:107) ∂ φ∂t ( s ) (cid:107) H ds + C ∆ t (cid:90) t n +1 t n ( d rdt ( s )) ds Then summing over n , n = 0 , , . . . , k −
1, and using Gronwall’s inequality, we have k − (cid:88) n =0 ∆ t (cid:107)∇ ¯ e n +1 / µ (cid:107) + (cid:107) (∆ + β )¯ e kφ (cid:107) + (cid:107) ¯ e kφ (cid:107) + ( e kr ) ≤ C ∆ t (cid:90) t k ( (cid:107) ( − ∆) − / ∂ φ∂t (cid:107) + (cid:107) ∂ φ∂t (cid:107) H + (cid:107) ∂ φ∂t (cid:107) H + | d rdt | ) ds + C (cid:107) φ (cid:107) L ∞ ( J ; H m +1 per (Ω)) N − m , (46)By the triangle inequality and (23), we obtain the desired results (28). Remark We observe that with
S > , i.e., with an active stabilization term, we need to assumeadditionally ∂ φ∂t ∈ L (0 , T ; H per (Ω)) . But this additional regularity requirement might be formally derivedfrom the assumption ∂ φ∂t ∈ L (0 , T ; H − (Ω)) since one time derivative is formally equivalent to six spatialderivatives.
6. Numerical results.
We present in this section several numerical examples to verify the accuracyof our S-SAV Fourier-spectral schemes and to illustrate the phase transition behaviors and crystal growthin a supercooled liquid. In the following simulations, we take M = 1, β = 1. This manuscript is for review purposes only. XIAOLI LI AND JIE SHEN
In this subsection, we test the accuracy of the first and second order fullydiscrete schemes (9) and (10), respectively. In this test, we take Ω = (0 , × (0 , T = 1, S = 5, λ = 0 . (cid:15) = 0 .
025 and the initial condition φ = sin( πx/
16) cos( πy/ N and 2 N is calculated by (cid:107) e φ (cid:107) = (cid:107) φ N − φ N (cid:107) and we use the similar procedure to compute the error between twodifferent time steps ∆ t and ∆ t/ N = 256 so the spatial error is negligible. In Tables 1 and2, we show the order of convergence in the temporal direction. To test the spectral accuracy, we choose thetime step to be sufficiently small so that the error is dominated by the spatial error, and compute the errorat T = 1, the results are plotted in 1 which shows that the error converges exponentially. These numericalresults are consistent with the error estimates in Theorem 5. Table 1
Errors and convergence rates in time for the first order scheme (9) . ∆ t (cid:107) e φ (cid:107) L ∞ ( J ; L (Ω)) Rate (cid:107) e r (cid:107) L ∞ ( J ) Rate1 / /
10 6.77E-2 0.90 4.64E-2 0.971 /
20 3.60E-2 0.91 2.33E-2 0.991 /
40 1.88E-2 0.94 1.17E-2 1.001 /
80 9.67E-3 0.96 5.84E-3 1.00
Table 2
Errors and convergence rates in time for the second order scheme (10) . ∆ t (cid:107) e φ (cid:107) L ∞ ( J ; L (Ω)) Rate (cid:107) e r (cid:107) L ∞ ( J ) Rate1 / /
10 1.02E-2 1.99 1.79E-3 1.971 /
20 2.12E-3 2.26 4.41E-4 2.021 /
40 4.99E-4 2.09 1.11E-4 1.991 /
80 1.22E-4 2.03 2.78E-5 1.99
Spectral order N -12 -10 -8 -6 -4 -2 E rr o r i n l og sc a l e e φ e r Fig. 1 . Spatial L errors at time T = 1 for the second order scheme (10) In this subsection, we simulatethe phase evolution behavior of the PFC model. The physical parameters are set as Ω = (0 , × (0 , This manuscript is for review purposes only. he SAV and Fourier-spectral method for the PFC Equation φ i,j = φ + η i,j , where φ = 0 .
06 and η i,j is a uniformly distributed randomnumber satisfying | η i,j | ≤ .
01. The other parameters are (cid:15) = 0 . λ = 0 . S = 0 .
01 and withoutstabilization. It can be observed that the energy decreases at all times, which indicates numerical evidencefor our method being unconditionally energy stable. But the modified SAV energy at large ∆ t is notconsistent with the original energy without stabilization, it only becomes close to the original energy with∆ t = 0 .
02. However, with the stabilization parameter S = 0 .
01, even the result with ∆ t = 1 producesreasonably accurate results. This example shows that while the stabilization is not needed for stability, it isessential for accuracy at larger time steps. So in the following simulations, we always add a stabilized termso that reasonable accuracy can be achieved without using exceedingly small time steps. Time E ne r g y Original Energy with ∆ t=1Modified SAV Energy with ∆ t=1Stabilized Original Energy with ∆ t=1Stabilized Modified SAV Energy with ∆ t=1Original Energy with ∆ t=0.02Modified SAV Energy with ∆ t=0.02Stabilized Original Energy with ∆ t=0.02Stabilized Modified SAV Energy with ∆ t=0.02 Fig. 2 . Time evolution of the free energy functional.
We present the evolution of the density field φ calculated using the second order scheme (10) with∆ t = 1, S = 0 .
01 and N = 256 in Figure 3. Fig. 3 . The evolution of the density field φ calculated using the second order scheme at t = 150 , , , , , , respectively. This manuscript is for review purposes only. XIAOLI LI AND JIE SHEN
In this subsection, we simulate the crystal growth ina supercooled liquid with the following expression to define the crystallites: φ ( x l , y l ) = φ ave + C (cid:18) cos ( C √ y l ) cos ( C x l ) − . cos ( 2 C √ y l ) (cid:19) , l = 1 , , , (47)where x l and y l define a local system of cartesian coordinates that is oriented with the crystallite lattice,and the constant parameters φ ave , C and C take the values φ ave = 0 . C = 0 .
446 and C = 0 .
66. Thenwe define the initial configuration by setting three perfect crystallites in three small square patches whichare located at (350 , , , θ = − π , , π respectively: x l ( x, y ) = x sin( θ ) + y cos( θ ) , y l ( x, y ) = − x cos( θ ) + y sin( θ ) . (48)In this simulation, we take the parameters (cid:15) = 0 . M = 1, S = 0 . λ = 0 .
001 and T = 800. 512 Fourier modes are used to discretize the space and ∆ t = 0 .
05. Figure 4 demonstrates the evolution of thephase transition behavior using the second-order scheme (10) at different times t = 0, 50, 100, 200, 300,400, 500, 600, 800, respectively. One can observe the growth of the crystalline phase and the motion ofwell-defined crystal-liquid interfaces. Besides, we can see that the different alignment of the crystallitescauses defects and dislocations. Fig. 4 . The dynamical behaviors of the crystal growth in a supercooled liquid. Snapshots of the numerical approximationof the density field φ are taken at t = 0 , 50, 100, 200, 300, 400, 500, 600, 800, respectively.
7. Conclusion.
We developed fully discrete, unconditionally energy stable schemes based on the SAVand stabilized SAV approaches in time and the Fourier-spectral method in space for the phase field crystal(PFC) equation. We also carried out a rigorous error analysis which provided optimal error estimates inboth time and space. To the best of our knowledge, this is the first such result for any fully discrete linearschemes for the PFC model.
This manuscript is for review purposes only. he SAV and Fourier-spectral method for the PFC Equation
REFERENCES[1]
M. Ainsworth and Z. Mao , Analysis and approximation of a fractional Cahn-Hilliard equation , SIAM Journal onNumerical Analysis, 55 (2017), pp. 1689–1718.[2]
K. Elder and M. Grant , Modeling elastic and plastic deformations in nonequilibrium processing using phase fieldcrystals , Physical Review E, 70 (2004), p. 051605.[3]
K. Elder, M. Katakowski, M. Haataja, and M. Grant , Modeling elasticity in crystal growth , Physical review letters,88 (2002), p. 245701.[4]
H. Gomez and X. Nogueira , An unconditionally energy-stable method for the phase field crystal equation , ComputerMethods in Applied Mechanics and Engineering, 249 (2012), pp. 52–61.[5]
R. Guo and Y. Xu , Local discontinuous Galerkin method and high order semi-implicit scheme for the phase field crystalequation , SIAM Journal on Scientific Computing, 38 (2016), pp. A105–A127.[6]
Q. Li, L. Mei, X. Yang, and Y. Li , Efficient numerical schemes with unconditional energy stabilities for the modifiedphase field crystal equation , Advances in Computational Mathematics, (2019), pp. 1–30.[7]
X. Li, J. Shen, and H. Rui , Energy stability and convergence of SAV block-centered finite difference method for gradientflows , Mathematics of Computation, (2019).[8]
Y. Li and J. Kim , An efficient and stable compact fourth-order finite difference scheme for the phase field crystal equation ,Computer Methods in Applied Mechanics and Engineering, 319 (2017), pp. 194–216.[9]
J. Ramos , C. canuto, my hussaini, a. quarteroni, ta zang, spectral methods in fluid dynamics, springer-verlag, new york(1988), dm 162 , 1991.[10]
J. Shen, T. Tang, and L.-L. Wang , Spectral methods: algorithms, analysis and applications , vol. 41, Springer Science& Business Media, 2011.[11]
J. Shen and J. Xu , Convergence and error analysis for the scalar auxiliary variable (SAV) schemes to gradient flows ,SIAM Journal on Numerical Analysis, 56 (2018), pp. 2895–2912.[12]
J. Shen, J. Xu, and J. Yang , The scalar auxiliary variable (SAV) approach for gradient flows , Journal of ComputationalPhysics, 353 (2018), pp. 407–416.[13]
J. Shen and X. Yang , Numerical approximations of Allen-Cahn and Cahn-Hilliard equations , Discrete Contin. Dyn.Syst, 28 (2010), pp. 1669–1691.[14]
J. Swift and P. C. Hohenberg , Hydrodynamic fluctuations at the convective instability , Physical Review A, 15 (1977),p. 319.[15]
S. M. Wise, C. Wang, and J. S. Lowengrub , An energy-stable and convergent finite-difference scheme for the phasefield crystal equation , SIAM Journal on Numerical Analysis, 47 (2009), pp. 2269–2288.[16]
X. Yang , Efficient schemes with unconditionally energy stability for the anisotropic Cahn-Hilliard equation using thestabilized-Scalar Augmented Variable (S-SAV) approach , arXiv preprint arXiv:1804.02619, (2018).[17]
X. Yang and D. Han , Linearly first-and second-order, unconditionally energy stable schemes for the phase field crystalmodel , Journal of Computational Physics, 330 (2017), pp. 1116–1134.[18]
Z. Zhang, Y. Ma, and Z. Qiao , An adaptive time-stepping strategy for solving the phase field crystal model , Journal ofComputational Physics, 249 (2013), pp. 204–215., Journal ofComputational Physics, 249 (2013), pp. 204–215.