ADInterfaceSorption / InterfaceSorption

Computes a sorption law at interface between solid and gas in isothermal conditions.

Description

The surface concentrations $var element = document.getElementById("moose-equation-ab261ce1-4ce9-44e0-8c58-5039846bc036");katex.render("C_s", element, {displayMode:false,throwOnError:false});$ in the solid are related to the partial pressure $var element = document.getElementById("moose-equation-39269801-4be6-4a2c-8db5-b06b868e8735");katex.render("P", element, {displayMode:false,throwOnError:false});$ (and corresponding gas concentration $var element = document.getElementById("moose-equation-ac97e8d8-7432-455c-9f7f-99cf14f09013");katex.render("C_g", element, {displayMode:false,throwOnError:false});$ assuming ideal gas behavior) of the neighboring gas by the interface sorption law: $var element = document.getElementById("moose-equation-8702108d-4c92-4a51-8dbd-dfa04095af96");katex.render("C_s = K P^n = K \\left(C_g R T \\right)^n,", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-b1e9d15f-0062-4c21-8429-d2d5721d470c");katex.render("R", element, {displayMode:false,throwOnError:false});$ is the ideal gas constant in J/mol/K, $var element = document.getElementById("moose-equation-6bbd7232-af20-46c0-b423-a4f9ee9d5e05");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the temperature in K, $var element = document.getElementById("moose-equation-7f7c27a3-d91e-4db1-b94a-d4787f76d4ee");katex.render("K", element, {displayMode:false,throwOnError:false});$ is the solubility, and $var element = document.getElementById("moose-equation-9b9bd2a6-e728-4ac9-a712-cf7d3a684db9");katex.render("n", element, {displayMode:false,throwOnError:false});$ is the exponent of the sorption law. $var element = document.getElementById("moose-equation-0afd25a9-63e6-4ab1-b786-5ffb5c2b50ca");katex.render("K", element, {displayMode:false,throwOnError:false});$ is described as an Arrhenius equation: $var element = document.getElementById("moose-equation-7ba53f7c-2a2c-4862-8742-8def5a64d72d");katex.render("K = K_0 \\exp \\left(\\frac{-E_a}{RT}\\right),", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-8e40a8c3-c6e3-4de9-ae46-4d3bc4c7fc47");katex.render("K_0", element, {displayMode:false,throwOnError:false});$ is the pre-exponential constant, and $var element = document.getElementById("moose-equation-f194f679-b26b-4da3-b553-99be13876e63");katex.render("E_a", element, {displayMode:false,throwOnError:false});$ is the activation energy in J/mol. The units of $var element = document.getElementById("moose-equation-3f8e23f9-65b1-4c0c-98c5-0a60224710fe");katex.render("K", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-2ef3e8cf-9426-4a2c-b091-af3e8262aeac");katex.render("K_0", element, {displayMode:false,throwOnError:false});$ depend on the units of $var element = document.getElementById("moose-equation-b3108041-467e-4d08-9368-56b9521a0cb1");katex.render("C_s", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-d9361256-7ccb-4c98-a60d-3221bfb5ff65");katex.render("P", element, {displayMode:false,throwOnError:false});$ , as well as the value of $var element = document.getElementById("moose-equation-e1150312-b30e-40db-8892-a97a37dee3b9");katex.render("n", element, {displayMode:false,throwOnError:false});$ .

The value of $var element = document.getElementById("moose-equation-5626dc0a-d261-4597-9730-1f7f023c646a");katex.render("n", element, {displayMode:false,throwOnError:false});$ determines the sorption law. For example:

$var element = document.getElementById("moose-equation-c8dfdbb8-784c-4651-a4a2-4920185d47bc");katex.render("n = 1", element, {displayMode:false,throwOnError:false});$ corresponds to Henry's law
$var element = document.getElementById("moose-equation-7ec77688-85ac-42e3-ae5f-55330b9b337a");katex.render("n = 1/2", element, {displayMode:false,throwOnError:false});$ corresponds to Sievert's law.

A penalty is applied to enforce the sorption law at the interface. Note that solubilities in the literature may only be compatible with concentrations expressed in certain units, and those units may differ from those used elsewhere in the simulation, e.g., $var element = document.getElementById("moose-equation-9e48fcc8-257f-4d90-b3a0-77523122f3c6");katex.render("mmol/kg", element, {displayMode:false,throwOnError:false});$ versus $var element = document.getElementById("moose-equation-1697a239-4995-4249-a4e6-80e59a03b25c");katex.render("mol/m^3", element, {displayMode:false,throwOnError:false});$ . The model accepts concentration unit conversion factors ("unit_scale" and "unit_scale_neighbor") to accommodate data from different sources. For example, unit conversion factors of $var element = document.getElementById("moose-equation-0d38b76e-840c-49d1-9acc-cdd9370e55cd");katex.render("1000", element, {displayMode:false,throwOnError:false});$ can be supplied to accommodate sets of solubilities for concentrations in $var element = document.getElementById("moose-equation-7ed40023-676d-4992-94d0-b79d6855b5f6");katex.render("mmol/m^3", element, {displayMode:false,throwOnError:false});$ . Note that for multi-atomic molecules, it is important to specify if the moles correspond to the moles of gas molecules, or the moles of atoms. The model assumes that the units associated with the partial pressures are equal and that all constants are compatible with absolute temperature.

To balance the mass flux across the gap, a second interfacial condition is given as:

$(1)var element = document.getElementById("moose-equation-449fadce-bafb-48ac-978a-aebb83ba91f0");katex.render("D_s \\nabla C_s \\cdot \\mathbf{n_s} + D_g \\nabla C_g \\cdot \\mathbf{n_g} = 0, ", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-08fe46d7-6178-418d-a329-21784eaec320");katex.render("D_s", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-f48d526f-a088-4af3-b01d-ea8ef944a8dd");katex.render("D_g", element, {displayMode:false,throwOnError:false});$ are the diffusivities in the solid and the gas, respectively, and $var element = document.getElementById("moose-equation-931bff56-b081-4039-9edf-03bb4dc67720");katex.render("\\mathbf{n_s}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-169376be-f5da-47a3-bf21-f40b3b5f81cc");katex.render("\\mathbf{n_g}", element, {displayMode:false,throwOnError:false});$ are the normals at the interface. Two methods are available to enforce flux conditions at the interface. By default, InterfaceSorption 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.

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

When using the default, non-penalty method to enforce 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 InterfaceSorption 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.

Note that the unit conversion factors ("unit_scale" and "unit_scale_neighbor") do not apply to the flux equation to ensure mass conservation.

Input File Usage Example

[InterfaceKernels<<<{"href": "../../syntax/InterfaceKernels/index.html"}>>>]
  [interface]
    type = InterfaceSorption<<<{"description": "Computes a sorption law at interface between solid and gas in isothermal conditions.", "href": "InterfaceSorption.html"}>>>
    K0<<<{"description": "The pre-exponential factor for the Arrhenius law for the solubility $K$ for the relationship $C_i = KP_i^n$"}>>> = ${solubility}
    Ea<<<{"description": "The activation energy for the Arrhenius law for the solubility $K$ for the relationship $C_i = KP_i^n$"}>>> = 0
    n_sorption<<<{"description": "The exponent $n$ for the relationship $C_i = KP_i^n$"}>>> = ${n_sorption}
    diffusivity<<<{"description": "The diffusion coefficient"}>>> = diffusivity
    unit_scale<<<{"description": "Unit conversion factor used to scale the concentration"}>>> = ${unit_scale}
    unit_scale_neighbor<<<{"description": "Unit conversion factor used to scale the neighbor (gas) concentration"}>>> = ${unit_scale_neighbor}
    temperature<<<{"description": "The coupled temperature variable"}>>> = temperature
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u2
    neighbor_var<<<{"description": "The variable on the other side of the interface."}>>> = u1
    sorption_penalty<<<{"description": "Penalty associated with the concentration imbalance"}>>> = 1e1
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = interface2
  []
[]

Input Parameters

K0The pre-exponential factor for the Arrhenius law for the solubility $K$ for the relationship $C_i = KP_i^n$
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The pre-exponential factor for the Arrhenius law for the solubility $K$ for the relationship $C_i = KP_i^n$
n_sorptionThe exponent $n$ for the relationship $C_i = KP_i^n$
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The exponent $n$ for the relationship $C_i = KP_i^n$
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.
temperatureThe coupled temperature variable
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The coupled temperature variable
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

Ea0The activation energy for the Arrhenius law for the solubility $K$ for the relationship $C_i = KP_i^n$
Default:0
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The activation energy for the Arrhenius law for the solubility $K$ for the relationship $C_i = KP_i^n$
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
matrix_onlyFalseWhether this object is only doing assembly to matrices (no vectors)
Default:False
C++ Type:bool
Controllable:No
Description:Whether this object is only doing assembly to matrices (no vectors)
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
unit_scale1Unit conversion factor used to scale the concentration
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Unit conversion factor used to scale the concentration
unit_scale_neighbor1Unit conversion factor used to scale the neighbor (gas) concentration
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Unit conversion factor used to scale the neighbor (gas) concentration
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/ver-1kc-2/ver-1kc-2.i)
(test/tests/val-2c/val-2c_base.i)
(test/tests/ver-1kb/ver-1kb.i)
(test/tests/interfacekernels/InterfaceSorption/interface_sorption.i)
(test/tests/interfacekernels/InterfaceSorption/interface_sorption_transient.i)
(test/tests/ver-1kd/ver-1kd.i)
(test/tests/ver-1kc-1/ver-1kc-1.i)

(test/tests/interfacekernels/InterfaceSorption/interface_sorption.i)

# This test is to verify the implementation of InterfaceSorption and its AD counterpart.
# It contains two 1D blocks separated by a continuous interface.
# InterfaceSorption is used to enforce the sorption law and preserve flux between the blocks.
# Checks are performed to verify concentration conservation, sorption behavior, and flux preservation.
# This input file uses BreakMeshByBlockGenerator, which is currently only supported for replicated
# meshes, so this file should not be run with the `parallel_type = DISTRIBUTED` flag

# In this input file, we apply the Sievert law with n_sorption=1/2.

# Physical Constants
R = 8.31446261815324 # J/mol/K, based on number used in include/utils/PhysicalConstants.h

solubility = 1.e-2
n_sorption = 0.5

unit_scale = 1
unit_scale_neighbor = 1


[GlobalParams]
  order = FIRST
  family = LAGRANGE
[]

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    nx = 10
    dim = 1
  []
  [block1]
    type = SubdomainBoundingBoxGenerator
    input = gen
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '0.5 1 0'
  []
  [block2]
    type = SubdomainBoundingBoxGenerator
    input = block1
    block_id = 2
    bottom_left = '0.5 0 0'
    top_right = '1 1 0'
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = block2
    primary_block = 1
    paired_block = 2
    new_boundary = interface
  []
  [interface2]
    type = SideSetsBetweenSubdomainsGenerator
    input = interface
    primary_block = 2
    paired_block = 1
    new_boundary = interface2
  []
[]

[Variables]
  [u1]
    block = 1
  []
  [u2]
    block = 2
  []
  [temperature]
  []
[]

[Kernels]
  [u1]
    type = MatDiffusion
    variable = u1
    diffusivity = diffusivity
    block = 1
  []
  [u2]
    type = MatDiffusion
    variable = u2
    diffusivity = diffusivity
    block = 2
  []
  [temperature]
    type = HeatConduction
    variable = temperature
  []
[]

[BCs]
  [left_u1]
    type = DirichletBC
    value = 1
    variable = u1
    boundary = left
  []
  [right_u2]
    type = DirichletBC
    value = 1
    variable = u2
    boundary = right
  []
  [left_temperature]
    type = DirichletBC
    value = 1100
    variable = temperature
    boundary = left
  []
  [right_temperature]
    type = DirichletBC
    value = 0
    variable = temperature
    boundary = right
  []
  [block1_2_temperature]
    type = DirichletBC
    value = 1000
    variable = temperature
    boundary = interface
  []
  [block2_1_temperature]
    type = DirichletBC
    value = 900
    variable = temperature
    boundary = interface2
  []
[]

[InterfaceKernels]
  [interface]
    type = InterfaceSorption
    K0 = ${solubility}
    Ea = 0
    n_sorption = ${n_sorption}
    diffusivity = diffusivity
    unit_scale = ${unit_scale}
    unit_scale_neighbor = ${unit_scale_neighbor}
    temperature = temperature
    variable = u2
    neighbor_var = u1
    sorption_penalty = 1e1
    boundary = interface2
  []
[]

[Materials]
  [properties_1]
    type = GenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity'
    prop_values = '1 1'
    block = 1
  []
  [properties_2]
    type = GenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity solubility'
    prop_values = '2 2 ${solubility}'
    block = 2
  []
[]

[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)'
  []
  [residual_concentration]
    type = ParsedFunction
    symbol_names = 'u_mid_inner u_mid_outer T'
    symbol_values = 'u_mid_inner u_mid_outer temperature_mid_inner'
    expression = 'u_mid_outer*${unit_scale} - ${solubility}*(u_mid_inner*${unit_scale_neighbor}*${R}*T)^${n_sorption}'
  []
  [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

  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
[]

[Postprocessors]
  [u_mid_inner]
    type = PointValue
    variable = u1
    point = '0.49999 0 0'
    outputs = 'csv console'
  []
  [u_mid_outer]
    type = PointValue
    variable = u2
    point = '0.50001 0 0'
    outputs = 'csv console'
  []
  [u_mid_diff]
    type = FunctionValuePostprocessor
    function = u_mid_diff
    outputs = 'csv console'
  []
  [temperature_mid_inner]
    type = PointValue
    variable = temperature
    point = '0.49999 0 0'
    outputs = csv
  []
  [temperature_mid_outer]
    type = PointValue
    variable = temperature
    point = '0.50001 0 0'
    outputs = csv
  []
  [residual_concentration]
    type = FunctionValuePostprocessor
    function = residual_concentration
    outputs = 'csv console'
  []

  [flux_inner] # verify flux preservation
    type = SideDiffusiveFluxIntegral
    variable = u1
    boundary = interface
    diffusivity = diffusivity
    outputs = 'csv console'
  []
  [flux_outer]
    type = SideDiffusiveFluxIntegral
    variable = u2
    boundary = interface2
    diffusivity = diffusivity
    outputs = 'csv console'
  []
  [flux_error]
    type = FunctionValuePostprocessor
    function = flux_error
    outputs = 'csv console'
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

[Debug]
  show_var_residual_norms = true
[]

(test/tests/ver-1kc-2/ver-1kc-2.i)

nb_segments_TMAP7 = 20
node_size_TMAP7 = '${units 1.25e-5 m}'
long_total = '${units ${fparse nb_segments_TMAP7 * node_size_TMAP7} m}'
nb_segments_TMAP8 = 1e2
simulation_time = '${units 0.25 s}'
temperature = '${units 500 K}'
R = '${units 8.31446261815324 J/mol/K}' # ideal gas constant from PhysicalConstants.h
initial_pressure_1 = '${units 1e5 Pa}'
initial_pressure_2 = '${units 1e-10 Pa}'
initial_concentration_1 = '${units ${fparse initial_pressure_1 / (R*temperature)} mol/m^3}'
initial_concentration_2 = '${units ${fparse initial_pressure_2 / (R*temperature)} mol/m^3}'
solubility = '${units ${fparse 10/sqrt(R*temperature)} mol/m^3/Pa^(1/2)}' # Sieverts' law solubility
diffusivity = '${units ${fparse 4.31e-6 * exp(-2818/temperature)} m^2/s}'
n_sorption = 0.5 # Sieverts' Law
K1 = '${units 4000 mol/m^3/s}' # reaction rate for H2+T2->2HT
equilibrium_constant = 2.0
K2 = '${units ${fparse (2*K1) / (equilibrium_constant)^2} mol/m^3/s}' # reaction rate for 2HT->H2+T2
unit_scale = 1
unit_scale_neighbor = 1

[Mesh]
  [generated]
    type = GeneratedMeshGenerator
    dim = 1
    nx = ${nb_segments_TMAP8}
    xmax = ${long_total}
  []
  [enclosure_1]
    type = SubdomainBoundingBoxGenerator
    input = generated
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '${fparse 1/3 * long_total} 0 0'
  []
  [enclosure_2]
    type = SubdomainBoundingBoxGenerator
    input = enclosure_1
    block_id = 2
    bottom_left = '${fparse 1/3 * long_total} 0 0'
    top_right = '${fparse long_total} 0 0'
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = enclosure_2
    primary_block = 1
    paired_block = 2
    new_boundary = interface
  []
  [interface2]
    type = SideSetsBetweenSubdomainsGenerator
    input = interface
    primary_block = 2
    paired_block = 1
    new_boundary = interface2
  []
[]

[Variables]

  # Variables for H2, T2, and HT in enclosure 1
  [concentration_H2_enclosure_1]
    block = 1
    initial_condition = '${initial_concentration_1}'
  []
  [concentration_T2_enclosure_1]
    block = 1
    initial_condition = '${initial_concentration_1}'
  []
  [concentration_HT_enclosure_1]
    block = 1
    initial_condition = '${initial_concentration_2}'
  []

  # Variables for H2, T2, and HT in enclosure 2
  [concentration_H2_enclosure_2]
    block = 2
    initial_condition = '${initial_concentration_2}'
  []
  [concentration_T2_enclosure_2]
    block = 2
    initial_condition = '${initial_concentration_2}'
  []
  [concentration_HT_enclosure_2]
    block = 2
    initial_condition = '${initial_concentration_2}'
  []
[]

[AuxVariables]

  # Auxiliary variables for H2, T2, and HT in enclosure 1
  [pressure_H2_enclosure_1]
    order = CONSTANT
    family = MONOMIAL
    block = 1
    initial_condition = '${initial_pressure_1}'
  []
  [pressure_T2_enclosure_1]
    order = CONSTANT
    family = MONOMIAL
    block = 1
    initial_condition = '${initial_pressure_1}'
  []
  [pressure_HT_enclosure_1]
    order = CONSTANT
    family = MONOMIAL
    block = 1
    initial_condition = '${initial_pressure_2}'
  []

  # Auxiliary variables for H2, T2, and HT in enclosure 2
  [pressure_H2_enclosure_2]
    order = CONSTANT
    family = MONOMIAL
    block = 2
    initial_condition = '${initial_pressure_2}'
  []
  [pressure_T2_enclosure_2]
    order = CONSTANT
    family = MONOMIAL
    block = 2
    initial_condition = '${initial_pressure_2}'
  []
  [pressure_HT_enclosure_2]
    order = CONSTANT
    family = MONOMIAL
    block = 2
    initial_condition = '${initial_pressure_2}'
  []
[]

[Kernels]
  # Diffusion equation for H2
  [H2_diffusion_enclosure_1]
    type = MatDiffusion
    variable = concentration_H2_enclosure_1
    diffusivity = ${diffusivity}
    block = '1'
  []
  [H2_time_derivative_enclosure_1]
    type = TimeDerivative
    variable = concentration_H2_enclosure_1
    block = '1'
  []
  [H2_diffusion_enclosure_2]
    type = MatDiffusion
    variable = concentration_H2_enclosure_2
    diffusivity = ${diffusivity}
    block = '2'
  []
  [H2_time_derivative_enclosure_2]
    type = TimeDerivative
    variable = concentration_H2_enclosure_2
    block = '2'
  []

  # Diffusion equation for T2
  [T2_diffusion_enclosure_1]
    type = MatDiffusion
    variable = concentration_T2_enclosure_1
    diffusivity = ${diffusivity}
    block = '1'
  []
  [T2_time_derivative_enclosure_1]
    type = TimeDerivative
    variable = concentration_T2_enclosure_1
    block = '1'
  []
  [T2_diffusion_enclosure_2]
    type = MatDiffusion
    variable = concentration_T2_enclosure_2
    diffusivity = ${diffusivity}
    block = '2'
  []
  [T2_time_derivative_enclosure_2]
    type = TimeDerivative
    variable = concentration_T2_enclosure_2
    block = '2'
  []

  # Diffusion equation for HT
  [HT_diffusion_enclosure_1]
    type = MatDiffusion
    variable = concentration_HT_enclosure_1
    diffusivity = ${diffusivity}
    block = '1'
  []
  [HT_time_derivative_enclosure_1]
    type = TimeDerivative
    variable = concentration_HT_enclosure_1
    block = '1'
  []
  [HT_diffusion_enclosure_2]
    type = MatDiffusion
    variable = concentration_HT_enclosure_2
    diffusivity = ${diffusivity}
    block = '2'
  []
  [HT_time_derivative_enclosure_2]
    type = TimeDerivative
    variable = concentration_HT_enclosure_2
    block = '2'
  []

  # Reaction H2+T2->2HT in enclosure 1
  [reaction_H2_encl_1_1]
    type = ADMatReactionFlexible
    variable = concentration_H2_enclosure_1
    vs = 'concentration_H2_enclosure_1 concentration_T2_enclosure_1'
    block = 1
    coeff = -1
    reaction_rate_name = ${K1}
  []
  [reaction_T2_encl_1_1]
    type = ADMatReactionFlexible
    variable = concentration_T2_enclosure_1
    vs = 'concentration_H2_enclosure_1 concentration_T2_enclosure_1'
    block = 1
    coeff = -1
    reaction_rate_name = ${K1}
  []
  [reaction_HT_encl_1_1]
    type = ADMatReactionFlexible
    variable = concentration_HT_enclosure_1
    vs = 'concentration_H2_enclosure_1 concentration_T2_enclosure_1'
    block = 1
    coeff = 2
    reaction_rate_name = ${K1}
  []

  # Reaction 2HT->H2+T2 in enclosure 1
  [reaction_H2_encl_1_2]
    type = ADMatReactionFlexible
    variable = concentration_H2_enclosure_1
    vs = 'concentration_HT_enclosure_1 concentration_HT_enclosure_1'
    block = 1
    coeff = 0.5
    reaction_rate_name = ${K2}
  []
  [reaction_T2_encl_1_2]
    type = ADMatReactionFlexible
    variable = concentration_T2_enclosure_1
    vs = 'concentration_HT_enclosure_1 concentration_HT_enclosure_1'
    block = 1
    coeff = 0.5
    reaction_rate_name = ${K2}
  []
  [reaction_HT_encl_1_2]
    type = ADMatReactionFlexible
    variable = concentration_HT_enclosure_1
    vs = 'concentration_HT_enclosure_1 concentration_HT_enclosure_1'
    block = 1
    coeff = -1
    reaction_rate_name = ${K2}
  []

  # Reaction H2+T2->2HT in enclosure 2
  [reaction_H2_encl_2_1]
    type = ADMatReactionFlexible
    variable = concentration_H2_enclosure_2
    vs = 'concentration_H2_enclosure_2 concentration_T2_enclosure_2'
    block = 2
    coeff = -1
    reaction_rate_name = ${K1}
  []
  [reaction_T2_encl_2_1]
    type = ADMatReactionFlexible
    variable = concentration_T2_enclosure_2
    vs = 'concentration_H2_enclosure_2 concentration_T2_enclosure_2'
    block = 2
    coeff = -1
    reaction_rate_name = ${K1}
  []
  [reaction_HT_encl_2_1]
    type = ADMatReactionFlexible
    variable = concentration_HT_enclosure_2
    vs = 'concentration_H2_enclosure_2 concentration_T2_enclosure_2'
    block = 2
    coeff = 2
    reaction_rate_name = ${K1}
  []

  # Reaction 2HT->H2+T2 in enclosure 2
  [reaction_H2_encl_2_2]
    type = ADMatReactionFlexible
    variable = concentration_H2_enclosure_2
    vs = 'concentration_HT_enclosure_2 concentration_HT_enclosure_2'
    block = 2
    coeff = 0.5
    reaction_rate_name = ${K2}
  []
  [reaction_T2_encl_2_2]
    type = ADMatReactionFlexible
    variable = concentration_T2_enclosure_2
    vs = 'concentration_HT_enclosure_2 concentration_HT_enclosure_2'
    block = 2
    coeff = 0.5
    reaction_rate_name = ${K2}
  []
  [reaction_HT_encl_2_2]
    type = ADMatReactionFlexible
    variable = concentration_HT_enclosure_2
    vs = 'concentration_HT_enclosure_2 concentration_HT_enclosure_2'
    block = 2
    coeff = -1
    reaction_rate_name = ${K2}
  []
[]

[AuxKernels]

  # Auxiliary kernels for H2, T2, and HT in enclosure 1
  [pressure_H2_enclosure_1]
    type = ParsedAux
    variable = pressure_H2_enclosure_1
    coupled_variables = 'concentration_H2_enclosure_1'
    expression = '${fparse R*temperature}*concentration_H2_enclosure_1'
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_T2_enclosure_1]
    type = ParsedAux
    variable = pressure_T2_enclosure_1
    coupled_variables = 'concentration_T2_enclosure_1'
    expression = '${fparse R*temperature}*concentration_T2_enclosure_1'
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_HT_enclosure_1]
    type = ParsedAux
    variable = pressure_HT_enclosure_1
    coupled_variables = 'concentration_HT_enclosure_1'
    expression = '${fparse R*temperature}*concentration_HT_enclosure_1'
    block = 1
    execute_on = 'initial timestep_end'
  []

  # Auxiliary kernels for H2, T2, and HT in enclosure 2
  [pressure_H2_enclosure_2]
    type = ParsedAux
    variable = pressure_H2_enclosure_2
    coupled_variables = 'concentration_H2_enclosure_2'
    expression = '${fparse R*temperature}*concentration_H2_enclosure_2'
    block = 2
    execute_on = 'initial timestep_end'
  []
  [pressure_T2_enclosure_2]
    type = ParsedAux
    variable = pressure_T2_enclosure_2
    coupled_variables = 'concentration_T2_enclosure_2'
    expression = '${fparse R*temperature}*concentration_T2_enclosure_2'
    block = 2
    execute_on = 'initial timestep_end'
  []
  [pressure_HT_enclosure_2]
    type = ParsedAux
    variable = pressure_HT_enclosure_2
    coupled_variables = 'concentration_HT_enclosure_2'
    expression = '${fparse R*temperature}*concentration_HT_enclosure_2'
    block = 2
    execute_on = 'initial timestep_end'
  []
[]

[InterfaceKernels]
  [interface_sorption_H2]
    type = InterfaceSorption
    K0 = ${solubility}
    Ea = 0
    n_sorption = ${n_sorption}
    diffusivity = ${diffusivity}
    unit_scale = ${unit_scale}
    unit_scale_neighbor = ${unit_scale_neighbor}
    temperature = ${temperature}
    variable = concentration_H2_enclosure_1
    neighbor_var = concentration_H2_enclosure_2
    sorption_penalty = 1e1
    boundary = interface
  []
  [interface_sorption_T2]
    type = InterfaceSorption
    K0 = ${solubility}
    Ea = 0
    n_sorption = ${n_sorption}
    diffusivity = ${diffusivity}
    unit_scale = ${unit_scale}
    unit_scale_neighbor = ${unit_scale_neighbor}
    temperature = ${temperature}
    variable = concentration_T2_enclosure_1
    neighbor_var = concentration_T2_enclosure_2
    sorption_penalty = 1e1
    boundary = interface
  []
  [interface_sorption_HT]
    type = InterfaceSorption
    K0 = ${solubility}
    Ea = 0
    n_sorption = ${n_sorption}
    diffusivity = ${diffusivity}
    unit_scale = ${unit_scale}
    unit_scale_neighbor = ${unit_scale_neighbor}
    temperature = ${temperature}
    variable = concentration_HT_enclosure_1
    neighbor_var = concentration_HT_enclosure_2
    sorption_penalty = 1e1
    boundary = interface
  []
[]

[Postprocessors]

  # postprocessors for H2
  [pressure_H2_enclosure_1]
    type = ElementAverageValue
    variable = pressure_H2_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_H2_enclosure_2]
    type = ElementAverageValue
    variable = pressure_H2_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []
  [concentration_H2_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = concentration_H2_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_H2_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface2
    variable = concentration_H2_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_ratio_H2]
    type = ParsedPostprocessor
    expression = 'concentration_H2_enclosure_1_at_interface / sqrt(concentration_H2_enclosure_2_at_interface)'
    pp_names = 'concentration_H2_enclosure_1_at_interface concentration_H2_enclosure_2_at_interface'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
  [pressure_H2_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = pressure_H2_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [pressure_H2_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface2
    variable = pressure_H2_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_H2_encl_1_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_H2_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [concentration_H2_encl_2_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_H2_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []

  # postprocessors for T2
  [pressure_T2_enclosure_1]
    type = ElementAverageValue
    variable = pressure_T2_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_T2_enclosure_2]
    type = ElementAverageValue
    variable = pressure_T2_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []
  [concentration_T2_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = concentration_T2_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_T2_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface2
    variable = concentration_T2_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_ratio_T2]
    type = ParsedPostprocessor
    expression = 'concentration_T2_enclosure_1_at_interface / sqrt(concentration_T2_enclosure_2_at_interface)'
    pp_names = 'concentration_T2_enclosure_1_at_interface concentration_T2_enclosure_2_at_interface'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
  [concentration_T2_encl_1_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_T2_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [concentration_T2_encl_2_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_T2_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []

  # postprocessors for HT
  [pressure_HT_enclosure_1]
    type = ElementAverageValue
    variable = pressure_HT_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_HT_enclosure_2]
    type = ElementAverageValue
    variable = pressure_HT_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []
  [concentration_HT_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = concentration_HT_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_HT_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface2
    variable = concentration_HT_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_ratio_HT]
    type = ParsedPostprocessor
    expression = 'concentration_HT_enclosure_1_at_interface / sqrt(concentration_HT_enclosure_2_at_interface)'
    pp_names = 'concentration_HT_enclosure_1_at_interface concentration_HT_enclosure_2_at_interface'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
  [concentration_HT_encl_1_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_HT_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [concentration_HT_encl_2_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_HT_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []

  # postprocessors for mass conservation
  [mass_conservation_sum_encl1_encl2]
    type = LinearCombinationPostprocessor
    pp_names = 'concentration_HT_encl_1_inventory concentration_HT_encl_2_inventory concentration_H2_encl_1_inventory concentration_H2_encl_2_inventory concentration_T2_encl_1_inventory concentration_T2_encl_2_inventory'
    pp_coefs = '1            1           1            1            1            1'
    execute_on = 'initial timestep_end'
  []

  # postprocessors for equilibrium constant
  [equilibrium_constant_encl_1]
    type = ParsedPostprocessor
    expression = 'pressure_HT_enclosure_1 / sqrt(pressure_H2_enclosure_1 * pressure_T2_enclosure_1)'
    pp_names = 'pressure_HT_enclosure_1 pressure_H2_enclosure_1 pressure_T2_enclosure_1'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
  [equilibrium_constant_encl_2]
    type = ParsedPostprocessor
    expression = 'pressure_HT_enclosure_2 / sqrt(pressure_H2_enclosure_2 * pressure_T2_enclosure_2)'
    pp_names = 'pressure_HT_enclosure_2 pressure_H2_enclosure_2 pressure_T2_enclosure_2'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
[]

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

[Executioner]
  type = Transient
  end_time = ${simulation_time}
  dtmax = 1e-2
  nl_max_its = 9
  l_max_its = 30
  scheme = 'bdf2'
  solve_type = PJFNK
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e-6
    optimal_iterations = 7
    iteration_window = 1
    growth_factor = 1.1
    cutback_factor = 0.9
    cutback_factor_at_failure = 0.9
  []
[]

[Outputs]
  file_base = 'ver-1kc-2_out_k10'
  csv = true
  exodus = true
  execute_on = 'initial timestep_end'
[]

(test/tests/val-2c/val-2c_base.i)

# Validation problem #2c from TMAP4/TMAP7 V&V document
# Test Cell Release Experiment based on
# D. F. Holland and R. A. Jalbert, "A Model for Tritium Concentration Following Tritium
# Release into a Test Cell and Subsequent Operation of an Atmospheric Cleanup Systen,"
# Proceedings, Eleventh Symposium of Fusion Engineering, November 18-22, 1985,. Austin,
# TX, Vol I, pp. 638-43, IEEE Cat. No. CH2251-7.

# Note that the approach to model this validation case is different in TMAP4 and TMAP7.

[Mesh]
  [base_mesh]
    type = GeneratedMeshGenerator
    dim = 1
    nx = ${fparse mesh_num_nodes_paint + 2}
    xmax = ${length_domain}
  []
  [subdomain_id]
    input = base_mesh
    type = SubdomainPerElementGenerator
    subdomain_ids = '0 0 0 0 0 0 0 0 0 0 0 0 1 1'
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = subdomain_id
    primary_block = '0' # paint
    paired_block = '1' # enclosure
    new_boundary = 'interface'
  []
  [interface_other_side]
    type = SideSetsBetweenSubdomainsGenerator
    input = interface
    primary_block = '1' # enclosure
    paired_block = '0' # paint
    new_boundary = 'interface_other'
  []
[]

[Variables]
  # T2 concentration in the enclosure in molecules/microns^3
  [t2_enclosure_concentration]
    block = 1
  []
  # HT concentration in the enclosure in molecules/microns^3
  [ht_enclosure_concentration]
    block = 1
  []
  # HTO concentration in the enclosure in molecules/microns^3
  [hto_enclosure_concentration]
    block = 1
  []
  # H2O concentration in the enclosure in molecules/microns^3
  [h2o_enclosure_concentration]
    block = 1
    initial_condition = ${initial_H2O_concentration}
  []
  # concentration of T2 in the paint in molecules/microns^3
  [t2_paint_concentration]
    block = 0
  []
  # concentration of HT in the paint in molecules/microns^3
  [ht_paint_concentration]
    block = 0
  []
  # concentration of HTO in the paint in molecules/microns^3
  [hto_paint_concentration]
    block = 0
  []
  # concentration of H2O in the paint in molecules/microns^3
  [h2o_paint_concentration]
    block = 0
  []
[]

[AuxVariables]
  # Used to prevent negative concentrations
  [bounds_dummy_t2_paint_concentration]
    order = FIRST
    family = LAGRANGE
  []
  [bounds_dummy_t2_enclosure_concentration]
    order = FIRST
    family = LAGRANGE
  []
  [bounds_dummy_hto_paint_concentration]
    order = FIRST
    family = LAGRANGE
  []
  [bounds_dummy_hto_enclosure_concentration]
    order = FIRST
    family = LAGRANGE
  []
  [bounds_dummy_h2o_paint_concentration]
    order = FIRST
    family = LAGRANGE
  []
  [bounds_dummy_h2o_enclosure_concentration]
    order = FIRST
    family = LAGRANGE
  []
[]

[Bounds]
  # To prevent negative concentrations
  [t2_paint_concentration_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_t2_paint_concentration
    bounded_variable = t2_paint_concentration
    bound_type = lower
    bound_value = ${lower_value_threshold}
  []
  [t2_enclosure_concentration_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_t2_enclosure_concentration
    bounded_variable = t2_enclosure_concentration
    bound_type = lower
    bound_value = ${lower_value_threshold}
  []
  [hto_paint_concentration_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_hto_paint_concentration
    bounded_variable = hto_paint_concentration
    bound_type = lower
    bound_value = ${lower_value_threshold}
  []
  [hto_enclosure_concentration_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_hto_enclosure_concentration
    bounded_variable = hto_enclosure_concentration
    bound_type = lower
    bound_value = ${lower_value_threshold}
  []
  [h2o_paint_concentration_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_h2o_paint_concentration
    bounded_variable = h2o_paint_concentration
    bound_type = lower
    bound_value = -1e-4
  []
  [h2o_enclosure_concentration_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_h2o_enclosure_concentration
    bounded_variable = h2o_enclosure_concentration
    bound_type = lower
    bound_value = -1e-4
  []
[]

[Kernels]
  # In the enclosure
  [t2_time_derivative]
    type = TimeDerivative
    variable = t2_enclosure_concentration
    block = 1
  []
  [t2_outflow]
    type = MaskedBodyForce
    variable = t2_enclosure_concentration
    value = '-1'
    mask = t2_enclosure_concentration_outflow
    block = 1
  []
  [t2_diffusion]
    type = MatDiffusion
    variable = t2_enclosure_concentration
    block = 1
    diffusivity = ${diffusivity_artificial_enclosure}
  []
  [ht_time_derivative]
    type = TimeDerivative
    variable = ht_enclosure_concentration
    block = 1
  []
  [ht_outflow]
    type = MaskedBodyForce
    variable = ht_enclosure_concentration
    value = '-1'
    mask = ht_enclosure_concentration_outflow
    block = 1
  []
  [ht_diffusion]
    type = MatDiffusion
    variable = ht_enclosure_concentration
    block = 1
    diffusivity = ${diffusivity_artificial_enclosure}
  []
  [hto_time_derivative]
    type = TimeDerivative
    variable = hto_enclosure_concentration
    block = 1
  []
  [hto_outflow]
    type = MaskedBodyForce
    variable = hto_enclosure_concentration
    value = '-1'
    mask = hto_enclosure_concentration_outflow
    block = 1
  []
  [hto_diffusion]
    type = MatDiffusion
    variable = hto_enclosure_concentration
    block = 1
    diffusivity = ${diffusivity_artificial_enclosure}
  []
  [h2o_time_derivative]
    type = TimeDerivative
    variable = h2o_enclosure_concentration
    block = 1
  []
  [h2o_inflow]
    type = MaskedBodyForce
    variable = h2o_enclosure_concentration
    mask = ${inflow_concentration}
    block = 1
  []
  [h2o_outflow]
    type = MaskedBodyForce
    variable = h2o_enclosure_concentration
    value = '-1'
    mask = h2o_enclosure_concentration_outflow
    block = 1
  []
  [h2o_diffusion]
    type = MatDiffusion
    variable = h2o_enclosure_concentration
    block = 1
    diffusivity = ${diffusivity_artificial_enclosure}
  []

  # reaction T2+H2O->HTO+HT
  [reaction_1_t2]
    type = ADMatReactionFlexible
    variable = t2_enclosure_concentration
    block = 1
    coeff = -1
    reaction_rate_name = reaction_rate_t2
  []
  [reaction_1_h2o]
    type = ADMatReactionFlexible
    variable = h2o_enclosure_concentration
    block = 1
    coeff = -1
    reaction_rate_name = reaction_rate_t2
  []
  [reaction_1_hto]
    type = ADMatReactionFlexible
    variable = hto_enclosure_concentration
    block = 1
    coeff = 1
    reaction_rate_name = reaction_rate_t2
  []
  [reaction_1_ht]
    type = ADMatReactionFlexible
    variable = ht_enclosure_concentration
    block = 1
    coeff = 1
    reaction_rate_name = reaction_rate_t2
  []
  # reaction HT+H2O->HTO+H2
  [reaction_2_HT]
    type = ADMatReactionFlexible
    variable = ht_enclosure_concentration
    block = 1
    coeff = -1
    reaction_rate_name = reaction_rate_ht
  []
  [reaction_2_h2o]
    type = ADMatReactionFlexible
    variable = h2o_enclosure_concentration
    block = 1
    coeff = -1
    reaction_rate_name = reaction_rate_ht
  []
  [reaction_2_hto]
    type = ADMatReactionFlexible
    variable = hto_enclosure_concentration
    block = 1
    coeff = 1
    reaction_rate_name = reaction_rate_ht
  []

  # In the paint
  [t2_paint_time]
    type = TimeDerivative
    variable = t2_paint_concentration
    block = 0
  []
  [t2_paint_diffusion]
    type = MatDiffusion
    variable = t2_paint_concentration
    block = 0
    diffusivity = '${diffusivity_elemental_tritium}'
  []
  [ht_paint_time]
    type = TimeDerivative
    variable = ht_paint_concentration
    block = 0
  []
  [ht_paint_diffusion]
    type = MatDiffusion
    variable = ht_paint_concentration
    block = 0
    diffusivity = '${diffusivity_elemental_tritium}'
  []
  [hto_paint_time]
    type = TimeDerivative
    variable = hto_paint_concentration
    block = 0
  []
  [hto_paint_diffusion]
    type = MatDiffusion
    variable = hto_paint_concentration
    block = 0
    diffusivity = '${diffusivity_tritiated_water}'
  []
  [h2o_paint_time]
    type = TimeDerivative
    variable = h2o_paint_concentration
    block = 0
  []
  [h2o_paint_diffusion]
    type = MatDiffusion
    variable = h2o_paint_concentration
    block = 0
    diffusivity = '${diffusivity_tritiated_water}'
  []
[]

[InterfaceKernels]
  # solubility at the surface of the paint
  [t2_solubility]
    type = ADInterfaceSorption
    variable = t2_paint_concentration
    neighbor_var = t2_enclosure_concentration # molecules/microns^3
    unit_scale_neighbor = ${fparse 1e18/NA} # to convert neighbor concentration to mol/m^3
    K0 = ${solubility_elemental_tritium}
    Ea = 0
    n_sorption = 1 # Henry's law
    temperature = ${temperature}
    diffusivity = diffusivity_elemental_tritium
    boundary = 'interface'
  []
  [ht_solubility]
    type = ADInterfaceSorption
    variable = ht_paint_concentration
    neighbor_var = ht_enclosure_concentration # molecules/microns^3
    unit_scale_neighbor = ${fparse 1e18/NA} # to convert neighbor concentration to mol/m^3
    K0 = ${solubility_elemental_tritium}
    Ea = 0
    n_sorption = 1 # Henry's law
    temperature = ${temperature}
    diffusivity = diffusivity_elemental_tritium
    boundary = 'interface'
  []
  [hto_solubility]
    type = ADInterfaceSorption
    variable = hto_paint_concentration
    neighbor_var = hto_enclosure_concentration # molecules/microns^3
    unit_scale_neighbor = ${fparse 1e18/NA} # to convert neighbor concentration to mol/m^3
    K0 = ${solubility_tritiated_water}
    Ea = 0
    n_sorption = 1 # Henry's law
    temperature = ${temperature}
    diffusivity = diffusivity_tritiated_water
    boundary = 'interface'
  []
  [h2o_solubility]
    type = ADInterfaceSorption
    variable = h2o_paint_concentration
    neighbor_var = h2o_enclosure_concentration # molecules/microns^3
    unit_scale_neighbor = ${fparse 1e18/NA} # to convert neighbor concentration to mol/m^3
    K0 = ${solubility_tritiated_water}
    Ea = 0
    n_sorption = 1 # Henry's law
    temperature = ${temperature}
    diffusivity = diffusivity_tritiated_water
    boundary = 'interface'
  []
[]

[Materials]
  [reaction_rate_t2]
    type = ADDerivativeParsedMaterial
    coupled_variables = 't2_enclosure_concentration ht_enclosure_concentration hto_enclosure_concentration'
    expression = '${reaction_rate} * 2 * t2_enclosure_concentration * (2*t2_enclosure_concentration + ht_enclosure_concentration + hto_enclosure_concentration)'
    property_name = reaction_rate_t2
    block = 1
  []
  [reaction_rate_ht]
    type = ADDerivativeParsedMaterial
    coupled_variables = 't2_enclosure_concentration ht_enclosure_concentration hto_enclosure_concentration'
    expression = '${reaction_rate} * ht_enclosure_concentration * (2*t2_enclosure_concentration + ht_enclosure_concentration + hto_enclosure_concentration)'
    property_name = reaction_rate_ht
    block = 1
  []
  [diffusivity_elemental_tritium_paint]
    type = ADConstantMaterial
    value = ${diffusivity_elemental_tritium}
    property_name = 'diffusivity_elemental_tritium'
    block = 0
  []
  [diffusivity_elemental_tritium_enclosure]
    type = ADConstantMaterial
    value = ${diffusivity_artificial_enclosure}
    property_name = 'diffusivity_elemental_tritium'
    block = 1
  []
  [diffusivity_tritiated_water_paint]
    type = ADConstantMaterial
    value = ${diffusivity_tritiated_water}
    property_name = 'diffusivity_tritiated_water'
    block = 0
  []
  [diffusivity_tritiated_water_enclosure]
    type = ADConstantMaterial
    value = ${diffusivity_artificial_enclosure}
    property_name = 'diffusivity_tritiated_water'
    block = 1
  []

  [t2_enclosure_concentration_outflow]
    type = DerivativeParsedMaterial
    coupled_variables = 't2_enclosure_concentration'
    expression = 't2_enclosure_concentration * ${outflow} / ${volume_enclosure}'
    property_name = 't2_enclosure_concentration_outflow'
    block = 1
  []
  [ht_enclosure_concentration_outflow]
    type = DerivativeParsedMaterial
    coupled_variables = 'ht_enclosure_concentration'
    expression = 'ht_enclosure_concentration * ${outflow} / ${volume_enclosure}'
    property_name = 'ht_enclosure_concentration_outflow'
    block = 1
  []
  [hto_enclosure_concentration_outflow]
    type = DerivativeParsedMaterial
    coupled_variables = 'hto_enclosure_concentration'
    expression = 'hto_enclosure_concentration * ${outflow} / ${volume_enclosure}'
    property_name = 'hto_enclosure_concentration_outflow'
    block = 1
  []
  [h2o_enclosure_concentration_outflow]
    type = DerivativeParsedMaterial
    coupled_variables = 'h2o_enclosure_concentration'
    expression = 'h2o_enclosure_concentration * ${outflow} / ${volume_enclosure}'
    property_name = 'h2o_enclosure_concentration_outflow'
    block = 1
  []
[]

[Postprocessors]
  # Pressures in enclosure
  [t2_enclosure_edge_concentration] # (molecules/microns^3)
    type = PointValue
    point = '${length_domain} 0 0' # on the far side of the enclosure
    variable = t2_enclosure_concentration
    execute_on = 'initial timestep_end'
  []
  [ht_enclosure_edge_concentration] # (molecules/microns^3)
    type = PointValue
    point = '${length_domain} 0 0' # on the far side of the enclosure
    variable = ht_enclosure_concentration
    execute_on = 'initial timestep_end'
  []
  [hto_enclosure_edge_concentration] # (molecules/microns^3)
    type = PointValue
    point = '${length_domain} 0 0' # on the far side of the enclosure
    variable = hto_enclosure_concentration
    execute_on = 'initial timestep_end'
  []
  [h2o_enclosure_edge_concentration] # (molecules/microns^3)
    type = PointValue
    point = '${length_domain} 0 0' # on the far side of the enclosure
    variable = h2o_enclosure_concentration
    execute_on = 'initial timestep_end'
  []

  # Inventory in enclosure
  [t2_enclosure_inventory] # (molecules/microns^2)
    type = ElementIntegralVariablePostprocessor
    variable = t2_enclosure_concentration
    execute_on = 'initial timestep_end'
    block = 1
  []
  [ht_enclosure_inventory] # (molecules/microns^2)
    type = ElementIntegralVariablePostprocessor
    variable = ht_enclosure_concentration
    execute_on = 'initial timestep_end'
    block = 1
  []
  [hto_enclosure_inventory] # (molecules/microns^2)
    type = ElementIntegralVariablePostprocessor
    variable = hto_enclosure_concentration
    execute_on = 'initial timestep_end'
    block = 1
  []
  [h2o_enclosure_inventory] # (molecules/microns^2)
    type = ElementIntegralVariablePostprocessor
    variable = h2o_enclosure_concentration
    execute_on = 'initial timestep_end'
    block = 1
  []

  # Inventory in paint
  [t2_paint_inventory] # (molecules/microns^2)
    type = ElementIntegralVariablePostprocessor
    variable = t2_paint_concentration
    execute_on = 'initial timestep_end'
    block = 0
  []
  [ht_paint_inventory] # (molecules/microns^2)
    type = ElementIntegralVariablePostprocessor
    variable = ht_paint_concentration
    execute_on = 'initial timestep_end'
    block = 0
  []
  [hto_paint_inventory] # (molecules/microns^2)
    type = ElementIntegralVariablePostprocessor
    variable = hto_paint_concentration
    execute_on = 'initial timestep_end'
    block = 0
  []
  [h2o_paint_inventory] # (molecules/microns^2)
    type = ElementIntegralVariablePostprocessor
    variable = h2o_paint_concentration
    execute_on = 'initial timestep_end'
    block = 0
  []

  [tritium_total_inventory]
    type = LinearCombinationPostprocessor
    pp_names = 't2_enclosure_inventory ht_enclosure_inventory hto_enclosure_inventory t2_paint_inventory ht_paint_inventory hto_paint_inventory'
    pp_coefs = '2                      1                      1                       2                  1                  1'
    execute_on = 'initial timestep_end'
  []
[]

[Debug]
  show_var_residual_norms = true
[]

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

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

  petsc_options_iname = '-pc_type -sub_pc_type -snes_type'
  petsc_options_value = 'asm      lu           vinewtonrsls' # This petsc option helps prevent negative concentrations with bounds'

  nl_rel_tol = 1.e-10
  automatic_scaling = true
  compute_scaling_once = false

  end_time = ${time_end}
  dtmax = ${dtmax}
  nl_max_its = 16
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = ${time_step}
    optimal_iterations = 12
    iteration_window = 1
    growth_factor = 1.1
    cutback_factor = 0.9
  []
[]

[Outputs]
  exodus = true
  perf_graph = true
  [csv]
    type = CSV
    execute_on = 'initial timestep_end'
  []
[]

(test/tests/ver-1kb/ver-1kb.i)

nb_segments_TMAP7 = 20
node_size_TMAP7 = '${units 1.25e-5 m}'
long_total = '${fparse nb_segments_TMAP7 * node_size_TMAP7}' # m
nb_segments_TMAP8 = 100
simulation_time = '${units 10 s}'
temperature = '${units 500 K}'
R = '${units 8.31446261815324 J/mol/K}' # ideal gas constant from PhysicalConstants.h
initial_pressure_1 = '${units 1e5 Pa}'
initial_pressure_2 = '${units 1e-10 Pa}'
initial_concentration_1 = '${units ${fparse initial_pressure_1 / (R*temperature)} mol/m^3}'
initial_concentration_2 = '${units ${fparse initial_pressure_2 / (R*temperature)} mol/m^3}'
solubility = '${units ${fparse 1/(R*temperature)} mol/m^3/Pa}' # Henry's law solubility
diffusivity = '${units ${fparse 4.31e-6 * exp(-2818/temperature)} m^2/s}'
n_sorption = 1 # Henry's Law
unit_scale = 1
unit_scale_neighbor = 1

[Mesh]
  [generated]
    type = GeneratedMeshGenerator
    dim = 1
    nx = ${nb_segments_TMAP8}
    xmax = ${long_total}
  []
  [enclosure_1]
    type = SubdomainBoundingBoxGenerator
    input = generated
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '${fparse 1/3 * long_total} 0 0'
  []
  [enclosure_2]
    type = SubdomainBoundingBoxGenerator
    input = enclosure_1
    block_id = 2
    bottom_left = '${fparse 1/3 * long_total} 0 0'
    top_right = '${fparse long_total} 0 0'
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = enclosure_2
    primary_block = 1
    paired_block = 2
    new_boundary = interface
  []
  [interface2]
    type = SideSetsBetweenSubdomainsGenerator
    input = interface
    primary_block = 2
    paired_block = 1
    new_boundary = interface2
  []
[]

[Variables]
  [concentration_enclosure_1]
    block = 1
    order = FIRST
    family = LAGRANGE
    initial_condition = '${initial_concentration_1}'
  []
  [concentration_enclosure_2]
    block = 2
    order = FIRST
    family = LAGRANGE
    initial_condition = '${initial_concentration_2}'
  []
[]

[AuxVariables]
  [pressure_enclosure_1]
    order = CONSTANT
    family = MONOMIAL
    block = 1
    initial_condition = '${initial_pressure_1}'
  []
  [pressure_enclosure_2]
    order = CONSTANT
    family = MONOMIAL
    block = 2
    initial_condition = '${initial_pressure_2}'
  []
[]

[Kernels]
  [diffusion_enclosure_1]
    type = MatDiffusion
    variable = concentration_enclosure_1
    diffusivity = ${diffusivity}
    block = '1'
  []
  [time_enclosure_1]
    type = TimeDerivative
    variable = concentration_enclosure_1
    block = '1'
  []
  [diffusion_enclosure_2]
    type = MatDiffusion
    variable = concentration_enclosure_2
    diffusivity = ${diffusivity}
    block = '2'
  []
  [time_enclosure_2]
    type = TimeDerivative
    variable = concentration_enclosure_2
    block = '2'
  []
[]

[AuxKernels]
  [pressure_enclosure_1]
    type = ParsedAux
    variable = pressure_enclosure_1
    coupled_variables = 'concentration_enclosure_1'
    expression = '${fparse R*temperature}*concentration_enclosure_1'
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_enclosure_2]
    type = ParsedAux
    variable = pressure_enclosure_2
    coupled_variables = 'concentration_enclosure_2'
    expression = '${fparse R*temperature}*concentration_enclosure_2'
    block = 2
    execute_on = 'initial timestep_end'
  []
[]

[InterfaceKernels]
  [interface_sorption]
    type = InterfaceSorption
    K0 = ${solubility}
    Ea = 0
    n_sorption = ${n_sorption}
    diffusivity = ${diffusivity}
    unit_scale = ${unit_scale}
    unit_scale_neighbor = ${unit_scale_neighbor}
    temperature = ${temperature}
    variable = concentration_enclosure_1
    neighbor_var = concentration_enclosure_2
    sorption_penalty = 1e1
    boundary = interface
  []
[]

[Postprocessors]
  [pressure_enclosure_1]
    type = ElementAverageValue
    variable = pressure_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_enclosure_2]
    type = ElementAverageValue
    variable = pressure_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []
  [concentration_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = concentration_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface2
    variable = concentration_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_ratio]
    type = ParsedPostprocessor
    expression = 'concentration_enclosure_1_at_interface / concentration_enclosure_2_at_interface'
    pp_names = 'concentration_enclosure_1_at_interface concentration_enclosure_2_at_interface'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
  [solubility]
    type = ParsedPostprocessor
    expression = 'concentration_enclosure_1_at_interface / pressure_enclosure_2_at_interface'
    pp_names = 'concentration_enclosure_1_at_interface pressure_enclosure_2_at_interface'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
  [pressure_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = pressure_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [pressure_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface2
    variable = pressure_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_encl_1_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [concentration_encl_2_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []
  [mass_conservation_sum_encl1_encl2]
    type = LinearCombinationPostprocessor
    pp_names = 'concentration_encl_1_inventory concentration_encl_2_inventory'
    pp_coefs = '1            1'
    execute_on = 'initial timestep_end'
  []
[]

[Executioner]
  type = Transient
  end_time = ${simulation_time}
  dtmax = 10
  nl_abs_tol = 1e-12
  nl_rel_tol = 1e-8
  scheme = 'bdf2'
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  nl_max_its = 6
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.01
    optimal_iterations = 5
    iteration_window = 3
    growth_factor = 1.05
    cutback_factor_at_failure = 0.8
  []
[]

[Outputs]
  file_base = 'ver-1kb_out_k1'
  csv = true
  exodus = true
  execute_on = 'initial timestep_end'
[]

(test/tests/interfacekernels/InterfaceSorption/interface_sorption.i)

# This test is to verify the implementation of InterfaceSorption and its AD counterpart.
# It contains two 1D blocks separated by a continuous interface.
# InterfaceSorption is used to enforce the sorption law and preserve flux between the blocks.
# Checks are performed to verify concentration conservation, sorption behavior, and flux preservation.
# This input file uses BreakMeshByBlockGenerator, which is currently only supported for replicated
# meshes, so this file should not be run with the `parallel_type = DISTRIBUTED` flag

# In this input file, we apply the Sievert law with n_sorption=1/2.

# Physical Constants
R = 8.31446261815324 # J/mol/K, based on number used in include/utils/PhysicalConstants.h

solubility = 1.e-2
n_sorption = 0.5

unit_scale = 1
unit_scale_neighbor = 1


[GlobalParams]
  order = FIRST
  family = LAGRANGE
[]

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    nx = 10
    dim = 1
  []
  [block1]
    type = SubdomainBoundingBoxGenerator
    input = gen
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '0.5 1 0'
  []
  [block2]
    type = SubdomainBoundingBoxGenerator
    input = block1
    block_id = 2
    bottom_left = '0.5 0 0'
    top_right = '1 1 0'
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = block2
    primary_block = 1
    paired_block = 2
    new_boundary = interface
  []
  [interface2]
    type = SideSetsBetweenSubdomainsGenerator
    input = interface
    primary_block = 2
    paired_block = 1
    new_boundary = interface2
  []
[]

[Variables]
  [u1]
    block = 1
  []
  [u2]
    block = 2
  []
  [temperature]
  []
[]

[Kernels]
  [u1]
    type = MatDiffusion
    variable = u1
    diffusivity = diffusivity
    block = 1
  []
  [u2]
    type = MatDiffusion
    variable = u2
    diffusivity = diffusivity
    block = 2
  []
  [temperature]
    type = HeatConduction
    variable = temperature
  []
[]

[BCs]
  [left_u1]
    type = DirichletBC
    value = 1
    variable = u1
    boundary = left
  []
  [right_u2]
    type = DirichletBC
    value = 1
    variable = u2
    boundary = right
  []
  [left_temperature]
    type = DirichletBC
    value = 1100
    variable = temperature
    boundary = left
  []
  [right_temperature]
    type = DirichletBC
    value = 0
    variable = temperature
    boundary = right
  []
  [block1_2_temperature]
    type = DirichletBC
    value = 1000
    variable = temperature
    boundary = interface
  []
  [block2_1_temperature]
    type = DirichletBC
    value = 900
    variable = temperature
    boundary = interface2
  []
[]

[InterfaceKernels]
  [interface]
    type = InterfaceSorption
    K0 = ${solubility}
    Ea = 0
    n_sorption = ${n_sorption}
    diffusivity = diffusivity
    unit_scale = ${unit_scale}
    unit_scale_neighbor = ${unit_scale_neighbor}
    temperature = temperature
    variable = u2
    neighbor_var = u1
    sorption_penalty = 1e1
    boundary = interface2
  []
[]

[Materials]
  [properties_1]
    type = GenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity'
    prop_values = '1 1'
    block = 1
  []
  [properties_2]
    type = GenericConstantMaterial
    prop_names = 'thermal_conductivity diffusivity solubility'
    prop_values = '2 2 ${solubility}'
    block = 2
  []
[]

[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)'
  []
  [residual_concentration]
    type = ParsedFunction
    symbol_names = 'u_mid_inner u_mid_outer T'
    symbol_values = 'u_mid_inner u_mid_outer temperature_mid_inner'
    expression = 'u_mid_outer*${unit_scale} - ${solubility}*(u_mid_inner*${unit_scale_neighbor}*${R}*T)^${n_sorption}'
  []
  [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

  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
[]

[Postprocessors]
  [u_mid_inner]
    type = PointValue
    variable = u1
    point = '0.49999 0 0'
    outputs = 'csv console'
  []
  [u_mid_outer]
    type = PointValue
    variable = u2
    point = '0.50001 0 0'
    outputs = 'csv console'
  []
  [u_mid_diff]
    type = FunctionValuePostprocessor
    function = u_mid_diff
    outputs = 'csv console'
  []
  [temperature_mid_inner]
    type = PointValue
    variable = temperature
    point = '0.49999 0 0'
    outputs = csv
  []
  [temperature_mid_outer]
    type = PointValue
    variable = temperature
    point = '0.50001 0 0'
    outputs = csv
  []
  [residual_concentration]
    type = FunctionValuePostprocessor
    function = residual_concentration
    outputs = 'csv console'
  []

  [flux_inner] # verify flux preservation
    type = SideDiffusiveFluxIntegral
    variable = u1
    boundary = interface
    diffusivity = diffusivity
    outputs = 'csv console'
  []
  [flux_outer]
    type = SideDiffusiveFluxIntegral
    variable = u2
    boundary = interface2
    diffusivity = diffusivity
    outputs = 'csv console'
  []
  [flux_error]
    type = FunctionValuePostprocessor
    function = flux_error
    outputs = 'csv console'
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

[Debug]
  show_var_residual_norms = true
[]

(test/tests/interfacekernels/InterfaceSorption/interface_sorption_transient.i)

# This test is to verify the implementation of InterfaceSorption and its AD counterpart in transient conditions.
# It contains two 1D blocks separated by a continuous interface.
# InterfaceSorption is used to enforce the sorption law and preserve flux between the blocks.
# Checks are performed to verify concentration conservation, sorption behavior, and flux preservation.
# This input file uses BreakMeshByBlockGenerator, which is currently only supported for replicated
# meshes, so this file should not be run with the `parallel_type = DISTRIBUTED` flag

# In this input file, we apply the Sievert law with n_sorption=1/2.

# Physical Constants
R = 8.31446261815324 # J/mol/K, based on number used in include/utils/PhysicalConstants.h

temperature = 1000 # K
initial_concentration = 1
n_sorption = 0.5
solubility = ${fparse 2/R^n_sorption/temperature^n_sorption}

unit_scale = 1
unit_scale_neighbor = 1


[GlobalParams]
  order = FIRST
  family = LAGRANGE
[]

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    nx = 10
    dim = 1
  []
  [block1]
    type = SubdomainBoundingBoxGenerator
    input = gen
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '0.5 1 0'
  []
  [block2]
    type = SubdomainBoundingBoxGenerator
    input = block1
    block_id = 2
    bottom_left = '0.5 0 0'
    top_right = '1 1 0'
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = block2
    primary_block = 1
    paired_block = 2
    new_boundary = interface
  []
  [interface2]
    type = SideSetsBetweenSubdomainsGenerator
    input = interface
    primary_block = 2
    paired_block = 1
    new_boundary = interface2
  []
[]

[Variables]
  [u1]
    block = 1
    initial_condition = ${initial_concentration}
  []
  [u2]
    block = 2
    initial_condition = ${initial_concentration}
  []
  [temperature]
    initial_condition = ${temperature}
  []
[]

[Kernels]
  [u1_time_derivative]
    type = TimeDerivative
    variable = u1
    block = 1
  []
  [u1_diffusion]
    type = MatDiffusion
    variable = u1
    diffusivity = diffusivity
    block = 1
  []
  [u2_time_derivative]
    type = TimeDerivative
    variable = u2
    block = 2
  []
  [u2_diffusion]
    type = MatDiffusion
    variable = u2
    diffusivity = diffusivity
    block = 2
  []
  [temperature]
    type = TimeDerivative
    variable = temperature
  []
[]

[BCs]
  [left_u1]
    type = NeumannBC
    value = 0
    variable = u1
    boundary = left
  []
  [right_u2]
    type = NeumannBC
    value = 0
    variable = u2
    boundary = right
  []
[]

[InterfaceKernels]
  [interface]
    type = InterfaceSorption
    K0 = ${solubility}
    Ea = 0
    n_sorption = ${n_sorption}
    diffusivity = diffusivity
    unit_scale = ${unit_scale}
    unit_scale_neighbor = ${unit_scale_neighbor}
    temperature = temperature
    variable = u2
    neighbor_var = u1
    sorption_penalty = 1e1
    boundary = interface2
  []
[]

[Materials]
  [properties_1]
    type = GenericConstantMaterial
    prop_names = 'diffusivity thermal_conductivity '
    prop_values = '0.2 1'
    block = '1'
  []
  [properties_2]
    type = GenericConstantMaterial
    prop_names = 'diffusivity solubility thermal_conductivity'
    prop_values = '0.2 ${solubility} 1'
    block = 2
  []
[]

[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)'
  []
  [residual_concentration]
    type = ParsedFunction
    symbol_names = 'u_mid_inner u_mid_outer'
    symbol_values = 'u_mid_inner u_mid_outer '
    expression = 'u_mid_outer*${unit_scale} - ${solubility}*(u_mid_inner*${unit_scale_neighbor}*${R}*${temperature})^${n_sorption}'
  []
  [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 = Transient
  solve_type = NEWTON

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  dt = 0.1
  end_time = 5
[]

[Postprocessors]
  [u_mid_inner]
    type = PointValue
    variable = u1
    point = '0.49999 0 0'
    outputs = 'csv console'
  []
  [u_mid_outer]
    type = PointValue
    variable = u2
    point = '0.50001 0 0'
    outputs = 'csv console'
  []
  [u_mid_diff]
    type = FunctionValuePostprocessor
    function = u_mid_diff
    outputs = 'csv console'
  []
  [u1_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = u1
    block = 1
  []
  [u2_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = u2
    block = 2
  []
  [mass_conservation_sum_u1_u2]
    type = LinearCombinationPostprocessor
    pp_names = 'u1_inventory u2_inventory'
    pp_coefs = '1            1'
  []
  [residual_concentration]
    type = FunctionValuePostprocessor
    function = residual_concentration
    outputs = 'csv console'
  []

  [flux_inner] # verify flux preservation
    type = SideDiffusiveFluxIntegral
    variable = u1
    boundary = interface
    diffusivity = diffusivity
    outputs = 'csv console'
  []
  [flux_outer]
    type = SideDiffusiveFluxIntegral
    variable = u2
    boundary = interface2
    diffusivity = diffusivity
    outputs = 'csv console'
  []
  [flux_error]
    type = FunctionValuePostprocessor
    function = flux_error
    outputs = 'csv console'
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

[Debug]
  show_var_residual_norms = true
[]

(test/tests/ver-1kd/ver-1kd.i)

nb_segments_TMAP7 = 20
node_size_TMAP7 = '${units 1.25e-5 m}'
long_total = '${units ${fparse nb_segments_TMAP7 * node_size_TMAP7} m}'
nb_segments_TMAP8 = 1e2
simulation_time = '${units 2 s}'
N_a = '${units 6.02214076e23 1/mol}' # Avogadro's number from PhysicalConstants.h
source = '${units ${fparse 1e23 / N_a} mol/m^3/s}' # Source term for T2 in enclosure 1
temperature = '${units 500 K}'
R = '${units 8.31446261815324 J/mol/K}' # ideal gas constant from PhysicalConstants.h
initial_pressure_1 = '${units 1e5 Pa}'
initial_pressure_2 = '${units 1e-10 Pa}'
initial_concentration_1 = '${units ${fparse initial_pressure_1 / (R*temperature)} mol/m^3}'
initial_concentration_2 = '${units ${fparse initial_pressure_2 / (R*temperature)} mol/m^3}'
solubility = '${units ${fparse 10/sqrt(R*temperature)} mol/m^3/Pa^(1/2)}' # Sieverts' law solubility
diffusivity = '${units ${fparse 4.31e-6 * exp(-2818/temperature)} m^2/s}'
n_sorption = 0.5 # Sieverts' Law
K1 = '${units 4000 mol/m^3/s}' # reaction rate for H2+T2->2HT
equilibrium_constant = 2.0
K2 = '${units ${fparse (2*K1) / (equilibrium_constant)^2} mol/m^3/s}' # reaction rate for 2HT->H2+T2
unit_scale = 1
unit_scale_neighbor = 1

[Mesh]
  [generated]
    type = GeneratedMeshGenerator
    dim = 1
    nx = ${nb_segments_TMAP8}
    xmax = ${long_total}
  []
  [enclosure_1]
    type = SubdomainBoundingBoxGenerator
    input = generated
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '${fparse 1/3 * long_total} 0 0'
  []
  [enclosure_2]
    type = SubdomainBoundingBoxGenerator
    input = enclosure_1
    block_id = 2
    bottom_left = '${fparse 1/3 * long_total} 0 0'
    top_right = '${fparse long_total} 0 0'
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = enclosure_2
    primary_block = 1
    paired_block = 2
    new_boundary = interface
  []
  [interface2]
    type = SideSetsBetweenSubdomainsGenerator
    input = interface
    primary_block = 2
    paired_block = 1
    new_boundary = interface2
  []
[]

[Variables]

  # Variables for H2, T2, and HT in enclosure 1
  [concentration_H2_enclosure_1]
    block = 1
    initial_condition = '${initial_concentration_1}'
  []
  [concentration_T2_enclosure_1]
    block = 1
    initial_condition = '${initial_concentration_1}'
  []
  [concentration_HT_enclosure_1]
    block = 1
    initial_condition = '${initial_concentration_2}'
  []

  # Variables for H2, T2, and HT in enclosure 2
  [concentration_H2_enclosure_2]
    block = 2
    initial_condition = '${initial_concentration_2}'
  []
  [concentration_T2_enclosure_2]
    block = 2
    initial_condition = '${initial_concentration_2}'
  []
  [concentration_HT_enclosure_2]
    block = 2
    initial_condition = '${initial_concentration_2}'
  []
[]

[AuxVariables]

  # Auxiliary variables for H2, T2, and HT in enclosure 1
  [pressure_H2_enclosure_1]
    order = CONSTANT
    family = MONOMIAL
    block = 1
    initial_condition = '${initial_pressure_1}'
  []
  [pressure_T2_enclosure_1]
    order = CONSTANT
    family = MONOMIAL
    block = 1
    initial_condition = '${initial_pressure_1}'
  []
  [pressure_HT_enclosure_1]
    order = CONSTANT
    family = MONOMIAL
    block = 1
    initial_condition = '${initial_pressure_2}'
  []

  # Auxiliary variables for H2, T2, and HT in enclosure 2
  [pressure_H2_enclosure_2]
    order = CONSTANT
    family = MONOMIAL
    block = 2
    initial_condition = '${initial_pressure_2}'
  []
  [pressure_T2_enclosure_2]
    order = CONSTANT
    family = MONOMIAL
    block = 2
    initial_condition = '${initial_pressure_2}'
  []
  [pressure_HT_enclosure_2]
    order = CONSTANT
    family = MONOMIAL
    block = 2
    initial_condition = '${initial_pressure_2}'
  []
[]

[Kernels]
  # Source term for T2 in enclosure 1
  [T2_source_term_enclosure_1]
    type = BodyForce
    variable = concentration_T2_enclosure_1
    block = '1'
    value = ${source}
  []

  # Diffusion equation for H2
  [H2_diffusion_enclosure_1]
    type = MatDiffusion
    variable = concentration_H2_enclosure_1
    diffusivity = ${diffusivity}
    block = '1'
  []
  [H2_time_derivative_enclosure_1]
    type = TimeDerivative
    variable = concentration_H2_enclosure_1
    block = '1'
  []
  [H2_diffusion_enclosure_2]
    type = MatDiffusion
    variable = concentration_H2_enclosure_2
    diffusivity = ${diffusivity}
    block = '2'
  []
  [H2_time_derivative_enclosure_2]
    type = TimeDerivative
    variable = concentration_H2_enclosure_2
    block = '2'
  []

  # Diffusion equation for T2
  [T2_diffusion_enclosure_1]
    type = MatDiffusion
    variable = concentration_T2_enclosure_1
    diffusivity = ${diffusivity}
    block = '1'
  []
  [T2_time_derivative_enclosure_1]
    type = TimeDerivative
    variable = concentration_T2_enclosure_1
    block = '1'
  []
  [T2_diffusion_enclosure_2]
    type = MatDiffusion
    variable = concentration_T2_enclosure_2
    diffusivity = ${diffusivity}
    block = '2'
  []
  [T2_time_derivative_enclosure_2]
    type = TimeDerivative
    variable = concentration_T2_enclosure_2
    block = '2'
  []

  # Diffusion equation for HT
  [HT_diffusion_enclosure_1]
    type = MatDiffusion
    variable = concentration_HT_enclosure_1
    diffusivity = ${diffusivity}
    block = '1'
  []
  [HT_time_derivative_enclosure_1]
    type = TimeDerivative
    variable = concentration_HT_enclosure_1
    block = '1'
  []
  [HT_diffusion_enclosure_2]
    type = MatDiffusion
    variable = concentration_HT_enclosure_2
    diffusivity = ${diffusivity}
    block = '2'
  []
  [HT_time_derivative_enclosure_2]
    type = TimeDerivative
    variable = concentration_HT_enclosure_2
    block = '2'
  []

  # Reaction H2+T2->2HT in enclosure 1
  [reaction_H2_encl_1_1]
    type = ADMatReactionFlexible
    variable = concentration_H2_enclosure_1
    vs = 'concentration_H2_enclosure_1 concentration_T2_enclosure_1'
    block = 1
    coeff = -1
    reaction_rate_name = ${K1}
  []
  [reaction_T2_encl_1_1]
    type = ADMatReactionFlexible
    variable = concentration_T2_enclosure_1
    vs = 'concentration_H2_enclosure_1 concentration_T2_enclosure_1'
    block = 1
    coeff = -1
    reaction_rate_name = ${K1}
  []
  [reaction_HT_encl_1_1]
    type = ADMatReactionFlexible
    variable = concentration_HT_enclosure_1
    vs = 'concentration_H2_enclosure_1 concentration_T2_enclosure_1'
    block = 1
    coeff = 2
    reaction_rate_name = ${K1}
  []

  # Reaction 2HT->H2+T2 in enclosure 1
  [reaction_H2_encl_1_2]
    type = ADMatReactionFlexible
    variable = concentration_H2_enclosure_1
    vs = 'concentration_HT_enclosure_1 concentration_HT_enclosure_1'
    block = 1
    coeff = 0.5
    reaction_rate_name = ${K2}
  []
  [reaction_T2_encl_1_2]
    type = ADMatReactionFlexible
    variable = concentration_T2_enclosure_1
    vs = 'concentration_HT_enclosure_1 concentration_HT_enclosure_1'
    block = 1
    coeff = 0.5
    reaction_rate_name = ${K2}
  []
  [reaction_HT_encl_1_2]
    type = ADMatReactionFlexible
    variable = concentration_HT_enclosure_1
    vs = 'concentration_HT_enclosure_1 concentration_HT_enclosure_1'
    block = 1
    coeff = -1
    reaction_rate_name = ${K2}
  []

  # Reaction H2+T2->2HT in enclosure 2
  [reaction_H2_encl_2_1]
    type = ADMatReactionFlexible
    variable = concentration_H2_enclosure_2
    vs = 'concentration_H2_enclosure_2 concentration_T2_enclosure_2'
    block = 2
    coeff = -1
    reaction_rate_name = ${K1}
  []
  [reaction_T2_encl_2_1]
    type = ADMatReactionFlexible
    variable = concentration_T2_enclosure_2
    vs = 'concentration_H2_enclosure_2 concentration_T2_enclosure_2'
    block = 2
    coeff = -1
    reaction_rate_name = ${K1}
  []
  [reaction_HT_encl_2_1]
    type = ADMatReactionFlexible
    variable = concentration_HT_enclosure_2
    vs = 'concentration_H2_enclosure_2 concentration_T2_enclosure_2'
    block = 2
    coeff = 2
    reaction_rate_name = ${K1}
  []

  # Reaction 2HT->H2+T2 in enclosure 2
  [reaction_H2_encl_2_2]
    type = ADMatReactionFlexible
    variable = concentration_H2_enclosure_2
    vs = 'concentration_HT_enclosure_2 concentration_HT_enclosure_2'
    block = 2
    coeff = 0.5
    reaction_rate_name = ${K2}
  []
  [reaction_T2_encl_2_2]
    type = ADMatReactionFlexible
    variable = concentration_T2_enclosure_2
    vs = 'concentration_HT_enclosure_2 concentration_HT_enclosure_2'
    block = 2
    coeff = 0.5
    reaction_rate_name = ${K2}
  []
  [reaction_HT_encl_2_2]
    type = ADMatReactionFlexible
    variable = concentration_HT_enclosure_2
    vs = 'concentration_HT_enclosure_2 concentration_HT_enclosure_2'
    block = 2
    coeff = -1
    reaction_rate_name = ${K2}
  []
[]

[AuxKernels]

  # Auxiliary kernels for H2, T2, and HT in enclosure 1
  [pressure_H2_enclosure_1]
    type = ParsedAux
    variable = pressure_H2_enclosure_1
    coupled_variables = 'concentration_H2_enclosure_1'
    expression = '${fparse R*temperature}*concentration_H2_enclosure_1'
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_T2_enclosure_1]
    type = ParsedAux
    variable = pressure_T2_enclosure_1
    coupled_variables = 'concentration_T2_enclosure_1'
    expression = '${fparse R*temperature}*concentration_T2_enclosure_1'
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_HT_enclosure_1]
    type = ParsedAux
    variable = pressure_HT_enclosure_1
    coupled_variables = 'concentration_HT_enclosure_1'
    expression = '${fparse R*temperature}*concentration_HT_enclosure_1'
    block = 1
    execute_on = 'initial timestep_end'
  []

  # Auxiliary kernels for H2, T2, and HT in enclosure 2
  [pressure_H2_enclosure_2]
    type = ParsedAux
    variable = pressure_H2_enclosure_2
    coupled_variables = 'concentration_H2_enclosure_2'
    expression = '${fparse R*temperature}*concentration_H2_enclosure_2'
    block = 2
    execute_on = 'initial timestep_end'
  []
  [pressure_T2_enclosure_2]
    type = ParsedAux
    variable = pressure_T2_enclosure_2
    coupled_variables = 'concentration_T2_enclosure_2'
    expression = '${fparse R*temperature}*concentration_T2_enclosure_2'
    block = 2
    execute_on = 'initial timestep_end'
  []
  [pressure_HT_enclosure_2]
    type = ParsedAux
    variable = pressure_HT_enclosure_2
    coupled_variables = 'concentration_HT_enclosure_2'
    expression = '${fparse R*temperature}*concentration_HT_enclosure_2'
    block = 2
    execute_on = 'initial timestep_end'
  []
[]

[InterfaceKernels]
  [interface_sorption_H2]
    type = InterfaceSorption
    K0 = ${solubility}
    Ea = 0
    n_sorption = ${n_sorption}
    diffusivity = ${diffusivity}
    unit_scale = ${unit_scale}
    unit_scale_neighbor = ${unit_scale_neighbor}
    temperature = ${temperature}
    variable = concentration_H2_enclosure_1
    neighbor_var = concentration_H2_enclosure_2
    sorption_penalty = 1e1
    boundary = interface
  []
  [interface_sorption_T2]
    type = InterfaceSorption
    K0 = ${solubility}
    Ea = 0
    n_sorption = ${n_sorption}
    diffusivity = ${diffusivity}
    unit_scale = ${unit_scale}
    unit_scale_neighbor = ${unit_scale_neighbor}
    temperature = ${temperature}
    variable = concentration_T2_enclosure_1
    neighbor_var = concentration_T2_enclosure_2
    sorption_penalty = 1e1
    boundary = interface
  []
  [interface_sorption_HT]
    type = InterfaceSorption
    K0 = ${solubility}
    Ea = 0
    n_sorption = ${n_sorption}
    diffusivity = ${diffusivity}
    unit_scale = ${unit_scale}
    unit_scale_neighbor = ${unit_scale_neighbor}
    temperature = ${temperature}
    variable = concentration_HT_enclosure_1
    neighbor_var = concentration_HT_enclosure_2
    sorption_penalty = 1e1
    boundary = interface
  []
[]

[Postprocessors]

  # postprocessors for H2
  [pressure_H2_enclosure_1]
    type = ElementAverageValue
    variable = pressure_H2_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_H2_enclosure_2]
    type = ElementAverageValue
    variable = pressure_H2_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []
  [concentration_H2_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = concentration_H2_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_H2_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface2
    variable = concentration_H2_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_ratio_H2]
    type = ParsedPostprocessor
    expression = 'concentration_H2_enclosure_1_at_interface / sqrt(concentration_H2_enclosure_2_at_interface)'
    pp_names = 'concentration_H2_enclosure_1_at_interface concentration_H2_enclosure_2_at_interface'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
  [pressure_H2_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = pressure_H2_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [pressure_H2_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface2
    variable = pressure_H2_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_H2_encl_1_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_H2_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [concentration_H2_encl_2_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_H2_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []

  # postprocessors for T2
  [pressure_T2_enclosure_1]
    type = ElementAverageValue
    variable = pressure_T2_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_T2_enclosure_2]
    type = ElementAverageValue
    variable = pressure_T2_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []
  [concentration_T2_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = concentration_T2_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_T2_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface2
    variable = concentration_T2_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_ratio_T2]
    type = ParsedPostprocessor
    expression = 'concentration_T2_enclosure_1_at_interface / sqrt(concentration_T2_enclosure_2_at_interface)'
    pp_names = 'concentration_T2_enclosure_1_at_interface concentration_T2_enclosure_2_at_interface'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
  [concentration_T2_encl_1_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_T2_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [concentration_T2_encl_2_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_T2_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []

  # postprocessors for HT
  [pressure_HT_enclosure_1]
    type = ElementAverageValue
    variable = pressure_HT_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_HT_enclosure_2]
    type = ElementAverageValue
    variable = pressure_HT_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []
  [concentration_HT_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = concentration_HT_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_HT_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface2
    variable = concentration_HT_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_ratio_HT]
    type = ParsedPostprocessor
    expression = 'concentration_HT_enclosure_1_at_interface / sqrt(concentration_HT_enclosure_2_at_interface)'
    pp_names = 'concentration_HT_enclosure_1_at_interface concentration_HT_enclosure_2_at_interface'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
  [concentration_HT_encl_1_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_HT_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [concentration_HT_encl_2_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_HT_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []

  # postprocessors for mass conservation
  [mass_conservation_sum_encl1_encl2]
    type = LinearCombinationPostprocessor
    pp_names = 'concentration_HT_encl_1_inventory concentration_HT_encl_2_inventory concentration_H2_encl_1_inventory concentration_H2_encl_2_inventory concentration_T2_encl_1_inventory concentration_T2_encl_2_inventory'
    pp_coefs = '1            1           1            1            1            1'
    execute_on = 'initial timestep_end'
  []

  # postprocessors for equilibrium constant
  [equilibrium_constant_encl_1]
    type = ParsedPostprocessor
    expression = 'pressure_HT_enclosure_1 / sqrt(pressure_H2_enclosure_1 * pressure_T2_enclosure_1)'
    pp_names = 'pressure_HT_enclosure_1 pressure_H2_enclosure_1 pressure_T2_enclosure_1'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
  [equilibrium_constant_encl_2]
    type = ParsedPostprocessor
    expression = 'pressure_HT_enclosure_2 / sqrt(pressure_H2_enclosure_2 * pressure_T2_enclosure_2)'
    pp_names = 'pressure_HT_enclosure_2 pressure_H2_enclosure_2 pressure_T2_enclosure_2'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
[]

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

[Executioner]
  type = Transient
  end_time = ${simulation_time}
  dtmax = 1e-2
  nl_max_its = 9
  l_max_its = 30
  scheme = 'bdf2'
  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e-6
    optimal_iterations = 7
    iteration_window = 1
    growth_factor = 1.1
    cutback_factor = 0.9
    cutback_factor_at_failure = 0.9
  []
[]

[Outputs]
  file_base = 'ver-1kd_out_k10'
  csv = true
  exodus = true
  execute_on = 'initial timestep_end'
[]

(test/tests/ver-1kc-1/ver-1kc-1.i)

nb_segments_TMAP7 = 20
node_size_TMAP7 = '${units 1.25e-5 m}'
long_total = '${units ${fparse nb_segments_TMAP7 * node_size_TMAP7} m}'
nb_segments_TMAP8 = 1e2
simulation_time = '${units 10 s}'
temperature = '${units 500 K}'
R = '${units 8.31446261815324 J/mol/K}' # ideal gas constant from PhysicalConstants.h
initial_pressure_1 = '${units 1e5 Pa}'
initial_pressure_2 = '${units 1e-10 Pa}'
initial_concentration_1 = '${units ${fparse initial_pressure_1 / (R*temperature)} mol/m^3}'
initial_concentration_2 = '${units ${fparse initial_pressure_2 / (R*temperature)} mol/m^3}'
solubility = '${units ${fparse 10/sqrt(R*temperature)} mol/m^3/Pa^(1/2)}' # Sieverts' law solubility
diffusivity = '${units ${fparse 4.31e-6 * exp(-2818/temperature)} m^2/s}'
n_sorption = 0.5 # Sieverts' Law
unit_scale = 1
unit_scale_neighbor = 1

[Mesh]
  [generated]
    type = GeneratedMeshGenerator
    dim = 1
    nx = ${nb_segments_TMAP8}
    xmax = ${long_total}
  []
  [enclosure_1]
    type = SubdomainBoundingBoxGenerator
    input = generated
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '${fparse 1/3 * long_total} 0 0'
  []
  [enclosure_2]
    type = SubdomainBoundingBoxGenerator
    input = enclosure_1
    block_id = 2
    bottom_left = '${fparse 1/3 * long_total} 0 0'
    top_right = '${fparse long_total} 0 0'
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = enclosure_2
    primary_block = 1
    paired_block = 2
    new_boundary = interface
  []
  [interface2]
    type = SideSetsBetweenSubdomainsGenerator
    input = interface
    primary_block = 2
    paired_block = 1
    new_boundary = interface2
  []
[]

[Variables]
  [concentration_enclosure_1]
    block = 1
    order = FIRST
    family = LAGRANGE
    initial_condition = '${initial_concentration_1}'
  []
  [concentration_enclosure_2]
    block = 2
    order = FIRST
    family = LAGRANGE
    initial_condition = '${initial_concentration_2}'
  []
[]

[AuxVariables]
  [pressure_enclosure_1]
    order = CONSTANT
    family = MONOMIAL
    block = 1
    initial_condition = '${initial_pressure_1}'
  []
  [pressure_enclosure_2]
    order = CONSTANT
    family = MONOMIAL
    block = 2
    initial_condition = '${initial_pressure_2}'
  []
[]

[Kernels]
  [diffusion_enclosure_1]
    type = MatDiffusion
    variable = concentration_enclosure_1
    diffusivity = ${diffusivity}
    block = '1'
  []
  [time_enclosure_1]
    type = TimeDerivative
    variable = concentration_enclosure_1
    block = '1'
  []
  [diffusion_enclosure_2]
    type = MatDiffusion
    variable = concentration_enclosure_2
    diffusivity = ${diffusivity}
    block = '2'
  []
  [time_enclosure_2]
    type = TimeDerivative
    variable = concentration_enclosure_2
    block = '2'
  []
[]

[AuxKernels]
  [pressure_enclosure_1]
    type = ParsedAux
    variable = pressure_enclosure_1
    coupled_variables = 'concentration_enclosure_1'
    expression = '${fparse R*temperature}*concentration_enclosure_1'
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_enclosure_2]
    type = ParsedAux
    variable = pressure_enclosure_2
    coupled_variables = 'concentration_enclosure_2'
    expression = '${fparse R*temperature}*concentration_enclosure_2'
    block = 2
    execute_on = 'initial timestep_end'
  []
[]

[InterfaceKernels]
  [interface_sorption]
    type = InterfaceSorption
    K0 = ${solubility}
    Ea = 0
    n_sorption = ${n_sorption}
    diffusivity = ${diffusivity}
    unit_scale = ${unit_scale}
    unit_scale_neighbor = ${unit_scale_neighbor}
    temperature = ${temperature}
    variable = concentration_enclosure_1
    neighbor_var = concentration_enclosure_2
    sorption_penalty = 1e1
    boundary = interface
  []
[]

[Postprocessors]
  [pressure_enclosure_1]
    type = ElementAverageValue
    variable = pressure_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_enclosure_2]
    type = ElementAverageValue
    variable = pressure_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []
  [concentration_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = concentration_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface2
    variable = concentration_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_ratio]
    type = ParsedPostprocessor
    expression = 'concentration_enclosure_1_at_interface / sqrt(concentration_enclosure_2_at_interface)'
    pp_names = 'concentration_enclosure_1_at_interface concentration_enclosure_2_at_interface'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
  [solubility]
    type = ParsedPostprocessor
    expression = 'concentration_enclosure_1_at_interface / pressure_enclosure_2_at_interface'
    pp_names = 'concentration_enclosure_1_at_interface pressure_enclosure_2_at_interface'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
  [pressure_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = pressure_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [pressure_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface2
    variable = pressure_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_encl_1_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [concentration_encl_2_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []
  [mass_conservation_sum_encl1_encl2]
    type = LinearCombinationPostprocessor
    pp_names = 'concentration_encl_1_inventory concentration_encl_2_inventory'
    pp_coefs = '1            1'
    execute_on = 'initial timestep_end'
  []
[]

[Executioner]
  type = Transient
  end_time = ${simulation_time}
  dtmax = 10
  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-6
  scheme = 'bdf2'
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  nl_max_its = 6
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.01
    optimal_iterations = 5
    iteration_window = 3
    growth_factor = 1.05
    cutback_factor_at_failure = 0.8
  []
[]

[Outputs]
  file_base = 'ver-1kc-1_out_k10'
  csv = true
  exodus = true
  execute_on = 'initial timestep_end'
[]

Description
Input File Usage Example
Input Parameters
Input Files