Sorption Isotherm Gap Interface

Computes sorption isotherm interfacial conditions.

Overview

The SorptionIsothermGapInterface class is used to enforce sorption as described by a concentration-pressure proportionality across a continuous interface. Users can select between two interface configurations using the interface_configuration parameter: (1) solid-solid and (2) solid-gas.

Solid-Solid Interface Configuration

The solid-solid interface configuration is designed to model sorption between two solid surfaces, which are separated by a gas-filled gap. It enforces the constraint: $(1)var element = document.getElementById("moose-equation-a0304c3f-6828-4ca6-aa01-71463dda0653");katex.render(" P_1 - P_2 = 0, ", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-dd5e4174-3cdd-47e0-8f4d-3cc9f6a943ed");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-30bd752d-b01b-4df6-bc29-9141c9735c40");katex.render("c_i", element, {displayMode:false,throwOnError:false});$ ) and temperature ( $var element = document.getElementById("moose-equation-bb0a9c20-4f64-4417-82bc-9fb4b640c0bf");katex.render("T_i", element, {displayMode:false,throwOnError:false});$ ) conditions at interface side $var element = document.getElementById("moose-equation-70ecafa7-15d6-470b-a6b8-ff882864fe35");katex.render("i", element, {displayMode:false,throwOnError:false});$ . Partial pressures can be calculated using SorptionPartialPressure. Constraining equal partial pressures at the two solid surfaces implicitly assumes that species transport instantaneously across the gas-filled gap.

Solid-Gas Interface Configuration

The solid-gas interface configuration is designed to model sorption between a solid surface and a meshed, gas-filled gap.

note:Boundary and interface kernel orientation conventions

By convention, the solid surface is assumed to be the primary side of the interface, and the gas is assumed to be the neighbor side of the interface. Mesh boundaries must be defined accordingly.

The following constraint, which is derived from the ideal gas law, is applied to solve for the concentration in the gas ( $var element = document.getElementById("moose-equation-201beb88-221a-499b-96f9-e12cb2f3d381");katex.render("c_2", element, {displayMode:false,throwOnError:false});$ ): $(2)var element = document.getElementById("moose-equation-b1949fef-bdbe-4877-9a07-45cb016ce21a");katex.render(" c_2 - \\frac{P_1}{RT_1} = 0, ", element, {displayMode:true,throwOnError:false});$ Note that this interface configuration requires users must supply an additional solid surface temperature ( $var element = document.getElementById("moose-equation-7cfd0a18-d9e7-49b5-b357-af8c006b853c");katex.render("T_1", element, {displayMode:false,throwOnError:false});$ ).The partial pressure can be calculated using SorptionPartialPressure as above.

Mass Flux Balance

Both of the above interface configurations balance the mass flux across the gap using the following constraint: $(3)var element = document.getElementById("moose-equation-aed8afcb-fec8-468c-8bf0-553bb027c8cc");katex.render(" \\mathbf{J_1} \\cdot \\mathbf{n_1} + \\mathbf{J_2} \\cdot \\mathbf{n_2} = 0, ", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-65a0f9cc-be8c-426b-b16e-90ea6621f6c7");katex.render("\\mathbf{J_i}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-d321e48f-afb3-4c0f-b545-a13d36831bbd");katex.render("\\mathbf{n_i}", element, {displayMode:false,throwOnError:false});$ are the mass flux and outward unit normal associated with interface side $var element = document.getElementById("moose-equation-95bbc0c6-c364-40d0-b69c-479f6960e37c");katex.render("i", element, {displayMode:false,throwOnError:false});$ , respectively. This condition implicitly assumes that the surface areas of the two sides of the interface are equal.

$var element = document.getElementById("moose-equation-fd27a54e-043d-4db7-9682-c59b39476b9b");katex.render("\\mathbf{J_i}", element, {displayMode:false,throwOnError:false});$ is further defined as the Fickian diffusion flux according to $(4)var element = document.getElementById("moose-equation-72dbd438-a68d-4bb6-b239-06305ae79c28");katex.render(" \\mathbf{J_i} = -D_i \\mathbf{\\nabla} c_i, ", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-aa3cb0a9-4244-456d-86d4-e88f6849945c");katex.render("D_i", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-1887ec3c-e0f7-473e-a9b7-590e1b2f9f24");katex.render("c_i", element, {displayMode:false,throwOnError:false});$ are the diffusion coefficient and concentration associated with species $var element = document.getElementById("moose-equation-ac8d9d81-f330-4735-b167-00950a8a111c");katex.render("i", element, {displayMode:false,throwOnError:false});$ , respectively.

Two methods are available to enforce flux conditions at the interface. By default, SorptionIsothermGapInterface applies the flux at the primary side of the interface to the residual on the neighbor side of the interface. Using a kernel such as MatDiffusion on the neighbor block adds the flux at the neighbor side of the interface to the residual, forming the full flux balance condition shown in the equation above.

note:Enforcing flux conditions using the default, non-penalty method

When using the default, non-penalty method to enforcing flux conditions at the interface, a kernel such as MatDiffusion must be applied to the neighbor block to form the correct residual.

Optionally, a penalty can be applied within SorptionIsothermGapInterface to directly enforce the full flux condition at the interface. The penalty method adds the correct residual to the neighbor side of the interface without requiring a second kernel, but it may over-constrain the problem under certain conditions.

Example Input Syntax

Example input syntax using the solid-solid interface configuration is shown below.

[InterfaceKernels<<<{"href": "../../syntax/InterfaceKernels/index.html"}>>>]
  [interface]
    type = SorptionIsothermGapInterface<<<{"description": "Computes sorption isotherm interfacial conditions.", "href": "SorptionIsothermGapInterface.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u
    neighbor_var<<<{"description": "The variable on the other side of the interface."}>>> = u
    interface_configuration<<<{"description": "The interface configuration solid_solid solid_gas"}>>> = solid_solid
    partial_pressure_name<<<{"description": "The name of the partial pressure if multiple partial pressure are defined on the same block"}>>> = partial_pressure
    sorption_penalty<<<{"description": "Penalty associated with the concentration imbalance"}>>> = 1e10
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = Block1_Block2
  []
[]

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_1]
    type = SorptionPartialPressure<<<{"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"}>>> = 1
    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_2]
    type = SorptionPartialPressure<<<{"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"}>>> = 2
    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

interface_configurationsolid_solidThe interface configuration solid_solid solid_gas
Default:solid_solid
C++ Type:MooseEnum
Options:solid_solid, solid_gas
Controllable:No
Description:The interface configuration solid_solid solid_gas
neighbor_varThe variable on the other side of the interface.
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The variable on the other side of the interface.
variableThe name of the variable that this residual object operates on
C++ Type:NonlinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The name of the variable that this residual object operates on

Required Parameters

boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
diffusivitydiffusivityThe diffusion coefficient
Default:diffusivity
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:The diffusion coefficient
flux_penalty1Penalty associated with the flux balance
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Penalty associated with the flux balance
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
sorption_penalty1Penalty associated with the concentration imbalance
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Penalty associated with the concentration imbalance
temperatureThe coupled temperature variable
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The coupled temperature variable
use_flux_penaltyFalseWhether to use the penalty formulation to enforce the flux balance
Default:False
C++ Type:bool
Controllable:No
Description:Whether to use the penalty formulation to enforce the flux balance

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

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.)
diag_save_in_var_sideThis parameter must exist if diag_save_in variables are specified and must have the same length as diag_save_in. This vector specifies whether the corresponding aux_var should save-in jacobian contributions from the primary ('p') or secondary side ('s').
C++ Type:MultiMooseEnum
Options:m, s
Controllable:No
Description:This parameter must exist if diag_save_in variables are specified and must have the same length as diag_save_in. This vector specifies whether the corresponding aux_var should save-in jacobian contributions from the primary ('p') or secondary side ('s').
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.)
save_in_var_sideThis parameter must exist if save_in variables are specified and must have the same length as save_in. This vector specifies whether the corresponding aux_var should save-in residual contributions from the primary ('p') or secondary side ('s').
C++ Type:MultiMooseEnum
Options:m, s
Controllable:No
Description:This parameter must exist if save_in variables are specified and must have the same length as save_in. This vector specifies whether the corresponding aux_var should save-in residual contributions from the primary ('p') or secondary side ('s').

Residual And Jacobian Debug Output 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_interface/sorption.i)
(test/tests/sorption_interface/desorption.i)
(examples/TRISO/accident_simulation/triso1D_accident.i)
(assessment/TRISO/validation/AGR-34/SharedFiles/capsule_Sr.i)
(test/tests/sorption_interface/test.i)
(test/tests/sorption_interface/adsorption.i)
(assessment/TRISO/validation/AGR-34/SharedFiles/capsule_Ag.i)
(assessment/TRISO/validation/AGR-34/SharedFiles/capsule_Cs.i)
(test/tests/sorption_interface/sorption_ad.i)

References

No citations exist within this document.

(test/tests/sorption_interface/test.i)

# This test is to verify the implementation of SorptionIsothermGapInterface and its AD counterpart.
# It contains two 2D blocks separated by a continuous interface.
# SorptionIsothermGapInterface is used 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]
  [gen]
    type = GeneratedMeshGenerator
    nx = 10
    ny = 10
    dim = 2
  []
  [block1]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '0.5 1 0'
    input = gen
  []
  [block2]
    type = SubdomainBoundingBoxGenerator
    block_id = 2
    bottom_left = '0.5 0 0'
    top_right = '1 1 0'
    input = block1
  []
  [breakmesh]
    input = block2
    type = BreakMeshByBlockGenerator
    block_pairs = '1 2'
    split_interface = true
    add_interface_on_two_sides = true
  []
[]

[Variables]
  [u]
  []
  [temperature]
  []
[]

[Kernels]
  [u]
    type = MatDiffusion
    variable = u
    diffusivity = diffusivity
  []
  [temp]
    type = HeatConduction
    variable = temperature
  []
[]

[BCs]
  [left_u]
    type = DirichletBC
    value = 1
    variable = u
    boundary = left
  []
  [right_u]
    type = DirichletBC
    value = 1
    variable = u
    boundary = right
  []
  [left_temp]
    type = DirichletBC
    value = 1100
    variable = temperature
    boundary = left
  []
  [right_temp]
    type = DirichletBC
    value = 0
    variable = temperature
    boundary = right
  []
  [block1_2_temp]
    type = DirichletBC
    value = 1000
    variable = temperature
    boundary = Block1_Block2
  []
  [block2_1_temp]
    type = DirichletBC
    value = 900
    variable = temperature
    boundary = Block2_Block1
  []
[]

[InterfaceKernels]
  [interface]
    type = SorptionIsothermGapInterface
    variable = u
    neighbor_var = u
    interface_configuration = solid_solid
    partial_pressure_name = partial_pressure
    sorption_penalty = 1e10
    boundary = Block1_Block2
  []
[]

[Materials]
  [properties_1]
    type = GenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity density'
    prop_values = '1 1 1'
    block = 1
  []
  [partial_pressure_1]
    type = SorptionPartialPressure
    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 = 1
    outputs = 'all'
    output_properties = partial_pressure
  []
  [properties_2]
    type = GenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity density'
    prop_values = '2 2 1'
    block = 2
  []
  [partial_pressure_2]
    type = SorptionPartialPressure
    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 = 2
    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'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  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.4999 0.5 0'
    outputs = csv
  []
  [u_mid_outer]
    type = PointValue
    variable = u
    point = '0.5001 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.4999 0.5 0'
    outputs = csv
  []
  [temperature_mid_outer]
    type = PointValue
    variable = temperature
    point = '0.5001 0.5 0'
    outputs = csv
  []
  [pressure_mid_inner]
    type = PointValue
    variable = partial_pressure
    point = '0.5001 0.5 0'
    outputs = csv
  []
  [pressure_mid_outer]
    type = PointValue
    variable = partial_pressure
    point = '0.5001 0.5 0'
    outputs = csv
  []
  [pressure_mid_error]
    type = FunctionValuePostprocessor
    function = pressure_mid_error
    outputs = console
  []
  [flux_inner] # verify flux preservation
    type = SideDiffusiveFluxIntegral
    variable = u
    boundary = Block1_Block2
    diffusivity = diffusivity
    outputs = csv
  []
  [flux_outer]
    type = SideDiffusiveFluxIntegral
    variable = u
    boundary = Block2_Block1
    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_interface/test.i)

# This test is to verify the implementation of SorptionIsothermGapInterface and its AD counterpart.
# It contains two 2D blocks separated by a continuous interface.
# SorptionIsothermGapInterface is used 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]
  [gen]
    type = GeneratedMeshGenerator
    nx = 10
    ny = 10
    dim = 2
  []
  [block1]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '0.5 1 0'
    input = gen
  []
  [block2]
    type = SubdomainBoundingBoxGenerator
    block_id = 2
    bottom_left = '0.5 0 0'
    top_right = '1 1 0'
    input = block1
  []
  [breakmesh]
    input = block2
    type = BreakMeshByBlockGenerator
    block_pairs = '1 2'
    split_interface = true
    add_interface_on_two_sides = true
  []
[]

[Variables]
  [u]
  []
  [temperature]
  []
[]

[Kernels]
  [u]
    type = MatDiffusion
    variable = u
    diffusivity = diffusivity
  []
  [temp]
    type = HeatConduction
    variable = temperature
  []
[]

[BCs]
  [left_u]
    type = DirichletBC
    value = 1
    variable = u
    boundary = left
  []
  [right_u]
    type = DirichletBC
    value = 1
    variable = u
    boundary = right
  []
  [left_temp]
    type = DirichletBC
    value = 1100
    variable = temperature
    boundary = left
  []
  [right_temp]
    type = DirichletBC
    value = 0
    variable = temperature
    boundary = right
  []
  [block1_2_temp]
    type = DirichletBC
    value = 1000
    variable = temperature
    boundary = Block1_Block2
  []
  [block2_1_temp]
    type = DirichletBC
    value = 900
    variable = temperature
    boundary = Block2_Block1
  []
[]

[InterfaceKernels]
  [interface]
    type = SorptionIsothermGapInterface
    variable = u
    neighbor_var = u
    interface_configuration = solid_solid
    partial_pressure_name = partial_pressure
    sorption_penalty = 1e10
    boundary = Block1_Block2
  []
[]

[Materials]
  [properties_1]
    type = GenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity density'
    prop_values = '1 1 1'
    block = 1
  []
  [partial_pressure_1]
    type = SorptionPartialPressure
    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 = 1
    outputs = 'all'
    output_properties = partial_pressure
  []
  [properties_2]
    type = GenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity density'
    prop_values = '2 2 1'
    block = 2
  []
  [partial_pressure_2]
    type = SorptionPartialPressure
    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 = 2
    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'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  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.4999 0.5 0'
    outputs = csv
  []
  [u_mid_outer]
    type = PointValue
    variable = u
    point = '0.5001 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.4999 0.5 0'
    outputs = csv
  []
  [temperature_mid_outer]
    type = PointValue
    variable = temperature
    point = '0.5001 0.5 0'
    outputs = csv
  []
  [pressure_mid_inner]
    type = PointValue
    variable = partial_pressure
    point = '0.5001 0.5 0'
    outputs = csv
  []
  [pressure_mid_outer]
    type = PointValue
    variable = partial_pressure
    point = '0.5001 0.5 0'
    outputs = csv
  []
  [pressure_mid_error]
    type = FunctionValuePostprocessor
    function = pressure_mid_error
    outputs = console
  []
  [flux_inner] # verify flux preservation
    type = SideDiffusiveFluxIntegral
    variable = u
    boundary = Block1_Block2
    diffusivity = diffusivity
    outputs = csv
  []
  [flux_outer]
    type = SideDiffusiveFluxIntegral
    variable = u
    boundary = Block2_Block1
    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_interface/test.i)

# This test is to verify the implementation of SorptionIsothermGapInterface and its AD counterpart.
# It contains two 2D blocks separated by a continuous interface.
# SorptionIsothermGapInterface is used 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]
  [gen]
    type = GeneratedMeshGenerator
    nx = 10
    ny = 10
    dim = 2
  []
  [block1]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '0.5 1 0'
    input = gen
  []
  [block2]
    type = SubdomainBoundingBoxGenerator
    block_id = 2
    bottom_left = '0.5 0 0'
    top_right = '1 1 0'
    input = block1
  []
  [breakmesh]
    input = block2
    type = BreakMeshByBlockGenerator
    block_pairs = '1 2'
    split_interface = true
    add_interface_on_two_sides = true
  []
[]

[Variables]
  [u]
  []
  [temperature]
  []
[]

[Kernels]
  [u]
    type = MatDiffusion
    variable = u
    diffusivity = diffusivity
  []
  [temp]
    type = HeatConduction
    variable = temperature
  []
[]

[BCs]
  [left_u]
    type = DirichletBC
    value = 1
    variable = u
    boundary = left
  []
  [right_u]
    type = DirichletBC
    value = 1
    variable = u
    boundary = right
  []
  [left_temp]
    type = DirichletBC
    value = 1100
    variable = temperature
    boundary = left
  []
  [right_temp]
    type = DirichletBC
    value = 0
    variable = temperature
    boundary = right
  []
  [block1_2_temp]
    type = DirichletBC
    value = 1000
    variable = temperature
    boundary = Block1_Block2
  []
  [block2_1_temp]
    type = DirichletBC
    value = 900
    variable = temperature
    boundary = Block2_Block1
  []
[]

[InterfaceKernels]
  [interface]
    type = SorptionIsothermGapInterface
    variable = u
    neighbor_var = u
    interface_configuration = solid_solid
    partial_pressure_name = partial_pressure
    sorption_penalty = 1e10
    boundary = Block1_Block2
  []
[]

[Materials]
  [properties_1]
    type = GenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity density'
    prop_values = '1 1 1'
    block = 1
  []
  [partial_pressure_1]
    type = SorptionPartialPressure
    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 = 1
    outputs = 'all'
    output_properties = partial_pressure
  []
  [properties_2]
    type = GenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity density'
    prop_values = '2 2 1'
    block = 2
  []
  [partial_pressure_2]
    type = SorptionPartialPressure
    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 = 2
    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'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  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.4999 0.5 0'
    outputs = csv
  []
  [u_mid_outer]
    type = PointValue
    variable = u
    point = '0.5001 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.4999 0.5 0'
    outputs = csv
  []
  [temperature_mid_outer]
    type = PointValue
    variable = temperature
    point = '0.5001 0.5 0'
    outputs = csv
  []
  [pressure_mid_inner]
    type = PointValue
    variable = partial_pressure
    point = '0.5001 0.5 0'
    outputs = csv
  []
  [pressure_mid_outer]
    type = PointValue
    variable = partial_pressure
    point = '0.5001 0.5 0'
    outputs = csv
  []
  [pressure_mid_error]
    type = FunctionValuePostprocessor
    function = pressure_mid_error
    outputs = console
  []
  [flux_inner] # verify flux preservation
    type = SideDiffusiveFluxIntegral
    variable = u
    boundary = Block1_Block2
    diffusivity = diffusivity
    outputs = csv
  []
  [flux_outer]
    type = SideDiffusiveFluxIntegral
    variable = u
    boundary = Block2_Block1
    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_interface/sorption.i)

# This test is verifies the implementation of SorptionIsothermGapInterface with
# interface_configuration = solid_gas. It contains two solids separated by a gas.
# One solid acts as a source of u. u desorbs into the gap and then onto the
# other solid, which acts as a sink. Temperature increases over time to verify
# sorption behavior and species conservation. For SorptionIsothermGapInterface,
# the solid is always the element and the gas is always the neighbor.

#  _________________________
# | source | gap    | sink   |
#  --------------------------

[GlobalParams]
  order = SECOND
[]

[Mesh]
  coord_type = RZ
  [gen]
    type = GeneratedMeshGenerator
    nx = 150
    xmax = 1.5
    dim = 1
    elem_type = EDGE3
  []
  [source_gap]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '0.75 0 0'
    input = gen
  []
  [gap_sink]
    type = SubdomainBoundingBoxGenerator
    block_id = 2
    bottom_left = '0.75 0 0'
    top_right = '1.5 0 0'
    input = source_gap
  []
  [source]
    type = SubdomainBoundingBoxGenerator
    block_id = 3
    bottom_left = '0 0 0'
    top_right = '0.5 0 0'
    input = gap_sink
  []
  [sink]
    type = SubdomainBoundingBoxGenerator
    block_id = 4
    bottom_left = '1 0 0'
    top_right = '1.5 0 0'
    input = source
  []
  [break]
    type = BreakMeshByBlockGenerator
    block_pairs = '1 3; 2 4'
    split_interface = true
    add_interface_on_two_sides = true
    input = sink
  []
[]

[Variables]
  [u]
  []
[]

[AuxVariables]
  [temperature]
  []
[]

[Kernels]
  [u_dt]
    type = TimeDerivative
    variable = u
  []
  [u]
    type = MatDiffusion
    variable = u
    diffusivity = diffusivity
  []
[]

[ICs]
  [u_source]
    type = FunctionIC
    variable = u
    function = 'if(x < 0.75, 1, 0)'
    block = 3
  []
[]

[AuxKernels]
  [temperature]
    type = FunctionAux
    variable = temperature
    function = '(4000 - 2000 * x) / 2.0 * (t + 0.1)'
    execute_on = 'INITIAL TIMESTEP_BEGIN LINEAR NONLINEAR TIMESTEP_END'
  []
[]

[InterfaceKernels]
  [u_solid_gas]
    type = SorptionIsothermGapInterface
    variable = u
    neighbor_var = u
    interface_configuration = solid_gas
    temperature = temperature
    partial_pressure_name = partial_pressure
    sorption_penalty = 1e10
    boundary = 'Block3_Block1 Block4_Block2'
  []
[]

[Materials]
  [properties_solid]
    type = GenericConstantMaterial
    prop_names = 'density diffusivity'
    prop_values = '1 0.1'
    block = '3 4'
  []
  [properties_gas]
    type = GenericConstantMaterial
    prop_names = diffusivity
    prop_values = 1
    block = '1 2'
  []
  [partial_pressure_source]
    type = SorptionPartialPressure
    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 = 3
    outputs = 'all'
    output_properties = partial_pressure
  []
  [partial_pressure_sink]
    type = SorptionPartialPressure
    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 = 4
    outputs = 'all'
    output_properties = partial_pressure
  []
[]

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

[Executioner]
  type = Transient
  solve_type = NEWTON
  automatic_scaling = true
  compute_scaling_once = false

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  line_search = none

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

  dt = 0.1
  end_time = 2
[]

[Functions]
  [source_surface_u_exact]
    type = ParsedFunction
    symbol_names = 'P T R'
    symbol_values = 'source_surface_pressure source_surface_temperature 8.31446261815324'
    expression = '1.0 / R / T * P'
  []
  [sink_surface_u_exact]
    type = ParsedFunction
    symbol_names = 'P T R'
    symbol_values = 'sink_surface_pressure sink_surface_temperature 8.31446261815324'
    expression = '1.0 / R / T * P'
  []
[]

[Postprocessors]
  [source_surface_pressure]
    type = SideAverageMaterialProperty
    property = partial_pressure
    boundary = Block3_Block1
  []
  [source_surface_temperature]
    type = SideAverageValue
    variable = temperature
    boundary = Block3_Block1
  []
  [source_surface_u]
    type = SideAverageValue
    variable = u
    boundary = Block1_Block3
  []
  [source_surface_u_exact]
    type = FunctionValuePostprocessor
    function = source_surface_u_exact
  []

  [sink_surface_pressure]
    type = SideAverageMaterialProperty
    property = partial_pressure
    boundary = Block4_Block2
  []
  [sink_surface_temperature]
    type = SideAverageValue
    variable = temperature
    boundary = Block4_Block2
  []
  [sink_surface_u]
    type = SideAverageValue
    variable = u
    boundary = Block2_Block4
  []
  [sink_surface_u_exact]
    type = FunctionValuePostprocessor
    function = sink_surface_u_exact
  []

  [u_source_tot]
    type = ElementIntegralVariablePostprocessor
    variable = u
    block = 3
  []
  [u_gap_tot]
    type = ElementIntegralVariablePostprocessor
    variable = u
    block = '1 2'
  []
  [u_sink_tot]
    type = ElementIntegralVariablePostprocessor
    variable = u
    block = 4
  []
  [u_tot]
    type = ElementIntegralVariablePostprocessor
    variable = u
    block = '1 2 3 4'
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

[Debug]
  show_var_residual_norms = true
[]

(test/tests/sorption_interface/desorption.i)

# This test is verifies the implementation of SorptionIsothermGapInterface with
# interface_configuration = solid_gas. It contains a solid having a single
# interface with a gas. The solid acts as a source of u. u desorbs into the gap.
# Temperature increases over time to verify desorption behavior and species
# conservation. For SorptionIsothermGapInterface, the solid is always the
# element and the gas is always the neighbor.

#  _________________
# | source | gap    |
#  -----------------

[GlobalParams]
  order = SECOND
[]

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    nx = 10
    ny = 2
    ymax = 0.2
    dim = 2
    elem_type = QUAD8
  []
  [source]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '0.5 0.2 0'
    input = gen
  []
  [gap]
    type = SubdomainBoundingBoxGenerator
    block_id = 2
    bottom_left = '0.5 0 0'
    top_right = '1 0.2 0'
    input = source
  []
  [break]
    type = BreakMeshByBlockGenerator
    block_pairs = '1 2'
    split_interface = true
    add_interface_on_two_sides = true
    input = gap
  []
[]

[Variables]
  [u]
  []
[]

[AuxVariables]
  [temperature]
  []
[]

[Kernels]
  [u_dt]
    type = TimeDerivative
    variable = u
  []
  [u]
    type = MatDiffusion
    variable = u
    diffusivity = diffusivity
  []
[]

[ICs]
  [u_source]
    type = ConstantIC
    variable = u
    value = 1
    block = 1
  []
[]

[AuxKernels]
  [temperature]
    type = FunctionAux
    variable = temperature
    function = '(4000 - 1000 * x) / 2.0 * (t + 0.1)'
    execute_on = 'INITIAL TIMESTEP_BEGIN LINEAR NONLINEAR TIMESTEP_END'
  []
[]

[InterfaceKernels]
  [u_source_gap]
    type = SorptionIsothermGapInterface
    variable = u
    neighbor_var = u
    interface_configuration = solid_gas
    temperature = temperature
    partial_pressure_name = partial_pressure
    sorption_penalty = 1e10
    boundary = Block1_Block2
  []
[]

[Materials]
  [properties_solid]
    type = GenericConstantMaterial
    prop_names = 'density diffusivity'
    prop_values = '1 0.1'
    block = 1
  []
  [properties_gas]
    type = GenericConstantMaterial
    prop_names = diffusivity
    prop_values = 1
    block = 2
  []
  [partial_pressure]
    type = SorptionPartialPressure
    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 = 1
    outputs = 'all'
    output_properties = partial_pressure
  []
[]

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

[Executioner]
  type = Transient
  solve_type = NEWTON
  automatic_scaling = true
  compute_scaling_once = false

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  line_search = none

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

  dt = 0.1
  end_time = 2
[]

[Functions]
  [surface_u_exact]
    type = ParsedFunction
    symbol_names = 'P T R'
    symbol_values = 'surface_pressure surface_temperature 8.31446261815324'
    expression = '1.0 / R / T * P'
  []
[]

[Postprocessors]
  [surface_pressure]
    type = SideAverageMaterialProperty
    property = partial_pressure
    boundary = Block1_Block2
  []
  [surface_temperature]
    type = SideAverageValue
    variable = temperature
    boundary = Block1_Block2
  []
  [surface_u]
    type = SideAverageValue
    variable = u
    boundary = Block2_Block1
  []
  [surface_u_exact]
    type = FunctionValuePostprocessor
    function = surface_u_exact
  []

  [u_source_tot]
    type = ElementIntegralVariablePostprocessor
    variable = u
    block = 1
  []
  [u_gap_tot]
    type = ElementIntegralVariablePostprocessor
    variable = u
    block = 2
  []
  [u_tot]
    type = ElementIntegralVariablePostprocessor
    variable = u
    block = '1 2'
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

[Debug]
  show_var_residual_norms = true
[]

(examples/TRISO/accident_simulation/triso1D_accident.i)

# This example is 1D spherical 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 an interface kernel to model gap mass transfer.
# 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

[GlobalParams]
  order = SECOND
  family = LAGRANGE
  displacements = disp_x
  flux_conversion_factor = 0.85
[]

[Mesh]
  coord_type = RSPHERICAL
  [gen] # exclude gap to establish buffer-IPyC neighbor relationships for the sorption interface kernel
    type = TRISO1DMeshGenerator
    elem_type = EDGE3
    coordinates = '0 2.125e-4 3.125e-4 3.525e-4 3.875e-4 4.275e-4'
    mesh_density = '10 5 2 2 2'
    block_names = 'fuel buffer IPyC SiC OPyC'
  []
  [break] # create gap between buffer and IPyC to model mechanical and thermal contact
    type = BreakMeshByBlockGenerator
    input = gen
    block_pairs = 'buffer IPyC'
    split_interface = true
    add_interface_on_two_sides = true
  []
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
[]

[Variables]
  [disp_x]
  []
  [temp]
    initial_condition = 1500.0
  []
  [conc]
    initial_condition = 0.0
  []
[]

[AuxVariables]
  [fission_rate]
    block = fuel
    order = CONSTANT
    family = MONOMIAL
  []
  [fluence]
    order = CONSTANT
    family = MONOMIAL
  []
  [burnup]
    block = fuel
    order = CONSTANT
    family = MONOMIAL
  []
  [gap_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_zz]
    order = CONSTANT
    family = MONOMIAL
  []
  [gap_HTC]
    order = CONSTANT
    family = MONOMIAL
  []
  [gap_distance]
    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 = HeatConductionTimeDerivative
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat]
    type = HeatConduction
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_source]
    type = NeutronHeatSource
    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 = TimeDerivative
    variable = conc
    extra_vector_tags = 'ref'
  []
  [mass]
    type = ArrheniusDiffusion
    variable = conc
    extra_vector_tags = 'ref'
  []
  [mass_source]
    type = BodyForce
    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
    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 = MaterialRealAux
    property = fast_neutron_fluence
    variable = fluence
  []
  [burnup]
    type = BurnupAux
    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 = RankTwoAux
    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 = RankTwoAux
    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 = RankTwoAux
    rank_two_tensor = creep_strain
    variable = creep_zz
    index_i = 2
    index_j = 2
    block = 'buffer IPyC SiC OPyC'
    execute_on = timestep_end
  []
  [conductance]
    type = MaterialRealAux
    property = gap_conductance
    variable = gap_cond
    boundary = buffer_IPyC
    execute_on = linear
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = IPyC_buffer
    secondary = buffer_IPyC
    penalty = 1e5
    model = frictionless
    formulation = penalty
  []
[]

[ThermalContact]
  [thermal_contact]
    type = GasGapHeatTransfer
    variable = temp
    primary = IPyC_buffer
    secondary = buffer_IPyC
    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 = SPHERE
    tangential_tolerance = 1e-6
    roughness_coef = 0.0
    quadrature = true
  []
[]

[InterfaceKernels]
  [cesium_gap]
    type = SorptionIsothermGapInterface
    variable = conc
    neighbor_var = conc
    partial_pressure_name = partial_pressure
    sorption_penalty = 1e5
    diffusivity = arrhenius_diffusion_coef
    use_flux_penalty = true
    flux_penalty = 1e3
    boundary = buffer_IPyC
    extra_vector_tags = 'ref'
  []
[]

[BCs]
  # pin particle along symmetry planes
  [no_disp_x]
    type = DirichletBC
    variable = disp_x
    boundary = xzero
    value = 0.0
    extra_vector_tags = 'ref'
  []
  # fix temperature on free surface
  [freesurf_temp]
    type = FunctionDirichletBC
    variable = temp
    boundary = exterior
    function = temp_bc
    extra_vector_tags = 'ref'
  []
  # fix concentration on free surface
  [freesurf_conc]
    type = DirichletBC
    variable = conc
    boundary = exterior
    value = 0.0
    extra_vector_tags = 'ref'
  []
  [PlenumPressure] #  apply plenum pressure on clad inner walls and pellet surfaces
    [plenumPressure]
      boundary = 'buffer_IPyC IPyC_buffer'
      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 = FastNeutronFlux
    calculate_fluence = true
    factor = 5e17
  []

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

  [fuel_thermal]
    type = UO2Thermal
    thermal_conductivity_model = FINK_LUCUTA
    block = fuel
    temperature = temp
    burnup = burnup
    initial_porosity = 0.0
  []
  [fuel_swelling]
    type = UO2VolumetricSwellingEigenstrain
    gas_swelling_model_type = MATPRO
    block = fuel
    temperature = temp
    burnup = burnup
    eigenstrain_name = 'swelling_strain'
    initial_fuel_density = ${initial_fuel_density}
  []
  [fuel_stress]
    type = ComputeFiniteStrainElasticStress
    block = 'fuel'
  []
  [fuel_elasticity]
    type = ComputeIsotropicElasticityTensor
    block = fuel
    youngs_modulus = 2.2e11
    poissons_ratio = .345
  []
  [fuel_thermal_strain]
    type = ComputeThermalExpansionEigenstrain
    block = fuel
    thermal_expansion_coeff = 10e-6
    stress_free_temperature = 1500.0
    eigenstrain_name = thermal_strain
    temperature = temp
  []
  [fuel_den]
    type = StrainAdjustedDensity
    block = fuel
    strain_free_density = ${initial_fuel_density} # kg/m^3
  []
  [fuel_conc]
    type = ArrheniusDiffusionCoef
    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 = PyCIrradiationEigenstrain
    block = buffer
    pyc_type = buffer
    eigenstrain_name = 'swelling_strain'
  []
  [buffer_thermal_strain]
    type = ComputeThermalExpansionEigenstrain
    block = buffer
    thermal_expansion_coeff = 5.65e-6
    stress_free_temperature = 1500.0
    eigenstrain_name = thermal_strain
    temperature = temp
  []
  [buffer_elasticity]
    type = ComputeIsotropicElasticityTensor
    block = buffer
    youngs_modulus = 2e10
    poissons_ratio = .23
  []
  [buffer_stress]
    type = PyCCreep
    block = buffer
    temperature = temp
  []
  [buffer_temp]
    type = HeatConductionMaterial
    block = buffer
    thermal_conductivity = 0.5 # J/m-s-K
    specific_heat = 720.0 # J/kg-K
  []
  [buffer_den]
    type = StrainAdjustedDensity
    strain_free_density = 1000.0 #kg/m^3
    block = buffer
  []
  [buffer_conc]
    type = ArrheniusDiffusionCoef
    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 = SorptionPartialPressure
    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 = density # convert from mmol/m^3 to mmol/kg
    concentration = conc
    temperature = temp
    block = 'buffer IPyC'
    outputs = 'all'
    output_properties = partial_pressure
  []
  [normal_vectors_triso]
    type = NormalVectorsTRISO
    block = 'IPyC OPyC buffer'
  []
  [IPyC_eigenstrain]
    type = PyCIrradiationEigenstrain
    block = IPyC
    pyc_type = dense
    eigenstrain_name = 'swelling_strain'
  []
  [IPyC_thermal_strain]
    type = ComputeThermalExpansionEigenstrain
    block = IPyC
    thermal_expansion_coeff = 5.65e-6
    stress_free_temperature = 1500.0
    eigenstrain_name = thermal_strain
    temperature = temp
  []
  [IPyC_elasticity]
    type = ComputeIsotropicElasticityTensor
    block = IPyC
    youngs_modulus = 4.74e10
    poissons_ratio = .23
  []
  [IPyC_disp]
    type = PyCCreep
    block = 'IPyC OPyC'
    temperature = temp
  []
  [IPyC_temp]
    type = HeatConductionMaterial
    block = 'IPyC OPyC'
    thermal_conductivity = 4.0
    specific_heat = 720.0
  []
  [IPyC_den]
    type = StrainAdjustedDensity
    block = 'IPyC OPyC'
    strain_free_density = 1900.0
  []
  [IPyC_conc]
    type = ArrheniusDiffusionCoef
    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
  []
  [SiC_thermal_strain]
    type = ComputeThermalExpansionEigenstrain
    block = SiC
    thermal_expansion_coeff = 4.9e-6
    stress_free_temperature = 1500.0
    eigenstrain_name = thermal_strain
    temperature = temp
  []
  [SiC_elasticity]
    type = ComputeIsotropicElasticityTensor
    block = SiC
    youngs_modulus = 3.4e11
    poissons_ratio = .13
  []
  [SiC_creep]
    type = MonolithicSiCCreepUpdate
    block = SiC
    temperature = temp
    k_function = k_function
  []
  [SiC_stress]
    type = ComputeMultipleInelasticStress
    block = SiC
    tangent_operator = elastic
    inelastic_models = 'SiC_creep'
  []
  [SiC_temp]
    type = HeatConductionMaterial
    block = SiC
    thermal_conductivity = 13.9 # J/m-s-K
    specific_heat = 620.0 # J/kg-K
  []
  [SiC_den]
    type = StrainAdjustedDensity
    strain_free_density = 3180.0 # kg/m^3
    block = SiC
  []
  [SiC_conc]
    type = ArrheniusDiffusionCoef
    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 = PyCIrradiationEigenstrain
    block = OPyC
    pyc_type = dense
    eigenstrain_name = 'swelling_strain'
  []
  [OPyC_thermal_strain]
    type = ComputeThermalExpansionEigenstrain
    block = OPyC
    thermal_expansion_coeff = 5.65e-6
    stress_free_temperature = 1500.0
    eigenstrain_name = thermal_strain
    temperature = temp
  []
  [OPyC_elasticity]
    type = ComputeIsotropicElasticityTensor
    block = OPyC
    youngs_modulus = 4.74e10
    poissons_ratio = .23
  []
  [OPyC_conc]
    type = ArrheniusDiffusionCoef
    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 = 'PJFNK'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'

  line_search = 'none'

  nl_rel_tol = 5e-4
  nl_abs_tol = 1e-9
  nl_max_its = 15

  l_tol = 1e-3
  l_max_its = 50

  start_time = 0.0
  end_time = 85.3682e6
  dt = 100

  dtmax = 2e6
  dtmin = 1

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 100
    optimal_iterations = 6
    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 = 1
    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 = SideIntegralMassFlux
    variable = conc
    boundary = exterior
    execute_on = timestep_end
  []
  [dt]
    type = TimestepSize
    execute_on = timestep_end
  []
  [fis_gas_produced] # fission gas produced (moles)
    type = ElementIntegralFisGasGeneratedSifgrs
    block = fuel
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  []
  [fis_gas_released] # fission gas released to plenum (moles)
    type = ElementIntegralFisGasReleasedSifgrs
    block = fuel
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  []
  [volumeTotal]
    type = InternalVolume
    boundary = exterior
    scale_factor = -1
    execute_on = 'initial timestep_end'
  []
  [volumeFuel]
    type = InternalVolume
    boundary = fuel_outer_boundary
    scale_factor = -1
    execute_on = 'initial timestep_end'
  []
  [volumeGas]
    type = InternalVolume
    boundary = 'buffer_IPyC IPyC_buffer'
    # 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 = 'buffer_IPyC IPyC_buffer'
    execute_on = 'initial timestep_end'
  []
  [ave_temp_interior]
    type = SideAverageValue
    boundary = 'buffer_IPyC IPyC_buffer'
    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 = SideDiffusiveFluxIntegral
    variable = conc
    boundary = buffer_IPyC
    diffusivity = arrhenius_diffusion_coef
  []
  [IPyC_integral_flux]
    type = SideDiffusiveFluxIntegral
    variable = conc
    boundary = IPyC_buffer
    diffusivity = arrhenius_diffusion_coef
  []
  [buffer_partial_pressure]
    type = SideAverageMaterialProperty
    property = partial_pressure
    boundary = buffer_IPyC
  []
  [IPyC_partial_pressure]
    type = SideAverageMaterialProperty
    property = partial_pressure
    boundary = IPyC_buffer
  []
  [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_outer_boundary
    scale_factor = -1
    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
  []
[]

(assessment/TRISO/validation/AGR-34/SharedFiles/capsule_Sr.i)

[InterfaceKernels]
  [fuel_matrix]
    type = SorptionIsothermGapInterface
    variable = conc_Sr
    neighbor_var = conc_Sr
    partial_pressure_name = partial_pressure
    sorption_penalty = 1e5
    diffusivity = arrhenius_diffusion_coef_Sr
    use_flux_penalty = true
    flux_penalty = 1e3
    boundary = 'DRV_IR'
  []
  [matrix_graphite]
    type = SorptionIsothermGapInterface
    variable = conc_Sr
    neighbor_var = conc_Sr
    partial_pressure_name = partial_pressure
    sorption_penalty = 1e5
    diffusivity = arrhenius_diffusion_coef_Sr
    use_flux_penalty = true
    flux_penalty = 1e3
    boundary = 'IR_OR'
  []
  [graphite_sink]
    type = SorptionIsothermGapInterface
    variable = conc_Sr
    neighbor_var = conc_Sr
    partial_pressure_name = partial_pressure
    sorption_penalty = 1e5
    diffusivity = arrhenius_diffusion_coef_Sr
    use_flux_penalty = true
    flux_penalty = 1e3
    boundary = 'OR_SR'
  []
[]

[Variables]
  [conc_Sr]
    scaling = 1e12 #1e18
  []
[]

[AuxVariables]
  [Sr_diff_coef]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [mass_Sr_dt]
    type = TimeDerivative
    variable = conc_Sr
    extra_vector_tags = 'ref'
  []
  [mass_Sr]
    type = ArrheniusDiffusion
    variable = conc_Sr
    arrhenius_prpty_name = arrhenius_diffusion_coef_Sr
    extra_vector_tags = 'ref'
  []
  [conc_source_dtf]
    type = BodyForce
    variable = conc_Sr
    block = 'DTF'
    function = DTF_MassIsotope
    extra_vector_tags = 'ref'
  []
  [conc_source_drv]
    type = BodyForce
    variable = conc_Sr
    block = 'DRV'
    function = DRV_MassIsotope
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [Sr_diff_coef]
    type = MaterialRealAux
    variable = Sr_diff_coef
    property = arrhenius_diffusion_coef_Sr
    execute_on = timestep_end
  []
[]

[BCs]
  [freesurf_conc_Sr]
    type = DirichletBC
    variable = conc_Sr
    boundary = 'right'
    value = 0.0
    extra_vector_tags = 'ref'
  []
[]

[Materials]
  [matrix_conc_Sr]
    type = ArrheniusDiffusionCoef
    d1 = 1.00e-2 # m^2/s
    q1 = 303000 # J/mol
    temperature = temperature
    arrhenius_prpty_name = arrhenius_diffusion_coef_Sr
  []
  [graphite_conc_Sr]
    type = ArrheniusDiffusionCoef
    d1 = 1.70e-2 # m^2/s
    q1 = 268000 # J/mol
    temperature = temperature
    arrhenius_prpty_name = arrhenius_diffusion_coef_Sr
  []
  [Sr_matrix_partial_pressure]
    type = SorptionPartialPressure
    A = 54.3
    B = -149000
    D = -8.52
    E = 28500
    d1 = 3.13
    d2 = 0
    unit_scale = 1e3 # convert from mol to mmol
    density = density # convert from mmol/m^3 to mmol/kg
    concentration = conc_Sr
    temperature = temperature
    output_properties = partial_pressure
  []
  [Sr_graphite_partial_pressure]
    type = SorptionPartialPressure
    A = 19.4
    B = -40100
    D = -0.32
    E = 4090
    d1 = -2.12
    d2 = 0
    unit_scale = 1e3 # convert from mol to mmol
    density = density # convert from mmol/m^3 to mmol/kg
    concentration = conc_Sr
    temperature = temperature
    output_properties = partial_pressure
  []
[]

[Postprocessors]
  [release_Sr_inc_pebble]
    type = SideIntegralMassFlux
    variable = conc_Sr
    boundary = right
    arrhenius_prpty_name = arrhenius_diffusion_coef_Sr
    execute_on = 'initial timestep_end'
  []
  [released_Sr]
    type = TimeIntegratedPostprocessor
    value = release_Sr_inc_pebble
    execute_on = 'initial timestep_end'
  []
  [retained_Sr]
    type = ElementIntegralVariablePostprocessor
    variable = conc_Sr
  []

  # check partial pressure equality and integral flux conservation from DRV to IR
  [flux_DRV_outer]
    type = SideDiffusiveFluxIntegral
    boundary = 'DRV_IR'
    variable = conc_Sr
    diffusivity = arrhenius_diffusion_coef_Sr
  []
  [flux_IR_inner]
    type = SideDiffusiveFluxIntegral
    boundary = 'IR_DRV'
    variable = conc_Sr
    diffusivity = arrhenius_diffusion_coef_Sr
  []
  [flux_DRV_IR_error]
    type = FunctionValuePostprocessor
    function = flux_DRV_IR_error
  []

  # check partial pressure equality and integral flux conservation from IR to OR
  [flux_IR_outer]
    type = SideDiffusiveFluxIntegral
    boundary = 'IR_OR'
    variable = conc_Sr
    diffusivity = arrhenius_diffusion_coef_Sr
  []
  [flux_OR_inner]
    type = SideDiffusiveFluxIntegral
    boundary = 'OR_IR'
    variable = conc_Sr
    diffusivity = arrhenius_diffusion_coef_Sr
  []
  [flux_IR_OR_error]
    type = FunctionValuePostprocessor
    function = flux_IR_OR_error
  []

  # check partial pressure equality and integral flux conservation from OR to SR
  [flux_OR_outer]
    type = SideDiffusiveFluxIntegral
    boundary = 'OR_SR'
    variable = conc_Sr
    diffusivity = arrhenius_diffusion_coef_Sr
  []
  [flux_SR_inner]
    type = SideDiffusiveFluxIntegral
    boundary = 'SR_OR'
    variable = conc_Sr
    diffusivity = arrhenius_diffusion_coef_Sr
  []
  [flux_OR_SR_error]
    type = FunctionValuePostprocessor
    function = flux_OR_SR_error
  []
[]

[VectorPostprocessors]
  [Sr_conc]
    type = LineValueSampler
    variable = conc_Sr
    start_point = '0.0 0 0.0'
    num_points = 1000
    use_displaced_mesh = false
    execute_on = final
    sort_by = x
    outputs = line_plot_Sr
  []
[]

[Outputs]
  [line_plot_Sr]
    type = CSV
    execute_on = 'FINAL'
    time_step_interval =  1
    create_final_symlink = true
    file_base = capsule_radial
  []
[]

(test/tests/sorption_interface/test.i)

# This test is to verify the implementation of SorptionIsothermGapInterface and its AD counterpart.
# It contains two 2D blocks separated by a continuous interface.
# SorptionIsothermGapInterface is used 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]
  [gen]
    type = GeneratedMeshGenerator
    nx = 10
    ny = 10
    dim = 2
  []
  [block1]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '0.5 1 0'
    input = gen
  []
  [block2]
    type = SubdomainBoundingBoxGenerator
    block_id = 2
    bottom_left = '0.5 0 0'
    top_right = '1 1 0'
    input = block1
  []
  [breakmesh]
    input = block2
    type = BreakMeshByBlockGenerator
    block_pairs = '1 2'
    split_interface = true
    add_interface_on_two_sides = true
  []
[]

[Variables]
  [u]
  []
  [temperature]
  []
[]

[Kernels]
  [u]
    type = MatDiffusion
    variable = u
    diffusivity = diffusivity
  []
  [temp]
    type = HeatConduction
    variable = temperature
  []
[]

[BCs]
  [left_u]
    type = DirichletBC
    value = 1
    variable = u
    boundary = left
  []
  [right_u]
    type = DirichletBC
    value = 1
    variable = u
    boundary = right
  []
  [left_temp]
    type = DirichletBC
    value = 1100
    variable = temperature
    boundary = left
  []
  [right_temp]
    type = DirichletBC
    value = 0
    variable = temperature
    boundary = right
  []
  [block1_2_temp]
    type = DirichletBC
    value = 1000
    variable = temperature
    boundary = Block1_Block2
  []
  [block2_1_temp]
    type = DirichletBC
    value = 900
    variable = temperature
    boundary = Block2_Block1
  []
[]

[InterfaceKernels]
  [interface]
    type = SorptionIsothermGapInterface
    variable = u
    neighbor_var = u
    interface_configuration = solid_solid
    partial_pressure_name = partial_pressure
    sorption_penalty = 1e10
    boundary = Block1_Block2
  []
[]

[Materials]
  [properties_1]
    type = GenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity density'
    prop_values = '1 1 1'
    block = 1
  []
  [partial_pressure_1]
    type = SorptionPartialPressure
    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 = 1
    outputs = 'all'
    output_properties = partial_pressure
  []
  [properties_2]
    type = GenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity density'
    prop_values = '2 2 1'
    block = 2
  []
  [partial_pressure_2]
    type = SorptionPartialPressure
    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 = 2
    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'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  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.4999 0.5 0'
    outputs = csv
  []
  [u_mid_outer]
    type = PointValue
    variable = u
    point = '0.5001 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.4999 0.5 0'
    outputs = csv
  []
  [temperature_mid_outer]
    type = PointValue
    variable = temperature
    point = '0.5001 0.5 0'
    outputs = csv
  []
  [pressure_mid_inner]
    type = PointValue
    variable = partial_pressure
    point = '0.5001 0.5 0'
    outputs = csv
  []
  [pressure_mid_outer]
    type = PointValue
    variable = partial_pressure
    point = '0.5001 0.5 0'
    outputs = csv
  []
  [pressure_mid_error]
    type = FunctionValuePostprocessor
    function = pressure_mid_error
    outputs = console
  []
  [flux_inner] # verify flux preservation
    type = SideDiffusiveFluxIntegral
    variable = u
    boundary = Block1_Block2
    diffusivity = diffusivity
    outputs = csv
  []
  [flux_outer]
    type = SideDiffusiveFluxIntegral
    variable = u
    boundary = Block2_Block1
    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_interface/adsorption.i)

# This test is verifies the implementation of SorptionIsothermGapInterface with
# interface_configuration = solid_gas. It contains a solid having a single
# interface with a gas. The gas acts as a source of u. u adsorbs onto the sink.
# Temperature increases over time to verify adsorption behavior and species
# conservation. For SorptionIsothermGapInterface, the solid is always the
# element and the gas is always the neighbor.

#  _________________
# | sink   | gap    |
#  -----------------

[GlobalParams]
  order = SECOND
[]

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    nx = 10
    ny = 2
    ymax = 0.2
    dim = 2
    elem_type = QUAD8
  []
  [sink]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '0.5 0.2 0'
    input = gen
  []
  [gap]
    type = SubdomainBoundingBoxGenerator
    block_id = 2
    bottom_left = '0.5 0 0'
    top_right = '1 0.2 0'
    input = sink
  []
  [break]
    type = BreakMeshByBlockGenerator
    block_pairs = '1 2'
    split_interface = true
    add_interface_on_two_sides = true
    input = gap
  []
[]

[Variables]
  [u]
  []
[]

[AuxVariables]
  [temperature]
  []
[]

[Kernels]
  [u_dt]
    type = TimeDerivative
    variable = u
  []
  [u]
    type = MatDiffusion
    variable = u
    diffusivity = diffusivity
  []
[]

[ICs]
  [u_gap]
    type = ConstantIC
    variable = u
    value = 1
    block = 2
  []
[]

[AuxKernels]
  [temperature]
    type = FunctionAux
    variable = temperature
    function = '(4000 - 1000 * x) / 2.0 * (2.0 - t + 0.1)'
    execute_on = 'INITIAL TIMESTEP_BEGIN LINEAR NONLINEAR TIMESTEP_END'
  []
[]

[InterfaceKernels]
  [u_sink_gap]
    type = SorptionIsothermGapInterface
    variable = u
    neighbor_var = u
    interface_configuration = solid_gas
    temperature = temperature
    partial_pressure_name = partial_pressure
    sorption_penalty = 1e10
    boundary = Block1_Block2
  []
[]

[Materials]
  [properties_solid]
    type = GenericConstantMaterial
    prop_names = 'density diffusivity'
    prop_values = '1 0.1'
    block = 1
  []
  [properties_gas]
    type = GenericConstantMaterial
    prop_names = diffusivity
    prop_values = 1
    block = 2
  []
  [partial_pressure]
    type = SorptionPartialPressure
    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 = 1
    outputs = 'all'
    output_properties = partial_pressure
  []
[]

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

[Executioner]
  type = Transient
  solve_type = NEWTON
  automatic_scaling = true
  compute_scaling_once = false

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  line_search = none

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

  dt = 0.1
  end_time = 2
[]

[Functions]
  [surface_u_exact]
    type = ParsedFunction
    symbol_names = 'P T R'
    symbol_values = 'surface_pressure surface_temperature 8.31446261815324'
    expression = '1.0 / R / T * P'
  []
[]

[Postprocessors]
  [surface_pressure]
    type = SideAverageMaterialProperty
    property = partial_pressure
    boundary = Block1_Block2
  []
  [surface_temperature]
    type = SideAverageValue
    variable = temperature
    boundary = Block1_Block2
  []
  [surface_u]
    type = SideAverageValue
    variable = u
    boundary = Block2_Block1
  []
  [surface_u_exact]
    type = FunctionValuePostprocessor
    function = surface_u_exact
  []

  [u_sink_tot]
    type = ElementIntegralVariablePostprocessor
    variable = u
    block = 1
  []
  [u_gap_tot]
    type = ElementIntegralVariablePostprocessor
    variable = u
    block = 2
  []
  [u_tot]
    type = ElementIntegralVariablePostprocessor
    variable = u
    block = '1 2'
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

[Debug]
  show_var_residual_norms = true
[]

(assessment/TRISO/validation/AGR-34/SharedFiles/capsule_Ag.i)

[InterfaceKernels]
  [fuel_matrix]
    type = SorptionIsothermGapInterface
    variable = conc_Ag
    neighbor_var = conc_Ag
    partial_pressure_name = partial_pressure
    sorption_penalty = 1e5
    diffusivity = arrhenius_diffusion_coef_Ag
    use_flux_penalty = true
    flux_penalty = 1e3
    boundary = 'DRV_IR'
  []
  [matrix_graphite]
    type = SorptionIsothermGapInterface
    variable = conc_Ag
    neighbor_var = conc_Ag
    partial_pressure_name = partial_pressure
    sorption_penalty = 1e5
    diffusivity = arrhenius_diffusion_coef_Ag
    use_flux_penalty = true
    flux_penalty = 1e3
    boundary = 'IR_OR'
  []
  [graphite_sink]
    type = SorptionIsothermGapInterface
    variable = conc_Ag
    neighbor_var = conc_Ag
    partial_pressure_name = partial_pressure
    sorption_penalty = 1e5
    diffusivity = arrhenius_diffusion_coef_Ag
    use_flux_penalty = true
    flux_penalty = 1e3
    boundary = 'OR_SR'
  []
[]

[Variables]
  [conc_Ag]
    scaling = 1e12 #1e18
  []
[]

[AuxVariables]
  [Ag_diff_coef]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [mass_Ag_dt]
    type = TimeDerivative
    variable = conc_Ag
    extra_vector_tags = 'ref'
  []
  [mass_Ag]
    type = ArrheniusDiffusion
    variable = conc_Ag
    arrhenius_prpty_name = arrhenius_diffusion_coef_Ag
    extra_vector_tags = 'ref'
  []
  [conc_source_dtf]
    type = BodyForce
    variable = conc_Ag
    block = 'DTF'
    function = DTF_MassIsotope
    extra_vector_tags = 'ref'
  []
  [conc_source_drv]
    type = BodyForce
    variable = conc_Ag
    block = 'DRV'
    function = DRV_MassIsotope
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [Ag_diff_coef]
    type = MaterialRealAux
    variable = Ag_diff_coef
    property = arrhenius_diffusion_coef_Ag
    execute_on = timestep_end
  []
[]

[BCs]
  [freesurf_conc_Ag]
    type = DirichletBC
    variable = conc_Ag
    boundary = 'right'
    value = 0.0
    extra_vector_tags = 'ref'
  []
[]

[Materials]
  [matrix_conc_Ag]
    type = ArrheniusDiffusionCoef
    d1 = 1.6 # m^2/s
    q1 = 258000 # J/mol
    temperature = temperature
    arrhenius_prpty_name = arrhenius_diffusion_coef_Ag
  []
  [graphite_conc_Ag]
    type = ArrheniusDiffusionCoef
    d1 = 1.38e-2 # m^2/s
    q1 = 226000 # J/mol
    temperature = temperature
    arrhenius_prpty_name = arrhenius_diffusion_coef_Ag
  []
  [Ag_matrix_partial_pressure]
    type = SorptionPartialPressure
    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 = density # convert from mmol/m^3 to mmol/kg
    concentration = conc_Ag
    temperature = temperature
    output_properties = partial_pressure
  []
  [Ag_graphite_partial_pressure]
    type = SorptionPartialPressure
    A = 24
    B = -35700
    D = -1.56
    E = 6120
    d1 = 2.04
    d2 = 1.79e-3
    unit_scale = 1e3 # convert from mol to mmol
    density = density # convert from mmol/m^3 to mmol/kg
    concentration = conc_Ag
    temperature = temperature
    output_properties = partial_pressure
  []
[]

[Postprocessors]
  [release_Ag_inc_pebble]
    type = SideIntegralMassFlux
    variable = conc_Ag
    boundary = right
    arrhenius_prpty_name = arrhenius_diffusion_coef_Ag
    execute_on = 'initial timestep_end'
  []
  [released_Ag]
    type = TimeIntegratedPostprocessor
    value = release_Ag_inc_pebble
    execute_on = 'initial timestep_end'
  []
  [retained_Ag]
    type = ElementIntegralVariablePostprocessor
    variable = conc_Ag
  []

  # check partial pressure equality and integral flux conservation from DRV to IR
  [flux_DRV_outer]
    type = SideDiffusiveFluxIntegral
    boundary = 'DRV_IR'
    variable = conc_Ag
    diffusivity = arrhenius_diffusion_coef_Ag
  []
  [flux_IR_inner]
    type = SideDiffusiveFluxIntegral
    boundary = 'IR_DRV'
    variable = conc_Ag
    diffusivity = arrhenius_diffusion_coef_Ag
  []
  [flux_DRV_IR_error]
    type = FunctionValuePostprocessor
    function = flux_DRV_IR_error
  []

  # check partial pressure equality and integral flux conservation from IR to OR
  [flux_IR_outer]
    type = SideDiffusiveFluxIntegral
    boundary = 'IR_OR'
    variable = conc_Ag
    diffusivity = arrhenius_diffusion_coef_Ag
  []
  [flux_OR_inner]
    type = SideDiffusiveFluxIntegral
    boundary = 'OR_IR'
    variable = conc_Ag
    diffusivity = arrhenius_diffusion_coef_Ag
  []
  [flux_IR_OR_error]
    type = FunctionValuePostprocessor
    function = flux_IR_OR_error
  []

  # check partial pressure equality and integral flux conservation from OR to SR
  [flux_OR_outer]
    type = SideDiffusiveFluxIntegral
    boundary = 'OR_SR'
    variable = conc_Ag
    diffusivity = arrhenius_diffusion_coef_Ag
  []
  [flux_SR_inner]
    type = SideDiffusiveFluxIntegral
    boundary = 'SR_OR'
    variable = conc_Ag
    diffusivity = arrhenius_diffusion_coef_Ag
  []
  [flux_OR_SR_error]
    type = FunctionValuePostprocessor
    function = flux_OR_SR_error
  []
[]

[VectorPostprocessors]
  [Ag_conc]
    type = LineValueSampler
    variable = conc_Ag
    start_point = '0.0 0 0.0'
    num_points = 1000
    use_displaced_mesh = false
    execute_on = final
    sort_by = x
    outputs = line_plot_Ag
  []
[]

[Outputs]
  [line_plot_Ag]
    type = CSV
    execute_on = 'FINAL'
    time_step_interval =  1
    create_final_symlink = true
    file_base = capsule_radial
  []
[]

(assessment/TRISO/validation/AGR-34/SharedFiles/capsule_Cs.i)

[InterfaceKernels]
  [fuel_matrix]
    type = SorptionIsothermGapInterface
    variable = conc_Cs
    neighbor_var = conc_Cs
    partial_pressure_name = partial_pressure
    sorption_penalty = 1e5
    diffusivity = arrhenius_diffusion_coef_Cs
    use_flux_penalty = true
    flux_penalty = 1e3
    boundary = 'DRV_IR'
  []
  [matrix_graphite]
    type = SorptionIsothermGapInterface
    variable = conc_Cs
    neighbor_var = conc_Cs
    partial_pressure_name = partial_pressure
    sorption_penalty = 1e5
    diffusivity = arrhenius_diffusion_coef_Cs
    use_flux_penalty = true
    flux_penalty = 1e3
    boundary = 'IR_OR'
  []
  [graphite_sink]
    type = SorptionIsothermGapInterface
    variable = conc_Cs
    neighbor_var = conc_Cs
    partial_pressure_name = partial_pressure
    sorption_penalty = 1e5
    diffusivity = arrhenius_diffusion_coef_Cs
    use_flux_penalty = true
    flux_penalty = 1e3
    boundary = 'OR_SR'
  []
[]

[Variables]
  [conc_Cs]
    scaling = 1e12 #1e18
  []
[]

[AuxVariables]
  [Cs_diff_coef]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [mass_Cs_dt]
    type = TimeDerivative
    variable = conc_Cs
    extra_vector_tags = 'ref'
  []
  [mass_Cs]
    type = ArrheniusDiffusion
    variable = conc_Cs
    arrhenius_prpty_name = arrhenius_diffusion_coef_Cs
    extra_vector_tags = 'ref'
  []
  [conc_source_dtf]
    type = BodyForce
    variable = conc_Cs
    block = 'DTF'
    function = DTF_MassIsotope
    extra_vector_tags = 'ref'
  []
  [conc_source_drv]
    type = BodyForce
    variable = conc_Cs
    block = 'DRV'
    function = DRV_MassIsotope
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [Cs_diff_coef]
    type = MaterialRealAux
    variable = Cs_diff_coef
    property = arrhenius_diffusion_coef_Cs
    execute_on = timestep_end
  []
[]

[BCs]
  [freesurf_conc_Cs]
    type = DirichletBC
    variable = conc_Cs
    boundary = 'right'
    value = 0.0
    extra_vector_tags = 'ref'
  []
[]

[Materials]
  [matrix_conc_Cs]
    type = ArrheniusDiffusionCoef
    d1 = 3.6e-4 #3.6e-4m^2/s
    q1 = 189000 # J/mol
    temperature = temperature
    arrhenius_prpty_name = arrhenius_diffusion_coef_Cs
  []
  [graphite_conc_Cs]
    type = ArrheniusDiffusionCoef
    d1 = 1.7e-6 # 1.7e-6 m^2/s
    q1 = 149000 # J/mol
    temperature = temperature
    arrhenius_prpty_name = arrhenius_diffusion_coef_Cs
  []
  [Cs_matrix_partial_pressure]
    type = SorptionPartialPressure
    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 = density # convert from mmol/m^3 to mmol/kg
    concentration = conc_Cs
    temperature = temperature
    output_properties = partial_pressure
  []
  [Cs_graphite_partial_pressure]
    type = SorptionPartialPressure
    A = 24.0
    B = -35730.
    D = -1.561
    E = 6123.
    d1 = -2.035
    d2 = 1.79e-3
    unit_scale = 1e3 # convert from mol to mmol
    density = density # convert from mmol/m^3 to mmol/kg
    concentration = conc_Cs
    temperature = temperature
    output_properties = partial_pressure
  []
[]

[Postprocessors]
  [release_Cs_inc_pebble]
    type = SideIntegralMassFlux
    variable = conc_Cs
    boundary = right
    arrhenius_prpty_name = arrhenius_diffusion_coef_Cs
    execute_on = 'initial timestep_end'
  []
  [released_Cs]
    type = TimeIntegratedPostprocessor
    value = release_Cs_inc_pebble
    execute_on = 'initial timestep_end'
  []
  [retained_Cs]
    type = ElementIntegralVariablePostprocessor
    variable = conc_Cs
  []

  # check partial pressure equality and integral flux conservation from DRV to IR
  [flux_DRV_outer]
    type = SideDiffusiveFluxIntegral
    boundary = 'DRV_IR'
    variable = conc_Cs
    diffusivity = arrhenius_diffusion_coef_Cs
  []
  [flux_IR_inner]
    type = SideDiffusiveFluxIntegral
    boundary = 'IR_DRV'
    variable = conc_Cs
    diffusivity = arrhenius_diffusion_coef_Cs
  []
  [flux_DRV_IR_error]
    type = FunctionValuePostprocessor
    function = flux_DRV_IR_error
  []

  # check partial pressure equality and integral flux conservation from IR to OR
  [flux_IR_outer]
    type = SideDiffusiveFluxIntegral
    boundary = 'IR_OR'
    variable = conc_Cs
    diffusivity = arrhenius_diffusion_coef_Cs
  []
  [flux_OR_inner]
    type = SideDiffusiveFluxIntegral
    boundary = 'OR_IR'
    variable = conc_Cs
    diffusivity = arrhenius_diffusion_coef_Cs
  []
  [flux_IR_OR_error]
    type = FunctionValuePostprocessor
    function = flux_IR_OR_error
  []

  # check partial pressure equality and integral flux conservation from OR to SR
  [flux_OR_outer]
    type = SideDiffusiveFluxIntegral
    boundary = 'OR_SR'
    variable = conc_Cs
    diffusivity = arrhenius_diffusion_coef_Cs
  []
  [flux_SR_inner]
    type = SideDiffusiveFluxIntegral
    boundary = 'SR_OR'
    variable = conc_Cs
    diffusivity = arrhenius_diffusion_coef_Cs
  []
  [flux_OR_SR_error]
    type = FunctionValuePostprocessor
    function = flux_OR_SR_error
  []
[]

[VectorPostprocessors]
  [Cs_conc]
    type = LineValueSampler
    variable = conc_Cs
    start_point = '0.0 0 0.0'
    num_points = 1000
    use_displaced_mesh = false
    execute_on = final
    sort_by = x
    outputs = line_plot_Cs
  []
[]

[Outputs]
  [line_plot_Cs]
    type = CSV
    execute_on = 'FINAL'
    time_step_interval =  1
    create_final_symlink = true
    file_base = capsule_radial
  []
[]

(test/tests/sorption_interface/sorption_ad.i)

# This test is verifies the implementation of SorptionIsothermGapInterface with
# interface_configuration = solid_gas. It contains two solids separated by a gas.
# One solid acts as a source of u. u desorbs into the gap and then onto the
# other solid, which acts as a sink. Temperature increases over time to verify
# sorption behavior and species conservation. For SorptionIsothermGapInterface,
# the solid is always the element and the gas is always the neighbor.

#  _________________________
# | source | gap    | sink   |
#  --------------------------

[GlobalParams]
  order = SECOND
[]

[Mesh]
  coord_type = RZ
  [gen]
    type = GeneratedMeshGenerator
    nx = 150
    xmax = 1.5
    dim = 1
    elem_type = EDGE3
  []
  [source_gap]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '0.75 0 0'
    input = gen
  []
  [gap_sink]
    type = SubdomainBoundingBoxGenerator
    block_id = 2
    bottom_left = '0.75 0 0'
    top_right = '1.5 0 0'
    input = source_gap
  []
  [source]
    type = SubdomainBoundingBoxGenerator
    block_id = 3
    bottom_left = '0 0 0'
    top_right = '0.5 0 0'
    input = gap_sink
  []
  [sink]
    type = SubdomainBoundingBoxGenerator
    block_id = 4
    bottom_left = '1 0 0'
    top_right = '1.5 0 0'
    input = source
  []
  [break]
    type = BreakMeshByBlockGenerator
    block_pairs = '1 3; 2 4'
    split_interface = true
    add_interface_on_two_sides = true
    input = sink
  []
[]

[Variables]
  [u]
  []
[]

[AuxVariables]
  [temperature]
  []
[]

[Kernels]
  [u_dt]
    type = ADTimeDerivative
    variable = u
  []
  [u]
    type = ADMatDiffusion
    variable = u
    diffusivity = diffusivity
  []
[]

[ICs]
  [u_source]
    type = FunctionIC
    variable = u
    function = 'if(x < 0.75, 1, 0)'
    block = 3
  []
[]

[AuxKernels]
  [temperature]
    type = FunctionAux
    variable = temperature
    function = '(4000 - 2000 * x) / 2.0 * (t + 0.1)'
    execute_on = 'INITIAL TIMESTEP_BEGIN LINEAR NONLINEAR TIMESTEP_END'
  []
[]

[InterfaceKernels]
  [u_solid_gas]
    type = ADSorptionIsothermGapInterface
    variable = u
    neighbor_var = u
    interface_configuration = solid_gas
    temperature = temperature
    partial_pressure_name = partial_pressure
    sorption_penalty = 1e10
    boundary = 'Block3_Block1 Block4_Block2'
  []
[]

[Materials]
  [properties_solid]
    type = ADGenericConstantMaterial
    prop_names = 'density diffusivity'
    prop_values = '1 0.1'
    block = '3 4'
  []
  [properties_gas]
    type = ADGenericConstantMaterial
    prop_names = diffusivity
    prop_values = 1
    block = '1 2'
  []
  [partial_pressure_source]
    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 = 3
    outputs = 'all'
    output_properties = partial_pressure
  []
  [partial_pressure_sink]
    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 = 4
    outputs = 'all'
    output_properties = partial_pressure
  []
[]

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

[Executioner]
  type = Transient
  solve_type = NEWTON
  automatic_scaling = true
  compute_scaling_once = false

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  line_search = none

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

  dt = 0.1
  end_time = 2
[]

[Functions]
  [source_surface_u_exact]
    type = ParsedFunction
    symbol_names = 'P T R'
    symbol_values = 'source_surface_pressure source_surface_temperature 8.31446261815324'
    expression = '1.0 / R / T * P'
  []
  [sink_surface_u_exact]
    type = ParsedFunction
    symbol_names = 'P T R'
    symbol_values = 'sink_surface_pressure sink_surface_temperature 8.31446261815324'
    expression = '1.0 / R / T * P'
  []
[]

[Postprocessors]
  [source_surface_pressure]
    type = ADSideAverageMaterialProperty
    property = partial_pressure
    boundary = Block3_Block1
  []
  [source_surface_temperature]
    type = SideAverageValue
    variable = temperature
    boundary = Block3_Block1
  []
  [source_surface_u]
    type = SideAverageValue
    variable = u
    boundary = Block1_Block3
  []
  [source_surface_u_exact]
    type = FunctionValuePostprocessor
    function = source_surface_u_exact
  []

  [sink_surface_pressure]
    type = ADSideAverageMaterialProperty
    property = partial_pressure
    boundary = Block4_Block2
  []
  [sink_surface_temperature]
    type = SideAverageValue
    variable = temperature
    boundary = Block4_Block2
  []
  [sink_surface_u]
    type = SideAverageValue
    variable = u
    boundary = Block2_Block4
  []
  [sink_surface_u_exact]
    type = FunctionValuePostprocessor
    function = sink_surface_u_exact
  []

  [u_source_tot]
    type = ElementIntegralVariablePostprocessor
    variable = u
    block = 3
  []
  [u_gap_tot]
    type = ElementIntegralVariablePostprocessor
    variable = u
    block = '1 2'
  []
  [u_sink_tot]
    type = ElementIntegralVariablePostprocessor
    variable = u
    block = 4
  []
  [u_tot]
    type = ElementIntegralVariablePostprocessor
    variable = u
    block = '1 2 3 4'
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

[Debug]
  show_var_residual_norms = true
[]

Overview
Solid - Solid Interface Configuration
Solid - Gas Interface Configuration
Mass Flux Balance
Example Input Syntax
Input Parameters
Input Files
References