Mass Sorption Constraint

Enforces sorption as described by a concentration-pressure proportionality between primary and secondary sides of a mortar interface using a Lagrange multiplier.

Description

The MassSorptionConstraint class is used to enforce sorption as described by a concentration-pressure proportionality between primary ( $var element = document.getElementById("moose-equation-ed881778-10b4-4525-9a33-4ec74ded7927");katex.render("1", element, {displayMode:false,throwOnError:false});$ ) and secondary ( $var element = document.getElementById("moose-equation-4a58b7ec-72c3-4d6e-8c8d-714b8b36e26c");katex.render("2", element, {displayMode:false,throwOnError:false});$ ) sides of a mortar interface using a Lagrange multiplier. It enforces the constraint $(1)var element = document.getElementById("moose-equation-241118fd-9a85-42c9-9581-28b2e056efb4");katex.render(" P_1 - P_2 = 0, ", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-c8ebd5a5-c11a-4e07-90fb-77367737c670");katex.render("P_i", element, {displayMode:false,throwOnError:false});$ is the partial pressure of a species in the gap that would result from concentration ( $var element = document.getElementById("moose-equation-88a7e10c-9c41-41bc-b3bc-33e86a303144");katex.render("c_i", element, {displayMode:false,throwOnError:false});$ ) and temperature ( $var element = document.getElementById("moose-equation-9f90d88b-aadb-4e54-bb7a-df709231d290");katex.render("T_i", element, {displayMode:false,throwOnError:false});$ ) conditions at interface side $var element = document.getElementById("moose-equation-3f1175ec-dff4-4092-b774-94a489e69b4b");katex.render("i", element, {displayMode:false,throwOnError:false});$ . Partial pressures can be calculated using SorptionPartialPressure.

A small non-zero shift factor, $var element = document.getElementById("moose-equation-c424ce4d-f9d2-4f70-b27d-b9c80e2691b7");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ , is applied to avoid a zero pivot when inverting the Jacobian. SorptionIsothermGapInterface implements similar physics across an interface that is not separated by a gap. Unlike SorptionIsothermGapInterface, MassSorptionConstraint does not necessarily preserve mass flux. This functionality is provided by pairing it with a class such as MassFluxConstraint.

Example Input Syntax

[Constraints<<<{"href": "../../syntax/Constraints/index.html"}>>>]
  [u]
    type = MassSorptionConstraint<<<{"description": "Enforces sorption as described by a concentration-pressure proportionality between primary and secondary sides of a mortar interface using a Lagrange multiplier.", "href": "MassSorptionConstraint.html"}>>>
    variable<<<{"description": "The name of the lagrange multiplier variable that this constraint is applied to. This parameter may not be supplied in the case of using penalty methods for example"}>>> = lm_u
    primary_variable<<<{"description": "Primal variable on primary surface. If this parameter is not provided then the primary variable will be initialized to the secondary variable"}>>> = u
    primary_boundary<<<{"description": "The name of the primary boundary sideset."}>>> = inner_right
    primary_subdomain<<<{"description": "The name of the primary subdomain."}>>> = primary
    secondary_variable<<<{"description": "Primal variable on secondary surface."}>>> = u
    secondary_boundary<<<{"description": "The name of the secondary boundary sideset."}>>> = outer_left
    secondary_subdomain<<<{"description": "The name of the secondary subdomain."}>>> = secondary
    partial_pressure_name<<<{"description": "The name of the partial pressure if multiple partial pressure are defined on the same block"}>>> = partial_pressure
    epsilon<<<{"description": "The non-zero shift factor applied to avoid a zero pivot"}>>> = 1e-9
    correct_edge_dropping<<<{"description": "Whether to enable correct edge dropping treatment for mortar constraints. When disabled any Lagrange Multiplier degree of freedom on a secondary element without full primary contributions will be set (strongly) to 0."}>>> = true
  []
[]

Here, the partial pressures on the two sides of the interface are calculated using SorptionPartialPressure as shown below.

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [partial_pressure_inner]
    type = ADSorptionPartialPressure<<<{"description": "Computes the sorption partial pressure of a solute in a material.", "href": "../materials/SorptionPartialPressure.html"}>>>
    A<<<{"description": "Temperature-independent coefficient of the intercept term in the concentration-pressure proportionality"}>>> = 19.3
    B<<<{"description": "Temperature-dependent coefficient of the intercept term in the concentration-pressure proportionality"}>>> = -47300
    D<<<{"description": "Temperature-independent coefficient of the order in the concentration-pressure proportionality"}>>> = 1.51
    E<<<{"description": "Temperature-dependent coefficient of the order in the concentration-pressure proportionality"}>>> = 4340
    d1<<<{"description": "Temperature-independent coefficient of the transition concentration"}>>> = 3.4
    d2<<<{"description": "Temperature-dependent coefficient of the transition concentration"}>>> = 6.15e-4
    unit_scale<<<{"description": "Unit conversion factor used to scale the concentration"}>>> = 1
    density<<<{"description": "Mass density of the material"}>>> = density
    concentration<<<{"description": "The coupled concentration variable"}>>> = u
    temperature<<<{"description": "The coupled temperature variable"}>>> = temperature
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = inner
    outputs<<<{"description": "Vector of output names where you would like to restrict the output of variables(s) associated with this object"}>>> = 'all'
    output_properties<<<{"description": "List of material properties, from this material, to output (outputs must also be defined to an output type)"}>>> = partial_pressure
  []
[]

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [partial_pressure_outer]
    type = ADSorptionPartialPressure<<<{"description": "Computes the sorption partial pressure of a solute in a material.", "href": "../materials/SorptionPartialPressure.html"}>>>
    A<<<{"description": "Temperature-independent coefficient of the intercept term in the concentration-pressure proportionality"}>>> = 24
    B<<<{"description": "Temperature-dependent coefficient of the intercept term in the concentration-pressure proportionality"}>>> = -35700
    D<<<{"description": "Temperature-independent coefficient of the order in the concentration-pressure proportionality"}>>> = -1.56
    E<<<{"description": "Temperature-dependent coefficient of the order in the concentration-pressure proportionality"}>>> = 6120
    d1<<<{"description": "Temperature-independent coefficient of the transition concentration"}>>> = 2.04
    d2<<<{"description": "Temperature-dependent coefficient of the transition concentration"}>>> = 1.79e-3
    unit_scale<<<{"description": "Unit conversion factor used to scale the concentration"}>>> = 1
    density<<<{"description": "Mass density of the material"}>>> = density
    concentration<<<{"description": "The coupled concentration variable"}>>> = u
    temperature<<<{"description": "The coupled temperature variable"}>>> = temperature
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = outer
    outputs<<<{"description": "Vector of output names where you would like to restrict the output of variables(s) associated with this object"}>>> = 'all'
    output_properties<<<{"description": "List of material properties, from this material, to output (outputs must also be defined to an output type)"}>>> = partial_pressure
  []
[]

Input Parameters

primary_boundaryThe name of the primary boundary sideset.
C++ Type:BoundaryName
Controllable:No
Description:The name of the primary boundary sideset.
primary_subdomainThe name of the primary subdomain.
C++ Type:SubdomainName
Controllable:No
Description:The name of the primary subdomain.
secondary_boundaryThe name of the secondary boundary sideset.
C++ Type:BoundaryName
Controllable:No
Description:The name of the secondary boundary sideset.
secondary_subdomainThe name of the secondary subdomain.
C++ Type:SubdomainName
Controllable:No
Description:The name of the secondary subdomain.

Required Parameters

aux_lmAuxiliary Lagrange multiplier variable that is utilized together with the Petrov-Galerkin approach.
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Auxiliary Lagrange multiplier variable that is utilized together with the Petrov-Galerkin approach.
compute_lm_residualsTrueWhether to compute Lagrange Multiplier residuals
Default:True
C++ Type:bool
Controllable:No
Description:Whether to compute Lagrange Multiplier residuals
compute_primal_residualsTrueWhether to compute residuals for the primal variable.
Default:True
C++ Type:bool
Controllable:No
Description:Whether to compute residuals for the primal variable.
correct_edge_droppingFalseWhether to enable correct edge dropping treatment for mortar constraints. When disabled any Lagrange Multiplier degree of freedom on a secondary element without full primary contributions will be set (strongly) to 0.
Default:False
C++ Type:bool
Controllable:No
Description:Whether to enable correct edge dropping treatment for mortar constraints. When disabled any Lagrange Multiplier degree of freedom on a secondary element without full primary contributions will be set (strongly) to 0.
debug_meshFalseWhether this constraint is going to enable mortar segment mesh debug information. An exodusfile will be generated if the user sets this flag to true
Default:False
C++ Type:bool
Controllable:No
Description:Whether this constraint is going to enable mortar segment mesh debug information. An exodusfile will be generated if the user sets this flag to true
epsilon1e-09The non-zero shift factor applied to avoid a zero pivot
Default:1e-09
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The non-zero shift factor applied to avoid a zero pivot
ghost_higher_d_neighborsFalseWhether we should ghost higher-dimensional neighbors. This is necessary when we are doing second order mortar with finite volume primal variables, because in order for the method to be second order we must use cell gradients, which couples in the neighbor cells.
Default:False
C++ Type:bool
Controllable:No
Description:Whether we should ghost higher-dimensional neighbors. This is necessary when we are doing second order mortar with finite volume primal variables, because in order for the method to be second order we must use cell gradients, which couples in the neighbor cells.
ghost_point_neighborsFalseWhether we should ghost point neighbors of secondary face elements, and consequently also their mortar interface couples.
Default:False
C++ Type:bool
Controllable:No
Description:Whether we should ghost point neighbors of secondary face elements, and consequently also their mortar interface couples.
interpolate_normalsTrueWhether to interpolate the nodal normals (e.g. classic idea of evaluating field at quadrature points). If this is set to false, then non-interpolated nodal normals will be used, and then the _normals member should be indexed with _i instead of _qp
Default:True
C++ Type:bool
Controllable:No
Description:Whether to interpolate the nodal normals (e.g. classic idea of evaluating field at quadrature points). If this is set to false, then non-interpolated nodal normals will be used, and then the _normals member should be indexed with _i instead of _qp
minimum_projection_angle40Parameter to control which angle (in degrees) is admissible for the creation of mortar segments. If set to a value close to zero, very oblique projections are allowed, which can result in mortar segments solving physics not meaningfully, and overprojection of primary nodes onto the mortar segment mesh in extreme cases. This parameter is mostly intended for mortar mesh debugging purposes in two dimensions.
Default:40
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Parameter to control which angle (in degrees) is admissible for the creation of mortar segments. If set to a value close to zero, very oblique projections are allowed, which can result in mortar segments solving physics not meaningfully, and overprojection of primary nodes onto the mortar segment mesh in extreme cases. This parameter is mostly intended for mortar mesh debugging purposes in two dimensions.
partial_pressure_namepartial_pressureThe name of the partial pressure if multiple partial pressure are defined on the same block
Default:partial_pressure
C++ Type:std::string
Controllable:No
Description:The name of the partial pressure if multiple partial pressure are defined on the same block
periodicFalseWhether this constraint is going to be used to enforce a periodic condition. This has the effect of changing the normals vector for projection from outward to inward facing
Default:False
C++ Type:bool
Controllable:No
Description:Whether this constraint is going to be used to enforce a periodic condition. This has the effect of changing the normals vector for projection from outward to inward facing
primary_variablePrimal variable on primary surface. If this parameter is not provided then the primary variable will be initialized to the secondary variable
C++ Type:VariableName
Unit:(no unit assumed)
Controllable:No
Description:Primal variable on primary surface. If this parameter is not provided then the primary variable will be initialized to the secondary variable
quadratureDEFAULTQuadrature rule to use on mortar segments. For 2D mortar DEFAULT is recommended. For 3D mortar, QUAD meshes are integrated using triangle mortar segments. While DEFAULT quadrature order is typically sufficiently accurate, exact integration of QUAD mortar faces requires SECOND order quadrature for FIRST variables and FOURTH order quadrature for SECOND order variables.
Default:DEFAULT
C++ Type:MooseEnum
Options:DEFAULT, FIRST, SECOND, THIRD, FOURTH
Controllable:No
Description:Quadrature rule to use on mortar segments. For 2D mortar DEFAULT is recommended. For 3D mortar, QUAD meshes are integrated using triangle mortar segments. While DEFAULT quadrature order is typically sufficiently accurate, exact integration of QUAD mortar faces requires SECOND order quadrature for FIRST variables and FOURTH order quadrature for SECOND order variables.
secondary_variablePrimal variable on secondary surface.
C++ Type:VariableName
Unit:(no unit assumed)
Controllable:No
Description:Primal variable on secondary surface.
use_petrov_galerkinFalseWhether to use the Petrov-Galerkin approach for the mortar-based constraints. If set to true, we use the standard basis as the test function and dual basis as the shape function for the interpolation of the Lagrange multiplier variable.
Default:False
C++ Type:bool
Controllable:No
Description:Whether to use the Petrov-Galerkin approach for the mortar-based constraints. If set to true, we use the standard basis as the test function and dual basis as the shape function for the interpolation of the Lagrange multiplier variable.
variableThe name of the lagrange multiplier variable that this constraint is applied to. This parameter may not be supplied in the case of using penalty methods for example
C++ Type:NonlinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The name of the lagrange multiplier variable that this constraint is applied to. This parameter may not be supplied in the case of using penalty methods for example

Optional Parameters

absolute_value_vector_tagsThe tags for the vectors this residual object should fill with the absolute value of the residual contribution
C++ 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 fill
C++ 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 fill
C++ 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 fill
Default: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 fill
Default: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.
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 form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
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/sorption_constraint/test.i)
(examples/TRISO/accident_simulation/triso2D_accident_ad.i)

References

No citations exist within this document.

(test/tests/sorption_constraint/test.i)

# This test is to verify the implementation of MassSorptionConstraint constraint.
# It contains two 2D blocks separated by an unmeshed flat gap.
# 1D mortar subdomains are set up to enforce sorption and preserve flux between the blocks.
# Checks are performed to verify concentration conservation, sorption behavior, and flux preservation.

[GlobalParams]
  order = FIRST
  family = LAGRANGE
[]

[Mesh]
  coord_type = XYZ
  patch_update_strategy = auto

  [inner] # create inner block
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = 0.5
    ymin = 0
    ymax = 1
    nx = 10
    ny = 10
    boundary_name_prefix = inner
    elem_type = QUAD4
  []
  [inner_id]
    type = SubdomainIDGenerator
    input = inner
    subdomain_id = 1
  []

  [outer] # create outer block
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0.51
    xmax = 1.01
    ymin = 0
    ymax = 1
    nx = 10
    ny = 10
    boundary_name_prefix = outer
    boundary_id_offset = 10
    elem_type = QUAD4
  []
  [outer_id]
    type = SubdomainIDGenerator
    input = outer
    subdomain_id = 2
  []

  [combined]
    type = MeshCollectionGenerator
    inputs = 'inner_id outer_id'
  []
  [block_rename]
    type = RenameBlockGenerator
    input = combined
    old_block = '1 2'
    new_block = 'inner outer'
  []

  [primary] # create left mortar subdomain between the inner and outer blocks
    type = LowerDBlockFromSidesetGenerator
    input = block_rename
    sidesets = 1
    new_block_id = 100
    new_block_name = primary
  []
  [secondary] # create right mortar subdomain between the inner and outer blocks
    type = LowerDBlockFromSidesetGenerator
    input = primary
    sidesets = 13
    new_block_id = 200
    new_block_name = secondary
  []
[]

[Variables]
  [u] # concentration
    block = 'inner outer'
    initial_condition = 2
  []
  [temperature]
    block = 'inner outer'
    initial_condition = 550
  []
  [lm_u] # constraint on the concentration in the gap
    block = secondary
  []
  [lm_dx] # constraint on the flux in the x direction in the gap
    block = secondary
  []
  [lm_dy] # constraint on the flux in the y direction in the gap
    block = secondary
  []
[]

[BCs]
  [u_left]
    type = ADDirichletBC
    value = 1
    variable = u
    boundary = inner_left
  []
  [u_right]
    type = ADDirichletBC
    value = 1
    variable = u
    boundary = outer_right
  []
  [temperature_left]
    type = ADDirichletBC
    value = 1100
    variable = temperature
    boundary = inner_left
  []
  [temperature_primary]
    type = ADDirichletBC
    value = 1000
    variable = temperature
    boundary = inner_right
  []
  [temperature_secondary]
    type = ADDirichletBC
    value = 900
    variable = temperature
    boundary = outer_left
  []
  [temperature_right]
    type = ADDirichletBC
    value = 0
    variable = temperature
    boundary = outer_right
  []
[]

[Kernels]
  [u]
    type = ADMatDiffusion
    variable = u
    diffusivity = diffusivity
    block = 'inner outer'
  []
  [temperature]
    type = ADHeatConduction
    variable = temperature
    thermal_conductivity = thermal_conductivity
    block = 'inner outer'
  []
[]

[Constraints]
  [u]
    type = MassSorptionConstraint
    variable = lm_u
    primary_variable = u
    primary_boundary = inner_right
    primary_subdomain = primary
    secondary_variable = u
    secondary_boundary = outer_left
    secondary_subdomain = secondary
    partial_pressure_name = partial_pressure
    epsilon = 1e-9
    correct_edge_dropping = true
  []
  [dx]
    type = MassFluxConstraint
    variable = lm_dx
    primary_variable = u
    diffusivity_primary = diffusivity
    primary_boundary = inner_right
    primary_subdomain = primary
    secondary_variable = u
    diffusivity_secondary = diffusivity
    secondary_boundary = outer_left
    secondary_subdomain = secondary
    component = 0
    epsilon = 1e-5
    correct_edge_dropping = true
  []
  [dy]
    type = MassFluxConstraint
    variable = lm_dy
    primary_variable = u
    diffusivity_primary = diffusivity
    primary_boundary = inner_right
    primary_subdomain = primary
    secondary_variable = u
    diffusivity_secondary = diffusivity
    secondary_boundary = outer_left
    secondary_subdomain = secondary
    component = 1
    epsilon = 1e-5
    correct_edge_dropping = true
  []
[]

[Materials]
  [properties_inner]
    type = ADGenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity density'
    prop_values = '1 1 1'
    block = inner
    outputs = all
  []
  [partial_pressure_inner]
    type = ADSorptionPartialPressure
    A = 19.3
    B = -47300
    D = 1.51
    E = 4340
    d1 = 3.4
    d2 = 6.15e-4
    unit_scale = 1
    density = density
    concentration = u
    temperature = temperature
    block = inner
    outputs = 'all'
    output_properties = partial_pressure
  []
  [properties_outer]
    type = ADGenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity density'
    prop_values = '2 2 1'
    block = outer
    outputs = all
  []
  [partial_pressure_outer]
    type = ADSorptionPartialPressure
    A = 24
    B = -35700
    D = -1.56
    E = 6120
    d1 = 2.04
    d2 = 1.79e-3
    unit_scale = 1
    density = density
    concentration = u
    temperature = temperature
    block = outer
    outputs = 'all'
    output_properties = partial_pressure
  []
[]

[Functions]
  [u_mid_diff]
    type = ParsedFunction
    symbol_names = 'u_mid_inner u_mid_outer'
    symbol_values = 'u_mid_inner u_mid_outer'
    expression = '(abs(u_mid_outer) - abs(u_mid_inner)) / abs(u_mid_inner)'
  []
  [pressure_mid_error]
    type = ParsedFunction
    symbol_names = 'pressure_mid_inner pressure_mid_outer'
    symbol_values = 'pressure_mid_inner pressure_mid_outer'
    expression = '(abs(pressure_mid_outer) - abs(pressure_mid_inner)) / abs(pressure_mid_inner)'
  []
  [flux_error]
    type = ParsedFunction
    symbol_names = 'flux_inner flux_outer'
    symbol_values = 'flux_inner flux_outer'
    expression = '(abs(flux_outer) - abs(flux_inner)) / abs(flux_inner)'
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Steady
  solve_type = NEWTON
  automatic_scaling = true
  scaling_group_variables = 'u; temperature; lm_u; lm_dx; lm_dy'

  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu       NONZERO               1e-15' # provide nonzero on-diagonals for lagrange multiplier
  line_search = none

  nl_rel_tol = 1e-15
  nl_abs_tol = 1e-9
  l_tol = 1e-3
  l_max_its = 50
  nl_max_its = 50
[]

[Postprocessors]
  [u_mid_inner] # verify sorption behavior
    type = PointValue
    variable = u
    point = '0.5 0.5 0'
    outputs = csv
  []
  [u_mid_outer]
    type = PointValue
    variable = u
    point = '0.51 0.5 0'
    outputs = csv
  []
  [u_mid_diff]
    type = FunctionValuePostprocessor
    function = u_mid_diff
    outputs = console
  []
  [temperature_mid_inner]
    type = PointValue
    variable = temperature
    point = '0.5 0.5 0'
    outputs = csv
  []
  [temperature_mid_outer]
    type = PointValue
    variable = temperature
    point = '0.51 0.5 0'
    outputs = csv
  []
  [pressure_mid_inner]
    type = PointValue
    variable = partial_pressure
    point = '0.51 0.5 0'
    outputs = csv
  []
  [pressure_mid_outer]
    type = PointValue
    variable = partial_pressure
    point = '0.51 0.5 0'
    outputs = csv
  []
  [pressure_mid_error]
    type = FunctionValuePostprocessor
    function = pressure_mid_error
    outputs = console
  []
  [flux_inner] # verify flux preservation
    type = ADSideDiffusiveFluxIntegral
    variable = u
    boundary = inner_right
    diffusivity = diffusivity
    outputs = csv
  []
  [flux_outer]
    type = ADSideDiffusiveFluxIntegral
    variable = u
    boundary = outer_left
    diffusivity = diffusivity
    outputs = csv
  []
  [flux_error]
    type = FunctionValuePostprocessor
    function = flux_error
    outputs = console
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

[Debug]
  show_var_residual_norms = true
[]

(test/tests/sorption_constraint/test.i)

# This test is to verify the implementation of MassSorptionConstraint constraint.
# It contains two 2D blocks separated by an unmeshed flat gap.
# 1D mortar subdomains are set up to enforce sorption and preserve flux between the blocks.
# Checks are performed to verify concentration conservation, sorption behavior, and flux preservation.

[GlobalParams]
  order = FIRST
  family = LAGRANGE
[]

[Mesh]
  coord_type = XYZ
  patch_update_strategy = auto

  [inner] # create inner block
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = 0.5
    ymin = 0
    ymax = 1
    nx = 10
    ny = 10
    boundary_name_prefix = inner
    elem_type = QUAD4
  []
  [inner_id]
    type = SubdomainIDGenerator
    input = inner
    subdomain_id = 1
  []

  [outer] # create outer block
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0.51
    xmax = 1.01
    ymin = 0
    ymax = 1
    nx = 10
    ny = 10
    boundary_name_prefix = outer
    boundary_id_offset = 10
    elem_type = QUAD4
  []
  [outer_id]
    type = SubdomainIDGenerator
    input = outer
    subdomain_id = 2
  []

  [combined]
    type = MeshCollectionGenerator
    inputs = 'inner_id outer_id'
  []
  [block_rename]
    type = RenameBlockGenerator
    input = combined
    old_block = '1 2'
    new_block = 'inner outer'
  []

  [primary] # create left mortar subdomain between the inner and outer blocks
    type = LowerDBlockFromSidesetGenerator
    input = block_rename
    sidesets = 1
    new_block_id = 100
    new_block_name = primary
  []
  [secondary] # create right mortar subdomain between the inner and outer blocks
    type = LowerDBlockFromSidesetGenerator
    input = primary
    sidesets = 13
    new_block_id = 200
    new_block_name = secondary
  []
[]

[Variables]
  [u] # concentration
    block = 'inner outer'
    initial_condition = 2
  []
  [temperature]
    block = 'inner outer'
    initial_condition = 550
  []
  [lm_u] # constraint on the concentration in the gap
    block = secondary
  []
  [lm_dx] # constraint on the flux in the x direction in the gap
    block = secondary
  []
  [lm_dy] # constraint on the flux in the y direction in the gap
    block = secondary
  []
[]

[BCs]
  [u_left]
    type = ADDirichletBC
    value = 1
    variable = u
    boundary = inner_left
  []
  [u_right]
    type = ADDirichletBC
    value = 1
    variable = u
    boundary = outer_right
  []
  [temperature_left]
    type = ADDirichletBC
    value = 1100
    variable = temperature
    boundary = inner_left
  []
  [temperature_primary]
    type = ADDirichletBC
    value = 1000
    variable = temperature
    boundary = inner_right
  []
  [temperature_secondary]
    type = ADDirichletBC
    value = 900
    variable = temperature
    boundary = outer_left
  []
  [temperature_right]
    type = ADDirichletBC
    value = 0
    variable = temperature
    boundary = outer_right
  []
[]

[Kernels]
  [u]
    type = ADMatDiffusion
    variable = u
    diffusivity = diffusivity
    block = 'inner outer'
  []
  [temperature]
    type = ADHeatConduction
    variable = temperature
    thermal_conductivity = thermal_conductivity
    block = 'inner outer'
  []
[]

[Constraints]
  [u]
    type = MassSorptionConstraint
    variable = lm_u
    primary_variable = u
    primary_boundary = inner_right
    primary_subdomain = primary
    secondary_variable = u
    secondary_boundary = outer_left
    secondary_subdomain = secondary
    partial_pressure_name = partial_pressure
    epsilon = 1e-9
    correct_edge_dropping = true
  []
  [dx]
    type = MassFluxConstraint
    variable = lm_dx
    primary_variable = u
    diffusivity_primary = diffusivity
    primary_boundary = inner_right
    primary_subdomain = primary
    secondary_variable = u
    diffusivity_secondary = diffusivity
    secondary_boundary = outer_left
    secondary_subdomain = secondary
    component = 0
    epsilon = 1e-5
    correct_edge_dropping = true
  []
  [dy]
    type = MassFluxConstraint
    variable = lm_dy
    primary_variable = u
    diffusivity_primary = diffusivity
    primary_boundary = inner_right
    primary_subdomain = primary
    secondary_variable = u
    diffusivity_secondary = diffusivity
    secondary_boundary = outer_left
    secondary_subdomain = secondary
    component = 1
    epsilon = 1e-5
    correct_edge_dropping = true
  []
[]

[Materials]
  [properties_inner]
    type = ADGenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity density'
    prop_values = '1 1 1'
    block = inner
    outputs = all
  []
  [partial_pressure_inner]
    type = ADSorptionPartialPressure
    A = 19.3
    B = -47300
    D = 1.51
    E = 4340
    d1 = 3.4
    d2 = 6.15e-4
    unit_scale = 1
    density = density
    concentration = u
    temperature = temperature
    block = inner
    outputs = 'all'
    output_properties = partial_pressure
  []
  [properties_outer]
    type = ADGenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity density'
    prop_values = '2 2 1'
    block = outer
    outputs = all
  []
  [partial_pressure_outer]
    type = ADSorptionPartialPressure
    A = 24
    B = -35700
    D = -1.56
    E = 6120
    d1 = 2.04
    d2 = 1.79e-3
    unit_scale = 1
    density = density
    concentration = u
    temperature = temperature
    block = outer
    outputs = 'all'
    output_properties = partial_pressure
  []
[]

[Functions]
  [u_mid_diff]
    type = ParsedFunction
    symbol_names = 'u_mid_inner u_mid_outer'
    symbol_values = 'u_mid_inner u_mid_outer'
    expression = '(abs(u_mid_outer) - abs(u_mid_inner)) / abs(u_mid_inner)'
  []
  [pressure_mid_error]
    type = ParsedFunction
    symbol_names = 'pressure_mid_inner pressure_mid_outer'
    symbol_values = 'pressure_mid_inner pressure_mid_outer'
    expression = '(abs(pressure_mid_outer) - abs(pressure_mid_inner)) / abs(pressure_mid_inner)'
  []
  [flux_error]
    type = ParsedFunction
    symbol_names = 'flux_inner flux_outer'
    symbol_values = 'flux_inner flux_outer'
    expression = '(abs(flux_outer) - abs(flux_inner)) / abs(flux_inner)'
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Steady
  solve_type = NEWTON
  automatic_scaling = true
  scaling_group_variables = 'u; temperature; lm_u; lm_dx; lm_dy'

  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu       NONZERO               1e-15' # provide nonzero on-diagonals for lagrange multiplier
  line_search = none

  nl_rel_tol = 1e-15
  nl_abs_tol = 1e-9
  l_tol = 1e-3
  l_max_its = 50
  nl_max_its = 50
[]

[Postprocessors]
  [u_mid_inner] # verify sorption behavior
    type = PointValue
    variable = u
    point = '0.5 0.5 0'
    outputs = csv
  []
  [u_mid_outer]
    type = PointValue
    variable = u
    point = '0.51 0.5 0'
    outputs = csv
  []
  [u_mid_diff]
    type = FunctionValuePostprocessor
    function = u_mid_diff
    outputs = console
  []
  [temperature_mid_inner]
    type = PointValue
    variable = temperature
    point = '0.5 0.5 0'
    outputs = csv
  []
  [temperature_mid_outer]
    type = PointValue
    variable = temperature
    point = '0.51 0.5 0'
    outputs = csv
  []
  [pressure_mid_inner]
    type = PointValue
    variable = partial_pressure
    point = '0.51 0.5 0'
    outputs = csv
  []
  [pressure_mid_outer]
    type = PointValue
    variable = partial_pressure
    point = '0.51 0.5 0'
    outputs = csv
  []
  [pressure_mid_error]
    type = FunctionValuePostprocessor
    function = pressure_mid_error
    outputs = console
  []
  [flux_inner] # verify flux preservation
    type = ADSideDiffusiveFluxIntegral
    variable = u
    boundary = inner_right
    diffusivity = diffusivity
    outputs = csv
  []
  [flux_outer]
    type = ADSideDiffusiveFluxIntegral
    variable = u
    boundary = outer_left
    diffusivity = diffusivity
    outputs = csv
  []
  [flux_error]
    type = FunctionValuePostprocessor
    function = flux_error
    outputs = console
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

[Debug]
  show_var_residual_norms = true
[]

(test/tests/sorption_constraint/test.i)

# This test is to verify the implementation of MassSorptionConstraint constraint.
# It contains two 2D blocks separated by an unmeshed flat gap.
# 1D mortar subdomains are set up to enforce sorption and preserve flux between the blocks.
# Checks are performed to verify concentration conservation, sorption behavior, and flux preservation.

[GlobalParams]
  order = FIRST
  family = LAGRANGE
[]

[Mesh]
  coord_type = XYZ
  patch_update_strategy = auto

  [inner] # create inner block
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = 0.5
    ymin = 0
    ymax = 1
    nx = 10
    ny = 10
    boundary_name_prefix = inner
    elem_type = QUAD4
  []
  [inner_id]
    type = SubdomainIDGenerator
    input = inner
    subdomain_id = 1
  []

  [outer] # create outer block
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0.51
    xmax = 1.01
    ymin = 0
    ymax = 1
    nx = 10
    ny = 10
    boundary_name_prefix = outer
    boundary_id_offset = 10
    elem_type = QUAD4
  []
  [outer_id]
    type = SubdomainIDGenerator
    input = outer
    subdomain_id = 2
  []

  [combined]
    type = MeshCollectionGenerator
    inputs = 'inner_id outer_id'
  []
  [block_rename]
    type = RenameBlockGenerator
    input = combined
    old_block = '1 2'
    new_block = 'inner outer'
  []

  [primary] # create left mortar subdomain between the inner and outer blocks
    type = LowerDBlockFromSidesetGenerator
    input = block_rename
    sidesets = 1
    new_block_id = 100
    new_block_name = primary
  []
  [secondary] # create right mortar subdomain between the inner and outer blocks
    type = LowerDBlockFromSidesetGenerator
    input = primary
    sidesets = 13
    new_block_id = 200
    new_block_name = secondary
  []
[]

[Variables]
  [u] # concentration
    block = 'inner outer'
    initial_condition = 2
  []
  [temperature]
    block = 'inner outer'
    initial_condition = 550
  []
  [lm_u] # constraint on the concentration in the gap
    block = secondary
  []
  [lm_dx] # constraint on the flux in the x direction in the gap
    block = secondary
  []
  [lm_dy] # constraint on the flux in the y direction in the gap
    block = secondary
  []
[]

[BCs]
  [u_left]
    type = ADDirichletBC
    value = 1
    variable = u
    boundary = inner_left
  []
  [u_right]
    type = ADDirichletBC
    value = 1
    variable = u
    boundary = outer_right
  []
  [temperature_left]
    type = ADDirichletBC
    value = 1100
    variable = temperature
    boundary = inner_left
  []
  [temperature_primary]
    type = ADDirichletBC
    value = 1000
    variable = temperature
    boundary = inner_right
  []
  [temperature_secondary]
    type = ADDirichletBC
    value = 900
    variable = temperature
    boundary = outer_left
  []
  [temperature_right]
    type = ADDirichletBC
    value = 0
    variable = temperature
    boundary = outer_right
  []
[]

[Kernels]
  [u]
    type = ADMatDiffusion
    variable = u
    diffusivity = diffusivity
    block = 'inner outer'
  []
  [temperature]
    type = ADHeatConduction
    variable = temperature
    thermal_conductivity = thermal_conductivity
    block = 'inner outer'
  []
[]

[Constraints]
  [u]
    type = MassSorptionConstraint
    variable = lm_u
    primary_variable = u
    primary_boundary = inner_right
    primary_subdomain = primary
    secondary_variable = u
    secondary_boundary = outer_left
    secondary_subdomain = secondary
    partial_pressure_name = partial_pressure
    epsilon = 1e-9
    correct_edge_dropping = true
  []
  [dx]
    type = MassFluxConstraint
    variable = lm_dx
    primary_variable = u
    diffusivity_primary = diffusivity
    primary_boundary = inner_right
    primary_subdomain = primary
    secondary_variable = u
    diffusivity_secondary = diffusivity
    secondary_boundary = outer_left
    secondary_subdomain = secondary
    component = 0
    epsilon = 1e-5
    correct_edge_dropping = true
  []
  [dy]
    type = MassFluxConstraint
    variable = lm_dy
    primary_variable = u
    diffusivity_primary = diffusivity
    primary_boundary = inner_right
    primary_subdomain = primary
    secondary_variable = u
    diffusivity_secondary = diffusivity
    secondary_boundary = outer_left
    secondary_subdomain = secondary
    component = 1
    epsilon = 1e-5
    correct_edge_dropping = true
  []
[]

[Materials]
  [properties_inner]
    type = ADGenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity density'
    prop_values = '1 1 1'
    block = inner
    outputs = all
  []
  [partial_pressure_inner]
    type = ADSorptionPartialPressure
    A = 19.3
    B = -47300
    D = 1.51
    E = 4340
    d1 = 3.4
    d2 = 6.15e-4
    unit_scale = 1
    density = density
    concentration = u
    temperature = temperature
    block = inner
    outputs = 'all'
    output_properties = partial_pressure
  []
  [properties_outer]
    type = ADGenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity density'
    prop_values = '2 2 1'
    block = outer
    outputs = all
  []
  [partial_pressure_outer]
    type = ADSorptionPartialPressure
    A = 24
    B = -35700
    D = -1.56
    E = 6120
    d1 = 2.04
    d2 = 1.79e-3
    unit_scale = 1
    density = density
    concentration = u
    temperature = temperature
    block = outer
    outputs = 'all'
    output_properties = partial_pressure
  []
[]

[Functions]
  [u_mid_diff]
    type = ParsedFunction
    symbol_names = 'u_mid_inner u_mid_outer'
    symbol_values = 'u_mid_inner u_mid_outer'
    expression = '(abs(u_mid_outer) - abs(u_mid_inner)) / abs(u_mid_inner)'
  []
  [pressure_mid_error]
    type = ParsedFunction
    symbol_names = 'pressure_mid_inner pressure_mid_outer'
    symbol_values = 'pressure_mid_inner pressure_mid_outer'
    expression = '(abs(pressure_mid_outer) - abs(pressure_mid_inner)) / abs(pressure_mid_inner)'
  []
  [flux_error]
    type = ParsedFunction
    symbol_names = 'flux_inner flux_outer'
    symbol_values = 'flux_inner flux_outer'
    expression = '(abs(flux_outer) - abs(flux_inner)) / abs(flux_inner)'
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Steady
  solve_type = NEWTON
  automatic_scaling = true
  scaling_group_variables = 'u; temperature; lm_u; lm_dx; lm_dy'

  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu       NONZERO               1e-15' # provide nonzero on-diagonals for lagrange multiplier
  line_search = none

  nl_rel_tol = 1e-15
  nl_abs_tol = 1e-9
  l_tol = 1e-3
  l_max_its = 50
  nl_max_its = 50
[]

[Postprocessors]
  [u_mid_inner] # verify sorption behavior
    type = PointValue
    variable = u
    point = '0.5 0.5 0'
    outputs = csv
  []
  [u_mid_outer]
    type = PointValue
    variable = u
    point = '0.51 0.5 0'
    outputs = csv
  []
  [u_mid_diff]
    type = FunctionValuePostprocessor
    function = u_mid_diff
    outputs = console
  []
  [temperature_mid_inner]
    type = PointValue
    variable = temperature
    point = '0.5 0.5 0'
    outputs = csv
  []
  [temperature_mid_outer]
    type = PointValue
    variable = temperature
    point = '0.51 0.5 0'
    outputs = csv
  []
  [pressure_mid_inner]
    type = PointValue
    variable = partial_pressure
    point = '0.51 0.5 0'
    outputs = csv
  []
  [pressure_mid_outer]
    type = PointValue
    variable = partial_pressure
    point = '0.51 0.5 0'
    outputs = csv
  []
  [pressure_mid_error]
    type = FunctionValuePostprocessor
    function = pressure_mid_error
    outputs = console
  []
  [flux_inner] # verify flux preservation
    type = ADSideDiffusiveFluxIntegral
    variable = u
    boundary = inner_right
    diffusivity = diffusivity
    outputs = csv
  []
  [flux_outer]
    type = ADSideDiffusiveFluxIntegral
    variable = u
    boundary = outer_left
    diffusivity = diffusivity
    outputs = csv
  []
  [flux_error]
    type = FunctionValuePostprocessor
    function = flux_error
    outputs = console
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

[Debug]
  show_var_residual_norms = true
[]

(test/tests/sorption_constraint/test.i)

# This test is to verify the implementation of MassSorptionConstraint constraint.
# It contains two 2D blocks separated by an unmeshed flat gap.
# 1D mortar subdomains are set up to enforce sorption and preserve flux between the blocks.
# Checks are performed to verify concentration conservation, sorption behavior, and flux preservation.

[GlobalParams]
  order = FIRST
  family = LAGRANGE
[]

[Mesh]
  coord_type = XYZ
  patch_update_strategy = auto

  [inner] # create inner block
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = 0.5
    ymin = 0
    ymax = 1
    nx = 10
    ny = 10
    boundary_name_prefix = inner
    elem_type = QUAD4
  []
  [inner_id]
    type = SubdomainIDGenerator
    input = inner
    subdomain_id = 1
  []

  [outer] # create outer block
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0.51
    xmax = 1.01
    ymin = 0
    ymax = 1
    nx = 10
    ny = 10
    boundary_name_prefix = outer
    boundary_id_offset = 10
    elem_type = QUAD4
  []
  [outer_id]
    type = SubdomainIDGenerator
    input = outer
    subdomain_id = 2
  []

  [combined]
    type = MeshCollectionGenerator
    inputs = 'inner_id outer_id'
  []
  [block_rename]
    type = RenameBlockGenerator
    input = combined
    old_block = '1 2'
    new_block = 'inner outer'
  []

  [primary] # create left mortar subdomain between the inner and outer blocks
    type = LowerDBlockFromSidesetGenerator
    input = block_rename
    sidesets = 1
    new_block_id = 100
    new_block_name = primary
  []
  [secondary] # create right mortar subdomain between the inner and outer blocks
    type = LowerDBlockFromSidesetGenerator
    input = primary
    sidesets = 13
    new_block_id = 200
    new_block_name = secondary
  []
[]

[Variables]
  [u] # concentration
    block = 'inner outer'
    initial_condition = 2
  []
  [temperature]
    block = 'inner outer'
    initial_condition = 550
  []
  [lm_u] # constraint on the concentration in the gap
    block = secondary
  []
  [lm_dx] # constraint on the flux in the x direction in the gap
    block = secondary
  []
  [lm_dy] # constraint on the flux in the y direction in the gap
    block = secondary
  []
[]

[BCs]
  [u_left]
    type = ADDirichletBC
    value = 1
    variable = u
    boundary = inner_left
  []
  [u_right]
    type = ADDirichletBC
    value = 1
    variable = u
    boundary = outer_right
  []
  [temperature_left]
    type = ADDirichletBC
    value = 1100
    variable = temperature
    boundary = inner_left
  []
  [temperature_primary]
    type = ADDirichletBC
    value = 1000
    variable = temperature
    boundary = inner_right
  []
  [temperature_secondary]
    type = ADDirichletBC
    value = 900
    variable = temperature
    boundary = outer_left
  []
  [temperature_right]
    type = ADDirichletBC
    value = 0
    variable = temperature
    boundary = outer_right
  []
[]

[Kernels]
  [u]
    type = ADMatDiffusion
    variable = u
    diffusivity = diffusivity
    block = 'inner outer'
  []
  [temperature]
    type = ADHeatConduction
    variable = temperature
    thermal_conductivity = thermal_conductivity
    block = 'inner outer'
  []
[]

[Constraints]
  [u]
    type = MassSorptionConstraint
    variable = lm_u
    primary_variable = u
    primary_boundary = inner_right
    primary_subdomain = primary
    secondary_variable = u
    secondary_boundary = outer_left
    secondary_subdomain = secondary
    partial_pressure_name = partial_pressure
    epsilon = 1e-9
    correct_edge_dropping = true
  []
  [dx]
    type = MassFluxConstraint
    variable = lm_dx
    primary_variable = u
    diffusivity_primary = diffusivity
    primary_boundary = inner_right
    primary_subdomain = primary
    secondary_variable = u
    diffusivity_secondary = diffusivity
    secondary_boundary = outer_left
    secondary_subdomain = secondary
    component = 0
    epsilon = 1e-5
    correct_edge_dropping = true
  []
  [dy]
    type = MassFluxConstraint
    variable = lm_dy
    primary_variable = u
    diffusivity_primary = diffusivity
    primary_boundary = inner_right
    primary_subdomain = primary
    secondary_variable = u
    diffusivity_secondary = diffusivity
    secondary_boundary = outer_left
    secondary_subdomain = secondary
    component = 1
    epsilon = 1e-5
    correct_edge_dropping = true
  []
[]

[Materials]
  [properties_inner]
    type = ADGenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity density'
    prop_values = '1 1 1'
    block = inner
    outputs = all
  []
  [partial_pressure_inner]
    type = ADSorptionPartialPressure
    A = 19.3
    B = -47300
    D = 1.51
    E = 4340
    d1 = 3.4
    d2 = 6.15e-4
    unit_scale = 1
    density = density
    concentration = u
    temperature = temperature
    block = inner
    outputs = 'all'
    output_properties = partial_pressure
  []
  [properties_outer]
    type = ADGenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity density'
    prop_values = '2 2 1'
    block = outer
    outputs = all
  []
  [partial_pressure_outer]
    type = ADSorptionPartialPressure
    A = 24
    B = -35700
    D = -1.56
    E = 6120
    d1 = 2.04
    d2 = 1.79e-3
    unit_scale = 1
    density = density
    concentration = u
    temperature = temperature
    block = outer
    outputs = 'all'
    output_properties = partial_pressure
  []
[]

[Functions]
  [u_mid_diff]
    type = ParsedFunction
    symbol_names = 'u_mid_inner u_mid_outer'
    symbol_values = 'u_mid_inner u_mid_outer'
    expression = '(abs(u_mid_outer) - abs(u_mid_inner)) / abs(u_mid_inner)'
  []
  [pressure_mid_error]
    type = ParsedFunction
    symbol_names = 'pressure_mid_inner pressure_mid_outer'
    symbol_values = 'pressure_mid_inner pressure_mid_outer'
    expression = '(abs(pressure_mid_outer) - abs(pressure_mid_inner)) / abs(pressure_mid_inner)'
  []
  [flux_error]
    type = ParsedFunction
    symbol_names = 'flux_inner flux_outer'
    symbol_values = 'flux_inner flux_outer'
    expression = '(abs(flux_outer) - abs(flux_inner)) / abs(flux_inner)'
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Steady
  solve_type = NEWTON
  automatic_scaling = true
  scaling_group_variables = 'u; temperature; lm_u; lm_dx; lm_dy'

  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu       NONZERO               1e-15' # provide nonzero on-diagonals for lagrange multiplier
  line_search = none

  nl_rel_tol = 1e-15
  nl_abs_tol = 1e-9
  l_tol = 1e-3
  l_max_its = 50
  nl_max_its = 50
[]

[Postprocessors]
  [u_mid_inner] # verify sorption behavior
    type = PointValue
    variable = u
    point = '0.5 0.5 0'
    outputs = csv
  []
  [u_mid_outer]
    type = PointValue
    variable = u
    point = '0.51 0.5 0'
    outputs = csv
  []
  [u_mid_diff]
    type = FunctionValuePostprocessor
    function = u_mid_diff
    outputs = console
  []
  [temperature_mid_inner]
    type = PointValue
    variable = temperature
    point = '0.5 0.5 0'
    outputs = csv
  []
  [temperature_mid_outer]
    type = PointValue
    variable = temperature
    point = '0.51 0.5 0'
    outputs = csv
  []
  [pressure_mid_inner]
    type = PointValue
    variable = partial_pressure
    point = '0.51 0.5 0'
    outputs = csv
  []
  [pressure_mid_outer]
    type = PointValue
    variable = partial_pressure
    point = '0.51 0.5 0'
    outputs = csv
  []
  [pressure_mid_error]
    type = FunctionValuePostprocessor
    function = pressure_mid_error
    outputs = console
  []
  [flux_inner] # verify flux preservation
    type = ADSideDiffusiveFluxIntegral
    variable = u
    boundary = inner_right
    diffusivity = diffusivity
    outputs = csv
  []
  [flux_outer]
    type = ADSideDiffusiveFluxIntegral
    variable = u
    boundary = outer_left
    diffusivity = diffusivity
    outputs = csv
  []
  [flux_error]
    type = FunctionValuePostprocessor
    function = flux_error
    outputs = console
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

[Debug]
  show_var_residual_norms = true
[]

(examples/TRISO/accident_simulation/triso2D_accident_ad.i)

# This example is 2D-RZ analysis of a TRISO fuel particle. Fully coupled
# heat transfer and solid mechanics, plus diffusion of the fission product
# species cesium (Cs) are simulated.  The mesh includes contact surfaces
# between the buffer and IPyC layers to facilitate a gap opening between
# these layers. These surfaces are initially in mechanical contact but
# are assumed to have no strength in tension. A coarse mesh is used to
# provide a short run time.

# The calculation simulates fuel-life in three steps. The first step is an
# irradiation period, where constant power and a fixed particle surface
# temperature (1500 K) are assumed over a lifetime of 76 Ms (2.4 yrs).
# For the second step, fuel removal and storage are simulated by setting
# the reactor power and Cs source terms to zero, reducing the particle
# surface temperature to ambient (300 K), and then holding it
# for 100 days. A third and final step simulates accident
# behavior by increasing the particle surface temperature from ambient
# to 2073 K over 2 hrs, and then holding it at this elevated temperature
# for an additional 200 hrs. At the particle outer boundary, the Cs
# concentration is held at zero and the pressure at ambient during the
# entire simulation. The particle is assumed to be stress-free at an
# initial temperature of 1500 K.
#
# Details about this simulation are given in Section 4 of the following
# article: J. D. Hales, R. L. Williamson, S. R. Novascone, D. M. Perez,
# B. W. Spencer and G. Pastore, "Multidimensional multiphysics simulation
# of TRISO particle fuel", Journal of Nuclear Materials, Vol. 443, p. 531,
# 2013.

# This is a version using a thermomechanical mortar approach. It uses
# Automatic Differentiation classes and models gap mass transfer using
# flux preserving and sorption mortar constraints. Sorption constants are
# given in Table 1 of the following article: A. Londono-Hurtado, I.
# Szlufarska, R. Bratton and D. Morgan, "A review of fission product
# sorption in carbon structures", Journal of Nuclear Materials, Vol. 426,
# p. 254, 2012.


initial_fuel_density = 11000.0

[GlobalParams]
  order = SECOND
  family = LAGRANGE
  displacements = 'disp_x disp_y'
  flux_conversion_factor = 0.85
  use_automatic_differentiation = true
[]

[Mesh]
  coord_type = RZ
  [file]
    type = FileMeshGenerator
    file = triso2Dmed.e
  []
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
  converge_on = 'disp_x disp_y temp conc'
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [temp]
    initial_condition = 1500.0
  []
  [conc]
    initial_condition = 0.0
  []
  [conc_lm]
    block = pellet_clad_mechanical_secondary_subdomain
  []
  [conc_dx_lm]
    block = pellet_clad_mechanical_secondary_subdomain
  []
  [conc_dy_lm]
    block = pellet_clad_mechanical_secondary_subdomain
  []
[]

[AuxVariables]
  [fission_rate]
    block = fuel
    order = CONSTANT
    family = MONOMIAL
  []
  [fluence]
    order = CONSTANT
    family = MONOMIAL
  []
  [burnup]
    block = fuel
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_zz]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 76e6 76.001e6'
    y = '1 1     0'
  []
  [temp_bc]
    type = PiecewiseLinear
    x = '0    76e6  76.001e6 84.641e6 84.6482e6'
    y = '1500 1500  300      300      2073'
  []
  [k_function]
    type = PiecewiseLinear
    x = '0     200e6'
    y = '4e-37 4e-37'
  []
  [d1_function]
    type = ParsedFunction
    expression = 'exp(t/4.5e25)'
  []
  [integral_flux_error]
    type = ParsedFunction
    symbol_names = 'buffer_integral_flux IPyC_integral_flux'
    symbol_values = 'buffer_integral_flux IPyC_integral_flux'
    expression = 'IPyC_integral_flux + buffer_integral_flux'
  []
  [partial_pressure_error]
    type = ParsedFunction
    symbol_names = 'buffer_partial_pressure IPyC_partial_pressure'
    symbol_values = 'buffer_partial_pressure IPyC_partial_pressure'
    expression = 'IPyC_partial_pressure - buffer_partial_pressure'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  generate_output = 'stress_xx stress_yy stress_zz stress_xy stress_yz stress_zx hydrostatic_stress'
  strain = FINITE
  incremental = true
  add_variables = false
  [default]
    block = 'fuel buffer IPyC OPyC'
    eigenstrain_names = 'thermal_strain swelling_strain'
    extra_vector_tags = 'ref'
  []
  [SiC]
    block = 'SiC'
    eigenstrain_names = 'thermal_strain'
    extra_vector_tags = 'ref'
  []
[]

[Kernels]
  [heat_ie]
    type = ADHeatConductionTimeDerivative
    variable = temp
    extra_vector_tags = 'ref'
    block = 'fuel buffer IPyC SiC OPyC'
  []
  [heat]
    type = ADHeatConduction
    variable = temp
    extra_vector_tags = 'ref'
    block = 'fuel buffer IPyC SiC OPyC'
  []
  [heat_source]
    type = ADNeutronHeatSource
    variable = temp
    block = fuel
    energy_per_fission = 3.2e-11 # units of J/fission
    fission_rate = fission_rate
    extra_vector_tags = 'ref'
  []
  [mass_ie]
    type = ADTimeDerivative
    variable = conc
    extra_vector_tags = 'ref'
    block = 'fuel buffer IPyC SiC OPyC'
  []
  [mass]
    type = ADArrheniusDiffusion
    variable = conc
    extra_vector_tags = 'ref'
    block = 'fuel buffer IPyC SiC OPyC'
  []
  [mass_source]
    type = ADBodyForce
    variable = conc
    function = power_history
    value = 1.22e-5 # units of moles/m**3-s
    block = fuel
    extra_vector_tags = 'ref'
  []
  [mass_decay]
    type = Decay
    variable = conc
    radioactive_decay_constant = 7.297e-10 # units:(1/sec)  The constant for Cesium
    block = 'fuel buffer IPyC SiC OPyC'
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [fission_rate]
    type = FissionRateGeneral
    fission_rate_formulation = GENERIC
    variable = fission_rate
    block = fuel
    fission_rate_function = power_history
    value = 3.89e19
    execute_on = timestep_begin
  []
  [fluence]
    type = ADMaterialRealAux
    property = fast_neutron_fluence
    variable = fluence
  []
  [burnup]
    type = ADBurnupAux
    variable = burnup
    block = fuel
    fission_rate = fission_rate
    molecular_weight = 0.270 # units of kg/mole
    execute_on = timestep_begin
    density = ${initial_fuel_density}
  []
  [creep_xx]
    type = ADRankTwoAux
    rank_two_tensor = creep_strain
    variable = creep_xx
    index_i = 0
    index_j = 0
    block = 'buffer IPyC SiC OPyC'
    execute_on = timestep_end
  []
  [creep_yy]
    type = ADRankTwoAux
    rank_two_tensor = creep_strain
    variable = creep_yy
    index_i = 1
    index_j = 1
    block = 'buffer IPyC SiC OPyC'
    execute_on = timestep_end
  []
  [creep_zz]
    type = ADRankTwoAux
    rank_two_tensor = creep_strain
    variable = creep_zz
    index_i = 2
    index_j = 2
    block = 'buffer IPyC SiC OPyC'
    execute_on = timestep_end
  []
[]

[ThermalContactMortar]
  [thermal]
    secondary_variable = temp
    primary_boundary = 15
    secondary_boundary = 17

    initial_moles = initial_moles # coupling to a postprocessor which supplies the initial plenum/gap gas mass
    gas_released = 'fis_gas_released co_production' # coupling to postprocessors which supply the fission gas addition, co addition
    released_gas_types = 'Kr Xe;
                          CO'
    released_fractions = '0.153 0.847;
                          1'
    gap_geometry_type = CYLINDER
    min_gap = 1e-7
    max_gap = 50e-6
    roughness_coef = 0.0
    correct_edge_dropping = true
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 15
    secondary = 17
    model = frictionless
    formulation = mortar
    c_normal = 1.0e8
    correct_edge_dropping = true
  []
[]

[Constraints]
    [cesium_gap_value]
      type = MassSorptionConstraint
      variable = conc_lm
      primary_variable = conc
      primary_boundary = 15
      primary_subdomain = pellet_clad_mechanical_primary_subdomain
      secondary_variable = conc
      secondary_boundary = 17
      secondary_subdomain = pellet_clad_mechanical_secondary_subdomain
      partial_pressure_name = partial_pressure
      epsilon = 1e-4
      correct_edge_dropping = true
    []
    [cesium_gap_flux_x]
      type = MassFluxConstraint
      variable = conc_dx_lm
      primary_variable = conc
      diffusivity_primary = arrhenius_diffusion_coef
      primary_boundary = 15
      primary_subdomain = pellet_clad_mechanical_primary_subdomain
      secondary_variable = conc
      diffusivity_secondary = arrhenius_diffusion_coef
      secondary_boundary = 17
      secondary_subdomain = pellet_clad_mechanical_secondary_subdomain
      component = 0
      epsilon = 1e-5
      correct_edge_dropping = true
    []
    [cesium_gap_flux_y]
      type = MassFluxConstraint
      variable = conc_dy_lm
      primary_variable = conc
      diffusivity_primary = arrhenius_diffusion_coef
      primary_boundary = 15
      primary_subdomain = pellet_clad_mechanical_primary_subdomain
      secondary_variable = conc
      diffusivity_secondary = arrhenius_diffusion_coef
      secondary_boundary = 17
      secondary_subdomain = pellet_clad_mechanical_secondary_subdomain
      component = 1
      epsilon = 1e-5
      correct_edge_dropping = true
    []
  []

[BCs]
  # pin particle along symmetry planes
  [no_disp_x]
    type = ADDirichletBC
    variable = disp_x
    boundary = xzero
    value = 0.0
    extra_vector_tags = 'ref'
  []
  [no_disp_y]
    type = ADDirichletBC
    variable = disp_y
    boundary = yzero
    value = 0.0
    extra_vector_tags = 'ref'
  []
  # fix temperature on free surface
  [freesurf_temp]
    type = ADFunctionDirichletBC
    variable = temp
    boundary = exterior
    function = temp_bc
    extra_vector_tags = 'ref'
  []
  # fix concentration on free surface
  [freesurf_conc]
    type = ADDirichletBC
    variable = conc
    boundary = exterior
    value = 0.0
    extra_vector_tags = 'ref'
  []
  [PlenumPressure] #  apply plenum pressure on clad inner walls and pellet surfaces
    [plenumPressure]
      boundary = BufferGapVol
      initial_pressure = 0
      startup_time = 1.0e4
      R = 8.3145
      output_initial_moles = initial_moles # coupling to post processor to get initial fill gas mass
      temperature = ave_temp_interior # coupling to post processor to get gas temperature approximation
      volume = volumeGas # coupling to post processor to get gas volume
      material_input = 'fis_gas_released co_production' # coupling to post processor to get fission gas added, co added
      output = plenum_pressure # coupling to post processor to output plenum/gap pressure
    []
  []
[]

[Materials]
  [flux]
    type = ADFastNeutronFlux
    calculate_fluence = true
    factor = 5e17
  []

  [fission_gas_release] # Sifgrs fission gas release mode
    type = ADUO2Sifgrs
    block = fuel
    temperature = temp
    fission_rate = fission_rate # coupling to fission_rate aux variable
    grain_radius_const = 5.0e-6
  []

  [fuel_thermal]
    type = ADUO2Thermal
    thermal_conductivity_model = FINK_LUCUTA
    block = fuel
    temperature = temp
    burnup = burnup
    initial_porosity = 0.0
  []
  [fuel_swelling]
    type = ADUO2VolumetricSwellingEigenstrain
    gas_swelling_model_type = MATPRO
    block = fuel
    temperature = temp
    burnup = burnup
    eigenstrain_name = 'swelling_strain'
    initial_fuel_density = ${initial_fuel_density}
  []
  [fuel_stress]
    type = ADComputeFiniteStrainElasticStress
    block = 'fuel'
  []
  [fuel_elasticity]
    type = ADComputeIsotropicElasticityTensor
    block = fuel
    youngs_modulus = 2.2e11
    poissons_ratio = .345
  []
  [fuel_thermal_strain]
    type = ADComputeThermalExpansionEigenstrain
    block = fuel
    thermal_expansion_coeff = 10e-6
    stress_free_temperature = 1500.0
    eigenstrain_name = thermal_strain
    temperature = temp
  []
  [fuel_den]
    type = ADStrainAdjustedDensity
    block = fuel
    strain_free_density = ${initial_fuel_density} # kg/m^3
  []
  [fuel_conc]
    type = ADArrheniusDiffusionCoef
    block = fuel
    d1 = 5.6e-8 # m^2/s
    q1 = 209.0e+3 # J/mol
    d2 = 5.2e-4 # m^2/s
    q2 = 362.0e+3 # J/mol
    gas_constant = 8.3143 # J/K-mol
    temperature = temp
  []

  [buffer_eigenstrain]
    type = ADPyCIrradiationEigenstrain
    block = buffer
    pyc_type = buffer
    eigenstrain_name = 'swelling_strain'
  []
  [buffer_thermal_strain]
    type = ADComputeThermalExpansionEigenstrain
    block = buffer
    thermal_expansion_coeff = 5.65e-6
    stress_free_temperature = 1500.0
    eigenstrain_name = thermal_strain
    temperature = temp
  []
  [buffer_elasticity]
    type = ADComputeIsotropicElasticityTensor
    block = buffer
    youngs_modulus = 2e10
    poissons_ratio = .23
  []
  [buffer_stress]
    type = ADPyCCreep
    block = buffer
    temperature = temp
  []
  [buffer_temp]
    type = ADHeatConductionMaterial
    block = buffer
    thermal_conductivity = 0.5 # J/m-s-K
    specific_heat = 720.0 # J/kg-K
  []
  [buffer_den]
    type = ADStrainAdjustedDensity
    strain_free_density = 1000.0 #kg/m^3
    block = buffer
  []
  [buffer_conc]
    type = ADArrheniusDiffusionCoef
    block = buffer
    d1 = 1.0e-12 # m^2/s
    q1 = 0.0
    d2 = 0.0
    q2 = 0.0
    gas_constant = 8.3143 # J/K-mol
    temperature = temp
  []
  [buffer_partial_pressure]
    type = ADSorptionPartialPressure
    A = 19.33
    B = -47290
    D = 1.518
    E = 4338
    d1 = 3.397
    d2 = 6.15e-4
    unit_scale = 1e3 # convert from mol to mmol
    density = 1000 # convert from mmol/m^3 to mmol/kg, using constant for compatibility with default AD derivative container size
    concentration = conc
    temperature = temp
    block = buffer
    outputs = 'all'
    output_properties = partial_pressure
  []
  [normal_vectors_triso]
    type = NormalVectorsTRISO
    block = 'IPyC OPyC buffer'
  []
  [IPyC_eigenstrain]
    type = ADPyCIrradiationEigenstrain
    block = IPyC
    pyc_type = dense
    eigenstrain_name = 'swelling_strain'
  []
  [IPyC_thermal_strain]
    type = ADComputeThermalExpansionEigenstrain
    block = IPyC
    thermal_expansion_coeff = 5.65e-6
    stress_free_temperature = 1500.0
    eigenstrain_name = thermal_strain
    temperature = temp
  []
  [IPyC_elasticity]
    type = ADComputeIsotropicElasticityTensor
    block = IPyC
    youngs_modulus = 4.74e10
    poissons_ratio = .23
  []
  [IPyC_disp]
    type = ADPyCCreep
    block = 'IPyC OPyC'
    temperature = temp
  []
  [IPyC_temp]
    type = ADHeatConductionMaterial
    block = 'IPyC OPyC'
    thermal_conductivity = 4.0
    specific_heat = 720.0
  []
  [IPyC_den]
    type = ADStrainAdjustedDensity
    block = 'IPyC OPyC'
    strain_free_density = 1900.0
  []
  [IPyC_conc]
    type = ADArrheniusDiffusionCoef
    block = IPyC
    d1 = 6.3e-8
    q1 = 222.0e+3
    d2 = 0.0
    q2 = 0.0
    gas_constant = 8.3143 # J/K-mol
    temperature = temp
  []
  [IPyC_partial_pressure]
    type = ADSorptionPartialPressure
    A = 19.33
    B = -47290
    D = 1.518
    E = 4338
    d1 = 3.397
    d2 = 6.15e-4
    unit_scale = 1e3 # convert from mol to mmol
    density = 1900 # convert from mmol/m^3 to mmol/kg, using constant for compatibility with default AD derivative container size
    concentration = conc
    temperature = temp
    block = IPyC
    outputs = 'all'
    output_properties = partial_pressure
  []
  [SiC_thermal_strain]
    type = ADComputeThermalExpansionEigenstrain
    block = SiC
    thermal_expansion_coeff = 4.9e-6
    stress_free_temperature = 1500.0
    eigenstrain_name = thermal_strain
    temperature = temp
  []
  [SiC_elasticity]
    type = ADComputeIsotropicElasticityTensor
    block = SiC
    youngs_modulus = 3.4e11
    poissons_ratio = .13
  []
  [SiC_creep]
    type = ADMonolithicSiCCreepUpdate
    block = SiC
    temperature = temp
    k_function = k_function
  []
  [SiC_stress]
    type = ADComputeMultipleInelasticStress
    block = SiC
    inelastic_models = 'SiC_creep'
  []
  [SiC_temp]
    type = ADHeatConductionMaterial
    block = SiC
    thermal_conductivity = 13.9 # J/m-s-K
    specific_heat = 620.0 # J/kg-K
  []
  [SiC_den]
    type = ADStrainAdjustedDensity
    strain_free_density = 3180.0 # kg/m^3
    block = SiC
  []
  [SiC_conc]
    type = ADArrheniusDiffusionCoef
    block = SiC
    d1 = 5.5e-14 # m^2/s
    d1_function = d1_function
    d1_function_variable = fluence
    q1 = 125.0e+3 # J/mol
    d2 = 1.6e-2 # m^2/s
    q2 = 514.0e+3 # J/mol
    gas_constant = 8.3143 # J/K-mol
    temperature = temp
  []

  [OPyC_eigenstrain]
    type = ADPyCIrradiationEigenstrain
    block = OPyC
    pyc_type = dense
    eigenstrain_name = 'swelling_strain'
  []
  [OPyC_thermal_strain]
    type = ADComputeThermalExpansionEigenstrain
    block = OPyC
    thermal_expansion_coeff = 5.65e-6
    stress_free_temperature = 1500.0
    eigenstrain_name = thermal_strain
    temperature = temp
  []
  [OPyC_elasticity]
    type = ADComputeIsotropicElasticityTensor
    block = OPyC
    youngs_modulus = 4.74e10
    poissons_ratio = .23
  []
  [OPyC_conc]
    type = ADArrheniusDiffusionCoef
    block = OPyC
    d1 = 6.3e-8 # m^2/s
    q1 = 222.0e+3 # J/mol
    d2 = 0.0
    q2 = 0.0
    gas_constant = 8.3143 # J/K-mol
    temperature = temp
  []
[]

[Dampers]
  [temp]
    type = MaxIncrement
    variable = temp
    max_increment = 50
  []
[]

[Debug]
  show_var_residual_norms = true
[]

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

[Executioner]
  type = Transient
  solve_type = 'NEWTON'
  petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -mat_mffd_err -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu    superlu_dist   1e-5          NONZERO               1e-14'

  snesmf_reuse_base = false

  line_search = 'none'
  nl_rel_tol = 5e-4
  nl_abs_tol = 1e-10
  nl_max_its = 20
  l_max_its = 8

  start_time = 0.0
  end_time = 85.3682e6
  dt = 100

  dtmax = 2e6
  dtmin = 1

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 100
    optimal_iterations = 10
    growth_factor = 1.5
    linear_iteration_ratio = 100
    time_t = '0    76e6  76.001e6 84.641e6 84.6482e6'
    time_dt = '20  20     20      20       20'
  []

  [Predictor]
    type = SimplePredictor
    scale = 0.5
    skip_times_old = '0    76e6  76.001e6 84.641e6 84.6482e6'
  []
[]

[Outputs]
  perf_graph = true
  exodus = true
  [console]
    type = Console
    max_rows = 25
  []
  [csv]
    type = CSV
    sync_times = '100 6308007 75696087'
    sync_only = true
  []
[]

[Postprocessors]
  [Cs_release]
    type = ADSideDiffusiveFluxIntegral
    variable = conc
    diffusivity = arrhenius_diffusion_coef
    boundary = exterior
    execute_on = timestep_end
  []
  [dt]
    type = TimestepSize
    execute_on = timestep_end
  []
  [fis_gas_produced] # fission gas produced (moles)
    type = ADElementIntegralFisGasGeneratedSifgrs
    block = fuel
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  []
  [fis_gas_released] # fission gas released to plenum (moles)
    type = ADElementIntegralFisGasReleasedSifgrs
    block = fuel
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  []
  [volumeTotal]
    type = InternalVolume
    boundary = exterior
    execute_on = 'initial timestep_end'
  []
  [volumeFuel]
    type = InternalVolume
    boundary = fuel
    execute_on = 'initial timestep_end'
  []
  [volumeGas]
    type = InternalVolume
    boundary = BufferGapVol
    # ro = 3.125e-4
    # ri = 2.125e-4
    # vb = 4/3*pi*(ro^3-ri^3) = 8.76e-11
    # buffer density = 1000
    # PyC density    = 1900
    # fill ratio = 10/19
    # vb*10/19 = 4.6e-11
    # Must remove 4.6e-11 m^3 from the volume
    addition = -4.6e-11
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  []
  [volumeBufferShell]
    type = InternalVolume
    boundary = BufferGapVol
    execute_on = 'initial timestep_end'
  []
  [ave_temp_interior]
    type = SideAverageValue
    boundary = BufferGapVol
    variable = temp
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  []

  # Postprocessors for CO production

  [total_fission_rate]
    type = ElementIntegralPower
    variable = temp
    fission_rate = fission_rate
    block = fuel
    energy_per_fission = 1.0
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  []
  [total_fissions]
    type = TimeIntegratedPostprocessor
    value = total_fission_rate
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  []
  [avg_surface_temp]
    type = SideAverageValue
    variable = temp
    boundary = exterior
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  []
  [time_int_surf_temp]
    type = TimeIntegratedPostprocessor
    value = avg_surface_temp
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  []
  [co_production]
    type = CarbonMonoxideProduction
    total_fissions = total_fissions
    time_integrated_triso_temperature = time_int_surf_temp
    initial_enrichment = 0.14029
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  []
  [num_lin_it]
    type = NumLinearIterations
  []
  [num_nonlin_it]
    type = NumNonlinearIterations
  []
  [tot_lin_it]
    type = CumulativeValuePostprocessor
    postprocessor = num_lin_it
  []
  [tot_nonlin_it]
    type = CumulativeValuePostprocessor
    postprocessor = num_nonlin_it
  []
  [alive_time]
    type = PerfGraphData
    section_name = Root
    data_type = TOTAL
  []

  [buffer_integral_flux]
    type = ADSideDiffusiveFluxIntegral
    variable = conc
    boundary = 17
    diffusivity = arrhenius_diffusion_coef
  []
  [IPyC_integral_flux]
    type = ADSideDiffusiveFluxIntegral
    variable = conc
    boundary = 15
    diffusivity = arrhenius_diffusion_coef
  []
  [buffer_partial_pressure]
    type = ADSideAverageMaterialProperty
    property = partial_pressure
    boundary = 17
  []
  [IPyC_partial_pressure]
    type = ADSideAverageMaterialProperty
    property = partial_pressure
    boundary = 15
  []
  [integral_flux_error]
    type = FunctionValuePostprocessor
    function = integral_flux_error
  []
  [partial_pressure_error]
    type = FunctionValuePostprocessor
    function = partial_pressure_error
  []
  [integral_Cs_release]
    type = TimeIntegratedPostprocessor
    value = Cs_release
  []
  [Cs_production]
    type = FunctionValuePostprocessor
    function = power_history
    scale_factor = 1.22e-5 # units of moles/m**3-s
  []
  [time_integral_Cs_production]
    type = TimeIntegratedPostprocessor
    value = Cs_production
  []
  [volumeFuel_initial]
    type = InternalVolume
    boundary = fuel
    execute_on = initial
  []
  [integral_Cs_production]
    type = ParsedPostprocessor
    pp_names = 'time_integral_Cs_production volumeFuel_initial'
    expression = 'time_integral_Cs_production * volumeFuel_initial'
  []
  [Cs_release_fraction]
    type = ParsedPostprocessor
    pp_names = 'integral_Cs_release integral_Cs_production'
    expression = 'integral_Cs_release / integral_Cs_production'
  []
[]

[VectorPostprocessors]
  [temperaturevpp]
    type = SideValueSampler
    boundary = 11
    variable = temp
    sort_by = x
    outputs = 'csv'
    use_displaced_mesh = true
  []
[]

Description
Example Input Syntax
Input Parameters
Input Files
References