# A crystal symmetry-invariant Kobayashi–Warren–Carter grain boundary model
and its implementation using a thresholding algorithm

Jaekwang Kim  Matt Jacobs  Stanley Osher  Nikhil Chandra Admal
<EMAIL_ADDRESS>Department of Mechanical Science and Engineering,
University of Illinois Urbana–Champaign  Department of Mathematics, University
of California Los Angeles

###### Abstract

One of the most important aims of grain boundary modeling is to predict the
evolution of a large collection of grains in phenomena such as abnormal grain
growth, coupled grain boundary motion, and recrystallization that occur under
extreme thermomechanical loads. A unified framework to study the coevolution
of grain boundaries with bulk plasticity has recently been developed by Admal
et al. (2018), which is based on modeling grain boundaries as continuum
dislocations governed by an energy based on the Kobayashi–Warren–Carter (KWC)
model (Kobayashi et al., 1998, 2000). While the resulting unified model
demonstrates coupled grain boundary motion and polygonization (seen in
recrystallization), it is restricted to grain boundary energies of the
Read–Shockley type, which applies only to small misorientation angles. In
addition, the implementation of the unified model using finite elements
inherits the computational challenges of the KWC model that originate from the
singular diffusive nature of its governing equations. The main goal of this
study is to generalize the KWC functional to grain boundary energies beyond
the Read—Shockley-type that respect the bicrystallography of grain boundaries.
The computational challenges of the KWC model are addressed by developing a
thresholding method that relies on a primal dual algorithm and the fast
marching method, resulting in an $\mathcal{O}(N\log N)$ algorithm, where $N$
is the number of grid points. We validate the model by demonstrating the
Herring angle and the von Neumann–Mullins relations, followed by a study of
the grain microstructure evolution in a two-dimensional face-centered cubic
copper polycrystal with crystal symmetry-invariant grain boundary energy data
obtained from the covariance grain boundary model of Runnels et al. (2016a,
b).

###### keywords:

A. phase field model , microstructures , B. motion by curvature , constitutive
behavior , polycrystalline materials

††journal: Journal of Computational Materials Science

##  1 Introduction

Most metals and ceramics exist as polycrystals, which are aggregates of single
crystal grains stacked together along grain boundaries. The microstructure of
a polycrystal is often characterized by the orientation distribution of its
grains commonly referred to as texture. The macroscopic properties of
polycrystals, which include yield strength, resistance to creep, fatigue,
thermal and magnetic properties, are strongly influenced by texture. Grain
boundary engineering refers to the strategy of enhancing the properties of
polycrystalline materials by transforming the grain boundary character
distribution to a desired state using thermomechanical processes (Watanabe,
2011). Mapping the microstructure-property relationship, and modeling the
evolution of microstructure under various manufacturing processes are
fundamental open problems relevant to grain boundary engineering.

The evolution of grain boundaries is driven by a long list of thermodynamic
forces, of which surface tension plays a central role. For instance, grain
boundaries in the isotropic Mullins’ model (Mullins, 1956) are driven by their
excess surface energy resulting in motion by curvature with velocity given by

$v=-m\gamma\kappa,$ (1.1)

where $\kappa$ and $\gamma$ denote curvature and misorientation-dependent
energy density of the grain boundary respectively, while $m$ represents a
constant mobility. More generally, an anisotropic grain boundary evolution
arises from the dependence of $m$ and $\gamma$ on the grain boundary character
defined by the five macroscopic degrees of freedom, which represent the
misorientation and the inclination of the grain boundary.111Under an
anisotropic energy density $\gamma$ that depends on the inclination $n$ of a
grain boundary, (1.1) transforms as
$v=-m\kappa(\gamma+\partial^{2}\gamma/\partial n^{2})$. Recent advances in the
development of accurate interatomic potentials have enabled us to build and
refine the grain boundary energy and mobility landscapes as functions of the
grain boundary character (Chen et al., 2020; Runnels et al., 2016a, b; Olmsted
et al., 2009; Bulatov et al., 2013). The grain boundary energy landscape
reflects the symmetry of a bicrystal, and understanding the role of crystal
symmetry contributes enormously towards characterizing grain microstructure
evolution. Moreover, precisely identifying the grain boundary character
distribution responsible for phenomena such as abnormal grain growth and
recrystallization remains an open problem in materials science. This motivates
us to undertake simulations of large ensembles of appropriately sampled
polycrystals to discover lower-order statistical models for grain
microstructure evolution. _The goal of this paper is to develop a lightweight
model for motion by curvature in the presence of a misorientation-dependent
grain boundary energy density, that can be implemented using an ultrafast
algorithm._

While motion by curvature is the simplest description of grain boundary
evolution, experiments (Rollett, 2018; Barmak et al., 2013) and atomistic
simulations (Upmanyu et al., 1998; Janssens et al., 2006) demonstrate that
surface tension alone is not the dominant force. Molecular dynamics (MD)
simulations have revealed that as grain boundaries evolve, they plastically
deform the underlying material resulting in lattice distortions that give rise
to additional forces on the grain boundary. In order to include the effect of
grain boundary plasticity, recent mesoscale models (Admal et al., 2018) have
focused on the coevolution grain boundaries and deformation. Evidently, such
models subsume motion by curvature as a special case, and are computationally
more expensive, which is another motivation for us to seek an ultrafast
algorithm. While the focus of this paper is motion by curvature, we ensure
that our model is amenable to generalizations that include grain boundary
plasticity.

Models for grain microstructure evolution can be broadly classified into three
categories: probabilistic, diffuse-interface, and sharp-interface models. An
example of a probabilistic grain growth model is the Monte-Carlo Potts model
(Anderson et al., 1958, 1989; Mendelev and Srolovitz, 2002; Upmanyu et al.,
2002; Yang et al., 2000), wherein a polycrystal is described using points in a
lattice, which are allocated to different grains. A grain boundary is
implicitly defined by adjacent lattice points that belong to different grains.
Evolution of the microstructure is carried out stochastically through random
jumps of boundaries in thermodynamically favorable directions. While the
advantage of the Mote-Carlo Potts model lies in the simplicity of its
implementation, it relies on heuristic rules that do not have a thermodynamic
basis.

In sharp interface models (Mullins, 1956; Hillert, 1965; Allen and Cahn,
1979), grain boundaries are modeled as surfaces that evolve according to
motion by curvature given in (1.1). Methods to implement (1.1) rely on either
implicitly or explicitly tracking the moving grain boundaries. For instance,
front tracking methods (Frost et al., 1990, 1988; Kinderlehrer et al., 2004,
2006) describe grain boundaries in two dimensions as line segments along with
their connectivity. Such a description breaks down at critical events
including disappearance of shrinking grains and topological changes due to
merging of grain boundaries. Therefore, front tracking methods are supplied
with additional rules to redefine the connectivity of line segments to
describe such critical events. The level set method (Zhao et al., 1996; Fausty
et al., 2018) addresses the above limitation using an implicit representation.
Each grain is described by a function that is positive within, and negative
outside the grain, which implies the zero-valued isosurface describes the
interface surrounding the grain. An implicit representation can describe
topological changes without any additional rules. Yet, the main disadvantage
of the level set method is that it does not extend to handle surfaces with
self intersection and junctions, which occur in polycrystals. More recently, a
thresholding method commonly referred to as the Merriman–Bence–Osher (MBO)
(Merriman et al., 1992) scheme and its generalization by Esedoḡlu and Otto
(2015) have been shown to simulate grain kinetics very efficiently in addition
to its ability to predict grain nucleation. The level set and the MBO methods
are memory intensive as they use as many functions as the number of grains in
order to describe a polycrystal. For example, a description of a 3D
polycrystal with $10,000$ grains on a 256$\times$256$\times$256 grid requires
$\sim 1$ TB of memory. Thus, for a large scale simulation, an additional
numerical technique is devised to employ a level set function that
approximates a large subset of spatially separated grains (Elsey et al., 2009,
2011). The thresholding method of Esedoḡlu and Otto (2015), which is first-
order accurate in time, has recently been extended by Zaitzeff et al. (2020)
to a second-order method that is unconditionally energy stable. In general,
all sharp-interface models can be incorporated with misorientation dependent
grain boundary energy densities and mobilities. However, including inclination
dependence is more challenging, and recent works by Basak and Gupta (2014);
Hallberg and Bulatov (2019); Joshi et al. (2020) have addressed this
challenge.

In diffuse interface models (Jokisaari et al., 2017), a polycrystal is defined
using functions called phase fields, which are constant in the interior of the
grains. The regions where the gradients of phase fields are non-zero are
identified as diffused grain boundaries, which have a characteristic width. A
numerical implementation of a diffuse interface model requires a grid that is
refined enough to resolve the width of the grain boundary. Therefore, diffuse
interface models are computationally more expensive than their sharp-interface
counterparts such as the level set and MBO methods. The multi phase field
(MPF) model (Chen, 2002; Hirouchi et al., 2012; Steinbach, 2009), and the
Kobayashi–Warren–Carter (KWC) model (Kobayashi et al., 1998, 2000; Warren et
al., 2003) are two examples of diffuse-interface models for grain boundaries.

The main advantage of the MPF model lies in the simplicity of its construction
to include misorientation dependent grain boundary energies and mobilities.
Similar to the level set and the MBO methods, a naive implementation of the
MPF method would use as many phase fields as the number of grains, resulting
in an excessive use of computational memory. Since, at any point in the
domain, only a few order parameters would be non-zero, recent implementations
(Fan et al., 2002; Permann et al., 2016) of the MPF model allow multiple
grains which do not share a common boundary to share the same order parameter.
A _grain remapping_ algorithm is used to strategically remap an order
parameter shared by two distant grains when they approach close to each other.
Recent advances (Ribot et al., 2019; Moelans et al., 2008; Kim et al., 2014)
in MPF models explore the full anisotropy of grain boundary energy which
consists of misorientation and inclination dependence.222Grain boundary energy
as a function of inclination is typically non-convex. For grain boundary
models that incorporate inclination dependence to be well-posed, they must
include curvature-dependent energy densities resulting in a higher-order model
which adds to the computationally intensive nature of the MPF model.

In contrast, the KWC model describes a two-dimensional polycrystal using only
two order parameters - one for structural order $\eta$ ranging from 0
(disordered phase) to 1 (crystalline state), and the other for crystal
orientation field $\theta$. The elegance of the KWC model is offset by the
severe restriction it imposes on the grain boundary energy. The energy
functional of the KWC model limits the dependence of grain boundary energy on
misorientation angle to a Read–Shockley-type (Read and Shockley, 1950) that
does not respect the crystal symmetry. In addition, the singular diffusive
nature of the KWC model results in stiff equations that are computationally
expensive to solve.

Recognizing the elegance of the KWC model, we formulate a generalization of
the KWC model that can incorporate arbitrary misorientation-dependent grain
boundary energies. In addition, we design a thresholding method that addresses
the challenge of solving the singular diffusive equation of the KWC model. The
resulting model inherits the memory efficiency of the original KWC model,
while having significantly more computational efficiency compared to
conventional numerical methods such as finite element and finite difference.

The paper is structured as follows. In Section 2, we discuss the role of
crystal symmetry on the anisotropy of grain boundary energies, and review the
original KWC model and its limitations. Then, we propose a generalization of
the KWC model to incorporate grain boundary energies beyond the Read–Shockley
type. In Section 3, we design a thresholding algorithm to implement the
generalized KWC model. Numerical experiments to validate and demonstrate the
computational efficiency of the thresholding method are discussed in Section
4. Finally, we summarize and conclude with a description of future directions.

##  2 Grain boundary energy and the Kobayashi–Warren–Carter model

Figure 1: A plot of grain boundary energy density as a function of
misorientation angle of a $[110]$ symmetric-tilt grain boundary in fcc copper,
computed using molecular dynamics Holm et al. (2010); Bulatov et al. (2014).
Misorientations corresponding to low energy $\Sigma$ boundaries are marked on
the upper axis.

A grain boundary is characterized by five macroscopic degrees of freedom where
three degrees represent a rotation associated with the misorientation between
the two grains, and the remaining two degrees correspond to the inclination of
the grain boundary. More precisely, the grain boundary character space is
given by the topological space $\mathcal{T}=SO(3)\times SO(3)/SO(2)$, where
$SO(n)$ is the special orthogonal group in $n$ dimensions. Grain boundaries
are equipped with a surface energy density, which is defined as a function on
$\mathcal{T}$. An energy density that is constant is referred to as an
_isotropic_ , and _anisotropic_ otherwise. Fig. 1 shows a plot of grain
boundary energy density as a function of misorientation angle for a $[110]$
symmetric-tilt grain boundary in face-centered cubic (fcc) copper, calculated
using molecular dynamics simulations (Holm et al., 2010; Bulatov et al.,
2014). Since the symmetry of an fcc lattice ensures that the energy of a
$[110]$ symmetric-tilt grain boundary is symmetric about the $180^{\circ}$
misorientation angle, Fig. 1 shows a plot of energy vs misorientation angles
up to $180^{\circ}$. In addition, $\gamma$ exhibits local minima at certain
misorientations, marked as $\Sigma_{3}$ and $\Sigma_{11}$ in Fig. 1, due to an
enhanced lattice matching (Runnels et al., 2016a, b; Wolf, 1990) between the
two adjoining grains. Recent efforts (Mason and Patala, 2019; Bulatov et al.,
2013; Runnels et al., 2016a, b; Olmsted et al., 2009; Kim et al., 2014) by
materials scientists in characterizing the grain boundary character space, and
parametrizing grain boundary energy using data from atomistic simulations and
experiments, bring us closer to developing an atomistically informed mesoscale
model for grain boundaries.

The motion of grain boundaries driven by surface tension to decrease the
interfacial energy is a defining characteristic of various grain
microstructure models such as the Mullins model (Mullins, 1956), and its
diffuse-interface counterparts such as the multiphase field and the KWC
models. The resulting grain boundary motion, commonly referred to as _motion
by curvature_ , in the presence of an anisotropic energy density has been
shown to have a considerable effect on grain statistics (Barmak et al., 2013)
leading to changes in the macroscopic properties of materials. In order to
explore the structure-property relationship, recent research efforts have
focused on developing ultrafast algorithms to simulate grain boundary
evolution in large polycrystals in the presence of atomistically derived
anisotropic grain boundary energies.

In this paper, we recognize the simplicity of the phase field model of
Kobayashi, Warren and Carter, and show that it can be generalized and
implemented using a thresholding method resulting in an ultrafast algorithm
for grain boundary evolution.

###  2.1 The KWC model

(a)

(b)

Figure 2: Results on the one-dimensional steady state solution of the original
KWC model describing a flat grain boundary. 2(a)) A steady state analytical
solution of the KWC model for a given misorientation. 2(b)) Variation of the
grain boundary energy as a function of misorientation in the KWC model.

The Kobayashi–Warren–Carter (KWC) model (Kobayashi et al., 1998, 2000; Warren
et al., 2003) is a dual-phase field model to study grain evolution in
polycrystalline materials. In this model, an arbitrary polycrystal in two
dimensions is described using only two order parameters $\eta$ and $\theta$.
This is one of the main advantages of employing the KWC model as opposed to
the multiphase field model, which uses as many order parameter as the number
of grains. The order parameter $\eta$ ranges from $0$, which signifies
disorder, to $1$ that describes crystalline order. On other hand, the order
parameter $\theta$ describes the orientation of the grains.333In three
dimensions, the order parameter $\theta$ is replaced by a rotation tensor.

The KWC model is governed by a free energy functional given by

$\mathcal{W}[\eta,\theta]=\int_{\Omega}\left[\frac{f(\eta)}{\epsilon}+\frac{\epsilon}{2}|\nabla\eta|^{2}+g(\eta)|\nabla\theta|+\frac{\epsilon}{2}|\nabla\theta|^{2}\right]\;dV,$
(2.1)

where

$f=\frac{(1-\eta)^{2}}{2}$ (2.2)

is a single-well potential with minimum at $\eta=1$, and

$g=-\ln(1-\eta)$ (2.3)

is an increasing function.444 The choice of the logarithmic function for $g$
is supported by the work of Alicandro et al. (1999), which shows that the KWC
functional converges (in the sense of $\Gamma$-convergence) to a surface
energy function when $g(1)=\infty$ (See Theorem 4.1 in Alicandro et al.
(1999)).  The functional in (2.1) is defined for all functions $\eta$ and
$\theta$ in the Hilbert space $H^{1}(\Omega)$.555$H^{1}(\Omega)$ denotes the
set of all functions on $\Omega$ whose first derivatives are square
integrable. The constant $\epsilon>0$ is a dimensionless scaling parameter
that determines the thickness of the grain boundary region (Lobkovsky and
Warren, 2001). The evolution equations for the order parameters, assuming a
gradient descent (with respect to the $L^{2}$-norm) of $\mathcal{W}$, are
obtained as

$\displaystyle\epsilon b_{\eta}\dot{\eta}$
$\displaystyle=\epsilon\Delta\eta-\frac{f^{\prime}(\eta)}{\epsilon}-g^{\prime}(\eta)|\nabla\theta|,$
(2.4a) $\displaystyle\epsilon b_{\theta}\dot{\theta}$
$\displaystyle=\nabla\cdot\left[\epsilon\nabla\theta+g(\eta)\frac{\nabla\theta}{|\nabla\theta|}\right],$
(2.4b)

where $b_{\eta}$ and $b_{\theta}$ are the inverse mobilities corresponding to
respective order parameters.666The KWC model was originally developed to
simultaneously model grain rotation _and_ grain boundary motion. The model can
be specialized to demonstrate _only_ grain boundary motion by enforcing zero
mobility for $\theta$ in the grain interior. This can be achieved by a
constant $b^{\phi}$, and a $\phi$-dependent $b^{\theta}$ (Dorr et al., 2010).

A one-dimensional steady state solution of (2.4) under Dirichlet boundary
conditions is plotted in Fig. 2(a). The value of $\eta<1$ in a neighborhood of
the grain boundary suggests a loss of crystalline order, and the orientation
$\theta$ is constant in the interior of the grains, and has a non-zero
gradient in a finite thickness around the grain boundary. In the limit
$\epsilon\to 0$, Lobkovsky and Warren (2001) have shown that the evolution
equations in (2.4) result in the shrinking of the grain boundary thickness
converging to the Mullins model.

Below, we briefly summarize the role of each term appearing in the KWC
functional. We refer the reader to Warren et al. (2003) for a more detailed
description, and Admal et al. (2019) for a generalization of the KWC model to
three dimensions. The function $f$ drives $\eta(\bm{x})$ towards $1$, while
the coupled term $g(\eta)|\nabla\theta|$ tends to decrease $\eta$ in a
neighborhood of the grain boundary. In addition, the coupled term tends to
localize the jump in $\theta$, while $|\nabla\theta|^{2}$ has a tendency to
diffuse it, resulting in a regularized step function for $\theta$. It is
interesting to note that in the absence of the $|\nabla\theta^{2}|$ term, the
steady state solution for $\theta$ is a pure step function, resulting in a
model with a blend of sharp- and diffuse-interface characteristics, i.e. while
$\theta$ is sharp, $\eta$ is diffused. Moreover, grain boundaries cease to
evolve in the absence of $|\nabla\theta|^{2}$ term (Lobkovsky and Warren,
2001). In other words, $|\nabla\theta|^{2}$ in the KWC model has a dual role
of not only regularizing $\theta$ but also rendering non-zero mobility to the
grain boundaries.

The grain boundary energy $\gamma$, as a function of misorientation angle,
predicted by the KWC model is of the Read–Shockley-type, as shown in Fig.
2(b). This is in contrast to the experimentally observed grain boundary
energies shown in Fig. 1. Despite the elegance of the KWC model in describing
polycrystals with only two order parameters, its restriction to Read–Shockley
type grain boundary energies is a major limitation compared to the flexibility
of incorporating arbitrary grain boundary energies into the multiphase field
model. The above-stated limitation is one of the main motivation for us to
seek a new formulation of the KWC model to incorporate arbitrary
misorientation-dependent grain boundary energies that respect the
bicrystallography of grain boundaries.

###  2.2 A crystal symmetry-invariant KWC model

In this section, we formulate a new KWC model that can incorporate arbitrary
misorientation-dependent grain boundary energies.

We begin with the KWC functional without the $|\nabla\theta|^{2}$ term. From
Section 2.1, recall that in the absence of the $|\nabla\theta|^{2}$ term, the
steady state solution for $\theta$ is a step function with the discontinuity
occurring at the grain boundary. In one-dimension, since a discontinuous
$\theta$ is not in $H^{1}(\Omega)$, the minimizer of $\mathcal{W}$, with
$|\nabla\theta|^{2}$ absent, is not attained. This observation motivates us to
redefine the domain of the modified KWC functional such that $\theta$ belongs
to the space of piecewise constant functions, as opposed to $H^{1}(\Omega)$,
and this enables us to simplify the functional as

$\mathcal{W}[\eta,\theta]=\int_{\Omega}\left[\frac{(1-\eta)^{2}}{2\epsilon}+\frac{\epsilon}{2}|\nabla\eta|^{2}\right]\,dV-\int_{\mathcal{S}}\ln{(1-\bar{\eta})}[\\![\theta]\\!]\,dS,$
(2.5)

where $\bar{\eta}:\mathcal{S}\to\mathbb{R}$ is the restriction of $\eta$ to
the jump set $\mathcal{S}$ of $\theta$, which represents the union of all
grain boundaries. The steady state solution, given in (A.8), corresponding to
a one-dimensional bicrystal governed by (2.5), and the resulting grain
boundary energy as a function of the misorientation,

$\gamma(\llbracket\theta\rrbracket)=\frac{\llbracket\theta\rrbracket}{2}\left(1-2\ln\left[\frac{\llbracket\theta\rrbracket}{2}\right]\right)$
(2.6)

are derived in A. From (2.6), it is clear that the grain boundary energy is of
a Read–Shockley-type, which does not respect the crystal symmetry.

The above observation leads us to the following generalization of the KWC
functional

$\mathcal{W}^{\rm
G}[\eta,\theta]=\int_{\Omega}\left[\frac{(1-\eta)^{2}}{2\epsilon}+\frac{\epsilon}{2}|\nabla\eta|^{2}\right]\,dV+\int_{\mathcal{S}}g(\bar{\eta})\mathcal{J}\left([\\![\theta]\\!]\right)\,dS,$
(2.7)

which is defined for all $\eta\in H^{1}(\Omega)$, and piecewise constant
functions $\theta$; and $\mathcal{J}$ is an even function of the jump in
orientation. Under this new formulation, the grain boundary energy function
modifies as

$\displaystyle\gamma^{\rm G}([\\![\theta]\\!])$
$\displaystyle=(1-\overline{\eta})^{2}-\ln{(1-\overline{\eta})}\mathcal{J}([\\![\theta]\\!])$
(2.8)
$\displaystyle=\frac{\mathcal{J}(\llbracket\theta\rrbracket)}{2}\left(1-2\ln\left[\frac{\mathcal{J}(\llbracket\theta\rrbracket)}{2}\right]\right),$

where $\bar{\eta}$ is the value of the stead-state solution on the grain
boundary, given implicitly in terms of
$\mathcal{J}(\llbracket\theta\rrbracket)$ as777 The analog of (2.9) in the
original KWC model is (A.7), whose derivation is shown in A.

$2(1-\overline{\eta})^{2}=\mathcal{J}([\\![\theta]\\!]).$ (2.9)

Inspired from the terminology in dislocations, we refer to
$\mathcal{J}(\llbracket\theta\rrbracket)$ as the core energy.

Figure 3: A plot (in green) of the core energy
$\mathcal{J}(\llbracket\theta\rrbracket)$ calculated to match $\gamma^{G}$
(2.8) to $\gamma^{\rm cov}$ (B.9) (in red) for $[110]$ symmetric-tilt grain
boundaries in fcc copper. $\gamma^{G}$ is identical to $\gamma^{\rm cov}$ by
construction.

From (2.8), it is clear that by appropriately constructing the core energy
$\mathcal{J}$, we can arrive at a $\gamma^{\rm G}$ that faithfully represents
the grain boundary energy and symmetry of the bicrystal. In other words, the
crystal symmetry of the new KWC model is inherited from the core energy. For
illustration, we consider the energy $\gamma^{\rm cov}$ of a $[110]$ symmetric
tilt grain boundary in face-centered cubic (fcc) copper, shown as red data
points in Fig. 3. $\gamma^{\rm cov}$ is computed using the covariance model
developed by Runnels et al. (2016a, b), wherein it is defined as the
covariance of the two lattices adjoining the grain boundary. For completeness,
in B, we describe the covariance model of grain boundary energy, the procedure
to compute it, and list the parameters used to arrive at the data plotted in
Fig. 3. We solve for $\mathcal{J}$ in (2.8) using the Newton’s method such
that the resulting $\gamma^{\rm G}=\gamma^{\rm cov}$ at each data point. In
other words, the grain boundary energy of the resulting KWC model is identical
to $\gamma^{\rm cov}$ by construction. Fig. 3 shows a plot of the solution
$\mathcal{J}$ in green, and highlights the common positions of the local
minimizers of $\mathcal{J}$ and $\gamma^{\rm cov}$. Typically, the number of
Newton iterations for convergence error (absolute value of the difference
between $\gamma^{\rm G}$ and $\gamma^{\rm cov}$) of $10^{-6}\rm[J/m^{2}]$ is
less than 20. While the modification of the KWC functional from (2.5) to (2.7)
allows us to model arbitrary misorientation-dependent grain boundary energies,
the absence of $|\nabla\theta|^{2}$ term in (2.7) renders the grain boundaries
immobile as mentioned in Section 2.1.888The evolution of $\theta$ by the
gradient descent of $\mathcal{W}^{\rm G}$ results in pure rotation while the
position of the grain boundaries remains fixed. In what follows, we address
this shortcoming by devising a thresholding method to move grain boundaries by
evolving the piecewise-constant $\theta$.

##  3 Grain boundary motion in the new KWC model

In this section, we present our approach in which we alternate between
evolving $\eta$ and $\theta$ to evolve a polycrystal governed by
$\mathcal{W}^{\rm G}$. The order parameter $\eta$ is solved in the following
minimization problem for a given $\theta$

$\eta^{*}=\operatornamewithlimits{arg\;min}_{\begin{subarray}{c}\eta\in
H^{1}(\Omega)\\\ \partial\eta/\partial
n|_{\partial\Omega}=0\end{subarray}}\mathcal{W}^{\rm G}[\eta,\theta].$ (3.1)

Next, the orientation field $\theta$ is evolved using a thresholding rule
described in the next section.

In order to solve for $\eta$ in (3.1), we note that since $g(\eta)\to\infty$
as $\eta\to 1$, the functional $W^{\rm G}$ is non-smooth in $\eta$, which
makes the Newton’s method not viable. Therefore, we use a primal-dual method
recently developed by Jacobs et al. (2019), which has a
$\mathcal{O}(\frac{1}{e}N\log N)$ complexity, where $e$ is the error in the
numerical solution to (3.1), and $N$ is the grid size. See C for a more
detailed description of the primal dual method.

Next, we develop a thresholding rule to evolve $\theta$ for a fixed $\eta^{*}$
obtained in (3.1). The alternate use of the primal dual method and the
thresholding rule at every time step constitutes our approach to evolving the
grain boundaries.

###  3.1 The thresholding rule

A thresholding method is a sequence of simple rules, executed every time step
to reinitialize the order parameter, such that its evolution describes the
motion of a grain boundary.

The original idea of using a thresholding method to evolve grain boundaries
goes back to the work (Merriman et al., 1992) of Merriman, Bence and Osher
(MBO) wherein, similar to the multiphase field model, grains in a polycrystal
are described using as many order parameters, with the caveat that the order
parameters are piecewise-constant implying a sharp interface. An order
parameter in the MBO method is evolved based on a two-step thresholding scheme
— a convolution of the order parameter with a Gaussian kernel followed by a
trivial thresholding — resulting in motion by curvature.

The MBO method has recently been generalized by Esedoḡlu and Otto (2015) to a
variational model, referred to as the _Gaussian kernel method_. In the MBO and
the Gaussian kernel methods, there are as many characteristic functions as the
number of distinct grains. While the end goal of the KWC model is also to
describe motion by curvature, it is markedly different from the Gaussian
kernel method as it uses only two order parameters to represent a polycrystal.
Therefore, it does not require additional techniques to address the memory
intensive nature of a naive implementation of the MBO/Gaussian kernel methods.
More importantly, the KWC model lends itself to further generalizations which
include the modeling of grain rotation. Therefore, the goal here is to seek a
thresholding algorithm to implement the KWC model.

In this section, we design a thresholding rule for $\theta$ that results in
motion by curvature. We first recall that $\theta$ is a piecewise-constant
field with a finite range of orientations. This implies, that a thresholding
rule for $\theta$ reassigns $\theta(\bm{x})$, for each point
$\bm{x}\in\Omega$, to one of the possible orientations. Our thresholding rule
originates from the observation that _the asymmetry of $\eta$ in the
neighborhood of a grain boundary characterizes its curvature._ Below, we
explicitly identify this asymmetry before describing our thresholding rule.

Figure 4: Level sets of the distance function $l(\bm{x})$ in a neighborhood of
$\bm{x}_{0}$. The grain boundary $\Gamma$ is depicted as a solid curve, and
the dashed curves correspond to the level sets $l(x)/\epsilon=\pm\tilde{l}$.

First, we note that the steady-state solution for $\eta$, derived in (A.8) for
a flat interface (zero curvature), is symmetric about the grain boundary.
Next, we derive an approximate form for $\eta$ in the presence of a non-zero
curvature, and show the dependence of asymmetry on the curvature. Let $\Gamma$
denote a grain boundary with a non-zero curvature that separates two grains
with orientations $\theta_{\rm L}$ and $\theta_{\rm R}$. We postulate that for
a $\Gamma$ with a small curvature such that $\epsilon\kappa\ll 1$, and for a
$\bm{x}_{0}\in\Gamma$ away from a triple junction, the solution $\eta^{*}$ to
(2.4a) is approximated in a small neighborhood of $\bm{x}_{0}$ as

$\eta^{*}(\bm{x})\approx u\left(\frac{l(\bm{x})}{\epsilon}\right),$ (3.2)

where $l(\bm{x})$ is the signed distance function from $\Gamma$ to $\bm{x}$,
as shown in Fig. 4. In other words, we assume that $\eta$ only depends on the
radial coordinate. Away from the grain boundary, the solution to the
minimization problem in (3.1) satisfies the equation

$\epsilon\Delta\eta^{*}-\frac{(\eta^{*}-1)}{\epsilon}=0.$ (3.3a)

The above equation can be simplified by using a local coordinate system
$\bm{x}=(\tilde{l},s)$, where $\tilde{l}=l(\bm{x})/\epsilon$ is the scaled
radial coordinate, and $s$ is the distance measured along $\Gamma$ between
$\bm{x}_{0}$ and the perpendicular projection of $\bm{x}$ on $\Gamma$. In this
coordinate system, we note that
$\triangle\eta(\bm{x})=u^{\prime\prime}/\epsilon^{2}+\kappa
u^{\prime}/\epsilon$, where $\kappa(\tilde{l})$ is the curvature of the
coordinate line $\\{\bm{x}\in
B_{r}(\bm{x}_{0}):l(\bm{x})=\tilde{l}\epsilon\\}$. Therefore, (3.3a)
simplifies as

$u^{\prime\prime}(\tilde{l})+\epsilon\kappa
u^{\prime}(\tilde{l})-u(\tilde{l})+1=0.$ (3.4)

Assuming $\kappa(\tilde{l})=\kappa(0)$, we obtain the following closed form
solution to (3.4):

$u(\tilde{l})=1+C_{1}\exp{\left[-\tilde{l}\left(\frac{\epsilon\kappa+\sqrt{4+\epsilon^{2}\kappa^{2}}}{2}\right)\right]}+C_{2}\exp{\left[-\tilde{l}\left(\frac{\epsilon\kappa-\sqrt{4+\epsilon^{2}\kappa^{2}}}{2}\right)\right]},$
(3.5)

where the constants $C_{1}$ and $C_{2}$ are determined using the boundary
conditions $u(\pm\infty)=1$.999The boundary conditions are interpreted in the
limit $\epsilon\to 0$, which results in the boundary conditions
$l/\epsilon=\pm\infty$ for the scaled radial coordinate. Subsequently, the
solution can be further approximated101010Here, we use the approximation
$\sqrt{4+\epsilon^{2}\kappa^{2}}\approx
2+\mathcal{O}(\epsilon^{2}\kappa^{2})$. under the assumption that both
$\epsilon$ and $\epsilon\kappa$ are small, resulting in

$u(\tilde{l})=\begin{cases}1+(u(0)-1)e^{-(1+0.5\epsilon\kappa)\tilde{l}}&\text{
if }\tilde{l}>0,\\\ 1+(u(0)-1)e^{(1-0.5\epsilon\kappa)\tilde{l}}&\text{
otherwise.}\end{cases}$ (3.6)

The asymmetry of $u$ is apparent from (3.6) by noting that in the presence of
a positive curvature, the rate at which $u\to 1$ as $\tilde{l}\to+\infty$ is
greater than when $\tilde{l}\to-\infty$.

Figure 5:  A plot of $\eta^{*}$ in a small neighborhood of $x_{0}$ (see Fig.
4) is shown in blue, while $(1-\eta^{*})^{2}$ is shown in red. The asymmetry
of $u$ around $x_{0}$ due to curvature $\kappa$ is characterized by the
position $x$ at which the two areas shown in yellow and green regions are
equal. The position $x$ is given in terms of $l_{0}$
($(x-x_{0})=l_{0}/\epsilon$), which is the solution of (3.7).

The asymmetry of $u$ forms the foundation of our thresholding scheme, which is
designed to reassign the values of $\theta$ in the neighborhood of the grain
boundaries resulting in a motion by curvature. To design a thresholding rule,
we identify a unique $l=l_{0}$, such that

$\int^{l_{0}}_{-\infty}(1-u(l/\epsilon))^{2}\,dl=\int^{+\infty}_{l_{0}}(1-u(l/\epsilon))^{2}\,dl.$
(3.7)

The two integrals in (3.7) are depicted as equal areas under the yellow and
green regions in Fig. 5, which clearly shows that in the presence of a non-
zero curvature, the asymmetry of $u$ results in $l_{0}\neq 0$. A
straightforward but tedious calculation (see D for details) shows that

$l_{0}=-\frac{\epsilon^{2}}{4}\kappa+\mathcal{O}(\epsilon^{3}).$ (3.8)

By reinitializing the orientations of all $\bm{x}$ with $l(\bm{x})<l_{0}$ to
$\theta_{\rm L}$, and to $\theta_{\rm R}$ when $l(\bm{x})>l_{0}$, we have a
thresholding rule that moves the grain boundary by $\epsilon^{2}\kappa$ in one
time step $dt=\mathsf{t}\epsilon^{2}/4$, where $\mathsf{t}=1$ is a unit
conversion factor. Alternating between the $\eta$-update using the primal-dual
method, and the $\theta$-update using the thresholding rule, results in a
grain boundary motion by curvature with mobility equal to the inverse of the
grain boundary energy.111111In this case, the _reduced mobility_ (Salvador and
Esedoḡlu, 2019; Martine La Boissonière et al., 2019), which is defined as the
product of grain boundary energy and mobility, is equal to $1$ for all grain
boundaries. Although this is a severe restriction on the mobility, we
postulate that this can be overcome by modifying the thresholding rule (3.7),
and this will be addressed in a future work. The efficiency of the
thresholding rule described above rests on the computation of $l_{0}$ in
(3.7). In the next section, we use the _fast marching method_ to not only
compute $l_{0}$ in an $\mathcal{O}(N\log N)$ algorithm, but also generalize
the above strategy to an arbitrary polycrystal.

###  3.2 Thresholding dynamics via the fast marching method

The fast marching method (FMM), developed by Tsitsiklis (1995), is an
algorithm to evolve a surface with a spatially varying normal velocity. A
general description of FMM with a stand-alone example is given in E. Here, we
focus on using the fast marching method to implement the thresholding
algorithm described in Section 3.1 by solving for $l_{0}$ in (3.7).

(a) Two interior regions of a bicrystal.

(b) The evolution of $\partial I_{1}$ and $\partial I_{2}$ solved using the
fast marching method.

Figure 6: Shrinking of a circular grain simulated using the thresholding
method. 6(a)) Two interior regions (red and blue) $I_{1}$ and $I_{2}$ are
grown towards the grain boundary with a speed $1/(1-\eta^{*})^{2}$ using the
fast marching method. 6(b)) A closeup of a rectangular region around the grain
boundary, depicted in (6(a)), shows the contour lines of the fast marching
method, which describe the time it takes for $\partial I_{1}$ or $\partial
I_{2}$ to arrive at a grid point. Therefore, the original grain boundary,
shown as a dashed black line in (6(b)), moves to a new position (solid green
line) where the two grain interiors meet.

We begin with a description of our implementation of the thresholding scheme
for a bicrystal consisting of a circular grain, followed by its generalization
to a polycrystal. We first recall that the boundary conditions
$u(\pm\infty)=1$ used to arrive at (3.6)–(3.7) apply only in the limit
$\epsilon\to 0$ as noted in footnote 9. In practice, we choose a finite limit
$l_{\rm b}>0$, and modify (3.7) as

$\text{Find $l_{0}$ such that }\int^{l_{0}}_{-l_{\rm
b}}(1-u(l/\epsilon))^{2}\,dl=\int^{l_{\rm
b}}_{l_{0}}(1-u(l/\epsilon))^{2}\,dl.$ (3.9)

In D, we show that the error in $l_{0}$ due to the introduction of $l_{\rm b}$
exponentially decreases as $\epsilon\to 0$. In order to use FMM to compute
$l_{0}$, we interpret the integrand $(1-u(l/\epsilon))^{2}$ in (3.9) as an
inverse of the normal velocity of a surface
$\mathcal{S}_{l}:=\\{\bm{x}:l(\bm{x})=l\\}$ traveling towards the grain
boundary. Under this interpretation, the integrals in (3.9) are a measure of
the time it takes for two initial surfaces $\mathcal{S}_{-l_{\rm{b}}}$ and
$\mathcal{S}_{l_{\rm b}}$ on either side of the grain boundary, to meet at
$l=l_{0}$. We use the fast marching method to evolve the surfaces
$\mathcal{S}_{-l_{\rm{b}}}$ and $\mathcal{S}_{l_{\rm b}}$, and implement the
thresholding rule described in Section 3.1 by reassigning the orientation of
any point $\bm{x}$ in the region $\\{\bm{x}\in\Omega:|l(\bm{x})|<l_{\rm b}\\}$
to $\theta_{\rm L}$ if it first encounters the evolving surface
$\mathcal{S}_{-l_{\rm b}}$, and to $\theta_{\rm R}$ otherwise.

In practice, however, we do not have access to the signed distance function
$l(\bm{x})$ to identify the surfaces $\mathcal{S}_{-l_{\rm b}}$ and
$\mathcal{S}_{l_{\rm b}}$. Instead, we first identify the grain interiors
$I_{p}$ defined as

$I_{p}=\\{\bm{x}\in\Omega:\theta(\bm{x})=\theta_{p},\;\eta(\bm{x})>1-\xi\\},$
(3.10)

where $\xi>0$ is some fixed small value. Fig. 6(a) shows the grain interiors
$I_{1}$ and $I_{2}$ in a bicrystal, and Fig. 6(b) is a closeup of a
rectangular region, marked in yellow, around the grain boundary. The original
grain boundary is marked as a black dashed line in Fig. 6(b). By construction,
the two surfaces $\partial I_{1}$ and $\partial I_{2}$ are equidistant, up to
$\mathcal{O}(\epsilon)$, from the grain boundary, and serve as substitutes for
$\mathcal{S}_{-l_{\rm b}}$ and $\mathcal{S}_{l_{\rm b}}$. The grain interiors
are grown in the outward direction with a velocity $(1-u(l/\epsilon))^{-2}$
using the fast marching method, and the surface where they meet is the new
grain boundary, shown as a green dashed line in Fig. 6(b).

We will now generalize the above implementation to an arbitrary polycrystal
consisting of $\mathcal{N}$ grains, described using a piecewise constant
$\theta$ with values in $\\{\theta_{1},\dots,\theta_{\mathcal{N}}\\}$. Using
(3.10), we identify the $\mathcal{N}$ grain interiors, and define $I$ as their
union. Next, we grow the grain interiors in their outward unit normal
directions until every point (in the almost everywhere sense) in the domain is
in precisely one grain. We implement this by first collecting all the
boundaries of the interior regions in $\partial I=\partial
I_{1}\cup\dots\cup\partial I_{n}$, and _simultaneously_ evolving them in the
outward normal direction with a speed of $1/(1-\eta^{*}(x))^{2}$ using the
fast marching method.121212Note that the fast marching method is used to
evolve _all_ grain interiors in unison as opposed to evolving them
individually. As the grain interiors grow, a point $\bm{x}\in\Omega-I$ is
reinitialized to an orientation $\theta_{q}$ if it encounters $\partial
I_{q}\subset\partial I$. At the end of the fast marching method, all points in
$\Omega-I$ have been reinitialized resulting in an updated polycrystal at the
end of a time step. Fig. 7 shows the implementation of the thresholding rule
in a tricrystal. Dirichlet boundary conditions on $\theta$ are imposed by
including all $\bm{x}\in\partial\Omega$ in the grain interiors. On the other
hand, periodic boundary conditions are achieved by periodically reinitializing
$\theta$ for $\bm{x}\in\partial\Omega$ during the fast marching step.

Figure 7: Movement of a triple junction according to the thresholding
algorithm. The triple junction initially at $(x,y)=(0.5,0.5)$ (black filled
circle) moves to a new position (yellow filled circle) where the three grain
interiors, evolved using the fast marching method, meet at the same time.

Input : A polycrystal with $\mathcal{N}$ grains with orientations
$\theta_{1},\dots,\theta_{\mathcal{N}}$, grain boundary core energies
$\mathcal{J}(\llbracket\theta\rrbracket)$; parameters: $\epsilon$, $\xi$,
total time $T$, and tolerance $\mathsf{e}$.

Output : Time evolution of the polycrystal

Construct the core energy function $\mathcal{J}([\\![\theta]\\!])$ from grain
boundary energy data

Initialize $t=0$, and the orientation field $\theta(\bm{x},0)$

while  _ $t<T$_ do

Compute the discrete jump fields $\llbracket\theta\rrbracket(\bm{x},t)$ and
$\bar{\mathcal{J}}:=\mathcal{J}(\llbracket\theta\rrbracket(\bm{x},t))$ on
$\Omega$

Regularize the jump field: $\mathcal{J}^{\star}=G*\bar{\mathcal{J}}$, where
$G(\bm{x})=(1/2\pi\epsilon^{2})e^{-\frac{\bm{|}\bm{x}|^{2}}{2\epsilon^{4}}}$

// Solve for $\eta(\bm{x},t)$ using the primal-dual algorithm

Initialize $\eta$ and the dual field $\psi$:
$\eta_{0}(\bm{x})=0,\;\psi_{0}(\bm{x})=0$, and $n=0$

do

$n=n+1$

Calculate $\eta_{n}$ using $\psi_{n-1}$ (C.6)

Calculate $\psi_{n}$ using $\eta_{n}$ (C.7)

while  _ $\|\eta_{n+1}-\eta_{n}\|_{\infty}\leq\mathsf{e}$ _;

$\eta(\bm{x},t)=\eta_{n+1}$

// Threshold/update the orientation field

Identify interiors of grains:
$I_{p}=\\{\bm{x}\in\Omega:\theta(\bm{x},t)=\theta_{p},\,\bar{\mathcal{J}}(\bm{x})<\xi\\}$,
and set $I=\cup_{p=1}^{n}I_{p}$

Evolve $I$ with speed $1/(1-\eta(\bm{x},t))^{2}$ using the fast marching
method and update/threshold the orientations at each point $\bm{x}\in\Omega-I$

$t=t+0.25\epsilon^{2}$ (3.8)

end while

Algorithm 1 Thresholding algorithm for the new KWC model

The primal-dual and the fast marching methods are implemented on a regular
grid of resolution, say $\delta x$. From (3.8), we know that the resolution of
the grid should be large enough to resolve a grain boundary movement of
$\epsilon^{2}\kappa$ in each time step, i.e.

$\delta x\ll\epsilon^{2}\kappa,$ (3.11)

which is a common requirement of other thresholding methods (Esedoḡlu and
Otto, 2015; Merriman et al., 1992). If this condition is not satisfied, grain
boundaries would stagnate. Since grain boundary evolution results in an
overall decrease in curvature, (3.11) may cease to hold as the simulation
progresses. Therefore, we adaptively increase $\epsilon$ when a grain boundary
stagnates, and as a consequence, we obtain a time adaptive algorithm since
$dt\propto\epsilon^{2}$. On the other hand, an extremely small $\epsilon$ will
increase the computational cost of the thresholding method.

| Grid Size
---|---
$\xi$ | $1024\times 1024$ | $2048\times 2048$
0.15 | 9.77 % | 2.32 %
0.10 | 7.74 % | 1.24 %
0.05 | 3.39 % | 0.71 %
0.02 | 2.58 % | 0.07 %
Table 1: The effect of parameter $\xi$ on deviations from the expected motion
by curvature. We note that for a $2048\times 2048$ grid, $\xi=0.05$ is small
enough to achieve an error less than 1%.

Finally, we explore the effect of $\xi$, introduced in (3.10), on the extent
to which (3.8) is satisfied. Recall that $\xi$ was introduced in (3.10) to
identify grain interiors. In the case of a circular grain (see Fig. 6(a)),
(3.8) implies the rate of change of radius is given by

$\dot{R}(t)=-\frac{\epsilon^{2}}{4R(t)}.$ (3.12)

To test if the above equation is satisfied, we executed the thresholding
algorithm using the $\eta$-solution from the primal dual algorithm with
$\epsilon=0.01$, and measured $\dot{R}(t)$. Relative % errors in shrinking-
rate $\dot{R}$ at different values of $\xi$ are summarized in Tab. 1. It is
confirmed that for a sufficiently small grid, $\xi=0.05$ is small enough to
achieve an error less than 1%.

Algorithm 1 summarizes our approach. The core energy data $\mathcal{J}$ (e.g.,
Fig. 3) is computed separately using the procedure described in Section 2.2,
and used as an input to our method. The algorithm alternates between the
primal-dual and the fast marching methods resulting in motion by curvature. A
C++ template library that implements Algorithm 1 is available at
https://github.com/admal-research-group/GBthresholding.

Finally, we remark on the computation of $\mathcal{J}(\bm{x})$ on a grid.
Since $\llbracket\theta\rrbracket$, calculated at a grid point $ij$ in either
$x$\- or $y$-directions using centered-difference, is shared between two grid
points, a factor of $1/2$ appears in the following expression used to compute
the total jump:

$\llbracket\theta\rrbracket_{ij}=\frac{1}{2}\sqrt{(\theta_{i+1,j}-\theta_{i-1,j})^{2}+(\theta_{i,j-1}-\theta_{i,j+1})^{2}}.$
(3.13)

##  4 Numerical results

In this section, we present examples that explore various features of grain
boundary evolution predicted by our model.

(a)

(b)

Figure 8: A comparison of the numerical solution resulting from the
thresholding algorithm, implemented with $\epsilon=0.1$ on a $512\times 512$
grid, with the analytical solution. Plots of 8(a)) the steady-state solution
$\eta$, and 8(b)) grain boundary energy as a function of misorientation.

We begin with a simulation of a one-dimensional bicrystal $\Omega=[0,1]$ with
a grain boundary at $x=0.5$, and
$\mathcal{J}(\llbracket\theta\rrbracket)=\llbracket\theta\rrbracket$. The
purpose of this simulation is to ensure that the results of the primal dual
algorithm are consistent with the analytical model described in A. A Neumann
boundary condition $d\eta/dx=0$ is enforced at the two ends. In the absence of
a curvature, we expect the grain boundary to remain at $x=0.5$, and $\eta$
reach its steady state. The tolerance $\mathsf{e}$ of the primal dual
algorithm (C.8) is set to $10^{-6}$. Fig. 8(a) confirms the agreement between
$\eta$ obtained from the primal dual algorithm and the analytical form given
in (A.8). Furthermore, Fig. 8(b) shows that the grain boundary energies
predicted by the primal-dual algorithm for various misorientation angles are
in agreement with the analytical result in (A.10).

###  4.1 Equilibrium of a triple junction

(a) Initial Condition

(b) Gaussian Kernel Method

(c) Current method

Figure 9: A comparison of the evolutions of a tricrystal under periodic
boundary conditions obtained using the Gaussian kernel method and the
generalized KWC model with $\epsilon=0.01$, implemented using our method. The
dihedral angles predicted by the Gaussian kernel method and our method are
$(93^{\circ},133.5^{\circ},133.5^{\circ})$ and
$(106^{\circ},127^{\circ},127^{\circ})$ respectively, while the theoretical
values are $(90.89^{\circ},134.56^{\circ},134.56^{\circ})$. In Fig. 10, we
demonstrate that the error in the dihedral angles predicted by the generalized
KWC model decreases as $\epsilon\to 0$.

(a) Initial Condition

(b) $\epsilon$=0.01

(c) $\epsilon$=0.006

(d) $\epsilon$=0.003

Figure 10: 10(a)) The orientation distribution in an initial tricrystal under
periodic boundary conditions with a triple junction at $(x,y)=(0.25.0.75)$,
and dihedral angles
$(\Theta_{1},\Theta_{2},\Theta_{3})=(90^{\circ},180^{\circ},90^{\circ})$. The
polycrystal is evolved using the thresholding algorithm with $\epsilon=0.01$,
$0.006$ and $0.003$. 10(b))-10(d)) Closeups of an evolving triple junction
(red box) clearly show that the dihedral angles converge to
$(\Theta_{1},\Theta_{2},\Theta_{3})=(134.56^{\circ},90.89^{\circ},134.56^{\circ})$
predicted by the Herring angle condition (4.1), as $\epsilon$ converges to
zero.

A triple junction is a line where three grains meet, and it is represented as
a point in two dimensions. The equilibrium of a triple junction is guaranteed
if it satisfies the Herring relation (Herring, 1951) given by

$\frac{\gamma^{12}}{\sin\Theta_{3}}=\frac{\gamma^{23}}{\sin\Theta_{1}}=\frac{\gamma^{31}}{\sin\Theta_{2}},$
(4.1)

where $\Theta_{i}$ is the dihedral angle of grain $i$, and $\gamma^{ij}$ is
the energy density of the grain boundary shared by grains $i$ and $j$. While
the Herring relation is derived in the sharp-interface framework, it is also
seen to hold for a triple junction governed by the original KWC model through
(2.4). This is not surprising since the KWC model converges to the Mullins
model in the sharp-interface limit _and_ the evolution in (2.4) has a
variational structure in the form of a gradient descent of the functional in
(2.1). On the other hand, it is not clear if our approach to evolve the
generalized KWC model arises from a variational formulation. Therefore, it is
necessary to examine the Herring relation using our thresholding algorithm. We
will now demonstrate that the Herring relation indeed holds provided the
parameter $\epsilon$ is chosen appropriately.

We study the evolution of a triple junction in a tricrystal with orientations
$\theta_{1}=0$, $\theta_{2}=\pi/6$, and $\theta_{3}=\pi/3$ in
$\Omega=[0,1]\times[0,1]$. Using the Read–Shockley core energy
$\mathcal{J}=\llbracket\theta\rrbracket$, we note from Fig. 8(b) that the
energy density of the three grain boundaries are $\gamma^{12}=0.62$,
$\gamma^{23}=0.62$, and $\gamma^{13}=0.87$. From the Herring relation in
(4.1), it follows that the steady state dihedral angles are
$\Theta_{1}=134.56^{\circ}$, $\Theta_{2}=90.89^{\circ}$, and
$\Theta_{3}=134.56^{\circ}$. In order to examine the Herring relation, we
consider a tricrystal under periodic boundary conditions, with an initial
orientation distribution given by

$\theta(\bm{x},t=0)=\begin{cases}\theta_{2}&\text{if }x_{2}\leq 0.25\text{ or
}x_{2}>0.75,\\\ \theta_{3}&\text{if }0.25<x_{2}\leq 0.75\text{ and }0.25\leq
x_{1}<0.75,\\\ \theta_{1}&\text{if }0.25<x_{2}\leq 0.75\text{ and
}x_{1}>0.25\text{ or }x_{1}>0.75,\end{cases}$ (4.2)

resulting in four triple junctions at
$(x_{1},x_{2})=(0.25,0.75),(0.75,0.75),(0.75,0.25)$, and $(0.25,0.25)$. The
initial dihedral angles of the triple junctions are $90^{\circ}$,
$180^{\circ}$, and $90^{\circ}$. Fig. 9(a) shows a plot of the initial
orientation distribution of the tricrystal.

We begin by comparing the evolution of a triple junction predicted by the KWC
model implemented using our thresholding scheme with that obtained using the
Gaussian kernel method (Esedoḡlu and Otto, 2015). The grain boundary energies
$(\gamma^{12},\gamma^{23},\gamma^{13})=(0.62,0.62,87)$, pre-computed using the
KWC model, are inputs to the Gaussian kernel method, and the respective
mobilities are set to the inverse of the grain boundary energies. The
parameter $\epsilon$ of the KWC model is taken as $0.01$. Both schemes are
simulated on a $1024\times 1024$ grid. As shown in Fig. 9, the evolution
dynamics of both schemes are qualitatively similar. As expected, the triple
junction adjusts at a faster time scale to satisfy the Herring angle condition
compared to the curvature-driven motion of grain boundaries (Esedoḡlu et al.,
2010). The motion of triple junctions induces a curvature in the grain
boundaries, which drives the shrinking of the embedded grains (blue and
yellow), while maintaining constant dihedral angles. The dihedral angles
predicted by the Gaussian kernel method are
$(93^{\circ},133.5^{\circ},133.5^{\circ})$, while the generalized KWC model
with $\epsilon=0.01$ yields $(106^{\circ},127^{\circ},127^{\circ})$.

To further investigate the dependence of the triple junction angles on
$\epsilon$, we implement the thresholding algorithm with $\epsilon=0.01$,
$0.006$ and $0.003$ on a $3000\times 3000$ grid. As shown in Fig. 10(b) to
Fig. 10(d), as $\epsilon$ decreases, the stabilized triple junction angles
converge to those predicted by the Herring relation. This test suggests that
the Herring relation is satisfied in the limit $\epsilon\to 0$.

###  4.2 The Von Neumann-Mullins Theory of Grain Growth

(a) Initial condition

(b) $t=5\times 10^{-3}$

Figure 11: The ideal grain growth simulate by the current scheme. The constant
core energy function $\mathcal{J}(\llbracket\theta\rrbracket)=0.5$ results an
isotropic grain boundary energy.  Figure 12:  The mean and standard deviation
of area change-rate for $n$-sided grain during the ideal grain growth shown in
Fig. 11. The black dashed line $\dot{A}(n)=(\pi/3)n-2\pi$ indicates the von
Neumann-Mullins equation (4.3) with the unit reduced mobility $m\gamma=1$.

In this section, we validate our thresholding scheme by testing the von
Neumann–Mullins relation for a polycrystal with uniform grain boundary
energies ($\gamma$) and mobilities ($m$). In two dimensions, the von Neumann-
Mullins law (von Neumann, 1952; Mullins, 1956) states

$\frac{dA}{dt}=\frac{\pi}{3}m\gamma(n-6),$ (4.3)

where $A$ is the area of a grain with $n$ sides, and $\gamma$ . In other
words, grains with more than six sides grow, while those with less than six
sides will shrink.

To test the relation given in (4.3), we select the core energy function as a
constant ($\mathcal{J}(\llbracket\theta\rrbracket)=0.5$), and consider an
initial polycrystal consisting of $\mathcal{N}=50$ grains. The initial
configuration, shown in Fig. 11(a), is generated using a Voronoi tessellation
of uniformly distributed random points. We implemented Algorithm 1 with
parameters $\epsilon=0.01$, $\mathsf{e}=10^{-6}$, and $\xi=0.05$, on a
$1024\times 1024$ grid.

As described in Section 4.1, the evolution of the polycrystal begins with the
motion of triple junctions to attain the dihedral angles
$(120^{\circ},120^{\circ},120^{\circ})$ predicted by (4.1) for constant
$\gamma$. Subsequently, grain boundary motion by curvature follows. A snapshot
of an evolving grain microstructure at $t=5\times 10^{-3}$, simulated using
our thresholding scheme, is shown in Fig. 11(b). In Fig. 12, we plot the mean
rate of change of area of grains — along with standard deviation — measured
during the time interval $[2.5\times 10^{-3},5.0\times 10^{-3}]$, as a
function of the number of sides. Noting that the plot in red is close to the
theoretically predicted black dashed line, we confirm that our thresholding
scheme accurately predicts the von Neumann-Mullins law.

###  4.3 Comparison with the finite element implementation of the KWC model

(a)

(b)

(c)

Figure 13: An initial polycrystal, shown in (13(a)), is evolved using the
thresholding and the finite element methods resulting in polycrystals shown in
(13(b)) and (13(c)) respectively. The two methods are consistent in predicting
the growth (e.g., \scriptsize{1}⃝, \scriptsize{2}⃝) and shrinkage (e.g.,
\scriptsize{3}⃝, \scriptsize{4}⃝) in various grains. The differences in the
evolution is attributed to the mobility function introduced in (4.4) to
prevent grain rotation.

In this section, we compare the evolutions of a polycrystal resulting form our
thresholding method and a finite element implementation of the KWC model,
which we will refer to as _FE-KWC_. An initial polycrystal consisting of
$\mathcal{N}=50$ grains, as shown in Fig. 13(a), is generated using a Voronoi
tessellation of uniformly distributed random points. The orientations of the
grains are randomly chosen from the interval $[0,\pi/2]$. The core energy is
chosen to be of the Read–Shockley-type, i.e.
$\mathcal{J}=\llbracket\theta\rrbracket$.

The thresholding algorithm is implemented on a $1024\times 1024$ grid, with
parameters $\epsilon=0.01$, $\mathsf{e}=10^{-6}$, and $\xi=0.05$. A snapshot
of an evolving grain microstructure at $t=5\times 10^{-3}$, simulated using
our thresholding scheme, is shown in Fig. 13(b).

(a)

(b)

(c)

Figure 14: An initial polycrystal, shown in (14(a)), is evolved using the
finite element method. (14(b)) and (14(c)) show the resulting polycrystals
with regularization parameters $\rho_{0}=2\times 10^{-4}$ and
$\rho_{0}=2\times 10^{-3}$ respectively. When $\rho_{0}$ is not sufficiently
small, grains with small misorientation (e.g., \scriptsize{1}⃝,
\scriptsize{2}⃝) blend out and grain boundaries easily become rounded.
However, the decrease in $\rho_{0}$ for simulating sharp interfaces, comes
with significant computational cost contrasts to the suggested scheme.

We note that FE-KWC, using continuous Lagrange finite elements, cannot be
carried out on our model since the solution for $\theta$ is discontinuous.
Therefore, we proceed with a finite element implementation of the regularized
KWC model given in (2.4). We use second-order quadrilateral Lagrange finite
elements to interpolate the order parameters. Since the regularized model
allows grain rotation, we inhibit rotation using the following
$\eta$-dependent mobility for $\theta$

$b_{\theta}^{-1}(\eta)=10^{-5}\epsilon+\left(1-\eta^{3}(10-15\eta+6\eta^{2})\right)(1-10^{-5})\epsilon,$
(4.4)

as suggested by Dorr et al. (2010). On the other hand,
$(b_{\eta})^{-1}=\epsilon$ is chosen to be constant. To address the
singularity due to the $|\nabla\theta|$ term in (2.4b), we use the
approximation

$g(\eta)|\nabla\theta|\approx g(\eta)\sqrt{\rho_{0}+|\nabla\theta|^{2}},$
(4.5)

where $\rho_{0}=2\times 10^{-3}$ is a constant. The manifestation of
$\rho_{0}$ on the solution will be discussed below. FE-KWC is performed on the
open-source computing platform Fenics Alnaes et al. (2015). We take an
implicit time step with $dt=0.012$. In order to compare the numerical
efficiency, we ensure that the number of degrees of freedom is the same in the
thresholding and the finite element simulations. The grain microstructure at
$t=1.71$, simulated using FE-KWC, is shown in Fig. 13(c).

Comparing Figs. 13(b)–13(c), we note that both the methods are consistent in
predicting growth (see \footnotesize{1}⃝, \footnotesize{2}⃝) and shrinkage
(see \footnotesize{3}⃝, \footnotesize{4}⃝) in various grains. It is observed
that grain boundaries become rounded in the finite element simulation, because
of the diffusive nature of orientation field. In addition, disparities are
more clear for small misorientation grain boundaries, e.g. \footnotesize{5}⃝,
which are highly diffused. This is a manifestation of $\rho_{0}$, which
results in a non-zero gradient in $\theta$ in the grain interiors. In Fig. 14,
we compare two finite element simulations with $\rho_{0}=2\times 10^{-3}$ and
$2\times 10^{-4}$, which shows that for a smaller $\rho_{0}$, the grain
boundaries retain their characteristic width.131313Recall that the
characteristic width of a grain boundary in the regularized KWC model is a
function of $\epsilon$. Thus, to simulate sharp grain interfaces comparable to
the our scheme, a small enough $\rho_{0}$ is required for FE-KWC. However, we
note that the decrease in $\rho_{0}$ increases the stiffness of the equations,
which significantly affects the computational time as discussed below.

Figure 15: A comparison of the complexity of the thresholding and the finite
element methods. The dashed lines represent reference slopes in the log-log
scale plot. Slope 1 and 2 represent $\mathcal{O}(N)$ and $\mathcal{O}(\log N)$
respectively. While both methods have a complexity of at most
$\mathcal{O}(N\log{N})$, the cost of the finite element method depends on the
choice of $\rho$.

Computational time study clearly highlights the advantage of the our
thresholding scheme. Performance tasks are executed on a single 1.6 GHz core
with 8 GB RAM, and we measured the wall-clock time to complete one-full time
step for the two methods. For our method, this includes solving for $\eta$
using the primal dual algorithm, and executing the fast marching based
thresholding algorithm to update $\theta$. In Fig. 15, we plot the dependence
of the wall-clock time, as a function of the number of degrees of freedom $N$.
The computational complexity of the current scheme is $\mathcal{O}(N\log{N})$,
with a dominant contribution from FFT used in the primal dual algorithm to
solve (C.7). On the other hand, the asymptotic computational cost of FE-KWC is
estimated to be in between $\mathcal{O}(N)$ and $\mathcal{O}(N^{2})$ as shown
in Fig. 15. The computational bottleneck of FE-KWC is in solving — using a
GMRES iterative solver (Saad and Schultz, 1986) — a linear system of equations
formed by an $N\times N$-sized sparse matrix. Though the asymptotic costs of
the two schemes are similar in terms of $N$, we note that the computational
cost of FE-KWC also depends on the choice of the regularization parameter
$\rho_{0}$, which increases the stiffness of the equations in the limit
$\rho\to 0$. Therefore, as demonstrated in Fig. 15, the current scheme can be
orders of magnitude faster than FE-KWC. Both, FE-KWC and the implementation of
our method, can be well-parallelized using the current generation of graphics
cards, which have the power, programmability and precision to implement FFT
and iterative matrix solvers (Li and Saad, 2013; Govindaraju et al., 2008)
respectively.

###  4.4 Grain growth in an fcc copper polycrystal

Figure 16:  Grain boundary energies used for the polycrystal simulation in
Section 4.4. Using a core energy $\mathcal{J}([\\![\theta]\\!])$ designed in
Section 2.2, we obtain a crystal symmetry-invariant KWC model with energy that
matches the covariance model. In order to compare the original and the new KWC
models, we scale the function $g$ of the original KWC model in (2.1) to
$g=-0.93\ln(1-\eta)$ such that the averages of the grain boundary energies
(with respect to misorientation) are identical in the two models. In other
words, the areas under the above plots are equal.

(a)

(b)

(c)

Figure 17: (17(a)) A polycrystal with $\mathcal{N}=50$ grains, and an initial
orientation distribution. (17(b)) and (17(c)) show evolved polycrystals using
the new and the original KWC models respectively. Grains $1$ and $2$ show
opposite growth/shrinkage trends in the two models due to the deviation of the
grain boundary energy from the Read–Shockley-type in the new formulation. The
blue and red colors represent the maximum and minimum orientation angles of
$0^{\circ}$ and $70.6^{\circ}$ respectively.

(a)

(b)

Figure 18: Initial distributions of grain boundary energies for the (18(a))
generalized, and (18(b)) the original KWC models. The grain boundary energy in
the original KWC model is of the Read–Shockley-type. On the other hand, the
grain boundary energy in the generalized KWC model reflects the crystal
symmetry of copper.

In this section, we examine grain growth in a two-dimensional fcc copper
polycrystal with $[110]$-type grain boundaries simulated using the generalized
KWC model,141414The $[110]$ direction of each grain is is aligned with the
$z$-axis (out of the plane). with crystal symmetry-invariant grain boundary
energy. We compare the results with the predictions of the original KWC model.

A two-dimensional polycrystal consisting of $\mathcal{N}=50$ grains, with
orientations in the range $[0,70.6^{\circ}]$ is generated using a Voronoi
tessellation of random points. Fig. 17 shows the initial orientation
distribution in the polycrystal. We assume that the grain boundary energy
density is independent of inclination. We use the core energy
$\mathcal{J}(\llbracket\theta\rrbracket)$ constructed in Section 2.2 (see Fig.
3). In order to compare the generalized model to the original KWC model, we
scale the function $g$ of the original KWC model in (2.1) to
$g=-0.93\ln(1-\eta)$ such that the mean of the grain boundary energies as
functions of misorientation in the range $[0,70.6^{\circ}]$ are identical for
the two models. Fig. 16 shows a comparison of the grain boundary energy
densities of the two models.

The orientation distributions of the polycrystal at the end of 200 time steps
for the generalized and the original KWC models are shown in Figs. 17(b)–17(c)
respectively. Comparing the resulting polycrystals with the initial
polycrystal in Fig. 17(a), we note that the generalized KWC model predicts a
growth for red grains while the original model results in their shrinkage.
This can be attributed to the difference in the grain boundary energies of the
two models, as shown in Fig. 18. For example, the grain boundary
\footnotesize{1}⃝, which has a misorientation of $\approx 70.6^{\circ}$ has a
relatively smaller energy in the generalized model due to crystal symmetry.

On the other hand, we note an opposite trend for light blue grains for which
the generalized model predicts shrinkage while the original model results in a
growth. This is a result of relatively larger energy of grain boundary
\footnotesize{2}⃝ in the generalized model compared to the original model. The
above observations suggest that the generalized model can result in the growth
of certain grains with large misorientation, highlighting the importance of
crystallography in grain growth.

##  5 Conclusion and future work

In this work, we generalized the two-dimensional KWC model for grain
boundaries such that it can incorporate arbitrary misorientation-dependent
grain boundary energies that respect the bicrystallography of grain
boundaries. In addition, we address the computational challenge of solving a
singular diffusive equation of the KWC model by developing an
$\mathcal{O}(N\log N)$ thresholding algorithm. Below, we summarize the
construction of our model, its implementation, and research directions for
future work.

First, we eliminate the $|\nabla\theta|^{2}$ term in the original KWC model,
which is responsible for regularizing the orientation order parameter, and
rendering a non-zero mobility to the grain boundaries. The lack of a
regularizing term results in a discontinuous orientation order parameter with
jump across the grain boundaries, and grain boundaries with no mobility. In
the presence of a piecewise-constant orientation order parameter, the modified
KWC functional can be separated into a bulk contribution that depends on
$\eta$, and a surface contribution, called the core energy, that depends
linearly on the jump $\llbracket\theta\rrbracket$ across the grain boundaries.
Next, we show that by generalizing the core energy from a linear function to
an arbitrary function $\mathcal{J}$ of $\llbracket\theta\rrbracket$, the model
can incorporate arbitrary dependence of grain boundary energies on
misorientation angles.

Since the absence of the regularizing $|\nabla\theta|^{2}$ term renders the
grain boundaries immobile, we design an $\mathcal{O}(N\log N)$ thresholding
algorithm to evolve grain boundaries by curvature, where $N$ is the number of
grid points. The algorithm, which employs a primal-dual and the fast marching
methods, is shown to be an order of magnitude faster than the finite element
implementation of the original KWC model. We validate our implementation by
predicting the Herring angle relation, and simulate a two-dimensional
polycrystal consisting of $[110]$ tilt grain boundaries. The computational
efficiency and flexibility of our approach opens the door to a number of
exiting directions for future work.

  * 1.

The present framework will enable us to carry out a statistical study of large
scale simulations of various ensembles of polycrystals to characterize
abnormal grain growth in terms of the grain boundary energy landscape and
crystal symmetry.

  * 2.

While arbitrary grain boundary energies can be incorporated into our model,
its implementation is restricted to grain boundary mobility equal to the
inverse of the energy. An extension of our algorithm to include mobilities
independently will be explored in a future work.

  * 3.

The present algorithm does not allow grain rotation, which is another
important phenomenon during recrystallization of polycrystalline
materials.151515We note that grain rotation may sometimes play an important
role during the transition from recovery to continuous dynamic
recrystallization. Dislocations agglomerate and form cell walls/subgrains at
the end of the recovery stage. In a phenomenon, commonly referred to as
_subgrain rotation recrystallization_ , few subgrains — aided by bulk
dislocations – increase their misorientation and transform to grains/nuclei
which grow (Li, 1962). From this perspective, grain rotation plays an
important role during the nucleation of recrystallized grains. We plan to
augment the current scheme with a step that models grain rotation.

  * 4.

A recent work by Admal et al. (2019) extended the two-dimensional KWC model to
a three-dimensional fully anisotropic (both misorientation and inclination
dependent) model, wherein the dependence of grain boundary energy on the
misorientation angle was restrictive to a Read–Shockley-type. Due to the high
computational cost of the finite element method, the implementation of the
three-dimensional model was restricted to simple bicrystals. It is envisaged
that the efficiency of our thresholding algorithm will enable us to explore
large three-dimensional polycrystals with fully anisotropic grain boundary
energy.

  * 5.

Finally, we recall from the introduction that surface tension is not the only
dominant driving force on a grain boundary due to grain boundary plasticity.
Adapting our thresholding algorithm into existing unified frameworks (Admal et
al., 2018), wherein grain microstructure and deformation evolve
contemporaneously, will enable us to quantify the role of grain boundary
plasticity, and study phenomena such as dynamic recrystallization,
superplasticity and severe plastic deformation (Thomas et al., 2017; Wei et
al., 2020; Runnels and Agrawal, 2020).

## Data availability

A C++ template library that implements Algorithm 1 is available at
https://github.com/admal-research-group/GBthresholding.

## CRediT author statement

Jaekwang Kim: Formal analysis, Investigation, Software, Validation, Writing-
Original draft, Visualization. Matt Jacobs: Conceptualization, Methodology,
Software, Formal analysis. Stanley Osher: Conceptualization. Nikhil Admal:
Conceptualization, Methodology, Writing-Review and Editing, Resources, Data
Curation, Supervision, Project Management.

##  Appendix A Results on the 1D KWC model

In this section, we collect results on the one-dimensional KWC model which
describes an infinite bicrystal with a grain boundary at the origin. In
particular, we present the derivation of the steady-state analytical solution
under Dirichlet boundary conditions, and the resulting grain boundary energy
as a function of misorientation.

Consider the following KWC energy functional without the $|\nabla\theta|^{2}$
regularizing term

$\mathcal{W}[\eta,\theta]=\int_{-\infty}^{\infty}\left[\frac{\epsilon}{2}|\nabla\eta|^{2}+\frac{(1-\eta)^{2}}{2\epsilon}+g(\eta)|\nabla\theta|\right]\,dV.$
(A.1)

The Euler–Lagrange equation associated with the above functional is

$\epsilon\triangle\eta-\frac{\eta-1}{\epsilon}-g,_{\eta}|\nabla\theta|=0,$
(A.2)

where $g,_{\eta}$ is used to denote $\partial g/\partial\eta$. In what
follows, we derive a steady-state solution of (A.2) under Dirichlet boundary
conditions

$\eta(\pm\infty)=1,\quad\theta(\infty)=-\theta(-\infty)=\theta/2.$ (A.3)

We begin with the ansatz that $\theta(x)$ is a step function satisfying (A.3)
with a discontinuity at the origin. Multiplying (A.2) by $\eta^{\prime}$, and
integrating with respect to $x$ in a region away from the origin, we obtain

$\frac{\epsilon}{2}\eta,^{2}_{x}-\frac{(1-\eta)^{2}}{2\epsilon}=0,\implies\eta,_{x}=\pm\frac{(1-\eta)}{\epsilon},$
(A.4)

On the other hand, multiplying (A.2) with $\eta^{\prime}$, and integrating
over an arbitrarily small neighborhood of $0$ results in the jump condition

$\epsilon\llbracket\eta,_{x}\rrbracket=g,_{\eta}(\bar{\eta})\llbracket\theta\rrbracket,$
(A.5)

where $\bar{\eta}:=\eta(0)$ is the value of $\eta$ at the grain boundary. From
(A.4) and (A.5), it follows that

$\epsilon\eta,_{x}=\begin{cases}1-\eta&\text{ if }x>0,\\\ -(1-\eta)&\text{
otherwise},\\\ \end{cases}$ (A.6)

and

$g,_{\eta}(\bar{\eta})\llbracket\theta\rrbracket=2(1-\eta),$ (A.7)

which relates $\bar{\eta}$ to $\llbracket\theta\rrbracket$. The analytical
solution for $\eta$ can be obtained by integrating (A.4). With our choice of
$g=-\ln(1-\eta)$, the result can be explicitly written as a function of
misorientation $\llbracket\theta\rrbracket$:

$\int_{\bar{\eta}}^{\eta}\frac{\epsilon}{1-\eta}d\eta=x,\implies\eta(x)=1-\sqrt{\frac{\llbracket\theta\rrbracket}{2}}\exp\left(-\frac{|x|}{\epsilon}\right).$
(A.8)

The grain boundary energy $\gamma$ as a function of misorientation is
calculated by evaluating $\mathcal{W}[\eta,\theta]$ using the steady state
solution for $\eta$ derived above. From (A.4), we have

$\displaystyle\gamma(\llbracket\theta\rrbracket)=\mathcal{W}[\eta,\theta]$
$\displaystyle=\int_{-\infty}^{\infty}\left[\frac{\epsilon}{2}\eta_{,x}^{2}+\frac{(1-\eta)^{2}}{2\epsilon}\right]\,dx+g(\bar{\eta})\llbracket\theta\rrbracket$
$\displaystyle=2\int_{0}^{\infty}\frac{(1-\eta)^{2}}{\epsilon}\,dx+g(\bar{\eta})\llbracket\theta\rrbracket$
$\displaystyle=2\int_{\bar{\eta}}^{1}(1-\eta)\,d\eta+g(\bar{\eta})\llbracket\theta\rrbracket=(1-\bar{\eta})^{2}+g(\bar{\eta})\llbracket\theta\rrbracket.$
(A.9)

Note that the grain boundary energy $\gamma$ and $\bar{\eta}$ are independent
of $\epsilon$, which reinforces that the model converges to its sharp
interface as $\epsilon\to 0$ while the energy remains unchanged. Again
specializing the analytical expression $\gamma$ with the choice of the
logarithmic $g$, we obtain

$\displaystyle\gamma(\llbracket\theta\rrbracket)$
$\displaystyle=(1-\bar{\eta})^{2}-\llbracket\theta\rrbracket\ln(1-\bar{\eta})$
$\displaystyle=\frac{\llbracket\theta\rrbracket}{2}-\llbracket\theta\rrbracket\ln\left(\sqrt{\frac{\llbracket\theta\rrbracket}{2}}\right).$
(A.10)

##  Appendix B The covariance model of grain boundary energy

The covariance model for grain boundary energy, developed by Runnels et al.
(2016a, b), estimates grain boundary energy using the covariance of atomic
densities of the two lattices adjoining a grain boundary.

In the covariance model, a lattice density measure $\bar{\rho}$ for a given
lattice161616A lattice $\mathcal{L}$ is defined using three lattice vectors
$\bm{l}_{1}$, $\bm{l}_{2}$, and $\bm{l}_{3}$ as
$\displaystyle\mathcal{L}=\\{n_{1}\bm{l}_{1}+n_{2}\bm{l}_{2}+n_{3}\bm{l}_{3}\>n_{i}\in\mathbb{Z}\\}.$
(B.1) $\mathcal{L}$, defined as an infinite sum of Dirac measures with support
at the lattice points points of $\mathcal{L}$:

$\bar{\rho}(\bm{x})=\sum_{\bm{d}\in\mathcal{L}}\delta(\bm{x}-\bm{d}).$ (B.2)

A lattice density field $\rho$ is introduced as the convolution of
$\bar{\rho}$ with a thermalization function $\xi$, i.e.

$\rho(\bm{x})=\overline{\rho}(\bm{x})*\xi(\bm{x}),$ (B.3)

where

$\xi(\bm{x})=\frac{1}{\sigma^{3}\pi^{3/2}}e^{-\|\bm{x}\|^{2}/\sigma^{2}},$
(B.4)

with $\sigma^{2}$ as the dimensionless temperature. The planar covariance of
two thermalized lattices $\mathcal{L}_{A}$ and $\mathcal{L}_{B}$ with their
respective density fields $\rho_{A}$ and $\rho_{B}$, measured on
$\mathbb{R}^{2}$, is defined as

$c[\rho_{A},\rho_{B}]=\int_{\bm{y}\in\mathbb{R}^{2}}\rho_{A}(P^{T}\bm{y})\rho_{B}(P^{T}\bm{y})\lambda(\bm{y})\;dA,$
(B.5)

where $\lambda(\bm{x})$ is an appropriately chosen window function (see (B.8),
$P:\mathbb{R}^{3}\to\mathbb{R}^{2}$ is the projection

$P=\left(\begin{array}[]{@{}ccc@{}}1&0&0\\\ 0&1&0\\\ \end{array}\right)$ (B.6)

on to the plane $\mathbb{R}^{2}$. Expressing the functions $\rho_{A}$ and
$\rho_{B}$ in Fourier series, the integral in (B.5) simplifies as

$c[\rho_{A},\rho_{B}]=\frac{1}{\hat{\lambda}(\bm{0})}\sum_{\bm{k}_{A}\in\mathcal{L}^{\prime}_{A}}\sum_{\bm{k}_{B}\in\mathcal{L}^{\prime}_{B}}\widehat{\rho}(\bm{k}_{A})\widehat{\rho}^{\;*}(\bm{k}_{B})\hat{\lambda}\left(P(\bm{k}_{B}-\bm{k}_{A})\right),$
(B.7)

where $\bm{k}_{A}$ and $\bm{k}_{B}$ are lattice vectors of the dual lattices
$\mathcal{L}^{\prime}_{A}$ and $\mathcal{L}^{\prime}_{B}$, and the window
function is defined in terms of its Fourier transform as

$\hat{\lambda}(\bm{k})=e^{-\|\bm{k}\|/\omega},$ (B.8)

with an adjustable parameter $\omega$. The grain boundary energy in the
covariance model is defined as

$\gamma^{\mathrm{cov}}=E_{0}\left(1-\frac{c[\rho^{\mathcal{A}},\rho^{\mathcal{B}}]}{c_{\mathrm{gs}}}\right),$
(B.9)

where $c_{\mathrm{gs}}$ is the ground state covariance defined as the
supremum, over all planes, of $c[\rho_{A},\rho_{A}]$. For example, in fcc,
$c_{\mathrm{gs}}$ corresponds to covariance measured with respect to the
$[111]$ plane. Finally, we note that the covariance model has three adjustable
parameters $\\{E_{0},\sigma,\omega\\}$ that can be used to fit
$\gamma^{\rm{cov}}$ to data from experiments or molecular dynamics
simulations. It is known that while (B.9) is a good indicator of grain
boundary energy, it over-predicts the energy for low angle grain boundaries as
the above model does not account for facet formation. Runnels et al. (2016a,
b) have shown that a further relaxation of the grain boundary energy, which
signifies the formation of facets, yields necessary corrections to the energy
predicted by the model.

Figure 19: A plot of the normalized grain boundary energy versus the
misorientation angle predicted by the covariance model for a $[110]$
symmetric-tilt grain boundary in fcc copper, computed using the relaxation
algorithm of Runnels et al. (2016a, b). For comparison, grain boundary
energies obtained from experiment (Miura et al., 1994), and MD simulations
(Wolf, 1990) are shown in blue and square points respectively.

Fig. 19 shows a plot of a relaxed $\gamma^{\rm{cov}}$ computed for $[110]$
symmetric-tilt grain boundaries in fcc copper using
$E_{0}=1.45\,\mathrm{J}/\mathrm{m}^{-2}$, $\omega=0.5$, and
$\sigma/a_{0}=0.175$, where $a_{0}=3.597$ is the lattice constant of copper.
From Fig. 19, it is clear that the grain boundary energy predicted by the
covariance model is in good overall agreement with data from molecular
dynamics simulations (Wolf, 1990; Miura et al., 1994).

##  Appendix C The primal-dual method

Primal-dual methods is a part of a class of first-order algorithms171717 An
algorithm that only requires the calculation of the gradient of a functional.
that have a long history in the context of optimization problems (Powell,
1978; Kuhn, 1955; Komodakis and Pesquet, 2015). As the name suggests, primal-
dual methods proceed by concurrently solving a primal problem and a dual
problem. The main benefit of primal-dual splitting is that it replaces an
original hard problem with a set of two easy sub-problems (primal- and dual-).
Because of this advantage, the method has been widely used in diverse fields
including compressed sensing, image processing, signal processing, and machine
learning (Donoho, 2006; Chambolle and Pock, 2011; Shalev-Shwartz and Singer,
2007; Combettes and Pesquet, 2011).

The motivation to use a primal-dual algorithm to solve the minimization
problem in (3.1) for $\eta$ arises due to the presence of a highly nonlinear
term $g(\eta)\mathcal{J}(\llbracket\theta\rrbracket)$ along with
$|\nabla\eta|^{2}$. Therefore, we adopt a primal-dual method by introducing an
auxiliary dual variable which enables us to cast (3.1) as an equivalent
optimization problem. The choice of the dual variable is based on the
observation that

$\displaystyle\frac{\epsilon}{2}\|\nabla\eta\|_{L^{2}(\Omega)}^{2}$
$\displaystyle=\epsilon\|\nabla\eta\|_{L^{2}(\Omega)}^{2}-\frac{\epsilon}{2}\|\nabla\eta\|_{L^{2}(\Omega)}^{2}$
$\displaystyle=-\epsilon\int_{\Omega}\eta\triangle\eta\,dV-\int_{\Omega}\frac{\epsilon}{2}\nabla\eta\cdot\nabla\eta\,dV,$
(C.1)

where we have used the divergence theorem, and the Neumann boundary condition
$\nabla\eta\cdot\bm{n}=\bm{0}$. Introducing an auxiliary variable $\psi$, and
identifying it with $-\epsilon\triangle\eta$, we have

$\displaystyle\frac{\epsilon}{2}\|\nabla\eta\|^{2}$
$\displaystyle=\sup\limits_{\psi\in(\dot{H}^{1}(\Omega))^{*}}\left[\int_{\Omega}\eta(x)\psi(x)\,dV-\frac{1}{2\epsilon}\|\Delta^{-1}\nabla\psi\|_{L^{2}(\Omega)}^{2}\right]$
$\displaystyle=\sup\limits_{\psi\in(\dot{H}^{1}(\Omega))^{*}}\left[\int_{\Omega}\eta(x)\psi(x)\,dV-\frac{1}{2\epsilon}\|\psi\|_{(\dot{H}^{1}(\Omega))^{*}}^{2}\right],$
(C.2)

where $\dot{H}^{1}(\Omega)$ denotes the set of all functions in
$H^{1}(\Omega)$ with zero average, and $(\dot{H}^{1}(\Omega))^{*}$ is its
dual. Substituting (C.2) into the KWC functional $\mathcal{W}^{\rm G}$, the
minimization problem in (3.1) transforms to the following saddle point
problem:

$\inf_{\eta\in
L^{2}(\Omega)}\,\sup_{\psi\in(\dot{H}^{1}(\Omega))^{*}}\Phi[\eta,\psi],$ (C.3)

where

$\Phi[\eta,\psi]=-\frac{1}{2\epsilon}\|\psi\|^{2}_{(\dot{H}^{1}(\Omega))^{*}}+\int_{\Omega}(\eta\psi+f(\eta))\,dV+\int_{\mathcal{S}}g(\eta)\mathcal{J}(\llbracket\theta\rrbracket)\,dS.$
(C.4)

The problems of minimizing $\Phi$ with respect to $\eta$, and maximizing it
with respect to $\psi$ are referred to as $\eta$ and $\psi$ sub-problems
respectively. The advantage of using a primal-dual algorithm is evident from
the observation that $\Psi$ _does not_ depend on the gradients of $\eta$,
which renders the $\eta$ sub-problem local, and the nonlinearity in $g(\eta)$
is no longer a concern. The existence and uniqueness of solutions to the sub-
problems follows from standard convex analysis.

We solve for the saddle point of $\Phi$ using the following primal-dual update
scheme (Algorithm 2 in (Chambolle and Pock, 2011)):

$\displaystyle\eta_{n+1}$
$\displaystyle=\operatornamewithlimits{arg\;min}_{\eta\in
L^{2}(\Omega)}\left[\Phi(\eta,\psi_{n})+\frac{1}{2\tau_{n}}\|\eta-\eta_{n}\|_{L^{2}(\Omega)}^{2}\right],$
(C.5a) $\displaystyle\psi_{n+1}$
$\displaystyle=\operatornamewithlimits{arg\;max}_{\psi\in(\dot{H}^{1}(\Omega))^{*}}\left[\Phi(\tilde{\eta}_{n+1},\psi)-\frac{1}{2\sigma_{n}}\|\psi-\psi_{n}\|^{2}_{L^{2}(\Omega)}\right],$
(C.5b)

where

$\tilde{\eta}_{n+1}=(1+\mu_{n})\eta_{n+1}-\mu_{n}\eta_{n}$

with

$\mu_{n}=1/\sqrt{1+2\tau_{n}/\epsilon},\quad\tau_{n+1}=\mu_{n}\tau_{n},\quad\sigma_{n+1}=\sigma_{n}/\mu_{n}.$

The scalars $\tau_{n}$ and $\sigma_{n}$ are the step sizes of the $\eta$\- and
$\psi$-update respectively. The stability (Jacobs et al., 2019; Chambolle and
Pock, 2011) of the update scheme in (C.5) is guaranteed if
$\tau_{n}\sigma_{n}\leq 1$. We select
$\tau_{0}=\epsilon,\sigma_{0}=1/\epsilon$. The solution to (C.5) is obtained
by solving the following Euler–Lagrange equations corresponding to gradient
flows of the two functionals in (C.5) 181818 In order to obtain (C.7), we note
that the constrained gradient in $(\dot{H}^{1}(\Omega))^{*}$ of
$\int_{\Omega}\tilde{\eta}_{n+1}\psi\,dV$ with respect to $\psi$ is
$-\triangle\tilde{\eta}_{n+1}$.

$\left(\frac{1}{\epsilon}+\frac{1}{\tau_{n}}\right)\eta^{2}(x)+\left(\psi_{n}(x)-\frac{2}{\epsilon}-(1+\eta_{n})\frac{1}{\tau_{n}}\right)\eta(x)-\mathcal{J}^{\star}\left([\\![\theta]\\!]\right)+\frac{1}{\epsilon}-\psi_{n}+\frac{1}{\tau_{n}}\eta_{n}=0,$
(C.6)
$(1/\epsilon-\Delta/\sigma_{n+1})\psi_{n+1}=-\Delta(\bar{\eta}_{n+1}+\psi_{n}/\sigma_{n+1}),$
(C.7)

where the surface measure $\mathcal{J}(\llbracket\theta\rrbracket)\,dS$ has
been replaced by a volume measure
$\mathcal{J}^{\star}\,dV=\mathcal{J}(\llbracket\theta\rrbracket)\exp(-x^{2}/2\epsilon^{4})\,dV$
that depends on the distance $x$ from the grain boundary. From (C.6), we note
that the primal dual algorithm along with the choice $g(\eta)=-\log(1-\eta)$
not only renders the $\eta$ sub-problem local but also analytically solvable.

We solve (C.6) and (C.7) on a uniform grid of size $N=N_{x}\times N_{y}$.
Since (C.6) is solved analytically at each grid point, its cost remains
$\mathcal{O}(N)$. We solve for $\psi_{n+1}$ in (C.7) using the fast Fourier
transform (FFT), resulting in an $\mathcal{O}(N\log{N})$ complexity for the
primal dual algorithm. We use the following stopping criterion for the update
scheme in (C.5),

$\|\eta_{n+1}-\eta_{n}\|_{\infty}=\max_{1\leq j\leq
N}|(\eta_{n+1})_{j}-(\eta_{n})_{j}|\leq\mathsf{e},$ (C.8)

where $\mathsf{e}$ is the tolerance of the iterative scheme. Finally, we note
that the use of FFT to solve (C.7) necessitates periodic boundary conditions
on $\eta$. On the other hand, for Neumann boundary conditions, we use the
discrete cosine transform given by

$\hat{\psi}_{pq}=\lambda_{p}\lambda_{q}\sum^{N_{x}-1}_{i=0}\sum^{N_{y}-1}_{j=0}\psi\left(\frac{i}{N_{x}},\frac{j}{N_{y}}\right)\cos{\left(\frac{\pi(2i+1)p}{2N_{x}}\right)}\cos{\left(\frac{\pi(2j+1)q}{2N_{y}}\right)},\quad\begin{aligned}
&0\leq p\leq N_{x}-1\\\ &0\leq q\leq N_{y}-1,\end{aligned}$ (C.9)

with

$\lambda_{p}=\left\\{\begin{aligned} &1/\sqrt{N_{x}},\quad p=0,\\\
&\sqrt{2/N_{x}},\quad 1\leq p\leq N_{x}-1,\end{aligned}\right.\text{ and
}\quad\lambda_{q}=\left\\{\begin{aligned} &1/\sqrt{N_{y}},\quad q=0,\\\
&\sqrt{2/N_{y}},\quad 1\leq q\leq N_{y}-1.\end{aligned}\right.$ (C.10)

##  Appendix D Note on the derivations of thresholding scheme

In this section, we describe the steps to obtain (3.8) from (3.7). We begin by
separating the domain of integration in (3.7) as

$\int^{l_{0}}_{-\infty}(1-u(l/\epsilon))^{2}\;dl=\int^{0}_{l_{0}}(1-u(l/\epsilon))^{2}\;dl+\int^{+\infty}_{0}(1-u(l/\epsilon))^{2}\;dl.$
(D.1)

Substituting the solution in (3.6) into (D.1), we have

$\frac{(u(0)-1)^{2}}{(\frac{2}{\epsilon}-\kappa)}\exp\left[\left(\frac{2}{\epsilon}-\kappa\right)l_{0}\right]=\frac{(u(0)-1)^{2}}{(\frac{2}{\epsilon}-\kappa)}\left(1-\exp\left[\left(\frac{2}{\epsilon}-\kappa\right)l_{0}\right]\right)+\frac{(u(0)-1)^{2}}{(\frac{2}{\epsilon}+\kappa)}.$
(D.2)

Dividing both sides by $(u(0)-1)^{2}$ and collecting the $l_{0}$ terms, we
obtain

$\frac{2\exp\left[\left(\frac{2}{\epsilon}-\kappa\right)l_{0}\right]}{(\frac{2}{\epsilon}-\kappa)}=\frac{1}{(\frac{2}{\epsilon}-\kappa)}+\frac{1}{(\frac{2}{\epsilon}+\kappa)}=\frac{\frac{4}{\epsilon}}{(\frac{4}{\epsilon^{2}}-\kappa^{2})}.$
(D.3)

Taking a logarithm, we have

$\left(\frac{2}{\epsilon}-\kappa\right)l_{0}=\log{\left(\frac{2/\epsilon}{2/\epsilon+\kappa}\right)}=\log{\left(\frac{1}{1+(\epsilon\kappa)/2}\right)}.$
(D.4)

A Taylor expansion of the right-hand-side of (D.4) with respect to
$\epsilon\kappa/2$ results in

$\left(\frac{2}{\epsilon}-\kappa\right)l_{0}=-\frac{\epsilon\kappa}{2}+\frac{\epsilon^{2}\kappa^{2}}{4}-\frac{\epsilon^{3}\kappa^{3}}{8}+O(\epsilon^{4}\kappa^{4}).$
(D.5)

Multiplying by $\epsilon$ on both sides, we have

$\left(2-\kappa\epsilon\right)l_{0}=-\frac{\epsilon^{2}\kappa}{2}+\frac{\epsilon^{3}\kappa^{2}}{4}+\frac{\epsilon^{4}\kappa^{3}}{8}+O(\epsilon^{4}\kappa^{4}).$
(D.6)

Finally, using the approximation $2-\kappa\epsilon\approx 2$, we get (3.8).

As mentioned in Section 3.2, in practice, the infinite bounds of the integral
in (D.1) are replaced by finite bounds of magnitude $l_{\rm b}$. Under this
change, (D.6) modifies as

$\left(\frac{2}{\epsilon}-\kappa\right)l_{0}=\log{\left(\frac{\frac{2}{\epsilon}}{\frac{2}{\epsilon}+\kappa}+\boxed{\frac{1}{2}\exp{\left(-\left(\frac{2}{\epsilon}-\kappa\right)l_{b}\right)}-\frac{\left(\frac{2}{\epsilon}-\kappa\right)}{2\left(\frac{2}{\epsilon}+\kappa\right)}\exp{\left(-\left(\frac{2}{\epsilon}+\kappa\right)l_{b}\right)}}\;\right)}.$
(D.7)

It can be easily shown that the boxed terms resulting from a finite value of
$d_{b}$ decay exponentially as $\epsilon\to 0$, which leaves (3.8) unchanged.

##  Appendix E Fast marching method

The fast marching method (FMM), developed by Tsitsiklis (1995) is used to
evolve a surface in the outward unit normal direction with a speed
$V(\bm{x})>0$. The fast marching method reformulates a time-dependent initial
value problem describing the evolution of a surface into an equivalent
boundary value formulation. In this section, we summarize the FMM algorithm as
described in Sethian (1996). For illustration, let $\bm{s}(t)$ describe a
surface evolving with speed $\mathcal{V}$ from a given initial surface
$\bm{s}(0)=\Gamma$. Instead of solving a time-dependent problem for
$\bm{s}(t)$, the fast marching method solves for a function $\zeta(\bm{x})$
which represents the time it takes for the surface to reach $\bm{x}$. By the
definition of $\zeta$, we have

$\zeta(\bm{s}(t))=t,$ (E.1)

with $\zeta=0$ on $\Gamma$. Differentiating (E.1) with respect to $t$, and
noting that $\nabla\zeta$ is normal to the surface, we arrive at the following
boundary value problem

$|\nabla\zeta|V=1,\quad\zeta=0\text{ on }\Gamma,$ (E.2)

commonly referred to as the Eikonal equation.

Next, we describe the algorithm to solve (E.2) on a two-dimensional grid. In
order to compute $|\nabla\zeta|$, an operator $D^{-x}_{ij}$, representing the
standard backward finite difference operation on the grid point $ij$, is
defined as

$D^{-x}_{ij}\zeta=\frac{\zeta_{ij}-\zeta_{(i-1)j}}{\Delta x}.$ (E.3)

Similarly, $D^{+x}$, $D^{-y}$, and $D^{+y}$ denote forward in $x$, backward
and forward in $y$ finite difference operators respectively. To guarantee a
unique viscosity solution191919See Sethian (1996) on the reason behind seeking
a viscosity solution. of the evolving surface, one should consider an upwind
finite difference scheme to compute the gradient, which is conveniently
written as

$\displaystyle|\nabla\zeta|$
$\displaystyle\approx\left[(\mathrm{max}(D^{-x}_{ij}\zeta,0)^{2}+\mathrm{min}(D^{+x}\zeta_{ij},0)^{2}+(\mathrm{max}(D^{-y}_{ij}\zeta,0)^{2}+\mathrm{min}(D^{+y}\zeta_{ij},0)^{2}\right]^{1/2}$
(E.4)
$\displaystyle=\left[(\mathrm{max}(D^{-x}_{ij}\zeta,0)^{2}+\mathrm{max}(-D^{+x}\zeta_{ij},0)^{2}+(\mathrm{max}(D^{-y}_{ij}\zeta,0)^{2}+\mathrm{max}(-D^{+y}\zeta_{ij},0)^{2}\right]^{1/2}.$

Using (E.4), we rewrite (E.2) in an algebraic form

$\left[(\mathrm{max}(D^{-x}_{ij}\zeta,D^{+x}\zeta_{ij},0)^{2}+(\mathrm{max}(D^{-y}_{ij}\zeta,-D^{+y}\zeta_{ij},0)^{2}\right]^{1/2}=\frac{1}{V(x,y)}.$
(E.5)

Note that if the neighboring values of $\zeta_{ij}$ are known, then (E.5) is a
quadratic equation for $\zeta_{ij}$ that can be solved analytically.

The fast marching method begins with the following initialization step

  1. 1.

Assign $\zeta(x)=0$ for grid points in the area enclosed by the initial
surface, and tag them as accepted.

  2. 2.

Assign $\zeta(x)=+\infty$ for the remaining grid points, and tag them as far.

  3. 3.

Among the accepted points, identify the points that are in the neighborhood of
points tagged as far, and tag them as considered.

The key step in the fast marching method is to update $\zeta$ with a trial
value using (E.5) for grid points tagged as considered , but only accept the
update with the smallest value. In order to identify the smallest value
efficiently, the grid points tagged as considered are stored in a min-
heap202020 A min-heap structure is a complete binary tree with a property that
the value at any given node is less than or equal to the values at its
children.  structure (Sedgewick and Wayne, 2008) borrowed from discrete
network algorithms. The fast marching method then proceeds as follows.

  1. 1.

Construct a min-heap structure for the considered points.

  2. 2.

Access the root (minimum value) of the heap.

  3. 3.

Find a trial solution $\tilde{\zeta}$ on the neighbors of the root using
(E.5). If the trial solution is smaller than the present values, then update
$\zeta(x)=\tilde{\zeta}$.

  4. 4.

If a point, previously tagged as far, is updated using a trial value, relabel
it as considered, and add it to the heap structure.

  5. 5.

Tag the root of the heap as accepted, and delete it from the heap.

  6. 6.

Repeat steps 2 to 5, until every grid point is tagged as accepted.

Figure 20: The level sets of the solution to the Eikonal equation (E.2),
computed using the fast marching method, describe a surface evolving with
outward normal velocity $V(x,y)=1$.

Fig. 20 demonstrates the fast marching method used to track an initial surface

$(9x^{2}-1)^{2}-(3y+1)(1-3y)^{3}=0,$ (E.6)

growing with a uniform outward normal velocity $V(x)=1$.

## References

## References

  * Admal et al. (2018) N. C. Admal, G. Po, J. Marian,  A unified framework for polycrystal plasticity with grain boundary evolution,  International Journal of Plasticity 106 (2018) 1–30.
  * Kobayashi et al. (1998) R. Kobayashi, J. A. Warren, W. C. Carter,  Vector-valued phase field model for crystallization and grain boundary formation,  Physica D: Nonlinear Phenomena 119 (1998) 415–423.
  * Kobayashi et al. (2000) R. Kobayashi, J. A. Warren, W. C. Carter,  A continuum model of grain boundaries,  Physica D: Nonlinear Phenomena 140 (2000) 141–150.
  * Runnels et al. (2016a) B. Runnels, I. J. Beyerlein, S. Conti, M. Ortiz,  An analytical model of interfacial energy based on a lattice-matching interatomic energy,  Journal of Mechanics and Physics of Solids 89 (2016a) 174–193.
  * Runnels et al. (2016b) B. Runnels, I. J. Beyerlein, S. Conti, M. Ortiz,  A relaxation method for the energy and morphology of grain boundaries and interfaces,  Journal of Mechanics and Physics of Solids 94 (2016b) 388–408.
  * Watanabe (2011) T. Watanabe,  Grain boundary engineering: historical perspective and future prospects,  Journal of Material Science 46 (2011) 4095–4115.
  * Mullins (1956) W. W. Mullins,  Two-dimensional motion of idealized grain boundaries,  Journal of Applied Physics 27 (1956) 900–904.
  * Chen et al. (2020) K. Chen, J. Han, X. Pan, D. J. Srolovitz,  The grain boundary mobility tensor,  Proceedings of the National Academy of Sciences of the United States of America 117 (2020) 4533–4538.
  * Olmsted et al. (2009) D. L. Olmsted, S. M. Foiles, E. A. Holm,  Survey of computed grain boundary properties in face-centered cubic metals: I. grain boundary energy,  Acta Materialla 57 (2009) 3694–3703.
  * Bulatov et al. (2013) V. V. Bulatov, B. W. Reed, M. Kumar,  Anisotropy of interfacial energy in five dimensions,  arXiv: Material Science (2013).
  * Rollett (2018) A. Rollett,  One crystal out of many,  Science 362 (2018) 996–998.
  * Barmak et al. (2013) K. Barmak, E. Eggeling, D. Kinderlehrer, R. Sharp, S. Ta’asan, A. D. Rollett, K. R. Coffey,  Grain growth and the puzzle of its stagnation in thin films: The curious tale of a tail and an ear,  Progress in Materials Science (2013) 987–1055.
  * Upmanyu et al. (1998) M. Upmanyu, R. W. Smith, D. J. Srolovitz,  Atomicstic simulation of curvature driven grain boundary migration,  Interface Science 6 (1998) 41–58.
  * Janssens et al. (2006) K. G. F. Janssens, D. Olmsted, E. A. Holm, S. M. Foiles, S. J. Plimton, P. M. Derlet,  Computing the mobility of grain boundaries,  Nature Materials 5 (2006) 124–127.
  * Anderson et al. (1958) M. P. Anderson, D. J. Srolovitz, G. S. Grest, P. S. Sahni,  Computer simulation of grain growth I,  Acta Metallurgica 32 (1958) 783–791.
  * Anderson et al. (1989) M. P. Anderson, G. S. Grest, D. J. Srolovitz,  Computer simulation of grain growth in three dimensions,  Philosophical Magazine B 59 (1989) 293–329.
  * Mendelev and Srolovitz (2002) M. I. Mendelev, D. J. Srolovitz,  Co-segregation effects on boundary migration,  Interface Science 10 (2002) 191–199.
  * Upmanyu et al. (2002) M. Upmanyu, G. N. Hassold, A. Kazaryan, E. A. Holm, Y. Wang, B. Patton, D. J. Srolovitz,  Boundary mobility and energy anisotropy effects on microstructural evolution during grain growth,  Interface Science 10 (2002) 201–216.
  * Yang et al. (2000) Z. Yang, S. Sista, J. W. Elmer, T. DebRoy,  Three dimensional Monte carlo simulation of grain growth during GTA welding of titanium,  Acta Materialla 48 (2000) 4813–4825.
  * Hillert (1965) M. Hillert,  On the theory of normal and abnormal grain growth,  Acta Metallurgica 13 (1965) 227–238.
  * Allen and Cahn (1979) S. M. Allen, J. W. Cahn,  A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening,  Acta Metallurgica 27 (1979) 1085–1095.
  * Frost et al. (1990) H. J. Frost, C. V. Thompson, D. T. Walton,  Simulation of thin film grain structures—I. grain growth stagnation,  Acta Metallugica et Materialia 38 (1990) 1455–1462.
  * Frost et al. (1988) H. J. Frost, C. V. Thompson, C. L. Howe, J. Whang,  A two dimensional computer simulation of capillarity-driven grain growth: Preliminary results,  Scripta Metallurgica 22 (1988) 65–70.
  * Kinderlehrer et al. (2004) D. Kinderlehrer, J. H. Lee, I. Livshits, A. D. Rollett, S. Ta’asan,  Mesoscale simulation of grain growth,  Materials Science Forum 467–470 (2004).
  * Kinderlehrer et al. (2006) D. Kinderlehrer, I. Livshits, S. Ta’asan,  A variational approach to modeling and simulation of grain growth,  SIAM Journal on Scientific Computing 28 (2006) 1694–1715.
  * Zhao et al. (1996) H.-K. Zhao, T. Chan, B. Merriman, S. Osher,  A variational level set approach to multiphase motion,  Journal of Computational Physics 127 (1996) 179–195.
  * Fausty et al. (2018) J. Fausty, N. Bozzolo, D. P. Munoz, M. Bernacki,  A novel level-set finite element formulation for grain growth with heterogeneous grain boundary energies,  Materials & Design 160 (2018) 578–590.
  * Merriman et al. (1992) B. Merriman, J. K. Bence, S. J. Osher,  Diffusion generated motion by mean curvature,  Proceedings of the Computational Crystal Growers Workshop (1992) 72–83.
  * Esedoḡlu and Otto (2015) S. Esedoḡlu, F. Otto,  Threshold dynamics for networks with arbitrary surface tensions,  Communications on Pure and Applied Mathematics 68 (2015) 808–864.
  * Elsey et al. (2009) M. Elsey, S. Esedoḡlu, P. Smereka,  Diffusion generated motion for grain growth in two and three dimensions,  Journal of Computational Physics 228 (2009) 8015–8033.
  * Elsey et al. (2011) M. Elsey, S. Esedoḡlu, P. Smereka,  Large scale simulations and parameter study for a simple recrystallization model,  Philosophical Magazine 91 (2011) 1607–1642.
  * Zaitzeff et al. (2020) A. Zaitzeff, S. Esedoḡlu, K. Garikipati,  Second order threshold dynamics schemes for two phase motion by mean curvature,  Journal of Computational Physics 410 (2020) 109404.
  * Basak and Gupta (2014) A. Basak, A. Gupta,  A two-dimensional study of coupled grain boundary motion using the level set method,  Modelling and Simulation in Materials Science and Engineering 22 (2014).
  * Hallberg and Bulatov (2019) H. Hallberg, V. V. Bulatov,  Modeling of grain growth under fully anisotropic grain boundary energy,  Modelling and Simulation in Materials Science and Engineering 27 (2019) 045002\.
  * Joshi et al. (2020) T. Joshi, R. Arora, A. Basak, A. Gupta,  Equilibrium shape of misfitting precipitates with anisotropic elasticity and anisotropic interfacial energy,  Modelling and Simulation in Materials Science and Engineering 28 (2020).
  * Jokisaari et al. (2017) A. M. Jokisaari, P. W. Voorhees, J. E. Guyer, J. Warren, O. G. Heinonen,  Benchmark problems for numerical implementations of phase field models,  Computational Materials Science 126 (2017) 139–151.
  * Chen (2002) L.-Q. Chen,  Phase-field models for microstructure evolution,  Annual Review of Materials Research 32 (2002) 113–140.
  * Hirouchi et al. (2012) T. Hirouchi, T. Tsuru, Y. Shibutani,  Grain growth prediction with inclination dependence of [110] tilt grain boundary using multi-phase-field model with penalty for multiple junctions,  Computational Materials Science 53 (2012) 474–482.
  * Steinbach (2009) I. Steinbach,  Phase-field models in materials science,  Modelling and Simulation in Materials Science and Engineering 17 (2009) 073001\.
  * Warren et al. (2003) J. A. Warren, R. Kobayashi, A. E. Lobkovsky, W. C. Carter,  Extending phase field models of solidification to polycrystalline materials,  Acta Materialla 51 (2003) 6035–6058.
  * Fan et al. (2002) D. Fan, S. P. Chen, L.-Q. Chen, P. W. Voorhees,  Phase-field simulation of 2-D ostwald ripening in the highvolume fraction regime,  Acta Materialla 50 (2002).
  * Permann et al. (2016) C. J. Permann, M. R. Tonks, B. Fromm, D. R. Gaston,  Order parameter re-mapping algorithm for 3D phase field model of grain growth using fem,  Computational Materials Science 115 (2016).
  * Ribot et al. (2019) J. G. Ribot, V. Agrawal, B. Runnels,  A new approach for phase field modeling of grain boundaries with strongly nonconvex energy,  Modelling and Simulation in Materials Science and Engineering 27 (2019).
  * Moelans et al. (2008) N. Moelans, B. Blanpain, P. Wollants,  Quantitative phase-field approach for simulating grain growth in anisotropic systems with arbitrary inclination and misorientation dependence,  Physical Review Letters 101 (2008) 025502.
  * Kim et al. (2014) H.-K. Kim, S. G. Kim, W. Dong, I. Steinbach, B.-J. Lee,  Phase-field modeling for 3d grain growth based on a grain boundary energy database,  Modelling and Simulation in Materials Science and Engineering 22 (2014) 034004\.
  * Read and Shockley (1950) W. T. Read, W. Shockley,  Dislocation models of crystal grain boundaries,  Physical Review 78 (1950) 275.
  * Holm et al. (2010) E. A. Holm, D. L. Olmsted, S. M. Foiles,  Comparing grain boundary energies in face-centered cubic metals: Al, Au, Cu and Ni,  Scripta Materialia 63 (2010) 905–908.
  * Bulatov et al. (2014) V. V. Bulatov, B. W. Reed, M. Kumar,  Grain boundary energy function for fcc metals,  Acta Materialla 65 (2014) 161–175.
  * Wolf (1990) D. Wolf,  Structure-energy correlation for grain boundaries in fcc metals—III. Symmetrical tilt boundaries,  Acta Metallugica et Materialia 38 (1990) 781–790.
  * Mason and Patala (2019) J. K. Mason, S. Patala,  Basis functions on the grain boundary space: Theory,  arXiv preprint arXiv:1909.11838 (2019).
  * Kim et al. (2014) H.-K. Kim, S. G. Kim, W. Dong, I. Steinbach, B.-J. Lee,  Phase-field modeling for 3D grain growth based on a grain boundary energy database,  Modelling and Simulation in Materials Science and Engineering 22 (2014) 034004\.
  * Alicandro et al. (1999) R. Alicandro, A. Braides, J. M. Shah,  Free-discontinuity problems via functionals involving the L1-norm of the gradient and their approximations,  Interface and Free Boundaries 1 (1999) 17–37.
  * Lobkovsky and Warren (2001) A. E. Lobkovsky, J. A. Warren,  Sharp interface limit of a phase-field model of crystal grains,  Physical Review E 63 (2001) 051605.
  * Dorr et al. (2010) M. R. Dorr, J.-L. Fattebert, M. E. Wickett, J. F. Belak, P. E. A. Turchi,  A numerical algorithm for the solution of a phase-field model of polycrystalline materials,  Journal of Computational Physics 229 (2010) 626–641.
  * Admal et al. (2019) N. C. Admal, J. Segurado, J. Marian,  A three-dimensional misorientation axis- and inclination-dependent Kobayashi–Warren–Carter grain boundary model,  Journal of the Mechanics and Physics of Solids 128 (2019) 32–53.
  * Jacobs et al. (2019) M. Jacobs, F. Leger, W. Li, S. Osher,  Solving large-scale optimization problems with a convergence rate independent of grid size,  SIAM Journal on Numerical Analysis 57 (2019) 1100–1123.
  * Salvador and Esedoḡlu (2019) T. Salvador, S. Esedoḡlu,  The role of surface tension and mobility model in simulations of grain growth,  arXiv:1907.11574 (2019).
  * Martine La Boissonière et al. (2019) G. Martine La Boissonière, R. Choksi, K. Barmak, S. Esedoḡlu,  Statistics of grain growth: Experiment versus the phase-field-crystal and mullins models,  Materialia 6 (2019) 100280.
  * Tsitsiklis (1995) J. N. Tsitsiklis,  Efficient algorithms for globally optimal trajectories,  IEEE Transactions on Automatic Control 40 (1995) 1528–1538.
  * Herring (1951) C. Herring, Surface Tension as a Motivation for Sintering, McGraw Hill, 1951. doi:10.1007/978-3-642-59938-5_2.
  * Esedoḡlu et al. (2010) S. Esedoḡlu, S. Ruuth, R. Tsai,  Diffusion generated motion using signed distance functions,  Journal of Computational Physics 229 (2010) 1017–1042.
  * von Neumann (1952) J. von Neumann, Metal Interfaces, American Society for Metals, Cleveland, 1952.
  * Alnaes et al. (2015) M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, G. N. Wells,  The FEniCS Project Version 1.5,  Archive of Numerical Software 3 (2015).
  * Saad and Schultz (1986) Y. Saad, M. H. Schultz,  GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems,  SIAM Journal on Scientific Computing 7 (1986) 856–869.
  * Li and Saad (2013) R. Li, Y. Saad,  GPU-accelerated preconditioned iterative linear solvers,  The Journal of Supercomputing 63 (2013) 443–466.
  * Govindaraju et al. (2008) N. K. Govindaraju, B. . Lloyd, Y. Dotsenko, B. Smith, J. Manferdelli,  High performance discrete fourier transforms on graphics processors,  Proceedings of the 2008 ACM/IEEE Conference on Supercomputing (2008).
  * Li (1962) J. C. Li,  Possibility of subgrain rotation during recrystallization,  Journal of Applied Physics 33 (1962) 2958–2965.
  * Thomas et al. (2017) S. Thomas, K. Chen, J. Han, P. K. Purohit, D. J. Srolovitz,  Reconciling grain growth and shear-coupled grain boundary migration,  Nature communications 8 (2017) 1–12.
  * Wei et al. (2020) C. Wei, L. Zhang, J. Han, D. J. Srolovitz, Y. Xiang,  Grain boundary triple junction dynamics: a continuum disconnection model,  SIAM Journal on Applied Mathematics 80 (2020) 1101–1122.
  * Runnels and Agrawal (2020) B. Runnels, V. Agrawal,  Phase field disconnections: A continuum method for disconnection-mediated grain boundary motion,  Scripta Materialia 186 (2020) 6–10.
  * Miura et al. (1994) H. Miura, M. Kato, T. Mori,  Temperature dependence of the energy of Cu [110] symmetrical tilt grain boundaries,  Journal of Material Science Letters 13 (1994) 46–48.
  * Powell (1978) M. J. D. Powell,  Algorithms for nonlinear constraints that use lagrangian functions,  Mathematical Programming 14 (1978) 224–248.
  * Kuhn (1955) H. W. Kuhn,  The hungarian method for the assignment problem,  Naval Research Logistics Quarterly 2 (1955) 83–97.
  * Komodakis and Pesquet (2015) N. Komodakis, J. Pesquet,  Playing with duality: An overview of recent Primal-dual approaches for solving large-scale optimization problems,  IEEE Signal Processing Magazine 32 (2015) 31–54.
  * Donoho (2006) D. L. Donoho,  Compressed sensing,  IEEE Transactions on Information Theory 52 (2006) 1289–1306.
  * Chambolle and Pock (2011) A. Chambolle, T. Pock,  A first-order primal-dual algorithm for convex problems with applications to imaging,  Journal of Mathematical Imaging and Vision 40 (2011) 120–145.
  * Shalev-Shwartz and Singer (2007) S. Shalev-Shwartz, Y. Singer,  A primal-dual perspective of online learning algorithms,  Machine Learning 69 (2007) 115–142.
  * Combettes and Pesquet (2011) P. L. Combettes, J.-C. Pesquet, Proximal Splitting Methods in Signal Processing, volume 49 of Springer Optimization and Its Applications, Springer, New York, 2011. doi:10.1007/978-1-4419-9569-8_10.
  * Sethian (1996) J. Sethian,  A fast marching level set method for monotonically advancing fronts,  Proceedings of the National Academy of Sciences of the United States of America 93 (1996) 1591–1595.
  * Sedgewick and Wayne (2008) R. Sedgewick, K. Wayne, Algorithms, Addison-Wesley, 2008\.