- reaction_rateThe reaction rate used with the kernelC++ Type:MaterialPropertyName Unit:(no unit assumed) Controllable:No Description:The reaction rate used with the kernel 
- variableThe name of the variable that this residual object operates onC++ Type:NonlinearVariableName Unit:(no unit assumed) Controllable:No Description:The name of the variable that this residual object operates on 
MatReaction
Kernel to add -L*v, where L=reaction rate, v=variable
Implements  where  (mob_name) is a reaction rate,  is either a coupled variable (v) or - if not explicitly specified - the nonlinear variable the kernel is operating on.
Note the negative sign, which does not appear in Reaction or CoefReaction.
Input Parameters
- argsVector of nonlinear variable arguments this object depends onC++ Type:std::vector<VariableName> Unit:(no unit assumed) Controllable:No Description:Vector of nonlinear variable arguments this object depends on 
- blockThe list of blocks (ids or names) that this object will be appliedC++ Type:std::vector<SubdomainName> Controllable:No Description:The list of blocks (ids or names) that this object will be applied 
- displacementsThe displacementsC++ Type:std::vector<VariableName> Unit:(no unit assumed) Controllable:No Description:The displacements 
- matrix_onlyFalseWhether this object is only doing assembly to matrices (no vectors)Default:False C++ Type:bool Controllable:No Description:Whether this object is only doing assembly to matrices (no vectors) 
- vSet this to make v a coupled variable, otherwise it will use the kernel's nonlinear variable for vC++ Type:std::vector<VariableName> Unit:(no unit assumed) Controllable:No Description:Set this to make v a coupled variable, otherwise it will use the kernel's nonlinear variable for v 
Optional Parameters
- absolute_value_vector_tagsThe tags for the vectors this residual object should fill with the absolute value of the residual contributionC++ Type:std::vector<TagName> Controllable:No Description:The tags for the vectors this residual object should fill with the absolute value of the residual contribution 
- extra_matrix_tagsThe extra tags for the matrices this Kernel should fillC++ Type:std::vector<TagName> Controllable:No Description:The extra tags for the matrices this Kernel should fill 
- extra_vector_tagsThe extra tags for the vectors this Kernel should fillC++ Type:std::vector<TagName> Controllable:No Description:The extra tags for the vectors this Kernel should fill 
- matrix_tagssystemThe tag for the matrices this Kernel should fillDefault:system C++ Type:MultiMooseEnum Options:nontime, system Controllable:No Description:The tag for the matrices this Kernel should fill 
- vector_tagsnontimeThe tag for the vectors this Kernel should fillDefault:nontime C++ Type:MultiMooseEnum Options:nontime, time Controllable:No Description:The tag for the vectors this Kernel should fill 
Contribution To Tagged Field Data Parameters
- control_tagsAdds user-defined labels for accessing object parameters via control logic.C++ Type:std::vector<std::string> Controllable:No 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> Unit:(no unit assumed) Controllable:No 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 Controllable:Yes Description:Set the enabled status of the MooseObject. 
- implicitTrueDetermines whether this object is calculated using an implicit or explicit formDefault:True C++ Type:bool Controllable:No 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> Unit:(no unit assumed) Controllable:No 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.) 
- search_methodnearest_node_connected_sidesChoice of search algorithm. All options begin by finding the nearest node in the primary boundary to a query point in the secondary boundary. In the default nearest_node_connected_sides algorithm, primary boundary elements are searched iff that nearest node is one of their nodes. This is fast to determine via a pregenerated node-to-elem map and is robust on conforming meshes. In the optional all_proximate_sides algorithm, primary boundary elements are searched iff they touch that nearest node, even if they are not topologically connected to it. This is more CPU-intensive but is necessary for robustness on any boundary surfaces which has disconnections (such as Flex IGA meshes) or non-conformity (such as hanging nodes in adaptively h-refined meshes).Default:nearest_node_connected_sides C++ Type:MooseEnum Options:nearest_node_connected_sides, all_proximate_sides Controllable:No Description:Choice of search algorithm. All options begin by finding the nearest node in the primary boundary to a query point in the secondary boundary. In the default nearest_node_connected_sides algorithm, primary boundary elements are searched iff that nearest node is one of their nodes. This is fast to determine via a pregenerated node-to-elem map and is robust on conforming meshes. In the optional all_proximate_sides algorithm, primary boundary elements are searched iff they touch that nearest node, even if they are not topologically connected to it. This is more CPU-intensive but is necessary for robustness on any boundary surfaces which has disconnections (such as Flex IGA meshes) or non-conformity (such as hanging nodes in adaptively h-refined meshes). 
- seed0The seed for the master random number generatorDefault:0 C++ Type:unsigned int Controllable:No 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 Controllable:No 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
- prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.C++ Type:MaterialPropertyName Unit:(no unit assumed) Controllable:No Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character. 
- use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.Default:False C++ Type:bool Controllable:No Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction. 
Material Property Retrieval Parameters
Input Files
- (test/tests/neml2/blocks_different_model.i)
- (modules/phase_field/test/tests/phase_field_kernels/CoupledAllenCahn.i)
- (modules/phase_field/test/tests/SimpleACInterface/SimpleCoupledACInterface.i)
- (test/tests/neml2/multiple_input_files.i)
- (test/tests/neml2/blocks_same_model.i)
- (modules/stochastic_tools/test/tests/actions/parameter_study_action/sub_eigen.i)
- (modules/phase_field/test/tests/phase_field_kernels/CoupledCoefAllenCahn.i)
- (modules/fsi/test/tests/2d-small-strain-transient/ad-fsi-flat-channel.i)
- (modules/phase_field/test/tests/phase_field_kernels/MatGradSquareCoupled.i)
- (modules/phase_field/test/tests/KKS_system/lagrange_multiplier.i)
- (modules/combined/examples/publications/rapid_dev/fig6.i)
- (test/tests/neml2/parameter.i)
- (test/tests/neml2/error.i)
- (modules/phase_field/test/tests/KKS_system/kks_multiphase.i)
- (modules/phase_field/test/tests/electrochem_sintering/ElectrochemicalSintering_test.i)
- (test/tests/neml2/mesh_change.i)
- (modules/fsi/test/tests/2d-small-strain-transient/fsi_flat_channel.i)
- (test/tests/neml2/simple_scheduler.i)
- (modules/phase_field/examples/multiphase/GrandPotential3Phase_masscons.i)
- (modules/phase_field/test/tests/mobility_derivative/coupledmatdiffusion.i)
- (test/tests/neml2/custom_model.i)
- (test/tests/neml2/simple_scheduler_async.i)
Child Objects
(test/tests/neml2/blocks_different_model.i)
[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 10
  []
  [A]
    type = SubdomainBoundingBoxGenerator
    input = gmg
    block_id = 1
    block_name = A
    bottom_left = '0 0 0'
    top_right = '0.5 1 0'
  []
  [B]
    type = SubdomainBoundingBoxGenerator
    input = A
    block_id = 2
    block_name = B
    bottom_left = '0.5 0 0'
    top_right = '1 1 0'
  []
[]
[NEML2]
  input = 'models/custom_model.i'
  verbose = true
  device = 'cpu'
  [A]
    model = 'model_A'
    block = 'A'
    moose_input_types = 'VARIABLE MATERIAL'
    moose_inputs = '     a        b'
    neml2_inputs = '     forces/A forces/B'
    moose_output_types = 'MATERIAL           MATERIAL'
    moose_outputs = '     neml2_sum          neml2_product'
    neml2_outputs = '     state/internal/sum state/internal/product'
    moose_derivative_types = 'MATERIAL'
    moose_derivatives = 'neml2_dproduct_da'
    neml2_derivatives = 'state/internal/product forces/A'
    export_outputs = 'neml2_sum neml2_product neml2_dproduct_da'
    export_output_targets = 'exodus; exodus; exodus'
  []
  [B]
    model = 'model_B'
    block = 'B'
    moose_input_types = 'VARIABLE MATERIAL'
    moose_inputs = '     a        b'
    neml2_inputs = '     forces/A forces/B'
    moose_output_types = 'MATERIAL           MATERIAL'
    moose_outputs = '     neml2_sum          neml2_product'
    neml2_outputs = '     state/internal/sum state/internal/product'
    moose_derivative_types = 'MATERIAL'
    moose_derivatives = 'neml2_dproduct_da'
    neml2_derivatives = 'state/internal/product forces/A'
    export_outputs = 'neml2_sum neml2_product neml2_dproduct_da'
    export_output_targets = 'exodus; exodus; exodus'
  []
[]
[Variables]
  [u]
  []
[]
[Kernels]
  [diffusion]
    type = Diffusion
    variable = u
  []
  [reaction_1]
    type = MatReaction
    variable = u
    reaction_rate = neml2_sum
  []
  [reaction_2]
    type = MatReaction
    variable = u
    reaction_rate = neml2_product
  []
  [reaction_3]
    type = MatReaction
    variable = u
    reaction_rate = neml2_dproduct_da
  []
[]
[AuxVariables]
  [a]
  []
[]
[ICs]
  [a]
    type = FunctionIC
    variable = a
    function = 'x'
  []
[]
[Materials]
  [b]
    type = GenericFunctionMaterial
    prop_names = 'b'
    prop_values = 'y+t'
    outputs = 'exodus'
  []
[]
[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  num_steps = 2
[]
[Outputs]
  exodus = true
[]
(modules/phase_field/test/tests/phase_field_kernels/CoupledAllenCahn.i)
#
# Test the coupled Allen-Cahn Bulk kernel
#
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 10
  ny = 10
  xmax = 12
  ymax = 12
  elem_type = QUAD4
[]
[Variables]
  [./w]
  [../]
  [./eta]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = SmoothCircleIC
      x1 = 0.0
      y1 = 0.0
      radius = 6.0
      invalue = 0.9
      outvalue = 0.1
      int_width = 3.0
    [../]
  [../]
[]
[Kernels]
  [./detadt]
    type = TimeDerivative
    variable = eta
  [../]
  [./ACBulk]
    type = CoupledAllenCahn
    variable = w
    v = eta
    f_name = F
  [../]
  [./W]
    type = Reaction
    variable = w
  [../]
  [./CoupledBulk]
    type = MatReaction
    variable = eta
    v = w
    reaction_rate = L
  [../]
  [./ACInterface]
    type = ACInterface
    variable = eta
    kappa_name = 1
  [../]
[]
[Materials]
  [./consts]
    type = GenericConstantMaterial
    prop_names  = 'L'
    prop_values = '1'
  [../]
  [./free_energy]
    type = DerivativeParsedMaterial
    property_name = F
    coupled_variables = 'eta'
    expression = '2 * eta^2 * (1-eta)^2 - 0.2*eta'
    derivative_order = 2
  [../]
[]
[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]
[Executioner]
  type = Transient
  scheme = 'bdf2'
  solve_type = 'PJFNK'
  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 = 2
  dt = 0.5
[]
[Debug]
  show_var_residual_norms = true
[]
[Outputs]
  hide = w
  file_base = AllenCahn_out
  exodus = true
[]
(modules/phase_field/test/tests/SimpleACInterface/SimpleCoupledACInterface.i)
#
# Test the coupled Allen-Cahn Bulk kernel
#
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 20
  ny = 20
  nz = 0
  xmin = 0
  xmax = 50
  ymin = 0
  ymax = 50
  zmin = 0
  zmax = 50
  elem_type = QUAD4
  uniform_refine = 1
[]
[Variables]
  [./w]
  [../]
  [./eta]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = SmoothCircleIC
      x1 = 25.0
      y1 = 25.0
      radius = 6.0
      invalue = 1.0
      outvalue = 0.0
      int_width = 5.0
    [../]
  [../]
[]
[Kernels]
  [./detadt]
    type = TimeDerivative
    variable = eta
  [../]
  [./ACBulk]
    type = AllenCahn
    variable = eta
    f_name = F
  [../]
  [./CoupledBulk]
    type = MatReaction
    variable = eta
    v = w
    reaction_rate = L
  [../]
  [./W]
    type = Reaction
    variable = w
  [../]
  [./CoupledACInterface]
    type = SimpleCoupledACInterface
    variable = w
    v = eta
    kappa_name = 1
  [../]
[]
[Materials]
  [./consts]
    type = GenericConstantMaterial
    prop_names  = 'L'
    prop_values = '1'
  [../]
  [./free_energy]
    type = DerivativeParsedMaterial
    property_name = F
    coupled_variables = 'eta'
    expression = 'eta^2 * (1-eta)^2'
    derivative_order = 2
  [../]
[]
[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]
[Executioner]
  type = Transient
  scheme = 'bdf2'
  solve_type = 'PJFNK'
  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 = 2
  dt = 2
[]
[Debug]
  show_var_residual_norms = true
[]
[Outputs]
  hide = w
  exodus = true
[]
(test/tests/neml2/multiple_input_files.i)
[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 10
  []
  [A]
    type = SubdomainBoundingBoxGenerator
    input = gmg
    block_id = 1
    block_name = A
    bottom_left = '0 0 0'
    top_right = '0.5 1 0'
  []
  [B]
    type = SubdomainBoundingBoxGenerator
    input = A
    block_id = 2
    block_name = B
    bottom_left = '0.5 0 0'
    top_right = '1 1 0'
  []
[]
[NEML2]
  verbose = true
  device = 'cpu'
  [A]
    input = 'models/custom_model.i'
    model = 'model_A'
    block = 'A'
    moose_input_types = 'VARIABLE MATERIAL'
    moose_inputs = '     a        b'
    neml2_inputs = '     forces/A forces/B'
    moose_output_types = 'MATERIAL           MATERIAL'
    moose_outputs = '     neml2_sum          neml2_product'
    neml2_outputs = '     state/internal/sum state/internal/product'
    moose_derivative_types = 'MATERIAL'
    moose_derivatives = 'neml2_dproduct_da'
    neml2_derivatives = 'state/internal/product forces/A'
    export_outputs = 'neml2_sum neml2_product neml2_dproduct_da'
    export_output_targets = 'exodus; exodus; exodus'
  []
  [B]
    input = 'models/custom_model_2.i'
    model = 'model_B'
    block = 'B'
    moose_input_types = 'VARIABLE MATERIAL'
    moose_inputs = '     a        b'
    neml2_inputs = '     forces/A forces/B'
    moose_output_types = 'MATERIAL           MATERIAL'
    moose_outputs = '     neml2_sum          neml2_product'
    neml2_outputs = '     state/internal/sum state/internal/product'
    moose_derivative_types = 'MATERIAL'
    moose_derivatives = 'neml2_dproduct_da'
    neml2_derivatives = 'state/internal/product forces/A'
    export_outputs = 'neml2_sum neml2_product neml2_dproduct_da'
    export_output_targets = 'exodus; exodus; exodus'
  []
[]
[Variables]
  [u]
  []
[]
[Kernels]
  [diffusion]
    type = Diffusion
    variable = u
  []
  [reaction_1]
    type = MatReaction
    variable = u
    reaction_rate = neml2_sum
  []
  [reaction_2]
    type = MatReaction
    variable = u
    reaction_rate = neml2_product
  []
  [reaction_3]
    type = MatReaction
    variable = u
    reaction_rate = neml2_dproduct_da
  []
[]
[AuxVariables]
  [a]
  []
[]
[ICs]
  [a]
    type = FunctionIC
    variable = a
    function = 'x'
  []
[]
[Materials]
  [b]
    type = GenericFunctionMaterial
    prop_names = 'b'
    prop_values = 'y+t'
    outputs = 'exodus'
  []
[]
[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  num_steps = 2
[]
[Outputs]
  exodus = true
[]
(test/tests/neml2/blocks_same_model.i)
[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 10
  []
  [A]
    type = SubdomainBoundingBoxGenerator
    input = gmg
    block_id = 1
    block_name = A
    bottom_left = '0 0 0'
    top_right = '0.5 1 0'
  []
  [B]
    type = SubdomainBoundingBoxGenerator
    input = A
    block_id = 2
    block_name = B
    bottom_left = '0.5 0 0'
    top_right = '1 1 0'
  []
[]
[NEML2]
  input = 'models/custom_model.i'
  verbose = true
  device = 'cpu'
  moose_input_types = 'VARIABLE MATERIAL'
  moose_inputs = '     a        b'
  neml2_inputs = '     forces/A forces/B'
  moose_output_types = 'MATERIAL           MATERIAL'
  moose_outputs = '     neml2_sum          neml2_product'
  neml2_outputs = '     state/internal/sum state/internal/product'
  moose_derivative_types = 'MATERIAL'
  moose_derivatives = 'neml2_dproduct_da'
  neml2_derivatives = 'state/internal/product forces/A'
  export_outputs = 'neml2_sum neml2_product neml2_dproduct_da'
  export_output_targets = 'exodus; exodus; exodus'
  [A]
    model = 'model'
    block = 'A'
  []
  [B]
    model = 'model'
    block = 'B'
  []
[]
[Variables]
  [u]
  []
[]
[Kernels]
  [diffusion]
    type = Diffusion
    variable = u
  []
  [reaction_1]
    type = MatReaction
    variable = u
    reaction_rate = neml2_sum
  []
  [reaction_2]
    type = MatReaction
    variable = u
    reaction_rate = neml2_product
  []
  [reaction_3]
    type = MatReaction
    variable = u
    reaction_rate = neml2_dproduct_da
  []
[]
[AuxVariables]
  [a]
  []
[]
[ICs]
  [a]
    type = FunctionIC
    variable = a
    function = 'x'
  []
[]
[Materials]
  [b]
    type = GenericFunctionMaterial
    prop_names = 'b'
    prop_values = 'y+t'
    outputs = 'exodus'
  []
[]
[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  num_steps = 2
[]
[Outputs]
  exodus = true
[]
(modules/stochastic_tools/test/tests/actions/parameter_study_action/sub_eigen.i)
[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = 10
    ymin = 0
    ymax = 10
    elem_type = QUAD4
    nx = 8
    ny = 8
  []
[]
[Variables]
  [u]
  []
[]
[Kernels]
  [diff]
    type = MatDiffusion
    variable = u
    diffusivity = D
  []
  [rhs]
    type = MatReaction
    variable = u
    reaction_rate = L
    extra_vector_tags = 'eigen'
  []
[]
[Materials]
  [mat]
    type = GenericFunctionMaterial
    prop_names = 'D L'
    prop_values = 'diff_fun react_fun'
  []
[]
[Functions]
  [diff_fun]
    type = ConstantFunction
    value = 1
  []
  [react_fun]
    type = ConstantFunction
    value = 1
  []
[]
[BCs]
  [homogeneous]
    type = DirichletBC
    variable = u
    boundary = '0 1 2 3'
    value = 0
  []
  [eigen]
    type = EigenDirichletBC
    variable = u
    boundary = '0 1 2 3'
  []
[]
[Executioner]
  type = Eigenvalue
[]
[VectorPostprocessors]
  [eigenvalues]
    type = Eigenvalues
  []
[]
[Postprocessors]
  [eigenvalue]
    type = VectorPostprocessorComponent
    vectorpostprocessor = eigenvalues
    vector_name = eigen_values_real
    index = 0
  []
[]
(modules/phase_field/test/tests/phase_field_kernels/CoupledCoefAllenCahn.i)
#
# Test the CoefReaction kernel (which adds -L*v to the residual) for the case
# where v is a coupled variable
#
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 15
  ny = 15
  nz = 0
  xmin = 0
  xmax = 50
  ymin = 0
  ymax = 50
  zmin = 0
  zmax = 50
  elem_type = QUAD4
[]
[Variables]
  [./w]
  [../]
  [./eta]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = SmoothCircleIC
      x1 = 25.0
      y1 = 25.0
      radius = 6.0
      invalue = 1.0
      outvalue = 0.0
      int_width = 3.0
    [../]
  [../]
[]
[Kernels]
  [./detadt]
    type = TimeDerivative
    variable = eta
  [../]
  [./ACBulk]
    type = CoupledAllenCahn
    variable = w
    v = eta
    f_name = F
    mob_name = 1
  [../]
  [./W]
    type = MatReaction
    variable = w
    reaction_rate = -1
  [../]
  [./CoupledBulk]
    type = MatReaction
    variable = eta
    v = w
    reaction_rate = L
  [../]
  [./ACInterface]
    type = ACInterface
    variable = eta
    kappa_name = 1
    mob_name = L
    coupled_variables = w
  [../]
[]
[Materials]
  [./mobility]
    type = DerivativeParsedMaterial
    property_name  = L
    coupled_variables = 'eta w'
    expression = '(1.5-eta)^2+(1.5-w)^2'
    derivative_order = 2
  [../]
  [./free_energy]
    type = DerivativeParsedMaterial
    property_name = F
    coupled_variables = 'eta'
    expression = 'eta^2 * (1-eta)^2'
    derivative_order = 2
  [../]
[]
[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]
[Executioner]
  type = Transient
  scheme = 'bdf2'
  solve_type = 'PJFNK'
  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 = 2
  dt = 0.5
[]
[Outputs]
  hide = w
  exodus = true
[]
(modules/fsi/test/tests/2d-small-strain-transient/ad-fsi-flat-channel.i)
[GlobalParams]
  displacements = 'disp_x disp_y'
  order = FIRST
  preset = false
  use_displaced_mesh = true
[]
[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = 3.0
    ymin = 0
    ymax = 1.0
    nx = 10
    ny = 15
    elem_type = QUAD4
  []
  [subdomain1]
    type = SubdomainBoundingBoxGenerator
    bottom_left = '0.0 0.5 0'
    block_id = 1
    top_right = '3.0 1.0 0'
    input = gmg
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '0'
    paired_block = '1'
    new_boundary = 'master0_interface'
    input = subdomain1
  []
  [break_boundary]
    type = BreakBoundaryOnSubdomainGenerator
    input = interface
  []
[]
[Variables]
  [vel]
    block = 0
    family = LAGRANGE_VEC
  []
  [p]
    block = 0
    order = FIRST
  []
  [disp_x]
  []
  [disp_y]
  []
  [vel_x_solid]
    block = 1
  []
  [vel_y_solid]
    block = 1
  []
[]
[Kernels]
  [mass]
    type = INSADMass
    variable = p
    block = 0
  []
  [mass_pspg]
    type = INSADMassPSPG
    variable = p
    block = 0
  []
  [momentum_time]
    type = INSADMomentumTimeDerivative
    variable = vel
    block = 0
  []
  [momentum_convection]
    type = INSADMomentumAdvection
    variable = vel
    block = 0
  []
  [momentum_viscous]
    type = INSADMomentumViscous
    variable = vel
    block = 0
  []
  [momentum_pressure]
    type = INSADMomentumPressure
    variable = vel
    pressure = p
    integrate_p_by_parts = true
    block = 0
  []
  [momentum_supg]
    type = INSADMomentumSUPG
    variable = vel
    material_velocity = relative_velocity
    block = 0
  []
  [momentum_mesh]
    type = INSADMomentumMeshAdvection
    variable = vel
    disp_x = 'disp_x'
    disp_y = 'disp_y'
    block = 0
  []
  [disp_x_fluid]
    type = Diffusion
    variable = disp_x
    block = 0
    use_displaced_mesh = false
  []
  [disp_y_fluid]
    type = Diffusion
    variable = disp_y
    block = 0
    use_displaced_mesh = false
  []
  [accel_tensor_x]
    type = CoupledTimeDerivative
    variable = disp_x
    v = vel_x_solid
    block = 1
    use_displaced_mesh = false
  []
  [accel_tensor_y]
    type = CoupledTimeDerivative
    variable = disp_y
    v = vel_y_solid
    block = 1
    use_displaced_mesh = false
  []
  [vxs_time_derivative_term]
    type = CoupledTimeDerivative
    variable = vel_x_solid
    v = disp_x
    block = 1
    use_displaced_mesh = false
  []
  [vys_time_derivative_term]
    type = CoupledTimeDerivative
    variable = vel_y_solid
    v = disp_y
    block = 1
    use_displaced_mesh = false
  []
  [source_vxs]
    type = MatReaction
    variable = vel_x_solid
    block = 1
    reaction_rate = 1
    use_displaced_mesh = false
  []
  [source_vys]
    type = MatReaction
    variable = vel_y_solid
    block = 1
    reaction_rate = 1
    use_displaced_mesh = false
  []
[]
[InterfaceKernels]
  [penalty]
    type = ADPenaltyVelocityContinuity
    variable = vel
    fluid_velocity = vel
    displacements = 'disp_x disp_y'
    solid_velocities = 'vel_x_solid vel_y_solid'
    boundary = master0_interface
    penalty = 1e6
  []
[]
[Physics/SolidMechanics/QuasiStatic]
  [solid_domain]
    strain = SMALL
    incremental = false
    # generate_output = 'strain_xx strain_yy strain_zz' ## Not at all necessary, but nice
    block = '1'
  []
[]
[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e2
    poissons_ratio = 0.3
    block = '1'
    use_displaced_mesh = false
  []
  [small_stress]
    type = ComputeLinearElasticStress
    block = 1
  []
  [const]
    type = ADGenericConstantMaterial
    block = 0
    prop_names = 'rho mu'
    prop_values = '1  1'
  []
  [ins_mat]
    type = INSADTauMaterial
    velocity = vel
    pressure = p
    block = 0
  []
[]
[BCs]
  [fluid_bottom]
    type = ADVectorFunctionDirichletBC
    variable = vel
    boundary = 'bottom'
    function_x = 0
    function_y = 0
  []
  [fluid_left]
    type = ADVectorFunctionDirichletBC
    variable = vel
    boundary = 'left_to_0'
    function_x = 'inlet_func'
    function_y = 0
    # The displacements actually affect the result of the function evaluation so in order to eliminate the impact
    # on the Jacobian we set 'use_displaced_mesh = false' here
    use_displaced_mesh = false
  []
  [no_disp_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'bottom top left_to_1 right_to_1 left_to_0 right_to_0'
    value = 0
  []
  [no_disp_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom top left_to_1 right_to_1 left_to_0 right_to_0'
    value = 0
  []
  [solid_x_no_slip]
    type = DirichletBC
    variable = vel_x_solid
    boundary = 'top left_to_1 right_to_1'
    value = 0.0
  []
  [solid_y_no_slip]
    type = DirichletBC
    variable = vel_y_solid
    boundary = 'top left_to_1 right_to_1'
    value = 0.0
  []
[]
[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]
[Executioner]
  type = Transient
  num_steps = 5
  # num_steps = 60
  dt = 0.1
  dtmin = 0.1
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu       NONZERO'
  line_search = none
  nl_rel_tol = 1e-50
  nl_abs_tol = 1e-10
[]
[Outputs]
  exodus = true
[]
[Functions]
  [inlet_func]
    type = ParsedFunction
    expression = '(-16 * (y - 0.25)^2 + 1) * (1 + cos(t))'
  []
[]
(modules/phase_field/test/tests/phase_field_kernels/MatGradSquareCoupled.i)
#
# Test the MatGradSquareCoupled kernel
#
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 15
  ny = 15
  nz = 0
  xmin = 0
  xmax = 50
  ymin = 0
  ymax = 50
  zmin = 0
  zmax = 50
  elem_type = QUAD4
[]
[Variables]
  [./w]
  [../]
  [./eta]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = SmoothCircleIC
      x1 = 25.0
      y1 = 25.0
      radius = 6.0
      invalue = 1.0
      outvalue = 0.0
      int_width = 3.0
    [../]
  [../]
[]
[Kernels]
  [./detadt]
    type = TimeDerivative
    variable = eta
  [../]
  [./ACBulk]
    type = CoupledAllenCahn
    variable = w
    v = eta
    f_name = F
    mob_name = 1
  [../]
  [./W]
    type = MatReaction
    variable = w
    reaction_rate = -1
  [../]
  [./CoupledBulk]
    type = MatReaction
    variable = eta
    v = w
    reaction_rate = L
  [../]
  [./ACInterface]
    type = ACInterface
    variable = eta
    kappa_name = 1
    mob_name = L
    coupled_variables = w
  [../]
# MatGradSquareCoupled kernel
  [./nabla_eta]
    type = MatGradSquareCoupled
    variable = w
    elec_potential = eta
    prefactor = 0.5
  [../]
[]
[Materials]
  [./mobility]
    type = DerivativeParsedMaterial
    property_name  = L
    coupled_variables = 'eta w'
    expression = '(1.5-eta)^2+(1.5-w)^2'
    derivative_order = 2
  [../]
  [./free_energy]
    type = DerivativeParsedMaterial
    property_name = F
    coupled_variables = 'eta'
    expression = 'eta^2 * (1-eta)^2'
    derivative_order = 2
  [../]
[]
[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]
[Executioner]
  type = Transient
  scheme = 'bdf2'
  solve_type = 'PJFNK'
  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 = 2
  dt = 0.5
[]
[Outputs]
  hide = w
  exodus = true
  console = true
[]
(modules/phase_field/test/tests/KKS_system/lagrange_multiplier.i)
#
# This test ensures that the equilibrium solution using two order parameters with a
# Lagrange multiplier constraint is identical to the dedicated two phase formulation
# in two_phase.i
#
[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 20
  xmax = 5
[]
[AuxVariables]
  [Fglobal]
    order = CONSTANT
    family = MONOMIAL
  []
[]
[Variables]
  # concentration
  [c]
    order = FIRST
    family = LAGRANGE
    [InitialCondition]
      type = FunctionIC
      function = x/5
    []
  []
  # order parameter 1
  [eta1]
    order = FIRST
    family = LAGRANGE
    initial_condition = 0.5
  []
  # order parameter 2
  [eta2]
    order = FIRST
    family = LAGRANGE
    initial_condition = 0.5
  []
  # phase concentration 1
  [c1]
    order = FIRST
    family = LAGRANGE
    initial_condition = 0.2
  []
  # phase concentration 2
  [c2]
    order = FIRST
    family = LAGRANGE
    initial_condition = 0.5
  []
  # Lagrange multiplier
  [lambda]
    order = FIRST
    family = LAGRANGE
    initial_condition = 0.0
  []
[]
[Materials]
  # simple toy free energies
  [f1] # = fd
    type = DerivativeParsedMaterial
    property_name = F1
    coupled_variables = 'c1'
    expression = '(0.9-c1)^2'
  []
  [f2] # = fm
    type = DerivativeParsedMaterial
    property_name = F2
    coupled_variables = 'c2'
    expression = '(0.1-c2)^2'
  []
  # Switching functions for each phase
  [h1_eta]
    type = SwitchingFunctionMaterial
    h_order = HIGH
    eta = eta1
    function_name = h1
  []
  [h2_eta]
    type = SwitchingFunctionMaterial
    h_order = HIGH
    eta = eta2
    function_name = h2
  []
  # Coefficients for diffusion equation
  [Dh1]
    type = DerivativeParsedMaterial
    material_property_names = 'D h1'
    expression = D*h1
    property_name = Dh1
  []
  [Dh2]
    type = DerivativeParsedMaterial
    material_property_names = 'D h2'
    expression = D*h2
    property_name = Dh2
  []
  # Barrier functions for each phase
  [g1]
    type = BarrierFunctionMaterial
    g_order = SIMPLE
    eta = eta1
    function_name = g1
  []
  [g2]
    type = BarrierFunctionMaterial
    g_order = SIMPLE
    eta = eta2
    function_name = g2
  []
  # constant properties
  [constants]
    type = GenericConstantMaterial
    prop_names  = 'D   L   kappa'
    prop_values = '0.7 0.7 0.2'
  []
[]
[Kernels]
  #Kernels for diffusion equation
  [diff_time]
    type = TimeDerivative
    variable = c
  []
  [diff_c1]
    type = MatDiffusion
    variable = c
    diffusivity = Dh1
    v = c1
  []
  [diff_c2]
    type = MatDiffusion
    variable = c
    diffusivity = Dh2
    v = c2
  []
  # Kernels for Allen-Cahn equation for eta1
  [deta1dt]
    type = TimeDerivative
    variable = eta1
  []
  [ACBulkF1]
    type = KKSMultiACBulkF
    variable = eta1
    Fj_names = 'F1 F2 '
    hj_names = 'h1 h2 '
    gi_name = g1
    eta_i = eta1
    wi = 0.2
    coupled_variables = 'c1 c2 eta2'
  []
  [ACBulkC1]
    type = KKSMultiACBulkC
    variable = eta1
    Fj_names = 'F1 F2'
    hj_names = 'h1 h2'
    cj_names = 'c1 c2'
    eta_i = eta1
    coupled_variables = 'eta2'
  []
  [ACInterface1]
    type = ACInterface
    variable = eta1
    kappa_name = kappa
  []
  [multipler1]
    type = MatReaction
    variable = eta1
    v = lambda
    reaction_rate = L
  []
  # Kernels for the Lagrange multiplier equation
  [mult_lambda]
    type = MatReaction
    variable = lambda
    reaction_rate = 2
  []
  [mult_ACBulkF_1]
    type = KKSMultiACBulkF
    variable = lambda
    Fj_names = 'F1 F2 '
    hj_names = 'h1 h2 '
    gi_name = g1
    eta_i = eta1
    wi = 0.2
    mob_name = 1
    coupled_variables = 'c1 c2 eta2 '
  []
  [mult_ACBulkC_1]
    type = KKSMultiACBulkC
    variable = lambda
    Fj_names = 'F1 F2'
    hj_names = 'h1 h2'
    cj_names = 'c1 c2'
    eta_i = eta1
    coupled_variables = 'eta2 '
    mob_name = 1
  []
  [mult_CoupledACint_1]
    type = SimpleCoupledACInterface
    variable = lambda
    v = eta1
    kappa_name = kappa
    mob_name = 1
  []
  [mult_ACBulkF_2]
    type = KKSMultiACBulkF
    variable = lambda
    Fj_names = 'F1 F2 '
    hj_names = 'h1 h2 '
    gi_name = g2
    eta_i = eta2
    wi = 0.2
    mob_name = 1
    coupled_variables = 'c1 c2 eta1 '
  []
  [mult_ACBulkC_2]
    type = KKSMultiACBulkC
    variable = lambda
    Fj_names = 'F1 F2'
    hj_names = 'h1 h2'
    cj_names = 'c1 c2'
    eta_i = eta2
    coupled_variables = 'eta1 '
    mob_name = 1
  []
  [mult_CoupledACint_2]
    type = SimpleCoupledACInterface
    variable = lambda
    v = eta2
    kappa_name = kappa
    mob_name = 1
  []
  # Kernels for constraint equation eta1 + eta2 = 1
  # eta2 is the nonlinear variable for the constraint equation
  [eta2reaction]
    type = MatReaction
    variable = eta2
    reaction_rate = 1
  []
  [eta1reaction]
    type = MatReaction
    variable = eta2
    v = eta1
    reaction_rate = 1
  []
  [one]
    type = BodyForce
    variable = eta2
    value = -1.0
  []
  # Phase concentration constraints
  [chempot12]
    type = KKSPhaseChemicalPotential
    variable = c1
    cb = c2
    fa_name = F1
    fb_name = F2
  []
  [phaseconcentration]
    type = KKSMultiPhaseConcentration
    variable = c2
    cj = 'c1 c2'
    hj_names = 'h1 h2'
    etas = 'eta1 eta2'
    c = c
  []
[]
[AuxKernels]
  [Fglobal_total]
    type = KKSMultiFreeEnergy
    Fj_names = 'F1 F2 '
    hj_names = 'h1 h2 '
    gj_names = 'g1 g2 '
    variable = Fglobal
    w = 0.2
    interfacial_vars = 'eta1  eta2 '
    kappa_names      = 'kappa kappa'
  []
[]
[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu       nonzero'
  l_max_its = 30
  nl_max_its = 10
  l_tol = 1.0e-4
  nl_rel_tol = 1.0e-10
  nl_abs_tol = 1.0e-11
  num_steps = 35
  dt = 10
[]
[VectorPostprocessors]
  [c]
    type = LineValueSampler
    variable = c
    start_point = '0 0 0'
    end_point = '5 0 0'
    num_points = 21
    sort_by = x
  []
[]
[Outputs]
  csv = true
  execute_on = FINAL
[]
(modules/combined/examples/publications/rapid_dev/fig6.i)
#
# Fig. 6 input for 10.1016/j.commatsci.2017.02.017
# D. Schwen et al./Computational Materials Science 132 (2017) 36-45
# Three phase interface simulation demonstrating the interfacial stability
# w.r.t. formation of a tspurious third phase
#
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 120
  ny = 120
  nz = 0
  xmin = 0
  xmax = 40
  ymin = 0
  ymax = 40
  zmin = 0
  zmax = 0
  elem_type = QUAD4
[]
[Variables]
  # concentration
  [./c]
  [../]
  # order parameter 1
  [./eta1]
  [../]
  # order parameter 2
  [./eta2]
  [../]
  # order parameter 3
  [./eta3]
  [../]
  # phase concentration 1
  [./c1]
    initial_condition = 0.4
  [../]
  # phase concentration 2
  [./c2]
    initial_condition = 0.5
  [../]
  # phase concentration 3
  [./c3]
    initial_condition = 0.8
  [../]
  # Lagrange multiplier
  [./lambda]
    initial_condition = 0.0
  [../]
[]
[AuxVariables]
  [./T]
    [./InitialCondition]
      type = FunctionIC
      function = 'x-10'
    [../]
  [../]
[]
[Functions]
  [./ic_func_eta1]
    type = ParsedFunction
    expression = '0.5*(1.0+tanh((x-10)/sqrt(2.0))) * 0.5*(1.0+tanh((y-10)/sqrt(2.0)))'
  [../]
  [./ic_func_eta2]
    type = ParsedFunction
    expression = '0.5*(1.0-tanh((x-10)/sqrt(2.0)))'
  [../]
  [./ic_func_eta3]
    type = ParsedFunction
    expression = '1 - 0.5*(1.0-tanh((x-10)/sqrt(2.0)))
              - 0.5*(1.0+tanh((x-10)/sqrt(2.0))) * 0.5*(1.0+tanh((y-10)/sqrt(2.0)))'
  [../]
  [./ic_func_c]
    type = ParsedFunction
    expression = '0.5 * 0.5*(1.0-tanh((x-10)/sqrt(2.0)))
              + 0.4 * 0.5*(1.0+tanh((x-10)/sqrt(2.0))) * 0.5*(1.0+tanh((y-10)/sqrt(2.0)))
              + 0.8 * (1 - 0.5*(1.0-tanh((x-10)/sqrt(2.0)))
                        - 0.5*(1.0+tanh((x-10)/sqrt(2.0))) * 0.5*(1.0+tanh((y-10)/sqrt(2.0))))'
  [../]
[]
[ICs]
  [./eta1]
    variable = eta1
    type = FunctionIC
    function = ic_func_eta1
  [../]
  [./eta2]
    variable = eta2
    type = FunctionIC
    function = ic_func_eta2
  [../]
  [./eta3]
    variable = eta3
    type = FunctionIC
    function = ic_func_eta3
  [../]
  [./c]
    variable = c
    type = FunctionIC
    function = ic_func_c
  [../]
[]
[Materials]
  # simple toy free energies
  [./f1]
    type = DerivativeParsedMaterial
    property_name = F1
    coupled_variables = 'c1'
    expression = '20*(c1-0.4)^2'
  [../]
  [./f2]
    type = DerivativeParsedMaterial
    property_name = F2
    coupled_variables = 'c2 T'
    expression = '20*(c2-0.5)^2 + 0.01*T'
  [../]
  [./f3]
    type = DerivativeParsedMaterial
    property_name = F3
    coupled_variables = 'c3'
    expression = '20*(c3-0.8)^2'
  [../]
  # Switching functions for each phase
  # h1(eta1, eta2, eta3)
  [./h1]
    type = SwitchingFunction3PhaseMaterial
    eta_i = eta1
    eta_j = eta2
    eta_k = eta3
    f_name = h1
  [../]
  # h2(eta1, eta2, eta3)
  [./h2]
    type = SwitchingFunction3PhaseMaterial
    eta_i = eta2
    eta_j = eta3
    eta_k = eta1
    f_name = h2
  [../]
  # h3(eta1, eta2, eta3)
  [./h3]
    type = SwitchingFunction3PhaseMaterial
    eta_i = eta3
    eta_j = eta1
    eta_k = eta2
    f_name = h3
  [../]
  # Coefficients for diffusion equation
  [./Dh1]
    type = DerivativeParsedMaterial
    material_property_names = 'D h1'
    expression = D*h1
    property_name = Dh1
  [../]
  [./Dh2]
    type = DerivativeParsedMaterial
    material_property_names = 'D h2'
    expression = D*h2
    property_name = Dh2
  [../]
  [./Dh3]
    type = DerivativeParsedMaterial
    material_property_names = 'D h3'
    expression = D*h3
    property_name = Dh3
  [../]
  # Barrier functions for each phase
  [./g1]
    type = BarrierFunctionMaterial
    g_order = SIMPLE
    eta = eta1
    function_name = g1
  [../]
  [./g2]
    type = BarrierFunctionMaterial
    g_order = SIMPLE
    eta = eta2
    function_name = g2
  [../]
  [./g3]
    type = BarrierFunctionMaterial
    g_order = SIMPLE
    eta = eta3
    function_name = g3
  [../]
  # constant properties
  [./constants]
    type = GenericConstantMaterial
    prop_names  = 'L   kappa  D'
    prop_values = '1.0 1.0    1'
  [../]
[]
[Kernels]
  #Kernels for diffusion equation
  [./diff_time]
    type = TimeDerivative
    variable = c
  [../]
  [./diff_c1]
    type = MatDiffusion
    variable = c
    diffusivity = Dh1
    v = c1
  [../]
  [./diff_c2]
    type = MatDiffusion
    variable = c
    diffusivity = Dh2
    v = c2
  [../]
  [./diff_c3]
    type = MatDiffusion
    variable = c
    diffusivity = Dh3
    v = c3
  [../]
  # Kernels for Allen-Cahn equation for eta1
  [./deta1dt]
    type = TimeDerivative
    variable = eta1
  [../]
  [./ACBulkF1]
    type = KKSMultiACBulkF
    variable  = eta1
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    gi_name   = g1
    eta_i     = eta1
    wi        = 1.0
    coupled_variables = 'c1 c2 c3 eta2 eta3'
  [../]
  [./ACBulkC1]
    type = KKSMultiACBulkC
    variable  = eta1
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    cj_names  = 'c1 c2 c3'
    eta_i     = eta1
    coupled_variables = 'eta2 eta3'
  [../]
  [./ACInterface1]
    type = ACInterface
    variable = eta1
    kappa_name = kappa
  [../]
  [./multipler1]
    type = MatReaction
    variable = eta1
    v = lambda
    reaction_rate = L
  [../]
  # Kernels for Allen-Cahn equation for eta2
  [./deta2dt]
    type = TimeDerivative
    variable = eta2
  [../]
  [./ACBulkF2]
    type = KKSMultiACBulkF
    variable  = eta2
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    gi_name   = g2
    eta_i     = eta2
    wi        = 1.0
    coupled_variables = 'c1 c2 c3 eta1 eta3'
  [../]
  [./ACBulkC2]
    type = KKSMultiACBulkC
    variable  = eta2
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    cj_names  = 'c1 c2 c3'
    eta_i     = eta2
    coupled_variables = 'eta1 eta3'
  [../]
  [./ACInterface2]
    type = ACInterface
    variable = eta2
    kappa_name = kappa
  [../]
  [./multipler2]
    type = MatReaction
    variable = eta2
    v = lambda
    reaction_rate = L
  [../]
  # Kernels for the Lagrange multiplier equation
  [./mult_lambda]
    type = MatReaction
    variable = lambda
    reaction_rate = 3
  [../]
  [./mult_ACBulkF_1]
    type = KKSMultiACBulkF
    variable  = lambda
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    gi_name   = g1
    eta_i     = eta1
    wi        = 1.0
    mob_name  = 1
    coupled_variables = 'c1 c2 c3 eta2 eta3'
  [../]
  [./mult_ACBulkC_1]
    type = KKSMultiACBulkC
    variable  = lambda
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    cj_names  = 'c1 c2 c3'
    eta_i     = eta1
    coupled_variables = 'eta2 eta3'
    mob_name  = 1
  [../]
  [./mult_CoupledACint_1]
    type = SimpleCoupledACInterface
    variable = lambda
    v = eta1
    kappa_name = kappa
    mob_name = 1
  [../]
  [./mult_ACBulkF_2]
    type = KKSMultiACBulkF
    variable  = lambda
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    gi_name   = g2
    eta_i     = eta2
    wi        = 1.0
    mob_name  = 1
    coupled_variables = 'c1 c2 c3 eta1 eta3'
  [../]
  [./mult_ACBulkC_2]
    type = KKSMultiACBulkC
    variable  = lambda
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    cj_names  = 'c1 c2 c3'
    eta_i     = eta2
    coupled_variables = 'eta1 eta3'
    mob_name  = 1
  [../]
  [./mult_CoupledACint_2]
    type = SimpleCoupledACInterface
    variable = lambda
    v = eta2
    kappa_name = kappa
    mob_name = 1
  [../]
  [./mult_ACBulkF_3]
    type = KKSMultiACBulkF
    variable  = lambda
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    gi_name   = g3
    eta_i     = eta3
    wi        = 1.0
    mob_name  = 1
    coupled_variables = 'c1 c2 c3 eta1 eta2'
  [../]
  [./mult_ACBulkC_3]
    type = KKSMultiACBulkC
    variable  = lambda
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    cj_names  = 'c1 c2 c3'
    eta_i     = eta3
    coupled_variables = 'eta1 eta2'
    mob_name  = 1
  [../]
  [./mult_CoupledACint_3]
    type = SimpleCoupledACInterface
    variable = lambda
    v = eta3
    kappa_name = kappa
    mob_name = 1
  [../]
  # Kernels for constraint equation eta1 + eta2 + eta3 = 1
  # eta3 is the nonlinear variable for the constraint equation
  [./eta3reaction]
    type = MatReaction
    variable = eta3
    reaction_rate = 1
  [../]
  [./eta1reaction]
    type = MatReaction
    variable = eta3
    v = eta1
    reaction_rate = 1
  [../]
  [./eta2reaction]
    type = MatReaction
    variable = eta3
    v = eta2
    reaction_rate = 1
  [../]
  [./one]
    type = BodyForce
    variable = eta3
    value = -1.0
  [../]
  # Phase concentration constraints
  [./chempot12]
    type = KKSPhaseChemicalPotential
    variable = c1
    cb       = c2
    fa_name  = F1
    fb_name  = F2
  [../]
  [./chempot23]
    type = KKSPhaseChemicalPotential
    variable = c2
    cb       = c3
    fa_name  = F2
    fb_name  = F3
  [../]
  [./phaseconcentration]
    type = KKSMultiPhaseConcentration
    variable = c3
    cj = 'c1 c2 c3'
    hj_names = 'h1 h2 h3'
    etas = 'eta1 eta2 eta3'
    c = c
  [../]
[]
[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-10
  nl_abs_tol = 1.0e-11
  num_steps = 1000
  [./TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.2
    optimal_iterations = 10
    iteration_window = 2
  [../]
[]
[Preconditioning]
  active = 'full'
  [./full]
    type = SMP
    full = true
  [../]
  [./mydebug]
    type = FDP
    full = true
  [../]
[]
[Outputs]
  exodus = true
  checkpoint = true
  print_linear_residuals = false
  [./csv]
    type = CSV
    execute_on = 'final'
  [../]
[]
#[VectorPostprocessors]
#  [./c]
#    type =  LineValueSampler
#    start_point = '-25 0 0'
#    end_point = '25 0 0'
#    variable = c
#    num_points = 151
#    sort_by =  id
#    execute_on = timestep_end
#  [../]
#  [./eta1]
#    type =  LineValueSampler
#    start_point = '-25 0 0'
#    end_point = '25 0 0'
#    variable = eta1
#    num_points = 151
#    sort_by =  id
#    execute_on = timestep_end
#  [../]
#  [./eta2]
#    type =  LineValueSampler
#    start_point = '-25 0 0'
#    end_point = '25 0 0'
#    variable = eta2
#    num_points = 151
#    sort_by =  id
#    execute_on = timestep_end
#  [../]
#  [./eta3]
#    type =  LineValueSampler
#    start_point = '-25 0 0'
#    end_point = '25 0 0'
#    variable = eta3
#    num_points = 151
#    sort_by =  id
#    execute_on = timestep_end
#  [../]
#[]
(test/tests/neml2/parameter.i)
[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 10
  []
[]
[NEML2]
  input = 'models/custom_model.i'
  [all]
    model = 'model'
    verbose = true
    device = 'cpu'
    moose_input_types = 'VARIABLE MATERIAL'
    moose_inputs = '     a        b'
    neml2_inputs = '     forces/A forces/B'
    moose_parameter_types = 'MATERIAL VARIABLE'
    moose_parameters = '     p1_mat   p2_var'
    neml2_parameters = '     p1       p2'
    moose_output_types = 'MATERIAL           MATERIAL'
    moose_outputs = '     neml2_sum          neml2_product'
    neml2_outputs = '     state/internal/sum state/internal/product'
    moose_derivative_types = 'MATERIAL'
    moose_derivatives = 'neml2_dproduct_da'
    neml2_derivatives = 'state/internal/product forces/A'
    export_outputs = 'neml2_sum neml2_product neml2_dproduct_da'
    export_output_targets = 'exodus; exodus; exodus'
  []
[]
[Variables]
  [u]
  []
[]
[Kernels]
  [diffusion]
    type = Diffusion
    variable = u
  []
  [reaction_1]
    type = MatReaction
    variable = u
    reaction_rate = neml2_sum
  []
  [reaction_2]
    type = MatReaction
    variable = u
    reaction_rate = neml2_product
  []
  [reaction_3]
    type = MatReaction
    variable = u
    reaction_rate = neml2_dproduct_da
  []
[]
[AuxVariables]
  [a]
  []
  [p2_var]
    initial_condition = 2
  []
[]
[ICs]
  [a]
    type = FunctionIC
    variable = a
    function = 'x'
  []
[]
[Materials]
  [b]
    type = GenericFunctionMaterial
    prop_names = 'b'
    prop_values = 'y+t'
    outputs = 'exodus'
  []
  [p1_mat]
    type = GenericConstantMaterial
    prop_names = 'p1_mat'
    prop_values = 3
  []
[]
[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  num_steps = 2
[]
[Outputs]
  exodus = true
[]
(test/tests/neml2/error.i)
[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 10
  []
[]
[NEML2]
  input = 'models/error.i'
  [all]
    model = 'model'
    verbose = true
    device = 'cpu'
    moose_input_types = 'VARIABLE MATERIAL'
    moose_inputs = '     a        b'
    neml2_inputs = '     forces/A forces/B'
    moose_output_types = 'MATERIAL           MATERIAL'
    moose_outputs = '     neml2_sum          neml2_product'
    neml2_outputs = '     state/internal/sum state/internal/product'
    moose_derivative_types = 'MATERIAL'
    moose_derivatives = 'neml2_dproduct_da'
    neml2_derivatives = 'state/internal/product forces/A'
  []
[]
[Variables]
  [u]
  []
[]
[Kernels]
  [diffusion]
    type = Diffusion
    variable = u
  []
  [reaction_1]
    type = MatReaction
    variable = u
    reaction_rate = neml2_sum
  []
  [reaction_2]
    type = MatReaction
    variable = u
    reaction_rate = neml2_product
  []
  [reaction_3]
    type = MatReaction
    variable = u
    reaction_rate = neml2_dproduct_da
  []
[]
[AuxVariables]
  [a]
  []
[]
[ICs]
  [a]
    type = FunctionIC
    variable = a
    function = 'x'
  []
[]
[Materials]
  [b]
    type = GenericFunctionMaterial
    prop_names = 'b'
    prop_values = 'y+t'
  []
[]
[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  num_steps = 1
  dtmin = 1
[]
(modules/phase_field/test/tests/KKS_system/kks_multiphase.i)
#
# This test is for the 3-phase KKS model
#
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 20
  ny = 20
  nz = 0
  xmin = 0
  xmax = 40
  ymin = 0
  ymax = 40
  zmin = 0
  zmax = 0
  elem_type = QUAD4
[]
[BCs]
  [./Periodic]
    [./all]
      auto_direction = 'x y'
    [../]
  [../]
[]
[AuxVariables]
  [./Energy]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]
[Variables]
  # concentration
  [./c]
    order = FIRST
    family = LAGRANGE
  [../]
  # order parameter 1
  [./eta1]
    order = FIRST
    family = LAGRANGE
  [../]
  # order parameter 2
  [./eta2]
    order = FIRST
    family = LAGRANGE
  [../]
  # order parameter 3
  [./eta3]
    order = FIRST
    family = LAGRANGE
    initial_condition = 0.0
  [../]
  # phase concentration 1
  [./c1]
    order = FIRST
    family = LAGRANGE
    initial_condition = 0.2
  [../]
  # phase concentration 2
  [./c2]
    order = FIRST
    family = LAGRANGE
    initial_condition = 0.5
  [../]
  # phase concentration 3
  [./c3]
    order = FIRST
    family = LAGRANGE
    initial_condition = 0.8
  [../]
  # Lagrange multiplier
  [./lambda]
    order = FIRST
    family = LAGRANGE
    initial_condition = 0.0
  [../]
[]
[ICs]
  [./eta1]
    variable = eta1
    type = SmoothCircleIC
    x1 = 20.0
    y1 = 20.0
    radius = 10
    invalue = 0.9
    outvalue = 0.1
    int_width = 4
  [../]
  [./eta2]
    variable = eta2
    type = SmoothCircleIC
    x1 = 20.0
    y1 = 20.0
    radius = 10
    invalue = 0.1
    outvalue = 0.9
    int_width = 4
  [../]
  [./c]
    variable = c
    type = SmoothCircleIC
    x1 = 20.0
    y1 = 20.0
    radius = 10
    invalue = 0.2
    outvalue = 0.5
    int_width = 2
  [../]
[]
[Materials]
  # simple toy free energies
  [./f1]
    type = DerivativeParsedMaterial
    property_name = F1
    coupled_variables = 'c1'
    expression = '20*(c1-0.2)^2'
  [../]
  [./f2]
    type = DerivativeParsedMaterial
    property_name = F2
    coupled_variables = 'c2'
    expression = '20*(c2-0.5)^2'
  [../]
  [./f3]
    type = DerivativeParsedMaterial
    property_name = F3
    coupled_variables = 'c3'
    expression = '20*(c3-0.8)^2'
  [../]
  # Switching functions for each phase
  # h1(eta1, eta2, eta3)
  [./h1]
    type = SwitchingFunction3PhaseMaterial
    eta_i = eta1
    eta_j = eta2
    eta_k = eta3
    property_name = h1
  [../]
  # h2(eta1, eta2, eta3)
  [./h2]
    type = SwitchingFunction3PhaseMaterial
    eta_i = eta2
    eta_j = eta3
    eta_k = eta1
    property_name = h2
  [../]
  # h3(eta1, eta2, eta3)
  [./h3]
    type = SwitchingFunction3PhaseMaterial
    eta_i = eta3
    eta_j = eta1
    eta_k = eta2
    property_name = h3
  [../]
  # Coefficients for diffusion equation
  [./Dh1]
    type = DerivativeParsedMaterial
    material_property_names = 'D h1'
    expression = D*h1
    property_name = Dh1
  [../]
  [./Dh2]
    type = DerivativeParsedMaterial
    material_property_names = 'D h2'
    expression = D*h2
    property_name = Dh2
  [../]
  [./Dh3]
    type = DerivativeParsedMaterial
    material_property_names = 'D h3'
    expression = D*h3
    property_name = Dh3
  [../]
  # Barrier functions for each phase
  [./g1]
    type = BarrierFunctionMaterial
    g_order = SIMPLE
    eta = eta1
    function_name = g1
  [../]
  [./g2]
    type = BarrierFunctionMaterial
    g_order = SIMPLE
    eta = eta2
    function_name = g2
  [../]
  [./g3]
    type = BarrierFunctionMaterial
    g_order = SIMPLE
    eta = eta3
    function_name = g3
  [../]
  # constant properties
  [./constants]
    type = GenericConstantMaterial
    prop_names  = 'L   kappa  D'
    prop_values = '0.7 1.0    1'
  [../]
[]
[Kernels]
  #Kernels for diffusion equation
  [./diff_time]
    type = TimeDerivative
    variable = c
  [../]
  [./diff_c1]
    type = MatDiffusion
    variable = c
    diffusivity = Dh1
    v = c1
  [../]
  [./diff_c2]
    type = MatDiffusion
    variable = c
    diffusivity = Dh2
    v = c2
  [../]
  [./diff_c3]
    type = MatDiffusion
    variable = c
    diffusivity = Dh3
    v = c3
  [../]
  # Kernels for Allen-Cahn equation for eta1
  [./deta1dt]
    type = TimeDerivative
    variable = eta1
  [../]
  [./ACBulkF1]
    type = KKSMultiACBulkF
    variable  = eta1
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    gi_name   = g1
    eta_i     = eta1
    wi        = 1.0
    coupled_variables      = 'c1 c2 c3 eta2 eta3'
  [../]
  [./ACBulkC1]
    type = KKSMultiACBulkC
    variable  = eta1
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    cj_names  = 'c1 c2 c3'
    eta_i     = eta1
    coupled_variables      = 'eta2 eta3'
  [../]
  [./ACInterface1]
    type = ACInterface
    variable = eta1
    kappa_name = kappa
  [../]
  [./multipler1]
    type = MatReaction
    variable = eta1
    v = lambda
    reaction_rate = L
  [../]
  # Kernels for Allen-Cahn equation for eta2
  [./deta2dt]
    type = TimeDerivative
    variable = eta2
  [../]
  [./ACBulkF2]
    type = KKSMultiACBulkF
    variable  = eta2
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    gi_name   = g2
    eta_i     = eta2
    wi        = 1.0
    coupled_variables      = 'c1 c2 c3 eta1 eta3'
  [../]
  [./ACBulkC2]
    type = KKSMultiACBulkC
    variable  = eta2
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    cj_names  = 'c1 c2 c3'
    eta_i     = eta2
    coupled_variables      = 'eta1 eta3'
  [../]
  [./ACInterface2]
    type = ACInterface
    variable = eta2
    kappa_name = kappa
  [../]
  [./multipler2]
    type = MatReaction
    variable = eta2
    v = lambda
    reaction_rate = L
  [../]
  # Kernels for the Lagrange multiplier equation
  [./mult_lambda]
    type = MatReaction
    variable = lambda
    reaction_rate = 3
  [../]
  [./mult_ACBulkF_1]
    type = KKSMultiACBulkF
    variable  = lambda
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    gi_name   = g1
    eta_i     = eta1
    wi        = 1.0
    mob_name  = 1
    coupled_variables      = 'c1 c2 c3 eta2 eta3'
  [../]
  [./mult_ACBulkC_1]
    type = KKSMultiACBulkC
    variable  = lambda
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    cj_names  = 'c1 c2 c3'
    eta_i     = eta1
    coupled_variables      = 'eta2 eta3'
    mob_name  = 1
  [../]
  [./mult_CoupledACint_1]
    type = SimpleCoupledACInterface
    variable = lambda
    v = eta1
    kappa_name = kappa
    mob_name = 1
  [../]
  [./mult_ACBulkF_2]
    type = KKSMultiACBulkF
    variable  = lambda
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    gi_name   = g2
    eta_i     = eta2
    wi        = 1.0
    mob_name  = 1
    coupled_variables      = 'c1 c2 c3 eta1 eta3'
  [../]
  [./mult_ACBulkC_2]
    type = KKSMultiACBulkC
    variable  = lambda
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    cj_names  = 'c1 c2 c3'
    eta_i     = eta2
    coupled_variables      = 'eta1 eta3'
    mob_name  = 1
  [../]
  [./mult_CoupledACint_2]
    type = SimpleCoupledACInterface
    variable = lambda
    v = eta2
    kappa_name = kappa
    mob_name = 1
  [../]
  [./mult_ACBulkF_3]
    type = KKSMultiACBulkF
    variable  = lambda
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    gi_name   = g3
    eta_i     = eta3
    wi        = 1.0
    mob_name  = 1
    coupled_variables      = 'c1 c2 c3 eta1 eta2'
  [../]
  [./mult_ACBulkC_3]
    type = KKSMultiACBulkC
    variable  = lambda
    Fj_names  = 'F1 F2 F3'
    hj_names  = 'h1 h2 h3'
    cj_names  = 'c1 c2 c3'
    eta_i     = eta3
    coupled_variables      = 'eta1 eta2'
    mob_name  = 1
  [../]
  [./mult_CoupledACint_3]
    type = SimpleCoupledACInterface
    variable = lambda
    v = eta3
    kappa_name = kappa
    mob_name = 1
  [../]
  # Kernels for constraint equation eta1 + eta2 + eta3 = 1
  # eta3 is the nonlinear variable for the constraint equation
  [./eta3reaction]
    type = MatReaction
    variable = eta3
    reaction_rate = 1
  [../]
  [./eta1reaction]
    type = MatReaction
    variable = eta3
    v = eta1
    reaction_rate = 1
  [../]
  [./eta2reaction]
    type = MatReaction
    variable = eta3
    v = eta2
    reaction_rate = 1
  [../]
  [./one]
    type = BodyForce
    variable = eta3
    value = -1.0
  [../]
  # Phase concentration constraints
  [./chempot12]
    type = KKSPhaseChemicalPotential
    variable = c1
    cb       = c2
    fa_name  = F1
    fb_name  = F2
  [../]
  [./chempot23]
    type = KKSPhaseChemicalPotential
    variable = c2
    cb       = c3
    fa_name  = F2
    fb_name  = F3
  [../]
  [./phaseconcentration]
    type = KKSMultiPhaseConcentration
    variable = c3
    cj = 'c1 c2 c3'
    hj_names = 'h1 h2 h3'
    etas = 'eta1 eta2 eta3'
    c = c
  [../]
[]
[AuxKernels]
  [./Energy_total]
    type = KKSMultiFreeEnergy
    Fj_names = 'F1 F2 F3'
    hj_names = 'h1 h2 h3'
    gj_names = 'g1 g2 g3'
    variable = Energy
    w = 1
    interfacial_vars =  'eta1  eta2  eta3'
    kappa_names =       'kappa kappa kappa'
  [../]
[]
[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-10
  nl_abs_tol = 1.0e-11
  num_steps = 2
  dt = 0.5
[]
[Preconditioning]
  active = 'full'
  [./full]
    type = SMP
    full = true
  [../]
  [./mydebug]
    type = FDP
    full = true
  [../]
[]
[Outputs]
  exodus = true
[]
(modules/phase_field/test/tests/electrochem_sintering/ElectrochemicalSintering_test.i)
[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 800
  xmin = 0
  xmax = 80
[]
[GlobalParams]
  op_num = 2
  var_name_base = gr
  int_width = 4
[]
[Variables]
  [wvy]
  []
  [wvo]
  []
  [phi]
  []
  [PolycrystalVariables]
  []
  [V]
  []
[]
[AuxVariables]
  [bnds]
  []
  [negative_V]
  []
  [E_x]
    order = CONSTANT
    family = MONOMIAL
  []
  [E_y]
    order = CONSTANT
    family = MONOMIAL
  []
  [ns_cat_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [ns_an_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [T]
  []
[]
[Functions]
  [ic_func_gr0]
    type = ParsedFunction
    expression = '0.5*(1.0-tanh((x)/sqrt(2.0*2.0)))'
  []
  [ic_func_gr1]
    type = ParsedFunction
    expression = '0.5*(1.0+tanh((x)/sqrt(2.0*2.0)))'
  []
[]
[ICs]
  [gr0_IC]
    type = FunctionIC
    variable = gr0
    function = ic_func_gr0
  []
  [gr1_IC]
    type = FunctionIC
    variable = gr1
    function = ic_func_gr1
  []
  [wvy_IC]
    type = ConstantIC
    variable = wvy
    value = 2.7827
  []
  [wvo_IC]
    type = ConstantIC
    variable = wvo
    value = 2.7827
  []
  [T_IC]
    type = ConstantIC
    variable = T
    value = 1600
  []
[]
[BCs]
  [v_left]
    type = DirichletBC
    preset = true
    variable = V
    boundary = left
    value = 1e-2
  []
  [v_right]
    type = DirichletBC
    preset = true
    variable = V
    boundary = right
    value = 0
  []
  [gr0_left]
    type = DirichletBC
    preset = true
    variable = gr0
    boundary = left
    value = 0.5 #Grain boundary at left hand side of domain
  []
  [gr1_left]
    type = DirichletBC
    preset = true
    variable = gr1
    boundary = left
    value = 0.5 #Grain boundary at left hand side of domain
  []
  [wvo_right]
    type = DirichletBC
    preset = true
    variable = wvo
    boundary = right
    value = 2.7827
  []
  [wvy_right]
    type = DirichletBC
    preset = true
    variable = wvy
    boundary = right
    value = 2.7827
  []
[]
[Materials]
  # Free energy coefficients for parabolic curves
  [ks_cat]
    type = ParsedMaterial
    property_name = ks_cat
    coupled_variables = 'T'
    constant_names = 'a b Va'
    constant_expressions = '-0.0017 140.44 0.03726'
    expression = '(a*T + b) * Va^2'
  []
  [ks_an]
    type = ParsedMaterial
    property_name = ks_an
    coupled_variables = 'T'
    constant_names = 'a b Va'
    constant_expressions = '-0.0017 140.44 0.03726'
    expression = '(a*T + b) * Va^2'
  []
  [kv_cat]
    type = ParsedMaterial
    property_name = kv_cat
    material_property_names = 'ks_cat'
    expression = '10*ks_cat'
  []
  [kv_an]
    type = ParsedMaterial
    property_name = kv_an
    material_property_names = 'ks_cat'
    expression = '10*ks_cat'
  []
  # Diffusivity and mobilities
  [chiDy]
    type = GrandPotentialTensorMaterial
    f_name = chiDy
    diffusivity_name = Dvy
    solid_mobility = L
    void_mobility = Lv
    chi = chi_cat
    surface_energy = 6.24
    c = phi
    T = T
    D0 = 5.9e11
    GBmob0 = 1.60e12
    Q = 4.14
    Em = 4.25
    bulkindex = 1
    gbindex = 1
    surfindex = 1
  []
  [chiDo]
    type = GrandPotentialTensorMaterial
    f_name = chiDo
    diffusivity_name = Dvo
    solid_mobility = Lo
    void_mobility = Lvo
    chi = chi_an
    surface_energy = 6.24
    c = phi
    T = T
    D0 = 5.9e11
    GBmob0 = 1.60e12
    Q = 4.14
    Em = 4.25
    bulkindex = 1
    gbindex = 1
    surfindex = 1
  []
  # Everything else
  [ns_y_min]
    type = DerivativeParsedMaterial
    property_name = ns_y_min
    coupled_variables = 'gr0 gr1 T'
    constant_names = 'Ef_B Ef_GB   kB          Va_Y'
    constant_expressions = '4.37 4.37    8.617343e-5 0.03726'
    derivative_order = 2
    expression = 'bnds:=gr0^2 + gr1^2; Ef:=Ef_B + 4.0 * (Ef_GB - Ef_B) * (1.0 - bnds)^2;
              '
               '  exp(-Ef/kB/T) / Va_Y'
  []
  [ns_o_min]
    type = DerivativeParsedMaterial
    property_name = ns_o_min
    coupled_variables = 'gr0 gr1 T'
    constant_names = 'Ef_B Ef_GB  kB          Va_O'
    constant_expressions = '4.37 4.37   8.617343e-5 0.02484'
    derivative_order = 2
    expression = 'bnds:=gr0^2 + gr1^2; Ef:=Ef_B + 4.0 * (Ef_GB - Ef_B) * (1.0 - bnds)^2;
              '
               '  exp(-Ef/kB/T) / Va_O'
  []
  [sintering]
    type = ElectrochemicalSinteringMaterial
    chemical_potentials = 'wvy wvo'
    electric_potential = V
    void_op = phi
    Temperature = T
    surface_energy = 6.24
    grainboundary_energy = 5.18
    solid_energy_coefficients = 'kv_cat kv_cat'
    void_energy_coefficients = 'kv_cat kv_an'
    min_vacancy_concentrations_solid = 'ns_y_min ns_o_min'
    min_vacancy_concentrations_void = '26.837 40.256'
    defect_charges = '-3 2'
    solid_relative_permittivity = 30
    solid_energy_model = DILUTE
  []
  [density_chi_y]
    type = ElectrochemicalDefectMaterial
    chemical_potential = wvy
    void_op = phi
    Temperature = T
    electric_potential = V
    void_density_name = nv_cat
    solid_density_name = ns_cat
    chi_name = chi_cat
    void_energy_coefficient = kv_cat
    min_vacancy_concentration_solid = ns_y_min
    min_vacancy_concentration_void = 26.837
    solid_energy_model = DILUTE
    defect_charge = -3
    solid_relative_permittivity = 30
  []
  [density_chi_o]
    type = ElectrochemicalDefectMaterial
    chemical_potential = wvo
    void_op = phi
    Temperature = T
    electric_potential = V
    void_density_name = nv_an
    solid_density_name = ns_an
    chi_name = chi_an
    void_energy_coefficient = kv_an
    min_vacancy_concentration_solid = ns_o_min
    min_vacancy_concentration_void = 40.256
    solid_energy_model = DILUTE
    defect_charge = 2
    solid_relative_permittivity = 30
  []
  [permittivity]
    type = DerivativeParsedMaterial
    property_name = permittivity
    coupled_variables = 'phi'
    material_property_names = 'hs hv'
    constant_names = 'eps_rel_solid   eps_void_over_e'
    constant_expressions = '30              5.52e-2' #eps_void_over_e in 1/V/nm
    derivative_order = 2
    expression = '-hs * eps_rel_solid * eps_void_over_e - hv * eps_void_over_e'
  []
  [void_pre]
    type = DerivativeParsedMaterial
    property_name = void_pre
    material_property_names = 'hv'
    constant_names = 'Z_cat   Z_an nv_y_min nv_o_min'
    constant_expressions = '-3      2    26.837   40.256'
    derivative_order = 2
    expression = '-hv * (Z_cat * nv_y_min + Z_an * nv_o_min)'
  []
  [cat_mu_pre]
    type = DerivativeParsedMaterial
    property_name = cat_mu_pre
    material_property_names = 'hv kv_cat'
    constant_names = 'Z_cat'
    constant_expressions = '-3'
    derivative_order = 2
    expression = '-hv * Z_cat / kv_cat'
  []
  [an_mu_pre]
    type = DerivativeParsedMaterial
    property_name = an_mu_pre
    material_property_names = 'hv kv_an'
    constant_names = 'Z_an'
    constant_expressions = '2'
    derivative_order = 2
    expression = '-hv * Z_an / kv_an'
  []
  [cat_V_pre]
    type = DerivativeParsedMaterial
    property_name = cat_V_pre
    material_property_names = 'hv kv_cat'
    constant_names = 'Z_cat   v_scale e '
    constant_expressions = '-3      1       1'
    derivative_order = 2
    expression = 'hv * Z_cat^2 * e * v_scale / kv_cat'
  []
  [an_V_pre]
    type = DerivativeParsedMaterial
    property_name = an_V_pre
    material_property_names = 'hv kv_an'
    constant_names = 'Z_an    v_scale e '
    constant_expressions = '2       1       1'
    derivative_order = 2
    expression = 'hv * Z_an^2 * e * v_scale / kv_an'
  []
[]
#This action adds most kernels needed for grand potential model
[Modules]
  [PhaseField]
    [GrandPotential]
      switching_function_names = 'hv hs'
      anisotropic = 'true true'
      chemical_potentials = 'wvy wvo'
      mobilities = 'chiDy chiDo'
      susceptibilities = 'chi_cat chi_an'
      free_energies_w = 'nv_cat ns_cat nv_an ns_an'
      gamma_gr = gamma
      mobility_name_gr = L
      kappa_gr = kappa
      free_energies_gr = 'omegav omegas'
      additional_ops = 'phi'
      gamma_grxop = gamma
      mobility_name_op = Lv
      kappa_op = kappa
      free_energies_op = 'omegav omegas'
    []
  []
[]
[Kernels]
  [barrier_phi]
    type = ACBarrierFunction
    variable = phi
    v = 'gr0 gr1'
    gamma = gamma
    mob_name = Lv
  []
  [kappa_phi]
    type = ACKappaFunction
    variable = phi
    mob_name = Lv
    kappa_name = kappa
  []
  [Laplace]
    type = MatDiffusion
    variable = V
    diffusivity = permittivity
    args = 'phi'
  []
  [potential_void_constants]
    type = MaskedBodyForce
    variable = V
    coupled_variables = 'phi'
    mask = void_pre
  []
  [potential_cat_mu]
    type = MatReaction
    variable = V
    v = wvy
    reaction_rate = cat_mu_pre
  []
  [potential_an_mu]
    type = MatReaction
    variable = V
    v = wvo
    reaction_rate = an_mu_pre
  []
  [potential_cat_V]
    type = MatReaction
    variable = V
    reaction_rate = cat_V_pre
  []
  [potential_an_V]
    type = MatReaction
    variable = V
    reaction_rate = an_V_pre
  []
  [potential_solid_cat]
    type = MaskedExponential
    variable = V
    w = wvy
    T = T
    coupled_variables = 'phi gr0 gr1'
    mask = hs
    species_charge = -3
    n_eq = ns_y_min
  []
  [potential_solid_an]
    type = MaskedExponential
    variable = V
    w = wvo
    T = T
    coupled_variables = 'phi gr0 gr1'
    mask = hs
    species_charge = 2
    n_eq = ns_o_min
  []
[]
[AuxKernels]
  [bnds_aux]
    type = BndsCalcAux
    variable = bnds
    execute_on = 'initial timestep_end'
  []
  [negative_V]
    type = ParsedAux
    variable = negative_V
    coupled_variables = V
    expression = '-V'
  []
  [E_x]
    type = VariableGradientComponent
    variable = E_x
    gradient_variable = negative_V
    component = x
  []
  [E_y]
    type = VariableGradientComponent
    variable = E_y
    gradient_variable = negative_V
    component = y
  []
  [ns_cat_aux]
    type = MaterialRealAux
    variable = ns_cat_aux
    property = ns_cat
  []
  [ns_an_aux]
    type = MaterialRealAux
    variable = ns_an_aux
    property = ns_an
  []
[]
[Postprocessors]
  [ns_cat_total]
    type = ElementIntegralMaterialProperty
    mat_prop = ns_cat
  []
  [ns_an_total]
    type = ElementIntegralMaterialProperty
    mat_prop = ns_an
  []
[]
[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]
[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = PJFNK
  petsc_options_iname = '-pc_type -sub_pc_type -pc_asm_overlap -ksp_gmres_restart -sub_ksp_type'
  petsc_options_value = ' asm      lu           1               31                 preonly'
  nl_max_its = 40
  l_max_its = 30
  l_tol = 1e-4
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-13
  start_time = 0
  num_steps = 2
  automatic_scaling = true
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1
    optimal_iterations = 8
    iteration_window = 2
  []
[]
[Outputs]
  exodus = true
[]
(test/tests/neml2/mesh_change.i)
[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 10
  []
  [A]
    type = SubdomainBoundingBoxGenerator
    input = gmg
    block_id = 1
    block_name = A
    bottom_left = '0 0 0'
    top_right = '0.5 1 0'
  []
  [B]
    type = SubdomainBoundingBoxGenerator
    input = A
    block_id = 2
    block_name = B
    bottom_left = '0.5 0 0'
    top_right = '1 1 0'
  []
[]
[MeshModifiers]
  [AB]
    type = CoupledVarThresholdElementSubdomainModifier
    coupled_var = u
    criterion_type = ABOVE
    threshold = 0.5
    subdomain_id = 2
    complement_subdomain_id = 1
    execute_on = TIMESTEP_BEGIN
  []
[]
[NEML2]
  input = 'models/custom_model.i'
  verbose = true
  device = 'cpu'
  [A]
    model = 'model_A_non_ad'
    block = 'A'
    moose_input_types = 'VARIABLE MATERIAL'
    moose_inputs = '     a        b'
    neml2_inputs = '     forces/A forces/B'
    moose_output_types = 'MATERIAL           MATERIAL'
    moose_outputs = '     neml2_sum          neml2_product'
    neml2_outputs = '     state/internal/sum state/internal/product'
    moose_derivative_types = 'MATERIAL'
    moose_derivatives = 'neml2_dproduct_da'
    neml2_derivatives = 'state/internal/product forces/A'
    export_outputs = 'neml2_sum neml2_product neml2_dproduct_da'
    export_output_targets = 'exodus; exodus; exodus'
  []
  [B]
    model = 'model_B_non_ad'
    block = 'B'
    moose_input_types = 'VARIABLE MATERIAL'
    moose_inputs = '     a        b'
    neml2_inputs = '     forces/A forces/B'
    moose_output_types = 'MATERIAL           MATERIAL'
    moose_outputs = '     neml2_sum          neml2_product'
    neml2_outputs = '     state/internal/sum state/internal/product'
    moose_derivative_types = 'MATERIAL'
    moose_derivatives = 'neml2_dproduct_da'
    neml2_derivatives = 'state/internal/product forces/A'
    export_outputs = 'neml2_sum neml2_product neml2_dproduct_da'
    export_output_targets = 'exodus; exodus; exodus'
  []
[]
[Variables]
  [u]
    [InitialCondition]
      type = FunctionIC
      function = 'x'
    []
  []
[]
[Kernels]
  [rate]
    type = TimeDerivative
    variable = u
  []
  [diffusion]
    type = Diffusion
    variable = u
  []
  [reaction_1]
    type = MatReaction
    variable = u
    reaction_rate = neml2_sum
  []
  [reaction_2]
    type = MatReaction
    variable = u
    reaction_rate = neml2_product
  []
  [reaction_3]
    type = MatReaction
    variable = u
    reaction_rate = neml2_dproduct_da
  []
[]
[AuxVariables]
  [a]
  []
[]
[ICs]
  [a]
    type = FunctionIC
    variable = a
    function = 'x'
  []
[]
[Materials]
  [b]
    type = GenericFunctionMaterial
    prop_names = 'b'
    prop_values = 'y+t'
    outputs = 'exodus'
  []
[]
[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  num_steps = 5
  dt = 0.1
[]
[Outputs]
  exodus = true
[]
(modules/fsi/test/tests/2d-small-strain-transient/fsi_flat_channel.i)
[GlobalParams]
  gravity = '0 0 0'
  integrate_p_by_parts = true
  laplace = true
  convective_term = true
  transient_term = true
  pspg = true
  supg = true
  displacements = 'disp_x disp_y'
  preset = false
  order = FIRST
  use_displaced_mesh = true
[]
[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = 3.0
    ymin = 0
    ymax = 1.0
    nx = 10
    ny = 15
    elem_type = QUAD4
  []
  [subdomain1]
    type = SubdomainBoundingBoxGenerator
    bottom_left = '0.0 0.5 0'
    block_id = 1
    top_right = '3.0 1.0 0'
    input = gmg
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '0'
    paired_block = '1'
    new_boundary = 'master0_interface'
    input = subdomain1
  []
  [break_boundary]
    type = BreakBoundaryOnSubdomainGenerator
    input = interface
  []
[]
[Variables]
  [./vel_x]
    block = 0
  [../]
  [./vel_y]
    block = 0
  [../]
  [./p]
    block = 0
  [../]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./vel_x_solid]
    block = 1
  [../]
  [./vel_y_solid]
    block = 1
  [../]
[]
[Kernels]
  [./vel_x_time]
    type = INSMomentumTimeDerivative
    variable = vel_x
    block = 0
  [../]
  [./vel_y_time]
    type = INSMomentumTimeDerivative
    variable = vel_y
    block = 0
  [../]
  [./mass]
    type = INSMass
    variable = p
    u = vel_x
    v = vel_y
    pressure = p
    block = 0
    disp_x = disp_x
    disp_y = disp_y
  [../]
  [./x_momentum_space]
    type = INSMomentumLaplaceForm
    variable = vel_x
    u = vel_x
    v = vel_y
    pressure = p
    component = 0
    block = 0
    disp_x = disp_x
    disp_y = disp_y
  [../]
  [./y_momentum_space]
    type = INSMomentumLaplaceForm
    variable = vel_y
    u = vel_x
    v = vel_y
    pressure = p
    component = 1
    block = 0
    disp_x = disp_x
    disp_y = disp_y
  [../]
  [./vel_x_mesh]
    type = ConvectedMesh
    disp_x = disp_x
    disp_y = disp_y
    variable = vel_x
    u = vel_x
    v = vel_y
    pressure = p
    block = 0
  [../]
  [./vel_y_mesh]
    type = ConvectedMesh
    disp_x = disp_x
    disp_y = disp_y
    variable = vel_y
    u = vel_x
    v = vel_y
    pressure = p
    block = 0
  [../]
  [./p_mesh]
    type = ConvectedMeshPSPG
    disp_x = disp_x
    disp_y = disp_y
    variable = p
    u = vel_x
    v = vel_y
    pressure = p
    block = 0
  [../]
  [./disp_x_fluid]
    type = Diffusion
    variable = disp_x
    block = 0
    use_displaced_mesh = false
  [../]
  [./disp_y_fluid]
    type = Diffusion
    variable = disp_y
    block = 0
    use_displaced_mesh = false
  [../]
  [./accel_tensor_x]
    type = CoupledTimeDerivative
    variable = disp_x
    v = vel_x_solid
    block = 1
    use_displaced_mesh = false
  [../]
  [./accel_tensor_y]
    type = CoupledTimeDerivative
    variable = disp_y
    v = vel_y_solid
    block = 1
    use_displaced_mesh = false
  [../]
  [./vxs_time_derivative_term]
    type = CoupledTimeDerivative
    variable = vel_x_solid
    v = disp_x
    block = 1
    use_displaced_mesh = false
  [../]
  [./vys_time_derivative_term]
    type = CoupledTimeDerivative
    variable = vel_y_solid
    v = disp_y
    block = 1
    use_displaced_mesh = false
  [../]
  [./source_vxs]
    type = MatReaction
    variable = vel_x_solid
    block = 1
    reaction_rate = 1
    use_displaced_mesh = false
  [../]
  [./source_vys]
    type = MatReaction
    variable = vel_y_solid
    block = 1
    reaction_rate = 1
    use_displaced_mesh = false
  [../]
[]
[InterfaceKernels]
  [./penalty_interface_x]
    type = CoupledPenaltyInterfaceDiffusion
    variable = vel_x
    neighbor_var = disp_x
    secondary_coupled_var = vel_x_solid
    boundary = master0_interface
    penalty = 1e6
  [../]
  [./penalty_interface_y]
    type = CoupledPenaltyInterfaceDiffusion
    variable = vel_y
    neighbor_var = disp_y
    secondary_coupled_var = vel_y_solid
    boundary = master0_interface
    penalty = 1e6
  [../]
[]
[Physics/SolidMechanics/QuasiStatic]
  [./solid_domain]
    strain = SMALL
    incremental = false
    # generate_output = 'strain_xx strain_yy strain_zz' ## Not at all necessary, but nice
    block = '1'
  [../]
[]
[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e2
    poissons_ratio = 0.3
    block = '1'
    use_displaced_mesh = false
  [../]
  [./small_stress]
    type = ComputeLinearElasticStress
    block = 1
  [../]
  [./const]
    type = GenericConstantMaterial
    block = 0
    prop_names = 'rho mu'
    prop_values = '1  1'
    use_displaced_mesh = false
  [../]
[]
[BCs]
  [./fluid_x_no_slip]
    type = DirichletBC
    variable = vel_x
    boundary = 'bottom'
    value = 0.0
  [../]
  [./fluid_y_no_slip]
    type = DirichletBC
    variable = vel_y
    boundary = 'bottom left_to_0'
    value = 0.0
  [../]
  [./x_inlet]
    type = FunctionDirichletBC
    variable = vel_x
    boundary = 'left_to_0'
    function = 'inlet_func'
  [../]
  [./no_disp_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'bottom top left_to_1 right_to_1 left_to_0 right_to_0'
    value = 0
  [../]
  [./no_disp_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom top left_to_1 right_to_1 left_to_0 right_to_0'
    value = 0
  [../]
  [./solid_x_no_slip]
    type = DirichletBC
    variable = vel_x_solid
    boundary = 'top left_to_1 right_to_1'
    value = 0.0
  [../]
  [./solid_y_no_slip]
    type = DirichletBC
    variable = vel_y_solid
    boundary = 'top left_to_1 right_to_1'
    value = 0.0
  [../]
[]
[Preconditioning]
  [./SMP]
    type = FDP
    full = true
  [../]
[]
[Executioner]
  type = Transient
  num_steps = 5
  # num_steps = 60
  dt = 0.1
  dtmin = 0.1
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu       NONZERO'
  line_search = none
  nl_rel_tol = 1e-50
  nl_abs_tol = 1e-10
[]
[Outputs]
  [./out]
    type = Exodus
  [../]
[]
[Functions]
  [./inlet_func]
    type = ParsedFunction
    expression = '(-16 * (y - 0.25)^2 + 1) * (1 + cos(t))'
  [../]
[]
(test/tests/neml2/simple_scheduler.i)
[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 10
  []
[]
[NEML2]
  input = 'models/custom_model.i'
  scheduler = 'simple'
  async_dispatch = false
  [all]
    model = 'model_non_ad'
    verbose = true
    device = 'cpu'
    moose_input_types = 'VARIABLE MATERIAL'
    moose_inputs = '     a        b'
    neml2_inputs = '     forces/A forces/B'
    moose_output_types = 'MATERIAL           MATERIAL'
    moose_outputs = '     neml2_sum          neml2_product'
    neml2_outputs = '     state/internal/sum state/internal/product'
    moose_derivative_types = 'MATERIAL'
    moose_derivatives = 'neml2_dproduct_da'
    neml2_derivatives = 'state/internal/product forces/A'
    export_outputs = 'neml2_sum neml2_product neml2_dproduct_da'
    export_output_targets = 'exodus; exodus; exodus'
  []
[]
[Variables]
  [u]
  []
[]
[Kernels]
  [diffusion]
    type = Diffusion
    variable = u
  []
  [reaction_1]
    type = MatReaction
    variable = u
    reaction_rate = neml2_sum
  []
  [reaction_2]
    type = MatReaction
    variable = u
    reaction_rate = neml2_product
  []
  [reaction_3]
    type = MatReaction
    variable = u
    reaction_rate = neml2_dproduct_da
  []
[]
[AuxVariables]
  [a]
  []
[]
[ICs]
  [a]
    type = FunctionIC
    variable = a
    function = 'x'
  []
[]
[Materials]
  [b]
    type = GenericFunctionMaterial
    prop_names = 'b'
    prop_values = 'y+t'
    outputs = 'exodus'
  []
[]
[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  num_steps = 2
[]
[Outputs]
  exodus = true
[]
(modules/phase_field/examples/multiphase/GrandPotential3Phase_masscons.i)
# This is an example of implementation of the multi-phase, multi-order parameter
# grand potential based phase-field model described in Phys. Rev. E, 98, 023309
# (2018). It includes 3 phases with 1 grain of each phase.
# This is a revised version of the model that eliminates small variations in mass
# that have been observed with the original formulation. In this version, rather
# than evolving the chemical potential as a field variable, we evolve the composition
# field using a normal Cahn-Hilliard equation, then relate chemical potential to
# composition using Eq. (22) from the paper (this relationship is derived from the
# grand potential functional and is valid only for parabolic free energies).
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 60
  ny = 60
  xmin = -15
  xmax = 15
  ymin = -15
  ymax = 15
[]
[Variables]
  [w]
  []
  [c]
  []
  [etaa0]
  []
  [etab0]
  []
  [etad0]
  []
[]
[ICs]
  [IC_etaa0]
    type = BoundingBoxIC
    variable = etaa0
    x1 = -10
    y1 = -10
    x2 = 10
    y2 = 10
    inside = 1.0
    outside = 0.0
  []
  [IC_etad0]
    type = BoundingBoxIC
    variable = etad0
    x1 = -10
    y1 = -10
    x2 = 10
    y2 = 10
    inside = 0.0
    outside = 1.0
  []
  [IC_c]
    type = BoundingBoxIC
    variable = c
    x1 = -10
    y1 = -10
    x2 = 10
    y2 = 10
    inside = 0.1
    outside = 0.5
  []
  [IC_w]
    type = FunctionIC
    variable = w
    function = ic_func_w
  []
[]
[Functions]
  [ic_func_w]
    type = ConstantFunction
    value = 0
  []
[]
[Kernels]
  # Order parameter eta_alpha0
  [ACa0_bulk]
    type = ACGrGrMulti
    variable = etaa0
    v = 'etab0 etad0'
    gamma_names = 'gab   gad'
  []
  [ACa0_sw]
    type = ACSwitching
    variable = etaa0
    Fj_names = 'omegaa omegab omegad'
    hj_names = 'ha     hb     hd'
    coupled_variables = 'etab0 etad0 w'
  []
  [ACa0_int]
    type = ACInterface
    variable = etaa0
    kappa_name = kappa
  []
  [ea0_dot]
    type = TimeDerivative
    variable = etaa0
  []
  # Order parameter eta_beta0
  [ACb0_bulk]
    type = ACGrGrMulti
    variable = etab0
    v = 'etaa0 etad0'
    gamma_names = 'gab   gbd'
  []
  [ACb0_sw]
    type = ACSwitching
    variable = etab0
    Fj_names = 'omegaa omegab omegad'
    hj_names = 'ha     hb     hd'
    coupled_variables = 'etaa0 etad0 w'
  []
  [ACb0_int]
    type = ACInterface
    variable = etab0
    kappa_name = kappa
  []
  [eb0_dot]
    type = TimeDerivative
    variable = etab0
  []
  # Order parameter eta_delta0
  [ACd0_bulk]
    type = ACGrGrMulti
    variable = etad0
    v = 'etaa0 etab0'
    gamma_names = 'gad   gbd'
  []
  [ACd0_sw]
    type = ACSwitching
    variable = etad0
    Fj_names = 'omegaa omegab omegad'
    hj_names = 'ha     hb     hd'
    coupled_variables = 'etaa0 etab0 w'
  []
  [ACd0_int]
    type = ACInterface
    variable = etad0
    kappa_name = kappa
  []
  [ed0_dot]
    type = TimeDerivative
    variable = etad0
  []
  #Concentration
  [c_dot]
    type = TimeDerivative
    variable = c
  []
  [Diffusion]
    type = MatDiffusion
    variable = c
    v = w
    diffusivity = DchiVm
    args = ''
  []
  #The following relate chemical potential to composition using Eq. (22)
  [w_rxn]
    type = MatReaction
    variable = w
    v = c
    reaction_rate = -1
  []
  [ca_rxn]
    type = MatReaction
    variable = w
    reaction_rate = 'hoverk_a'
    args = 'etaa0 etab0 etad0'
  []
  [ca_bodyforce]
    type = MaskedBodyForce
    variable = w
    mask = ha
    coupled_variables = 'etaa0 etab0 etad0'
    value = 0.1 #caeq
  []
  [cb_rxn]
    type = MatReaction
    variable = w
    reaction_rate = 'hoverk_b'
    args = 'etaa0 etab0 etad0'
  []
  [cb_bodyforce]
    type = MaskedBodyForce
    variable = w
    mask = hb
    coupled_variables = 'etaa0 etab0 etad0'
    value = 0.9 #cbeq
  []
  [cd_rxn]
    type = MatReaction
    variable = w
    reaction_rate = 'hoverk_d'
    args = 'etaa0 etab0 etad0'
  []
  [cd_bodyforce]
    type = MaskedBodyForce
    variable = w
    mask = hd
    coupled_variables = 'etaa0 etab0 etad0'
    value = 0.5 #cdeq
  []
[]
[Materials]
  [ha_test]
    type = SwitchingFunctionMultiPhaseMaterial
    h_name = ha
    all_etas = 'etaa0 etab0 etad0'
    phase_etas = 'etaa0'
  []
  [hb_test]
    type = SwitchingFunctionMultiPhaseMaterial
    h_name = hb
    all_etas = 'etaa0 etab0 etad0'
    phase_etas = 'etab0'
  []
  [hd_test]
    type = SwitchingFunctionMultiPhaseMaterial
    h_name = hd
    all_etas = 'etaa0 etab0 etad0'
    phase_etas = 'etad0'
  []
  [omegaa]
    type = DerivativeParsedMaterial
    coupled_variables = 'w'
    property_name = omegaa
    material_property_names = 'Vm ka caeq'
    expression = '-0.5*w^2/Vm^2/ka-w/Vm*caeq'
    derivative_order = 2
  []
  [omegab]
    type = DerivativeParsedMaterial
    coupled_variables = 'w'
    property_name = omegab
    material_property_names = 'Vm kb cbeq'
    expression = '-0.5*w^2/Vm^2/kb-w/Vm*cbeq'
    derivative_order = 2
  []
  [omegad]
    type = DerivativeParsedMaterial
    coupled_variables = 'w'
    property_name = omegad
    material_property_names = 'Vm kd cdeq'
    expression = '-0.5*w^2/Vm^2/kd-w/Vm*cdeq'
    derivative_order = 2
  []
  [rhoa]
    type = DerivativeParsedMaterial
    coupled_variables = 'w'
    property_name = rhoa
    material_property_names = 'Vm ka caeq'
    expression = 'w/Vm^2/ka + caeq/Vm'
    derivative_order = 2
  []
  [rhob]
    type = DerivativeParsedMaterial
    coupled_variables = 'w'
    property_name = rhob
    material_property_names = 'Vm kb cbeq'
    expression = 'w/Vm^2/kb + cbeq/Vm'
    derivative_order = 2
  []
  [rhod]
    type = DerivativeParsedMaterial
    coupled_variables = 'w'
    property_name = rhod
    material_property_names = 'Vm kd cdeq'
    expression = 'w/Vm^2/kd + cdeq/Vm'
    derivative_order = 2
  []
  [const]
    type = GenericConstantMaterial
    prop_names = 'kappa_c  kappa   L   D    Vm   ka    caeq kb    cbeq  kd    cdeq  gab gad gbd  mu  tgrad_corr_mult'
    prop_values = '0        1       1.0 1.0  1.0  10.0  0.1  10.0  0.9   10.0  0.5   1.5 1.5 1.5  1.0 0.0'
  []
  [Mobility]
    type = DerivativeParsedMaterial
    property_name = DchiVm
    material_property_names = 'D chi Vm' #Factor of Vm is needed to evolve c instead of rho
    expression = 'D*chi*Vm'
    derivative_order = 2
  []
  [chi]
    type = DerivativeParsedMaterial
    property_name = chi
    material_property_names = 'Vm ha(etaa0,etab0,etad0) ka hb(etaa0,etab0,etad0) kb hd(etaa0,etab0,etad0) kd'
    expression = '(ha/ka + hb/kb + hd/kd) / Vm^2'
    coupled_variables = 'etaa0 etab0 etad0'
    derivative_order = 2
  []
  [hoverk_a]
    type = DerivativeParsedMaterial
    material_property_names = 'ha(etaa0,etab0,etad0) Vm ka'
    property_name = hoverk_a
    expression = 'ha / Vm / ka'
  []
  [hoverk_b]
    type = DerivativeParsedMaterial
    material_property_names = 'hb(etaa0,etab0,etad0) Vm kb'
    property_name = hoverk_b
    expression = 'hb / Vm / kb'
  []
  [hoverk_d]
    type = DerivativeParsedMaterial
    material_property_names = 'hd(etaa0,etab0,etad0) Vm kd'
    property_name = hoverk_d
    expression = 'hd / Vm / kd'
  []
[]
[Postprocessors]
  [c_total]
    type = ElementIntegralVariablePostprocessor
    variable = c
  []
[]
[Executioner]
  type = Transient
  nl_max_its = 15
  scheme = bdf2
  solve_type = NEWTON
  petsc_options_iname = -pc_type
  petsc_options_value = asm
  l_max_its = 15
  l_tol = 1.0e-3
  nl_rel_tol = 1.0e-8
  start_time = 0.0
  num_steps = 20
  nl_abs_tol = 1e-10
  dt = 1.0
[]
[Outputs]
  csv = true
  exodus = true
[]
(modules/phase_field/test/tests/mobility_derivative/coupledmatdiffusion.i)
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 15
  ny = 15
  xmax = 15.0
  ymax = 15.0
  elem_type = QUAD4
[]
[Variables]
  [./c]
    [./InitialCondition]
      type = CrossIC
      x1 = 0.0
      x2 = 30.0
      y1 = 0.0
      y2 = 30.0
    [../]
  [../]
  [./d]
    [./InitialCondition]
      type = SmoothCircleIC
      x1 = 15
      y1 = 15
      radius = 8
      int_width = 3
      invalue = 2
      outvalue = 0
    [../]
  [../]
  [./u]
  [../]
  [./w]
  [../]
[]
[Kernels]
  [./ctime]
    type = TimeDerivative
    variable = c
  [../]
  [./umat]
    type = MatReaction
    variable = c
    v = u
    reaction_rate = 1
  [../]
  [./urxn]
    type = Reaction
    variable = u
  [../]
  [./cres]
    type = MatDiffusion
    variable = u
    diffusivity = Dc
    args = d
    v = c
  [../]
  [./dtime]
    type = TimeDerivative
    variable = d
  [../]
  [./wmat]
    type = MatReaction
    variable = d
    v = w
    reaction_rate = 1
  [../]
  [./wrxn]
    type = Reaction
    variable = w
  [../]
  [./dres]
    type = MatDiffusion
    variable = w
    diffusivity = Dd
    args = c
    v = d
  [../]
[]
[Materials]
  [./Dc]
    type = DerivativeParsedMaterial
    property_name = Dc
    expression = '0.01+c^2+d'
    coupled_variables = 'c d'
    derivative_order = 1
  [../]
  [./Dd]
    type = DerivativeParsedMaterial
    property_name = Dd
    expression = 'd^2+c+1.5'
    coupled_variables = 'c d'
    derivative_order = 1
  [../]
[]
[Preconditioning]
  [./SMP]
    type = SMP
    full = true
  [../]
[]
[Executioner]
  type = Transient
  scheme = 'BDF2'
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -ksp_gmres_restart -sub_pc_type -pc_asm_overlap'
  petsc_options_value = 'asm      31                  lu           1'
  dt = 1
  num_steps = 2
[]
[Outputs]
  exodus = true
[]
(test/tests/neml2/custom_model.i)
[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 10
  []
[]
[NEML2]
  input = 'models/custom_model.i'
  [all]
    model = 'model'
    verbose = true
    device = 'cpu'
    moose_input_types = 'VARIABLE MATERIAL'
    moose_inputs = '     a        b'
    neml2_inputs = '     forces/A forces/B'
    moose_output_types = 'MATERIAL           MATERIAL'
    moose_outputs = '     neml2_sum          neml2_product'
    neml2_outputs = '     state/internal/sum state/internal/product'
    moose_derivative_types = 'MATERIAL'
    moose_derivatives = 'neml2_dproduct_da'
    neml2_derivatives = 'state/internal/product forces/A'
    export_outputs = 'neml2_sum neml2_product neml2_dproduct_da'
    export_output_targets = 'exodus; exodus; exodus'
  []
[]
[Variables]
  [u]
  []
[]
[Kernels]
  [diffusion]
    type = Diffusion
    variable = u
  []
  [reaction_1]
    type = MatReaction
    variable = u
    reaction_rate = neml2_sum
  []
  [reaction_2]
    type = MatReaction
    variable = u
    reaction_rate = neml2_product
  []
  [reaction_3]
    type = MatReaction
    variable = u
    reaction_rate = neml2_dproduct_da
  []
[]
[AuxVariables]
  [a]
  []
[]
[ICs]
  [a]
    type = FunctionIC
    variable = a
    function = 'x'
  []
[]
[Materials]
  [b]
    type = GenericFunctionMaterial
    prop_names = 'b'
    prop_values = 'y+t'
    outputs = 'exodus'
  []
[]
[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  num_steps = 2
[]
[Outputs]
  exodus = true
[]
(test/tests/neml2/simple_scheduler_async.i)
[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 10
  []
[]
[NEML2]
  input = 'models/custom_model.i'
  scheduler = 'simple'
  async_dispatch = true
  [all]
    model = 'model_non_ad'
    verbose = true
    device = 'cpu'
    moose_input_types = 'VARIABLE MATERIAL'
    moose_inputs = '     a        b'
    neml2_inputs = '     forces/A forces/B'
    moose_output_types = 'MATERIAL           MATERIAL'
    moose_outputs = '     neml2_sum          neml2_product'
    neml2_outputs = '     state/internal/sum state/internal/product'
    moose_derivative_types = 'MATERIAL'
    moose_derivatives = 'neml2_dproduct_da'
    neml2_derivatives = 'state/internal/product forces/A'
    export_outputs = 'neml2_sum neml2_product neml2_dproduct_da'
    export_output_targets = 'exodus; exodus; exodus'
  []
[]
[Variables]
  [u]
  []
[]
[Kernels]
  [diffusion]
    type = Diffusion
    variable = u
  []
  [reaction_1]
    type = MatReaction
    variable = u
    reaction_rate = neml2_sum
  []
  [reaction_2]
    type = MatReaction
    variable = u
    reaction_rate = neml2_product
  []
  [reaction_3]
    type = MatReaction
    variable = u
    reaction_rate = neml2_dproduct_da
  []
[]
[AuxVariables]
  [a]
  []
[]
[ICs]
  [a]
    type = FunctionIC
    variable = a
    function = 'x'
  []
[]
[Materials]
  [b]
    type = GenericFunctionMaterial
    prop_names = 'b'
    prop_values = 'y+t'
    outputs = 'exodus'
  []
[]
[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  num_steps = 2
[]
[Outputs]
  exodus = true
[]
(test/include/kernels/RenamedMatReaction.h)
// This file is part of the MOOSE framework
// https://mooseframework.inl.gov
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html
#pragma once
#include "MatReaction.h"
class RenamedMatReaction : public MatReaction
{
public:
  static InputParameters validParams();
  RenamedMatReaction(const InputParameters & parameters);
};