We now consider a number of time integration schemes taking part in the numerical comparisons presented below. All these schemes are well known and described in numerous literature and, in particular, in [64]. All schemes considered in this subsection approximate equation (2.1) and are conservative.

Throughout this paper space derivatives are discretized by finite differences. Assume that a uniform grid with grid size $h=1/N$ and nodes $x_{i}=(i-1/2)h$, $i=\overline{0,N+1}$, with $N$ being the number of grid cells, is used to discretize the equation in the domain $\Omega$. The domain boundaries $x=0$ and $x=1$ correspond to points $x_{1/2}=(x_{0}+x_{1})/2=0$ and $x_{N+1/2}=(x_{N}+x_{N+1})/2=Nh$ and the grid nodes $x_{1},x_{2},\ldots,x_{N}$ lie within the domain. For time discretization a uniform grid with time step size $\tauup$ and nodes $t_{n}=n\tauup$, $n=0,1,2,\ldots$, is used. Thus, the solution of the discretized problem is defined in points $(x_{i},t_{j})$, $i=\overline{1,N}$, $j=0,1,\ldots$. The corresponding values of the grid functions $c_{h}$ are denoted by $c_{i}^{n}=c(x_{i},t_{n})$.

Let $\Delta_{h}$ be a standard three-point finite-difference approximation of the Laplacian. In the 1D case under consideration, its values in the nodes $i=\overline{1,N-1}$ of the space grid are

$$\Delta_{h}c_{i}=\frac{1}{h^{2}}\left(c_{i+1}-2c_{i}+c_{i-1}\right).$$

The finite difference equations

$$c_{0}-c_{1}=0,\quad c_{N}-c_{N-1}=0$$

(3.1)

approximating the homogeneous Neumann boundary conditions (2.2) are related to the boundary nodes $i=0,N$. Then the resulting discrete Laplacian, with incorporated discretized boundary conditions, is denoted by $\Delta_{h}$.

Explicit Euler (EE) scheme, the simplest discretization of (2.1), reads

$$\displaystyle\frac{1}{M}\frac{c^{n+1}_{i}-c^{n}_{i}}{\tauup}=\Delta_{h}\muup^{n}_{i},\quad\muup^{n}_{i}=F^{\prime}(c^{n}_{i})-\epsilonup^{2}\Delta_{h}c^{n}_{i}.$$

The scheme is first order accurate in time and second order accurate in space, is stable for $\tauup\sim h^{4}$ and is not energy stable.

Implicit Euler (IE) scheme reads

$$\displaystyle\frac{1}{M}\frac{c^{n+1}_{i}-c^{n}_{i}}{\tauup}=\Delta_{h}\muup^{n+1}_{i},\quad\muup^{n+1}_{i}=F^{\prime}(c^{n+1}_{i})-\epsilonup^{2}\Delta_{h}c^{n+1}_{i}.$$

The scheme is first order accurate in time and second order accurate in space and is nonlinear (in the sense that a nonlinear equation has to be solved at each time step to determine $c^{n+1}_{i}$). It is conditionally (i.e., for a sufficiently small $\tauup$) gradient stable and conditionally uniquely solvable, see Table 1.

Crank-Nicolson (CN) scheme can be written as

$$\displaystyle\frac{1}{M}\frac{c^{n+1}_{i}-c^{n}_{i}}{\tauup}=\Delta_{h}\muup^{n+1/2}_{i},\quad\muup^{n+1/2}_{i}=\frac{1}{2}(\muup^{n}_{i}+\muup^{n+1}_{i}),$$

with $\muup_{i}^{n}$, $\muup_{i}^{n+1}$ defined above. The scheme is second order accurate in time and second order accurate in space and is nonlinear. It is conditionally gradient stable and conditionally uniquely solvable, see Table 1.

Semi-implicit Euler (SIE) reads

$$\displaystyle\frac{1}{M}\frac{c^{n+1}_{i}-c^{n}_{i}}{\tauup}=\Delta_{h}\muup^{n+1}_{i},\quad\muup^{n+1}_{i}=F^{\prime}(c^{n}_{i})-\epsilonup^{2}\Delta_{h}c^{n+1}_{i}.$$

The scheme is first order accurate in time and second order accurate in space and is linear. It is unconditionally uniquely solvable and is not energy stable, see Table 1.

The time integration schemes considered up to now are not energy stable and/or not uniquely solvable for arbitrary large time step sizes [12, 13]. An efficient and general way to design time integration schemes, which are both uniquely solvable and energy stable for arbitrary time step sizes, is proposed in papers of David J. Eyre [14, 15]. Since we essentially employ this approach to construct new schemes, we discuss it in more details here.

The essence of the Eyre approach is as follows. The system free energy density $\psiup$, c.f. (2.4), is not a convex function. This is because the dividing part of the free energy is a double-well potential. Assume that $\psiup$ can be split as

$$\psiup=\psiup_{\text{c}}+\psiup_{\text{e}},$$

(3.2)

where $\psiup_{\text{c}}$ and $-\psiup_{\text{e}}$ are convex functions. The subindices ‘‘c’’ (meaning ‘‘contraction’’) and ‘‘e’’ (meaning ‘‘expansion’’) reflect their physical meaning. Here the convexity of a function is understood as a positive semi-definiteness of its Hessian, see [13], i.e., a function $f(x)$ is convex provided that $f^{\prime\prime}(x)\equiv d^{2}f/dx^{2}\geqslant 0$. Then, the semi-discrete (discrete in time) approximations

$$\displaystyle\frac{1}{M}\frac{c^{n+1}-c^{n}}{\tauup}=\Delta\left[\muup_{\text{c}}(c^{n+1})+\muup_{\text{e}}(c^{n})\right],$$

(3.3)

with

$$\muup_{\text{c,e}}=\frac{\deltaup\psiup_{\text{c,e}}}{\deltaup c},\quad\muup=\muup_{\text{c}}+\muup_{\text{e}},$$

are gradient stable and

$$\psiup(c^{n+1})\leqslant\psiup(c^{n})$$

holds for any time step size $\tauup>0$. In addition, the convexity of the energy part$\psiup_{\text{c}}$ corresponding to the implicitly treated term, ensures the unique solvability of (3.3).

Consider the specific form of the energy functional (2.4). Taking into account that the gradient term in (2.4) belongs to the convex part of the energy, we can take in (3.2)

$$\psiup_{\text{c}}=\frac{1}{2}|\nabla c|^{2}+F_{\text{c}},\quad\psiup_{\text{e}}=F_{\text{e}},\quad F_{\text{c}}+F_{\text{e}}=F.$$

Thus, constructing a suitable splitting of the system energy is reduced to splitting the dividing part of the energy $F$ in such a way that

$$F=F_{\text{c}}+F_{\text{e}},\quad F^{\prime\prime}_{\text{c}}\geqslant 0,\quad-F^{\prime\prime}_{\text{e}}\geqslant 0.$$

(3.4)

Then, the chemical potential takes form

$$\muup=\muup_{\text{c}}+\muup_{\text{e}},\quad\muup_{\text{c}}=-\epsilonup^{2}\Delta c+F^{\prime}_{\text{c}},\quad\muup_{\text{e}}=F^{\prime}_{\text{e}},$$

whereas the semi-discrete system (3.3) can be written as

$$\frac{1}{M}\frac{c^{n+1}-c^{n}}{\tauup}=-\epsilonup^{2}\Delta^{2}c^{n+1}+\Delta F^{\prime}_{\text{c}}(c^{n+1})+\Delta F^{\prime}_{\text{e}}(c^{n}).$$

(3.5)

Splitting of the form (3.4) can be constructed in a number of ways. One of them is to write the dividing part $F$ of the free energy as a sum of two terms according to (3.2). Note that the gradient part of the free energy (2.5) should be handled together with the convex part $\psiup_{\text{c}}$. Another way is to choose a regularizing function $F_{\text{r}}$ such that

$$F_{\text{c}}=F+F_{\text{r}},\quad F_{\text{e}}=F-F_{\text{r}}.$$

Although both ways are essentially quite similar, it is convenient to distinguish them due to various possible generalizations.

Consider the specific form (2.5) of $F$ used in this work. We have

$$F(c)=c^{2}(1-c^{2}),\quad F^{\prime}(c)=2c(2c^{2}-3c+1),\quad F^{\prime\prime}(c)=2(6c^{2}-6c+1).$$

Equation $F^{\prime\prime}(c)=0$ has roots $c_{1,2}=\left(3\pm\sqrt{3}\right)/6$, and, hence, $F^{\prime\prime}(c)<0$ for $\left(3-\sqrt{3}\right)/6<c<\left(3+\sqrt{3}\right)/6$ ; $F^{\prime\prime}(c)\geqslant 0$ for $c\leqslant\left(3-\sqrt{3}\right)/6$ and $c\geqslant\left(3+\sqrt{3}\right)/6$. Function $F^{\prime\prime}(c)$ attains its minimum, equal to $-1$, at $c=1/2$. The maxima of $F^{\prime\prime}(c)$ for $c\in[0,1]$ are attained at $c=0$ and $c=1$ and equal $2$. Thus, to ensure that $F_{\text{c}}$ and $-F_{\text{e}}$ are convex, it suffices to set

$$F_{\text{c}}=F(c)+\frac{1}{2}c^{2},\quad F_{\text{e}}=-\frac{1}{2}c^{2}.$$

This leads to non-linearly stabilized splitting scheme (NLSS),

$$\displaystyle\frac{1}{M}\frac{c^{n+1}_{i}-c^{n}_{i}}{\tauup}=\Delta_{h}\muup^{n+1},\quad\muup^{n+1}_{i}=\left[F^{\prime}(c^{n+1}_{i})+c^{n+1}_{i}\right]-c^{n}_{i}-\epsilonup^{2}\Delta_{h}c^{n+1}_{i}.$$

(3.6)

The scheme can also be written as

$$\frac{1}{M}\frac{c^{n+1}_{i}-c^{n}_{i}}{\tauup}=\left(\Delta_{h}-\epsilonup^{2}\Delta_{h}^{2}\right)c^{n+1}_{i}+\Delta_{h}F^{\prime}(c^{n+1}_{i})-\Delta c^{n}_{i}.$$

It is first order accurate in time, second order in space, implicit, nonlinear in $c^{n+1}$ and energy stable for arbitrary $\tauup>0$.

Another option is to set

$$F_{\text{c}}=c^{2},\quad F_{\text{e}}=-c^{2}+F(c),$$

which leads to linearly stabilized splitting scheme (LSS),

$$\displaystyle\frac{1}{M}\frac{c^{n+1}_{i}-c^{n}_{i}}{\tauup}=\Delta_{h}\muup^{n+1}_{i},\quad\muup^{n+1}_{i}=2c^{n+1}_{i}+\left[F^{\prime}(c^{n}_{i})-2c^{n}_{i}\right]-\epsilonup^{2}\Delta_{h}c^{n+1}_{i}.$$

It may be convenient to write it as

$$\frac{1}{M}\frac{c^{n+1}_{i}-c^{n}_{i}}{\tauup}=\left(2\Delta_{h}-\epsilonup^{2}\Delta_{h}^{2}\right)c^{n+1}_{i}+\Delta_{h}F^{\prime}(c^{n}_{i})-2\Delta c^{n}_{i}.$$

The scheme is first order accurate in time, second order in space, implicit and, provided that splitting parameters are chosen properly [40], energy stable for arbitrary $\tauup>0$. Unlike the NLSS scheme (3.6), it is linear with respect to $c^{n+1}$. A three-parameter family of linearly stabilized schemes, including the LSS scheme, is analyzed in [13].

Properties of considered time integration schemes are summarized in Table 1.

For a fully discrete scheme (i.e., when discretization in time and in space is applied), its gradient stability implies that the free energy functional (2.4) is approximated in a certain way. In this work such an approximation, determined by the space disretization, reads

$$\psiup_{h}(\mathbf{c})=h\sum\limits_{i=1}^{N}F^{\prime}(c_{i})+\frac{1}{2}\epsilonup^{2}h\sum\limits_{i=1}^{N}\frac{(c_{i+1}-c_{i})^{2}}{h^{2}}.$$

(3.7)

Let a $N\times N$ matrix $\mathbf{A}$ be the matrix of the discussed above discretization $-\Delta_{h}$ of the operator $-\Delta$, with incorporated homogeneous Neumann boundary conditions (3.1). The space discretization reduces the initial boundary-value problem (2.1)–(2.3) to an initial value problem

$$\frac{1}{M}\mathbf{c}^{\prime}=-\mathbf{A}\left(F^{\prime}(\mathbf{c})+\epsilonup^{2}\mathbf{A}\mathbf{c}\right),\quad\mathbf{c}(0)=\mathbf{c}^{0},$$

(3.8)

where the entries of the vector function $\mathbf{c}(t):\mathbb{R}\rightarrow\mathbb{R}^{N}$ are the sought after solution values on the space grid and the vector $\mathbf{c}^{0}$ contains the grid values of the given function $c^{0}(x)$.

Let $\mathbf{c}^{n}$ be an $N$ vector containing the numerical solution values for time level $n$. Then, the time integration schemes discussed above can be written in a compact form as follows.

• •

  Explicit Euler (EE):

  $$\frac{\mathbf{c}^{n+1}-\mathbf{c}^{n}}{\tauup}=-\mathbf{A}\left(F^{\prime}(\mathbf{c}^{n})+\epsilonup^{2}\mathbf{A}\mathbf{c}^{n}\right);$$

  (3.9)

• •

  implicit Euler (IE):

  $$\frac{\mathbf{c}^{n+1}-\mathbf{c}^{n}}{\tauup}=-\mathbf{A}\left(F^{\prime}(\mathbf{c}^{n+1})+\epsilonup^{2}\mathbf{A}\mathbf{c}^{n+1}\right);$$

  (3.10)

• •

  Crank-Nicolson (CN):

  $$\frac{\mathbf{c}^{n+1}-\mathbf{c}^{n}}{\tauup}=-\frac{1}{2}\mathbf{A}\left(F^{\prime}(\mathbf{c}^{n+1})+\epsilonup^{2}\mathbf{A}\mathbf{c}^{n+1}+F^{\prime}(\mathbf{c}^{n})+\epsilonup^{2}\mathbf{A}\mathbf{c}^{n}\right);$$

  (3.11)

• •

  Semi-implicit Euler (SIE):

  $$\frac{\mathbf{c}^{n+1}-\mathbf{c}^{n}}{\tauup}=-\mathbf{A}\left(F^{\prime}(\mathbf{c}^{n})+\epsilonup^{2}\mathbf{A}\mathbf{c}^{n+1}\right);$$

  (3.12)

• •

  nonlinearly stabilized splitting (NLSS)

  $$\frac{\mathbf{c}^{n+1}-\mathbf{c}^{n}}{\tauup}=-\mathbf{A}\left(F^{\prime}(\mathbf{c}^{n+1})+\mathbf{c}^{n+1}-\mathbf{c}^{n}+\epsilonup^{2}\mathbf{A}\mathbf{c}^{n+1}\right);$$

  (3.13)

• •

  linearly stabilized splitting (LSS)

  $$\frac{\mathbf{c}^{n+1}-\mathbf{c}^{n}}{\tauup}=-\mathbf{A}\left(F^{\prime}(\mathbf{c}^{n})-2\mathbf{c}^{n}+2\mathbf{c}^{n+1}+\epsilonup^{2}\mathbf{A}\mathbf{c}^{n+1}\right).$$

  (3.14)

We now briefly discuss computational complexity of these schemes. In implicit nonlinear schemes (i.e., IE, CN, NLSS), each time step a nonlinear system of equations has to be solved, e.g., by a Newton iteration. In turn, each Newton iteration involves solving the Jacobian linear system. Implicit linear schemes SIE and LSS require solving a single linear system per time step. Explicit schemes involve neither nonlinear nor linear system solution and have lowest costs per time step. Although implicit time integration schemes allow to use significantly larger time step sizes than explicit schemes, solving large linear systems, within Newton iterations or independently, may become a computationally expensive task. This is because efficiency in sparse direct solvers, as well as in preconditioned iterative solvers, is difficult to retain on modern high-performance hybrid CPU/GPU platforms.

The aim of this work is to develop an efficient numerical algorithm which, on one hand, is gradient stable and allows large time step sizes and, on the other hand, can be efficiently implemented on modern high-performance hybrid platforms. Our approach is based on the combination of the Eyre splitting, in particular the LSS scheme, with the local iteration modified (LIM) scheme. Then, the Eyre splitting allows to use large time step sizes while keeping gradient stability, whereas the LIM mechanism provides a computationally efficient alternative to solving linear systems in implicit schemes.

The local iteration modified (LIM) scheme can be viewed as a special explicit scheme where a Chebyshev polynomial of a sufficiently high order $p$ is employed to ensure stability and monotonicity for a given time step size $\tauup>0$. We now describe a LIM variant based on the linearly stabilized splitting (LSS) scheme (3.14), which we write as

$$\mathbf{c}^{n+1}=\left[\mathbf{I}+\tauup\widehat{\mathbf{A}}\right]^{-1}\widehat{\mathbf{f}}^{n},\quad\widehat{\mathbf{A}}=\mathbf{A}(2\mathbf{I}+\epsilonup^{2}\mathbf{A}),\quad\widehat{\mathbf{f}}^{n}=\mathbf{c}^{n}+\tauup\mathbf{A}(2\mathbf{c}^{n}-F^{\prime}(\mathbf{c}^{n})),$$

(3.15)

where $\mathbf{I}$ is the $N\times N$ identity matrix. We emphasize that the inverse matrix here is, of course, never computed, rather a linear system with this matrix is solved. Within the LIM approach the inverse matrix is replaced by a specially chosen Chebyshev polynomial. Let $\lambdaup_{\infty}$ be an upper bound for the largest eigenvalue of $\widehat{\mathbf{A}}$ (in practice one can set $\lambdaup_{\infty}:=\|\widehat{\mathbf{A}}\|_{1}=\max_{j}\sum_{i}|\widehat{a}_{ij}|$). The polynomial order $p$ for which stability takes place is chosen as

$$p=\left\lceil\frac{\piup}{4}\sqrt{\tauup\lambdaup_{\infty}+1}\,\right\rceil,$$

(3.16)

where $\lceil x\rceil$, for $x\in\mathbb{R}$, denotes the smallest integer greater than or equal to $x$. Let the Chebyshev polynomial roots

$$\{\betaup_{m},\;m=1,\ldots,p\}=\left\{\cos\piup\frac{2i-1}{2p},\;i=1,\dots,p\right\}$$

be ordered in such a way that no numerical instability occurs for arbitrarily large $p$ and $\betaup_{1}=\cos({\piup}/{2p})$ is the first root. We set $z_{1}=\betaup_{1}$ and define Chebyshev polynomial parameters

$$a_{m}=\frac{\lambdaup_{\infty}}{1+z_{1}}(z_{1}-\betaup_{m}),\quad m=1,\dots,p.$$

Denoting $\mathbf{y}^{(0)}=\mathbf{c}^{n}$, we obtain the LIM solution $\mathbf{c}^{n+1}$ on the next time level by carrying out $2p-1$ Chebyshev iterations as

	$\displaystyle\mathbf{y}^{(m)}$	$\displaystyle=\frac{1}{1+\tauup a_{m}}\left(\mathbf{c}^{n}+\tauup a_{m}\mathbf{y}^{(m-1)}+\tauup(\widehat{\mathbf{f}}^{n}-\widehat{\mathbf{A}}\mathbf{y}^{(m-1)})\right),\quad m=1,\dots,p,$		(3.17)
	$\displaystyle\mathbf{y}^{(p+m-1)}$	$\displaystyle=\frac{1}{1+\tauup a_{m}}\left(\mathbf{c}^{n}+\tauup a_{m}\mathbf{y}^{(m-1)}+\tauup(\widehat{\mathbf{f}}^{n}-\widehat{\mathbf{A}}\mathbf{y}^{(m-1)})\right),\quad m=2,\dots,p,$
	$\displaystyle\mathbf{c}^{n+1}$	$\displaystyle=\mathbf{y}^{(2p-1)}.$

Note that $a_{1}=0$ and it is easy to check that the first iteration solution $\mathbf{y}^{(1)}$ is the explicit Euler solution on the next time level. In this sense, the Chebyshev iterations actually start at the second step by computing $\mathbf{y}^{(2)}$. If only the first sweep of iterations is carried out in (3.17), with $\mathbf{y}^{(1)}$, …, $\mathbf{y}^{(p)}$ computed, then setting $\mathbf{c}^{n+1}=\mathbf{y}^{(p)}$ we get a time integration scheme which is known as the regular (non-modified) local iteration scheme (the LI scheme). As already mentioned, in the regular LI scheme stability of the numerical solution is guaranteed for parabolic problems but the solution nonnegativity is not. As one can see, in the LIM scheme (3.17) the second Chebyshev iteration sweep, with $\mathbf{y}^{(p+1)}$, …, $\mathbf{y}^{(2p-1)}$ being computed, the same parameters $a_{m}$ are employed, except that the first iteration with $a_{1}=0$ (which gives the explicit Euler solution) is not carried out. One can show [58] that formulas (3.17) can be written in an operator form as

$$\mathbf{c}^{n+1}=(\mathbf{I}-\mathbf{F}_{p}^{2})\left[\mathbf{I}+\tauup\widehat{\mathbf{A}}\right]^{-1}\widehat{\mathbf{f}}^{n},$$

(3.18)

where $\widehat{\mathbf{A}}$ and $\widehat{\mathbf{f}}^{n}$ are defined in (3.15) and $\mathbf{F}_{p}$ is the Chebyshev polynomial operator,

$$\mathbf{F}_{p}=\prod_{m=p}^{m=1}\left(\mathbf{I}-\frac{1}{1+\tauup a_{m}}(\mathbf{I}+\tauup\widehat{\mathbf{A}})\right).$$

By replacing $\mathbf{F}_{p}^{2}$ in (3.18) by $\mathbf{F}_{p}$ we obtain an operator form of the (regular) LI scheme. Thus, as we see, squaring the Chebyshev polynomial leads to monotonicity in the LIM scheme.