SAV Method Applied to Fractional Allen-Cahn Equation

Fractional differential equations are becoming increasingly used as a modeling tool for processes associated with anomalous diffusion or spatial heterogeneity. We consider in this paper a nonlocal phase transition model, in particular described by the Allen-Cahn equation. Precisely, we deal with a fractional Allen-Cahn equation (FACE) with spectral definition of the fractional Laplacian. For the space discretization, we first use spectral Galerkin method to solve the eigenvalue problem (EVP) for the Laplacian operator with homogeneous Neumann boundary condition, then we employ the computed eigenfunctions as trial and test bases to solve the FACE. For the time discretization, based on the scalar auxiliary variable (SAV) approach, we construct unconditionally second-order energy stable BDF scheme (SAV/BDF2). Finally, the developed algorithm is used to simulate the phase evolution of a benchmark problem.


Introduction
The Allen-Cahn equation was originally introduced to describe the motion of antiphase boundaries in crystalline solids [1]. There have been a large body of work on numerical analysis of Allen-Cahn equations (cf. [2][3][4][5] and the references therein). We aim in this paper to use the SAV scheme, recently introduced and analyzed by a number of researchers; see, e.g., [5]  (1.1) In the above, γ is a positive kinetic coefficient, s ∈ (0, 1), ⊂ R d is a bounded domain, n is the outward normal, f (φ) = F (φ) with a given function F (φ) = 1 4ε 2 (φ 2 − 1) 2 being the Ginzburg-Landau double-well potential. The phase field φ is such that φ = ⎧ ⎨ ⎩ 1, phase 1, −1, phase 2, and ε represents the thickness of the smooth transition layer connecting the two phases, which is small compared to the characteristic length of the system scale. The homogeneous Neumann boundary condition implies that no mass loss occurs across the boundary walls.
Among the different definitions of fractional Laplacians (see [6,7] for a quantitative assessment of new numerical methods as well as available state-of-theart methods for discretizing the fractional Laplacians problems), we choose in this paper to focus on the fractional spectral definition. It is defined by where λ i , e i are the eigenvalues and eigenfunctions of the Laplace operator − in with homogeneous Neumann boundary condition, i.e., they satisfy While, a i represents the projection of u on the direction e i , a i = (u, e i ) L 2 . The spectral fractional Laplacian is nonlocal on the interior for noninteger s ∈ (0, 1). We see that to compute the inner product a i = (u, e i ) L 2 , it suffices for u to be defined on the interior of . No information about u on the exterior R d \ is required. Thus, from a conceptual viewpoint, in boundary value problems the spectral fractional Laplacian can admit the same type of boundary conditions as the standard, local Laplacian − . In this paper, we let  The rest of this paper in organized as follows. In Sect. 2, we present briefly the spectral method by giving some notations and reminders. The fractional Laplace operator and its possible applications is discussed in Sect. 3. To demonstrate the applicability of the approximative fractional Laplacian for real applications, we consider a fractional Allen-Cahn equation (FACE). Based on the scalar auxiliary variable (SAV) approach, we construct an unconditionally second-order energy stable BDF scheme (SAV/BDF2) for FACE. We present numerical results for a test case as well as a benchmark example in Sect. 4.

Spatial Discretizations
We limit here the description of the spectral approximation to the introduction of some notations and reminders (see [8,9]). For complex domain, we can use spectral element method [10]. Let = {(ξ i , ρ i ); 0 ≤ i ≤ N } denote the sets of Gauss-Lobatto-Legendre quadrature nodes and weights associated to polynomials of degree N . These quantities are such that on where P N ( ) denotes the space of polynomials of degree ≤ N . We recall that the nodes The canonical polynomial interpolation basis h i (x) ∈ P N ( ) built on is given by the relationships: with the elementary cardinality property where δ ij is Kronecker's delta symbol.
In the sequel the phase field φ will be approximated in space variable by suitable polynomial functions φ N as follows The L 2 -inner products involved in the calculation will be achieved using Gauss-Lobatto-Legendre quadrature, which reads: for all continuous functions ϕ and ψ in ,

Scalar Auxiliary Variable (SAV) Approach for FACE
SAV approach was introduced in [4,5] to solve gradient flows. The main purpose of this section is to construct efficient unconditionally stable scheme based on this approach for (1.1). Throughout the paper, we assume there exists a constant C 0 such that F (φ)dx + C 0 > 0. We first introduce a scalar auxiliary variable Then, we rewrite the phase-field equation (1.1) under an equivalent form as: find φ : (0, T ] × → R and r : (3.1) , then we have the following energy dissipation law Proof Taking the inner product of the first two equations with μ, ∂φ ∂t respectively, and multiplying the third equation with 2r(t), then adding them together, We obtain Given initial conditions φ 0 = φ 0 , and let r 0 = F (φ 0 )dx + C 0 , find φ n+1 ∈ H s ( ) and r n+1 ∈ R, n = 1, . . . , such that In the above,φ n+1 can be any explicit approximation of φ(t n+1 ) with an error of O( t 2 ). For instance, we may choose the following onē φ n+1 = 2φ n − φ n−1 .

Theorem 3.2 The scheme (3.5)-(3.7) is unconditionally stable in the sense that
(3.9) Proof The result can be directly deduced from taking the inner product of the first two equations (3.5) and (3.6) with μ n+1 and 3φ n+1 −4φ n +φ n−1 2 t respectively, and multiplying the third equation (3.7) with 2r n+1 , then using the following identity:

Implementation
Besides its unconditional stability, a most remarkable feature of the above scheme is that it can be solved very efficiently. Indeed, by inserting (3.6) and (3.7) into (3.5), and letF n+1 := F (φ n+1 )dx + C 0 , we obtain where We shall first determine (f (φ n+1 ), φ n+1 ) from (3.10). To this end we multiply 2F n+1 (f (φ n+1 ), β n+1 ) = (f (φ n+1 ), α n+1 ), (3.13) Then, we have . (3.14) Thus we obtain an expression to compute φ n+1 by bringing back (3.14) into (3.10): Finally we compute r n+1 through We now summarize the algorithm of the Scalar Auxiliary Variable approach/Semi-Implicit Second-Order Scheme (3.5)-(3.7) as follows: In this section, we first present a numerical example to illustrate the efficiency of the SAV scheme in terms of stability and accuracy. We then use the proposed scheme to simulate a benchmark problem.

Test of the Convergence Order
In order to validate the proposed SAV/BDF2 scheme for the fractional phase-field equation, we consider a fabricated forcing term so that the exact solution to (1.1) is φ(x, t) = sin(t) cos(π x) cos(πy). In this test we set γ = 1, =] − 1, 1[ 2 , and the nonlinear term is given by In the calculation we use polynomial degree 32×32 for the spatial discretization, which is large enough so that the spatial discretization error is negligible compared to the temporal error. Figure 1 shows the L 2 -errors at T = 1.0 in log-log scale as a function of the time step size for several fractional orders. It is observed from this figure that the convergence rate of the time stepping scheme is exactly second order as expected for all tested values of s. It is worthy to mention that no numerical instability was observed for all time step sizes used in the calculation. This implies that the proposed scheme is unconditionally stable.

Benchmark Test
In this subsection, we apply the SAV/BDF2 scheme to the fractional version of a classical benchmark problem (cf. [11]) that we describe below. Our main purpose in this test is to demonstrate the applicability of the constructed method for the FACE. We are particularly interested in numerically investigating the impact of the fractional order on the evolution of the phase interface.
At the initial state, there is a circular phase interface of the radius R 0 = 100 in the rectangular domain ] − 128, 128[ 2 . In other words, the initial condition is given by Such a circular interface is unstable and the driving force will make it shrink and eventually disappear. It has been shown that in the limit that the radius of the circle is much larger than the interfacial thickness, the velocity and the radius of the moving interface are given (see [1]) by In the implementation we map the computational domain ] − 128, 128[ 2 to ] − 1, 1[ 2 . Therefore actually we are led to solve the fractional Allen-Cahn equation (1.1) with the coefficients γ = 1/128 2 and ε = 0.0078. In the simulation, the space resolution is set to N = 512, and the time step size is t = 0.1. The computed radius R(t) for s = 1 using the SAV/BDF2 scheme is plotted in Fig. 2. We observe that R(t) keeps monotonously decreasing and very close to the sharp interface limit value. This confirms the accuracy of the proposed method, at least in the case s = 1.
Next we apply the proposed scheme to investigate the impact of the fractional order on the radius behavior. In Fig. 3 we present the numerical radius evolution for a number of the fractional orders. Specifically, Fig. 4 shows the circle shrinking for fractional orders s = 1.0, 0.9, 0.8. It is clearly indicated that the radius decay rate slow down when the fractional order decreases. However, for the time being the physical meaning and mathematical explanation of this phenomena remain unknown. We plan to address this issue in future work.  The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.