ACMultiInterface

Gradient energy Allen-Cahn Kernel with cross terms

Implements Allen-Cahn interface terms for a multiphase system. This includes cross terms of the form

$var element = document.getElementById("moose-equation-6ee7acc6-b195-418c-bb89-26642fba1ea0");katex.render("\\sum_j\\frac12\\kappa_{ij}\\left( \\eta_i\\nabla\\eta_j - \\eta_j\\nabla\\eta_i \\right)^2", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

where $var element = document.getElementById("moose-equation-76aac4ea-7106-447e-aff6-0cd7aa390fae");katex.render("\\eta_i", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is the non-linear variable the kernel is acting on, $var element = document.getElementById("moose-equation-830e6ddc-1fbd-4743-a9cc-aa8d6b34ece2");katex.render("\\eta_j", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ (etas are all non-conserved order parameters in the system, $var element = document.getElementById("moose-equation-436aef61-8ea8-4482-b566-de3696806423");katex.render("\\kappa_{ij}", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ (kappa_name) are the gradient energy coefficents, and $var element = document.getElementById("moose-equation-5998b879-cc59-401c-b69f-83c1c709b226");katex.render("L_i", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ (mob_name) is the scalar (isotropic) mobility associated with the $var element = document.getElementById("moose-equation-891d3ea8-f65e-451e-af0f-b9be1906a706");katex.render("\\eta_i", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ order parameter.

Derivation

The interfacial free energy density $var element = document.getElementById("moose-equation-1900ee28-5343-4c19-b1b1-73b0d66281af");katex.render("f_{int}", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is implemented following Nestler and Wheeler (1998) equations (7) and (8) (also see footnote 1)

$var element = document.getElementById("moose-equation-b77d99e2-9dfb-415e-b9b9-ff8edcceb0d0");katex.render("f_{int} = \\sum_{\\substack{a,b \\\\ b\\neq a}} \\frac12 \\kappa_{ab} \\left| \\eta_a\\nabla\\eta_b - \\eta_b\\nabla\\eta_a\\right|^2,", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

Where the sum is taken over unique tuples _a,b_ (i.e. without the permutations _b,a_). We take the functional derivative taken using the lemma

$var element = document.getElementById("moose-equation-b879ed64-e245-4fcf-bd4c-d5945309443d");katex.render("\\frac{\\delta f}{\\delta\\eta} = \\frac{\\partial f}{\\partial\\eta} - \\nabla\\frac{\\partial f}{\\partial\\nabla\\eta}.", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

We obtain a one dimensional sum for each of the $var element = document.getElementById("moose-equation-46bf3f47-0135-4f99-9d06-6902634b7699");katex.render("\\eta", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ -derivatives.

$var element = document.getElementById("moose-equation-2207a195-56b8-4199-849e-23d2d987d0f2");katex.render("\\begin{aligned} \\frac{\\delta f_{int}}{\\delta\\eta_a} & = \\sum_b \\kappa_{ab} \\left[ (\\eta_a\\nabla\\eta_b - \\eta_b\\nabla\\eta_a)\\nabla\\eta_b + \\nabla\\left((\\eta_a\\nabla\\eta_b - \\eta_b\\nabla\\eta_a)\\eta_b\\right) \\right] \\\\ &= \\sum_{\\substack{b\\\\b\\neq a}} \\kappa_{ab} \\left[ \\underbrace{2(\\eta_a\\nabla\\eta_b - \\eta_b\\nabla\\eta_a)\\nabla\\eta_b}_{\\text{order 1}} + \\underbrace{\\eta_b(\\eta_a\\nabla^2\\eta_b - \\eta_b\\nabla^2\\eta_a)}_{\\text{order 2}} \\right] \\end{aligned}", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

We transform this expression into the weak form and see that the derivative order on the _order 2_ term has to be reduced by shifting a gradient onto the test function by applying the product rule

$var element = document.getElementById("moose-equation-8a487bf2-a6e7-40b1-9a9f-04deb6b492b1");katex.render("v \\nabla\\cdot\\mathbf{w} = -\\nabla v \\cdot \\mathbf{w} + \\nabla\\cdot (v\\mathbf{w}),", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

after multiplying with the test function $var element = document.getElementById("moose-equation-632572f1-dc3a-4494-b8ba-a6a9f9018692");katex.render("\\psi", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ and integrating over the volume $var element = document.getElementById("moose-equation-c5cf4f32-5370-4d5c-9e90-744774cfb784");katex.render("\\Omega", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ . We identify $var element = document.getElementById("moose-equation-8eecaa76-48df-4383-8c84-16bd665ce888");katex.render("v", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ and $var element = document.getElementById("moose-equation-2d4740ec-26a4-4759-aa29-73116b96c943");katex.render("w", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ as follows

$var element = document.getElementById("moose-equation-50012e75-1af5-4df8-b7a0-16e4696e51cd");katex.render("\\int_\\Omega\\psi\\eta_b(\\eta_a\\nabla^2\\eta_b - \\eta_b\\nabla^2\\eta_a) = \\int_\\Omega\\underbrace{\\psi\\eta_a\\eta_b}_{v_1}\\nabla\\cdot\\underbrace{\\nabla\\eta_b}_{\\mathbf{w}_1} - \\int_\\Omega\\underbrace{\\psi\\eta_b^2}_{v_2}\\nabla\\cdot\\underbrace{\\nabla\\eta_a}_{\\mathbf{w}_2}", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

$var element = document.getElementById("moose-equation-8c3becd7-ee7b-47ee-a79f-b9cf62c68af3");katex.render(" \\int_\\Omega\\left[ -\\nabla(\\psi\\eta_a\\eta_b)\\cdot\\nabla\\eta_b\\right] - \\int_\\omega\\left[-\\nabla(\\psi\\eta_b^2)\\cdot\\nabla\\eta_a\\right] + \\int_\\Omega \\nabla\\cdot(\\psi\\eta_a\\eta_b\\nabla\\eta_b) - \\int_\\Omega \\nabla\\cdot(\\psi\\eta_b^2\\nabla\\eta_a).", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

We get rid of the last two terms by applying the divergence theorem and obtain

$var element = document.getElementById("moose-equation-623deea8-f8f2-40a6-a035-1aec327e41ad");katex.render("\\underbrace{ \\int_\\Omega\\left[ -\\nabla(\\psi\\eta_a\\eta_b)\\cdot\\nabla\\eta_b\\right] - \\int_\\Omega\\left[-\\nabla(\\psi\\eta_b^2)\\cdot\\nabla\\eta_a\\right] }_{\\text{volume terms}} + \\underbrace{ \\int_{\\partial\\Omega} \\mathbf{n}\\cdot(\\psi\\eta_a\\eta_b\\nabla\\eta_b) - \\int_{\\partial\\Omega} \\mathbf{n}\\cdot(\\psi\\eta_b^2\\nabla\\eta_a). }_{\\text{boundary terms}}", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

to convert them from volume to surface/boundary integrals. We again apply the product rule to expand the gradient of the product in the _volume terms_ and obtain

$var element = document.getElementById("moose-equation-6225d5b7-c59e-4636-8f4d-50ed88941761");katex.render("\\int_\\Omega\\left[ -\\left( \\eta_a\\eta_b\\nabla\\psi + \\psi\\eta_b\\nabla\\eta_a + \\psi\\eta_a\\nabla\\eta_b \\right) \\cdot\\nabla\\eta_b\\right] - \\int_\\Omega\\left[ -\\left( \\eta_b^2\\nabla\\psi + 2\\psi\\eta_b\\nabla\\eta_b \\right)\\cdot\\nabla\\eta_a\\right].", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

Residual

The total residual $var element = document.getElementById("moose-equation-696d5b9a-cf54-48ae-9ccc-72a56f22b8ea");katex.render("R_a", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is then

$var element = document.getElementById("moose-equation-a573aec4-d421-4604-822b-48a2fa8877cc");katex.render("\\begin{aligned} R_a = L_a\\sum_{\\substack{b\\\\b\\neq a}}\\kappa_{ab} \\Bigg[ &\\,& \\int_\\Omega\\left[ 2\\psi(\\eta_a\\nabla\\eta_b - \\eta_b\\nabla\\eta_a)\\nabla\\eta_b \\right]& \\\\ &+& \\int_\\Omega\\left[ -\\left( \\eta_a\\eta_b\\nabla\\psi + \\psi\\eta_b\\nabla\\eta_a + \\psi\\eta_a\\nabla\\eta_b \\right) \\cdot\\nabla\\eta_b\\right] & \\\\ &-& \\int_\\Omega\\left[ -\\left( \\eta_b^2\\nabla\\psi + 2\\psi\\eta_b\\nabla\\eta_b \\right)\\cdot\\nabla\\eta_a\\right]& \\Bigg]. \\end{aligned}", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

On-diagonal Jacobian

The on-diagonal jacobian $var element = document.getElementById("moose-equation-435953a3-dcfa-4d0d-a049-bd2df364faff");katex.render("J_a", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is obtained by taking the derivative with respect to $var element = document.getElementById("moose-equation-fc395d22-3819-428f-b354-44c98d7ef5c4");katex.render("\\eta_{aj}", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ , where $var element = document.getElementById("moose-equation-c788f9d8-2e8f-4d15-b813-d74f83bd46c2");katex.render("\\frac{\\partial \\eta_a}{\\partial \\eta_{aj}} = \\phi_j", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ and $var element = document.getElementById("moose-equation-c0fa739a-e140-4c42-af1f-007f64e038c3");katex.render("\\frac{\\partial \\nabla\\eta_{aj}}{\\partial \\eta_{aj}} = \\nabla\\phi_j", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

$var element = document.getElementById("moose-equation-270969ca-9d6a-4e0e-8d68-f538d35c6de2");katex.render("\\begin{aligned} J_a = L_a\\sum_{\\substack{b\\\\b\\neq a}}\\kappa_{ab} \\Bigg[ &\\,&\\int_\\Omega2\\psi\\left[ (\\phi_j\\nabla\\eta_b - \\eta_b\\nabla\\phi_j)\\nabla\\eta_b \\right]& \\\\ &+& \\int_\\Omega\\left[ -\\left( \\phi_j\\eta_b\\nabla\\psi + \\psi\\eta_b\\nabla\\phi_j + \\psi\\phi_j\\nabla\\eta_b \\right)\\cdot\\nabla\\eta_b \\right]& \\\\ &-& \\int_\\Omega\\left[ -\\left( \\eta_b^2\\nabla\\psi + 2\\psi\\eta_b\\nabla\\eta_b \\right)\\cdot\\nabla\\phi_j \\right]& \\Bigg]. \\end{aligned}", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

Off-diagonal jacobian

For the off diagonal Jacobian entry $var element = document.getElementById("moose-equation-3062fe62-7bb4-410e-9f90-a1dbe812b69d");katex.render("J_{ab}", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ we take the derivative $var element = document.getElementById("moose-equation-50464386-077d-4fde-9307-2015008e24f5");katex.render("\\frac{\\partial R_a}{\\partial \\eta_{bj}}", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ and obtain

\begin{equation} \begin{aligned} J_{ab} = &\,& L_a\kappa_{ab}\int_\Omega2\psi\left[ (\eta_a\nabla\phi_j - \phi_j\nabla\eta_a)\nabla\eta_b + (\eta_a\nabla\eta_b - \eta_b\nabla\eta_a)\nabla\phi_j \right]
&+& \int_\Omega\left[ -\left( \eta_a\phi_j\nabla\psi + \psi\phi_j\nabla\eta_a + \psi\eta_a\nabla\phi_j \right) \cdot\nabla\eta_b -\left( \eta_a\eta_b\nabla\psi + \psi\eta_b\nabla\eta_a + \psi\eta_a\nabla\eta_b \right) \cdot\nabla\phi_j \right]
&-& \int_\Omega\left[ %-\left( \eta_b^2\nabla\psi + 2\psi\eta_b\nabla\eta_b \right)\cdot\nabla\eta_a -\left( 2\eta_b\phi_j\nabla\psi + 2\psi(\phi_j\nabla\eta_b + \eta_b\nabla\phi_j) \right)\cdot\nabla\eta_a \right] \end{aligned}\end{equation}

----

1) Note, that in the two-phase case with $var element = document.getElementById("moose-equation-1e88cb28-1d54-4d52-a1af-0706c7e2bbc6");katex.render("\\eta_b=1-\\eta_a", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ this reduces to

$var element = document.getElementById("moose-equation-1dafaf23-ec51-4d7b-a2cb-4ec61094a431");katex.render("\\begin{aligned} f_{int} & = \\frac12 \\kappa_{ab} \\left| \\eta_a\\nabla(1-\\eta_a) - (1-\\eta_a)\\nabla\\eta_a\\right|^2 \\\\ & = \\frac12 \\kappa_{ab} \\left| -\\eta_a\\nabla\\eta_a - \\nabla\\eta_a + \\eta_a\\nabla\\eta_a\\right|^2 \\\\ & = \\frac12 \\kappa_{ab} \\left| \\nabla\\eta_a \\right|^2, \\end{aligned}", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

which is the familiar form implemented by ACInterface.

Input Parameters

etasAll eta_i order parameters of the multiphase problem
C++ Type:std::vector<VariableName>
Options:
Description:All eta_i order parameters of the multiphase problem
kappa_namesThe kappa used with the kernel
C++ Type:std::vector<MaterialPropertyName>
Options:
Description:The kappa used with the kernel
variableThe name of the variable that this residual object operates on
C++ Type:NonlinearVariableName
Options:
Description:The name of the variable that this residual object operates on

Required Parameters

blockThe list of block ids (SubdomainID) that this object will be applied
C++ Type:std::vector<SubdomainName>
Options:
Description:The list of block ids (SubdomainID) that this object will be applied
displacementsThe displacements
C++ Type:std::vector<VariableName>
Options:
Description:The displacements
mob_nameLThe mobility used with the kernel
Default:L
C++ Type:MaterialPropertyName
Options:
Description:The mobility used with the kernel

Optional Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Options:
Description:Adds user-defined labels for accessing object parameters via control logic.
diag_save_inThe name of auxiliary variables to save this Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Options:
Description:The name of auxiliary variables to save this Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Options:
Description:Set the enabled status of the MooseObject.
implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Options:
Description:Determines whether this object is calculated using an implicit or explicit form
save_inThe name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Options:
Description:The name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Options:
Description:The seed for the master random number generator
use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Options:
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters

extra_matrix_tagsThe extra tags for the matrices this Kernel should fill
C++ Type:std::vector<TagName>
Options:
Description:The extra tags for the matrices this Kernel should fill
extra_vector_tagsThe extra tags for the vectors this Kernel should fill
C++ Type:std::vector<TagName>
Options:
Description:The extra tags for the vectors this Kernel should fill
matrix_tagssystemThe tag for the matrices this Kernel should fill
Default:system
C++ Type:MultiMooseEnum
Options:nontime, system
Description:The tag for the matrices this Kernel should fill
vector_tagsnontimeThe tag for the vectors this Kernel should fill
Default:nontime
C++ Type:MultiMooseEnum
Options:nontime, time
Description:The tag for the vectors this Kernel should fill

Tagging Parameters

Input Files

(modules/combined/examples/phase_field-mechanics/Pattern1.i)
(modules/phase_field/test/tests/MultiPhase/acmultiinterface.i)
(modules/phase_field/examples/multiphase/DerivativeMultiPhaseMaterial.i)
(modules/combined/examples/publications/rapid_dev/fig8.i)
(modules/phase_field/test/tests/MultiPhase/acmultiinterface_aux.i)

References

B. Nestler and A. A. Wheeler. Anisotropic multi-phase-field model: interfaces and junctions. Phys. Rev. E, 57:2602–2609, Mar 1998. URL: https://link.aps.org/doi/10.1103/PhysRevE.57.2602, doi:10.1103/PhysRevE.57.2602.

@article{Nest98,
    author = "Nestler, B. and Wheeler, A. A.",
    title = "Anisotropic multi-phase-field model: Interfaces and junctions",
    journal = "Phys. Rev. E",
    volume = "57",
    issue = "3",
    pages = "2602--2609",
    numpages = "0",
    year = "1998",
    month = "Mar",
    publisher = "American Physical Society",
    doi = "10.1103/PhysRevE.57.2602",
    url = "https://link.aps.org/doi/10.1103/PhysRevE.57.2602"
}

(modules/combined/examples/phase_field-mechanics/Pattern1.i)

#
# Pattern example 1
#
# Phase changes driven by a combination mechanical (elastic) and chemical
# driving forces. In this three phase system a matrix phase, an oversized and
# an undersized precipitate phase compete. The chemical free energy favors a
# phase separation into either precipitate phase. A mix of both precipitate
# emerges to balance lattice expansion and contraction.
#
# This example demonstrates the use of
# * ACMultiInterface
# * SwitchingFunctionConstraintEta and SwitchingFunctionConstraintLagrange
# * DerivativeParsedMaterial
# * ElasticEnergyMaterial
# * DerivativeMultiPhaseMaterial
# * MultiPhaseStressMaterial
# which are the components to se up a phase field model with an arbitrary number
# of phases
#

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 80
  ny = 80
  nz = 0
  xmin = -20
  xmax = 20
  ymin = -20
  ymax = 20
  zmin = 0
  zmax = 0
  elem_type = QUAD4
[]

[GlobalParams]
  # CahnHilliard needs the third derivatives
  derivative_order = 3
  enable_jit = true
  displacements = 'disp_x disp_y'
[]

# AuxVars to compute the free energy density for outputting
[AuxVariables]
  [./local_energy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./cross_energy]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[AuxKernels]
  [./local_free_energy]
    type = TotalFreeEnergy
    variable = local_energy
    interfacial_vars = 'c'
    kappa_names = 'kappa_c'
    additional_free_energy = cross_energy
  [../]
  [./cross_terms]
    type = CrossTermGradientFreeEnergy
    variable = cross_energy
    interfacial_vars = 'eta1 eta2 eta3'
    kappa_names = 'kappa11 kappa12 kappa13
                   kappa21 kappa22 kappa23
                   kappa31 kappa32 kappa33'
  [../]
[]

[Variables]
  # Solute concentration variable
  [./c]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = RandomIC
      min = 0
      max = 0.8
      seed = 1235
    [../]
  [../]

  # Order parameter for the Matrix
  [./eta1]
    order = FIRST
    family = LAGRANGE
    initial_condition = 0.5
  [../]
  # Order parameters for the 2 different inclusion orientations
  [./eta2]
    order = FIRST
    family = LAGRANGE
    initial_condition = 0.1
  [../]
  [./eta3]
    order = FIRST
    family = LAGRANGE
    initial_condition = 0.1
  [../]

  # Mesh displacement
  [./disp_x]
    order = FIRST
    family = LAGRANGE
  [../]
  [./disp_y]
    order = FIRST
    family = LAGRANGE
  [../]

  # Lagrange-multiplier
  [./lambda]
    order = FIRST
    family = LAGRANGE
    initial_condition = 1.0
  [../]
[]

[Kernels]
  # Set up stress divergence kernels
  [./TensorMechanics]
  [../]

  # Cahn-Hilliard kernels
  [./c_res]
    type = CahnHilliard
    variable = c
    f_name = F
    args = 'eta1 eta2 eta3'
  [../]
  [./time]
    type = TimeDerivative
    variable = c
  [../]

  # Allen-Cahn and Lagrange-multiplier constraint kernels for order parameter 1
  [./deta1dt]
    type = TimeDerivative
    variable = eta1
  [../]
  [./ACBulk1]
    type = AllenCahn
    variable = eta1
    args = 'eta2 eta3 c'
    mob_name = L1
    f_name = F
  [../]
  [./ACInterface1]
    type = ACMultiInterface
    variable = eta1
    etas = 'eta1 eta2 eta3'
    mob_name = L1
    kappa_names = 'kappa11 kappa12 kappa13'
  [../]
  [./lagrange1]
    type = SwitchingFunctionConstraintEta
    variable = eta1
    h_name   = h1
    lambda = lambda
  [../]

  # Allen-Cahn and Lagrange-multiplier constraint kernels for order parameter 2
  [./deta2dt]
    type = TimeDerivative
    variable = eta2
  [../]
  [./ACBulk2]
    type = AllenCahn
    variable = eta2
    args = 'eta1 eta3 c'
    mob_name = L2
    f_name = F
  [../]
  [./ACInterface2]
    type = ACMultiInterface
    variable = eta2
    etas = 'eta1 eta2 eta3'
    mob_name = L2
    kappa_names = 'kappa21 kappa22 kappa23'
  [../]
  [./lagrange2]
    type = SwitchingFunctionConstraintEta
    variable = eta2
    h_name   = h2
    lambda = lambda
  [../]

  # Allen-Cahn and Lagrange-multiplier constraint kernels for order parameter 3
  [./deta3dt]
    type = TimeDerivative
    variable = eta3
  [../]
  [./ACBulk3]
    type = AllenCahn
    variable = eta3
    args = 'eta1 eta2 c'
    mob_name = L3
    f_name = F
  [../]
  [./ACInterface3]
    type = ACMultiInterface
    variable = eta3
    etas = 'eta1 eta2 eta3'
    mob_name = L3
    kappa_names = 'kappa31 kappa32 kappa33'
  [../]
  [./lagrange3]
    type = SwitchingFunctionConstraintEta
    variable = eta3
    h_name   = h3
    lambda = lambda
  [../]

  # Lagrange-multiplier constraint kernel for lambda
  [./lagrange]
    type = SwitchingFunctionConstraintLagrange
    variable = lambda
    etas    = 'eta1 eta2 eta3'
    h_names = 'h1   h2   h3'
    epsilon = 1e-6
  [../]
[]

[Materials]
  # declare a few constants, such as mobilities (L,M) and interface gradient prefactors (kappa*)
  [./consts]
    type = GenericConstantMaterial
    prop_names  = 'M   kappa_c  L1 L2 L3  kappa11 kappa12 kappa13 kappa21 kappa22 kappa23 kappa31 kappa32 kappa33'
    prop_values = '0.2 0        1  1  1   2.00    2.00    2.00    2.00    2.00    2.00    2.00    2.00    2.00   '
  [../]

  # We use this to output the level of constraint enforcement
  # ideally it should be 0 everywhere, if the constraint is fully enforced
  [./etasummat]
    type = ParsedMaterial
    f_name = etasum
    args = 'eta1 eta2 eta3'
    material_property_names = 'h1 h2 h3'
    function = 'h1+h2+h3-1'
    outputs = exodus
  [../]

  # This parsed material creates a single property for visualization purposes.
  # It will be 0 for phase 1, -1 for phase 2, and 1 for phase 3
  [./phasemap]
    type = ParsedMaterial
    f_name = phase
    args = 'eta2 eta3'
    function = 'if(eta3>0.5,1,0)-if(eta2>0.5,1,0)'
    outputs = exodus
  [../]

  # matrix phase
  [./elasticity_tensor_1]
    type = ComputeElasticityTensor
    base_name = phase1
    C_ijkl = '3 3'
    fill_method = symmetric_isotropic
  [../]
  [./strain_1]
    type = ComputeSmallStrain
    base_name = phase1
    displacements = 'disp_x disp_y'
  [../]
  [./stress_1]
    type = ComputeLinearElasticStress
    base_name = phase1
  [../]

  # oversized phase
  [./elasticity_tensor_2]
    type = ComputeElasticityTensor
    base_name = phase2
    C_ijkl = '7 7'
    fill_method = symmetric_isotropic
  [../]
  [./strain_2]
    type = ComputeSmallStrain
    base_name = phase2
    displacements = 'disp_x disp_y'
    eigenstrain_names = eigenstrain
  [../]
  [./stress_2]
    type = ComputeLinearElasticStress
    base_name = phase2
  [../]
  [./eigenstrain_2]
    type = ComputeEigenstrain
    base_name = phase2
    eigen_base = '0.02'
    eigenstrain_name = eigenstrain
  [../]

  # undersized phase
  [./elasticity_tensor_3]
    type = ComputeElasticityTensor
    base_name = phase3
    C_ijkl = '7 7'
    fill_method = symmetric_isotropic
  [../]
  [./strain_3]
    type = ComputeSmallStrain
    base_name = phase3
    displacements = 'disp_x disp_y'
    eigenstrain_names = eigenstrain
  [../]
  [./stress_3]
    type = ComputeLinearElasticStress
    base_name = phase3
  [../]
  [./eigenstrain_3]
    type = ComputeEigenstrain
    base_name = phase3
    eigen_base = '-0.05'
    eigenstrain_name = eigenstrain
  [../]

  # switching functions
  [./switching1]
    type = SwitchingFunctionMaterial
    function_name = h1
    eta = eta1
    h_order = SIMPLE
  [../]
  [./switching2]
    type = SwitchingFunctionMaterial
    function_name = h2
    eta = eta2
    h_order = SIMPLE
  [../]
  [./switching3]
    type = SwitchingFunctionMaterial
    function_name = h3
    eta = eta3
    h_order = SIMPLE
  [../]

  [./barrier]
    type = MultiBarrierFunctionMaterial
    etas = 'eta1 eta2 eta3'
  [../]

  # chemical free energies
  [./chemical_free_energy_1]
    type = DerivativeParsedMaterial
    f_name = Fc1
    function = '4*c^2'
    args = 'c'
    derivative_order = 2
  [../]
  [./chemical_free_energy_2]
    type = DerivativeParsedMaterial
    f_name = Fc2
    function = '(c-0.9)^2-0.4'
    args = 'c'
    derivative_order = 2
  [../]
  [./chemical_free_energy_3]
    type = DerivativeParsedMaterial
    f_name = Fc3
    function = '(c-0.9)^2-0.5'
    args = 'c'
    derivative_order = 2
  [../]

  # elastic free energies
  [./elastic_free_energy_1]
    type = ElasticEnergyMaterial
    base_name = phase1
    f_name = Fe1
    derivative_order = 2
    args = 'c' # should be empty
  [../]
  [./elastic_free_energy_2]
    type = ElasticEnergyMaterial
    base_name = phase2
    f_name = Fe2
    derivative_order = 2
    args = 'c' # should be empty
  [../]
  [./elastic_free_energy_3]
    type = ElasticEnergyMaterial
    base_name = phase3
    f_name = Fe3
    derivative_order = 2
    args = 'c' # should be empty
  [../]

  # phase free energies (chemical + elastic)
  [./phase_free_energy_1]
    type = DerivativeSumMaterial
    f_name = F1
    sum_materials = 'Fc1 Fe1'
    args = 'c'
    derivative_order = 2
  [../]
  [./phase_free_energy_2]
    type = DerivativeSumMaterial
    f_name = F2
    sum_materials = 'Fc2 Fe2'
    args = 'c'
    derivative_order = 2
  [../]
  [./phase_free_energy_3]
    type = DerivativeSumMaterial
    f_name = F3
    sum_materials = 'Fc3 Fe3'
    args = 'c'
    derivative_order = 2
  [../]

  # global free energy
  [./free_energy]
    type = DerivativeMultiPhaseMaterial
    f_name = F
    fi_names = 'F1  F2  F3'
    hi_names = 'h1  h2  h3'
    etas     = 'eta1 eta2 eta3'
    args = 'c'
    W = 3
  [../]

  # Generate the global stress from the phase stresses
  [./global_stress]
    type = MultiPhaseStressMaterial
    phase_base = 'phase1 phase2 phase3'
    h          = 'h1     h2     h3'
  [../]
[]

[BCs]
  # the boundary conditions on the displacement enforce periodicity
  # at zero total shear and constant volume
  [./bottom_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom'
    value = 0
  [../]
  [./top_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'top'
    value = 0
  [../]
  [./left_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0
  [../]
  [./right_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'right'
    value = 0
  [../]

  [./Periodic]
    [./disp_x]
      auto_direction = 'y'
    [../]
    [./disp_y]
      auto_direction = 'x'
    [../]

    # all other phase field variables are fully periodic
    [./c]
      auto_direction = 'x y'
    [../]
    [./eta1]
      auto_direction = 'x y'
    [../]
    [./eta2]
      auto_direction = 'x y'
    [../]
    [./eta3]
      auto_direction = 'x y'
    [../]
    [./lambda]
      auto_direction = 'x y'
    [../]
  [../]
[]

[Preconditioning]
  [./SMP]
    type = SMP
    full = true
  [../]
[]

# We monitor the total free energy and the total solute concentration (should be constant)
[Postprocessors]
  [./total_free_energy]
    type = ElementIntegralVariablePostprocessor
    variable = local_energy
  [../]
  [./total_solute]
    type = ElementIntegralVariablePostprocessor
    variable = c
  [../]
[]

[Executioner]
  type = Transient
  scheme = bdf2

  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -sub_pc_type'
  petsc_options_value = 'asm      ilu'

  l_max_its = 30
  nl_max_its = 10
  l_tol = 1.0e-4
  nl_rel_tol = 1.0e-8
  nl_abs_tol = 1.0e-10
  start_time = 0.0
  num_steps = 200

  [./TimeStepper]
    type = SolutionTimeAdaptiveDT
    dt = 0.1
  [../]
[]

[Outputs]
  execute_on = 'timestep_end'
  exodus = true
  [./table]
    type = CSV
    delimiter = ' '
  [../]
[]

[Debug]
  # show_var_residual_norms = true
[]

(modules/phase_field/test/tests/MultiPhase/acmultiinterface.i)

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 20
  ny = 10
  nz = 0
  xmin = -10
  xmax = 10
  ymin = -5
  ymax = 5
  elem_type = QUAD4
[]

[Variables]
  [./eta1]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = SmoothCircleIC
      x1 = -3.5
      y1 =  0.0
      radius = 4.0
      invalue = 0.9
      outvalue = 0.1
      int_width = 2.0
    [../]
  [../]
  [./eta2]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = SmoothCircleIC
      x1 =  3.5
      y1 =  0.0
      radius = 4.0
      invalue = 0.9
      outvalue = 0.1
      int_width = 2.0
    [../]
  [../]
  [./eta3]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = SpecifiedSmoothCircleIC
      x_positions = '-4.0 4.0'
      y_positions = ' 0.0 0.0'
      z_positions = ' 0.0 0.0'
      radii = '4.0 4.0'
      invalue = 0.1
      outvalue = 0.9
      int_width = 2.0
    [../]
  [../]

  [./lambda]
    order = FIRST
    family = LAGRANGE
    initial_condition = 1.0
  [../]
[]

[Kernels]
  [./deta1dt]
    type = TimeDerivative
    variable = eta1
  [../]
  [./ACBulk1]
    type = AllenCahn
    variable = eta1
    args = 'eta2 eta3'
    mob_name = L1
    f_name = F
  [../]
  [./ACInterface1]
    type = ACMultiInterface
    variable = eta1
    etas = 'eta1 eta2 eta3'
    mob_name = L1
    kappa_names = 'kappa11 kappa12 kappa13'
  [../]
  [./lagrange1]
    type = SwitchingFunctionConstraintEta
    variable = eta1
    h_name   = h1
    lambda = lambda
  [../]

  [./deta2dt]
    type = TimeDerivative
    variable = eta2
  [../]
  [./ACBulk2]
    type = AllenCahn
    variable = eta2
    args = 'eta1 eta3'
    mob_name = L2
    f_name = F
  [../]
  [./ACInterface2]
    type = ACMultiInterface
    variable = eta2
    etas = 'eta1 eta2 eta3'
    mob_name = L2
    kappa_names = 'kappa21 kappa22 kappa23'
  [../]
  [./lagrange2]
    type = SwitchingFunctionConstraintEta
    variable = eta2
    h_name   = h2
    lambda = lambda
  [../]

  [./deta3dt]
    type = TimeDerivative
    variable = eta3
  [../]
  [./ACBulk3]
    type = AllenCahn
    variable = eta3
    args = 'eta1 eta2'
    mob_name = L3
    f_name = F
  [../]
  [./ACInterface3]
    type = ACMultiInterface
    variable = eta3
    etas = 'eta1 eta2 eta3'
    mob_name = L3
    kappa_names = 'kappa31 kappa32 kappa33'
  [../]
  [./lagrange3]
    type = SwitchingFunctionConstraintEta
    variable = eta3
    h_name   = h3
    lambda = lambda
  [../]

  [./lagrange]
    type = SwitchingFunctionConstraintLagrange
    variable = lambda
    etas    = 'eta1 eta2 eta3'
    h_names = 'h1   h2   h3'
    epsilon = 0
  [../]
[]

[BCs]
  [./Periodic]
    [./All]
      auto_direction = 'x y'
    [../]
  [../]
[]

[Materials]
  [./consts]
    type = GenericConstantMaterial
    prop_names  = 'Fx  L1 L2 L3  kappa11 kappa12 kappa13 kappa21 kappa22 kappa23 kappa31 kappa32 kappa33'
    prop_values = '0   1  1  1   1       1       1       1       1       1       1       1       1      '
  [../]

  [./etasummat]
    type = ParsedMaterial
    f_name = etasum
    args = 'eta1 eta2 eta3'
    material_property_names = 'h1 h2 h3'
    function = 'h1+h2+h3'
  [../]

  [./switching1]
    type = SwitchingFunctionMaterial
    function_name = h1
    eta = eta1
    h_order = SIMPLE
  [../]
  [./switching2]
    type = SwitchingFunctionMaterial
    function_name = h2
    eta = eta2
    h_order = SIMPLE
  [../]
  [./switching3]
    type = SwitchingFunctionMaterial
    function_name = h3
    eta = eta3
    h_order = SIMPLE
  [../]

  [./barrier]
    type = MultiBarrierFunctionMaterial
    etas = 'eta1 eta2 eta3'
  [../]

  [./free_energy]
    type = DerivativeMultiPhaseMaterial
    f_name = F
    # we use a constant free energy (GeneriConstantmaterial property Fx)
    fi_names = 'Fx  Fx  Fx'
    hi_names = 'h1  h2  h3'
    etas     = 'eta1 eta2 eta3'
    # the free energy is given by the MultiBarrierFunctionMaterial only
    W = 1
    derivative_order = 2
  [../]
[]

[Preconditioning]
  [./SMP]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient
  scheme = 'bdf2'

  solve_type = 'PJFNK'
  #petsc_options = '-snes_ksp -snes_ksp_ew'
  #petsc_options = '-ksp_monitor_snes_lg-snes_ksp_ew'
  #petsc_options_iname = '-ksp_gmres_restart'
  #petsc_options_value = '1000              '

  l_max_its = 15
  l_tol = 1.0e-6

  nl_max_its = 50
  nl_rel_tol = 1.0e-8
  nl_abs_tol = 1.0e-10

  start_time = 0.0
  num_steps = 2
  dt = 0.2
[]

[Outputs]
  execute_on = 'timestep_end'
  exodus = true
[]

(modules/phase_field/examples/multiphase/DerivativeMultiPhaseMaterial.i)

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 40
  ny = 40
  nz = 0
  xmin = -12
  xmax = 12
  ymin = -12
  ymax = 12
  elem_type = QUAD4
[]

[GlobalParams]
  # let's output all material properties for demonstration purposes
  outputs = exodus

  # prefactor on the penalty function kernels. The higher this value is, the
  # more rigorously the constraint is enforced
  penalty = 1e3
[]

#
# These AuxVariables hold the directly calculated free energy density in the
# simulation cell. They are provided for visualization purposes.
#
[AuxVariables]
  [./local_energy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./cross_energy]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[AuxKernels]
  [./local_free_energy]
    type = TotalFreeEnergy
    variable = local_energy
    interfacial_vars = 'c'
    kappa_names = 'kappa_c'
    additional_free_energy = cross_energy
  [../]

  #
  # Helper kernel to cpompute the gradient contribution from interfaces of order
  # parameters evolved using the ACMultiInterface kernel
  #
  [./cross_terms]
    type = CrossTermGradientFreeEnergy
    variable = cross_energy
    interfacial_vars = 'eta1 eta2 eta3'
    #
    # The interface coefficient matrix. This should be symmetrical!
    #
    kappa_names = 'kappa11 kappa12 kappa13
                   kappa21 kappa22 kappa23
                   kappa31 kappa32 kappa33'
  [../]
[]

[Variables]
  [./c]
    order = FIRST
    family = LAGRANGE
    #
    # We set up a smooth cradial concentrtaion gradient
    # The concentration will quickly change to adapt to the preset order
    # parameters eta1, eta2, and eta3
    #
    [./InitialCondition]
      type = SmoothCircleIC
      x1 = 0.0
      y1 = 0.0
      radius = 5.0
      invalue = 1.0
      outvalue = 0.01
      int_width = 10.0
    [../]
  [../]

  [./eta1]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = FunctionIC
      #
      # Note: this initial conditions sets up a _sharp_ interface. Ideally
      # we should start with a smooth interface with a width consistent
      # with the kappa parameter supplied for the given interface.
      #
      function = 'r:=sqrt(x^2+y^2);if(r<=4,1,0)'
    [../]
  [../]
  [./eta2]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = FunctionIC
      function = 'r:=sqrt(x^2+y^2);if(r>4&r<=7,1,0)'
    [../]
  [../]
  [./eta3]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = FunctionIC
      function = 'r:=sqrt(x^2+y^2);if(r>7,1,0)'
    [../]
  [../]
[]

[Kernels]
  #
  # Cahn-Hilliard kernel for the concentration variable.
  # Note that we are not using an interfcae kernel on this variable, but rather
  # rely on the interface width enforced on the order parameters. This allows us
  # to use a direct solve using the CahnHilliard kernel _despite_ only using first
  # order elements.
  #
  [./c_res]
    type = CahnHilliard
    variable = c
    f_name = F
    args = 'eta1 eta2 eta3'
  [../]
  [./time]
    type = TimeDerivative
    variable = c
  [../]

  #
  # Order parameter eta1
  # Each order parameter is acted on by 4 kernels:
  #  1. The stock time derivative deta_i/dt kernel
  #  2. The Allen-Cahn kernel that takes a Dervative Material for the free energy
  #  3. A gradient interface kernel that includes cross terms
  #     see http://mooseframework.org/wiki/PhysicsModules/PhaseField/DevelopingModels/MultiPhaseModels/ACMultiInterface/
  #  4. A penalty contribution that forces the interface contributions h(eta)
  #     to sum up to unity
  #
  [./deta1dt]
    type = TimeDerivative
    variable = eta1
  [../]
  [./ACBulk1]
    type = AllenCahn
    variable = eta1
    args = 'eta2 eta3 c'
    mob_name = L1
    f_name = F
  [../]
  [./ACInterface1]
    type = ACMultiInterface
    variable = eta1
    etas = 'eta1 eta2 eta3'
    mob_name = L1
    kappa_names = 'kappa11 kappa12 kappa13'
  [../]
  [./penalty1]
    type = SwitchingFunctionPenalty
    variable = eta1
    etas    = 'eta1 eta2 eta3'
    h_names = 'h1   h2   h3'
  [../]

  #
  # Order parameter eta2
  #
  [./deta2dt]
    type = TimeDerivative
    variable = eta2
  [../]
  [./ACBulk2]
    type = AllenCahn
    variable = eta2
    args = 'eta1 eta3 c'
    mob_name = L2
    f_name = F
  [../]
  [./ACInterface2]
    type = ACMultiInterface
    variable = eta2
    etas = 'eta1 eta2 eta3'
    mob_name = L2
    kappa_names = 'kappa21 kappa22 kappa23'
  [../]
  [./penalty2]
    type = SwitchingFunctionPenalty
    variable = eta2
    etas    = 'eta1 eta2 eta3'
    h_names = 'h1   h2   h3'
  [../]

  #
  # Order parameter eta3
  #
  [./deta3dt]
    type = TimeDerivative
    variable = eta3
  [../]
  [./ACBulk3]
    type = AllenCahn
    variable = eta3
    args = 'eta1 eta2 c'
    mob_name = L3
    f_name = F
  [../]
  [./ACInterface3]
    type = ACMultiInterface
    variable = eta3
    etas = 'eta1 eta2 eta3'
    mob_name = L3
    kappa_names = 'kappa31 kappa32 kappa33'
  [../]
  [./penalty3]
    type = SwitchingFunctionPenalty
    variable = eta3
    etas    = 'eta1 eta2 eta3'
    h_names = 'h1   h2   h3'
  [../]
[]

[BCs]
  [./Periodic]
    [./All]
      auto_direction = 'x y'
    [../]
  [../]
[]

[Materials]
  # here we declare some of the model parameters: the mobilities and interface
  # gradient prefactors. For this example we use arbitrary numbers. In an actual simulation
  # physical mobilities would be used, and the interface gradient prefactors would
  # be readjusted to the free energy magnitudes.
  [./consts]
    type = GenericConstantMaterial
    prop_names  = 'M   kappa_c L1 L2 L3  kappa11 kappa12 kappa13 kappa21 kappa22 kappa23 kappa31 kappa32 kappa33'
    prop_values = '0.2 0.75    1  1  1   0.75    0.75    0.75    0.75    0.75    0.75    0.75    0.75    0.75   '
  [../]

  # This material sums up the individual phase contributions. It is written to the output file
  # (see GlobalParams section above) and can be used to check the constraint enforcement.
  [./etasummat]
    type = ParsedMaterial
    f_name = etasum
    args = 'eta1 eta2 eta3'
    material_property_names = 'h1 h2 h3'
    function = 'h1+h2+h3'
  [../]

  # The phase contribution factors for each material point are computed using the
  # SwitchingFunctionMaterials. Each phase with an order parameter eta contributes h(eta)
  # to the global free energy density. h is a function that switches smoothly from 0 to 1
  [./switching1]
    type = SwitchingFunctionMaterial
    function_name = h1
    eta = eta1
    h_order = SIMPLE
  [../]
  [./switching2]
    type = SwitchingFunctionMaterial
    function_name = h2
    eta = eta2
    h_order = SIMPLE
  [../]
  [./switching3]
    type = SwitchingFunctionMaterial
    function_name = h3
    eta = eta3
    h_order = SIMPLE
  [../]

  # The barrier function adds a phase transformation energy barrier. It also
  # Drives order parameters toward the [0:1] interval to avoid negative or larger than 1
  # order parameters (these are set to 0 and 1 contribution by the switching functions
  # above)
  [./barrier]
    type = MultiBarrierFunctionMaterial
    etas = 'eta1 eta2 eta3'
  [../]

  # We use DerivativeParsedMaterials to specify three (very) simple free energy
  # expressions for the three phases. All necessary derivatives are built automatically.
  # In a real problem these expressions can be arbitrarily complex (or even provided
  # by custom kernels).
  [./phase_free_energy_1]
    type = DerivativeParsedMaterial
    f_name = F1
    function = '(c-1)^2'
    args = 'c'
  [../]
  [./phase_free_energy_2]
    type = DerivativeParsedMaterial
    f_name = F2
    function = '(c-0.5)^2'
    args = 'c'
  [../]
  [./phase_free_energy_3]
    type = DerivativeParsedMaterial
    f_name = F3
    function = 'c^2'
    args = 'c'
  [../]

  # The DerivativeMultiPhaseMaterial ties the phase free energies together into a global free energy.
  # http://mooseframework.org/wiki/PhysicsModules/PhaseField/DevelopingModels/MultiPhaseModels/
  [./free_energy]
    type = DerivativeMultiPhaseMaterial
    f_name = F
    # we use a constant free energy (GeneriConstantmaterial property Fx)
    fi_names = 'F1  F2  F3'
    hi_names = 'h1  h2  h3'
    etas     = 'eta1 eta2 eta3'
    args = 'c'
    W = 1
  [../]
[]

[Postprocessors]
  # The total free energy of the simulation cell to observe the energy reduction.
  [./total_free_energy]
    type = ElementIntegralVariablePostprocessor
    variable = local_energy
  [../]

  # for testing we also monitor the total solute amount, which should be conserved.
  [./total_solute]
    type = ElementIntegralVariablePostprocessor
    variable = c
  [../]
[]

[Preconditioning]
  # This preconditioner makes sure the Jacobian Matrix is fully populated. Our
  # kernels compute all Jacobian matrix entries.
  # This allows us to use the Newton solver below.
  [./SMP]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient
  scheme = 'bdf2'

  # Automatic differentiation provedes a _full_ Jacobian in this example
  # so we can safely use NEWTON for a fast solve
  solve_type = 'NEWTON'

  l_max_its = 15
  l_tol = 1.0e-6

  nl_max_its = 50
  nl_rel_tol = 1.0e-6
  nl_abs_tol = 1.0e-6

  start_time = 0.0
  end_time   = 150.0

  [./TimeStepper]
    type = SolutionTimeAdaptiveDT
    dt = 0.1
  [../]
[]

[Debug]
  # show_var_residual_norms = true
[]

[Outputs]
  execute_on = 'timestep_end'
  exodus = true
  [./table]
    type = CSV
    delimiter = ' '
  [../]
[]

(modules/combined/examples/publications/rapid_dev/fig8.i)

#
# Fig. 8 input for 10.1016/j.commatsci.2017.02.017
# D. Schwen et al./Computational Materials Science 132 (2017) 36-45
# Two growing particles with differnet anisotropic Eigenstrains
#

[Mesh]
  [./gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 80
    ny = 40
    xmin = -20
    xmax = 20
    ymin = 0
    ymax = 20
    elem_type = QUAD4
  [../]
  [./cnode]
    type = ExtraNodesetGenerator
    input = gen
    coord = '0.0 0.0'
    new_boundary = 100
    tolerance = 0.1
  [../]
[]

[GlobalParams]
  # CahnHilliard needs the third derivatives
  derivative_order = 3
  enable_jit = true
  displacements = 'disp_x disp_y'
  int_width = 1
[]

# AuxVars to compute the free energy density for outputting
[AuxVariables]
  [./local_energy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./cross_energy]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[AuxKernels]
  [./local_free_energy]
    type = TotalFreeEnergy
    variable = local_energy
    interfacial_vars = 'c'
    kappa_names = 'kappa_c'
    additional_free_energy = cross_energy
    execute_on = 'INITIAL TIMESTEP_END'
  [../]
  [./cross_terms]
    type = CrossTermGradientFreeEnergy
    variable = cross_energy
    interfacial_vars = 'eta1 eta2 eta3'
    kappa_names = 'kappa11 kappa12 kappa13
                   kappa21 kappa22 kappa23
                   kappa31 kappa32 kappa33'
    execute_on = 'INITIAL TIMESTEP_END'
  [../]
[]

# particle x positions and radius
P1X=8
P2X=-4
PR=2

[Variables]
  # Solute concentration variable
  [./c]
    [./InitialCondition]
      type = SpecifiedSmoothCircleIC
      x_positions = '${P1X} ${P2X}'
      y_positions = '0 0'
      z_positions = '0 0'
      radii = '${PR} ${PR}'
      outvalue = 0.5
      invalue = 0.9
    [../]
  [../]
  [./w]
  [../]

  # Order parameter for the Matrix
  [./eta1]
    [./InitialCondition]
      type = SpecifiedSmoothCircleIC
      x_positions = '${P1X} ${P2X}'
      y_positions = '0 0'
      z_positions = '0 0'
      radii = '${PR} ${PR}'
      outvalue = 1.0
      invalue = 0.0
    [../]
  [../]
  # Order parameters for the 2 different inclusion orientations
  [./eta2]
    [./InitialCondition]
      type = SmoothCircleIC
      x1 = ${P2X}
      y1 = 0
      radius = ${PR}
      invalue = 1.0
      outvalue = 0.0
    [../]
  [../]
  [./eta3]
    [./InitialCondition]
      type = SmoothCircleIC
      x1 = ${P1X}
      y1 = 0
      radius = ${PR}
      invalue = 1.0
      outvalue = 0.0
    [../]
  [../]

  # Lagrange-multiplier
  [./lambda]
    initial_condition = 1.0
  [../]
[]

[Modules]
  [./TensorMechanics]
    [./Master]
      [./all]
        add_variables = true
        strain = SMALL
        eigenstrain_names = eigenstrain
      [../]
    [../]
  [../]
[]

[Kernels]
  # Split Cahn-Hilliard kernels
  [./c_res]
    type = SplitCHParsed
    variable = c
    f_name = F
    args = 'eta1 eta2 eta3'
    kappa_name = kappa_c
    w = w
  [../]
  [./wres]
    type = SplitCHWRes
    variable = w
    mob_name = M
  [../]
  [./time]
    type = CoupledTimeDerivative
    variable = w
    v = c
  [../]

  # Allen-Cahn and Lagrange-multiplier constraint kernels for order parameter 1
  [./deta1dt]
    type = TimeDerivative
    variable = eta1
  [../]
  [./ACBulk1]
    type = AllenCahn
    variable = eta1
    args = 'eta2 eta3 c'
    mob_name = L1
    f_name = F
  [../]
  [./ACInterface1]
    type = ACMultiInterface
    variable = eta1
    etas = 'eta1 eta2 eta3'
    mob_name = L1
    kappa_names = 'kappa11 kappa12 kappa13'
  [../]
  [./lagrange1]
    type = SwitchingFunctionConstraintEta
    variable = eta1
    h_name   = h1
    lambda = lambda
  [../]

  # Allen-Cahn and Lagrange-multiplier constraint kernels for order parameter 2
  [./deta2dt]
    type = TimeDerivative
    variable = eta2
  [../]
  [./ACBulk2]
    type = AllenCahn
    variable = eta2
    args = 'eta1 eta3 c'
    mob_name = L2
    f_name = F
  [../]
  [./ACInterface2]
    type = ACMultiInterface
    variable = eta2
    etas = 'eta1 eta2 eta3'
    mob_name = L2
    kappa_names = 'kappa21 kappa22 kappa23'
  [../]
  [./lagrange2]
    type = SwitchingFunctionConstraintEta
    variable = eta2
    h_name   = h2
    lambda = lambda
  [../]

  # Allen-Cahn and Lagrange-multiplier constraint kernels for order parameter 3
  [./deta3dt]
    type = TimeDerivative
    variable = eta3
  [../]
  [./ACBulk3]
    type = AllenCahn
    variable = eta3
    args = 'eta1 eta2 c'
    mob_name = L3
    f_name = F
  [../]
  [./ACInterface3]
    type = ACMultiInterface
    variable = eta3
    etas = 'eta1 eta2 eta3'
    mob_name = L3
    kappa_names = 'kappa31 kappa32 kappa33'
  [../]
  [./lagrange3]
    type = SwitchingFunctionConstraintEta
    variable = eta3
    h_name   = h3
    lambda = lambda
  [../]

  # Lagrange-multiplier constraint kernel for lambda
  [./lagrange]
    type = SwitchingFunctionConstraintLagrange
    variable = lambda
    etas    = 'eta1 eta2 eta3'
    h_names = 'h1   h2   h3'
    epsilon = 1e-6
  [../]
[]

[Materials]
  # declare a few constants, such as mobilities (L,M) and interface gradient prefactors (kappa*)
  [./consts]
    type = GenericConstantMaterial
    block = 0
    prop_names  = 'M   kappa_c  L1 L2 L3  kappa11 kappa12 kappa13 kappa21 kappa22 kappa23 kappa31 kappa32 kappa33'
    prop_values = '0.2 0.5      1  1  1   2.00    2.00    2.00    2.00    2.00    2.00    2.00    2.00    2.00   '
  [../]

  # We use this to output the level of constraint enforcement
  # ideally it should be 0 everywhere, if the constraint is fully enforced
  [./etasummat]
    type = ParsedMaterial
    f_name = etasum
    args = 'eta1 eta2 eta3'
    material_property_names = 'h1 h2 h3'
    function = 'h1+h2+h3-1'
    outputs = exodus
  [../]

  # This parsed material creates a single property for visualization purposes.
  # It will be 0 for phase 1, -1 for phase 2, and 1 for phase 3
  [./phasemap]
    type = ParsedMaterial
    f_name = phase
    args = 'eta2 eta3'
    function = 'if(eta3>0.5,1,0)-if(eta2>0.5,1,0)'
    outputs = exodus
  [../]

  # global mechanical properties
  [./elasticity_tensor]
    type = ComputeElasticityTensor
    C_ijkl = '400 400'
    fill_method = symmetric_isotropic
  [../]
  [./stress]
    type = ComputeLinearElasticStress
  [../]

  # eigenstrain
  [./eigenstrain_2]
    type = GenericConstantRankTwoTensor
    tensor_name = s2
    tensor_values = '0 -0.05 0  0 0 0'
  [../]
  [./eigenstrain_3]
    type = GenericConstantRankTwoTensor
    tensor_name = s3
    tensor_values =  '-0.05 0 0  0 0 0'
  [../]
  [./eigenstrain]
    type = CompositeEigenstrain
    weights = 'h2 h3'
    tensors = 's2 s3'
    args = 'eta2 eta3'
    eigenstrain_name = eigenstrain
  [../]

  # switching functions
  [./switching1]
    type = SwitchingFunctionMaterial
    function_name = h1
    eta = eta1
    h_order = SIMPLE
  [../]
  [./switching2]
    type = SwitchingFunctionMaterial
    function_name = h2
    eta = eta2
    h_order = SIMPLE
  [../]
  [./switching3]
    type = SwitchingFunctionMaterial
    function_name = h3
    eta = eta3
    h_order = SIMPLE
  [../]

  [./barrier]
    type = MultiBarrierFunctionMaterial
    etas = 'eta1 eta2 eta3'
  [../]

  # chemical free energies
  [./chemical_free_energy_1]
    type = DerivativeParsedMaterial
    f_name = Fc1
    function = '4*c^2'
    args = 'c'
    derivative_order = 2
  [../]
  [./chemical_free_energy_2]
    type = DerivativeParsedMaterial
    f_name = Fc2
    function = '(c-0.9)^2-0.4'
    args = 'c'
    derivative_order = 2
  [../]
  [./chemical_free_energy_3]
    type = DerivativeParsedMaterial
    f_name = Fc3
    function = '(c-0.9)^2-0.5'
    args = 'c'
    derivative_order = 2
  [../]

  # global chemical free energy
  [./chemical_free_energy]
    type = DerivativeMultiPhaseMaterial
    f_name = Fc
    fi_names = 'Fc1  Fc2  Fc3'
    hi_names = 'h1  h2  h3'
    etas     = 'eta1 eta2 eta3'
    args = 'c'
    W = 3
  [../]

  # global elastic free energy
  [./elastic_free_energy]
    type = ElasticEnergyMaterial
    f_name = Fe
    args = 'eta2 eta3'
    outputs = exodus
    output_properties = Fe
    derivative_order = 2
  [../]

  # Penalize phase 2 and 3 coexistence
  [./multi_phase_penalty]
    type = DerivativeParsedMaterial
    f_name = Fp
    function = '50*(eta2*eta3)^2'
    args = 'eta2 eta3'
    derivative_order = 2
    outputs = exodus
    output_properties = Fp
  [../]

  # free energy
  [./free_energy]
    type = DerivativeSumMaterial
    f_name = F
    sum_materials = 'Fc Fe Fp'
    args = 'c eta1 eta2 eta3'
    derivative_order = 2
  [../]
[]

[BCs]
  # fix center point location
  [./centerfix_x]
    type = DirichletBC
    boundary = 100
    variable = disp_x
    value = 0
  [../]

  # fix side point x coordinate to inhibit rotation
  [./angularfix]
    type = DirichletBC
    boundary = bottom
    variable = disp_y
    value = 0
  [../]
[]

[Preconditioning]
  [./SMP]
    type = SMP
    full = true
  [../]
[]

# We monitor the total free energy and the total solute concentration (should be constant)
[Postprocessors]
  [./total_free_energy]
    type = ElementIntegralVariablePostprocessor
    variable = local_energy
    execute_on = 'INITIAL TIMESTEP_END'
  [../]
  [./total_solute]
    type = ElementIntegralVariablePostprocessor
    variable = c
    execute_on = 'INITIAL TIMESTEP_END'
  [../]
[]

[Executioner]
  type = Transient
  scheme = bdf2

  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -sub_pc_type'
  petsc_options_value = 'asm      lu'

  l_max_its = 30
  nl_max_its = 10
  l_tol = 1.0e-4
  nl_rel_tol = 1.0e-8
  nl_abs_tol = 1.0e-10
  start_time = 0.0
  end_time = 12.0

  [./TimeStepper]
    type = IterationAdaptiveDT
    optimal_iterations = 8
    iteration_window = 1
    dt = 0.01
  [../]
[]

[Outputs]
  print_linear_residuals = false
  execute_on = 'INITIAL TIMESTEP_END'
  exodus = true
  [./table]
    type = CSV
    delimiter = ' '
  [../]
[]

[Debug]
  # show_var_residual_norms = true
[]

(modules/phase_field/test/tests/MultiPhase/acmultiinterface_aux.i)

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 20
  ny = 10
  nz = 0
  xmin = -10
  xmax = 10
  ymin = -5
  ymax = 5
  elem_type = QUAD4
[]

[AuxVariables]
  [./eta1]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = SmoothCircleIC
      x1 = -3.5
      y1 =  0.0
      radius = 4.0
      invalue = 0.9
      outvalue = 0.1
      int_width = 2.0
    [../]
  [../]
[]

[Variables]
  [./eta2]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = SmoothCircleIC
      x1 =  3.5
      y1 =  0.0
      radius = 4.0
      invalue = 0.9
      outvalue = 0.1
      int_width = 2.0
    [../]
  [../]
  [./eta3]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = SpecifiedSmoothCircleIC
      x_positions = '-4.0 4.0'
      y_positions = ' 0.0 0.0'
      z_positions = ' 0.0 0.0'
      radii = '4.0 4.0'
      invalue = 0.1
      outvalue = 0.9
      int_width = 2.0
    [../]
  [../]

  [./lambda]
    order = FIRST
    family = LAGRANGE
    initial_condition = 1.0
  [../]
[]

[Kernels]
  [./deta2dt]
    type = TimeDerivative
    variable = eta2
  [../]
  [./ACBulk2]
    type = AllenCahn
    variable = eta2
    args = 'eta1 eta3'
    mob_name = L2
    f_name = F
  [../]
  [./ACInterface2]
    type = ACMultiInterface
    variable = eta2
    etas = 'eta1 eta2 eta3'
    mob_name = L2
    kappa_names = 'kappa21 kappa22 kappa23'
  [../]
  [./lagrange2]
    type = SwitchingFunctionConstraintEta
    variable = eta2
    h_name   = h2
    lambda = lambda
  [../]

  [./deta3dt]
    type = TimeDerivative
    variable = eta3
  [../]
  [./ACBulk3]
    type = AllenCahn
    variable = eta3
    args = 'eta1 eta2'
    mob_name = L3
    f_name = F
  [../]
  [./ACInterface3]
    type = ACMultiInterface
    variable = eta3
    etas = 'eta1 eta2 eta3'
    mob_name = L3
    kappa_names = 'kappa31 kappa32 kappa33'
  [../]
  [./lagrange3]
    type = SwitchingFunctionConstraintEta
    variable = eta3
    h_name   = h3
    lambda = lambda
  [../]

  [./lagrange]
    type = SwitchingFunctionConstraintLagrange
    variable = lambda
    etas    = 'eta1 eta2 eta3'
    h_names = 'h1   h2   h3'
    epsilon = 0
  [../]
[]

[BCs]
  [./Periodic]
    [./All]
      auto_direction = 'x y'
    [../]
  [../]
[]

[Materials]
  [./consts]
    type = GenericConstantMaterial
    prop_names  = 'Fx  L1 L2 L3  kappa11 kappa12 kappa13 kappa21 kappa22 kappa23 kappa31 kappa32 kappa33'
    prop_values = '0   1  1  1   1       1       1       1       1       1       1       1       1      '
  [../]

  [./switching1]
    type = SwitchingFunctionMaterial
    function_name = h1
    eta = eta1
    h_order = SIMPLE
  [../]
  [./switching2]
    type = SwitchingFunctionMaterial
    function_name = h2
    eta = eta2
    h_order = SIMPLE
  [../]
  [./switching3]
    type = SwitchingFunctionMaterial
    function_name = h3
    eta = eta3
    h_order = SIMPLE
  [../]

  [./barrier]
    type = MultiBarrierFunctionMaterial
    etas = 'eta1 eta2 eta3'
  [../]

  [./free_energy]
    type = DerivativeMultiPhaseMaterial
    f_name = F
    # we use a constant free energy (GeneriConstantmaterial property Fx)
    fi_names = 'Fx  Fx  Fx'
    hi_names = 'h1  h2  h3'
    etas     = 'eta1 eta2 eta3'
    # the free energy is given by the MultiBarrierFunctionMaterial only
    W = 1
    derivative_order = 2
  [../]
[]

[Preconditioning]
  [./SMP]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient
  scheme = 'bdf2'

  solve_type = 'PJFNK'
  #petsc_options = '-snes_ksp -snes_ksp_ew'
  #petsc_options = '-ksp_monitor_snes_lg-snes_ksp_ew'
  #petsc_options_iname = '-ksp_gmres_restart'
  #petsc_options_value = '1000              '

  l_max_its = 15
  l_tol = 1.0e-6

  nl_max_its = 50
  nl_rel_tol = 1.0e-8
  nl_abs_tol = 1.0e-10

  start_time = 0.0
  num_steps = 2
  dt = 0.2
[]

[Outputs]
  execute_on = 'timestep_end'
  exodus = true
[]

Derivation
Input Parameters
Input Files
References

Install MOOSE

New Users

Examples and Tutorials

Application Usage

Physics and Syntax

Application Development

Framework Development

MOOSEDocs

Infrastructure

Questions

Information and Tools

INL Applications and Remote Access