Improving the accuracy of convexity splitting methods for gradient flow equations

Abstract

This paper introduces numerical time discretization methods which significantly improve the accuracy of the convexity-splitting approach of Eyre (1998) [7], while retaining the same numerical cost and stability properties.

A first order method is constructed by iteration of a semi-implicit method based upon decomposing the energy into convex and concave parts. A second order method is also presented based on backwards differentiation formulas. Several extrapolation procedures for iteration initialization are proposed. We show that, under broad circumstances, these methods have an energy decreasing property, leading to good numerical stability.

The new schemes are tested using two evolution equations commonly used in materials science: the Cahn–Hilliard equation and the phase field crystal equation. We find that our methods can increase accuracy by many orders of magnitude in comparison to the original convexity-splitting algorithm. In addition, the optimal methods require little or no iteration, making their computation cost similar to the original algorithm.

Introduction

High order parabolic equations which derive from gradient descent of an energy functional (“gradient flows”) are commonly encountered in the study of material microstructure evolution and related phase field models [1], [2]. Two extensively studied models in this class are the Cahn–Hilliard (CH) equation [3]

$$u_t=\mathrm{\Delta }(-{ϵ}^2\mathrm{\Delta }u+u^3-u),$$

and phase field crystal (PFC) equation [4]

$$u_t=\mathrm{\Delta }({(\mathrm{\Delta }+1)}^{2}u-ru+u^3).$$

Numerical methods with explicit time discretization are impractical for these equations because of limitations imposed by numerical stability. This has motivated the development of a variety of implicit and semi-implicit numerical methods [5], [6], [7], [8], [9], [10], [11]. A simple and elegant semi-implicit approach was formulated by Eyre [7]. Although it was originally introduced for the Cahn–Hilliard equation, it has been successfully implemented in many different contexts [11], [12], [13], [14]. Various improvements and extensions of the method have also been made [15], [16], [17].

In the algorithm's most basic form, an energy function $F(u)$ is decomposed into convex and concave parts as $F=F_{+}+F_{-}$. For a gradient flow of the form $u_t=-\mathrm{\nabla }F(u)$, the method treats the convex part implicitly, and the concave part explicitly, leading to the formally first order method

$$\frac{u_{n+1}-u_n}{h}=-\mathrm{\nabla }{F}_{+}(u_{n+1})-\mathrm{\nabla }{F}_{-}(u_n),$$

where $u_n$ is an approximation to $u(nh)$ and h is the timestep. For this scheme, it can be shown that the energy decreasing property holds in a discrete way: $F(u_{n+1})\le F(u_n)$. An algorithm with this property is said to be energetically stable. Under some circumstances, this notion of stability coincides with the usual idea of unconditional numerical stability.

The main advantage to the convexity splitting method is that $\mathrm{\nabla }{F}_{+}$ can often be chosen to be linear, and moreover $I+h\mathrm{\nabla }{F}_{+}$ is often easy to invert. This means that (3) can be accomplished very efficiently, without the need for Newton-iteration steps and various sorts of linear algebra typical of fully implicit methods. In addition, practical implementations in phase field modeling appears to produce qualitatively correct features, such as the development and propagation of interfaces. On the other hand, it has been widely observed that the convexity-splitting scheme can have large temporal errors [9], [18], [19]. The purpose of this paper is to improve the accuracy of the convexity splitting approach while retaining all of its other benefits.

Our new methods fall into a broad class of semi-implicit algorithms for high order parabolic equations [9], [20], [21], [22]. There have been several recent advances for these types of methods in the context of materials science. The approach proposed by Hu et al. [22] used finite differences and a multi-grid approach to solve the PFC equation. Gomez and Hughes [9] introduced a second order accurate in time variational method for the CH equation, using mixed finite elements in space. Second order accurate algorithms have also been used for the PFC equation by Vignal et al. [21]. The work in [21] presented a computational framework that relies on the convex-concave splitting approach. Christlieb et al. [20] suggested that high accuracy can be achieved using fully implicit methods with a pseudo-spectral spatial discretization, employing Newton iterations and a conjugate gradient solver.

In this paper, we utilize the convexity splitting of the free energy as a basis for new algorithms, based upon various combinations of iteration, extrapolation, and higher order discretization. In section 2, we review the variational formulation of abstract gradient flow equations and explain stability of convexity splitting methods in this context. Iterative and higher order algorithms are introduced in section 3, which are shown to be energetically stable under broad circumstances. These methods are tested in section 4 for the Cahn–Hilliard and phase field crystal equations.