DerivativeTwoPhaseMaterial

Two phase material that combines two single phase materials using a switching function.

The simplified two-phase model** uses a single order parameter to switch between the two phases. A global free energy is constructed using a meta material class that combines the phase free energies.

For two phase models the DerivativeTwoPhaseMaterial can be used to combine two phase free energies into a global free energy (which the AllenCahn and Cahn-Hilliard kernels use to evolve the system) as

$var element = document.getElementById("moose-equation-9a733fc0-54f3-4b43-b9da-2a084ab09e51");katex.render("F = \\left(1-h(\\eta)\\right) F_a + h(\\eta)F_b + Wg(\\eta)", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

Input Parameters

etaOrder parameter
C++ Type:std::vector<VariableName>
Options:
Description:Order parameter
fa_namePhase A material (at eta=0)
C++ Type:MaterialPropertyName
Options:
Description:Phase A material (at eta=0)
fb_namePhase A material (at eta=1)
C++ Type:MaterialPropertyName
Options:
Description:Phase A material (at eta=1)

Required Parameters

W0Energy barrier for the phase transformation from A to B
Default:0
C++ Type:double
Options:
Description:Energy barrier for the phase transformation from A to B
argsArguments of fa and fb - use vector coupling
C++ Type:std::vector<VariableName>
Options:
Description:Arguments of fa and fb - use vector coupling
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
boundaryThe list of boundary IDs from the mesh where this boundary condition applies
C++ Type:std::vector<BoundaryName>
Options:
Description:The list of boundary IDs from the mesh where this boundary condition applies
computeTrueWhen false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
Default:True
C++ Type:bool
Options:
Description:When false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
constant_onNONEWhen ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
Default:NONE
C++ Type:MooseEnum
Options:NONE, ELEMENT, SUBDOMAIN
Description:When ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
derivative_order3Maximum order of derivatives taken (2 or 3)
Default:3
C++ Type:unsigned int
Options:
Description:Maximum order of derivatives taken (2 or 3)
displacement_gradientsVector of displacement gradient variables (see Modules/PhaseField/DisplacementGradients action)
C++ Type:std::vector<VariableName>
Options:
Description:Vector of displacement gradient variables (see Modules/PhaseField/DisplacementGradients action)
f_nameFBase name of the free energy function (used to name the material properties)
Default:F
C++ Type:std::string
Options:
Description:Base name of the free energy function (used to name the material properties)
ggBarrier Function Material that provides g(eta)
Default:g
C++ Type:MaterialPropertyName
Options:
Description:Barrier Function Material that provides g(eta)
hhSwitching Function Material that provides h(eta)
Default:h
C++ Type:MaterialPropertyName
Options:
Description:Switching Function Material that provides h(eta)

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.
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
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

output_propertiesList of material properties, from this material, to output (outputs must also be defined to an output type)
C++ Type:std::vector<std::string>
Options:
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
outputsnone Vector of output names were you would like to restrict the output of variables(s) associated with this object
Default:none
C++ Type:std::vector<OutputName>
Options:
Description:Vector of output names were you would like to restrict the output of variables(s) associated with this object

Outputs Parameters

Input Files

(modules/phase_field/test/tests/MultiPhase/mixedswitchingfunctionmaterial.i)
(modules/combined/examples/phase_field-mechanics/Nonconserved.i)
(modules/combined/test/tests/surface_tension_KKS/surface_tension_KKS.i)
(modules/combined/examples/publications/rapid_dev/fig7a.i)
(modules/combined/examples/phase_field-mechanics/kks_mechanics_VTS.i)
(modules/phase_field/test/tests/MultiPhase/derivativetwophasematerial.i)
(modules/phase_field/test/tests/GrandPotentialPFM/GrandPotentialPFM.i)
(modules/combined/examples/publications/rapid_dev/fig7b.i)

Note that the phase free energies usually have single well character. The global free energy landscape will however have a double well character in the examples above.

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

# This is a test of the MixedSwitchingfunctionmaterial
# Several mixed type of switching function with ajustable weight parameter

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

[Variables]
  [./eta]
  [../]
[]

[ICs]
  [./IC_eta]
    type = SmoothCircleIC
    variable = eta
    x1 = 10
    y1 = 10
    radius = 5
    invalue = 1
    outvalue = 0
    int_width = 1
  [../]
[]

[Kernels]
  [./eta_bulk]
    type = AllenCahn
    variable = eta
    f_name = F
  [../]
  [./eta_interface]
    type = ACInterface
    variable = eta
    kappa_name = kappa_eta
  [../]

  [./detadt]
    type = TimeDerivative
    variable = eta
  [../]
[]

[Materials]
  [./consts]
    type = GenericConstantMaterial
    prop_names  = 'L kappa_eta'
    prop_values = '1.0 1.0'
  [../]

  [./switching]
    type = MixedSwitchingFunctionMaterial
    function_name = h
    eta = eta
    h_order = MIX234
    weight = 1.0
  [../]

  [./barrier]
    type = BarrierFunctionMaterial
    eta = eta
    g_order = SIMPLE
  [../]

# Total free energy: F = Fa*(1-h) + Fb*h
  [./free_energy]
    type = DerivativeTwoPhaseMaterial
    f_name = F
    fa_name = '0'
    fb_name = '-1'
    eta = eta
    W = 3.1
    derivative_order = 2
    outputs = exodus
  [../]
[]

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

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = 'PJFNK'

  l_max_its = 30
  nl_max_its = 10
  l_tol = 1.0e-4
  nl_rel_tol = 1.0e-10
  nl_abs_tol = 1.0e-12
  start_time = 0.0
  num_steps = 2

  [./TimeStepper]
    type = IterationAdaptiveDT
    optimal_iterations = 9
    iteration_window = 2
    growth_factor = 1.1
    cutback_factor = 0.75
    dt = 0.3
  [../]
[]

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

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

#
# Example 2
# Phase change driven by a mechanical (elastic) driving force.
# An oversized phase inclusion grows under a uniaxial tensile stress.
# Check the file below for comments and suggestions for parameter modifications.
#

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

[Variables]
  [./eta]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = SmoothCircleIC
      x1 = 0
      y1 = 0
      radius = 30.0
      invalue = 1.0
      outvalue = 0.0
      int_width = 10.0
    [../]
  [../]
  [./disp_x]
    order = FIRST
    family = LAGRANGE
  [../]
  [./disp_y]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[Kernels]
  [./TensorMechanics]
    displacements = 'disp_x disp_y'
  [../]

  [./eta_bulk]
    type = AllenCahn
    variable = eta
    f_name = F
  [../]
  [./eta_interface]
    type = ACInterface
    variable = eta
    kappa_name = 1
  [../]
  [./time]
    type = TimeDerivative
    variable = eta
  [../]
[]

#
# Try visualizing the stress tensor components as done in Conserved.i
#

[Materials]
  [./consts]
    type = GenericConstantMaterial
    block = 0
    prop_names  = 'L'
    prop_values = '1'
  [../]

  # matrix phase
  [./stiffness_a]
    type = ComputeElasticityTensor
    base_name = phasea
    block = 0
    # lambda, mu values
    C_ijkl = '7 7'
    # Stiffness tensor is created from lambda=7, mu=7 for symmetric_isotropic fill method
    fill_method = symmetric_isotropic
    # See RankFourTensor.h for details on fill methods
  [../]
  [./strain_a]
    type = ComputeSmallStrain
    block = 0
    displacements = 'disp_x disp_y'
    base_name = phasea
  [../]
  [./stress_a]
    type = ComputeLinearElasticStress
    block = 0
    base_name = phasea
  [../]
  [./elastic_free_energy_a]
    type = ElasticEnergyMaterial
    base_name = phasea
    f_name = Fea
    block = 0
    args = ''
  [../]

  # oversized precipitate phase (simulated using thermal expansion)
  [./stiffness_b]
    type = ComputeElasticityTensor
    base_name = phaseb
    block = 0
    # Stiffness tensor lambda, mu values
    # Note that the two phases could have different stiffnesses.
    # Try reducing the precipitate stiffness (to '1 1') rather than making it oversized
    C_ijkl = '7 7'
    fill_method = symmetric_isotropic
  [../]
  [./strain_b]
    type = ComputeSmallStrain
    block = 0
    displacements = 'disp_x disp_y'
    base_name = phaseb
    eigenstrain_names = eigenstrain
  [../]
  [./eigenstrain_b]
    type = ComputeEigenstrain
    base_name = phaseb
    eigen_base = '0.1 0.1 0.1'
    eigenstrain_name = eigenstrain
  [../]
  [./stress_b]
    type = ComputeLinearElasticStress
    block = 0
    base_name = phaseb
  [../]
  [./elastic_free_energy_b]
    type = ElasticEnergyMaterial
    base_name = phaseb
    f_name = Feb
    block = 0
    args = ''
  [../]

  # Generate the global free energy from the phase free energies
  [./switching]
    type = SwitchingFunctionMaterial
    block = 0
    eta = eta
    h_order = SIMPLE
  [../]
  [./barrier]
    type = BarrierFunctionMaterial
    block = 0
    eta = eta
    g_order = SIMPLE
  [../]
  [./free_energy]
    type = DerivativeTwoPhaseMaterial
    block = 0
    f_name = F
    fa_name = Fea
    fb_name = Feb
    eta = eta
    args = ''
    W = 0.1
    derivative_order = 2
  [../]

  # Generate the global stress from the phase stresses
  [./global_stress]
    type = TwoPhaseStressMaterial
    block = 0
    base_A = phasea
    base_B = phaseb
  [../]
[]

[BCs]
  [./bottom_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom'
    value = 0
  [../]
  [./top_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'top'
    value = 5
  [../]
  [./left_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0
  [../]
[]

[Preconditioning]
  # active = ' '
  [./SMP]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient
  scheme = bdf2

  # this gives best performance on 4 cores
  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
  num_steps = 200

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

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

(modules/combined/test/tests/surface_tension_KKS/surface_tension_KKS.i)

#
# KKS coupled with elasticity. Physical parameters for matrix and precipitate phases
# are gamma and gamma-prime phases, respectively, in the Ni-Al system.
# Parameterization is as described in L.K. Aagesen et al., Computational Materials
# Science, 140, 10-21 (2017), with isotropic elastic properties in both phases
# and without eigenstrain.
#

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 200
  xmax = 200
[]

[Problem]
  coord_type = RSPHERICAL
[]

[GlobalParams]
  displacements = 'disp_x'
[]

[Variables]
  # order parameter
  [./eta]
    order = FIRST
    family = LAGRANGE
  [../]

  # solute concentration
  [./c]
    order = FIRST
    family = LAGRANGE
  [../]

  # chemical potential
  [./w]
    order = FIRST
    family = LAGRANGE
  [../]

  # solute phase concentration (matrix)
  [./cm]
    order = FIRST
    family = LAGRANGE
    initial_condition = 0.13
  [../]
  # solute phase concentration (precipitate)
  [./cp]
    order = FIRST
    family = LAGRANGE
    initial_condition = 0.235
  [../]
[]

[AuxVariables]
  [./energy_density]
    family = MONOMIAL
  [../]
  [./extra_xx]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./extra_yy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./extra_zz]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./strain_xx]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./strain_yy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./strain_zz]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[ICs]
  [./eta_ic]
    variable = eta
    type = FunctionIC
    function = ic_func_eta
  [../]
  [./c_ic]
    variable = c
    type = FunctionIC
    function = ic_func_c
  [../]
[]

[Functions]
  [./ic_func_eta]
    type = ParsedFunction
    value = 'r:=sqrt(x^2+y^2+z^2);0.5*(1.0-tanh((r-r0)/delta_eta/sqrt(2.0)))'
    vars = 'delta_eta r0'
    vals = '6.431     100'
  [../]
  [./ic_func_c]
    type = ParsedFunction
    value = 'r:=sqrt(x^2+y^2+z^2);eta_an:=0.5*(1.0-tanh((r-r0)/delta/sqrt(2.0)));0.235*eta_an^3*(6*eta_an^2-15*eta_an+10)+0.13*(1-eta_an^3*(6*eta_an^2-15*eta_an+10))'
    vars = 'delta r0'
    vals = '6.431 100'
  [../]
[]

[Modules/TensorMechanics/Master]
  [./all]
    add_variables = true
    generate_output = 'hydrostatic_stress stress_xx stress_yy stress_zz'
  [../]
[]


[Kernels]

  # enforce c = (1-h(eta))*cm + h(eta)*cp
  [./PhaseConc]
    type = KKSPhaseConcentration
    ca       = cm
    variable = cp
    c        = c
    eta      = eta
  [../]

  # enforce pointwise equality of chemical potentials
  [./ChemPotVacancies]
    type = KKSPhaseChemicalPotential
    variable = cm
    cb       = cp
    fa_name  = f_total_matrix
    fb_name  = f_total_ppt
  [../]

  #
  # Cahn-Hilliard Equation
  #
  [./CHBulk]
    type = KKSSplitCHCRes
    variable = c
    ca       = cm
    fa_name  = f_total_matrix
    w        = w
  [../]

  [./dcdt]
    type = CoupledTimeDerivative
    variable = w
    v = c
  [../]
  [./ckernel]
    type = SplitCHWRes
    mob_name = M
    variable = w
  [../]

  #
  # Allen-Cahn Equation
  #
  [./ACBulkF]
    type = KKSACBulkF
    variable = eta
    fa_name  = f_total_matrix
    fb_name  = f_total_ppt
    w        = 0.0033
    args = 'cp cm'
  [../]
  [./ACBulkC]
    type = KKSACBulkC
    variable = eta
    ca       = cm
    cb       = cp
    fa_name  = f_total_matrix
  [../]
  [./ACInterface]
    type = ACInterface
    variable = eta
    kappa_name = kappa
  [../]
  [./detadt]
    type = TimeDerivative
    variable = eta
  [../]
[]

[AuxKernels]
  [./extra_xx]
    type = RankTwoAux
    rank_two_tensor = extra_stress
    index_i = 0
    index_j = 0
    variable = extra_xx
  [../]
  [./extra_yy]
    type = RankTwoAux
    rank_two_tensor = extra_stress
    index_i = 1
    index_j = 1
    variable = extra_yy
  [../]
  [./extra_zz]
    type = RankTwoAux
    rank_two_tensor = extra_stress
    index_i = 2
    index_j = 2
    variable = extra_zz
  [../]
  [./strain_xx]
    type = RankTwoAux
    rank_two_tensor = mechanical_strain
    index_i = 0
    index_j = 0
    variable = strain_xx
  [../]
  [./strain_yy]
    type = RankTwoAux
    rank_two_tensor = mechanical_strain
    index_i = 1
    index_j = 1
    variable = strain_yy
  [../]
  [./strain_zz]
    type = RankTwoAux
    rank_two_tensor = mechanical_strain
    index_i = 2
    index_j = 2
    variable = strain_zz
  [../]
[]

[Materials]
  # Chemical free energy of the matrix
  [./fm]
    type = DerivativeParsedMaterial
    f_name = fm
    args = 'cm'
    function = '6.55*(cm-0.13)^2'
  [../]
# Elastic energy of the matrix
  [./elastic_free_energy_m]
    type = ElasticEnergyMaterial
    base_name = matrix
    f_name = fe_m
    args = ' '
  [../]
# Total free energy of the matrix
  [./Total_energy_matrix]
    type = DerivativeSumMaterial
    f_name = f_total_matrix
    sum_materials = 'fm fe_m'
    args = 'cm'
  [../]

  # Free energy of the precipitate phase
  [./fp]
    type = DerivativeParsedMaterial
    f_name = fp
    args = 'cp'
    function = '6.55*(cp-0.235)^2'
  [../]

# Elastic energy of the precipitate
  [./elastic_free_energy_p]
    type = ElasticEnergyMaterial
    base_name = ppt
    f_name = fe_p
    args = ' '
  [../]

# Total free energy of the precipitate
  [./Total_energy_ppt]
    type = DerivativeSumMaterial
    f_name = f_total_ppt
    sum_materials = 'fp fe_p'
    args = 'cp'
  [../]

# Total elastic energy
  [./Total_elastic_energy]
    type = DerivativeTwoPhaseMaterial
    eta = eta
    f_name = f_el_mat
    fa_name = fe_m
    fb_name = fe_p
    outputs = exodus
    W = 0
  [../]

  # h(eta)
  [./h_eta]
    type = SwitchingFunctionMaterial
    h_order = HIGH
    eta = eta
  [../]

  # g(eta)
  [./g_eta]
    type = BarrierFunctionMaterial
    g_order = SIMPLE
    eta = eta
    outputs = exodus
  [../]

  # constant properties
  [./constants]
    type = GenericConstantMaterial
    prop_names  = 'M   L   kappa'
    prop_values = '0.7 0.7 0.1365'
  [../]

  #Mechanical properties
  [./Stiffness_matrix]
    type = ComputeElasticityTensor
    C_ijkl = '74.25 14.525'
    base_name = matrix
    fill_method = symmetric_isotropic
  [../]
  [./Stiffness_ppt]
    type = ComputeElasticityTensor
    C_ijkl = '74.25 14.525'
    base_name = ppt
    fill_method = symmetric_isotropic
  [../]
  [./strain_matrix]
    type = ComputeRSphericalSmallStrain
    base_name = matrix
  [../]
  [./strain_ppt]
    type = ComputeRSphericalSmallStrain
    base_name = ppt
  [../]
  [./stress_matrix]
    type = ComputeLinearElasticStress
    base_name = matrix
  [../]
  [./stress_ppt]
    type = ComputeLinearElasticStress
    base_name = ppt
  [../]
  [./global_stress]
    type = TwoPhaseStressMaterial
    base_A = matrix
    base_B = ppt
  [../]
  [./interface_stress]
    type = ComputeSurfaceTensionKKS
    v = eta
    kappa_name = kappa
    w = 0.0033
  [../]
[]

[BCs]
  [./left_r]
    type = DirichletBC
    variable = disp_x
    boundary = left
    value = 0
  [../]
[]

#
# Precondition using handcoded off-diagonal terms
#
[Preconditioning]
  [./full]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -sub_pc_type   -sub_pc_factor_shift_type'
  petsc_options_value = 'asm       lu            nonzero'

  l_max_its = 30
  nl_max_its = 10
  l_tol = 1.0e-4
  nl_rel_tol = 1.0e-9
  nl_abs_tol = 1.0e-10

  num_steps = 2
  dt = 0.5
[]

[Outputs]
  exodus = true
  [./csv]
    type = CSV
    execute_on = 'final'
  [../]
[]

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

#
# Fig. 7 input for 10.1016/j.commatsci.2017.02.017
# D. Schwen et al./Computational Materials Science 132 (2017) 36-45
# Solid gray curve (1)
# Eigenstrain and elastic energies ar computed per phase and then interpolated.
# Supply the RADIUS parameter (10-35) on the command line to generate data
# for all curves in the plot.
#

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 32
  xmin = 0
  xmax = 100
  second_order = true
[]

[Problem]
  coord_type = RSPHERICAL
[]

[GlobalParams]
  displacements = 'disp_r'
[]

[Functions]
  [./diff]
    type = ParsedFunction
    value = '${RADIUS}-pos_c'
    vars = pos_c
    vals = pos_c
  [../]
[]

# 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'
    execute_on = 'INITIAL TIMESTEP_END'
  [../]
[]

[Variables]
  # Solute concentration variable
  [./c]
    [./InitialCondition]
      type = SmoothCircleIC
      invalue = 1
      outvalue = 0
      x1 = 0
      y1 = 0
      radius = ${RADIUS}
      int_width = 3
    [../]
  [../]
  [./w]
  [../]

  # Phase order parameter
  [./eta]
    [./InitialCondition]
      type = SmoothCircleIC
      invalue = 1
      outvalue = 0
      x1 = 0
      y1 = 0
      radius = ${RADIUS}
      int_width = 3
    [../]
  [../]

  # Mesh displacement
  [./disp_r]
    order = SECOND
  [../]

  [./Fe_fit]
    order = SECOND
  [../]
[]

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

  # Split Cahn-Hilliard kernels
  [./c_res]
    type = SplitCHParsed
    variable = c
    f_name = F
    args = 'eta'
    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
  [./detadt]
    type = TimeDerivative
    variable = eta
  [../]
  [./ACBulk1]
    type = AllenCahn
    variable = eta
    args = 'c'
    mob_name = L
    f_name = F
  [../]
  [./ACInterface]
    type = ACInterface
    variable = eta
    mob_name = L
    kappa_name = kappa_eta
  [../]

  [./Fe]
    type = MaterialPropertyValue
    prop_name = Fe
    variable = Fe_fit
  [../]

  [./autoadjust]
    type = MaskedBodyForce
    variable = w
    function = diff
    mask = mask
  [../]
[]

[Materials]
  # declare a few constants, such as mobilities (L,M) and interface gradient prefactors (kappa*)
  [./consts]
    type = GenericConstantMaterial
    prop_names  = 'M   L   kappa_c kappa_eta'
    prop_values = '1.0 1.0 0.5     1'
  [../]

  # forcing function mask
  [./mask]
    type = ParsedMaterial
    f_name = mask
    function = grad/dt
    material_property_names = 'grad dt'
  [../]
  [./grad]
    type = VariableGradientMaterial
    variable = c
    prop = grad
  [../]
  [./time]
    type = TimeStepMaterial
  [../]

  # global mechanical properties
  [./elasticity_tensor_1]
    type = ComputeElasticityTensor
    C_ijkl = '1 1'
    base_name = phase1
    fill_method = symmetric_isotropic
  [../]
  [./elasticity_tensor_2]
    type = ComputeElasticityTensor
    C_ijkl = '1 1'
    base_name = phase2
    fill_method = symmetric_isotropic
  [../]
  [./strain_1]
    type = ComputeRSphericalSmallStrain
    base_name = phase1
  [../]
  [./strain_2]
    type = ComputeRSphericalSmallStrain
    base_name = phase2
    eigenstrain_names = eigenstrain
  [../]
  [./stress_1]
    type = ComputeLinearElasticStress
    base_name = phase1
  [../]
  [./stress_2]
    type = ComputeLinearElasticStress
    base_name = phase2
  [../]

  # eigenstrain per phase
  [./eigenstrain2]
    type = ComputeEigenstrain
    eigen_base = '0.05 0.05 0.05 0 0 0'
    base_name = phase2
    eigenstrain_name = eigenstrain
  [../]

  # switching functions
  [./switching]
    type = SwitchingFunctionMaterial
    function_name = h
    eta = eta
    h_order = SIMPLE
  [../]
  [./barrier]
    type = BarrierFunctionMaterial
    eta = eta
  [../]

  # chemical free energies
  [./chemical_free_energy_1]
    type = DerivativeParsedMaterial
    f_name = Fc1
    function = 'c^2'
    args = 'c'
    derivative_order = 2
  [../]
  [./chemical_free_energy_2]
    type = DerivativeParsedMaterial
    f_name = Fc2
    function = '(1-c)^2'
    args = 'c'
    derivative_order = 2
  [../]

  # elastic free energies
  [./elastic_free_energy_1]
    type = ElasticEnergyMaterial
    f_name = Fe1
    args = ''
    base_name = phase1
    derivative_order = 2
  [../]
  [./elastic_free_energy_2]
    type = ElasticEnergyMaterial
    f_name = Fe2
    args = ''
    base_name = phase2
    derivative_order = 2
  [../]

  # per phase free energies
  [./free_energy_1]
    type = DerivativeSumMaterial
    f_name = F1
    sum_materials = 'Fc1 Fe1'
    args = 'c'
    derivative_order = 2
  [../]
  [./free_energy_2]
    type = DerivativeSumMaterial
    f_name = F2
    sum_materials = 'Fc2 Fe2'
    args = 'c'
    derivative_order = 2
  [../]

  # global chemical free energy
  [./global_free_energy]
    type = DerivativeTwoPhaseMaterial
    f_name = F
    fa_name = F1
    fb_name = F2
    eta = eta
    args = 'c'
    W = 4
  [../]

  # global stress
  [./global_stress]
    type = TwoPhaseStressMaterial
    base_A = phase1
    base_B = phase2
  [../]

  [./elastic_free_energy]
    type = DerivativeTwoPhaseMaterial
    f_name = Fe
    fa_name = Fe1
    fb_name = Fe2
    eta = eta
    args = 'c'
    W = 0
  [../]
[]

[BCs]
  [./left_r]
    type = DirichletBC
    variable = disp_r
    boundary = 'left'
    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'
    outputs = 'table console'
  [../]
  [./total_solute]
    type = ElementIntegralVariablePostprocessor
    variable = c
    execute_on = 'INITIAL TIMESTEP_END'
    outputs = 'table console'
  [../]
  [./pos_c]
    type = FindValueOnLine
    start_point = '0 0 0'
    end_point = '100 0 0'
    v = c
    target = 0.582
    tol = 1e-8
    execute_on = 'INITIAL TIMESTEP_END'
    outputs = 'table console'
  [../]
  [./pos_eta]
    type = FindValueOnLine
    start_point = '0 0 0'
    end_point = '100 0 0'
    v = eta
    target = 0.5
    tol = 1e-8
    execute_on = 'INITIAL TIMESTEP_END'
    outputs = 'table console'
  [../]
  [./c_min]
    type = ElementExtremeValue
    value_type = min
    variable = c
    execute_on = 'INITIAL TIMESTEP_END'
    outputs = 'table console'
  [../]
[]

[VectorPostprocessors]
  [./line]
    type = LineValueSampler
    variable = 'Fe_fit c w'
    start_point = '0 0 0'
    end_point =   '100 0 0'
    num_points = 5000
    sort_by = x
    outputs = vpp
    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 = 15
  l_tol = 1.0e-4
  nl_rel_tol = 1.0e-8
  nl_abs_tol = 2.0e-9
  start_time = 0.0
  end_time = 100000.0

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

  [./Adaptivity]
    initial_adaptivity = 5
    interval = 10
    max_h_level = 5
    refine_fraction = 0.9
    coarsen_fraction = 0.1
  [../]
[]

[Outputs]
  print_linear_residuals = false
  perf_graph = true
  execute_on = 'INITIAL TIMESTEP_END'
  [./table]
    type = CSV
    delimiter = ' '
    file_base = radius_${RADIUS}/energy_pp
  [../]
  [./vpp]
    type = CSV
    delimiter = ' '
    sync_times = '10 50 100 500 1000 5000 10000 50000 100000'
    sync_only = true
    time_data = true
    file_base = radius_${RADIUS}/energy_vpp
  [../]
[]

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

# KKS phase-field model coupled with elasticity using the Voigt-Taylor scheme as
# described in L.K. Aagesen et al., Computational Materials Science, 140, 10-21 (2017)
# Original run #170329e

[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 640
  ny = 1
  nz = 1
  xmin = -10
  xmax = 10
  ymin = 0
  ymax = 0.03125
  zmin = 0
  zmax = 0.03125
  elem_type = HEX8
[]


[Variables]
  # order parameter
  [./eta]
    order = FIRST
    family = LAGRANGE
  [../]

  # solute concentration
  [./c]
    order = FIRST
    family = LAGRANGE
  [../]

  # chemical potential
  [./w]
    order = FIRST
    family = LAGRANGE
  [../]

  # solute phase concentration (matrix)
  [./cm]
    order = FIRST
    family = LAGRANGE
  [../]
  # solute phase concentration (precipitate)
  [./cp]
    order = FIRST
    family = LAGRANGE
  [../]
  [./disp_x]
    order = FIRST
    family = LAGRANGE
  [../]
  [./disp_y]
    order = FIRST
    family = LAGRANGE
  [../]
  [./disp_z]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[ICs]
  [./eta_ic]
    variable = eta
    type = FunctionIC
    function = ic_func_eta
    block = 0
  [../]
  [./c_ic]
    variable = c
    type = FunctionIC
    function = ic_func_c
    block = 0
  [../]
  [./w_ic]
    variable = w
    type = ConstantIC
    value = 0.00991
    block = 0
  [../]
  [./cm_ic]
    variable = cm
    type = ConstantIC
    value = 0.131
    block = 0
  [../]
  [./cp_ic]
    variable = cp
    type = ConstantIC
    value = 0.236
    block = 0
  [../]
[]

[Functions]
  [./ic_func_eta]
    type = ParsedFunction
    value = '0.5*(1.0+tanh((x)/delta_eta/sqrt(2.0)))'
    vars = 'delta_eta'
    vals = '0.8034'
  [../]
  [./ic_func_c]
    type = ParsedFunction
    value = '0.2388*(0.5*(1.0+tanh(x/delta/sqrt(2.0))))^3*(6*(0.5*(1.0+tanh(x/delta/sqrt(2.0))))^2-15*(0.5*(1.0+tanh(x/delta/sqrt(2.0))))+10)+0.1338*(1-(0.5*(1.0+tanh(x/delta/sqrt(2.0))))^3*(6*(0.5*(1.0+tanh(x/delta/sqrt(2.0))))^2-15*(0.5*(1.0+tanh(x/delta/sqrt(2.0))))+10))'
    vars = 'delta'
    vals = '0.8034'
  [../]
  [./psi_eq_int]
    type = ParsedFunction
    value = 'volume*psi_alpha'
    vars = 'volume psi_alpha'
    vals = 'volume psi_alpha'
  [../]
  [./gamma]
    type = ParsedFunction
    value = '(psi_int - psi_eq_int) / dy / dz'
    vars = 'psi_int psi_eq_int dy       dz'
    vals = 'psi_int psi_eq_int 0.03125  0.03125'
  [../]
[]

[AuxVariables]
  [./sigma11]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./sigma22]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./sigma33]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./e11]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./e12]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./e22]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./e33]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./e_el11]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./e_el12]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./e_el22]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./f_el]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./eigen_strain00]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./Fglobal]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./psi]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[AuxKernels]
  [./matl_sigma11]
    type = RankTwoAux
    rank_two_tensor = stress
    index_i = 0
    index_j = 0
    variable = sigma11
  [../]
  [./matl_sigma22]
    type = RankTwoAux
    rank_two_tensor = stress
    index_i = 1
    index_j = 1
    variable = sigma22
  [../]
  [./matl_sigma33]
    type = RankTwoAux
    rank_two_tensor = stress
    index_i = 2
    index_j = 2
    variable = sigma33
  [../]
  [./matl_e11]
    type = RankTwoAux
    rank_two_tensor = total_strain
    index_i = 0
    index_j = 0
    variable = e11
  [../]
  [./matl_e12]
    type = RankTwoAux
    rank_two_tensor = total_strain
    index_i = 0
    index_j = 1
    variable = e12
  [../]
  [./matl_e22]
    type = RankTwoAux
    rank_two_tensor = total_strain
    index_i = 1
    index_j = 1
    variable = e22
  [../]
  [./matl_e33]
    type = RankTwoAux
    rank_two_tensor = total_strain
    index_i = 2
    index_j = 2
    variable = e33
  [../]
  [./f_el]
    type = MaterialRealAux
    variable = f_el
    property = f_el_mat
    execute_on = timestep_end
  [../]
  [./GlobalFreeEnergy]
    variable = Fglobal
    type = KKSGlobalFreeEnergy
    fa_name = fm
    fb_name = fp
    w = 0.0264
    kappa_names = kappa
    interfacial_vars = eta
  [../]
  [./psi_potential]
    variable = psi
    type = ParsedAux
    args = 'Fglobal w c f_el sigma11 e11'
    function = 'Fglobal - w*c + f_el - sigma11*e11'
  [../]
[]


[BCs]
  [./left_x]
    type = DirichletBC
    variable = disp_x
    boundary = left
    value = 0
  [../]
  [./right_x]
    type = DirichletBC
    variable = disp_x
    boundary = right
    value = 0
  [../]
  [./front_y]
    type = DirichletBC
    variable = disp_y
    boundary = front
    value = 0
  [../]
  [./back_y]
    type = DirichletBC
    variable = disp_y
    boundary = back
    value = 0
  [../]
  [./top_z]
    type = DirichletBC
    variable = disp_z
    boundary = top
    value = 0
  [../]
  [./bottom_z]
    type = DirichletBC
    variable = disp_z
    boundary = bottom
    value = 0
  [../]
[]

[Materials]
  # Chemical free energy of the matrix
  [./fm]
    type = DerivativeParsedMaterial
    f_name = fm
    args = 'cm'
    function = '6.55*(cm-0.13)^2'
  [../]
# Elastic energy of the matrix
  [./elastic_free_energy_m]
    type = ElasticEnergyMaterial
    base_name = matrix
    f_name = fe_m
    args = ' '
    outputs = exodus
  [../]
# Total free energy of the matrix
  [./Total_energy_matrix]
    type = DerivativeSumMaterial
    f_name = f_total_matrix
    sum_materials = 'fm fe_m'
    args = 'cm'
  [../]

  # Free energy of the precipitate phase
  [./fp]
    type = DerivativeParsedMaterial
    f_name = fp
    args = 'cp'
    function = '6.55*(cp-0.235)^2'
  [../]

# Elastic energy of the precipitate
  [./elastic_free_energy_p]
    type = ElasticEnergyMaterial
    base_name = ppt
    f_name = fe_p
    args = ' '
    outputs = exodus
  [../]

# Total free energy of the precipitate
  [./Total_energy_ppt]
    type = DerivativeSumMaterial
    f_name = f_total_ppt
    sum_materials = 'fp fe_p'
    args = 'cp'
  [../]

  # Total elastic energy
    [./Total_elastic_energy]
      type = DerivativeTwoPhaseMaterial
      eta = eta
      f_name = f_el_mat
      fa_name = fe_m
      fb_name = fe_p
      outputs = exodus
      W = 0
    [../]

  # h(eta)
  [./h_eta]
    type = SwitchingFunctionMaterial
    h_order = HIGH
    eta = eta
  [../]

  # g(eta)
  [./g_eta]
    type = BarrierFunctionMaterial
    g_order = SIMPLE
    eta = eta
  [../]

  # constant properties
  [./constants]
    type = GenericConstantMaterial
    prop_names  = 'M   L   kappa     misfit'
    prop_values = '0.7 0.7 0.01704   0.00377'
  [../]

  #Mechanical properties
  [./Stiffness_matrix]
    type = ComputeElasticityTensor
    C_ijkl = '103.3 74.25 74.25 103.3 74.25 103.3 46.75 46.75 46.75'
    base_name = matrix
    fill_method = symmetric9
  [../]
  [./Stiffness_ppt]
    type = ComputeElasticityTensor
    C_ijkl = '100.7 71.45 71.45 100.7 71.45 100.7 50.10 50.10 50.10'
    base_name = ppt
    fill_method = symmetric9
  [../]
  [./stress_matrix]
    type = ComputeLinearElasticStress
    base_name = matrix
  [../]
  [./stress_ppt]
    type = ComputeLinearElasticStress
    base_name = ppt
  [../]
  [./strain_matrix]
    type = ComputeSmallStrain
    displacements = 'disp_x disp_y disp_z'
    base_name = matrix
  [../]
  [./strain_ppt]
    type = ComputeSmallStrain
    displacements = 'disp_x disp_y disp_z'
    base_name = ppt
    eigenstrain_names = 'eigenstrain_ppt'
  [../]
  [./eigen_strain]
    type = ComputeEigenstrain
    base_name = ppt
    eigen_base = '1 1 1 0 0 0'
    prefactor = misfit
    eigenstrain_name = 'eigenstrain_ppt'
  [../]
  [./global_stress]
    type = TwoPhaseStressMaterial
    base_A = matrix
    base_B = ppt
  [../]
  [./global_strain]
    type = ComputeSmallStrain
    displacements = 'disp_x disp_y disp_z'
  [../]
[]

[Kernels]
  [./TensorMechanics]
    displacements = 'disp_x disp_y disp_z'
  [../]

  # enforce c = (1-h(eta))*cm + h(eta)*cp
  [./PhaseConc]
    type = KKSPhaseConcentration
    ca       = cm
    variable = cp
    c        = c
    eta      = eta
  [../]

  # enforce pointwise equality of chemical potentials
  [./ChemPotVacancies]
    type = KKSPhaseChemicalPotential
    variable = cm
    cb       = cp
    fa_name  = f_total_matrix
    fb_name  = f_total_ppt
  [../]

  #
  # Cahn-Hilliard Equation
  #
  [./CHBulk]
    type = KKSSplitCHCRes
    variable = c
    ca       = cm
    fa_name  = f_total_matrix
    w        = w
  [../]

  [./dcdt]
    type = CoupledTimeDerivative
    variable = w
    v = c
  [../]
  [./ckernel]
    type = SplitCHWRes
    mob_name = M
    variable = w
  [../]

  #
  # Allen-Cahn Equation
  #
  [./ACBulkF]
    type = KKSACBulkF
    variable = eta
    fa_name  = f_total_matrix
    fb_name  = f_total_ppt
    w        = 0.0264
    args = 'cp cm'
  [../]
  [./ACBulkC]
    type = KKSACBulkC
    variable = eta
    ca       = cm
    cb       = cp
    fa_name  = f_total_matrix
  [../]
  [./ACInterface]
    type = ACInterface
    variable = eta
    kappa_name = kappa
  [../]
  [./detadt]
    type = TimeDerivative
    variable = eta
  [../]
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -sub_pc_type   -sub_pc_factor_shift_type'
  petsc_options_value = 'asm       ilu            nonzero'

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

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

[VectorPostprocessors]
  #[./eta]
  #  type =  LineValueSampler
  #  start_point = '-10 0 0'
  #  end_point = '10 0 0'
  #  variable = eta
  #  num_points = 321
  #  sort_by =  id
  #[../]
  #[./eta_position]
  #  type = FindValueOnLineSample
  #  vectorpostprocessor = eta
  #  variable_name = eta
  #  search_value = 0.5
  #[../]
#  [./f_el]
#    type =  LineMaterialRealSampler
#    start = '-20 0 0'
#    end   = '20 0 0'
#    sort_by = id
#    property = f_el
#  [../]
#  [./f_el_a]
#    type =  LineMaterialRealSampler
#    start = '-20 0 0'
#    end   = '20 0 0'
#    sort_by = id
#    property = fe_m
#  [../]
#  [./f_el_b]
#    type =  LineMaterialRealSampler
#    start = '-20 0 0'
#    end   = '20 0 0'
#    sort_by = id
#    property = fe_p
#  [../]
#  [./h_out]
#    type =  LineMaterialRealSampler
#    start = '-20 0 0'
#    end   = '20 0 0'
#    sort_by = id
#    property = h
#  [../]
#  [./fm_out]
#    type =  LineMaterialRealSampler
#    start = '-20 0 0'
#    end   = '20 0 0'
#    sort_by = id
#    property = fm
#  [../]
[]

[Postprocessors]
  [./f_el_int]
    type = ElementIntegralMaterialProperty
    mat_prop = f_el_mat
  [../]
  [./c_alpha]
    type =  SideAverageValue
    boundary = left
    variable = c
  [../]
  [./c_beta]
    type =  SideAverageValue
    boundary = right
    variable = c
  [../]
  [./e11_alpha]
    type =  SideAverageValue
    boundary = left
    variable = e11
  [../]
  [./e11_beta]
    type =  SideAverageValue
    boundary = right
    variable = e11
  [../]
  [./s11_alpha]
    type =  SideAverageValue
    boundary = left
    variable = sigma11
  [../]
  [./s22_alpha]
    type =  SideAverageValue
    boundary = left
    variable = sigma22
  [../]
  [./s33_alpha]
    type =  SideAverageValue
    boundary = left
    variable = sigma33
  [../]
  [./s11_beta]
    type =  SideAverageValue
    boundary = right
    variable = sigma11
  [../]
  [./s22_beta]
    type =  SideAverageValue
    boundary = right
    variable = sigma22
  [../]
  [./s33_beta]
    type =  SideAverageValue
    boundary = right
    variable = sigma33
  [../]
  [./f_el_alpha]
    type =  SideAverageValue
    boundary = left
    variable = f_el
  [../]
  [./f_el_beta]
    type =  SideAverageValue
    boundary = right
    variable = f_el
  [../]
  [./f_c_alpha]
    type =  SideAverageValue
    boundary = left
    variable = Fglobal
  [../]
  [./f_c_beta]
    type =  SideAverageValue
    boundary = right
    variable = Fglobal
  [../]
  [./chem_pot_alpha]
    type =  SideAverageValue
    boundary = left
    variable = w
  [../]
  [./chem_pot_beta]
    type =  SideAverageValue
    boundary = right
    variable = w
  [../]
  [./psi_alpha]
    type =  SideAverageValue
    boundary = left
    variable = psi
  [../]
  [./psi_beta]
    type =  SideAverageValue
    boundary = right
    variable = psi
  [../]
  [./total_energy]
    type = ElementIntegralVariablePostprocessor
    variable = Fglobal
  [../]
  # Get simulation cell size from postprocessor
  [./volume]
    type = ElementIntegralMaterialProperty
    mat_prop = 1
  [../]
  [./psi_eq_int]
    type = FunctionValuePostprocessor
    function = psi_eq_int
  [../]
  [./psi_int]
    type = ElementIntegralVariablePostprocessor
    variable = psi
  [../]
  [./gamma]
    type = FunctionValuePostprocessor
    function = gamma
  [../]
[]

#
# Precondition using handcoded off-diagonal terms
#
[Preconditioning]
  [./full]
    type = SMP
    full = true
  [../]
[]


[Outputs]
  [./exodus]
    type = Exodus
    interval = 20
  [../]
  [./csv]
    type = CSV
    execute_on = 'final'
  [../]
#[./console]
#    type = Console
#    output_file = true
#  [../]
[]

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

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 14
  ny = 10
  nz = 0
  xmin = 10
  xmax = 40
  ymin = 15
  ymax = 35
  elem_type = QUAD4
[]

[Variables]
  [./c]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = SmoothCircleIC
      x1 = 25.0
      y1 = 25.0
      radius = 6.0
      invalue = 0.9
      outvalue = 0.1
      int_width = 3.0
    [../]
  [../]
  [./w]
    order = FIRST
    family = LAGRANGE
  [../]
  [./eta]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = SmoothCircleIC
      x1 = 30.0
      y1 = 25.0
      radius = 4.0
      invalue = 0.9
      outvalue = 0.1
      int_width = 2.0
    [../]
  [../]
[]

[Kernels]
  [./detadt]
    type = TimeDerivative
    variable = eta
  [../]
  [./ACBulk]
    type = AllenCahn
    variable = eta
    args = c
    f_name = F
  [../]
  [./ACInterface]
    type = ACInterface
    variable = eta
    kappa_name = kappa_eta
  [../]

  [./c_res]
    type = SplitCHParsed
    variable = c
    f_name = F
    kappa_name = kappa_c
    w = w
    args = 'eta'
  [../]
  [./w_res]
    type = SplitCHWRes
    variable = w
    mob_name = M
  [../]
  [./time]
    type = CoupledTimeDerivative
    variable = w
    v = c
  [../]
[]

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

[Materials]
  [./consts]
    type = GenericConstantMaterial
    prop_names  = 'L kappa_eta'
    prop_values = '1 1        '
  [../]
  [./consts2]
    type = GenericConstantMaterial
    prop_names  = 'M kappa_c'
    prop_values = '1 1'
  [../]

  [./switching]
    type = SwitchingFunctionMaterial
    eta = eta
    h_order = SIMPLE
  [../]
  [./barrier]
    type = BarrierFunctionMaterial
    eta = eta
    g_order = SIMPLE
  [../]

  [./free_energy_A]
    type = DerivativeParsedMaterial
    f_name = Fa
    args = 'c'
    function = '(c-0.1)^2*(c-1)^2 + c*0.01'
    derivative_order = 2
    enable_jit = true
  [../]
  [./free_energy_B]
    type = DerivativeParsedMaterial
    f_name = Fb
    args = 'c'
    function = 'c^2*(c-0.9)^2 + (1-c)*0.01'
    derivative_order = 2
    enable_jit = true
  [../]

  [./free_energy]
    type = DerivativeTwoPhaseMaterial
    f_name = F
    fa_name = Fa
    fb_name = Fb
    args = 'c'
    eta = eta
    derivative_order = 2
    outputs = exodus
    output_properties = 'F dF/dc dF/deta d^2F/dc^2 d^2F/dcdeta d^2F/deta^2'
  [../]
[]

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

[Executioner]
  type = Transient
  scheme = 'bdf2'
  solve_type = 'NEWTON'

  l_max_its = 15
  l_tol = 1.0e-4

  nl_max_its = 10
  nl_rel_tol = 1.0e-11

  start_time = 0.0
  num_steps = 1
  dt = 0.1
[]

[Outputs]
  exodus = true
[]

(modules/phase_field/test/tests/GrandPotentialPFM/GrandPotentialPFM.i)

# this input file test the implementation of the grand-potential phase-field model based on M.Plapp PRE 84,031601(2011)
# in this simple example, the liquid and solid free energies are parabola with the same curvature and the material properties are constant
# Note that this example also test The SusceptibilityTimeDerivative kernels
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 16
  ny = 16
  xmax = 32
  ymax = 32
[]

[GlobalParams]
  radius = 20.0
  int_width = 4.0
  x1 = 0
  y1 = 0
[]

[Variables]
  [./w]
  [../]
  [./eta]
  [../]
[]

[ICs]
  [./w]
    type = SmoothCircleIC
    variable = w
    # note w = A*(c-cleq), A = 1.0, cleq = 0.0 ,i.e., w = c (in the matrix/liquid phase)
    outvalue = -0.2
    invalue = 0.2
  [../]
  [./eta]
    type = SmoothCircleIC
    variable = eta
    outvalue = 0.0
    invalue = 1.0
  [../]
[]

[Kernels]
  [./w_dot]
    type = SusceptibilityTimeDerivative
    variable = w
    f_name = chi
    args = '' # in this case chi (the susceptibility) is simply a constant
  [../]
  [./Diffusion]
    type = MatDiffusion
    variable = w
    diffusivity = D
    args = ''
  [../]
  [./coupled_etadot]
    type = CoupledSusceptibilityTimeDerivative
    variable = w
    v = eta
    f_name = ft
    args = 'eta'
  [../]
  [./AC_bulk]
    type = AllenCahn
    variable = eta
    f_name = F
    args = 'w'
  [../]
  [./AC_int]
    type = ACInterface
    variable = eta
  [../]
  [./e_dot]
    type = TimeDerivative
    variable = eta
  [../]
[]

[Materials]
  [./constants]
    type = GenericConstantMaterial
    prop_names  = 'kappa_op  D    L    chi  cseq   cleq   A'
    prop_values = '4.0       1.0  1.0  1.0  0.0  1.0  1.0'
  [../]

  [./liquid_GrandPotential]
    type = DerivativeParsedMaterial
    function = '-0.5 * w^2/A - cleq * w'
    args = 'w'
    f_name = f1
    material_property_names = 'cleq A'
  [../]
  [./solid_GrandPotential]
    type = DerivativeParsedMaterial
    function = '-0.5 * w^2/A - cseq * w'
    args = 'w'
    f_name = f2
    material_property_names = 'cseq A'
  [../]
  [./switching_function]
    type = SwitchingFunctionMaterial
    eta = eta
    h_order = HIGH
  [../]
  [./barrier_function]
    type = BarrierFunctionMaterial
    eta = eta
  [../]
  [./cs]
    type = DerivativeParsedMaterial
    args = 'w'
    f_name = cs
    material_property_names = 'A cseq'
    function = 'w/A + cseq' # since w = A*(c-cseq)
    derivative_order = 2
  [../]
  [./cl]
    type = DerivativeParsedMaterial
    args = 'w'
    f_name = cl
    material_property_names = 'A cleq'
    function = 'w/A + cleq' # since w = A*(c-cleq)
    derivative_order = 2
  [../]
  [./total_GrandPotential]
    type = DerivativeTwoPhaseMaterial
    args = 'w'
    eta = eta
    fa_name = f1
    fb_name = f2
    derivative_order = 2
    W = 1.0
  [../]
  [./coupled_eta_function]
    type = DerivativeParsedMaterial
    function = '(cs - cl) * dh'
    args = 'eta w'
    f_name = ft
    material_property_names = 'cs cl dh:=D[h,eta]'
    derivative_order = 1
    outputs = exodus
  [../]

  [./concentration]
    type = ParsedMaterial
    f_name = c
    material_property_names = 'dF:=D[F,w]'
    function = '-dF'
    outputs = exodus
  [../]
[]

[Postprocessors]
  [./C]
    type = ElementIntegralMaterialProperty
    mat_prop = c
    execute_on = 'initial timestep_end'
  [../]
[]

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

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = NEWTON

  l_max_its = 15
  l_tol = 1e-3
  nl_max_its = 15
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-8

  num_steps = 5
  dt = 10.0
[]

[Outputs]
  exodus = true
  csv = true
  execute_on = 'TIMESTEP_END'
[]

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

#
# Fig. 7 input for 10.1016/j.commatsci.2017.02.017
# D. Schwen et al./Computational Materials Science 132 (2017) 36-45
# Dashed black curve (2)
# Eigenstrain is globally applied. Single global elastic free energies.
# Supply the RADIUS parameter (10-35) on the command line to generate data
# for all curves in the plot.
#

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 32
  xmin = 0
  xmax = 100
  second_order = true
[]

[Problem]
  coord_type = RSPHERICAL
[]

[GlobalParams]
  displacements = 'disp_r'
[]

[Functions]
  [./diff]
    type = ParsedFunction
    value = '${RADIUS}-pos_c'
    vars = pos_c
    vals = pos_c
  [../]
[]

# 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'
    execute_on = 'INITIAL TIMESTEP_END'
  [../]
[]

[Variables]
  # Solute concentration variable
  [./c]
    [./InitialCondition]
      type = SmoothCircleIC
      invalue = 1
      outvalue = 0
      x1 = 0
      y1 = 0
      radius = ${RADIUS}
      int_width = 3
    [../]
  [../]
  [./w]
  [../]

  # Phase order parameter
  [./eta]
    [./InitialCondition]
      type = SmoothCircleIC
      invalue = 1
      outvalue = 0
      x1 = 0
      y1 = 0
      radius = ${RADIUS}
      int_width = 3
    [../]
  [../]

  [./Fe_fit]
    order = SECOND
  [../]
[]

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

[Kernels]
  # Split Cahn-Hilliard kernels
  [./c_res]
    type = SplitCHParsed
    variable = c
    f_name = F
    args = 'eta'
    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
  [./detadt]
    type = TimeDerivative
    variable = eta
  [../]
  [./ACBulk1]
    type = AllenCahn
    variable = eta
    args = 'c'
    mob_name = L
    f_name = F
  [../]
  [./ACInterface]
    type = ACInterface
    variable = eta
    mob_name = L
    kappa_name = kappa_eta
  [../]

  [./Fe]
    type = MaterialPropertyValue
    prop_name = Fe
    variable = Fe_fit
  [../]

  [./autoadjust]
    type = MaskedBodyForce
    variable = w
    function = diff
    mask = mask
  [../]
[]

[Materials]
  # declare a few constants, such as mobilities (L,M) and interface gradient prefactors (kappa*)
  [./consts]
    type = GenericConstantMaterial
    prop_names  = 'M   L   kappa_c kappa_eta'
    prop_values = '1.0 1.0 0.5     1'
  [../]

  # forcing function mask
  [./mask]
    type = ParsedMaterial
    f_name = mask
    function = grad/dt
    material_property_names = 'grad dt'
  [../]
  [./grad]
    type = VariableGradientMaterial
    variable = c
    prop = grad
  [../]
  [./time]
    type = TimeStepMaterial
  [../]

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

  # eigenstrain as a function of phase
  [./eigenstrain]
    type = ComputeVariableEigenstrain
    eigen_base = '0.05 0.05 0.05 0 0 0'
    prefactor = h
    args = eta
    eigenstrain_name = eigenstrain
  [../]

  # switching functions
  [./switching]
    type = SwitchingFunctionMaterial
    function_name = h
    eta = eta
    h_order = SIMPLE
  [../]
  [./barrier]
    type = BarrierFunctionMaterial
    eta = eta
  [../]

  # chemical free energies
  [./chemical_free_energy_1]
    type = DerivativeParsedMaterial
    f_name = Fc1
    function = 'c^2'
    args = 'c'
    derivative_order = 2
  [../]
  [./chemical_free_energy_2]
    type = DerivativeParsedMaterial
    f_name = Fc2
    function = '(1-c)^2'
    args = 'c'
    derivative_order = 2
  [../]

  # global chemical free energy
  [./chemical_free_energy]
    type = DerivativeTwoPhaseMaterial
    f_name = Fc
    fa_name = Fc1
    fb_name = Fc2
    eta = eta
    args = 'c'
    W = 4
  [../]

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

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

[BCs]
  [./left_r]
    type = DirichletBC
    variable = disp_r
    boundary = 'left'
    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'
    outputs = 'table console'
  [../]
  [./total_solute]
    type = ElementIntegralVariablePostprocessor
    variable = c
    execute_on = 'INITIAL TIMESTEP_END'
    outputs = 'table console'
  [../]
  [./pos_c]
    type = FindValueOnLine
    start_point = '0 0 0'
    end_point = '100 0 0'
    v = c
    target = 0.582
    tol = 1e-8
    execute_on = 'INITIAL TIMESTEP_END'
    outputs = 'table console'
  [../]
  [./pos_eta]
    type = FindValueOnLine
    start_point = '0 0 0'
    end_point = '100 0 0'
    v = eta
    target = 0.5
    tol = 1e-8
    execute_on = 'INITIAL TIMESTEP_END'
    outputs = 'table console'
  [../]
  [./c_min]
    type = ElementExtremeValue
    value_type = min
    variable = c
    execute_on = 'INITIAL TIMESTEP_END'
    outputs = 'table console'
  [../]
[]

[VectorPostprocessors]
  [./line]
    type = LineValueSampler
    variable = 'Fe_fit c w'
    start_point = '0 0 0'
    end_point =   '100 0 0'
    num_points = 5000
    sort_by = x
    outputs = vpp
    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 = 15
  l_tol = 1.0e-4
  nl_rel_tol = 1.0e-8
  nl_abs_tol = 2.0e-9
  start_time = 0.0
  end_time = 100000.0

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

  [./Adaptivity]
    initial_adaptivity = 5
    interval = 10
    max_h_level = 5
    refine_fraction = 0.9
    coarsen_fraction = 0.1
  [../]
[]

[Outputs]
  print_linear_residuals = false
  perf_graph = true
  execute_on = 'INITIAL TIMESTEP_END'
  [./table]
    type = CSV
    delimiter = ' '
    file_base = radius_${RADIUS}/eigenstrain_pp
  [../]
  [./vpp]
    type = CSV
    delimiter = ' '
    sync_times = '10 50 100 500 1000 5000 10000 50000 100000'
    sync_only = true
    time_data = true
    file_base = radius_${RADIUS}/eigenstrain_vpp
  [../]
[]

Input Parameters
Input Files

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