Rapid multiphasefield model development using a modular free energy based approach with automatic differentiation in MOOSE/MARMOT
Abstract
We present a novel phasefield model development capability in the open source MOOSE finite element framework. This facility is based on the “modular free energy” approach in which the phasefield equations are implemented in a general form that is logically separated from modelspecific data such as the thermodynamic free energy density and mobility functions. Free energy terms contributing to a phasefield model are abstracted into selfcontained objects that can be dynamically combined at simulation run time. Combining multiple chemical and mechanical free energy contributions expedites the construction of coupled phasefield, mechanics, and multiphase models. This approach allows computational material scientists to focus on implementing new material models, and to reuse existing solution algorithms and data processing routines. A key new aspect of the rapid phasefield development approach that we discuss in detail is the automatic symbolic differentiation capability. Automatic symbolic differentiation is used to compute derivatives of the free energy density functionals, and removes potential sources of human error while guaranteeing that the nonlinear system Jacobians are accurately approximated. Through justintime compilation, we greatly reduce the computational expense of evaluating the differentiated expressions. The new capability is demonstrated for a variety of representative applications.
keywords:
phasefield, finite element, automatic differentiationPacs:
, 46.15.x, 05.10.a, 02.70.DhMsc:
[2010] 6504, 65Z05[correspondingauthor]Corresponding author
1 Introduction
The phasefield method is a wellestablished tool for simulating the coevolution of microstructure and physical properties at the mesoscale (1); (2). In the phasefield method, the microstructure is described by a system of continuous variables, also called order parameters. Microstructure interfaces are approximated using a finite width, and the order parameters vary smoothly over the interfaces. In isolated systems, the evolution of these variables leads to a monotonically decreasing free energy as a function of time. The phasefield method has been used to model a large range of physical phenomena, including solidification (3); (4), phase transformation (5); (6), and grain growth (7); (8).
The microstructure evolution is governed by partial differential equations (PDEs) that depend on the derivatives of a functional defining the free energy of the system in terms of the phasefield variables. The phasefield PDEs have been solved using the finite difference (9), finite volume (10), and finite element (FEM) methods (11); (12), as well as spectral methods based on the fast Fourier transform (13). The finite difference method is the easiest to implement, while the spectral methods are computationally efficient for small systems and offer better convergence properties. The finite volume and finite element methods are the most flexible, allowing a large range of boundary conditions, domain shapes, and coupling to other physics.
No matter which method is used to solve the phasefield equations, their basic form remains the same for most models: only the free energy functional and the mobility expressions change when modeling different materials and physical phenomena. In this work, we continue the efforts described in (12) where a free energy based approach to phasefield model development is employed. This strategy takes advantage of the “fixed” form of the governing phasefield equations, and restricts the task of programming new phasefield models to the task of programming new free energy functionals. To the authors’ knowledge, this specific approach is not currently followed by other open source phasefield software development efforts (14); (15); (10); (16), however it is also a fairly natural one that could be implemented retroactively in existing code bases.
The free energy based approach employs three capabilities which simplify and accelerate the development of new multiphasefield models: (i) a symbolic differentiation module which eliminates the need to compute derivatives (up to thirdorder) of the free energy, (ii) a system to modularize, recombine, and reuse the various free energy contributions across models, and (iii) a generic framework for multiphasefield simulations in which new free energy modules can be employed. The sum of these components represents an important step toward a modular and general phasefield approach.
2 PhaseField method summary
In the phasefield method, the evolution of nonconserved order parameters (e.g. phase regions and grains) is governed by the Allen–Cahn (17) equation (1) and the evolution of conserved order parameters (e.g. concentrations) is governed by the Cahn–Hilliard (18) equation (2):
(1)  
(2) 
Here, is the total free energy of the system, which can be formulated as a volume integral
(3) 
where is the simulation domain,
(4) 
is the local free energy density, and
(5) 
is the gradient energy contribution. is the free energy contribution due to external driving forces. In general, the total free energy depends on all of the conserved and nonconserved order parameters and their gradients. Computing the variational derivatives in (1) and (2) yields terms which depend on the derivatives of the local free energy density, , with respect to all order parameters, i.e.
(6)  
(7) 
Note that the thermodynamics of the modeled system are determined by the local free energy density , while the gradient contribution produces diffusewidth interfaces and contributes only to the interfacial energy. The gradient contribution is often a known functional that exposes only scalar parameters for tuning the interfacial width, as discussed in (2). is therefore the primary input needed to formulate a new phasefield material model.
In this work, we solve the resulting PDEs using FEM and implicit time integration with the open source Multiphysics ObjectOriented Simulation Environment (MOOSE) (19); (12). MOOSE is a finiteelement framework primarily developed at Idaho National Laboratory (INL) which includes several physics modules that assist users in developing phasefield, thermal transport, solid mechanics, and chemistry models. All the techniques and applications discussed in this work are currently available in the MOOSE phase_field module. MOOSE solves systems of PDEs in a tightlycoupled manner by forming a single residual vector comprising all the unknowns. To form the residual vector, the equations have to be transformed into their weak form via multiplication by a suitably defined test function and integration over . We subsequently utilize the Gauss divergence theorem to reduce the order of the spatial derivatives in the resulting residual equations (see (12) for more detail on the development of the weak form).
Due to the use of implicit time integration, the system of phasefield equations requires a nonlinear solve at each time step. To solve this system of nonlinear equations, MOOSE typically employs the preconditioned Jacobianfree Newton Krylov (20); (21) method (PJFNK), via interfaces provided by the libMesh (22) and PETSc (23) libraries. To improve the convergence properties of the nonlinear solve, the preconditioning matrix should be as close as possible to the actual Jacobian of the nonlinear system of equations. Computing the Jacobian matrix entries requires derivatives of the residual vector entries with respect to each of the nonlinear variables in the system. This includes, for instance, crossderivatives of the free energy density functional with respect to all phasefield variables used in the model.
The Cahn–Hilliard Eq. (7) involves a fourthorder spatial derivative on the concentration variable. There are two standard formulations which are commonly used to solve (7), see, for instance, the discussion in (24). The first is to directly solve the timedependent fourthorder PDE (using a continuous finite element discretization) and the second is to split the equation into two secondorder PDEs which can be solved with more traditional finite elements. The split solve of the Cahn–Hilliard equation depends on , while the the direct solve involves
(8) 
If there are phasefield variables, computing the Jacobian requires additional second derivatives for the Allen–Cahn equation/split solve of the Cahn–Hilliard equation, and third derivatives for the direct solve of the Cahn–Hilliard equation. The MOOSE phase_field module contains FEM discretizations for both the split and nonsplit Cahn–Hilliard formulations, since there are advantages and disadvantages to using both (24).
3 Free energy based approach
The development of numerical simulation tools for solving phasefield mathematical models can be simplified by employing an approach based on the free energy functional. In such an approach, the Cahn–Hilliard equation (both the nonsplit and split forms) and the Allen–Cahn equation obtain all the required free energy density derivatives from free energy objects (in the sense of objectoriented programming) which are specialized to the application in question. When a new phasefield model is developed, the existing implementations of the residual equations are reused—only the code required to define the specialized free energy object itself must be written.
In the MOOSE phase_field module, Cahn–Hilliard and Allen–Cahn residuals and Jacobians have been implemented using small program units called “Kernels”. A Kernel is a C++ class representing a term in the weak form of a PDE. The computation of the free energy density and its derivatives is performed in MOOSE by “Material” objects. A Material computes values, which may depend on nonlinear variables, at specified points in the simulation domain. All simulationspecific physics are implemented through Material objects which calculate the free energy density. The free energy density derivatives are then used in a generic way by the Kernel objects.
4 Automatic differentiation of free energies
While the free energy based approach to phasefield model implementation simplifies the overall process, manually programming the free energy derivatives can be both burdensome and a potential source of error. Therefore, automatic differentiation is used to simplify the process even further, requiring the user to implement only the free energy density expression itself. All required derivatives are computed analytically in a fullyautomated way.
To facilitate automatic differentiation, and to allow userdefined functions to be supplied via input files, MOOSE employs the Function Parser library (25) that is included as a thirdparty plugin in the underlying libMesh finite element library. The Function Parser library accepts a mathematical function definition given as a plain text string. The expression string is lexically parsed into an intermediate tree representation and then transformed into stack machine bytecode. This bytecode can then be executed by the Function Parser bytecode interpreter module as often as necessary without further transformation. This intermediate tree representation of the Function Parser expressions, illustrated in Fig. 1, readily lends itself to algorithmic transformations such as automatic differentiation.
The automatic differentiation system, which is a major new contribution in this work, operates on the tree representation of the free energy function. In this tree structure, leaf nodes can correspond to constants or variables, and internal nodes correspond to mathematical operators and functions, the arguments of which are contained in child nodes or subtrees. The derivatives of the leaf nodes are for all nodes that do not represent the variable the derivative is taken with respect to, and for all the nodes that do. The derivatives of the internal nodes are constructed recursively according to a set of elementary derivative rules.
Construction of the derivative starts at the root node of the expression tree. For the example expression tree in Fig. 1, which represents the expression , the root node holds the multiplication . To obtain the derivative with respect to , we need to calculate the derivative of the root node, . We set according to the product rule. This expression contains derivatives of the nodes and . These derivatives are recursively constructed until leaf nodes with a zero derivative value are reached. This happens in all cases except , which evaluates to one. The full derivative expression that is constructed in this manner is .
Our method is a member of the class of source transformation or symbolic differentiation algorithms (26); (27). In contrast to the commonly used forward and backward accumulation automatic differentiation algorithms, the derivative evaluation is not tied to the evaluation of the undifferentiated function. Each derivative is algebraically optimized, compiled, and evaluated exactly when needed. In a finite element code, the weak form residual evaluation (which, in a split form Cahn–Hilliard problem, contains the first derivative of the free energy only) is the most common operation by far, occurring at every linear iteration. The Jacobian, containing the second derivatives, is only evaluated once per nonlinear iteration. The undifferentiated form is not evaluated for the purpose of solving the phasefield equations, but only when needed for output or postprocessing stages, which is generally once per time step. MOOSE furthermore detects which derivatives are needed by the residual and Jacobian evaluations and skips evaluation of unused derivatives, which reduces the computational burden when constructing approximate Jacobian matrices for preconditioning purposes.
The Function Parser library provides a comprehensive algebraic optimizer that groups, reorders, and transforms the function expression into an equivalent but fastertoevaluate form. The algebraic simplifications are essential for removing the trivial leaf node derivatives which may lead to evaluation errors, such as division by zero, and can be avoided by simple term cancellations. In the above example, the simplifications reduce the derivative expression to .
To further improve the performance of the parsed and runtime interpreted functions, we have also developed a justintime (JIT) compilation module. At runtime, the generated bytecode sequences are automatically transformed into small C source code files. A compiler is dispatched to compile each function file into a dynamically linkable library, which is then loaded “on the fly” with a dlopen (28) POSIX system call. This occurs once during simulation initialization. If the JIT compilation fails, the function evaluation falls back on the bytecode interpreter, otherwise the generated machine code is called. The average time overhead of the additional compilation step is on the order of 0.1 s per function expression or less, depending on the system the simulation is executed on. This is further mitigated by a caching system. A unique SHA1 hash (29) is computed from the function bytecode, and the compiled functions are stored permanently using the hash as a filename. Recompilation will only occur if the bytecode, and thus the function expression, changes. Trivial function changes, namely the modification of constants, will in most cases not trigger a recompilation.
In Fig. 2, the run times for an iron chromium phasefield simulation (described in detail below) are compared for different implementations of the free energy. Results (a) through (e) are obtained using our symbolic differentiation system, while results (f) and (g) are obtained using hand coded derivatives. The outputs generated by those runs have been confirmed to be identical within the limits of the convergence tolerance. The blue (light gray) bars were obtained using the LLVMbased Clang compiler (30); (31) (version 3.9.0), and the red (dark gray) results were obtained using the GNU Compiler Collection (32) (GCC) C++ compiler (version 6.2.0). Results (a) and (b) employ the standard Function Parser bytecode compiler. Result (b) shows the effect of applying the builtin algebraic optimization to the function byte code. Results (c), (d) and (e) use the JIT compilation module developed for this work. JIT compilation turns out to be the crucial step needed to ensure performance on par with handcoded derivatives as shown in result (g).
The JIT module has been incorporated into the Function Parser library fork in libMesh. Through this automatic symbolic differentiation system, we achieve a significant reduction in developer time and remove a source of developer error that is difficult to track down and debug. We note that a naive implementation of the free energy obtained via manual differentiation can exhibit significantly worse performance. One factor is the use of the C++ library function std::pow, which for integer powers is significantly slower than repeated multiplication operators. The Function Parser library performs this optimization automatically.
To demonstrate the correctness of the implementation, we compare the numerical solution of a Cahn–Hilliard equation for a simple double well free energy given by (interface parameter equal to 1) to the analytical solution, for the equilibrium interface in a one dimensional domain. Fig. 3 shows good agreement of the interface profiles obtained numerically (red circles) and the analytical solution (solid line). The inset shows the expected quadratic convergence with respect to the mesh resolution for linear Lagrange elements. The split formulation of the Cahn–Hilliard equation was used along with automatic differentiation of the free energy. To the authors’ best knowledge this integrated approach of utilizing run time symbolic differentiation and JIT compiling of free energies in a phasefield framework is unique.
4.1 Smoothlyextrapolated logarithm
Free energies that contain a term for the configurational entropy derived from ideal or regular (33) solution models will contain terms of the form
where is a conserved order parameter. As the natural logarithm is only defined for strictly positive numbers, this expression restricts the domain of the free energy to numbers on the open interval . This poses numerical challenges for systems with equilibrium concentrations near 0 or 1.
To improve the convergence behavior, we have developed a smoothlyextrapolated logarithm surrogate function. For input arguments , we compute a Taylorseries expansion of the logarithm function centered around instead. For , we evaluate the standard logarithm function. This surrogate logarithm function also extends to negative arguments, and is twice differentiable everywhere. In the resulting free energy expressions, the extension to negative arguments manifests as a free energy penalty which drives the solution back to physicallyallowable concentrations without incurring numerical errors. In previous work (34) piecewise constructions using Taylor expansions of the full free energy expressions outside its domain have been suggested. This requires knowledge of the domain of the free energy to set the interval boundaries over which the Taylor expansion is active. Our approach retains a single free energy expression with the Taylor expansion only occurring at the level of the logarithm functions, which are the underlying cause of the limitation of the domain of the free energy. The domain of the log function is well defined and thus our model requires fewer user inputs.
Care has to be taken to choose small enough to not adversely impact the thermodynamic properties of the simulated system. In particular, large values of can change the phase diagram by moving the common tangent points to larger concentrations.
As an example, we used the techniques presented above to implement a phasefield model based on the published (35) free energy surface of an iron chromium binary alloy. We have encoded the full free energy expression from the publication into a MOOSE input file as a parsed function Material with automatic differentiation and Taylor expansion substitution for small logarithm arguments. A running phasefield model was obtained with little additional effort. Fig. 4 shows a simulation snapshot obtained using this free energy. The system is in the particle coarsening stage, having previously undergone spinodal decomposition. A line scan performed on the center precipitate is shown in the inset. In agreement with the published free energy surface and resulting phase diagram, we observe practically no solubility of iron in the chromium precipitates, while the chromium solubility in the iron matrix is approximately 6.7% at the simulation temperature of 500 K.
In addition to the free energy, the user has to provide a mobility model, which for this example we assumed to be concentration independent.
We used experimental data on chromium tracer diffusivity in iron as a guide to set the actual value of the mobility . For the chosen simulation length scale of 20 nm 20nm and mesh of quadrilateral elements, an appropriate interfacial energy parameter eV nm was chosen. The initial choice of is somewhat arbitrary, as its order of magnitude rarely changes at a given length scale, and common excess free energies are of similar magnitude for many alloy systems.
Further refinement of the value to obtain a specific surface energy may require a few simulation trial runs, which on onedimensional test systems take only seconds.
The example simulation was run with an adaptive time stepper to a simulation time of about 24 h.
All input files required to rerun the simulations in this paper are available from the MOOSE GitHub repository
5 Multiphase phasefield models
The tools presented above allow the rapid implementation of single phase material systems, however, many material systems of interest exhibit complex phase diagrams with multiple phases potentially coexisting in a simulation volume. These phases have separate free energies and can potentially have different mechanical and thermal properties. To model systems with multiple coexisting phases, the construction of a global free energy functional spanning the entire phase space of the system is necessary. We now show that the materialbased modular free energy system presented here lends itself to the convenient construction of such multiphase free energies. The free energy of each individual phase can be provided using a Material object that encapsulates all required free energy derivatives and leverages the capabilities of the symbolic differentiation module. The global free energy is then constructed as a combination of the phase free energies using nonconserved order parameters to indicate the phase distribution. Just as the free energy values are combined, the derivatives for the phase free energies are used to construct the derivatives of the global free energy. The MOOSE phasefield module has a selection of premade Material objects for global free energy construction which users can provide single phase free energies for.
One common approach is to use a linear combination of the free energy densities of each phase in the system based on the WBM model (5), e.g.
(9) 
The switching function varies smoothly from to as goes from to . The barrier function (multiplied by the barrier height ) penalizes phase mixtures over pure phase regions. The total weight of all phase free energy contributions at each point in the simulation volume is exactly unity, which can be expressed as the constraint
(10) 
Two phase systems can be modeled using a single order parameter with the explicit constraint . The symmetric switching function then satisfies the constraint (10). For phase systems with , it is advantageous to use order parameters. In this case, the constraint is not automatically satisfied and needs to be enforced by other means. In the MOOSE phasefield module we provide two methods of enforcing the switching function constraint: a “hard” constraint utilizing the Lagrange multiplier technique, and a “soft” constraint implemented via a free energy penalty term.
The soft constraint is implemented by adding the contribution,
(11) 
to the free energy, where is a usertunable penalty coefficient. In contrast, the hard constraint is imposed by introducing a Lagrange multiplier as a field variable. The variational statement of the problem is then: find such that the boundary conditions are satisfied, and
(12)  
(13) 
hold for all admissible test functions , where is the weak form (Allen–Cahn) residual for the th nonconserved order parameter. We note that these equations alter the character of the Jacobian matrix of the nonlinear problem by introducing a zero block on the diagonal. This can complicate the task of preconditioning and iteratively solving the system substantially. By replacing the constraint with the modified constraint
(14) 
the term introduces an dependence in the constraint. This results in a nonzero ondiagonal Jacobian contribution of for Eq. (13), avoiding “zero pivot” errors arising from PETSc preconditioners (such as pc_type lu, which implements only partial pivoting). The value of should be chosen slightly larger than the linear solver tolerance, and results in a trade off between accuracy and solver performance. This approach does result in a violation of the constraint of , however it was found that this discrepancy did not adversely affect the overall quality of the solution, and improved the convergence characteristics of the solver.
5.1 KKS models
An additional multiphase model implemented in the phase_field module is the KimKimSuzuki (KKS) model (6). KKS addresses the issue of systems with large phase free energy differences in the interfacial regions. One relevant example is the xenon gas bubble problem shown in Fig. 5, computed using the phase free energies from Li et al. (36). Here the gas solubility in the solid UO matrix is very low, with large free energy penalties for large gas concentrations. In the bubble phase, the equilibrium gas concentration is near unity. In the bubble/matrix interfacial region, both the order parameter and the concentration change from the bubble equilibrium values to the matrix equilibrium values over a finite distance due to the soft interface approximation. In that interfacial region, the free energy of both phases is computed for the intermediate concentration range, which results in large free energy densities from the solid phase contribution. This increases the interfacial free energy of the bubbles to a nonphysical value.
The KKS model solves this problem by introducing the concept of phase concentrations, which are effectively the fractions of the total concentration held in a given phase. In this model, the gas concentration is largely shifted to the gas phase to avoid the free energy penalty. In the KKS model, the interfacial free energy is decoupled from the diffuse interface width, allowing the use of wider interfaces and therefore coarser spatial discretizations requiring fewer computational resources. However, solving for these new variables requires additional differential equations. In the KKS model, these are given by requiring mass conservation and equality of the component chemical potentials across all phases:
(15)  
(16)  
(17) 
The MOOSE phasefield module currently implements a twophase version of the KKS model that uses Kernels for the phasefield equations as well as the KKS constraint equations. The free energy is supplied to those Kernels using the derivative Material system outlined above in the same manner as for the other multiphase approaches discussed above.
In addition to the twophase KKS model, a threephase model based on (37) has been implemented in the MOOSE phasefield module. One limitation of many multiphasefield models is that a binary interface between two phases is unstable with respect to the spurious formation of additional phases in the diffuse interfacial region (38); (39). The formation of the spurious third phase occurs when the free energy is not convex with respect to that third phase. Its formation distorts the interfacial energy of the binary interface, and can lead to nucleation of the third phase in unphysical locations. This can be mitigated by requiring the switching functions to have the additional property of having zero slope and positive curvature perpendicular to the transformation path between two phases. Functions with that property are called tilting functions. In the threephase model implemented in MOOSE, the tilting functions developed in (38) and applied to the KKS approach in (37) are used to prevent spurious thirdphase formation at a twophase interface.
In this model, the three phases are represented by order parameters , , and . They are constrained such that
(18) 
The free energy of the system is given by
(19) 
The local energy is given by
(20) 
where is the free energy density of each phase, is the phase concentration, is the potential barrier height, and the are the tilting functions (38); (37)
(21) 
for , and where
(22)  
(23) 
i.e. form a cyclic permutation. The gradient energy is given by
(24) 
The concentration is defined as a function of the phase concentrations as
(25) 
Because the tilting function reduces to the commonly used twophase interpolation function along the twophase interfaces (38), this constraint on the concentrations reduces to (16), the constraint in the twophase KKS model, for interfaces. The phase concentrations are constrained such that the chemical potentials of each phase are equal:
(26) 
To enforce the constraint of (18), a term is added to the free energy functional, where is a Lagrange multiplier. This results in the Lagrangian :
(27) 
The time evolution of the order parameters is given by the Allen–Cahn equation
(28) 
where using the variational derivative of the Lagrangian enforces the constraint (18). Using the fact that and assuming the mobilities for each phase are equal at each position (though they may be dependent on the local values of the order parameters), the Lagrange multiplier can be eliminated and the Lagrangian can be written in terms of the nonconstrained variational derivatives as
(29) 
Substituting for ,
(30) 
where
(31) 
and
(32) 
The time evolution of the concentration field is given by the Cahn–Hilliard equation:
(33) 
where is the Cahn–Hilliard mobility. Using the chain rule, it can be shown following Kim et al. (40) that this is equivalent to
(34) 
where is the solute diffusion coefficient. This demonstrates that the evolution equation can be solved in terms of the phase concentrations.
The threephase KKS model was used to simulate a trijunction in a simple model system, as shown in Fig. 6. The three phases have parabolic free energies, supplied using the derivative Material system, with minima at , , and . The free energy of phase 2 was made temperaturedependent, and a fixed temperature gradient was imposed to stabilize the trijunction. As seen in Fig. 6, at equilibrium a stable trijunction is formed. The boundaries between phases are indicated with contour lines at and . No thirdphase formation is observed at any twophase interface. It was also verified separately that the order parameter and composition profiles through a 1D interface match the analytical solutions, , , , and (where is the midpoint of the interface and ), and that the interfacial energy matched the analytical solution, . Thus, the threephase KKS model implemented in the MOOSE phasefield model allows simulation of threephase systems with arbitrary free energies (specified as before with the derivative Material system), prevents thirdphase formation at a twophase interface, decouples interfacial energy from local (chemical) energies, and allows the interfacial energy and thickness to be set independently for optimal efficiency.
6 Mechanics coupling
Heterogeneous material properties in multiphase simulations, such as lattice mismatches and variations in the elasticity tensor, will introduce an interplay between chemical and mechanical driving forces. Thus, it is critical to couple the phasefield model equations to equations defining the mechanical behavior of the material. Mechanics simulations are available in the MOOSE tensor_mechanics module, which can be easily coupled to the phase_field module. The local displacement vector is determined by solving the stress divergence equation
(35) 
where is the stress and is an applied body force. The system is solved given suitable boundary conditions and a constitutive law defining the relationship between stress and the strain, . A stress free strain (or eigenstrain) accounts for lattice mismatches, thermal expansion, etc. The elastic energy of the system
(36) 
is added to the total free energy to account for its impact on the microstructure evolution.
In multiphase models, two approaches exist to model the elastic free energy of the total system:

In the VoightTaylor scheme (41) each phase has its own mechanical properties, including the elastic constants, constitutive model (and stress), and eigenstrain, which can depend on composition variables. Each phase free energy contains an elastic free energy contribution. The global elastic free energy is computed analogously to the total chemical free energy, i.e.
(37) and the total stress is calculated in a similar manner
(38)
In both cases, it is necessary for the mechanical property materials to provide derivatives with respect to their dependent variables. A solute concentration dependent eigenstrain Material, for example, provides the derivative. The Function Parser system currently only operates with scalar quantities. Free energy contributions requiring vector or tensor terms need to be implemented separately. MOOSE provides a component that computes linear elastic energy contributions to the free energy. To maintain a modular free energy approach, we have created materials that define the elastic contributions to the free energy and its derivatives. In the first approach, this elastic contribution is computed for each phase individually, while in the second approach only one elastic free energy computation is performed in the system. Additional material classes combine the energy contributions according to Eq. (37). A utility Material is then used to add the chemical and mechanical free energies into a total phase free energy. Mechanical material properties with a direct dependence on the phasefield variables will automatically yield coupling terms between the mechanics and phasefield equations.
The two approaches result in different elastic energy densities in interfacial regions. This is demonstrated in Fig. 7, which shows the elastic energy density of six spherical particles with varying radii as a function of distance from the particle center. The particle phase has a 5% eigenstrain. Each phase has a parabolic free energy with an equilibrium solute concentration of 1 in the precipitate phase and 0 in the matrix phase. The multiphase model is assembled using (9), with , , and . A single order parameter is used to model the twophase system. The stiffness tensor is symmetric and isotropic with the nonzero entries being 1 eV/nm.
Note that the VTS perphase elastic energy simulations result in a pronounced elastic energy excess at the interface. In both cases, the elastic strain state varies smoothly across the interface. In the Khachaturyan scheme, the effective eigenstrain also varies smoothly across the interface. In the VTS scheme each phase has a fixed eigenstrain value, and in the interface region the matrix is strongly stressed for both phases. This effect, and its consequences for simulating systems with strong elastic anisotropies, will be explored in a future publication.
As a further demonstration of the mechanics coupling capability, we consider an immiscible threephase system consisting of a matrix phase and two precipitate phases with anisotropic eigenstrains simulated with the second (global mechanical properties) approach. Both particles have 5% lattice contraction along their minor axis direction. The elastic properties of precipitates and matrix are set to a bulk modulus of 20 GPa and a shear modulus of 7 GPa. Noflux boundary conditions are applied for the phasefield variables and the null space for the displacement variables was eliminated by pinning select nodes. Fig. 8 shows the simulation state after the precipitate growth has progressed. The mesh displacement is plotted with an amplification factor of 5. The bottom half of the plot shows the local strain energy density in eV/nm. The long range stress field in the minor axis direction enforces the lenticular shape of the particles. Both precipitate phases have a simple harmonic free energy with a minimum at , while the matrix has its chemical free energy minimum at . As an initial condition, two spherical nuclei with and a radius of 2 nm were inserted in a super saturated matrix with to provide solute for particle growth.
7 Conclusions
In this work, we have summarized a novel capability for the rapid development of multiphasefield models using automatic symbolic differentiation in the open source MOOSE framework. A modular free energy based approach allows researchers to focus on material model development without the need to touch the underlying numerical details of the coupled partial differential equation system solves. Encapsulating free energies together with their derivatives, which are needed by the phase field evolution equations, allows them to be recombined like building blocks at runtime to set up simulation scenarios. Both the automatic symbolic differentiation capability and the free energy based approach to the solution of phasefield models have been available in MOOSE for some time.
Symbolic differentiation enhances developer productivity in three important ways: it lowers the bar of entry required to investigate new/experimental phasefield models, it reduces the amount of time computational scientists must spend developing code, allowing that time to instead be spent on analysis, and it prevents an extremely common class of errors, i.e. incorrect Jacobians, from negatively impacting the efficiency and accuracy of results. We have shown that the performance of the automatically generated symbolic derivatives is at least on par with carefully handcrafted code when using the provided justintime compilation capability.
Leveraging the modular free energy system, we have implemented a set of multicomponent, multiphase models such as WBM and KKS that allow the user to combine arbitrary single phase free energies into multiphase free energies. Tight coupling to linear elasticity is enabled though free energy modules that provide the strain energy and its derivatives to the modular free energy system. The various free energy contributions are combined at runtime, alleviating the need to modify and recompile code. Together, these features offer increased flexibility when implementing multiphase models and solution methods, and when coupling to mechanics.
Acknowledgements
This work was funded by the Department of Energy Nuclear Energy Advanced Modeling and Simulation program and the Light Water Reactor Sustainability program. This manuscript has been authored by Battelle Energy Alliance, LLC under Contract No. DEAC0705ID14517 with the US Department of Energy. The publisher, by accepting the article for publication, acknowledges that the United States Government retains a nonexclusive, paidup, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.
References
Footnotes
 journal: Computational Materials Science
 https://github.com/idaholab/moose/tree/devel/modules/combined/examples/publications/rapid_dev
References
 L.Q. Chen, Phasefield models for microstructure evolution, Annual Review of Materials Research 32 (1) (2002) 113–140. doi:10.1146/annurev.matsci.32.112001.132041.
 N. Moelans, B. Blanpain, P. Wollants, An introduction to phasefield modeling of microstructure evolution, CALPHAD 32 (2) (2008) 268–294. doi:10.1016/j.calphad.2007.11.003.
 J. A. Warren, W. J. Boettinger, Prediction of dendritic growth and microsegregation patterns in a binary alloy using the phasefield method, Acta Metallurgica et Materialia 43 (2) (1995) 689–703. doi:10.1016/09567151(94)00285P.
 A. Karma, W.J. Rappel, Phasefield method for computationally efficient modeling of solidification with arbitrary interface kinetics, Physical Review E 53 (4) (1996) R3017 (4 pages). doi:10.1103/PhysRevE.53.R3017.
 A. A. Wheeler, W. J. Boettinger, G. B. McFadden, Phasefield model for isothermal phase transitions in binary alloys, Physical Review A 45 (10) (1992) 7424–7439. doi:10.1103/PhysRevA.45.7424.
 S. G. Kim, W. T. Kim, T. Suzuki, Phasefield model for binary alloys, Physical Review E 60 (6) (1999) 7186–7197. doi:10.1103/PhysRevE.60.7186.
 D. Fan, L.Q. Chen, Diffusioncontrolled grain growth in twophase solids, Acta Materialia 45 (8) (1997) 3297–3310. doi:10.1016/S13596454(97)000220.
 N. Moelans, B. Blanpain, P. Wollants, Quantitative analysis of grain boundary properties in a generalized phase field model for grain growth in anisotropic systems, Physical Review B 78 (2) (2008) 024113 (23 pages). doi:10.1103/PhysRevB.78.024113.
 A. Wheeler, B. Murray, R. Schaefer, Computation of dendrites using a phase field model, Physica D: Nonlinear Phenomena 66 (1) (1993) 243–262. doi:10.1016/01672789(93)90242S.
 J. E. Guyer, D. Wheeler, J. A. Warren, FiPy: Partial differential equations with Python, Computing in Science & Engineering 11 (3) (2009) 6–15. doi:10.1109/MCSE.2009.52.
 T. Takaki, T. Fukuoka, Y. Tomita, Phasefield simulation during directional solidification of a binary alloy using adaptive finite element method, Journal of Crystal Growth 283 (1) (2005) 263–278.
 M. R. Tonks, D. Gaston, P. C. Millett, D. Andrs, P. Talbot, An objectoriented finite element framework for multiphysics phase field simulations, Computational Materials Science 51 (1) (2012) 20–29. doi:10.1016/j.commatsci.2011.07.028.
 L.Q. Chen, J. Shen, Applications of semiimplicit fourierspectral method to phase field equations, Computer Physics Communications 108 (2) (1998) 147–158. doi:10.1016/S00104655(97)00115X.
 B. Puchala, G. Tarcea, E. A. Marquis, M. Hedstrom, H. V. Jagadish, J. E. Allison, The Materials Commons: A collaboration platform and information repository for the global materials community, JOM 68 (8) (2016) 2035–2044. doi:10.1007/s1183701619987.
 A. Logg, K.A. Mardal, G. Wells, Automated solution of differential equations by the finite element method: The FEniCS book, Vol. 84, Springer Science & Business Media, 2012.
 T. Keller, et al., The mesoscale microstructure simulation project, https://github.com/mesoscale/mmsp (2017).
 S. M. Allen, J. W. Cahn, Ground state structures in ordered binary alloys with second neighbor interactions, Acta Metallurgica 20 (3) (1972) 423–433. doi:10.1016/00016160(72)900375.
 J. W. Cahn, J. E. Hilliard, Free Energy of a Nonuniform System. I. Interfacial Free Energy, The Journal of Chemical Physics 28 (2) (1958) 258–267. doi:10.1063/1.1744102.
 D. R. Gaston, C. J. Permann, J. W. Peterson, A. E. Slaughter, D. Andrš, Y. Wang, M. P. Short, D. M. Perez, M. R. Tonks, J. Ortensi, R. C. Martineau, Physicsbased multiscale coupling for full core nuclear reactor simulation, Annals of Nuclear Energy, Special Issue on MultiPhysics Modelling of LWR Static and Transient Behaviour 84 (2015) 45–54. doi:10.1016/j.anucene.2014.09.060.
 P. N. Brown, A. C. Hindmarsh, Matrixfree methods for stiff systems of ODEs, SIAM Journal on Numerical Analysis 23 (3) (1986) 610–638. doi:10.1137/0723039.
 D. A. Knoll, D. E. Keyes, Jacobianfree NewtonKrylov methods: A survey of approaches and applications, Journal of Computational Physics 193 (2) (2004) 357–397. doi:10.1016/j.jcp.2003.08.010.
 B. S. Kirk, J. W. Peterson, R. H. Stogner, G. F. Carey, libMesh: A C++ Library for Parallel Adaptive Mesh Refinement/Coarsening Simulations, Engineering with Computers 22 (3–4) (2006) 237–254. doi:10.1007/s0036600600493.

S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman,
L. Dalcin, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C.
McInnes, K. Rupp, B. F. Smith, S. Zampini, H. Zhang, H. Zhang,
PETSc users manual, Tech. Rep.
ANL95/11  Revision 3.7, Argonne National Laboratory (2016).
URL http://www.mcs.anl.gov/petsc  L. Zhang, M. R. Tonks, D. Gaston, J. W. Peterson, D. Andrs, P. C. Millett, B. S. Biner, A quantitative comparison between and elements for solving the cahnhilliard equation, Journal of Computational Physics 236 (2013) 74–80. doi:10.1016/j.jcp.2012.12.001.

J. Nieminen, J. Yliluoma,
Function Parser Web
page (2011).
URL http://warp.povusers.org/FunctionParser  J. E. Tolsma, P. I. Barton, On computational differentiation, Computers & Chemical Engineering 22 (4–5) (1998) 475–490. doi:10.1016/S00981354(97)002640.
 G. Kedem, Automatic differentiation of computer programs, ACM Transactions on Mathematical Software (TOMS) 6 (2) (1980) 150–165. doi:10.1145/355887.355890.

The Open Group Base
Specifications Issue 6 – dlopen (2004).
URL http://www.opengroup.org/susv3xsh/dlopen.html  Q. H. Dang, Secure Hash Standard (SHS), Tech. Rep. FIPS1804, Information Technology Laboratory, National Institute of Standards and Technology, our implementation is from: http://www.tamale.net/sha1 (Mar. 2012). doi:10.6028/NIST.FIPS.1804.

C. Lattner, V. Adve,
LLVM: A
compilation framework for lifelong program analysis & transformation, in:
Proceedings of the International Symposium on Code Generation and
Optimization: Feedbackdirected and Runtime Optimization (CGO), 2004, pp.
75–87.
URL http://dl.acm.org/citation.cfm?id=977395.977673 
C. Lattner, et al., clang: A C language family
frontend for LLVM (2014).
URL http://clang.llvm.org/ 
R. Stallman, et al., GCC, the GNU compiler
collection (2014).
URL https://gcc.gnu.org/  J. H. Hildebrand, The term ‘regular solution’, Nature 168 (1951) 868. doi:10.1038/168868a0.
 A. Jokisaari, K. Thornton, General method for incorporating CALPHAD free energies of mixing into phase field models: Application to the zirconium/hydride system, CALPHAD 51 (2015) 334–343. doi:10.1016/j.calphad.2015.10.011.
 D. Schwen, E. Martinez, A. Caro, On the analytic calculation of critical size for alpha prime precipitation in FeCr, Journal of Nuclear Materials 439 (13) (2013) 180–184. doi:10.1016/j.jnucmat.2013.03.057.
 Y. Li, S. Hu, R. Montgomery, F. Gao, X. Sun, Phasefield simulations of intragranular fission gas bubble evolution in UO under postirradiation thermal annealing, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 303 (2013) 62–67. doi:10.1016/j.nimb.2012.11.028.
 M. Ohno, K. Matsuura, Quantitative phasefield modeling for twophase solidification process involving diffusion in the solid, Acta Materialia 58 (17) (2010) 5749–5758. doi:10.1016/j.actamat.2010.06.050.
 R. Folch, M. Plapp, Quantitative phasefield modeling of twophase growth, Physical Review E 72 (1) (2005) 011602 (27 pages). doi:10.1103/PhysRevE.72.011602.
 G. I. Tóth, T. Pusztai, L. Gránásy, Consistent multiphasefield theory for interface driven multidomain dynamics, Physical Review B 92 (18) (2015) 184105 (19 pages). doi:10.1103/PhysRevB.92.184105.
 S. G. Kim, W. T. Kim, T. Suzuki, M. Ode, Phasefield modeling of eutectic solidification, Journal of Crystal Growth 261 (1) (2004) 135–158. doi:10.1016/j.jcrysgro.2003.08.078.
 K. Ammar, B. Appolaire, G. Cailletaud, S. Forest, Combining phase field approach and homogenization methods for modelling phase transformation in elastoplastic media, European Journal of Computational Mechanics 18 (56) (2009) 485–523. doi:10.3166/ejcm.18.485523.
 A. G. Khachaturyan, Theory of Structural Transformations in Solids, Wiley, New York, 1983.
 V. Vaithyanathan, C. Wolverton, L.Q. Chen, Multiscale modeling of precipitate microstructure evolution, Physical Review Letters 88 (12) (2002) 125503 (4 pages). doi:10.1103/PhysRevLett.88.125503.
 V. Vaithyanathan, C. Wolverton, L.Q. Chen, Multiscale modeling of precipitation in AlCu binary alloys, Acta Materialia 52 (10) (2004) 2973–2987. doi:10.1016/j.actamat.2004.03.001.