Parallel energy-stable phase field crystal simulations based on domain decomposition methods

TL;DR

This paper introduces a parallel, energy-stable phase field crystal simulation algorithm that employs domain decomposition and adaptive time stepping, achieving high scalability and second-order accuracy in 2D and 3D cases.

Contribution

It develops a novel parallel Newton–Krylov–Schwarz method with improved boundary conditions for efficient phase field crystal simulations.

Findings

Second-order accuracy in space and time.

Energy stability with large time steps.

Scalability to over ten thousand cores.

Abstract

In this paper, we present a parallel numerical algorithm for solving the phase field crystal equation. In the algorithm, a semi-implicit finite difference scheme is derived based on the discrete variational derivative method. Theoretical analysis is provided to show that the scheme is unconditionally energy stable and can achieve second-order accuracy in both space and time. An adaptive time step strategy is adopted such that the time step size can be flexibly controlled based on the dynamical evolution of the problem. At each time step, a nonlinear algebraic system is constructed from the discretization of the phase field crystal equation and solved by a domain decomposition based, parallel Newton--Krylov--Schwarz method with improved boundary conditions for subdomain problems. Numerical experiments with several two and three dimensional test cases show that the proposed algorithm is second-order accurate in both space and time, energy stable with large time steps, and highly scalable to over ten thousands processor cores on the Sunway TaihuLight supercomputer.

Tables

Table 1. Table 1: Performance of the NKS with different subdomain solvers. Here n/c means the divergence of the GMRES solver.

Subdomain solver		ILU(0)	ILU(1)	ILU(2)	LU	ILU-reuse
Test A	Total Newton	30	30	30	30	30
	GMRES/Newton	6.0	3.93	3.87	3.87	6.0
	Total Time (s)	3.17	3.83	5.39	40.43	2.01
Test B	Total Newton	n/c	30	30	30	30
	GMRES/Newton	n/c	101.03	8.17	7.87	8.4
	Total Time (s)	n/c	30.47	8.29	56.75	5.12
Test C	Total Newton	34	33	33	33	34
	GMRES/Newton	5.88	3.94	3.91	3.91	6.56
	Total Time (s)	9.99	11.28	14.86	191.75	6.02
Test D	Total Newton	31	31	31	31	31
	GMRES/Newton	7.29	7.23	7.23	7.23	7.42
	Total Time (s)	24.01	284.68	982.22	683.94	10.41

Table 2. Table 2: Performance of NKS with different types of preconditioner and different overlapping size.

Preconditioner		classical-AS			left-RAS		right-RAS
$\delta$		$0$	$1$	$2$	$1$	$2$	$1$	$2$
Test A	Total Newton	30	30	30	30	30	30	30
	GMRES/Newton	10.0	6.0	6.1	2.67	3.0	2.67	3.0
	Total Time (s)	2.20	2.01	2.28	1.57	1.78	1.55	1.79
Test B	Total Newton	30	30	30	30	30	30	30
	GMRES/Newton	113.23	8.4	84.73	8.7	55.13	7.53	37.7
	Total Time (s)	31.31	5.12	33.18	5.21	22.13	4.86	16.07
Test C	Total Newton	34	34	34	34	34	34	34
	GMRES/Newton	13.82	6.56	6.56	2.91	2.88	2.91	2.88
	Total Time (s)	8.25	6.02	6.66	4.62	5.11	4.63	5.10
Test D	Total Newton	32	31	32	31	31	31	31
	GMRES/Newton	12.88	7.42	5.0	2.06	2.06	2.06	2.06
	Total Time (s)	3.69	10.41	35.06	8.27	30.49	8.24	30.51

Table 3. Table 3: Performance results with different number of processor cores, here NP means the number of processor cores.

NP	216	432	864	1,728	3,456
Total Newton	34	34	34	34	34
GMRES/Newton	2.88	2.88	2.88	2.88	2.88

Table 4. Table 4: Performance results with different number of processor cores.

NP	1,536	3,072	6,144	12,288	24,576
Total Newton	32	32	31	31	31
GMRES/Newton	12.63	12.66	11.97	12.65	12.65

Equations

$$F(\varphi )={\int }_{\Omega }\left\{\frac{1-\gamma }{2}{\varphi }^2+\frac{1}{4}{\varphi }^4-|\nabla \varphi {|}^2+\frac{1}{2}(\Delta \varphi {)}^2\right\}\mathrm{d}\mathbf{x},$$

$$\mathbf{n}\cdot \nabla \varphi {|}_{\partial \Omega }=\mathbf{n}\cdot \nabla (\Delta \varphi ){|}_{\partial \Omega }=\mathbf{n}\cdot \nabla ({\Delta }^2\varphi ){|}_{\partial \Omega }=0$$

$$G(\varphi )=\frac{1-\gamma }{2}{\varphi }^2+\frac{1}{4}{\varphi }^4-|\nabla \varphi {|}^2+\frac{1}{2}(\Delta \varphi {)}^2.$$

$${\varphi }_t=\nabla \cdot \left(M(\varphi )\nabla \frac{\delta G}{\delta \varphi }\right),$$

$$\frac{\delta G}{\delta \varphi }$$

$$=(1-\gamma )\varphi +{\varphi }^3+2\Delta \varphi +{\Delta }^2\varphi .$$

$$\frac{d}{dt}F=\frac{d}{dt}{\int }_{\Omega }G\mathrm{d}\mathbf{x}$$

$$={\int }_{\Omega }\left(\frac{\partial G}{\partial \varphi }\frac{\partial \varphi }{\partial t}+\frac{\partial G}{\partial {\varphi }_x}\frac{\partial {\varphi }_x}{\partial t}+\frac{\partial G}{\partial {\varphi }_y}\frac{\partial {\varphi }_y}{\partial t}+\frac{\partial G}{\partial {\varphi }_{xx}}\frac{\partial {\varphi }_{xx}}{\partial t}+\frac{\partial G}{\partial {\varphi }_{yy}}\frac{\partial {\varphi }_{yy}}{\partial t}\right)\mathrm{d}\mathbf{x}$$

$$={\int }_{\Omega }\frac{\delta G}{\delta \varphi }{\varphi }_t\mathrm{d}\mathbf{x}+B_1={\int }_{\Omega }\frac{\delta G}{\delta \varphi }\nabla \cdot \left(M(\varphi )\nabla \frac{\delta G}{\delta \varphi }\right)\mathrm{d}\mathbf{x}+B_1$$

$$=-{\int }_{\Omega }M(\varphi )\nabla \frac{\delta G}{\delta \varphi }\cdot \nabla \frac{\delta G}{\delta \varphi }\mathrm{d}\mathbf{x}+B_1+B_2\le 0,$$

$${\int }_{\Omega }\left(G_{\delta }-G(\varphi ,{\varphi }_x,{\varphi }_y,{\varphi }_{xx},{\varphi }_{yy})\right)\mathrm{d}\mathbf{x}$$

$$={\int }_{\Omega }\left(\frac{\partial G}{\partial \varphi }\delta \varphi +\frac{\partial G}{\partial {\varphi }_x}\delta {\varphi }_x+\frac{\partial G}{\partial {\varphi }_y}\delta {\varphi }_y+\frac{\partial G}{\partial {\varphi }_{xx}}\delta {\varphi }_{xx}+\frac{\partial G}{\partial {\varphi }_{yy}}\delta {\varphi }_{yy}\right)\mathrm{d}\mathbf{x}$$

$$+O\left((\delta \varphi {)}^2\right)$$

$$={\int }_{\Omega }\frac{\delta G}{\delta \varphi }\delta \varphi \mathrm{d}\mathbf{x}+B_3+O\left((\delta \varphi {)}^2\right)$$

$$D_x^{+}{\varphi }_{i,j}=\frac{{\varphi }_{i+1,j}-{\varphi }_{i,j}}{\Delta x},D_x^{-}{\varphi }_{i,j}=\frac{{\varphi }_{i,j}-{\varphi }_{i-1,j}}{\Delta x},$$

$$D_x{\varphi }_{i,j}=\frac{{\varphi }_{i+\frac{1}{2},j}-{\varphi }_{i-\frac{1}{2},j}}{\Delta x},D_x^{⟨2⟩}{\varphi }_{i,j}=\frac{{\varphi }_{i+1,j}-2{\varphi }_{i,j}+{\varphi }_{i-1,j}}{(\Delta x{)}^2},$$

$${\nabla }_d=(D_x,D_y)$$

$${\Delta }_d={\nabla }_d\cdot {\nabla }_d=D_x^{⟨2⟩}+D_y^{⟨2⟩},$$

$$(\nabla \cdot M(\varphi )\nabla {)}_d={\nabla }_d\cdot M({\varphi }_{i,j}){\nabla }_d=D_{x}M({\varphi }_{i,j})D_x+D_{y}M({\varphi }_{i,j})D_y,$$

$$G_d({\varphi }_{i,j}^{(n)})$$

$$-\left(\frac{(D_x^{+}{\varphi }_{i,j}^{(n)}{)}^2+(D_x^{-}{\varphi }_{i,j}^{(n)}{)}^2}{2}+\frac{(D_y^{+}{\varphi }_{i,j}^{(n)}{)}^2+(D_y^{-}{\varphi }_{i,j}^{(n)}{)}^2}{2}\right).$$

$$F_d^{(n)}=\sum _{i=1}^{N_x}\sum _{j=1}^{N_y}{G}_d({\varphi }_{i,j}^{(n)})\Delta x\Delta y.$$

$$G_{\delta }-G(\varphi ,{\varphi }_x,{\varphi }_y,{\varphi }_{xx},{\varphi }_{yy}):=G_d({\varphi }_{i,j}^{(n+1)})-G_d({\varphi }_{i,j}^{(n)}).$$

$$\sum _{i,j}\left(G_d({\varphi }_{i,j}^{(n+1)})-G_d({\varphi }_{i,j}^{(n)})\right)\Delta x\Delta y$$

$$=\sum _{i,j}\frac{\delta G_d}{\delta ({\varphi }^{(n+1)},{\varphi }^{(n)}{)}_{i,j}}\left({\varphi }_{i,j}^{(n+1)}-{\varphi }_{i,j}^{(n)}\right)\Delta x\Delta y,$$

$$\frac{\delta G_d}{\delta ({\varphi }^{(n+1)},{\varphi }^{(n)}{)}_{i,j}}=\left[(1-\gamma )+2{\Delta }_d+{\Delta }_d^2\right]\left(\frac{{\varphi }_{i,j}^{(n+1)}+{\varphi }_{i,j}^{(n)}}{2}\right)$$

$$+\frac{({\varphi }_{i,j}^{(n+1)}{)}^3+({\varphi }_{i,j}^{(n+1)}{)}^2{\varphi }_{i,j}^{(n)}+{\varphi }_{i,j}^{(n+1)}({\varphi }_{i,j}^{(n)}{)}^2+({\varphi }_{i,j}^{(n)}{)}^3}{4}.$$

$$\frac{{\varphi }_{i,j}^{(n+1)}-{\varphi }_{i,j}^{(n)}}{\Delta t_n}={\nabla }_d\cdot M_d{\nabla }_d\frac{\delta G_d}{\delta ({\varphi }^{(n+1)},{\varphi }^{(n)}{)}_{i,j}}.$$

$$F_d({\varphi }^{(n+1)})\le F_d({\varphi }^{(n)}).$$

$$\sum _{i=1}^{N_x}{g}_{i,j}(D_x{M}_{i,j}{D}_x)g_{i,j}=-\sum _{i=1}^{N_x}\frac{M_{i,j}}{2}\left[(D_x^{+}{g}_{i,j}{)}^2+(D_x^{-}{g}_{i,j}{)}^2\right],$$

$$\sum _{j=1}^{N_y}{g}_{i,j}(D_y{M}_{i,j}{D}_y)g_{i,j}=-\sum _{j=1}^{N_y}\frac{M_{i,j}}{2}\left[(D_y^{+}{g}_{i,j}{)}^2+(D_y^{-}{g}_{i,j}{)}^2\right].$$

$$\frac{F_d({\varphi }^{(n+1)})-F_d({\varphi }^{(n)})}{\Delta t_n}$$

$$=\sum _{i=1}^{N_x}\sum _{j=1}^{N_y}\frac{\delta G_d}{\delta ({\varphi }^{(n+1)},{\varphi }^{(n)}{)}_{i,j}}\left(\frac{{\varphi }_{i,j}^{(n+1)}-{\varphi }_{i,j}^{(n)}}{\Delta t_n}\right)\Delta x\Delta y$$

$$=\sum _{i=1}^{N_x}\sum _{j=1}^{N_y}\left(\frac{\delta G_d}{\delta ({\varphi }^{(n+1)},{\varphi }^{(n)}{)}_{i,j}}\right){\nabla }_d\cdot M_{i,j}{\nabla }_d\left(\frac{\delta G_d}{\delta ({\varphi }^{(n+1)},{\varphi }^{(n)}{)}_{i,j}}\right)\Delta x\Delta y$$

$$=-\sum _{i=1}^{N_x}\sum _{j=1}^{N_y}\frac{M_{i,j}}{2}\left[(D_x^{+}\mathcal{A}{)}^2+(D_x^{-}\mathcal{A}{)}^2+(D_y^{+}\mathcal{A}{)}^2+(D_y^{-}\mathcal{A}{)}^2\right]\Delta x\Delta y\le 0.$$

$$\Delta t_n=\max \left(\Delta t_{min},\frac{\Delta t_{max}}{\sqrt{1+\eta |F^{\prime }(t){|}^2}}\right),$$

$$\mathcal{F}(\mathbf{\Phi })=0$$

Full text

††journal: Elsevier

1 Introduction

The phase field crystal (PFC) equation is a popular model for simulating microstructural evolution in material sciences. As an atomic description of crystalline materials on the diffusive time scale, the PFC equation was originally proposed to model the dynamics of crystal growth by Elder et. al. [1, 2]. Since then, during the past decade, it has been applied with significant successes for the simulation of phenomena found in various solid-liquid systems such as the crystal growth in a supercooled liquid [3, 4, 5], the crack propagation in a ductile material [6, 2], the dendritic and eutectic solidification [7, 6], and the epitaxial growth [7, 8]. The PFC equation is usually derived from the free-energy functional that originates from the more advanced density functional theory of Hohenberg and Kohn [9], resulting in a six-order nonlinear partial differential equation of the atomic density function. Due to the existence of the nonlinearity and the high-order terms, the PFC equation requires high-fidelity simulations with advanced numerical methods.

Together with the introduction of the PFC model, an explicit Euler method was employed to solve the PFC equation successfully by Elder et. al. [2, 7], in which the time step size is in proportion to the sixth order of grid sizes. The computational cost of the explicit method is high due to the small time step for maintaining the stability. To relax the restriction of the small time step, some implicit methods were proposed, see, e.g., [10, 11, 5, 12]. Although relatively large time step size can be employed in these implicit methods, no energy stability analyses were presented. Recently, many researches were done to design energy stable time-stepping algorithms that are fully free from the time step constraint due to the stability condition; examples include the energy stable finite difference methods [13, 14, 15, 3, 4] and the energy stable finite element methods [6, 16, 17]. But because of the lack of the supportive theory, the design of these methods are often problem-dependent and not easy to generalize. In this paper, we design a second-order discretization scheme for the PFC equation based on the discrete variational derivative (DVD) method [18] to naturally achieve the energy stability. And by exploiting the DVD method, the numerical scheme can be designed in a general way, which is suitable for the PFC equation with different types of the mobility and boundary conditions. Further more, we employ an adaptive time stepping as a companion of the unconditionally energy stable method so as to adjust the time step size without losing the accuracy.

When an implicit discretization of the PFC equation is applied, a sparse linear or nonlinear algebraic system arises at each time step. And the solution of the discretized system could be time consuming especially when a fine mesh is used. It is therefore of great importance to study highly efficient solvers to accelerate the simulation at large scale. Despite the fact that some effective approaches, such as the multigrid method, have been applied in solving the discretized PFC equations [17, 12, 14], dedicated studies on efficient parallel solution algorithms are less to be seen. In this paper, we propose a highly scalable parallel solver based on the Newton–Krylov–Schwarz (NKS) algorithm [19] with modified subdomain boundary conditions for solving the nonlinear algebraic system arising at each implicit time step. Several key parameters in the solver, including the type of the Schwarz preconditioner, the size of the overlap, and the solver for subdomain problems, are discussed and tested to achieve the optimal performance. We show by experiments that the proposed solver can scale well to over ten thousands processor cores.

The remainder of this paper is organized as follows. The PFC equation is introduced in Sec. 2. In Sec. 3, we employ the DVD method to obtain an unconditionally stable scheme for the PFC equation. In Sec. 4, we introduce the NKS algorithm to solve the nonlinear system at each time step. Several numerical simulations are reported in Sec. 5 and concluding remarks are given in Sec. 6.

2 Phase Field Crystal Equation

The free-energy functional in the PFC model takes the following dimensionless form [1, 20]:

[TABLE]

where $\gamma>0$ is the quench depth for supercooling the material and $\phi:\Omega\rightarrow\mathbb{R}$ is the density distribution function to approximate the number density of atoms. In the PFC model, the quench depth is proportional to the deviation of the temperature from the melting temperature and the density distribution function is conserved during the non-equilibrium process. For simplicity, we focus the discussion on a two dimensional rectangle domain $\Omega=[0,L_{x}]\times[0,L_{y}]$. The three dimesional cases will be studied in the numerical experiments. We impose the doubly periodic boundary conditions or Neumann-type boundary conditions

[TABLE]

for $\phi$, with $\mathbf{n}$ being the outward normal of $\partial\Omega$.

Based on (1), a local energy functional $G(\phi)$ can be defined as

[TABLE]

Then the PFC equation takes the following form [18]

[TABLE]

where $M(\phi)\geq 0$ is the mobility, and

[TABLE]

is the variational derivative of the local energy $G(\phi)$. For the PFC equation, it can be proved that the time derivative of the free-energy $F(\phi)$ satisfies

[TABLE]

which demonstrates the energy dissipative of the system. Here $B_{1}$ and $B_{2}$ are the boundary terms coming from the integration-by-parts formula. In particular, we have $B_{1}=B_{2}=0$ for the periodic and Neumann-type boundary conditions.

3 The discretization of the PFC equation

In the section, we use the DVD method to construct a numerical scheme of the PFC equation. By denoting $G_{\delta}$ as $G(\phi+\delta\phi,...,\phi_{yy}+\delta\phi_{yy})$, the following integral relationship between the variational derivative and the local energy holds

[TABLE]

The first equality comes from the Taylor expansion formula. The second equality comes from the integration-by-parts formula, and the boundary term $B_{3}=0$ under the given boundary conditions. Eq. (6) plays an important role in the DVD method due to the fact that it shows the connection between the discrete variational derivative and the discrete local energy. 􏰙

We use a uniform mesh of $N_{x}\times N_{y}$ elements with mesh sizes $\Delta x=L_{x}/N_{x}$ and $\Delta y=L_{y}/N_{y}$ to cover the computational domain $\Omega$. The solution $\phi$ is approximated as $\phi_{i,j}\approx\phi(x_{i},y_{j})$, in which $(x_{i},y_{j})=\left((i-\frac{1}{2})\Delta x,(j-\frac{1}{2})\Delta y\right)$, $1\leq i\leq N_{x},1\leq j\leq N_{y}$. $\phi_{i,j}^{(n)}$ denotes the numerical solution at $n$-th time step corresponding to time $t_{n}$. We introduce some useful notations as following

[TABLE]

[TABLE]

and $D_{y}^{+},D_{y}^{-},D_{y},D_{y}^{\langle 2\rangle}$ are defined similarly. In this article, we use the subscript $d$ to indicate the corresponding discrete forms of operators and functions. And we denote

[TABLE]

as the discrete gradient operator. The Laplacian operator $\Delta$ is discretized by

[TABLE]

and operator $\nabla\cdot M(\phi)\nabla$ is discretized by

[TABLE]

where the value of $M$ on the midpoint of an edge is approximated by the averaged value of $M$ on the two adjacent nodes of the edge; for instance, $M(\phi_{i+\frac{1}{2},j})\approx\frac{M(\phi_{i+1,j})+M(\phi_{i,j})}{2}$.

To compute the discrete variational derivative, we first define the discrete local energy and discrete free-energy for the PFC equation. The discrete approximation of the local free-energy $G$ in Eq. (2) is

[TABLE]

We define the discrete free-energy $F_{d}$ at time $t_{n}$ as

[TABLE]

Taking $\delta\phi:=\phi^{(n+1)}_{i,j}-\phi^{(n)}_{i,j}$, we can obtain

[TABLE]

Then the discretization of Eq. (6) is derived as a summation formula

[TABLE]

where $\frac{\delta G_{d}}{\delta(\phi^{(n+1)},\phi^{(n)})_{i,j}}$ expresses the discrete variational derivative obtained by solving Eq. (10). According to Eq. (7) and (10), we can obtain the discrete variational derivative as

[TABLE]

When a first-order finite difference formula is applied, we may obtain a semi-implicit scheme for Eq. (3) as

[TABLE]

Here, $\Delta t_{n}=t_{n+1}-t_{n}$ is the time step size, and $M_{d}=M\left(\frac{\phi_{i,j}^{(n+1)}+\phi_{i,j}^{(n)}}{2}\right)$. The scheme can be viewed as the standard Crank–Nicolson scheme, which is trivial for the linear terms, but is non-trivial when dealing with the nonlinear terms. Numerical simulations presented later will show that the scheme (12) has numerical second-order accuracy in time and space [21]. It’s worth mentioning that the discrete variational derivative in Eq. (11) derived by DVD method happens to have a similar form with the numerical scheme for PFC equation in [3]. However, the numerical scheme presented in [3] employ a first-order forward difference instead of a second-order centered one to discrete the gradient operator, therefore the scheme is formally first-order accurate in space in order to maintain the energy stability.

The scheme (12) has a crucial feature that the dissipation property is kept no matter what time step size is adopted. The proof can be derived analogously to the continuous case. In addition, it can be easily proved that the scheme (12) keeps the conservation property of the density distribution function $\phi$ during the whole time evolution process.

Proposition. Under the periodic or Neumann-type boundary conditions, the numerical scheme (12) is unconditionally energy stable. More precisely, for any time step size $\Delta t>0$, the solution of the scheme (12) satisfies the energy dissipation

[TABLE]

Proof. With the periodic boundary conditions, we present two vital formulas that will be used in the proof

[TABLE]

The Eq. (14) can be regarded as a discrete analogue of the integration-by-parts formula and their demonstration is omitted here. According to Eq. (14), the dissipation property of the discrete free-energy $F_{d}$ is given by

[TABLE]

In the third equality, we simplify the notation $\frac{\delta G_{d}}{\delta(\phi^{(n+1)},\phi^{(n)})_{i,j}}$ as $\mathcal{A}$ for brevity. The demonstration can be derived analogously for the Neumann-type boundary conditions.

The unconditional stability of the semi-implicit scheme (12) is obtained from the dissipation property. But a blind increase of the time step size is adverse for keeping the computational accuracy. Numerical experiments show that using a large constant time step may produce nonphysical solutions [3]. This is because the PFC equation, similar to other phase-field equations, contains multiple time scales that may vary in orders of magnitude during the coarsening and phase separation processes. Therefore an adaptive control of the time step size is necessary in the numerical simulation, in which the time step size is selected based on the desired solution accuracy and the dynamic features of the system.

In this paper, we exploit the adaptive time step strategies described in [3], the time step size is computed by

[TABLE]

where $\lvert F^{\prime}(t)\rvert=|\frac{F_{d}^{(n)}-F_{d}^{(n-1)}}{\Delta t_{n-1}}|$ corresponds to the change rate of free-energy on the two previous time steps, and $\eta$ is chosen to adjust the level of adaptivity in which large (small) $\eta$ indicates strengthening (easing) the restriction to the time step size. $\Delta t_{min}$ and $\Delta t_{max}$ are defined as the lower and upper bounds of the time step size, namely $\Delta t_{min}\leq\Delta t_{n}\leq\Delta t_{max}$. In the PFC model, the free-energy $F$ decays quickly at the early stage of dynamics because of the nonlinear interaction, and then decays rather slowly until it reaches a steady state. By using the adaptive time step strategy in Eq. (16), the time step size is adjusted promptly based on the change rate of the free-energy $F$, in which the large $\lvert F^{\prime}(t)\rvert$ leads to the small time step size, and the small $\lvert F^{\prime}(t)\rvert$ yields the large time step size.

4 Newton–Krylov–Schwarz Algorithm

In the semi-implicit scheme (12), the PFC equation is discretized into a nonlinear system

[TABLE]

at each time step. In the paper, we solve the nonlinear system (17) on a parallel supercomputers by adopting an NKS type algorithm [19]. The NKS algorithm consists of three important components: (1) an inexact Newton method as the outer iteration; (2) a Krylov method as an inner iteration for the linear Jacobian system at each Newton iteration; and (3) a Schwarz preconditioner to improve the Krylov method.

At each time step, the nonlinear system (17) is solved by an inexact Newton method. The solution of the previous time step is set to be the initial guess which has a great impact on the convergence of the iteration. At the $(m+1)$-th step of the inexact Newton iteration, the new approximate solution $\mathbf{\Phi}_{m+1}$ is obtained from the current approximate solution $\mathbf{\Phi}_{m}$ through

[TABLE]

Here $\lambda_{m}$ is the step length determined by a line search procedure [22], and $\mathbf{S}_{m}$ is the search direction obtained by approximately solving the following linear Jacobian system

[TABLE]

where $J_{m}=\frac{\partial{\cal F}(\mathbf{\Phi}_{m})}{\partial\mathbf{\Phi}_{m}}$ is the Jacobian matrix. The stopping condition for the Newton iteration (18) is

[TABLE]

where $\varepsilon_{r},\varepsilon_{a}\geq 0$ are the relative and absolute tolerances for the nonlinear iterations, respectively. Compared to the classical Newton method, the inexact method is superior especially when the number of unknowns is large (e.g., of the order of millions or larger) due to the reason that the linear Jacobian system is solved approximately instead of exactly, leading to a substantial reduction of the computational cost.

To accelerate the convergence of the linear solver, a right-preconditioned linear system

[TABLE]

is solved instead of the original Jacobian system (19). In our study, the Generalized Minimal Residual (GMRES) method [23] is applied to approximately solve the right-preconditioned linear system (21) until the linear residual $\mathbf{r}_{m}=J_{m}\mathbf{S}_{m}+{\cal F}(\mathbf{\Phi}_{m})$ satisfies the stopping condition

[TABLE]

where $\xi_{r},\xi_{\alpha}\geq 0$ are the relative and absolute tolerances for the linear iterations, respectively. In the GMRES method, the additive Schwarz type preconditioner $H_{m}^{-1}$ is the key to the success of the linear solver. To define the additive Schwarz type preconditioner $H_{m}^{-1}$ , we first partition the computational domain $\Omega$ into $np$ non-overlapping subdomains $\Omega_{k},\ (k=1,2,\cdots,np)$, then extend each subdomain by $3\delta$ mesh layers to form an overlapping decomposition $\Omega=\cup_{k=1}^{np}\Omega_{k}^{\delta}$, in which $3$ is the stencil width of the finite difference scheme for the PFC equation. The classical additive Schwarz preconditioner [24] is defined as

[TABLE]

Here $R_{k}^{\delta}$ serves as a restriction operator that maps a vector to a new one that is defined in the subdomain $\Omega_{k}^{\delta}$, by discarding the components outside $\Omega_{k}^{\delta}$; $(R_{k}^{\delta})^{T}$ represents an extension operator that maps a vector defined in the subdomain $\Omega_{k}^{\delta}$ to a new one that is defined in the whole domain, by putting zeros at the components outside $\Omega_{k}^{\delta}$.

There are two modified versions of the AS preconditioner that may have some potential advantages. The first one is the left restricted additive Schwarz (left-RAS, [25]) preconditioner that reads

[TABLE]

The only difference between the left-RAS preconditioner and the AS preconditioner is the extension operator. Instead of $(R_{k}^{\delta})^{T}$, the left-RAS preconditioner uses $(R_{k}^{0})^{T}$ which puts zeros not only outside $\Omega_{k}^{\delta}$ but also outside $\Omega_{k}$. The other modification to the AS preconditioner is the right restricted additive Schwarz (right-RAS, [26]) preconditioner that is given by

[TABLE]

The only difference between the right-RAS preconditioner and the AS preconditioner is the restriction operator. Instead of $R_{k}^{\delta}$, the right-RAS preconditioner uses $R_{k}^{0}$ which ignores the entries outside $\Omega_{k}$ when doing the extension.

In Eq. (22)-(24), the subdomain matrix $A_{k}$ can be directly generated as

[TABLE]

But people usually do not employ the method in Eq. (25) due to the high computational cost caused by the usage of the global matrix $J_{m}$ in the formula. Likewise, we use the discretization of the subdomain problem to generate $A_{k}$. Therefore, we need to consider what boundary conditions should be imposed on the subdomain boundaries. In our approach, we employ the boundary conditions as follows

[TABLE]

where $\Omega_{k}^{\delta+1}$ is the $3$ mesh layers expansion region of $\Omega_{k}^{\delta}$ that ensures all mesh points in $\Omega_{k}^{\delta}$ can acquire sufficient information to perform the stencil calculations when solving the subdomain problem. After defining suitable interface conditions for the subdomain problems, we then solve them either directly by using a sparse LU factorization or approximately by using a sparse incomplete LU (ILU) factorization.

A great advantage of the additive Schwarz preconditioners is that communication only occurs between neighboring subdomains during the restriction and extension processes. The major cost of the additive Schwarz preconditioners is the subdomain solves which are done sequentially without any inter-process communication. Therefore the additive Schwarz preconditioners is naturally suitable to parallel computing as long as the number of iterations is kept low. We further remark that compared to the classical AS preconditioner, the communication in the two restricted versions is reduced approximately by half because only the restriction or the extension step requires communication. And experiment results show that the restricted Schwarz preconditioners also have similar levels of convergence rate to the classic case, if not superior to. This may further improve the performance of the preconditioner.

5 Numerical experiments

In this section, we investigate the numerical behavior and parallel performance of the proposed algorithm. We begin with several numerical tests for the PFC equation to validate the discretization of the proposed algorithm. In addition, we investigate different performance-related parameters in the NKS algorithm to obtain the best performance. We consider four test cases, including three two dimensional (2D) test cases and a three dimensional (3D) test case. We mainly focus on (1) the verification of the numerical accuracy of the semi-implicit method, (2) the parallel performance of the semi-implicit method for various parameters, including the subdomain solvers and preconditioners, and (3) the parallel scalability of the proposed algorithm.

We perform our numerical experiments on the Sunway TaihuLight supercomputer, which tops the TOP–500 list as of June, 2016, with a peak performance greater than 100 PFlops. The computing power of TaihuLight is provided by a homegrown many–core SW26010 CPU [27], in which we only enable one core per CPU socket for the current study. The algorithm for the PFC equation is implemented on top of the Portable, Extensible Toolkits for Scientific computations (PETSc, [28]) library. The stopping conditions for the nonlinear and linear iterations are as follows.

The relative tolerance for the nonlinear iteration: $\varepsilon_{r}=1\times 10^{-8}$.

• 2.

The absolute tolerance for the nonlinear iteration: $\varepsilon_{a}=1\times 10^{-10}$.

• 3.

The relative tolerance for the linear iteration: $\xi_{r}=1\times 10^{-3}$.

• 4.

The absolute tolerance for the linear iteration: $\xi_{a}=1\times 10^{-11}$.

To test the accuracy of our scheme, we define the relative $l_{2}$ error as follows:

[TABLE]

where $\phi_{i,j}$ is the numerical solution, $\tilde{\phi}_{i,j}$ is the analytical solution.

5.1 Validation of the proposed method

A. Crystal growth in a 2D supercooled liquid

First, we use the PFC equation to model the crystal growth in a 2D supercooled homogeneous liquid. The simulation is conducted on a periodic square domain $\Omega=[0,128]^{2}$ with a random initial value $\phi_{i,j}^{(0)}=\bar{\phi}+\sigma_{i,j}$ in which $\bar{\phi}=0.07$ is the average density of the liquid-crystal system and $\sigma_{i,j}\in[-0.07,0.07]$ is chosen randomly. The positive parameter $\gamma$ is set to be $0.025$, and the mobility $M(\phi)$ is set to be $1$ for convenience.

Similar simulations were reported in [5, 14, 4]. We perform the simulation on a $256\times 256$ uniform mesh. The time step size is adaptively controlled by using the adaptive time step strategy with $\Delta t_{min}=0.01,\Delta t_{max}=20$ and the parameter $\eta=400,000$.

Fig. 1 shows the pseudocolor plots of the density distribution at times $t=1$, $500$, $1,200$ and $3,200$, from which we observe: (1) the fluid quickly crystallizes under the supercooling before $t=1,200$; (2) after $t=1,200$ the crystallized material gradually stabilizes as a solid lattice with periodic hexagonal pattern; (3) the domain is filled perfectly with periodic regular hexagonal pattern at $t=3,200$. We present the scaled total free-energy $F/(L_{x}L_{y})$ and the history of the time step size in Fig. 2, respectively. As shown in Fig. 2 (a), the free-energy decreases monotonically to the minimal as the solution evolves to the steady state. We observe that the free-energy decays quickly at the early stage and then decays rather slowly. From Fig. 2 (b), we observe that the time step size is rapidly adjusted from $\Delta t_{min}$ to $\Delta t_{max}$ due to the change of the free-energy. The time step is controlled by the variation of the free-energy on the two previous time step, in which the small time step means that the free-energy varies quickly and the large time step indicates that the free-energy varies slowly, which is consistent with the variation of the free-energy in Fig. 2 (a). Compared to the other publications of the PFC equation [14, 3, 17], we obtain the same solution that ascertains the rightness of the semi-implicit scheme.

Next, we study the accuracy of the proposed scheme. Because the random initial value is inappropriate for the check of the grid convergence, we consider the following smooth initial value

[TABLE]

where $\Omega=[0,32]^{2}$ is periodic. The positive parameter $\gamma$ and mobility $M(\phi)$ is unchanged. To understand the accuracy of the spatial discretization, we run the test on gradually refined meshes. Since the exact density distribution function $phi$ is unknown, the numerical solution on a very fine mesh $1,024\times 1,024$ and small time step size $\Delta t=0.01$ is taken as the analytical solution $\tilde{\phi}_{i,j}$. The $l_{2}$ errors at time $t=5$ for several mesh sizes ($N_{x}=N_{y}=N$) are plotted in Fig. 3 (a) which clearly shows the pesented semi-implicit scheme has a second-order accuracy in space.

Next we fix the spatial mesh as $512\times 512$ to study the accuracy of the semi-implicit scheme in time. The numerical solution at a fixed time step size $\Delta t=0.0005$ is regarded as the analytical solution. The $l_{2}$ errors at time $t=5$ are illustrated in Fig. 3 (b), which validates the second-order accuracy of the semi-implicit scheme in time.

B. Polycrystalline growth in a 2D supercooled liquid

The second test case is the polycrystalline growth in a 2D supercooled liquid in which the growth of different orientated crystallites is studied. In the simulation, three initial crystallites with hexagonal pattern oriented in different directions are seeded in the liquid. Similar numerical experiments were reported in [1, 14, 6, 17]. The simulation is conducted on a square domain $\Omega=[0,400]^{2}$, in which a uniform mesh with $1,024\times 1,024$ elements is applied. The positive parameter $\gamma$ is set to be $0.25$. The time step size is adaptively controlled with $\eta=5,000,\Delta t_{min}=0.02$, and $\Delta t_{max}=10$. To define the initial value, we first set the density function $\bar{\phi}$ to be a constant $0.285$ in the computational domain, and then replace three hexagonal lattice crystallites in three $25\times 25$ square patches of the domain.

We use the following periodic function to describe the two-dimensional hexagonal lattice structure

[TABLE]

where $A=0.446$ is a constant which represents an amplitude of the fluctuations in density, and $q_{x}=q_{y}=0.66$ are the period-related parameters. To obtain the orientated crystallites, we define a local Cartesian coordinate system $\tilde{\mathbf{x}}$ that is obtained by rotating the original Cartesian coordinate $\mathbf{x}=(x,y)$ with a certain angle $\omega$. Here $\tilde{\mathbf{x}}$ is defined as

[TABLE]

We imply rotation angles $\omega=-\pi/4,0,\pi/4$ on three crystal patches, respectively.

We consider four test cases, which are (a) periodic boundary conditions with $M(\phi)=1$; (b) Neumann-type boundary conditions with $M(\phi)=1$; (c) periodic boundary conditions with $M(\phi)=1-\phi^{2}$; and (d) Neumann-type boundary conditions with $M(\phi)=1-\phi^{2}$. Fig. 4 shows the distributions of the density function under these situations. The solutions $\phi$ are displayed at times $t=0,\,200,\,430,\,10,000$, from which we observe the growth of the crystallines. We observe a common phenomena that the three hexagonal lattice crystallites grow separately and then fuse with each other as time goes on. The defects and the dislocations caused by the different alignment of the crystallites can be clearly observed in Fig. 4. In the cases with Neumann-type boundary conditions, the “X” shape defects and the dislocations finally appear in the domain. In the cases with periodic boundary conditions, the hexagonal lattice crystallites can go through the boundaries resulting to the appearance of defects and dislocations near the boundaries. It is worth to note that the results of $M(\phi)=1$ and $M(\phi)=1-\phi^{2}$ are almost the same due to $\phi<1$, which agree well with the results reported in [17].

We show the scaled total free-energy and the history of the time step size of the tests in Fig. 5. For the four tested scenarios, the total free-energy decreases monotonically at a similar pathway.

Meanwhile, the time step sizes are successfully adjusted from $\Delta t_{min}$ to $\Delta t_{max}$, and have almost same line graph. In this sense, we summarize that these two kinds of boundary conditions and mobilities have little impact on the free-energy in test case B.

C. Crack propagation in a 2D ductile material

In the test case, we employ the PFC equation to model the crack propagation in a periodic rectangle domain with ductile material. Similar simulations can be found in [2, 6]. In order to define the initial value, we first set a crystal lattice given by the expression

[TABLE]

in the computational domain. Here, $q_{x}$ and $q_{y}$ are the parameters which determine the crystal period. We take $q_{x}=\sqrt{3}/2$, and $q_{y}=0.9q_{x}$, which means that the initial crystal has no stretching in the $x$ direction and a 1/9 stretching in the $y$ direction.

The numerical simulation is conducted on a maximum periodic system $591\frac{2\pi}{q_{x}}\times 76\frac{2\pi\sqrt{3}}{q_{y}}$ contained in domain $4,096\pi/3\times 1,024\pi/3$. The computational mesh is composed of $4,096\times 1,024$ elements. In the centre of computational domain, a notch of size $40􏰰\Delta x\times􏰯20\Delta y$ is cut out and replaced with a coexisting liquid ($\phi=0.79$) [2]􏱍. The notch provides a nucleating cite for a crack to start propagating. In this simulation, the parameter $\gamma$, $M(\phi)$ equal $1$, and the time step size is adaptively controlled with $\Delta t_{min}=1,\Delta t_{max}=5$ and $\eta=100$.

Fig. 6 depicts the pseudocolor plots of the density distribution at three times, $t=0$ (the initial shape), and $t=8,000,\,16,000$. We observe that the crack grows from a little rectangle and keeps growing outward like a tree on the endings. We show the scaled total free-energy and the history of the time step size in Fig. 7. It can be seen that there is no increase in free-energy, and the time step size increases from $\Delta t_{min}$ to $\Delta t_{max}$.

D. Polycrystalline growth in a 3D supercooled liquid

The experiment is conducted to simulate the polycrystalline growth in three dimensional space, which can be regarded as the 3D version of test B. Similar simulation was reported in [16]. The computational domain is 􏰕$[0,80\pi]^{3}$, and periodic boundary conditions are imposed in all directions. An uniform mesh comprised of $256\times 256\times 256$ elements is used, and the time step size is controlled adaptively with $\Delta t_{min}=0.5,\Delta t_{max}=10$, and $\eta=500$. Three initial crystallite spheres with BCC configuration oriented in different directions are placed in the domain. We can predict that the grain boundaries emerge eventually when the crystallites meet due to the inconformity of orientations.

Analytically, the BCC configuration is defined as [29, 30]

[TABLE]

where $\mathbf{x}=(x,y,z)$ represents the point in the three-dimensional Cartesian coordinate system and $q$ represents a wavelength related to the BCC crystalline structure. In our simulation, a regular crystallite has the form as follows

[TABLE]

where $\bar{\phi}$ represents the average density of the liquid-crystal system, and $A$ represents an amplitude of the fluctuations in density. The scaling function $\rho(\mathbf{x})$ is defined as

[TABLE]

where $\mathbf{x}_{0}=(x_{0},y_{0},z_{0})$ is the center of the crystallite sphere, and $d_{0}$ is the radius of the crystallite sphere. In order to define the crystallites oriented in different direction, we replace $\mathbf{x}$ with $\tilde{\mathbf{\mathbf{x}}}$, in which a system of local Cartesian coordinates $\tilde{x},\tilde{y},\tilde{z}$ is used to generate the crystallites in different directions. We make an affine transformation of the global coordinates to produce a rotation with an angle $\omega$ along the $z$-axis, in which the local Cartesian coordinates are given as follows

[TABLE]

In the experiment, we situate the three crystallite spheres on the $z=L/2$ plane with the same radius $d_{0}=5\pi$, where $L=80\pi$. The centers of the three crystallite spheres are $(L/2,3L/4,L/2)$, $(L/4,L/4,L/2)$, $(3L/4,L/4,L/2)$, respectively. The rotation angles $\omega$ for the three crystallite spheres are $0,-\pi/8,\pi/8$, respectively. In the simulation, $\bar{\phi}=-0.35$, $q=1/\sqrt{2}$, $A=1$, $\gamma=0.35$ and $M(\phi)=1$ are used.

The numerical solutions at different times are shown in Fig. 8. Similarly to the two dimensional case, the crystallite spheres grow in the liquid, and the grain boundaries appear when the crystals meet due to the orientation mismatch.

The free-energy evolution for the simulation is shown in Fig. 9, in which we can see that the free-energy has no increase. The time step size keeps the lower bound $\Delta t_{min}$ because the change of free-energy is rapid in the whole time interval $[0,650]$. To observe the grain boundaries clearly, our initial crystallites are fixed on the same plane and only have rotations along the $z$-axis, if the readers are interested in other situations, the positions and the rotation angles both can be changed freely and similar simulations can be done easily.

5.2 Performance tuning

It is well known that the performance of the Schwarz preconditioner depends on the choice of the subdomain solver. To investigate it, we run the above four test cases for the first $10$ time steps, respectively, in which the periodic boundary conditions are imposed on computational, and the mobility $M(\phi)$ is fixed to be $1$ for convenience. To check the influence of the subdomain solvers, we limit the test to the classical AS preconditioner and fix the overlapping size to $\delta=1$. The ILU factorizations with [math], $1$, and $2$ levels of fill-in and LU factorizations are considered. The number of processor cores, the mesh size, and the time step size for the four tests are listed as follows.

In the test $A$, $24$ processor cores with a $256\times 256$ mesh and a fixed time step size $\Delta t=0.01$ are used.

• 2.

In the test $B$, $288$ processor cores with a $1,024\times 1,024$ mesh and a fixed time step size $\Delta t=1$ are applied.

• 3.

In the test $C$, $576$ processor cores with a $4,096\times 1,024$ mesh and a fixed time step size $\Delta t=1$ are used.

• 4.

In the test $D$, $4,096$ processor cores with a $128\times 128\times 128$ mesh and a fixed time step size $\Delta t=0.1$ are applied. And the computational domain is scaled down to $[0,40\pi]^{3}.$

The numbers of Newton and GMRES iterations together with the total compute time are provided in Table 1.

For all test cases, the number of Newton iterations is insensitive to the subdomain solver. It is clear that by increasing the fill-in level, the number of GMRES iterations decreases, but the compute time does not necessarily reduce due to the increased cost of the subdomain solver. In summary, we find that the optimal choice in terms of the total compute time is the ILU(0) or ILU(2) subdomain solver. In the Newton method, the Jacobian matrices of the each Newton iteration have very similar structures, so it is possible to save the compute time by only performing the factorization once and reusing the preconditioner matrices within the all Newton iteration. Here, we apply the reuse strategy to the optimum subdomain solver which is ILU(0) for test cases $A$, $C$, $D$ and is ILU(2) for test case $B$. The results are listed in the last column of Table 1, which indicates that the reuse strategy can save nearly 50% of the compute time.

We then investigate the performance of the NKS solver by changing the type of the AS preconditioner and the overlapping factor $\delta$. The number of processor cores, the mesh size and the time step size for the four test cases are the same with the previous simulation. Based on the previous report, we take the optimal choice of subdomain solver with the reuse strategy throughout the test cases. The classical-AS (22), the left-RAS (23), and the right-RAS (24) preconditioners with overlapping size $\delta$ gradually increasing from 0 to 2 are considered. The numbers of Newton and GMRES iterations together with the total compute time are listed in Table 2.

From the Table, we conclude that the left-RAS and right-RAS preconditioner are superior to the classical-AS preconditioner. For the test cases $A$ and $B$, the minimal compute time and the least number of GMRES iterations are obtained when the overlapping size $\delta=1$. For the test case $C$, the minimal compute time is achieved when the overlapping size $\delta=1$ but the least number of GMRES iterations is obtained at $\delta=2$. For the test $D$, the minimal compute time and the greatest number of GMRES iterations are obtained when the overlapping size $\delta=0$. The observations reflect that the number of GMRES iterations does not necessarily reduce as the overlapping size increase.

5.3 Large-scale scalability

To study the parallel scalability, we run the test case $C$ and $D$ for the first $10$ time steps with different number of processor cores, respectively. In the test case $C$, a $24,576\times 6,144$ mesh is considered and time step size $\Delta t$ is fixed to be $1$. Based on the observations from the above subsection, we use the left-RAS preconditioner with the overlapping $\delta=1$ and employ the sparse ILU(0) factorization with the reuse strategy as the subdomain solver. The numbers of nonlinear and linear iterations are reported in Table 3,

which displays that the total number of nonlinear iterations and the average number of linear iterations keep constant during the increase process of the number of processor cores. Fig. 10 shows the results on the total compute time and the parallel scalability.

It can be seen from Fig. 10 that the total compute time decreases almost linearly, as the number of processor cores increases. The overall speedup from $216$ to $3,456$ cores is around $11.7$, which indicates that the proposed algorithm has a good parallel efficiency for the 2D test case.

We run the test case $D$ on $512\times 512\times 512$ mesh with a fixed time step size $\Delta t=0.1$. We use the classical-AS preconditioner with overlapping size $\delta=0$ and employ the sparse ILU(0) factorization with reuse strategy as the subdomain solver. The numbers of nonlinear and linear iterations are reported in Table 4.

We notices that the number of nonlinear iterations and the average number of linear iterations are almost unchanged during the increase of the number of processors, which implies that the total number of nonlinear iterations and the average number of linear iterations are insensitive to the number of processor cores. Fig. 11 shows the results on the total compute time and the parallel scalability.

The total compute time decreases almost linearly as the number of processor cores increases. The overall speedup from $1,536$ to $24,576$ cores is around $13.9$, which indicates an almost ideal parallel efficiency of the proposed algorithm for the 3D test case.

6 Conclusion

In this paper, a semi-implicit finite difference scheme and a highly parallel domain decomposition algorithm are proposed for solving the PFC equation. The semi-implicit finite difference scheme is derived based on the DVD method and is proved to be unconditionally stable and satisfies the second order accuracy in time and space. For the steady state calculation, an adaptive time stepping strategy is successfully incorporated into the semi-implicit time integration scheme so that the time step size is controlled based on the state of solution. The nonlinear system constructed by the discretization of PFC equation at each time step is solved by the NKS method with modified boundary conditions for the subdomain solves. The accuracy and applicability of the proposed method are validated by several two and three dimensional test cases. The performance of the NKS method is tuned by changing the subdomain solver, the type of the Schwarz preconditioner and the overlapping size. Large scale numerical experiments show that the proposed algorithm can scale well to over ten thousands processor cores on the Sunway TaihuLight supercomputer.

Acknowledgments

This work was supported in part by Natural Science Foundation of China (grant# 91530323, 11501554), National Key R&D Plan of China (grant# 2016YFB0200603), and Key Research Program of Frontier Sciences from CAS (grant# QYZDB-SSW-SYS006).