# Parallel energy stable phase field simulations of Ni-based alloys system
††thanks: This work was supported in part by NSFC (grant # 11871069 and
12131002), and the Strategic Priority Research Program of the Chinese Academy
of Sciences (grant # XDB22020100).

Jizu Huang LSEC, Academy of Mathematics and Systems Science, Chinese Academy
of Sciences, Beijing 100190, China, and School of Mathematical Sciences,
University of Chinese Academy of Sciences, Beijing 100049, China
(huangjz@lsec.cc.ac.cn). Chao Yang Corresponding author: Chao Yang
<EMAIL_ADDRESS>School of Mathematical Sciences, Peking University,
Beijing 100871, China, and Institute for Computing and Digital Economy, Peking
University, Changsha 410205, China. (chao_yang@pku.edu.cn).

###### Abstract

In this paper, we investigate numerical methods for solving Nickel-based phase
field system related to free energy, including the elastic energy and
logarithmic type functionals. To address the challenge posed by the particular
free energy functional, we propose a semi-implicit scheme based on the
discrete variational derivative method, which is unconditionally energy stable
and maintains the energy dissipation law and the mass conservation law. Due to
the good stability of the semi-implicit scheme, the adaptive time step
strategy is adopted, which can flexibly control the time step according to the
dynamic evolution of the problem. A domain decomposition based, parallel
Newton–Krylov–Schwarz method is introduced to solve the nonlinear algebraic
system constructed by the discretization at each time step. Numerical
experiments show that the proposed algorithm is energy stable with large time
steps, and highly scalable to six thousand processor cores.

Key words. Ni-based alloys phase field system, Allen–Cahn/Ginzburg–Landau
equations, discrete variational derivative method, unconditionally energy
stability, Newton–Krylov–Schwarz

AMS Subject Classifications: 74S20, 65Y05

##  1 Introduction

Nickel-based (Ni-based) alloys are often used as high-temperature structural
materials, such as working blades of aero-engine and turbine disk [5, 35] in
the aerospace fields, due to their excellent toughness, high temperature
creeps strength, oxidation resistance, corrosion resistance, and anti-fatigue
properties. The precipitation strengthening of ordered precipitates
($\gamma^{\prime}$ phase) plays an important role in maintaining the high-
temperature mechanical properties of Ni-based alloys, which are embedded in
the matrix ($\gamma$ phase) through the heat treatment process. And the
properties of Ni-based alloys are closely related to the volume fraction of
the ordered precipitates [41]. It is therefore important to study the effect
of volume fraction on the precipitation kinetics of the $\gamma^{\prime}$
phase, including the particle size, number density, and spatial distribution
[36].

The phase field model has been widely used to study microstructure evolutions
in Ni-based alloys [24, 25, 39, 40, 49, 50], such as phase transitions and
dynamic descriptions of the $\gamma^{\prime}$ phase growth and coarsening. In
these phase field models, the $\gamma/\gamma^{\prime}$ microstructure can be
described by a composition field variable and several long-range order
parameters. The total free energy of a Ni-based alloy system, including the
local chemical free energy, the elastic energy, and the interfacial energy, is
defined as a functional of all field variables. The local chemical free energy
is usually constructed by a Landau-type polynomial [39] or using the CALPHAD
method [49]. The elastic energy makes an important contribution to the
morphology evolution and the kinetics of growth and coarsening of the
$\gamma^{\prime}$ phase, which is defined by using the micro-elasticity theory
proposed in [19, 21]. Driven by the total free energy, a gradient flow system
for the Ni-based alloys, such as Cahn–Hilliard (CH) and Ginzburg–Landau (GL)
coupled equations, can be derived.

Over the past two decades, a series of efforts [24, 25, 39, 40, 42, 45, 49,
50] have been made to numerically simulate the morphological evolution of the
Ni-based alloys. Most of these works rely on the semi-implicit Fourier
spectral method [9, 48] to solve the CH and GL coupled equations together with
the linear elasticity equations in order to maintain stability and achieve
high accurate solutions. It is well known that there are a plethora of
efficient solution algorithms for gradient flow problems, including the convex
splitting method [4, 13, 14, 32], the stabilization method [48, 34], the
exponential time differencing [18], the invariant energy quadratization [43,
47], the scalar auxiliary variable [33], the discrete variational derivative
(DVD) method [12, 17, 15, 44], etc. However, due to the existence of the
linear elasticity equations and the logarithmic functions in the free energy,
it is a non-trivial task to extend the above methods to handle the Ni-based
alloys phase field model. On the other hand, the interfacial thickness between
the $\gamma$ and $\gamma^{\prime}$ phase is less than 5 nm [1] in a typical
Ni-based alloy system. To maintain numerical accuracy and stability, the mesh
size should be smaller than the interfacial thickness, which makes the
simulation of Ni-based alloy systems very expensive. As a result, there are
only a few works on three dimensional simulations, usually at the scale of
tens of nanometers. Such a small computational system is usually far from
adequate to study the coarsening process in real three dimensional Ni-based
alloys, where a sufficient number of particles is required for proper
statistical measurements [50].

To overcome the above difficulties, we propose an unconditional energy stable
semi-implicit scheme based on the DVD method to efficiently solve the Ni-based
alloys phase field system, which preserves many properties of the original
phase-field system, for instance, the free energy dissipation, the mass
conservation, etc. Thanks to the unconditional energy stability of the semi-
implicit scheme, an adaptive time stepping strategy [17] can be introduced to
dynamically select the time step size according to the accuracy requirement of
the solutions and the dynamic features of the system, which can improve the
efficiency of the newly proposed scheme. In the semi-implicit method, the most
expensive step is to solve a large sparse nonlinear algebraic system at each
time step. A parallel, highly scalable, Newton– Krylov–Schwarz (NKS) algorithm
[7] is applied to solve such a system, where the nonlinear system arisen at
each time step is solved by an inexact Newton method. And the Jacobian problem
of the NKS algorithm is solved with a preconditioned Krylov subspace method,
within which an LU or ILU factorization is used for solving the subdomain
problem. To achieve the optimal performance of the NKS algorithm, several key
parameters, including the type of the Schwarz preconditioner, the overlap
size, and the subdomain solver, are discussed and tested.

The remainder of this paper is organized as follows. In Sec. 2, we introduce a
phase field model for Ni-based alloys and study the mass conservation law and
energy dissipation law. In Sec. 3, an unconditional energy stable semi-
implicit scheme is proposed based on the DVD method, which is proved to enjoy
the mass conservation law and energy dissipation law. In Sec. 4, we introduce
the NKS algorithm together with the adaptive time step strategy to solve the
nonlinear system. Several three dimensional simulations are reported in Sec. 5
and the paper is concluded in Sec. 6.

##  2 Phase field model for Ni-based alloys

In this section, we introduce a phase field model for the Ni-based alloys to
describe the evolution of the $\gamma/\gamma^{\prime}$ phase driven by total
free energy, which includes the Gibbs energy, the interface energy, and the
elastic energy. The total free energy for the Ni-based alloys has the
following form:

(1) ${\cal
E}=\int\limits_{\Omega}\left[{E}_{\textnormal{G}}(c,\eta)+\frac{\gamma_{c}}{2}|\nabla
c|^{2}+\frac{3\gamma_{\eta}}{2}|\nabla\eta|^{2}+E_{\textnormal{el}}(c)\right]\hbox{d}\mathbf{x},$

where $\Omega\in\mathbb{R}^{d}$ with $d=1,\,2,$ or $3$ is the computational
domain, ${E}_{\textnormal{G}}(c,\eta)$ is the Gibbs energy, and
$E_{\textnormal{el}}(c)$ is the elastic energy. In (1), the conserved
concentration field $c:=c(\mathbf{x},t)\in(0,0.25)$ denotes the composition of
Aluminum (Al) element, $\eta:=\eta(\mathbf{x},t)\in(0,1)$ is a long-range
order parameter field, and $\gamma_{c},\,\gamma_{\eta}$ are the interface
energy coefficients of the composition and long-range order parameter,
respectively. The Gibbs energy ${E}_{\textnormal{G}}(c,\eta)$ is defined as:

(2)
${E}_{\textnormal{G}}(c,\eta)=\Phi(c,\eta)+\theta\left\\{\Psi(c)+\frac{3}{4}\Psi(c\eta)+\frac{1}{4}\Psi(c\eta^{\prime})\right\\},$

where $\Phi(c,\eta)$ is a polynomial functional of $c$ and $\eta$ given in
Appendix A, $\Psi(z)=z\ln z+(1-z)\ln(1-z)$, $\theta$ is a dimensionless
parameter, and $\eta^{\prime}=4-3\eta$. The elastic strain energy density
$E_{\textnormal{el}}(c)$ is given by:

(3)
$E_{\textnormal{el}}(c)=\frac{1}{2}\epsilon^{\textnormal{el}}:C:\epsilon^{\textnormal{el}},$

where $C$ is the Hooke’s tensor and
$\epsilon^{\textnormal{el}}=(\epsilon^{\textnormal{el}}_{IJ})$ is the elastic
strain. The elastic strain can be reformulated as
$\epsilon^{\textnormal{el}}=\epsilon-\epsilon^{0}$, where $\epsilon$ is the
total strain related to the displacement and the eigenstrain $\epsilon^{0}$
expressed as $\epsilon^{0}=\epsilon_{0}(c-\bar{c}){\cal I}$. Here
$\epsilon_{0}$ is a lattice parameter, ${\cal I}$ is the identity matrix, and
$\bar{c}=\frac{1}{|\Omega|}\int_{\Omega}c(\mathbf{x},t)\mathrm{d}\mathbf{x}$
is the average composition. According to the relationship between the strain
and the displacement, the total strain $\epsilon$ can be expressed as:

$\epsilon_{IJ}(\mathbf{x},t)=\frac{1}{2}\left[\frac{\partial u_{I}}{\partial
x_{J}}+\frac{\partial u_{J}}{\partial x_{I}}\right],$

where $\mathbf{u}(\mathbf{x},t):=(u_{1},\cdots,u_{d})^{T}$ denotes the
displacement. Since the mechanical equilibrium for the elastic displacement is
established much faster than composition diffusional process, the system is
always in the mechanical equilibrium

(4) $\nabla\cdot\sigma=\nabla\cdot\sigma^{0},$

where $\sigma=C:\epsilon$ and $\sigma^{0}=C:\epsilon^{0}$.

In the Ni-based alloys, the morphological evolution of the composition field
variable $c(\mathbf{x},t)$ and the long-range order parameter
$\eta(\mathbf{x},t)$ is described by the CH–GL coupled equations in the
following dimensionless form:

(5) $\left\\{\begin{aligned} &\frac{\partial c}{\partial
t}=\nabla\left[M\nabla\frac{\delta{\cal E}}{\delta
c}\right]\quad\hbox{in}~{}\Omega\times(0,\mathcal{T}],\\\
&\frac{\partial\eta}{\partial t}=-\frac{\delta{\cal
E}}{\delta\eta}\quad\hbox{in}~{}\Omega\times(0,\mathcal{T}],\end{aligned}\right.$

where $M=\varkappa c(1-c)$ with $\varkappa$ being a mobility parameter. The
variational derivatives of the free energy functional $\frac{\delta{\cal
E}}{\delta c}$ and $\frac{\delta{\cal E}}{\delta\eta}$ in (5) are defined as:

(6) $\left\\{\begin{aligned} &\frac{\delta{\cal E}}{\delta c}=\frac{\partial
E_{\textnormal{G}}(c,\eta)}{\partial c}-\gamma_{c}\Delta c+\frac{\partial
E_{\textnormal{el}}(c)}{\partial c},\\\ &\frac{\delta{\cal
E}}{\delta\eta}=\frac{\partial
E_{\textnormal{G}}(c,\eta)}{\partial\eta}-3\gamma_{\eta}\Delta\eta.\end{aligned}\right.$

More details about the Ni-based alloys phase field system are discussed in
Appendix A. The initial conditions for CH–GL coupled equations Eq. 5 are given
as $c^{0}=c(\mathbf{x},0)$ and $\eta^{0}=\eta(\mathbf{x},0)$. We consider
periodic boundary conditions or the following homogeneous Neumann boundary
conditions

(7) $\mathbf{n}\cdot\nabla
c|_{\partial\Omega}=\mathbf{n}\cdot\nabla\eta|_{\partial\Omega}=\mathbf{n}\cdot
M\nabla\frac{\delta{\cal E}}{\delta
c}\Big{|}_{\partial\Omega}=\mathbf{n}\cdot{\sigma^{\textnormal{el}}}|_{\partial\Omega}=0,$

where $\sigma^{\textnormal{el}}=C:\epsilon^{\textnormal{el}}$.

Several important properties of the phase field model for Ni-based alloys are
described in the following theorem.

######  Theorem 2.1.

The Ni-based alloys phase field model given by Eq. 5 and Eq. 4 together with
periodic boundary conditions or homogeneous Neumann boundary conditions Eq. 7
satisfies the following properties.

  * 1.

Mass conservation law:
$\int_{\Omega}c(\mathbf{x},t)\mathrm{d}\mathbf{x}=\int_{\Omega}c(\mathbf{x},0)\mathrm{d}\mathbf{x}$.

  * 2.

Energy dissipation law: $\frac{\mathrm{d}}{\mathrm{d}t}{\cal E}\leq 0$.

  * 3.

The mechanical equilibrium equation: $\frac{\delta{\cal
E}}{\delta\mathbf{u}}=0$.

######  Proof 2.2.

Integrating the first equation of Eq. 5 over domain $\Omega$, we obtain the
conservation law
$\int_{\Omega}c(\mathbf{x},t)\mathrm{d}\mathbf{x}=\int_{\Omega}c(\mathbf{x},0)\mathrm{d}\mathbf{x}$
immediately. According to Eq. 3 and Eq. 4, we have

(8) $\frac{\delta{\cal E}}{\delta\mathbf{u}}=-\nabla\cdot\frac{\partial
E_{\textnormal{el}}}{\partial(\nabla\mathbf{u})}=-\nabla\cdot\frac{\partial
E_{\textnormal{el}}}{\partial\epsilon}=-\nabla\cdot\sigma^{\textnormal{el}}=0.$

The free energy functional of the Ni-based alloys phase field system satisfies
the following equation

(9) $\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}{\cal E}$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}t}\int\limits_{\Omega}\left\\{E_{\textnormal{G}}(c,\eta)+\frac{\gamma_{c}}{2}|\nabla
c|^{2}+\frac{3\gamma_{\eta}}{2}|\nabla\eta|^{2}+E_{\textnormal{el}}(c)\right\\}\hbox{d}\mathbf{x}$
$\displaystyle=\int\limits_{\Omega}\Big{\\{}\frac{\partial
E_{\textnormal{G}}}{\partial c}\frac{\partial c}{\partial t}+\gamma_{c}\nabla
c\cdot\frac{\partial\nabla c}{\partial t}+\frac{\partial
E_{\textnormal{el}}}{\partial c}\frac{\partial c}{\partial t}+\frac{\partial
E_{\textnormal{el}}}{\partial(\nabla\mathbf{u})}\cdot\frac{\partial(\nabla\mathbf{u})}{\partial
t}$ $\displaystyle\qquad+\frac{\partial
E_{\textnormal{G}}}{\partial\eta}\frac{\partial\eta}{\partial
t}+3\gamma_{\eta}\nabla\eta\cdot\frac{\partial\nabla\eta}{\partial
t}\Big{\\}}\mathrm{d}\mathbf{x}$
$\displaystyle=\int\limits_{\Omega}\left\\{\frac{\delta{\cal E}}{\delta
c}\frac{\partial c}{\partial t}+\frac{\delta{\cal
E}}{\delta\eta}\frac{\partial\eta}{\partial
t}\right\\}\mathrm{d}\mathbf{x}+B_{1}$
$\displaystyle=-\int\limits_{\Omega}\left\\{M\nabla\frac{\delta{\cal
E}}{\delta c}\cdot\nabla\frac{\delta{\cal E}}{\delta
c}+\left(\frac{\delta{\cal
E}}{\delta\eta}\right)^{2}\right\\}\mathrm{d}\mathbf{x}+B_{2},$

where $B_{1}$ and $B_{2}$ are the boundary terms due to integration-by-parts,
which vanishes with the given boundary conditions. Thus we have

(10) $\frac{\mathrm{d}}{\mathrm{d}t}{\cal E}\leq 0,$

which demonstrates the energy dissipation of the Ni-based alloys phase field
system.

By denoting ${\cal E}_{\delta}:={\cal E}(c+\delta
c,\eta+\delta\eta,\mathbf{u}+\delta\mathbf{u})$, the following integral
relationship between the variational derivative and the free energy holds

(11) $\displaystyle{\cal E}_{\delta}-{\cal E}(c,\eta,\mathbf{u})$
$\displaystyle=\int\limits_{\Omega}\Big{\\{}\frac{\partial
E_{\textnormal{G}}}{\partial c}\delta c+\gamma_{c}\nabla c\cdot\nabla\delta
c+\frac{\partial E_{\textnormal{el}}}{\partial c}\delta c+\frac{\partial
E_{\textnormal{el}}}{\partial(\nabla\mathbf{u})}\cdot\nabla\delta\mathbf{u}$
$\displaystyle\qquad+\frac{\partial
E_{\textnormal{G}}}{\partial\eta}\delta\eta+3\gamma_{\eta}\nabla\eta\cdot\nabla\delta\eta\Big{\\}}\mathrm{d}\mathbf{x}$
$\displaystyle\approx\int\limits_{\Omega}\left\\{\frac{\delta{\cal E}}{\delta
c}\delta c+\frac{\delta{\cal
E}}{\delta\eta}\delta\eta\right\\}\mathrm{d}\mathbf{x}+B.$

The first equality is obtained by using Taylor expansion, and the second
equality comes from the integration-by-parts formula and $\frac{\delta{\cal
E}}{\delta\mathbf{u}}=0$. Here $B$ represents the boundary terms from the
integration-by-parts formula and equals zero due to the given boundary
conditions. Equation (11) shows the connection between the variational
derivative and the total free energy functional, which plays an important role
in the construction of the energy-stable numerical scheme for the Ni-based
alloys phase field system.

##  3 Discretization for the Ni-based alloys phase field model

We rewrite the Ni-based alloys phase field model given by Eq. 5 and Eq. 4 in a
vector form

(12) $\left\\{\begin{aligned} &\frac{\partial\mathbf{V}}{\partial t}=-{\cal
A}\mathbf{G},\\\
&\nabla\cdot\sigma=\nabla\cdot\sigma^{0},\end{aligned}\right.$

where $\mathbf{V}=(c,\,\eta)^{T}$, ${\cal A}=\textnormal{diag}(-\nabla\cdot
M\nabla,1)$ with the mobility $M=\varkappa c(1-c)$, and
$\mathbf{G}:=\frac{\delta{\cal E}}{\delta\mathbf{V}}$ represents the
variational derivative of the free energy functional ${\cal E}$ with respect
to $\mathbf{V}$. The local free energy is decomposed into four parts, which
are

(13) $\displaystyle E_{1}(\mathbf{V})=\Phi(c,\eta),$ $\displaystyle
E_{2}(\mathbf{V})=\theta\left\\{\Psi(c)+\frac{3}{4}\Psi(c\eta)+\frac{1}{4}\Psi(c\eta^{\prime})\right\\},$
$\displaystyle E_{3}(\mathbf{V},\nabla\mathbf{u})=E_{\textnormal{el}}(c),$
$\displaystyle E_{4}(\nabla\mathbf{V}^{T})=\frac{\gamma_{c}}{2}|\nabla
c|^{2}+\frac{3\gamma_{\eta}}{2}|\nabla\eta|^{2}.$

It is easy to check that ${\cal A}$ is a semi-positive operator for $c\in(0,\
0.25)$ and $\eta\in[0,\ 1]$.

Without loss of generality, we consider discretizing the Ni-based alloys phase
field system on a one dimensional domain $\Omega=[0,L_{x}]$, which is covered
by a uniform mesh with the mesh size $\Delta x=L_{x}/N_{x}$. The temporal
interval $[0,{\cal T}]$ is split by a set of nonuniform points
$\\{t_{n}\\}_{0}^{N_{\cal T}}$ with the $n$-th time step size $\Delta
t_{n}=t_{n+1}-t_{n}$. Let us denote
$\tilde{\mathbf{V}}^{n}_{i}\approx\mathbf{V}(x_{i},t_{n})$ and
$\tilde{{u}}^{n}_{i}\approx{u}(x_{i},t_{n})$ as the approximate solutions of
the Ni-based alloys phase field system at $t_{n}$, where
$x_{i}=(i-\frac{1}{2})\Delta x$ with $1\leq i\leq N_{x}$. Throughout the
paper, notations with a tilde are the corresponding approximate solutions or
functions at the discrete level. We denote
$[\tilde{\phi}]_{i}^{n+\frac{1}{2}}:=\frac{1}{2}([\tilde{\phi}]_{i}^{n}+[\tilde{\phi}]_{i}^{n+1})$
as the approximation of any function $\phi$ at $t_{n+\frac{1}{2}}$. Let us
introduce some useful notations as following:

$[D_{x}^{+}\tilde{\phi}]^{n}_{i}=\frac{\tilde{\phi}^{n}_{i+1}-\tilde{\phi}^{n}_{i}}{\Delta
x},\quad\
[D_{x}^{-}\tilde{\phi}]^{n}_{i}=\frac{\tilde{\phi}^{n}_{i}-\tilde{\phi}^{n}_{i-1}}{\Delta
x},\quad\
[D_{x}\tilde{\phi}]^{n}_{i}=\frac{\tilde{\phi}^{n}_{i+\frac{1}{2}}-\tilde{\phi}^{n}_{i-\frac{1}{2}}}{\Delta
x},$
$[\hat{D}_{x}\tilde{\phi}]^{n}_{i}=\frac{\tilde{\phi}^{n}_{i+1}-\tilde{\phi}^{n}_{i-1}}{2\Delta
x},\quad\
\left[(D_{x}^{\pm}\tilde{\phi})^{2}\right]^{n}_{i}=\frac{\left([D_{x}^{+}\tilde{\phi}]^{n}_{i}\right)^{2}+\left([D_{x}^{-}\tilde{\phi}]^{n}_{i}\right)^{2}}{2}.$

The operator $\nabla\cdot M\nabla$ in $\cal A$ is discretized by
$[(\nabla\cdot\tilde{M}\nabla)\tilde{\phi}]_{i}^{n}=[D_{x}\tilde{M}D_{x}\tilde{\phi}]_{i}^{n},$
where the approximate value of the mobility $M$ is calculated as

$\tilde{M}_{i+\frac{1}{2}}=\varkappa\frac{\tilde{c}^{n}_{i}(1-\tilde{c}^{n}_{i})+\tilde{c}^{n}_{i+1}(1-\tilde{c}^{n}_{i+1})}{2}.$

The discrete operator $\tilde{{\cal A}}$ for the Ni-based alloy phase field
system is then defined as $\tilde{{\cal
A}}(\tilde{\varphi}_{i}^{n},\tilde{\phi}_{i}^{n})^{T}=(-[\nabla\cdot\tilde{M}\nabla\tilde{\varphi}]_{i}^{n},\tilde{\phi}_{i}^{n})^{T}$.
With the periodic boundary conditions or the homogeneous Neumann boundary
conditions, we present vital formula

(14)
$\displaystyle-\sum_{i=1}^{N_{x}}\tilde{\varphi}^{n}_{i}\left[(D_{x}\tilde{M}D_{x})\tilde{\phi}\right]^{n}_{i}\color[rgb]{0,0,0}=\frac{1}{2}\sum_{i=1}^{N_{x}}\left(\tilde{M}^{n}_{i+\frac{1}{2}}[D_{x}^{+}\tilde{\varphi}]^{n}_{i}[D_{x}^{+}\tilde{\phi}]^{n}_{i}\right).$

It follows from the Cauchy–Schwarz inequality and (14) that

(15) $\left\\{\begin{aligned}
&\left<\tilde{\phi}^{n}_{i}[{A}_{11}\tilde{\phi}]_{i}^{n}\right>\geq 0,\\\
&\left<\tilde{\varphi}^{n}_{i}[{A}_{11}\tilde{\phi}]_{i}^{n}\right>=\left<\tilde{\phi}^{n}_{i}[{A}_{11}\tilde{\varphi}]_{i}^{n}\right>,\\\
&\left<\tilde{\varphi}^{n}_{i}[{A}_{11}\tilde{\phi}]_{i}^{n}\right>\leq\left<\tilde{\phi}^{n}_{i}[{A}_{11}\tilde{\phi}]_{i}^{n}\right>^{1/2}\left<\tilde{\varphi}^{n}_{i}[{A}_{11}\tilde{\varphi}]_{i}^{n}\right>^{1/2},\end{aligned}\right.$

where $A_{11}=-\nabla\cdot M\nabla$ and $\left<\Box_{i}\right>$ is defined as
$\left<\Box_{i}\right>:=\sum_{i=1}^{N_{x}}\Box_{i}\Delta x$. Thus, we conclude
that the discrete operator $\tilde{\cal A}$ for the Ni-based alloys phase
field system is semi-positive.

With the aforementioned notations, we construct a semi-implicit scheme for the
Ni-based alloys phase field system Eq. 12 by the following steps. In the first
step, the discretizations of the local free energies $E_{i}$ (i=1, 2, 3, 4)
and the total free energy ${\cal E}$ at time $t_{n}$ are respectively defined
as

(16) $\left\\{\begin{aligned}
\left[\tilde{E}_{1}\right]_{i}^{n}=&E_{1}(\tilde{c}_{i}^{n},\tilde{\eta}_{i}^{n}),\\\
\left[\tilde{E}_{2}\right]_{i}^{n}=&E_{2}(\tilde{c}_{i}^{n},\tilde{\eta}_{i}^{n}),\\\
\left[\tilde{E}_{3}\right]_{i}^{n}=&\frac{1}{2}\left[\tilde{\epsilon}^{\textnormal{el}}\right]_{i}^{n}:C:\left[\tilde{\epsilon}^{\textnormal{el}}\right]_{i}^{n},\\\
\left[\tilde{E}_{4}\right]_{i}^{n}=&\frac{\gamma_{c}}{2}\left[(D_{x}^{\pm}\tilde{c})^{2}\right]^{n}_{i}+\frac{3\gamma_{\eta}}{2}\left[(D_{x}^{\pm}\tilde{\eta})^{2}\right]^{n}_{i},\end{aligned}\right.$

and

(17) $\displaystyle\tilde{{\cal
E}}^{n}:=\sum_{i=1}^{N_{x}}\left([\tilde{E}_{1}]^{n}_{i}+[\tilde{E}_{2}]^{n}_{i}+[\tilde{E}_{3}]^{n}_{i}+[\tilde{E}_{4}]^{n}_{i}\right)\Delta
x.$

In Eq. 16, the discrete elastic strain
$\left[\tilde{\epsilon}^{\textnormal{el}}\right]_{i}^{n}$ is defined as

$\left[\tilde{\epsilon}^{\textnormal{el}}\right]_{i}^{n}=\left[{\tilde{\epsilon}}\right]^{n}_{i}-\epsilon_{0}\left(\tilde{c}^{n}_{i}-\bar{c}^{0}\right),$

where $\bar{c}^{0}=\frac{1}{|\Omega|}\int_{\Omega}c(x,0)\mathrm{d}x$ and the
heterogeneous strain $\left[\tilde{\epsilon}\right]^{n}_{i}$ is discretized as

$\left[\tilde{\epsilon}\right]^{n}_{i}=[\hat{D}_{x}\tilde{u}]^{n}_{i}.$

The discrete elastic strain energy is obtained by solving the following
discretization system for the mechanical equilibrium equation

(18) $\left[\hat{D}_{x}\tilde{\sigma}^{\textnormal{el}}\right]_{i}^{n}=0,$

where
$\left[\tilde{\sigma}^{\textnormal{el}}\right]_{i}^{n}=C\left[\tilde{\epsilon}^{\textnormal{el}}\right]_{i}^{n}$.

In the second step, a discrete form of the variational derivative
${\mathbf{G}}$ is constructed. By taking
$\mathbf{V}(x_{i},t_{n}):=\tilde{\mathbf{V}}^{n}_{i}$ and
$\delta\mathbf{V}\big{|}_{x_{i}}:=\tilde{\mathbf{V}}^{n+1}_{i}-\tilde{\mathbf{V}}^{n}_{i}$,
we choose the discrete variational derivative $\tilde{\mathbf{G}}_{i}$ such
that the following summation formula exactly holds

(19) $\displaystyle\begin{aligned} \tilde{{\cal E}}^{n+1}-\tilde{{\cal
E}}^{n}&=\left<\tilde{\mathbf{G}}_{i}\cdot\left(\tilde{\mathbf{V}}_{i}^{n+1}-\tilde{\mathbf{V}}_{i}^{n}\right)\right>.\end{aligned}$

Equation (19) can be viewed as a discrete form of (11). According to (13) and
(19), we can derive the discrete variational derivative to the following
formula

$\tilde{\mathbf{G}}_{i}:=\sum_{K=1}^{4}\tilde{\mathbf{G}}^{K}_{i}=\sum_{K=1}^{3}\mathbf{G}_{K}(\tilde{\mathbf{V}}_{i}^{n+1},\tilde{\mathbf{V}}_{i}^{n})+\mathbf{G}_{4}([\nabla\tilde{\mathbf{V}}^{T}]_{i}^{n+1},[\nabla\tilde{\mathbf{V}}^{T}]_{i}^{n}),$

where $\tilde{\mathbf{G}}^{K}_{i}$ are derived from $E_{K}$ with $K=1,\ 2,\
3,$ and $4$ in local free energy, respectively. Following the framework
proposed in [17], we have that $\tilde{\mathbf{G}}^{1}_{i}$ is also in
polynomial form. The cumbersome definition of $\tilde{\mathbf{G}}^{1}_{i}$ is
not given here but instead in Appendix B. $\tilde{\mathbf{G}}^{3}_{i}$ and
$\tilde{\mathbf{G}}^{4}_{i}$ are respectively defined as follows

(20) $\tilde{\mathbf{G}}^{3}_{i}:=\left(\begin{array}[]{l}\begin{aligned}
-\epsilon_{0}C[\tilde{\epsilon}^{\textnormal{el}}]^{n+\frac{1}{2}}_{i}\\\
0\qquad\qquad\end{aligned}\end{array}\right),$

and

(21) $\tilde{\mathbf{G}}^{4}_{i}:=\left(\begin{array}[]{l}\begin{aligned}
-{\gamma_{c}}\left[D_{x}^{2}\tilde{c}\right]_{i}^{n+\frac{1}{2}}\\\
-{3\gamma_{\eta}}\left[D_{x}^{2}\tilde{\eta}\right]_{i}^{n+\frac{1}{2}}\end{aligned}\end{array}\right).$

However, since $E_{2}$ is not a polynomial of $\mathbf{V}$, there does not
exist a $\tilde{\mathbf{G}}^{2}_{i}$ with polynomial form such that (19)
exactly holds. Consider an alternative form of $\tilde{\mathbf{G}}^{2}_{i}$ as

(22) $\tilde{\mathbf{G}}^{2}_{i}:=\frac{\theta}{4}\left(\begin{aligned}
4\Psi[\tilde{c}_{i}^{n+1},\tilde{c}_{i}^{n}]+3\eta^{n+\frac{1}{2}}_{i}\Psi[p_{i}^{n+1},p_{i}^{n}]+[\tilde{\eta}^{\prime}]^{n+\frac{1}{2}}_{i}\Psi[q_{i}^{n+1},q_{i}^{n}]\\\
3\tilde{c}^{n+\frac{1}{2}}_{i}\left(\Psi[p_{i}^{n+1},p_{i}^{n}]-\Psi[q_{i}^{n+1},q_{i}^{n}]\right)\qquad\qquad\end{aligned}\right)$

such that
$\tilde{\mathbf{G}}^{2}_{i}\cdot\left(\tilde{\mathbf{V}}_{i}^{n+1}-\tilde{\mathbf{V}}_{i}^{n}\right)=[\tilde{E}_{2}]_{i}^{n+1}-[\tilde{E}_{2}]_{i}^{n}$.
Here $p_{i}^{n}=\tilde{c}_{i}^{n}\tilde{\eta}_{i}^{n}$,
$q_{i}^{n}=\tilde{c}_{i}^{n}(4-3\tilde{\eta}_{i}^{n})$, and
$\Psi[y_{1},y_{2}]:=(\Psi(y_{1})-\Psi(y_{2}))/(y_{1}-y_{2})$ is the first
order finite quotient of function $\Psi$. Using the trapezoidal rule at the
half-time level, we obtain a fully discretized scheme for the Ni-based alloys
phase field system as:

(23) $\left\\{\begin{aligned}
&\frac{\tilde{\mathbf{V}}_{i}^{n+1}-\tilde{\mathbf{V}}_{i}^{n}}{\Delta
t_{n}}=-\tilde{{\cal A}}\tilde{\mathbf{G}}_{i},\\\
&\left[\hat{D}_{x}\tilde{\sigma}^{\textnormal{el}}\right]_{i}^{n+1}=0.\end{aligned}\right.$

The stability of the proposed scheme (23) is given by the following theorem.

######  Theorem 3.1.

For any given time step $\Delta t_{n}>0$, the numerical scheme (23) is
unconditionally energy stable and satisfies the following properties.

  * 1.

Mass conservation law:
$\left<\tilde{c}^{n+1}_{i}\right>=\left<\tilde{c}^{n}_{i}\right>$.

  * 2.

Energy dissipation law: $\tilde{\cal E}^{n+1}\leq\tilde{\cal E}^{n}$.

######  Proof 3.2.

According to (14) and (23), by setting $\varphi=1$ and $\phi=\frac{\delta{\cal
E}}{\delta c}$, we have the following mass conservation law

(24) $\left<\tilde{c}^{n+1}_{i}\right>-\left<\tilde{c}^{n}_{i}\right>=\Delta
t_{n}\left<\left[(D_{x}\tilde{M}D_{x})\tilde{\phi}\right]^{n}_{i}\right>=0.$

To prove the energy dissipation law, we shall show that the summation formula
(19) is exactly held with the definitions of discrete variational derivatives
$\tilde{\mathbf{G}}_{i}$. Similar to the proof of Eq. 14, we obtain the
following equation from the discrete mechanical equilibrium scheme Eq. 18

(25)
$\left<[\tilde{\epsilon}]^{n}_{i}[\tilde{\epsilon}^{\textnormal{el}}]^{n+1}_{i}\right>=\left<[\tilde{\epsilon}]^{n+1}_{i}[\tilde{\epsilon}^{\textnormal{el}}]^{n+1}_{i}\right>=0.$

Then, we get the following formula

(26)
$\displaystyle\left<\left[\tilde{E}_{\textnormal{el}}\right]_{i}^{n+1}-\left[\tilde{E}_{\textnormal{el}}\right]_{i}^{n}\right>$
$\displaystyle\qquad=-C\left<\left[\tilde{\epsilon}^{\textnormal{el}}\right]_{i}^{n+\frac{1}{2}}(\left[\tilde{\epsilon}^{0}\right]_{i}^{n+1}-\left[\tilde{\epsilon}^{0}\right]_{i}^{n})\right>+C\left<(\left[\tilde{\epsilon}\right]_{i}^{n+1}-\left[\tilde{\epsilon}\right]_{i}^{n})\left[\tilde{\epsilon}^{\textnormal{el}}\right]_{i}^{n+\frac{1}{2}}\right>$
$\displaystyle\qquad=-\left<\epsilon_{0}C\left[\tilde{\epsilon}^{\textnormal{el}}\right]_{i}^{n+\frac{1}{2}}(\tilde{c}_{i}^{n+1}-\tilde{c}_{i}^{n})\right>=\left<\tilde{\mathbf{G}}^{3}_{i}\cdot\left({\tilde{\mathbf{V}}_{i}^{n+1}-\tilde{\mathbf{V}}_{i}^{n}}\right)\right>.$

Using the definitions of $\tilde{\mathbf{G}}^{1}_{i}$,
$\tilde{\mathbf{G}}^{2}_{i}$, $\tilde{\mathbf{G}}^{4}_{i}$ and Eq. 14, we have

(27) $\left\\{\begin{aligned}
&\left<\left[\tilde{E}_{1}\right]_{i}^{n+1}-\left[\tilde{E}_{1}\right]_{i}^{n}\right>=\left<\tilde{\mathbf{G}}^{1}_{i}\cdot\left({\tilde{\mathbf{V}}_{i}^{n+1}-\tilde{\mathbf{V}}_{i}^{n}}\right)\right>,\\\
&\left<\left[\tilde{E}_{2}\right]_{i}^{n+1}-\left[\tilde{E}_{2}\right]_{i}^{n}\right>=\left<\tilde{\mathbf{G}}^{2}_{i}\cdot\left({\tilde{\mathbf{V}}_{i}^{n+1}-\tilde{\mathbf{V}}_{i}^{n}}\right)\right>,\\\
&\left<\left[\tilde{E}_{4}\right]_{i}^{n+1}-\left[\tilde{E}_{4}\right]_{i}^{n}\right>=\left<\tilde{\mathbf{G}}^{4}_{i}\cdot\left({\tilde{\mathbf{V}}_{i}^{n+1}-\tilde{\mathbf{V}}_{i}^{n}}\right)\right>.\\\
\end{aligned}\right.$

By combining (15), (19), (23), Eq. 26, and Eq. 27, the discrete free energy
induced by scheme (23) satisfies

(28) $\displaystyle\frac{\tilde{\cal E}^{n+1}-\tilde{\cal E}^{n}}{\Delta
t_{n}}=\left<\tilde{\mathbf{G}}_{i}\cdot\frac{\tilde{\mathbf{V}}_{i}^{n+1}-\tilde{\mathbf{V}}_{i}^{n}}{\Delta
t_{n}}\right>=-\left<\tilde{\mathbf{G}}_{i}\cdot\tilde{{\cal
A}}\tilde{\mathbf{G}}_{i}\right>\leq 0.$

Since the total free energy functional $\cal E$ is non-convex, the existence
and uniqueness of the solution for system (23) (especially for large values of
$\Delta t$) are not immediate. Nonetheless, the proof of Theorem 2.1 does not
require a unique solution to (23). Even if one could prove that scheme (23) is
unconditionally energy stable, it is numerically unstable when applied to
solve the Ni-based alloys phase field system. The reason is that the numerical
computation of the first order finite quotient $\Psi[y_{1},y_{2}]$ is unstable
and inaccurate when $y_{2}-y_{1}$ is close to zero. Unfortunately, in the Ni-
based alloys phase field system, $\delta p_{i}=p_{i}^{n+1}-p_{i}^{n}$ or
$\delta q_{i}=q_{i}^{n+1}-q_{i}^{n}$ is often close to zero, which indicates
that the numerical calculation of $\tilde{\mathbf{G}}^{2}_{i}$ by (22) is
numerically unstable and inaccurate in this situation. As a result, the scheme
(23) may lead to an inaccurate or non-physical solution with the existence of
the complicated function in $E_{2}$.

To overcome the difficulty, we propose a stable and highly accurate
approximation for $\tilde{\mathbf{G}}^{2}_{i}$ in the third step, which was
first introduced by us in [17] for solving a coupled Allen–Cahn and
Cahn–Hilliard system. Let us assume
$\Psi_{S}(y,\delta)=\Psi(y)+\sum_{s=1}^{2S}\frac{\delta^{s}\textnormal{d}^{s}\Psi(y)}{s!(\textnormal{d}y)^{s}}$,
which is a polynomial of the small parameter $\delta$. The first order
quotient $\Psi[y_{1},y_{2}]$ can be approximately computed by

$\Psi[y_{1},y_{2}]\approx\Psi_{S}[y_{1},y_{2}]=\frac{\Psi_{S}((y_{1}+y_{2})/2,\delta)-\Psi_{S}((y_{1}+y_{2})/2,-\delta)}{2\delta},$

which is a polynomial of the small parameter $\delta$ with
$\delta=(y_{2}-y_{1})/2$. Then, $\tilde{\mathbf{G}}^{2}_{i}$ is approximated
as

(29) $\displaystyle\tilde{\mathbf{G}}^{2,S}_{i}:$
$\displaystyle=\frac{\theta}{4}\left(\begin{aligned}
4\Psi_{S}[\tilde{c}_{i}^{n+1},\tilde{c}_{i}^{n}]+3\eta^{n+\frac{1}{2}}_{i}\Psi_{S}[p_{i}^{n+1},p_{i}^{n}]+[\tilde{\eta}^{\prime}]^{n+\frac{1}{2}}_{i}\Psi_{S}[q_{i}^{n+1},q_{i}^{n}]\\\
3\tilde{c}^{n+\frac{1}{2}}_{i}\left(\Psi_{S}[p_{i}^{n+1},p_{i}^{n}]-\Psi_{S}[q_{i}^{n+1},q_{i}^{n}]\right)\qquad\qquad\end{aligned}\right)$
$\displaystyle=\tilde{\mathbf{G}}^{2}_{i}-{\mathbf{R}}^{S}_{i}.$

Based on the Taylor expansion of $\Psi$, we have the truncation error
${{\mathbf{R}}}^{S}_{i}\rightarrow 0$ as $S\rightarrow\infty$.

With the approximation, the scheme (23) for the Ni-bases phase field system is
replaced with

(30) $\left\\{\begin{aligned}
&\frac{\tilde{\mathbf{V}}_{i}^{n+1}-\tilde{\mathbf{V}}_{i}^{n}}{\Delta
t_{n}}=-\tilde{{\cal A}}\tilde{\mathbf{G}}_{i}^{S},\\\
&\left[\hat{D}_{x}\tilde{\sigma}^{\textnormal{el}}\right]_{i}^{n+1}=0,\end{aligned}\right.$

where the approximate discrete variational derivative
$\tilde{\mathbf{G}}_{i}^{S}$ is defined as

(31)
$\tilde{\mathbf{G}}_{i}^{S}:=\tilde{\mathbf{G}}^{1}_{i}+\tilde{\mathbf{G}}_{i}^{2,S}+\tilde{\mathbf{G}}_{i}^{3}+\tilde{\mathbf{G}}_{i}^{4}.$

It should be noted that the solution obtained by using scheme (30) is usually
different from the one by scheme (23) due to the different approximations made
by the two schemes. For simplicity, we still denote the solution of scheme
(30) as $\tilde{\mathbf{V}}_{i}^{n+1}$. According to (19) and (29), we have

(32) $\displaystyle\qquad\quad\tilde{{\cal E}}^{n+1}-\tilde{{\cal
E}}^{n}=\left<\tilde{\mathbf{G}}_{i}\cdot\left(\tilde{\mathbf{V}}_{i}^{n+1}-\tilde{\mathbf{V}}_{i}^{n}\right)\right>=\left<\left(\tilde{\mathbf{G}}_{i}^{S}+{\mathbf{R}}^{S}_{i}\right)\cdot\left(\tilde{\mathbf{V}}_{i}^{n+1}-\tilde{\mathbf{V}}_{i}^{n}\right)\right>.~{}~{}~{}~{}$

The stability of the proposed scheme (30) is given by the following theorem.

######  Theorem 3.3.

Given $\Delta t_{n}>0$, there exists an integer $S_{0}$ such that the scheme
(30) is unconditionally energy stable and satisfies the following properties.

  * 1.

Mass conservation law:
$\left<\tilde{c}^{n+1}_{i}\right>=\left<\tilde{c}^{n}_{i}\right>$.

  * 2.

Energy dissipation law: $\tilde{\cal E}^{n+1}\leq\tilde{\cal E}^{n}$ for any
$S>S_{0}$.

######  Proof 3.4.

Similar to Theorem 3.1, one can complete the proof of mass conservation law.
According to (15), (30), and (32), the discrete free energy decided by scheme
(30) satisfies

(33) $\displaystyle\frac{\tilde{\cal E}^{n+1}-\tilde{\cal E}^{n}}{\Delta
t_{n}}=\left<\tilde{\mathbf{G}}_{i}\cdot\frac{\tilde{\mathbf{V}}_{i}^{n+1}-\tilde{\mathbf{V}}_{i}^{n}}{\Delta
t_{n}}\right>=-\left<\tilde{\mathbf{G}}_{i}\cdot\tilde{{\cal
A}}\tilde{\mathbf{G}}^{S}_{i}\right>$
$\displaystyle\qquad=-\left<\tilde{\mathbf{G}}_{i}\cdot\tilde{{\cal
A}}\tilde{\mathbf{G}}_{i}-\tilde{\mathbf{G}}_{i}\cdot\tilde{{\cal
A}}{\mathbf{R}}_{i}^{S}\right>$
$\displaystyle\qquad\leq-\left<\tilde{\mathbf{G}}_{i}\cdot\tilde{{\cal
A}}\tilde{\mathbf{G}}_{i}\right>+\left<\tilde{\mathbf{G}}_{i}\cdot\tilde{{\cal
A}}\tilde{\mathbf{G}}_{i}\right>^{1/2}\left<{\mathbf{R}}^{S}_{i}\cdot\tilde{{\cal
A}}{\mathbf{R}}^{S}_{i}\right>^{1/2}.$

If $\left<\tilde{\mathbf{G}}_{i}\cdot\tilde{{\cal
A}}\tilde{\mathbf{G}}_{i}\right>=0$, we have $\tilde{\cal
E}^{n+1}\leq\tilde{\cal E}^{n}$. Otherwise, we set
$\left<\tilde{\mathbf{G}}_{i}\cdot\tilde{{\cal
A}}\tilde{\mathbf{G}}_{i}\right>=\varrho>0$. Since
${\mathbf{R}}^{S}_{i}\rightarrow 0$ as $S\rightarrow\infty$, there exists an
integer $S_{0}>0$ such that $\left<{\mathbf{R}}^{S}_{i}\cdot\tilde{{\cal
A}}{\mathbf{R}}^{S}_{i}\right>\leq\varrho/2$ holds for any $S>S_{0}$. Thus, we
have $\tilde{\cal E}^{n+1}\leq\tilde{\cal E}^{n}$, which completes the proof
of energy dissipation law.

##  4 Newton–Krylov–Schwarz solver with adaptive time stepping

By discretizing the Ni-based phase field system with the proposed energy
stable scheme (30), a discrete nonlinear equation system ${\cal
F}(\mathbf{X}^{n})=0$ is constructed and required to be solved at each time
step. We omit the superscript $n$ in the remainder of the section. The unknown
$\mathbf{X}$ is organized by points, i.e.
$\mathbf{X}=(\tilde{c}_{1},\,\tilde{\eta}_{1}$, $\tilde{u}_{1},$
$\tilde{c}_{2},\,\tilde{\eta}_{2}$, $\tilde{u}_{2},\,\cdots)^{T}$, to maintain
the strong coupling of the concentration field and the order parameter. We
solve the nonlinear system on a parallel supercomputer with the Netwon-Krylov-
Schwarz (NKS) type algorithm [7]. Assume $\mathbf{X}_{m}$ is the current
approximate solution, then a new solution $\mathbf{X}_{m+1}$ can be obtained
by the following steps.

  * •

Compute an inexact Newton direction $\mathbf{S}_{m}$ by solving the following
Jacobian system:

(34) $\begin{array}[]{ll}J_{m}\mathbf{S}_{m}=-{\cal
F}(\mathbf{X}_{m}),\end{array}$

until $\|J_{m}\mathbf{S}_{m}+{\cal
F}(\mathbf{X}_{m})\|\leq\textnormal{max}\\{\xi_{r}\|{\cal
F}(\mathbf{X}_{m})\|,\xi_{a}\\}$. Here $J_{m}:=\nabla{\cal F}(\mathbf{X}_{m})$
is the Jacobian matrix at the solution $\mathbf{X}_{m}$, and
$\xi_{r},\xi_{\alpha}\geq 0$ are the relative and absolute tolerances for the
linear iteration, respectively.

  * •

Calculate a new solution through
$\mathbf{X}_{m+1}=\mathbf{X}_{m}+\lambda_{m}\mathbf{S}_{m}$, where the step
length $\lambda_{m}\in(0,1]$ is determined by a line search procedure [10].

The stopping condition for the nonlinear iteration is set as

(35) $\|{\cal F}(\mathbf{X}_{m+1})\|\leq\max\\{\varepsilon_{r}\|{\cal
F}(\mathbf{X}_{0})\|,\varepsilon_{a}\\},$

where $\varepsilon_{r},\varepsilon_{a}\geq 0$ are the relative and absolute
tolerances for the nonlinear iteration, respectively.

In the NKS algorithm, the Jacobian system (34) is solved by using an iterative
method based on Krylov subspace together with an additive Schwarz type
preconditioner. In particular, we employ the restarted Generalized Minimal
Residual (GMRES) method [30] as the iterative solver and consider three
preconditioners including the classical additive Schwarz (AS) preconditioner
[11], the left restricted additive Schwarz (left-RAS, [8]) preconditioner, and
the right restricted additive Schwarz (right-RAS, [6]) preconditioner, for
comparison purpose. We refer the readers to [6, 11, 8] for more details about
the three preconditioners.

The Ni-based alloys phase field system usually contains multiple time scales
that may vary in orders of magnitude during the nucleation, the growth, and
the coarsening of $\gamma^{\prime}$ particles. Therefore an adaptive control
of the time step size is necessary, in which the time step size is selected
based on the desired solution accuracy and the dynamic features of the system.
According to Theorem 3.3, the implicit scheme (30) is unconditional energy
stable and the time step size can be chosen from a wide range. Following the
idea of [17, 29, 44, 46], the time step size at the $n$-th time step is
predicted to

(36) $\tilde{\Delta}t_{n}=\max\left(\Delta t_{\min},\frac{\Delta
t_{\max}}{\sqrt{1+\zeta\mathbf{X}_{d}^{\prime}(t_{n})^{2}}}\right),$

where
$\mathbf{X}_{d}^{\prime}(t_{n})=\frac{\|\mathbf{X}^{n}-\mathbf{X}^{n-1}\|}{\Delta
t_{n-1}}$ corresponds to the change rate of numerical solutions on the two
previous time steps, and $\zeta$ is a positive pre-chosen parameter. In (36),
$\Delta t_{\min}$ and $\Delta t_{\max}$ are defined as the lower and upper
bounds of the time step size. Then we use the NKS algorithm to solve the
discrete system (30) with the predicted time step size $\tilde{\Delta}t_{n}$.
If the NKS solver diverges with the currently predicted time step size
$\tilde{\Delta}t_{n}$, a smaller predicted time step size
$\tilde{\Delta}t_{n}:=\tilde{\Delta}t_{n}/\sqrt{2}$ is chosen to restart the
NKS algorithm and $\zeta:=2\zeta$. The loop is broken down and the time step
size $\Delta t_{n}$ is set to be $\tilde{\Delta}t_{n}$ until the NKS solver
converges.

##  5 Numerical experiments

In this section, we investigate the numerical behavior and parallel
performance of the proposed algorithm for the Ni-based alloys phase field
system Eq. 12. We carry out several three dimensional tests to validate the
discretization of the proposed algorithm and study the morphology evolution of
$\gamma^{\prime}$ precipitates in Ni-based alloys during aging process. We
mainly focus on (1) the morphology evolution of $\gamma^{\prime}$ particles
during the nucleation, the growth, and the coarsening; (2) the performance of
the NKS algorithm with different preconditioners; and (3) the parallel
scalability of the proposed algorithm.

The algorithm for the Ni-based alloys phase field system is implemented on top
of the Portable, Extensible Toolkits for Scientific computations (PETSc)
library [3]. The numerical experiments are carried out on the BSCC-A6
supercomputer. The computing nodes of BSCC-A6 are comprised of two 32-core AMD
CPUs with 256 GB local memory, and are interconnected via a proprietary high
performance network. In the numerical experiments we use all 64 CPU cores in
each node and assign one subdomain to each core. In the discrete scheme Eq.
30, we set $S=10$ so that the error coming from the Taylor approximation can
be ignored. The stopping conditions for the nonlinear and linear iterations
are set as $\varepsilon_{r}=10^{2}\xi_{r}=1\times 10^{-8}$ and
$\varepsilon_{a}=10^{-3}\xi_{a}=1\times 10^{-6}$, respectively.

###  5.1 The shape evolution of a single $\gamma^{\prime}$ particle

The shape of the $\gamma^{\prime}$ precipitate mainly depends on the balance
of the interfacial energy and the elastic energy. To understand the elastic
energy contribution to the shape evolution of the $\gamma^{\prime}$
precipitate, we run a set of simulations by multiplying Hooke’s tensor by a
dimensionless parameter ${\cal L}$, which changes the ratio of the elastic
energy and the interfacial energy. According to [28], doubling the elastic
energy $E_{\textnormal{el}}$ achieves a similar effect of the double average
characteristic length of the particles. In all simulations, the computational
domain $\Omega$ is covered by an $80\times 80\times 80$ uniform mesh with
$h=\Delta x=\Delta y=\Delta z=0.25$. The time step size is initially set as
$\Delta t_{1}=0.01$ and then adaptively controlled with $\zeta=100$, $\Delta
t_{\textnormal{min}}=0.01$, and $\Delta t_{\textnormal{max}}=2$. A spherical
$\gamma^{\prime}$ particle is embedded in the $\gamma$ matrix and located at
the center of computational domain $\Omega$ with radius $r=7.5$. The initial
values are set as $c=0.238,\ \eta=0.01$ for $\gamma^{\prime}$ phase and
$c=0.1375,\ \eta=0.99$ for $\gamma$ phase, respectively.

Four simulations with ${\cal L}=0,\,1,\,3,$ and 5 are considered. The
simulation with ${\cal L}=0$ is corresponding to a standard phase field model
without the elastic energy, which is taken as a comparison. The simulations
with ${\cal L}=1,\,3,\,5$ are equivalent to the simulations with initial
radius $r=7.5,\,22.5,$ and $37.5$, respectively. The shape evolutions of a
single $\gamma^{\prime}$ particle for the four simulations are given in Fig.
1. As the particle size increases, the equilibrium shape of the
$\gamma^{\prime}$ particle transitions from spherical to approximately cubic.
We plot the evolutions of the total free-energy, the interfacial energy, and
the elastic energy for ${\cal L}=1,\,3,\,5$ in Fig. 2, which clearly indicates
that the proposed scheme satisfies the energy dissipation law even with a
relatively large time step size. By combining Fig. 1 and Fig. 2, we found that
the interfacial energy will dominate in the decay of the total energy, and the
equilibrium particle shape will be near-spherical as the particle size is
small. On the contrary, the elastic energy gradually dominates in the decay of
the total energy, and the equilibrium particle shape will approach a cube as
the particle size becomes larger. The history of the time step size is also
displayed in 2, which shows the time step size is successfully adjusted from
$\Delta t_{\min}$ to $\Delta t_{\max}$ by two orders of magnitude.

(a) $t=10$ (b) $t=50$

(c) $t=100$ (d) $t=300$

Figure 1: The shape evolutions of a single $\gamma^{\prime}$ particle with
${\cal L}=0,\,1,\,3,$ and 5. The contour plots for $c=0.22$ at $x=10$ are
displayed.

(a) ${\cal L}=1$ (b) ${\cal L}=3$ (c) ${\cal L}=5$

Figure 2: The evolutions of the relative total free energy $\tilde{\cal
E}^{n}-\tilde{\cal E}^{1}$, the relative interfacial energy, the relative
elastic energy (left axis), and the history of the time step size (right
axis).

###  5.2 Coarsening rate constant of a single $\gamma^{\prime}$ particle

The precipitation process of the $\gamma^{\prime}$ phase can be roughly
divided into three stages: the nucleation, the growth, and the coarsening. In
this subsection, we focus on the coarsening rate constants of a single
$\gamma^{\prime}$ particle. According to Lifshitz–Slyozov–Wagner theory [23],
the increase in the mean particle size as a function of time will follow a
cubic growth law

(37) $\left<R\right>^{3}-\left<R_{0}\right>^{3}=K(t-t_{0}),$

where $\left<R\right>$ is the mean particle size, $\left<R_{0}\right>$ and
$t_{0}$ are the mean particle size and time at the beginning of the
coarsening. $K$ is known as the coarsening rate constant, which depends on the
surface free energy of the particle-matrix interface and the diffusion
coefficient of the solute in the matrix. To verify the growth law (37), we run
simulations with different volume fractions of $\gamma^{\prime}$ phase, which
are $V_{f}=10\%,$ $17.5\%,$ $25\%,$ $32.5\%,$ and $40\%,$ respectively. The
computational domain is covered by a $60\times 60\times 60$ uniform mesh size.
The mesh size and parameters used in the adaptive time stepping strategy are
the same as in the previous simulations, except that the initial time step
size is set to $\Delta t_{\textnormal{min}}=0.001$. Initially, a spherical
$\gamma^{\prime}$ particle with $c=0.234$ is embedded in a matrix with
$c=0.147+0.087V_{f}$. The particle is located at the center of $\Omega$ with
radius $r=1.5$.

The relationship between $R^{3}$ and $t$ is illustrated in Fig. 3-(b). Through
data fitting, the dimensionless coarsening rate constants of a single
$\gamma^{\prime}$ particle are 0.24, 0.83, 1.75, 3.45, and 5.25, which are
corresponding to $V_{f}=10\%,$ $17.5\%,$ $25\%,$ $32.5\%,$ and $40\%,$
respectively. It is obvious that the dimensionless coarsening rate constant
increases as the volume fraction of $\gamma^{\prime}$ increasing. The
evolutions of the total free energy, the interfacial energy, the elastic
energy, and the average composition for the case with $V_{f}=0.4$ are shown in
Fig. 3-(a). The results in other cases are very similar and will not be
repeated here. From Fig. 3-(a), we found that the interfacial energy and the
elastic energy increase as the particle gets larger, but the total free keeps
decaying due to the energy dissipation law and unconditional energy stability
of the newly presented scheme. We also plot the average composition $\bar{c}$
in Fig. 3-(a), which verifies that the newly proposed scheme fully complies
with the mass conservation law.

(a) $V_{f}=0.4$ (b)

Figure 3: (a) The evolutions of the total free energy, the interfacial energy,
the elastic energy (left axis), and the average composition $\bar{c}$ (right
axis); (b) Evolutions of $R^{3}$ for different volume fractions.

###  5.3 Morphological evolution of the $\gamma^{\prime}$ phase and particles
number density change with increased volume fraction

In this subsection, we investigate the morphological evolution and particle
number density changes of the $\gamma^{\prime}$ phase by simulating several
tests for different volume fractions. The computational domain
$\Omega=[0,\,75]\times[0,\,75]\times[0,\times 10]$ is covered by a $300\times
300\times 40$ uniform mesh with $h=0.25$. The time step size is initially set
as $\Delta t_{1}=0.01$ and then adaptively controlled with $\zeta=100$,
$\Delta t_{\textnormal{min}}=0.01$, and $\Delta t_{\textnormal{max}}=2$. A
homogeneous solution with small composition fluctuations around the average
composition is taken as the initial value, which is set to be a randomly
distributed state $(c^{0},\,\eta^{0})=(\bar{c}+\delta_{c},0.1+\delta_{\eta})$,
where $\delta_{c}$ or $\delta_{\eta}$ is a uniform random distribution
function of $-0.05$ to 0.05. The value of the average composition $\bar{c}$ is
determined by the volume fraction of the $\gamma^{\prime}$ phase. For the
reason of comparison, we run four simulations with $\bar{c}\approx 0.1622$,
0.1687, 0.1753, and 0.1818, which are corresponding to the volume fraction
$V_{f}=17.5\%$, $25\%$, $32.5\%$, and $40\%$, respectively. All simulations
are stopped at $t=200$.

The distributions of the $\gamma^{\prime}$ phase at $t=200$ for the four
simulations are given in Fig. 4, from which we found a remarkable morphology
change of $\gamma^{\prime}$ phase with increased volume fraction. The
particles size, shape, and number density change as the volume fraction of the
$\gamma^{\prime}$ phase increases. Next, we take the simulation with
$V_{f}=17.5\%$ as an example to study the morphological evolution and particle
number density of the $\gamma^{\prime}$ phase. The pseudo-color plots of the
composition $c$ at $t=12,\,39$, and $171$ are displayed in Fig. 5 (a)-(c),
which are corresponding to three stages: the nucleation, the growth, and the
coarsening, respectively. The nucleation of new $\gamma^{\prime}$ phase forms
randomly due to small composition fluctuations. The growth of these nuclei is
caused by diffusive transportation, and the coarsening involves the growth of
large particles at the expense of the dissolution of small ones, driven by an
overall reduction in the interfacial energy. Contour plots of the third
component of displacement, $u_{3}$, on surface $z=0$ are also shown in Fig. 5
(d)-(f), from which we conclude that the elastic interaction is long range and
can lead to the spatial correlation between the $\gamma^{\prime}$ particles.

(a) $V_{f}=17.5\%$

(b) $V_{f}=25\%$

(c) $V_{f}=32.5\%$

(d) $V_{f}=40\%$

(e) $V_{f}=17.5\%$

(f) $V_{f}=25\%$

(g) $V_{f}=32.5\%$

(h) $V_{f}=40\%$

Figure 4: (a)-(d): Pseudo-color plots of composition $c$ at $t=200$ with
different volume fractions of the $\gamma^{\prime}$ phase. (e)-(h): Contour
plots of long range order parameter $\eta$ on surface $z=0$ at $t=200$.

(a) $t=12$, nucleation

(b) $t=39$, growth

(c) $t=171$, coarsening

(d) $t=12$

(e) $t=39$

(f) $t=171$

Figure 5: Volume fraction $V_{f}=17.5\%$. (a)-(c): Pseudo-color plots of
composition $c$ at different times. (d)-(f): Contour plots of the third
component of displacement $u_{3}$ on surface $z=0$ at different times.

The evolutions of the total free energy, the interfacial energy, and the
elastic energy for the simulation with $V_{f}=17.5\%$ are shown in Fig. 6-(a).
During the nucleation stage, the total free energy decreases rapidly as the
interfacial energy and the elastic energy increase. At the growth stage, the
total energy decay rate becomes smaller as the interfacial energy and the
elastic energy increase rapidly. During the coarsening stage, the total energy
decays continuously with the shape evolution of $\gamma^{\prime}$ particles,
which is mainly determined by the balance between the interfacial energy and
the elastic energy. We also plot the average composition $\bar{c}$ in Fig.
6-(a), which clearly verifies that the newly proposed scheme fully complies
with the mass conservation law. The evolution of $\left<R^{3}\right>$ is given
in Fig. 6-(b), which verifies the cubic growth law (37). During the growth
stage and early coarsening stage, the slope $K$ in the cubic growth law (37)
is near to 0.83, which is the coarsening rate constant of a single
$\gamma^{\prime}$ particle for volume fraction $V_{f}=17.5\%$. As coarsening
proceeds, the coarsening rate constant decreases and eventually approaches
zero when a steady state is reached. The number of $\gamma^{\prime}$ particles
as a function of time $t$ is also plotted in Fig. 6-(b), from which we found
that the particle number decreases rapidly during the growth stage. The
temporal evolution of the normalized particle size distribution (PSD) of
$R/\left<R\right>$ for the $\gamma^{\prime}$ phase are shown in Fig. 7, which
shows a good agreement with the Lifshitz–Slyozov–Wagner (LSW) theory [23, 38]
and the experiment results reported in [2].

(a)  (b)

Figure 6: Volume fraction $V_{f}=17.5\%$. (a) The evolutions of the total free
energy, the interfacial energy, the elastic energy (left axis), and the
average composition $\bar{c}$ (right axis). (b) The evolutions of
$\left<R^{3}\right>$ and number of $\gamma^{\prime}$ particles. For
visualization convenience, the total free energy is shifted by subtracting
$(\tilde{\cal E}^{1}-4.6\times 10^{4})$.

(a) $t=66$ (b) $t=127$ (c) $t=171$ (d) $t=200$

Figure 7: Volume fraction $V_{f}=17.5\%$. The temporal evolutions of the
normalized PSD of $R/\left<R\right>$ for the $\gamma^{\prime}$ phase.

###  5.4 Parallel scalability tests

In this subsection, we study the parallel scalability, in both the weak and
strong sense, of the proposed method. To begin, we discuss several issues
including the influence of different overlaps, subdomain solvers and so on. In
all simulations, we run the code for the first $10$ time steps with a fixed
time step size $\Delta t=0.01$ and an $80\times 80\times 80$ mesh by using
$256$ processor cores. To check the influence of the subdomain solvers, we
limit the test to the classical AS preconditioner and fix the overlapping size
to $2\delta=2$. The ILU factorizations with $0$, $1$, and $2$ levels of fill-
in and LU factorization are considered. The number of Newton iterations and
the averaged number of GMRES iterations together with the total compute times
are provided in Table 1. From Table 1, we can see the number of Newton
iterations is invariable, which means the number of Newton iterations is
insensitive to the subdomain solver. In addition, the number of GMRES
iterations decreases by increasing the fill-in level, but the total compute
time keeps growing 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) subdomain solver and the reuse strategy can further save nearly
20-30% of the compute time.

Table 1:  Performance of the NKS algorithm with different subdomain solvers.  Subdomain solver | ILU(0) | ILU(1) | ILU(2) | LU | ILU(0)-reuse
---|---|---|---|---|---
Total Newton | 25 | 24 | 25 | 25 | 25
GMRES/Newton | 23.9 | 14.5 | 14.9 | 13.5 | 24
Total Time (s) | 126.3 | 881.4 | 3,815.5 | 8,325.2 | 95.1

Based on the above tests, we take the ILU(0)-reuse as the subdomain solver and
investigate the performance of the NKS solver by changing the type of the AS
preconditioner and the overlapping factor $\delta$. The classical-AS, the
left-RAS, and the right-RAS preconditioners with overlapping size
$2\delta=0,\,2,\,4$ are considered. The number of Newton iterations, the
averaged number of GMRES iterations, and the total compute times are listed in
Table 2. From Table 2, we can see that the number of Newton iterations is
insensitive to the type of the AS preconditioner and the overlapping size. We
also conclude that the classical-AS preconditioner with $\delta=0$ is superior
to the other AS preconditioners with $\delta=1$ or $2$ based on the averaged
number of GMRES iterations and the compute times.

Table 2: Performance of NKS solver with respect to the types of preconditioner and overlapping sizes.  Type of preconditioner | classical-AS | left-RAS | right-RAS
---|---|---|---
$\delta$ | $0$ | $1$ | $2$ | $1$ | $2$ | $1$ | $2$
Total Newton | 26 | 25 | 25 | 25 | 25 | 25 | 25
GMRES/Newton | 26.2 | 24 | 19.4 | 14.9 | 15.9 | 14.8 | 17.0
Total Time (s) | 69.1 | 95.1 | 127.0 | 79.1 | 121.4 | 78.7 | 192.4

Then, we focus on parallel scalability tests to examine the weak and strong
scaling performance. Based on the above observations, we use the classical-AS
preconditioner with the overlapping $\delta=0$ and employ the ILU(0)
factorization with the reuse strategy as the subdomain solver in all
simulations. In the weak scalability tests, a subdomain with a fixed number of
mesh points is handled by each processor core. The weak scaling performance is
testified by increasing the number of the subdomains and the processor cores
simultaneously. The number of processors increases from $384$ to $3,072$, and
the number of physical grid points increases from $100^{3}$ to $200^{3}$. The
numbers of Newton and GMRES iterations together with the total compute times
are provided in Table 3, which shows that the averaged number of GMRES
iterations and the total compute time increase slowly as more processors are
used. The good weak scalability of our method is validated by the simulation.
To study the strong scalability of our method, we run the tests with a
$250^{3}$ mesh. The numbers of nonlinear and linear iterations are reported in
Table 4. We notice that the number of nonlinear iterations is almost unchanged
and the average number of linear iterations increases slightly as the number
of used processor core increases. The total compute time decreases nearly half
as the number of processor cores doubles. The overall speedup from $384$ to
$6,144$ cores is around $6.79$, which indicates a good strong parallel
efficiency of the proposed algorithm.

Table 3: The weak scalability of the proposed method. Here “NP” denotes the number of processor cores.  Mesh size | $100^{3}$ | $2\times 100^{3}$ | $4\times 100^{3}$ | $200^{3}$
---|---|---|---|---
NP | 384 | 768 | 1,536 | 3,072
Total Newton | 26 | 26 | 26 | 26
GMRES/Newton | 28.0 | 37.5 | 42.8 | 42.1
Total Time (s) | 99.8 | 115.8 | 122.4 | 124.4
Table 4:  Performance of the NKS solver with different numbers of processor cores. Here “NP” denotes the number of processor cores.  NP | 384 | 768 | 1,536 | 3,072 | 6,144
---|---|---|---|---|---
Total Newton | 25 | 25 | 25 | 26 | 26
GMRES/Newton | 40.1 | 39.6 | 40.6 | 45.7 | 48.8
Total Time (s) | 1251.2 | 634.2 | 423.3 | 236.4 | 184.4
Speed up | 1 | 1.97 | 2.96 | 5.29 | 6.79

##  6 Conclusion

The Ni-based alloys phase field system consists of CH/GL equations and linear
elastic equations, in which the total free energy includes the elastic energy
and logarithmic type functionals. By using the DVD method, we propose a semi-
implicit finite difference scheme, which solves the challenge posed by the
special free energy functional. The newly presented semi-implicit scheme is
proved to be unconditionally energy stable and enjoys the energy dissipative
law and mass conservation law. Thanks to the unconditional energy stability of
the scheme, based on the desired solution accuracy and the dynamic features of
the nucleation, growth, and coarsening procedures, the time step size is
adaptively selected by using an adaptive time stepping strategy. At each time
step, a large sparse nonlinear algebraic system is constructed by the semi-
implicit method. We introduce a parallel NKS algorithm to efficiently solve
the nonlinear system. Several three dimensional test cases are used to
validate the energy stability of the semi-implicit scheme and study the
morphological evolution of the $\gamma^{\prime}$ phase through the heat
treatment process. Large scale numerical experiments show that the proposed
algorithm can achieve high scalability with up to six thousand processor cores
in both strong and weak senses.

## Appendix A

In this appendix, we introduce more details on the Ni-based alloys phase field
system, i.e., the definition of the total free energy (1). According to [1,
24, 49], the $\gamma/\gamma^{\prime}$ two-phase microstructure in Ni-based
alloys can be described by a composition field variable $c$ of the Al element
and four long-range order parameter fields $\eta_{i}:=\eta_{i}(\mathbf{x},t)$
with $i=1,2,3$, and $4$. As shown in Fig. 8, a face center cubic (FCC)
superlattice is divided into four sub-lattices. Each sub-lattice includes four
atoms and the site fractions of Al element in sub-lattices A, B, C, and D are
denoted as $c\eta_{i}$. The overall composition $c$ can be computed by
$c=\sum\limits_{i=1}^{4}c\eta_{i}/4$, which implies constraints
$\sum\limits_{i=1}^{4}\eta_{i}=4$ and $\eta_{i}\in[0,4]$. The composition
field variable for the Ni element is set as $1-c$. For Ni-based alloys
including $\gamma$/$\gamma^{\prime}$ phases, the composition field variable
$c$ is in the range of $(0,0.25)$. In the disordered phase ($\gamma$ phase),
the four long-range order parameter fields $\eta_{i}$ are the same, i.e.,
$\eta_{i}=1$. If the first three long-range order parameter fields have the
same value and different from the fourth, there is the ordering phase
($\gamma^{\prime}$ phase). The total free energy of a $\gamma/\gamma^{\prime}$
two-phase system, including the local chemical (Gibbs) free energy, the
interfacial energy, and the elastic energy, can be written as the following
function of $c$ and $\eta_{i}$ with $i=1,2,3$ [1, 24, 49]

(38) ${\cal
E}=\int\limits_{\Omega}\left[\frac{1}{V_{m}}E_{\textnormal{G}}(c,\eta_{1},\eta_{2},\eta_{3})+\frac{\gamma_{c}}{2}|\nabla
c|^{2}+\frac{\gamma_{\eta}}{2}\sum\limits_{i=1}^{3}|\nabla\eta_{i}|^{2}+E_{\textnormal{el}}(c)\right]\hbox{d}\mathbf{x},$

where $\Omega\in\mathbb{R}^{d}$ is the computational domain, $V_{m}=1.48\times
10^{-5}{\hbox{m}}^{3}/\hbox{mol}$ is the mole volume,
$E_{\textnormal{G}}(c,\eta_{1},\eta_{2},\eta_{3})$ is the local chemical free
energy, and $E_{\textnormal{el}}(c)$ is the elastic energy. Here,
$\gamma_{c}=2.5\times 10^{-9}$J/m and $\gamma_{\eta}=6.0\times 10^{-12}$J/m
are the gradient energy coefficients of the composition and long-range order
parameters, respectively.

Figure 8: A FCC lattice is divided into four cubic sub-lattices, labeled A, B,
C, and D, each with lattice constant $a$. The D site is nonmagnetic and the
others are magnetic.

The Gibbs free energy defines the basic thermodynamic properties of the Ni-
based alloys system, which can be written in terms of the composition and the
long-range order parameter fields as [1, 49]

(39)
$E_{\textnormal{G}}(c,\eta_{1},\eta_{2},\eta_{3})=E^{\textnormal{dis}}(c)+E^{\textnormal{ord}}(c,\eta_{1},\eta_{2},\eta_{3})-E^{\textnormal{ord}}(c,1,1,1),$

where $E^{\textnormal{dis}}$ and $E^{\textnormal{ord}}$ are the Gibbs free
energies of the disordered phase and ordered phase, respectively. The Gibbs
free energy of a pure disorder phase is defined as

(40) $E^{\textnormal{dis}}(c)=E^{\textnormal{dis}}_{0}(c)+\Delta
E^{\textnormal{dis}}_{\textnormal{ideal}}(c)+\Delta
E^{\textnormal{dis}}_{\textnormal{ex}}(c),$

where the reference Gibbs energy
$E^{\textnormal{dis}}_{0}(c)=E_{0}^{\textnormal{Al}}c+E_{0}^{\textnormal{Ni}}(1-c)$,
the entropy $\Delta E^{\textnormal{dis}}_{\textnormal{ideal}}(c)=RT\Psi(c)$,
and the excess Gibbs energy $\Delta
E^{\textnormal{dis}}_{\textnormal{ex}}(c)=c(1-c)\sum\limits_{i=0}^{3}[L_{i}(2c-1)^{i}]$.
The Gibbs free energy $E^{\textnormal{ord}}$ of a pure order phase is given by

(41)
$E^{\textnormal{ord}}(c,\eta_{1},\eta_{2},\eta_{3})=E^{\textnormal{ord}}_{0}(c,\eta_{1},\eta_{2},\eta_{3})+\Delta
E^{\textnormal{ord}}_{\textnormal{ideal}}(c,\eta_{1},\eta_{2},\eta_{3})+\Delta
E^{\textnormal{ord}}_{\textnormal{ex}}(c,\eta_{1},\eta_{2},\eta_{3}),$

where the reference Gibbs energy $E^{\textnormal{ord}}_{0}$, the entropy
$\Delta E^{\textnormal{ord}}_{\textnormal{ideal}}$, and the excess Gibbs
energy $\Delta E^{\textnormal{ord}}_{\textnormal{ex}}(c)$ are respectively
defined as

(42)
$E^{\textnormal{ord}}_{0}(c,\eta_{1},\eta_{2},\eta_{3})=6U_{1}c^{2}\sum_{i=1}^{3}\phi_{i}^{2},$
(43) $\Delta
E^{\textnormal{ord}}_{\textnormal{ideal}}(c,\eta_{1},\eta_{2},\eta_{3})=\frac{RT}{4}\sum_{i=1}^{4}\Psi(c\eta_{i}),$

and

(44) $\Delta
E^{\textnormal{ord}}_{\textnormal{ex}}(c,\eta_{1},\eta_{2},\eta_{3})=-2U_{1}c^{2}\sum_{i=1}^{3}\phi_{i}^{2}+12U_{4}(1-2c)c^{2}\sum_{i=1}^{3}\phi_{i}^{2}-48U_{4}c^{3}\phi_{1}\phi_{2}\phi_{3}.$

Here $\phi_{1}=1-\frac{\eta_{1}+\eta_{2}}{2}$,
$\phi_{2}=1-\frac{\eta_{1}+\eta_{3}}{2}$, and
$\phi_{3}=1-\frac{\eta_{2}+\eta_{3}}{2}$. All the parameters in the above free
energies such as $E_{0}^{\textnormal{Al}}$, $E_{0}^{\textnormal{Ni}}$, $L_{i}$
with $i=0,1,2,3$, $U_{1}$, and $U_{4}$ are functions of temperature $T$ and
the gas constant $R$, which can be obtained from CALPHAD database [1, 31].

Since the three long-range order parameter fields $\\{\eta_{i}\\}_{i=1}^{3}$
take the same value in the $\gamma$ and $\gamma^{\prime}$ phase, the phase
field system can be approximately described by a composition field variable
$c$ and one long-range order parameter $\eta$ by taken
$\eta(\mathbf{x},t)=\eta_{1}=\eta_{2}=\eta_{3}$ [24]. The total free energy is
then rewritten as

(45) ${\cal
E}=\int\limits_{\Omega}\left[\frac{1}{V_{m}}E_{\textnormal{G}}(c,\eta)+\frac{\gamma_{c}}{2}|\nabla
c|^{2}+\frac{3\gamma_{\eta}}{2}|\nabla\eta|^{2}+E_{\textnormal{el}}(c)\right]\hbox{d}\mathbf{x}.$

The contour plot of the Gibbs free energy $E_{\textnormal{G}}(c,\eta)$ at
$T=1,037$ K is plotted in Fig. 9-(a). There are two wells on the surface and
they correspond to the equilibrium states of the disorder and order phases.
The values of the composition and the order parameter for these two phases are
$\eta=1,c=0.137$, and $\eta=0.02,c=0.229$, respectively. Fig. 9-(b) displays
the relative chemical free energy of the $\gamma$ phase and $\gamma^{\prime}$
phase of the Ni-based alloys as a function of composition $c$. In the stress-
free state, the maximum driving force is approximately 292.5 J/mol. Fig. 9-(c)
shows the phase diagram for $\gamma$/$\gamma^{\prime}$ phase with $T\in[800,\
1,300]$, which is consistent with the experiment results [26] and calculated
results by using EAM potential [27].

(a)  (b)  (c)

Figure 9: (a) Contour plot of the Gibbs free energy $E(c,\eta)$ at $T=1,073$
K. (b) The Gibbs free energy of the $\gamma$ and $\gamma^{\prime}$ phases. (c)
Phase boundary of the $\gamma$ and $\gamma^{\prime}$ phases calculated by
$E(c,\eta)$ for temperatures ranging from 800K to 1300K.

The elastic energy has important contributions to the morphological evolution
and coarsening kinetics of $\gamma^{\prime}$ particles, which can be obtained
by using the micro-elasticity theory [20, 37]. The main factors of the elastic
energy fluctuation in the computational domain are: (1) the lattice mismatch
between the different phases; (2) the elastic constant inhomogeneity; (3)
externally applied stress. In this paper, the phase field system is
constructed without external stress. According to [50], the elastic constant
difference between the $\gamma$ and $\gamma^{\prime}$ phases (usually less
than 7%) can be ignored, if no external stress is applied. Therefore, the
elastic strain energy [16, 22] is approximately calculated by:

(46)
$E_{\textnormal{el}}(c)=\frac{1}{2}\epsilon^{\textnormal{el}}:C:\epsilon^{\textnormal{el}},$

where $C$ is the Hooke’s tensor,
$\epsilon^{\textnormal{el}}=(\epsilon^{\textnormal{el}}_{IJ})$ is the elastic
strain. The elastic strain is reformulated as
$\epsilon^{\textnormal{el}}=\epsilon-\epsilon^{0}$, where $\epsilon$ is the
total strain related to the displacement and the eigenstrain $\epsilon^{0}$
expressed as $\epsilon^{0}=\epsilon_{0}(c-\bar{c}){\cal I}$. Here
$\epsilon_{0}=\frac{1}{a}\frac{\hbox{d}a}{\hbox{d}c}\approx\frac{a_{\gamma^{\prime}}-a_{\gamma}}{a_{\gamma}(c_{\gamma^{\prime}}-c_{\gamma})}=0.049$
is the composition expansion coefficient of the lattice parameter, ${\cal I}$
is the identity matrix, and
$\bar{c}=\frac{1}{|\Omega|}\int_{\Omega}c(\mathbf{x},t)\mathrm{d}\mathbf{x}$
is the average composition. The lattice constants for the $\gamma$ and
$\gamma^{\prime}$ phases are denoted as $a_{\gamma}$ and
$a_{\gamma^{\prime}}$, respectively. The Ni-based alloys can be taken as a
hexagonal crystal, which has three independent elastic constants
$C_{11},\,C_{12},$ and $C_{44}$. The elastic strain energy is finally written
as

(47) $\displaystyle
E_{\textnormal{el}}(c):=\frac{1}{2}\left(C_{12}\left(\sum\limits_{I=1}^{d}\epsilon_{II}^{\textnormal{el}}\right)^{2}+(C_{11}-C_{12})\sum\limits_{I=1}^{d}\left(\epsilon_{II}^{\textnormal{el}}\right)^{2}+4C_{44}\sum\limits_{I<J}(\epsilon_{IJ}^{\textnormal{el}})^{2}\right).$

The morphological evolution of the composition field variable
$c(\mathbf{x},t)$ and long-range order parameter $\eta(\mathbf{x},t)$ are
described by the CH–GL coupled equations:

(48) $\left\\{\begin{aligned} &\frac{\partial c}{\partial
t}=\nabla\left[M_{c}\nabla\frac{\delta{\cal E}}{\delta
c}\right]\quad\hbox{in}~{}\Omega\times(0,\mathcal{T}],\\\
&\frac{\partial\eta}{\partial t}=-M_{\eta}\frac{\delta{\cal
E}}{\delta\eta}\quad\hbox{in}~{}\Omega\times(0,\mathcal{T}],\end{aligned}\right.$

where $M_{c}$ and $M_{\eta}$ are two mobility parameters. By taken
dimensionless units for time $t$ as $M_{\eta}\Delta E$ with $\Delta
E=3.3\times 10^{7}$ $\hbox{J/m}^{3}$, for $\mathbf{x}$ as $1/{\mathbb{L}}$
with ${\mathbb{L}}=1.5$ nm, we obtain the dimensionless equations Eq. 5, in
which $\varkappa=0.008$.

## Appendix B

In this appendix, we present the definition of $\tilde{\mathbf{G}}^{1}_{i}$.
In the definition of the variational derivatives Eq. 6, $\frac{\partial
E_{\textnormal{G}}(c,\eta)}{\partial c}$, $\frac{\partial
E_{\textnormal{G}}(c,\eta)}{\partial\eta}$, and $\frac{\partial
E_{\textnormal{el}}(c)}{\partial c}$ are given as follows

(49) $\left\\{\begin{aligned} &\frac{\partial
E_{\textnormal{G}}(c,\eta)}{\partial
c}=M_{0}+M_{1}(\phi)c+M_{2}(\phi)c^{2}+M_{3}c^{3}+M_{4}c^{4}+\frac{RT}{4}\varphi^{\prime}_{c}(c),\\\
&\frac{\partial
E_{\textnormal{G}}(c,\eta)}{\partial\eta}=M_{5}(c)\phi+M_{6}(c)\phi^{2}+\frac{RT}{4}\varphi^{\prime}_{\eta},\\\
&\frac{\partial E_{\textnormal{el}}(c)}{\partial c}=-\epsilon_{0}{\cal
I}:C:\epsilon^{\textnormal{el}},\end{aligned}\right.$

where the dimensionless parameters (functions) are defined by

(50) $\left\\{\begin{aligned} &M_{0}=\frac{1}{V_{m}\Delta
E}\left[E_{0}^{\textnormal{Al}}-E_{0}^{\textnormal{Ni}}+(L_{0}-L1+L2-L3)\right],\\\
&M_{1}(\phi)=\frac{1}{V_{m}\Delta
E}\left[-(2L_{0}-6L_{1}+10L_{2}-14L_{3})+24U_{1}\phi^{2}+72U_{4}\phi^{2}\right],\\\
&M_{2}(\phi)=\frac{1}{V_{m}\Delta
E}\left[-(6L_{1}-24L_{2}+54L_{3})-216U_{4}\phi^{2}-144U_{4}\phi^{3}\right],\\\
&M_{3}=-\frac{16L_{2}-80L_{3}}{V_{m}\Delta E},\\\
&M_{4}=-\frac{40L_{3}}{V_{m}\Delta E},\\\
&M_{5}(c)=-\frac{24U_{1}c^{2}+72U_{4}(1-2c)c^{2}}{V_{m}\Delta E},\\\
&M_{6}(c)=\frac{144U_{4}c^{3}}{V_{m}\Delta E},\\\
&\varphi^{\prime}_{c}=\frac{1}{V_{m}\Delta
E}\left[3\eta\Phi^{\prime}(c\eta)+(4-3\eta)\Phi^{\prime}(c(4-3\eta))\right],\\\
&\varphi^{\prime}_{\eta}=\frac{1}{V_{m}\Delta
E}\left[3c\Phi^{\prime}(c\eta)-3c\Phi^{\prime}(c(4-3\eta))\right].\end{aligned}\right.$

The definition of $\tilde{\mathbf{G}}^{1}_{i}$ is then given by

(51) $\tilde{\mathbf{G}}^{1}_{i}=\left(\begin{aligned}
&\frac{M_{1}(\tilde{\phi}_{i}^{n})+M_{1}(\tilde{\phi}_{i}^{n+1})}{2}\tilde{c}_{i}^{n+1/2}+\frac{M_{2}(\tilde{\phi}_{i}^{n})+M_{2}(\tilde{\phi}_{i}^{n+1})}{6}\frac{(\tilde{c}_{i}^{n+1})^{3}-(\tilde{c}_{i}^{n})^{3}}{\tilde{c}_{i}^{n+1}-\tilde{c}_{i}^{n}}\\\
&\qquad\qquad+\frac{M_{3}}{4}\frac{(\tilde{c}_{i}^{n+1})^{4}-(\tilde{c}_{i}^{n})^{4}}{\tilde{c}_{i}^{n+1}-\tilde{c}_{i}^{n}}+\frac{M_{4}}{5}\frac{(\tilde{c}_{i}^{n+1})^{5}-(\tilde{c}_{i}^{n})^{5}}{\tilde{c}_{i}^{n+1}-\tilde{c}_{i}^{n}}+M_{0},\\\
&\frac{M_{5}(\tilde{c}_{i}^{n})+M_{5}(\tilde{c}_{i}^{n+1})}{2}\tilde{\phi}_{i}^{n+1/2}+\frac{M_{6}(\tilde{c}_{i}^{n})+M_{6}(\tilde{c}_{i}^{n+1})}{6}\frac{(\tilde{\phi}_{i}^{n+1})^{3}-(\tilde{\phi}_{i}^{n})^{3}}{\tilde{\phi}_{i}^{n+1}-\tilde{\phi}_{i}^{n}}\end{aligned}\right).$

## References

  * [1] I. Ansara, N. Dupin, H. L. Lukas, and B. Sundman, Thermodynamic assessment of the al-ni system, Journal of Alloys and Compounds, 247 (1997), pp. 20–30.
  * [2] A. Ardell and R. Nicholson, The coarsening of $\gamma$’in ni-al alloys, Journal of Physics and Chemistry of Solids, 27 (1966), pp. 1793–1794.
  * [3] S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. Gropp, et al., PETSc users manual, (2019).
  * [4] A. Baskaran, J. S. Lowengrub, C. Wang, and S. M. Wise, Convergence analysis of a second order convex splitting scheme for the modified phase field crystal equation, SIAM Journal on Numerical Analysis, 51 (2013), pp. 2851–2873.
  * [5] N. Bozzolo, N. Souaï, and R. E. Logé, Evolution of microstructure and twin density during thermomechanical processing in a $\gamma$-$\gamma$’nickel-based superalloy, Acta materialia, 60 (2012), pp. 5056–5066.
  * [6] X.-C. Cai, M. Dryja, and M. Sarkis, Restricted additive Schwarz preconditioners with harmonic overlap for symmetric positive definite linear systems, 41 (2003), pp. 1209–1231.
  * [7] X.-C. Cai, W. D. Gropp, D. E. Keyes, and M. D. Tidriri, Newton-Krylov-Schwarz methods in CFD, in Proceedings of the International Workshop on the Navier-Stokes Equations, Notes in Numerical Fluid Mechanics, R. Rannacher, ed., Braunschweig, 1994, Vieweg Verlag, pp. 123–135.
  * [8] X.-C. Cai and M. Sarkis, A restricted additive Schwarz preconditioner for general sparse linear systems, 21 (1999), pp. 792–797.
  * [9] L. Q. Chen and J. Shen, Applications of semi-implicit fourier-spectral method to phase field equations, Computer Physics Communications, 108 (1998), pp. 147–158.
  * [10] J. E. Dennis and R. B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Society for Industrial and Applied Mathematics, Philadelphia, 1996.
  * [11] M. Dryja and O. B. Widlund, Domain decomposition algorithms with small overlap, 15 (1994), pp. 604–620.
  * [12] Q. Du and R. A. Nicolaides, Numerical analysis of a continuum model of phase transition, SIAM Journal on Numerical Analysis, 28 (1991), pp. 1310–1322.
  * [13] C. M. Elliott and A. Stuart, The global dynamics of discrete semilinear parabolic equations, SIAM journal on numerical analysis, 30 (1993), pp. 1622–1663.
  * [14] D. J. Eyre, Unconditionally gradient stable time marching the cahn-hilliard equation, MRS Online Proceedings Library (OPL), 529 (1998).
  * [15] D. Furihata and T. Matsuo, Discrete variational derivative method : a structure-preserving numerical method for partial differential equations, Chapman and Hall/CRC, 2011.
  * [16] S. Hu and L. Chen, A phase-field model for evolving microstructures with strong elastic inhomogeneity, Acta materialia, 49 (2001), pp. 1879–1890.
  * [17] J. Huang, C. Yang, and Y. Wei, Parallel energy-stable solver for a coupled allen–cahn and cahn–hilliard system, SIAM Journal on Scientific Computing, 42 (2020), pp. C294–C312.
  * [18] L. Ju, J. Zhang, L. Zhu, and Q. Du, Fast explicit integration factor methods for semilinear parabolic equations, Journal of Scientific Computing, 62 (2015), pp. 431–455.
  * [19] A. Khachaturyan, S. Semenovskaya, and T. Tsakalakos, Elastic strain energy of inhomogeneous solids, Physical Review B, 52 (1995), p. 15909.
  * [20] A. G. Khachaturyan, Theory of structural transformation in solids, (1983).
  * [21] A. G. Khachaturyan, Theory of structural transformations in solids, Courier Corporation, 2013.
  * [22] Y.-S. Li, H. Zhu, L. Zhang, and X.-L. Cheng, Phase decomposition and morphology characteristic in thermal aging fe–cr alloys under applied strain: A phase-field simulation, Journal of nuclear materials, 429 (2012), pp. 13–18.
  * [23] I. M. Lifshitz and V. V. Slyozov, The kinetics of precipitation from supersaturated solid solutions, Journal of physics and chemistry of solids, 19 (1961), pp. 35–50.
  * [24] C. Liu, Y. Li, L. Zhu, S. Shi, S. Chen, and S. Jin, Morphology and kinetics evolution of $\gamma$’ phase with increased volume fraction in ni–al alloys, Materials Chemistry and Physics, 217 (2018), pp. 23–30.
  * [25] Y. Lu, D. Jia, T. Hu, Z. Chen, and L. Zhang, Phase-field study the effects of elastic strain energy on the occupation probability of cr atom in ni–al–cr alloy, Superlattices and Microstructures, 66 (2014), pp. 105–111.
  * [26] T. Massalski, Phase diagrams in materials science, Metallurgical Transactions A, 20 (1989), pp. 1295–1323.
  * [27] Y. Mishin, Atomistic modeling of the $\gamma$ and $\gamma$’-phases of the ni–al system, Acta Materialia, 52 (2004), pp. 1451–1467.
  * [28] R. Mueller and D. Gross, 3d simulation of equilibrium morphologies of precipitates, Computational materials science, 11 (1998), pp. 35–44.
  * [29] Z. Qiao, Z. Zhang, and T. Tang, An adaptive time-stepping strategy for the molecular beam epitaxy models, SIAM Journal on Scientific Computing, 33 (2011), pp. 1395–1414.
  * [30] Y. Saad and M. H. Schultz, Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput., 7 (1986), pp. 856–869.
  * [31] N. Saunders and A. P. Miodownik, CALPHAD (calculation of phase diagrams): a comprehensive guide, Elsevier, 1998.
  * [32] J. Shen, C. Wang, X. Wang, and S. M. Wise, Second-order convex splitting schemes for gradient flows with ehrlich–schwoebel type energy: application to thin film epitaxy, SIAM Journal on Numerical Analysis, 50 (2012), pp. 105–125.
  * [33] J. Shen, J. Xu, and J. Yang, A new class of efficient and robust energy stable schemes for gradient flows, SIAM Review, 61 (2019), pp. 474–506.
  * [34] J. Shen and X. Yang, Numerical approximations of allen-cahn and cahn-hilliard equations, Discrete & Continuous Dynamical Systems, 28 (2010), p. 1669.
  * [35] T. Smith, Y. Rao, Y. Wang, M. Ghazisaeidi, and M. Mills, Diffusion processes during creep at intermediate temperatures in a ni-based superalloy, Acta Materialia, 141 (2017), pp. 261–272.
  * [36] X. Su, Q. Xu, R. Wang, Z. Xu, S. Liu, and B. Liu, Microstructural evolution and compositional homogenization of a low re-bearing ni-based single crystal superalloy during through progression of heat treatment, Materials & Design, 141 (2018), pp. 296–322.
  * [37] P. Turchi, Phase transformations and evolution in materials, (2000).
  * [38] C. Wagner, Theory of precipitate change by redissolution, Z. Elektrochem, 65 (1961), pp. 581–591.
  * [39] Y. Wang, D. Banerjee, C.-C. Su, and A. Khachaturyan, Field kinetic model and computer simulation of precipitation of l12 ordered intermetallics from fcc solid solution, Acta materialia, 46 (1998), pp. 2983–3001.
  * [40] X. Wu, Y. Li, M. Huang, W. Liu, and Z. Hou, Precipitation kinetics of ordered $\gamma$’ phase and microstructure evolution in a nial alloy, Materials Chemistry and Physics, 182 (2016), pp. 125–132.
  * [41] V. Yadav, L. Vanherpe, and N. Moelans, Effect of volume fractions on microstructure evolution in isotropic volume-conserved two-phase alloys: A phase-field study, Computational Materials Science, 125 (2016), pp. 297–308.
  * [42] M. Yang, H. Wei, J. Zhang, Y. Zhao, T. Jin, L. Liu, and X. Sun, Phase-field study on effects of antiphase domain and elastic energy on evolution of $\gamma$’ precipitates in nickel-based superalloys, Computational Materials Science, 129 (2017), pp. 211–219.
  * [43] X. Yang, Linear, first and second-order, unconditionally energy stable numerical schemes for the phase field model of homopolymer blends, Journal of Computational Physics, 327 (2016), pp. 294–316.
  * [44] C. Y. Ying Wei and J. Huang, Parallel energy-stable phase field crystal simulations based on domain decomposition methods, arXiv:1703.01202v1, (2017).
  * [45] Y. Zhang, C. Yang, and Q. Xu, Numerical simulation of microstructure evolution in ni-based superalloys during p-type rafting using multiphase-field model and crystal plasticity, Computational Materials Science, 172 (2020), p. 109331.
  * [46] Z. Zhang, Y. Ma, and Z. Qiao, An adaptive time-stepping strategy for solving the phase field crystal model, J. Comput. Phys., 249 (2013), pp. 204–215.
  * [47] J. Zhao, Q. Wang, and X. Yang, Numerical approximations for a phase field dendritic crystal growth model based on the invariant energy quadratization approach, International Journal for Numerical Methods in Engineering, 110 (2017), pp. 279–300.
  * [48] J. Zhu, L.-Q. Chen, J. Shen, and V. Tikare, Coarsening kinetics from a variable-mobility cahn-hilliard equation: Application of a semi-implicit fourier spectral method, Physical Review E, 60 (1999), p. 3564.
  * [49] J. Zhu, Z. Liu, V. Vaithyanathan, and L. Chen, Linking phase-field model to calphad: application to precipitate shape evolution in ni-base alloys, Scripta Materialia, 46 (2002), pp. 401–406.
  * [50] J. Zhu, T. Wang, A. Ardell, S. Zhou, Z. Liu, and L. Chen, Three-dimensional phase-field simulations of coarsening kinetics of $\gamma$’ particles in binary ni–al alloys, Acta materialia, 52 (2004), pp. 2837–2845.