An efficient numerical framework for the amplitude expansion of the phase-field crystal model

TL;DR

This paper introduces an efficient finite element-based numerical framework for the amplitude expansion of the phase-field crystal model, enabling large-scale 3D simulations of polycrystalline microstructures with dislocation networks.

Contribution

It develops a real space finite element approach with an optimized preconditioner and mesh adaptivity for the APFC model, enhancing simulation efficiency and scale.

Findings

Enables simulation of large 3D polycrystalline systems

Captures microstructure evolution and dislocation dynamics

Improves convergence with a new preconditioner

Abstract

The study of polycrystalline materials requires theoretical and computational techniques enabling multiscale investigations. The amplitude expansion of the phase field crystal model (APFC) allows for describing crystal lattice properties on diffusive timescales by focusing on continuous fields varying on length scales larger than the atomic spacing. Thus, it allows for the simulation of large systems still retaining details of the crystal lattice. Fostered by the applications of this approach, we present here an efficient numerical framework to solve its equations. In particular, we consider a real space approach exploiting the finite element method. An optimized preconditioner is developed in order to improve the convergence of the linear solver. Moreover, a mesh adaptivity criterion based on the local rotation of the polycrystal is used. This results in an unprecedented capability of simulating large, three-dimensional systems including the dynamical description of the microstructures in polycrystalline materials together with their dislocation networks.

Equations

$$n(\mathbf{r})=n_0+\sum _{j=1}^J\left[{\eta }_j(\mathbf{r})e^{\text{i}{\mathbf{k}}_j\cdot \mathbf{r}}+{\eta }_j^{\ast }(\mathbf{r})e^{-\text{i}{\mathbf{k}}_j\cdot \mathbf{r}}\right],$$

$$F(\vec{\mathbf{\eta }}):={\int }_{\Omega }{f}^{\text{s}}(\vec{\mathbf{\eta }},{\vec{\mathbf{\eta }}}^{\ast })+\frac{1}{2}\sum _{j=1}^J|{\mathcal{G}}_j{\eta }_j{|}^2\text{d}\mathbf{r},$$

$$\frac{\partial {\eta }_j}{\partial t}=-{\kappa }_j\frac{\delta F}{\delta {\eta }_j^{\ast }},j=1,\dots ,J,$$

$$\frac{\delta F}{\delta {\eta }_j^{\ast }}={\mathcal{G}}_j^2{\eta }_j+\frac{\delta f^{\text{s}}(\vec{\mathbf{\eta }},{\vec{\mathbf{\eta }}}^{\ast })}{\delta {\eta }_j^{\ast }}.$$

$$f^{\text{s}}(\vec{\mathbf{\eta }},{\vec{\mathbf{\eta }}}^{\ast })=\frac{a_1}{a_0}A+\frac{a_2}{a_0}\left(A^2-2\sum _{j=1}^J|{\eta }_j{|}^4\right)+g^{\text{s}}(\vec{\mathbf{\eta }},{\vec{\mathbf{\eta }}}^{\ast }),$$

$$g^{\text{s}}(\vec{\mathbf{\eta }},{\vec{\mathbf{\eta }}}^{\ast })=-2\frac{a_3}{a_0}\left({\eta }_1{\eta }_2{\eta }_3+{\eta }_1^{\ast }{\eta }_2^{\ast }{\eta }_3^{\ast }\right),$$

$${\mathbf{k}}_1=\left(-\frac{\sqrt{3}}{2},-\frac{1}{2}\right),{\mathbf{k}}_2=(0,1),{\mathbf{k}}_3=\left(\frac{\sqrt{3}}{2},-\frac{1}{2}\right).$$

$$\begin{split}(\zeta_{j},\,\vartheta_{1})+(\nabla\eta_{j},\,\nabla\vartheta_{1})-2\textrm{i}(\mathbf{k}_{j}\cdot\nabla\eta_{j},\,\vartheta_{1})&=0\\ \left(\frac{\partial\eta_{j}}{\partial t},\,\vartheta_{2}\right)-\kappa_{j}\big{(}(\nabla\zeta_{j},\,\nabla\vartheta_{2})-2\textrm{i}(\mathbf{k}_{j}\cdot\nabla\zeta_{j},\,\vartheta_{2})\big{)}&=\left(\frac{\delta f^{\textrm{s}}(\vec{\boldsymbol{\eta}},\vec{\boldsymbol{\eta}}^{*})}{\delta\eta_{j}^{*}},\,\vartheta_{2}\right)\,,\qquad\forall\vartheta_{1},\vartheta_{2}\in H^{1}(\Omega,\mathbb{C})\end{split}$$

$$\begin{array}{rl}{\mathcal{V}}_h& =\{\vartheta \in C(\Omega ,\mathbb{C}):\vartheta {|}_S\in {\mathbb{P}}_1(S,\mathbb{C}),S\in {\mathcal{T}}_h\},\end{array}$$

$$\begin{split}(\zeta^{h}_{j},\,\vartheta^{h}_{1})+(\nabla\eta^{h}_{j},\,\nabla\vartheta^{h}_{1})-2\textrm{i}(\mathbf{k}_{j}\cdot\nabla\eta^{h}_{j},\,\vartheta^{h}_{1})&=0\\ \left(\frac{\partial\eta^{h}_{j}}{\partial t},\,\vartheta^{h}_{2}\right)-\kappa_{j}\big{(}(\nabla\zeta^{h}_{j},\,\nabla\vartheta^{h}_{2})-2\textrm{i}(\mathbf{k}_{j}\cdot\nabla\zeta^{h}_{j},\,\vartheta^{h}_{2})\big{)}&=\left(\frac{\delta f^{\textrm{s}}(\vec{\boldsymbol{\eta}}^{h},\vec{\boldsymbol{\eta}}^{h*})}{\delta\eta_{j}^{*}},\,\vartheta^{h}_{2}\right)\,,\qquad\forall\vartheta^{h}_{1},\vartheta^{h}_{2}\in\mathcal{V}_{h}\,,\end{split}$$

$$d{F}_j^{\text{s}}(\vec{\mathbf{\eta }}):=\frac{\delta F_j^{\text{s}}}{\delta {\eta }_j}(\vec{\mathbf{\eta }})=\frac{{\delta }^2{f}^s}{\delta {\eta }_j\delta {\eta }_j^{\ast }}(\vec{\mathbf{\eta }}).$$

$$\left[\begin{array}{cc}\mathbf{M}& \mathbf{K}-\text{i}{\mathbf{T}}^j\\ -{\kappa }_j(\mathbf{K}-\text{i}{\mathbf{T}}^j)& \frac{1}{{\tau }_m}\mathbf{M}+{\mathbf{J}}^j({\vec{\mathbf{\eta }}}^{(k-1)})\end{array}\right]\left(\begin{array}{c}{\mathbf{\zeta }}_j^{(k)}\\ {\mathbf{\eta }}_j^{(k)}\end{array}\right)=\left(\begin{array}{c}0\\ \frac{1}{{\tau }_m}\mathbf{M}{\mathbf{\eta }}_j^{m-1}+{\mathbf{h}}^j({\vec{\mathbf{\eta }}}^{(k-1)})\end{array}\right).$$

$$\underset{\mathbf{A}}{\underbrace{\left[\begin{array}{cc}\mathbf{M}& -\mathbf{G}\\ \kappa \mathbf{G}& \frac{1}{\tau }\mathbf{M}+\mathbf{J}\end{array}\right]}}\underset{\mathbf{x}}{\underbrace{\left(\begin{array}{c}\mathbf{\zeta }\\ \mathbf{\eta }\end{array}\right)}}=\left(\begin{array}{c}0\\ \frac{1}{\tau }\mathbf{M}{\mathbf{\eta }}^m+\mathbf{h}\end{array}\right)$$

$$\mathbf{P}:=\left[\begin{array}{cc}\mathbf{M}& 0\\ \kappa \mathbf{G}& \frac{1}{\sqrt{\tau }}\mathbf{M}+\sqrt{\kappa }\mathbf{K}\end{array}\right]\left[\begin{array}{cc}\mathbf{I}& -{\mathbf{M}}^{-1}\mathbf{K}\\ 0& {\mathbf{M}}^{-1}(\frac{1}{\sqrt{\tau }}\mathbf{M}+\sqrt{\kappa }\mathbf{K})\end{array}\right]=\left[\begin{array}{cc}\mathbf{M}& -\mathbf{K}\\ \kappa \mathbf{G}& \frac{1}{\tau }\mathbf{M}+2\sqrt{\frac{\kappa }{\tau }}\mathbf{K}+\kappa (\mathbf{K}{\mathbf{M}}^{-1}\mathbf{K}-\mathbf{G}{\mathbf{M}}^{-1}\mathbf{K})\end{array}\right]$$

$$\hat{\eta }(\mathbf{q}):=\mathcal{F}\eta ={\int }_{\Omega }\eta (\mathbf{r})e^{-i\mathbf{q}\cdot \mathbf{r}}\text{d}\mathbf{r},$$

$$\mathcal{A}=\left[\begin{array}{cc}1& {\mathbf{q}}^2+2\mathbf{q}\cdot \mathbf{k}\\ -\kappa ({\mathbf{q}}^2+2\mathbf{q}\cdot \mathbf{k})& \frac{1}{\tau }\end{array}\right],\mathcal{P}=\left[\begin{array}{cc}1& 0\\ -\kappa ({\mathbf{q}}^2+2\mathbf{q}\cdot \mathbf{k})& \frac{1}{\sqrt{\tau }}+\sqrt{\kappa }{\mathbf{q}}^2\end{array}\right]\left[\begin{array}{cc}1& -{\mathbf{q}}^2\\ 0& \frac{1}{\sqrt{\tau }}+\sqrt{\kappa }{\mathbf{q}}^2\end{array}\right].$$

$$\sigma(\mathcal{Q})=\left\{1,\lambda(\mathbf{q})\coloneqq\frac{\mathbf{q}^{4}\mp 4\mathbf{q}^{3}+4\mathbf{q}^{2}+1}{\mathbf{q}^{4}+2\mathbf{q}^{2}+1}\;\big{|}\;\mathbf{q}\in\mathbb{R}\right\}\,,$$

$$\rho :=\frac{\sqrt{{\lambda }_{\text{max}}/{\lambda }_{\text{min}}}-1}{\sqrt{{\lambda }_{\text{max}}/{\lambda }_{\text{min}}}+1}=0.796,$$

$${\eta }_j={\varphi }_j{e}^{i\delta {\mathbf{k}}_j(\theta )\cdot \mathbf{r}},$$

$${\omega }_{ij}=\frac{1}{2}\left[\frac{\partial u_i}{\partial x_j}-\frac{\partial u_j}{\partial x_i}\right],$$

$${\mathbf{k}}_j\cdot \mathbf{u}=\arctan \left[\frac{\text{Im}({\eta }_j)}{\text{Re}({\eta }_j)}\right].$$

$$h_{\mathrm{amp}}=\frac{1}{n}\underset{j}{\min }\left(\left|{\lambda }_j[\underset{d}{\max }({\omega }_d)]\right|\right).$$

$$h=\{\begin{array}{ll}{h}_{\mathrm{int}},& \text{mbox}if|\nabla A|\ge ϵ\\ \min (\max (h_{\mathrm{amp}},h_{\mathrm{int}}),h_{\max }),& \text{mbox}ifA>0\text{mbox}and|\nabla A|<ϵ\\ h_{\max },& \text{mbox}elsewhere.\end{array}$$

Full text

I Introduction

The study of polycrystalline materials requires modeling over various length scales. While the nucleation and structure of defects has to be considered on an atomistic resolution, the size, shapes and arrangements of grains can be considered on mesoscopic length scales and mechanical behavior can be tackled on a macroscopic scale. However, all these phenomena strongly affect material properties, such as fracture or yield stress and thus require multiscale investigations Rollett et al. (2015).

Coarse-grained, mesoscale approaches are therefore viable tools to provide bridging-scale information of polycrystalline materials. Among the many different models reported in the literature, the so-called phase-field crystal (PFC) model has been developed in order to filter out atom vibrations and enabling the investigation of relatively long dynamics Elder et al. (2002); Elder and Grant (2004); Emmerich et al. (2012). It describes the crystal lattice by means of a continuous field that is the atomic probability density, $n$, averaged over vibrational length scales. Although the original model has been further extended, as e.g. in Refs Greenwood et al. (2010); Huang et al. (2010); Kocher and Provatas (2015), and applied to describe several different properties and mechanisms for crystals and quasicrystals, as e.g. in Refs. Rottler et al. (2012); Berry et al. (2014); Achim et al. (2014); Hirvonen et al. (2016); Wang et al. (2016); Yamanaka et al. (2017), it is restricted to relatively small sizes as the continuous density has still to be resolved at the atomic length scale.

This limit has been overcome by the so-called amplitude expansion of the PFC model (APFC), providing coarse-graining in both time and space in a single framework Goldenfeld et al. (2005); Athreya et al. (2006); Goldenfeld et al. (2006); Yeon et al. (2010). This model focuses on the amplitudes of the atomic probability density which vary on a larger length scale than the atomic spacing. In its standard formulation, it is restricted to relatively small deformations Spatschek and Karma (2010) and approximate atomic rearrangements such as dislocation-core structures. However, it encodes a detailed description of deformations matching continuum-elasticity theories and allows for large simulations approaching the macroscopic length scales, still retaining details of the atomic length scale Salvalaglio et al. (2019). The original model can be further extended as it has been done, for instance, to account for binary systems Elder et al. (2010a) and to control the energy of defects and interfaces Salvalaglio et al. (2017). Moreover, it has been recently used to study the anisotropic shrinkage in 3D of small-angle spherical grain boundaries (GBs) regardless of the crystal lattice symmetry Salvalaglio et al. (2018). In addition, the APFC framework has been proven suitable to allow for the description of hydrodynamics Heinonen et al. (2016), dislocation dynamics Skaugen et al. (2018a) and surface-energy anisotropy Ofori-Opoku et al. (2018) within the more general PFC framework.

In this work, we report on the development of an efficient numerical framework to integrate the partial differential equations of the APFC model, enabling unprecendented simulations of large polycrystalline systems in two- and three-dimensions. We provide an optimized discretization of the APFC equations based on the one reported in Ref. Salvalaglio et al. (2017). Then, we propose a new preconditioner allowing for faster convergence of iterative linear solvers. We also illustrate a mesh-adaptivity strategy exploiting continuous fields, namely local crystal rotations, which can be derived directly by the complex amplitudes to solve for in the APFC model. Performance studies showing comparisons with a standard solver and assessment of the numerical parameter entering the method are reported. Then, we use the developed framework to simulate the large scale growth of a polycrystal in 2D and 3D accounting for both the evolution of the microstructures and the defects forming between grains having different orientations.

II Model and Implementation

The APFC model Goldenfeld et al. (2005); Athreya et al. (2006); Goldenfeld et al. (2006); Yeon et al. (2010) is based on the representation of the atomic probability density, $n$, as sum of plane waves

[TABLE]

with $n_{0}$ the average density, here set to zero without loss of generality, $\eta_{j}(\mathbf{r})$ the amplitude of each plane wave, and $\mathbf{k}_{j}$ the reciprocal lattice vector representing a specific crystal symmetry. Complex amplitude functions $\eta_{j}:\Omega\rightarrow\mathbb{C}$ with $\Omega$ a rectangular (2D) or parallelepiped (3D) domain are considered, allowing for describing distortions and rotations of the crystal lattice with respect to a reference state described by the $\mathbf{k}_{j}$ vectors. The model is based on the definition of a free energy in terms of the $\eta_{j}$’s in the approximation of slowly varying amplitudes, i.e., varying over a length scale significantly larger than the lattice spacing,

[TABLE]

with $\vec{\boldsymbol{\eta}}=\{\eta_{j}\}$, $\mathcal{G}_{j}\coloneqq\Delta+2\textrm{i}\mathbf{k}_{j}\cdot\nabla$, and $f^{\textrm{s}}$ a polynomial bulk free energy density. For the sake of convenience, the energy is here normalized by a positive, non-zero parameter $a_{0}$ which would multiply the second term at the right-hand side of (2) Elder et al. (2010b); Salvalaglio et al. (2017). The evolution laws for amplitudes $\eta_{j}$ are given by the $L^{2}$-gradient flow of the energy $F$,

[TABLE]

with $\kappa_{j}=a_{0}|\mathbf{k}_{j}|^{2}$ a mobility coefficient and

[TABLE]

The bulk energy term $f^{\textrm{s}}(\vec{\boldsymbol{\eta}},\vec{\boldsymbol{\eta}}^{*})$ entering Eq. (2) can be written as

[TABLE]

with $a_{i}$ positive parameters, $g^{\textrm{s}}(\vec{\boldsymbol{\eta}},\vec{\boldsymbol{\eta}}^{*})$ a polynomial in $\vec{\boldsymbol{\eta}}$ and $\vec{\boldsymbol{\eta}}^{*}$ encoding the lattice symmetry, and $A=2\sum_{j=1}^{J}|\eta_{j}|^{2}$. In this work we consider a triangular lattice in 2D and a face-centered cubic (FCC) lattice in 3D. The former is obtained by setting $J=3$, $g^{\rm s}$ entering Eq. (4) as

[TABLE]

and the set of $\mathbf{k}_{j}$ vectors as

[TABLE]

For the FCC lattice $J=7$, thus leading to a larger number of amplitude functions to be considered. For the sake of brevity we refer to Ref. Salvalaglio et al. (2017) for the corresponding $g^{\rm s}$ polynomial and $\mathbf{k}_{j}$ vectors. Following Refs. Elder et al. (2010b); Salvalaglio et al. (2017, 2018), the parameters entering the energy are set to favor the growth of the solid phase as $a_{0}=0.98$, $a_{1}=0.01$, $a_{2}=0.25$, $a_{3}=0.5$.

II.1 Discretization and Nonlinear solver

Throughout the following, we use the $L^{2}(\Omega,\mathbb{C})$ scalar product $(u,\,v)\coloneqq\int_{\Omega}u^{*}(\mathbf{r})\cdot v(\mathbf{r})\,\text{d}\mathbf{r}$, with $u,v$ complex scalar or vector valued functions and $u^{*}$ the complex conjugate of $u$. By introducing a set of auxiliary variables $\zeta_{j}\coloneqq\mathcal{G}_{j}\eta_{j}$, we derive the weak form of the evolution equations (3) using a natural splitting of the operator $\mathcal{G}_{j}^{2}$:

[TABLE]

This set of coupled second order equations can then be discretized using conforming finite-dimensional approximations of the functions space $H^{1}(\Omega,\mathbb{C})$. Therefore, let us now consider a triangulation of the domain $\Omega$, $\mathcal{T}_{h}$, with $h$ the minimal diameter of the grid elements, and the finite-element space $\mathcal{V}_{h}\subset H^{1}(\Omega,\mathbb{C})$,

[TABLE]

the continuous space of complex valued local linear polynomials on elements $S\in\mathcal{T}_{h}$ of the triangulation. The space discretization method of (8) corresponds to find $\eta^{h}_{j},\zeta^{h}_{j}\in L^{1}([0,T],\mathcal{V}_{h})$ such that

[TABLE]

for $j=1,\ldots,J$ in the time interval $[0,T]$ provided that $\eta_{j}^{h}(t=0)=\eta_{j}^{h,0}$ is given. For better readability, we drop the superscript $h$ in the following.

Due to the bulk polynomial $f^{\textrm{s}}$, equations (9) are a set of nonlinear equations. We combine a simplified Newton method with the idea of an operator splitting method, in order to separate the equation for each amplitude. This will allow us to efficiently solve the set of equations also if the number of amplitudes is large, as in the case of FCC symmetry. We introduce the Jacobian of $F^{\textrm{s}}_{j}\coloneqq\delta f^{\textrm{s}}/\delta\eta_{j}^{*}$, but only in direction $\eta_{j}$, i.e.,

[TABLE]

Let $\{\phi_{i}\}$ be a basis of $\mathcal{V}_{h}$. Moreover, let us introduce the mass matrix \mathbf{M}\coloneqq\big{[}{(\phi_{i},\,\phi_{j})}\big{]}, the stiffness matrix \mathbf{K}\coloneqq\big{[}{(\nabla\phi_{i},\nabla\phi_{j})}\big{]}, \mathbf{T}^{k}\coloneqq\big{[}{(2\mathbf{k}_{k}\cdot\nabla\phi_{i},\phi_{j})}\big{]}, \mathbf{J}^{k}(\vec{\boldsymbol{\eta}})\coloneqq\big{[}{(-dF^{\textrm{s}}_{k}(\vec{\boldsymbol{\eta}})\phi_{i},\,\phi_{j})}\big{]}, and the vector \mathbf{h}^{k}(\vec{\boldsymbol{\eta}})\coloneqq\big{[}{(F^{\textrm{s}}_{k}(\vec{\boldsymbol{\eta}})-dF^{\textrm{s}}_{k}(\vec{\boldsymbol{\eta}})\eta_{k},\,\phi_{j})}\big{]}. Note, that these matrices are assembled with the real valued basis functions $\phi_{i}$.

For the time discretization we choose a backward Euler method. In the following let $0=t_{0}<t_{1}<\ldots<t_{M}=T$ be a discretization of the time interval $[0,T]$ with $\tau_{m}:=t_{m}-t_{m-1}$ and the timestep solution $\eta^{m}_{j}=\sum_{i}\eta^{m}_{j,i}\phi_{i}\equiv\eta_{j}(t_{m})$ with \boldsymbol{\eta}_{j}^{m}=\big{[}{\eta^{m}_{j,i}}\big{]} the vector of the coefficients for the $j$-th amplitude at the timestep $m$. Analogously, we define $\zeta^{m}_{j}=\sum_{i}\zeta^{m}_{j,i}\phi_{i}\equiv\zeta_{j}(t_{m})$ and \boldsymbol{\zeta}_{j}^{m}=\big{[}{\zeta^{m}_{j,i}}\big{]}.

The timestep iteration with an inner Picard iterative process for the linearized $F^{\textrm{s}}_{j}$ can be found in Alg. missing 1. The discretization is implemented in the finite-element framework AMDiS Vey and Voigt (2007); Witkowski et al. (2015).

II.2 Solving the Linear System

The skew symmetric linear system (10) needs to be solved in each timestep and in each Newton iteration. For large-scale numerical setups in 3D, the widely used direct spare LU factorizations are out of applicability due to large memory requirements. Therefore, we concentrate on preconditioned iterative methods, belonging in particular to the Krylov subspace, to solve (10). Two strategies will be compared: domain decomposition with local sparse LU factorization also known as block Jacobi (bjacobi) preconditioner Saad (2003) and a dedicated preconditioner based on the schur-complement method (apfc), as detailed in the following.

For the construction of the schur-complement preconditioner, we follow the lines of Boyanova et al. (2012); Axelsson and Neytcheva (2013), and first simplify the notation by neglecting the amplitude index $j$. For each amplitude a structurally similar system has to be solved, so we can construct one preconditioner that can be applied to all amplitudes, with adapted coefficients. Also, we drop the timestep index $m$ and Newton iterate index $(k)$ and just write $\tau\equiv\tau_{m}$. In spite of the operator $\mathcal{G}_{j}$ we introduce $\mathbf{G}^{j}\coloneqq-\mathbf{K}+\textrm{i}\mathbf{T}^{j}$ and write the linear system as

[TABLE]

with schur-complement $\mathbf{S}=\frac{1}{\tau}\mathbf{M}+\mathbf{J}+\kappa\mathbf{G}\mathbf{M}^{-1}\mathbf{G}$.

While $\mathbf{M}$ and $\mathbf{K}$ are real matrices, $\mathbf{G}$, $\mathbf{J}$, and $\mathbf{h}$ may be complex valued. Following the ideas of Boyanova et al. (2012); Axelsson and Neytcheva (2013) to approximate the schur-complement matrix $\mathbf{S}$ in a way that allows for a simple factorization, we propose the preconditioner:

[TABLE]

with schur complement $\mathbf{S}_{P}\coloneqq(\frac{1}{\sqrt{\tau}}\mathbf{M}+\sqrt{\kappa}\mathbf{K})\mathbf{M}^{-1}(\frac{1}{\sqrt{\tau}}\mathbf{M}+\sqrt{\kappa}\mathbf{K})$. Since the highest order term in $\mathbf{G}\mathbf{M}^{-1}\mathbf{G}$ corresponds to $\mathbf{K}\mathbf{M}^{-1}\mathbf{K}$, the matrix $\mathbf{S}_{P}$ is expected to be a good approximation of $\mathbf{S}$.

The application of $\mathbf{P}$ to a vector $\mathbf{b}=(\boldsymbol{\beta}_{1},\boldsymbol{\beta}_{2})^{\top}\in\mathbb{C}^{2\cdot\operatorname{dim}(V_{h})}$, i.e., $\mathbf{x}=\mathbf{P}_{2}^{-1}\mathbf{b}$, with $\mathbf{x}=(\boldsymbol{\chi}_{1},\boldsymbol{\chi}_{2})^{\top}$, can be performed in four steps:

solve $\mathbf{M}\boldsymbol{\beta}^{\prime}_{1}=\boldsymbol{\beta}_{1}$ 2. 2.

solve $(\frac{1}{\sqrt{\tau}}\mathbf{M}+\sqrt{\kappa}\mathbf{K})\boldsymbol{\beta}^{\prime}_{2}=\boldsymbol{\beta}_{2}-\kappa\mathbf{G}\boldsymbol{\beta}^{\prime}_{1}$ 3. 3.

solve $(\frac{1}{\sqrt{\tau}}\mathbf{M}+\sqrt{\kappa}\mathbf{K})\boldsymbol{\chi}_{2}=\mathbf{M}\boldsymbol{\beta}^{\prime}_{2}$ 4. 4.

$\boldsymbol{\chi}_{1}=\boldsymbol{\beta}^{\prime}_{1}-\frac{1}{\sqrt{\kappa}}(\boldsymbol{\beta}^{\prime}_{2}-\frac{1}{\sqrt{\tau}}\boldsymbol{\chi}_{2})$.

Thus, only one mass-matrix and two diffusion-matrix systems need to be solved, plus two matrix-vector products need to be performed in the application of the preconditioner. Since efficient (iterative) linear solvers exist to approximate the solution of $\mathbf{M}^{-1}\mathbf{b}$ and $(\frac{1}{\sqrt{\tau}}\mathbf{M}+\sqrt{\kappa}\mathbf{K})^{-1}\mathbf{b}$, the preconditioner is cheap to apply.

In order to analyze the quality of the preconditioner $\mathbf{P}$, we look at the eigenvalue spectrum of the right-preconditioned linear system $\mathbf{A}\mathbf{P}^{-1}$. Instead of evaluating the eigenvalues directly, we consider the symbols of the differential operators and the resulting spectrum for a simplified system. Therefore, let $\Omega$ be a periodic one-dimensional domain of length $2\pi$. We set the nonlinear term to $\mathbf{J}=0$, and consider only lattice vectors $\mathbf{k}\in\{-1,+1\}$. Let $\mathcal{F}$ denote the non-unitary Fourier transform in angular frequency $\mathbf{q}$, i.e.,

[TABLE]

with inverse transform $\mathcal{F}^{-1}$ and the linear operators $\mathbf{A}$ and $\mathbf{P}$ expressed as $\mathbf{A}\mathbf{x}=\mathcal{F}^{-1}(\mathcal{A}\widehat{\mathbf{x}})$ and $\mathbf{P}\mathbf{x}=\mathcal{F}^{-1}(\mathcal{P}\widehat{\mathbf{x}})$, respectively, with symbols $\mathcal{A}$ and $\mathcal{P}$. In terms of the frequency $\mathbf{q}$, these symbols can be written as

[TABLE]

For $\mathcal{Q}:=\mathcal{A}\mathcal{P}^{-1}$, $\kappa=1$, and $\tau=1$, we obtain the spectrum

[TABLE]

where the sign of the second term in the nominator depends on the choice of the k-vector, see Fig. missing 1 for a visualization. The eigenfunction $\lambda$ approaches the value 1 for zero and $\pm\infty$ frequencies and has a minimum $\min\{\lambda(\mathbf{q})\,|\,\mathbf{q}\in\mathbb{R}\}\approx 0.035$ and maximum $\max\{\lambda(\mathbf{q})\,|\,\mathbf{q}\in\mathbb{R}\}\approx 2.70$. The operator symbol $\mathcal{A}$ forms a normal matrix for $|\mathbf{k}_{j}|=1$ and the preconditioned operator spectrum is in the positive and real interval $[0.035,2.7]$, leading to an asymptotic convergence factor (see Ref. Saad (2003))

[TABLE]

and thus rapid convergence of an iterative Krylov subspace method, like the Flexible Generalized Minimal Residual Method Saad (1993) (known as FGMRES). See Praetorius and Voigt (2015) for a similar analysis on a preconditioner for the discretized PFC equation.

Note, if the operator $\mathcal{G}^{j}$ contains higher order derivatives, a similar procedure for the construction of the preconditioner can be performed that may lead to a slightly different sequence of mass and diffusion-like systems but follows the same structure. A construction of a preconditioner for a sixth order PFC, an eighth order PFC, and a Lifschitz-Petrich type energy is discussed in Praetorius (2015), eventually leading to a sequence of diffusion equations to solve in the application of the preconditioner. The state potential $f^{s}$ may be extended in different directions and we expect the same procedure to work as above, provided that $f^{s}$ contains no derivative terms and the problem is well defined.

II.3 Mesh adaptivity

The variables to solve for within the APFC model, namely $\eta_{j}$’s, are constant for relaxed crystals, oscillate with different periodicity according to the local distortion of the crystal with respect to the reference one, and exhibit significant variation at defects and solid-liquid interfaces. This allows for exploiting mesh adaptivity in order to optimize the numerical approach Athreya et al. (2007); Berčič and Kugler (2018). We set a local grid refinement in order to resolve the oscillation of amplitudes within grains and ensuring a proper resolution at defects and interfaces.

Amplitude functions describing a rotated crystal of an angle $\theta$ can be computed as

[TABLE]

with $\delta\mathbf{k}_{j}(\theta)=\mathbf{k}_{j}\cdot\mathbf{R}(\theta)-\mathbf{k}_{j}$, $\mathbf{R}(\theta)$ the counterclockwise rotation matrix and $\phi_{j}$ the real amplitudes describing a relaxed, unrotated crystal. Therefore, the wavelength of $\eta_{j}$’s in the presence of a certain rotation $\theta$ is $\lambda_{j}(\theta)=2\pi/|\delta\mathbf{k}_{j}(\theta)|$ and the spatial discretization is set as a fraction of the smallest $\lambda_{j}(\theta)$ i.e. h=$\min_{j}{(|\lambda_{j}(\theta)|)}/n$. We determined that a proper discretization is achieved with $n\geq 10$. In order to practically use this idea we then need to evaluate the local rotation field during the evolution. Notice that Eq. (13) cannot be just inverted due to its functional form. In order to compute the local rotation we then use the approach described in Salvalaglio et al. (2019) where the rotational tensor $\boldsymbol{\omega}$ is determined by considering the curl of the local displacement with respect to the relaxed crystal. In practice, its components $\omega_{ij}$ representing the rotation in the $x_{i}$-$x_{j}$-plane, are given by

[TABLE]

with $u_{i}$ the result of inverting the system of equations

[TABLE]

Explicit expressions for $u_{i}$, both in 2D and 3D can be found in Ref. Salvalaglio et al. (2019). In 2D, Eq. (14) delivers just a scalar field $\omega_{ij}\equiv\omega$. In 3D, $\boldsymbol{\omega}$ has three independent components, $\omega_{d}$ with $d=1,2,3$, and the largest rotation component is chosen in order to ensure the proper resolution for the fastest oscillation in the system. Within grains the discretization, $h_{\rm amp}$, then reads

[TABLE]

The evaluation of rotations considered here is well posed in the solid phase and also at defects Salvalaglio et al. (2019). However, it is not well posed for the liquid, disordered phase as $\eta_{j}$’s vanish. Therefore, in order to properly discretize interface regions we define an additional refinement for the interfaces controlled by $h_{\rm int}$ where $|\nabla A|$ is significantly larger than a relatively small threshold $\epsilon$ as done in Salvalaglio et al. (2017). As a result, $h_{\rm int}$ is ensured also at defects and we set it as the smallest resolution imposed in the system. In addition, a large discretization bound $h_{\rm max}$ is defined for region where $A\sim 0$ or where the local rotation vanishes. Summarizing these ideas, the local discretization, $h$, is set as

[TABLE]

Although not addressed here, this criterion can be extended in order to account for amplitude oscillations due to strain fields, exploiting continuous strain-field components, as derived in Ref. Salvalaglio et al. (2019), instead of $\boldsymbol{\omega}$.

The model and the mesh discretization are illustrated in Fig. missing 2 where a simulation reproducing the growth of 20 crystal seeds in 2D with $\theta\in(-15^{\circ},15^{\circ})$ is illustrated. A square simulation domain with a side length of $200\pi$ is considered (more details are reported in the following section). The initial rotation of the grain is set by initializing the amplitudes by means of Eq. (missing) 13. $A$, $\text{Re}(\eta_{1})$ and $\omega$ at different stages during the evolution over time are shown superposed to the mesh used for the simulation. The refinement of the mesh where $A$ changes can be recognized at the solid-liquid interface and at defects, see in particular Fig. missing 2(a). The variables to solve for are the amplitudes $\eta_{j}$ as illustrated in Fig. missing 2(b) by $\text{Re}(\eta_{1})$. Notice that the correct resolution of these oscillating functions is ensured by the discretization (17), exploiting the local rotations shown in Fig. missing 2(c).

The mesh adaptivity criterion illustrated in Fig. missing 2 leads to larger refined regions and, in principle, to larger computational costs than the one reported in Ref. Athreya et al. (2007) where an efficient mesh refinement strategy exploiting a polar representation of amplitudes has been proposed. However, as discussed therein, a robust regularization scheme for phase equations and special care to treat high-order derivatives were needed. This practically restricts the applicability of that method and prevents its extension to 3D calculations. An improvement of this description has been recently proposed in Ref. Berčič and Kugler (2018), exploiting spatially-dependent rotation of the $\mathbf{k}_{j}$ vectors, although it has been demonstrated only in 2D. The approach we use here allows for an optimized spatial discretization by using the standard APFC model. As it will be illustrated in the following it can be readily exploited to simulate different systems regardless of system dimensionality.

III Performance studies

For the numerical evaluation of the linear solvers and the mesh adaption strategy, we concentrate on two types of setups, a triangular symmetry in 2D (TRI) and a FCC symmetry in 3D. The domain $\Omega$ is of size $[0,s\cdot 10\pi]^{d}$ where $d$ is the space dimension and $s$ a scaling factor. The unit of lenght and time are dimensionless as they enter in the equations of the APFC model reported above. For the setups (TRI, FCC 1, FCC 2, FCC 3), we have chosen $s$ as ($20,7,14,28$), respectively. To test the numerical framework we consider configurations already known and discussed in literature. The 2D triangular setup is initialized with 20 randomly positioned and oriented grains with orientation in the range $(-15^{\circ},15^{\circ})$, corresponding to the configuration for the growth of a polycrystal as shown in Fig. missing 2 (see, e.g., Refs. Athreya et al. (2007); Spatschek and Karma (2010); Salvalaglio et al. (2019)). The 3D, FCC configurations, are initialized with one grain of size $30\pi$, $60\pi$, or $120\pi$, for the setup FCC 1, FCC 2, or FCC 3, respectively, rotated by $10^{\circ}$ about the [111] direction with respect to a surrounding, unrotated crystal (see Ref. Salvalaglio et al. (2018)). The discretization bounds are $h_{\rm int}\approx 2.0$ and $h_{\rm max}\approx 40.0$ in 2D and $h_{\rm int}\approx 3.0$ and $h_{\rm max}\approx 60.0$ in 3D.

In a first comparison, we evaluate the linear-solver performance for the preconditioner $\mathbf{P}$ defined in Eq. (12) (apfc) and a block Jacobi preconditioner (bjacobi) applied to the blocks resulting from an element-wise domain decomposition, with a local sparse direct solver UMFPACK Davis (2004). The domain decomposition is combined with a distributed memory parallelization. In Fig. missing 3 the timings for solving the linear system per timestep are summarized for TRI, FCC 1, and FCC 2 setup. Timings for the two preconditioners are evaluated by averaging over 100 timesteps over all amplitudes. Depending on the size of the setup different improvements of the apfc preconditioner compared to the bjacobi preconditioner can be found, ranging from a factor $10$ for small number of cores and thus large local partitions, to $2$ for large number of cores and thus small local partitions. In 2D the improvement is nearly constant over all subdomain sizes, but in 3D the local linear system size matters for the bjacobi solver significantly. The three setups show a similar behaviour: the time to solve one linear system goes down with increasing number of processors, while it stagnates after some threshold in the number of cores. This is due to the fact that the local domain sizes become too small and additionally the load balancing gets more and more complicated. Increasing the number of cores further would increase also the solver time. Then, communication cost imbalance of the local problems would dominate the solution procedure.

While the bjacobi preconditioner uses a direct solver on each local subdomain, that results in large memory costs, the apfc preconditioner requires the solution of mass matrix and diffusion like equations, using an iterative solver optimized for the specific type of linear system. For the mass matrix, we use three iterations of a diagonally preconditioned conjugate gradient (CG) method. The diffusion system, on the other hand, is solved using either a diagonally preconditioned CG method with 5 iterations (parameter=1), or an algebraic multigrid method (AMG, parameter=2), based on the so-called Hypre BoomerAMG Henson and Yang (2002), with one V-cycle and symmetric relaxation. In Fig. missing 4(a) the difference in the solver time per timestep is plotted, depending on the timestep parameter $\tau$ that is part of the preconditioner definition (12). Increasing the timestep size results in slower convergence of all the linear solvers, but the effect is more pronounced with the apfc preconditioner. In the same figure also the difference in the solver time for different solver parameters is shown. For this setup, a simple CG iteration outperforms the AMG subsolver. Concerning just the performance of the linear solver the optimal timestep would be as large as possible, but the accuracy of the scheme is just first order in $\tau$ so that a balance between accuracy and performance has to be found for a concrete setup.

A third benchmark considers the parallel scalability of the apfc preconditioner. Therefore, we run setup FCC 3 with an increasing number of subdomains, that is, an increasing number of compute cores. These benchmarks are run on the JUWELS cluster of the Jülich Supercomputing Centre on regular $2\times 24$ core nodes. Since an allocation of a different number of running processes per node results in different node performance, e.g., due to parallel usage of memory pipelines, we have assigned a fixed number of cores per node for each scaling benchmark. In Fig. missing 4(b) the relative speedup for an allocation of 24 and 48 cores per node is plotted with respect to the smallest possible number of cores, i.e., one node with 24 or 48 cores, respectively. Using less cores per node leads to a slightly better speedup that might even exceed the ideal scaling line. Although the parallel efficiency drops down if the local problems get too small, in the intermediate range of up to 300 cores, we see very good speedup and thus parallel efficiency of the apfc preconditioner.

IV Application: Growth of polycrystals with dislocations

In this section we illustrate some applications of the numerical approach illustrated above. In order to allow for estimates and comparisons in terms of number of atoms, e.g. to other methods, simulations or real systems, the length scale is here reported in terms of the lattice constant for the considered symmetries as set by the corresponding set of $\mathbf{k}_{j}$ vectors. They read $a_{\rm tri}$ and $a_{\rm fcc}$ for the triangular and FCC lattice symmetry, respectively.

Fig. missing 5 shows the growth of 200 crystal seeds in 2D having triangular symmetry. A rotation of crystals $\theta\in(-15^{\circ},15^{\circ})$ set as initial condition by means of Eq. (missing) 13 as in Fig. missing 2 is considered. We set here a square domain of side length $\sim 10^{3}a_{\rm tri}$. The initial crystal seeds are generated in a square region at the center (see Fig. missing 5(a), $t=0$). As shown by two representative steps during the evolution, $t=2000$ and $t=5500$, the growth of these seeds results in a polycrystal with several dislocations at the center while almost straight GBs form between the peripheral grains which are free to grow towards the liquid phase. The formation of such GBs is also highlighted at $t=12000$ by means of a portion of the entire simulation domain. Insets of Fig. missing 5(a) illustrate the arrangement of defects on a smaller length scale in a portion of the crystal. The motion and eventual annihilation of dislocations is accounted for by the approach as illustrated in Fig. missing 5(b) where the arrangement of defects in the central region of the simulation (green shaded square superposed to the crystal at $t=5500$) is compared at $t=600$ (blue, filled) and $t=12000$ (red, empty). The coarsening dynamics results faster in the early stages than at later times as can be noticed by comparing the main features of the dislocation networks at different times in Fig. missing 5. Indeed, small grains are present at the beginning leading to the formation of curved grain boundaries formed by a few dislocations, which, in turn, move fast Mullins (1956); Doherty et al. (1997). Later, the grains at the center of the polycrystal are larger and the resulting GBs are more stable having smaller curvatures. However, the coarsening dynamics with motion and annihilation of defects continues as can be notice in the insets of Fig. missing 5(a), where significant changes in the arrangement of defects are observed at later times. The straight grain boundaries forming due to the growth of the polycrystal can be considered as long-lasting defects as their curvature is negligible, while the spacing between dislocations depends on the relative tilts of grains. Notice that the growth velocity of the polycrystal is a function of the parameters entering the free energy that control the energy difference between the solid and the liquid phase Elder et al. (2007). Therefore, faster or slower grain growth compared to defect motion can be inspected.

The main features of the dynamics obtained by the standard APFC approach considered here qualitatively correspond to predictions of classical theories and simulations Mullins (1956); Doherty et al. (1997). However, a quantitative description can be achieved by accounting for extensions of the PFC and APFC approach, including a proper description of elastic and plastic relaxation during the motion of defects Heinonen et al. (2016); Skaugen et al. (2018b).

A 3D example is shown in Fig. missing 6. Two FCC crystals, having a rotation of $\pm 5^{\circ}$ about the [111] direction are considered. The simulation size is set as in the setup FCC 3, corresponding to a square domain of side length $\sim 80a_{\rm fcc}$. They are also aligned along the [111] direction, therefore their growth is expected to form a (111) twisted GB with a typical hexagonal arrangement of defects Scott and Goodhew (1981); De Hosson and Vitek (1990); Salvalaglio et al. (2018). Fig. missing 6 reports four stages during the evolution in terms of $A>0$ along with the dislocation networks at the resulting planar GB illustrated as insets. They correspond to regions where $0<A<0.8\max_{\Omega}(A)$ at the (111) blue shaded plane. The study of these kind of planar GBs, in particular for what concern the morphology of the emerging dislocation network, are typically accounted for by assuming ideally infinite crystals by means of periodic boundary conditions or even just by 2D approaches. Here, the GB is obtained together with the explicit description of crystal growth which enable more general investigations tackling the simultaneous presence of different GBs with different orientations. A more general case, illustrating the general capability of the approach, is reported in Fig. missing 7. Therein 30 crystals having random rotation about the [111] direction are considered. A simulation domain that is double the size of the setup FCC 3, namely corresponding to a square domain of side length $\sim 160a_{\rm fcc}$, is considered here. Fig. missing 7(a) shows the morphologies of the seeds and of the resulting, growing polycrystal in terms of $A>0$ regions. Fig. missing 7(b) shows half of the crystals reported in Fig. missing 7(a) by the isosurface $A=0.8\max_{\Omega}(A)$ revealing also the dislocation network forming at the internal GBs. In this case the initial seeds are distributed randomly and are not aligned along a specific direction. Therefore, together with the twist GBs as shown in Fig. missing 6, other orientations for the boundaries between grains are present thus leading to different morphologies for the dislocation networks. This can be observed in more detail in Fig. missing 7(b)–(d). Hexagonal patterns mostly lying on (111) planes, can be recognized along with elongated defects typical of pure tilt GBs. Moreover, similar patterns having different spacing between dislocations are present due to different relative rotations between grains.

The simulations reported in this section are unprecedented in terms of sizes for what concern APFC (and PFC) approaches. Just by focusing on the largest system of Fig. 7, we were able to simulate a crystalline system with FCC lattice symmetry including $\sim 8\cdot 10^{6}$ atoms. For specific materials exhibiting such a lattice symmetry as, e.g, Cu, Ag/Au, Pb this would mean a volume of $\sim(45\ {\rm nm})^{3}$, $\sim(50\ {\rm nm})^{3}$, $\sim(65\ {\rm nm})^{3}$ respectively Davey (1925), which lie in the typical size range of nanoparticles and nanostructures. The simulation in Fig. 7 is done on 720 cores. Increasing the number of cores will also allow to consider even larger samples and thus realize the envisioned multiscale approach, ranging from atomistic details to micrometer sizes.

V Conclusions

We illustrated a numerical approach to solve the equations of the APFC model efficiently. A specific discretization scheme combined with a nonlinear solver has been proved to allow for unprecedented size and performances for both 2D and 3D APFC simulations. In particular, we have constructed a schur-complement preconditioner for iterative Krylov-subspace methods that outperforms a block Jacobi solver for the linear systems arising from the discretization with adaptive finite elements. The schur-complement solver requires much less memory, converges fast in terms of wall-clock time, and scales well in parallel setups. On the other hand, it is more sensitive to an increase in timestep size, compared to the bjacobi solver.

Moreover, an optimized criterion for mesh adaptivity has been proposed and used, exploiting features of the complex amplitude functions $\eta_{j}$ and derived physical quantities such as the local rotation field. This can be applied in both 2D and 3D and allows for a significant reduction of the overall number of DOFs.

Some applications involving the growth of polycrystals in 2D and 3D as well as the simultaneous description of dislocations forming at GBs have been shown. They set new limits for APFC and then, more in general, PFC approaches, enabling the investigation of large systems matching the size of real nanostructures. This is a crucial step to address large, mesoscale problems still retaining details of the atomic length scale. Future work will be devoted to explicitly include further details, compatible with the APFC model, in the numerical framework presented here as, for instance, an improved description of interface-energy anisotropy Ofori-Opoku et al. (2018); Salvalaglio et al. (2015), binary systems Elder et al. (2010b), an improved description of the dynamics Heinonen et al. (2016); Skaugen et al. (2018b), the coupling with magnetic fields Backofen et al. (2019), as well as improvements on scaling properties of the numerical approach to enable larger systems addressing grain growth in 3D.

VI Acknowledgements

The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this work by providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC), under the grant No. HDR06. A.V. acknowledge the financial support from the German Research Foundation (DFG) under Grant No. SPP 1959.