ADMatInterfaceReactionYHxPCT

Implements a reaction to establish ReactionRate=k_f*u-k_b*v to compute the surface H concentration in YHx from the temperature and partial pressure based on the PCT curves with u the concentration in the solid and v (neighbor) the concentration in the gas in mol/m^3.

Description

The composition of metal hydrides like yttrium hydride (YHx) is described by pressure-composition-temperature (PCT) data. In TMAP8, the PCT curve can be imposed by an interface kernel that then dictates the material composition. Hence, the ADMatInterfaceReactionYHxPCT interface kernel imposes the surface concentration of H in YHx based on the input pressure (Pa) and temperature (K). ADMatInterfaceReactionYHxPCT is related to ADMatInterfaceReaction, but includes the YHx PCT curves. At the interface between a solid (main) and a gas (neighbor), it imposes: $(1)var element = document.getElementById("moose-equation-0f4398a4-443b-4b7e-9d69-b4df10708626");katex.render(" \\frac{d C_s}{dt} = 0 = K_b f_{at}(T,P) \\rho - K_f C_s ,", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-f4c0745f-79e2-4e63-89a2-905deb7dfc5f");katex.render("C_s", element, {displayMode:false,throwOnError:false});$ is the surface H concentration in mol/m $var element = document.getElementById("moose-equation-d80fcee5-b603-40dd-aff4-955fe684cbdb");katex.render("^3", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-ff333863-9497-4ac3-b502-fa60869e3999");katex.render("K_b", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-be7344a7-5381-4cb8-bfe4-212597d5bc1c");katex.render("K_f", element, {displayMode:false,throwOnError:false});$ are the backward and forward surface reaction rate in 1/s, respectively, $var element = document.getElementById("moose-equation-c98d533b-e3b7-4427-8cf5-33873253236a");katex.render("f_{at}(T,P)", element, {displayMode:false,throwOnError:false});$ is the composition in atomic fraction of H in YHx given a gas temperature $var element = document.getElementById("moose-equation-f933a5da-0d6e-4fdb-a46e-8fbb3a927f46");katex.render("T", element, {displayMode:false,throwOnError:false});$ and gas pressure $var element = document.getElementById("moose-equation-8674c137-02d4-4cf9-805a-3e33674676e3");katex.render("P", element, {displayMode:false,throwOnError:false});$ , and $var element = document.getElementById("moose-equation-54973409-3f48-4cd8-aa73-2c8c0480195a");katex.render("\\rho", element, {displayMode:false,throwOnError:false});$ is the yttrium atomic density in mol/m $var element = document.getElementById("moose-equation-7f71aa69-5714-4e6e-b9ca-499b7fea0e1b");katex.render("^3", element, {displayMode:false,throwOnError:false});$ .

Figure 1 shows the data used in this interface kernel. The experimental data originates from Lundin and Blackledge (1962), and the fit is from Matthews et al. (2021).

Figure 1: PCT data for YHx from Lundin and Blackledge (1962) with fits from Matthews et al. (2021) implemented in TMAP8, and test results.

The plateau representing phase transition is captured as (Matthews et al., 2021): $(2)var element = document.getElementById("moose-equation-78b32f0e-f314-4aa1-999c-a02b675a4f29");katex.render(" P_{lim} = \\exp\\left(-26.1+3.88 \\times 10^{-2} T - 9.7 \\times 10^{-6} T^2 \\right),", element, {displayMode:true,throwOnError:false});$ with $var element = document.getElementById("moose-equation-601bbafc-6da9-4ccb-9d44-3162bd1458a9");katex.render("P_{lim}", element, {displayMode:false,throwOnError:false});$ being the hydrogen partial pressure limit delineating the plateau in Pa and $var element = document.getElementById("moose-equation-4b5e7270-888a-4c11-9294-0a943e4e550d");katex.render("T", element, {displayMode:false,throwOnError:false});$ being the temperature in K. This fit is shown in Figure 2.

Figure 2: Fit phase transition region as pressure as a function of temperature.

The model currently only captures the high pressure region of the data, i.e., for pressure and temperatures above the plateau. In this region, the atomic fraction $var element = document.getElementById("moose-equation-c7a221fd-20d5-4c8b-a087-bc22a73caa20");katex.render("f_{at}", element, {displayMode:false,throwOnError:false});$ is determined as (Matthews et al., 2021): $(3)var element = document.getElementById("moose-equation-0692a7d9-646d-4c92-ae49-91fbc6764f93");katex.render(" f_{at}(T,P) = 2-\\left[1+\\exp(21.6 -0.0225 T + (-0.0445 + 7.18 \\times 10^{-4} T) (\\log\\left(P - P_{lim}(T)\\right)))\\right]^{-1},", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-0b66bf27-aba9-4499-88ec-35cb6ba3959d");katex.render("P", element, {displayMode:false,throwOnError:false});$ is the hydrogen partial pressure in Pa. This fit is plotted in Figure 1.

The "silence_warnings" option can be used to dictate how TMAP8 reacts when the pressure gets out of bounds. If silence_warnings = false, which is the default behavior, then TMAP8 will print a warning stating that the pressure and temperature are outside the bounds of the atomic fraction correlation. If silence_warnings = true, then TMAP8 will let the simulation continue without issuing any warnings.

Test

(test/tests/yttrium_hydrogen_system/YHx_PCT.i) tests the implantation of the YHx PCT curves in TMAP8. The domain contains two blocks: gas (left) and YHx (right) with an interface between the two blocks. The diffusion is given by (Majer et al., 2002) and the surface reaction rate $var element = document.getElementById("moose-equation-8e14ebe9-6d20-4b49-9d3e-ed178ede3ddf");katex.render("K", element, {displayMode:false,throwOnError:false});$ is taken from (Fisher and Tanase, 1984) ( $var element = document.getElementById("moose-equation-82656fee-818d-4d48-a1f2-293d6439ddc0");katex.render("K_f=K_b=K", element, {displayMode:false,throwOnError:false});$ ). To model the interface, the input file employs the InterfaceDiffusion object to model the flux of hydrogen at the surface, and ADMatInterfaceReactionYHxPCT to model the steady-state condition for the hydrogen concentration at the surface $var element = document.getElementById("moose-equation-faf81b58-604b-4f0e-a158-2d22c87acb1b");katex.render("C_s", element, {displayMode:false,throwOnError:false});$ defined by: $(4)var element = document.getElementById("moose-equation-2e944d88-c471-43a0-99af-14fcd29de79d");katex.render(" \\frac{d C_s}{dt} = 0 = K (f_{at}(T,P) \\rho - C_s),", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-c879621a-21a2-4aa6-8a96-3eccb1a4bd4c");katex.render("\\rho", element, {displayMode:false,throwOnError:false});$ is the yttrium atomic density.

The results of the test for ( $var element = document.getElementById("moose-equation-57172ce1-33fb-4f10-98d5-73872ba514a1");katex.render("T", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-2e5e841e-fb64-422a-9eb5-c2144de001c1");katex.render("P", element, {displayMode:false,throwOnError:false});$ ) = (1173.15 K, $var element = document.getElementById("moose-equation-bbcdbaab-0482-4f71-bd6c-e51d5af7a1b6");katex.render("1 \\times 10^{3}", element, {displayMode:false,throwOnError:false});$ Pa), (1173.15 K, $var element = document.getElementById("moose-equation-8563ed70-d04a-4f88-b0b0-19cc0f58651e");katex.render("1 \\times 10^{4}", element, {displayMode:false,throwOnError:false});$ Pa), (1173.15 K, $var element = document.getElementById("moose-equation-ed5f961c-7572-4050-8ca1-06fef2887435");katex.render("5 \\times 10^{4}", element, {displayMode:false,throwOnError:false});$ Pa), and(1473.15 K, $var element = document.getElementById("moose-equation-7bdefaef-d349-463d-9fa1-3cb4e1a6d237");katex.render("5 \\times 10^{4}", element, {displayMode:false,throwOnError:false});$ Pa) are shown in Figure 1 and show good agreement with Eq. (3).

Example Input File Syntax

[InterfaceKernels<<<{"href": "../../syntax/InterfaceKernels/index.html"}>>>]
  [interface_reaction_YHx_PCT]
    type = ADMatInterfaceReactionYHxPCT<<<{"description": "Implements a reaction to establish ReactionRate=k_f*u-k_b*v to compute the surface H concentration in YHx from the temperature and partial pressure based on the PCT curves with u the concentration in the solid and v (neighbor) the concentration in the gas in mol/m^3.", "href": "ADMatInterfaceReactionYHxPCT.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = concentration_H_enclosure_2
    neighbor_var<<<{"description": "The variable on the other side of the interface."}>>> = concentration_H_enclosure_1
    neighbor_temperature<<<{"description": "The variable on the other side of the interface for temperature (K)."}>>> = temperature
    density<<<{"description": "Density of the solid in (mol/m^3)."}>>> = ${density_Y}
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = interface_2
    forward_rate<<<{"description": "Forward reaction rate coefficient (1/s)."}>>> = 'reaction_rate_surface_YHx'
    backward_rate<<<{"description": "Backward reaction rate coefficient (1/s)."}>>> = 'reaction_rate_surface_YHx'
  []
[]

Input Parameters

neighbor_temperatureThe variable on the other side of the interface for temperature (K).
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The variable on the other side of the interface for temperature (K).
neighbor_varThe variable on the other side of the interface.
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The variable on the other side of the interface.
variableThe name of the variable that this residual object operates on
C++ Type:NonlinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The name of the variable that this residual object operates on

Required Parameters

backward_ratekbBackward reaction rate coefficient (1/s).
Default:kb
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Backward reaction rate coefficient (1/s).
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
densitydensityDensity of the solid in (mol/m^3).
Default:density
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Density of the solid in (mol/m^3).
forward_ratekfForward reaction rate coefficient (1/s).
Default:kf
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Forward reaction rate coefficient (1/s).
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)
silence_warningsFalseWhether to silence correlation out of bound warnings
Default:False
C++ Type:bool
Controllable:No
Description:Whether to silence correlation out of bound warnings

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/yttrium_hydrogen_system/YHx_PCT.i)

References

P. W. Fisher and M. Tanase. Diffusivities of hydrogen in yttrium and yttrium alloys. Journal of Nuclear Materials, 123(1):1536–1540, 1984. URL: https://www.sciencedirect.com/science/article/pii/0022311584902976, doi:https://doi.org/10.1016/0022-3115(84)90297-6.

@article{FISHER19841536,
    author = "Fisher, P. W. and Tanase, M.",
    title = "Diffusivities of hydrogen in yttrium and yttrium alloys",
    journal = "Journal of Nuclear Materials",
    volume = "123",
    number = "1",
    pages = "1536-1540",
    year = "1984",
    issn = "0022-3115",
    doi = "https://doi.org/10.1016/0022-3115(84)90297-6",
    url = "https://www.sciencedirect.com/science/article/pii/0022311584902976"
}

C. E. Lundin and J. P. Blackledge. Pressure‐temperature‐composition relationships of the yttrium‐hydrogen system. Journal of The Electrochemical Society, 109(9):838, sep 1962. URL: https://dx.doi.org/10.1149/1.2425565, doi:10.1149/1.2425565.

@article{Lundin_1962,
    author = "Lundin, C. E. and Blackledge, J. P.",
    doi = "10.1149/1.2425565",
    url = "https://dx.doi.org/10.1149/1.2425565",
    year = "1962",
    month = "sep",
    publisher = "The Electrochemical Society, Inc.",
    volume = "109",
    number = "9",
    pages = "838",
    title = "Pressure‐Temperature‐Composition Relationships of the Yttrium‐Hydrogen System",
    journal = "Journal of The Electrochemical Society"
}

G. Majer, J. Gottwald, D. T. Peterson, and R. G. Barnes. Model-independent measurements of hydrogen diffusivity in the yttrium dihydrides. Journal of Alloys and Compounds, 330-332:438–442, 2002. Proceedings of the International Symposium on Metal-Hydrogen Systems, Fundamentals and Applications (MH2000). URL: https://www.sciencedirect.com/science/article/pii/S0925838801014529, doi:https://doi.org/10.1016/S0925-8388(01)01452-9.

@article{MAJER2002438,
    author = "Majer, G. and Gottwald, J. and Peterson, D. T. and Barnes, R. G.",
    title = "Model-independent measurements of hydrogen diffusivity in the yttrium dihydrides",
    journal = "Journal of Alloys and Compounds",
    volume = "330-332",
    pages = "438-442",
    year = "2002",
    note = "Proceedings of the International Symposium on Metal-Hydrogen Systems, Fundamentals and Applications (MH2000)",
    issn = "0925-8388",
    doi = "https://doi.org/10.1016/S0925-8388(01)01452-9",
    url = "https://www.sciencedirect.com/science/article/pii/S0925838801014529"
}

C. Matthews, Shivprasad A. P., and Cooper M. W. D. Metal Hydride Simulations Using SWIFT. Technical Report LA-UR-21-27538, Los Alamos National Laboratory, Los Alamos, NM (United States), 2021.

@techreport{Matthews2021SWIFT,
    author = "Matthews, C. and P., Shivprasad A. and D., Cooper M. W.",
    title = "{Metal Hydride Simulations Using SWIFT}",
    year = "2021",
    institution = "Los Alamos National Laboratory, Los Alamos, NM (United States)",
    number = "LA-UR-21-27538"
}

(test/tests/yttrium_hydrogen_system/YHx_PCT.i)

# Physical constants
R = '${units 8.31446261815324 J/mol/K}' # ideal gas constant from PhysicalConstants.h
N_a = '${units 6.02214076e23 1/mol}' # Avogadro's number from PhysicalConstants.h
boltzmann_constant = '${units 1.380649e-23 J/K}' # Boltzmann constant from PhysicalConstants.h

# Simulation conditions and material properties
temperature = '${units 1200 K}'
density_Y = '${units 48605 mol/m^3}'
initial_pressure_H2_enclosure_1 = '${units 1e4 Pa}'
initial_concentration_H_enclosure_1 = '${units ${fparse 2*initial_pressure_H2_enclosure_1 / (R*temperature)} mol/m^3}'
initial_atomic_fraction = 1.8 # (-)
initial_concentration_H_enclosure_2 = '${units ${fparse initial_atomic_fraction*density_Y} mol/m^3}'

# diffusivity from Majer et al., Journal of Alloys and Compounds 330-332 (2002) 438-442.
diffusivity_Do = '${units 1.e-8 m^2/s}'
diffusivity_Ea = '${units 0.38 eV -> J}'
diffusivity_ratio_gas_YHx = '${fparse initial_concentration_H_enclosure_2 / initial_concentration_H_enclosure_1 * 10}' # this ratio is large and helps InterfaceDiffusion due to the ratio of concentrations
# Surface reaction rate from P. W. Fisher, M. Tanase, Journal of Nuclear Materials 122-123 (1984) 1536-1540.
reaction_rate_0 = '${units 4.95e5 1/s}'
reaction_rate_Ea = '${units 1.52 eV -> J}'

# Domain size and mesh parameters
domain_length = '${units 1 m}'
num_nodes = 8

# time
simulation_time = '${units 1e9 s}'
dt_max = '${fparse simulation_time/100}'
dt_init = '${units 1e-3 s}'
tau_constant_BC = '${fparse dt_init*2e-2}' # the smaller, the faster the up-ramp for the pressure BC

# convergence parameters
lower_value_threshold_concentration_enclosure_1 = -1e-20
lower_value_threshold_concentration_enclosure_2 = -1e-20

# file base
output_file_base = 'YHx_PCT_out'

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

[Variables]
  [concentration_H_enclosure_1]
    block = 1
    initial_condition = '${initial_concentration_H_enclosure_1}'
  []
  [concentration_H_enclosure_2]
    block = 2
    initial_condition = '${initial_concentration_H_enclosure_2}'
  []
[]

[Bounds]
  # To prevent negative concentrations
  [concentration_H_enclosure_1_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_concentration_H_enclosure_1
    bounded_variable = concentration_H_enclosure_1
    bound_type = lower
    bound_value = ${lower_value_threshold_concentration_enclosure_1}
  []
  [concentration_H_enclosure_2_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_concentration_H_enclosure_2
    bounded_variable = concentration_H_enclosure_2
    bound_type = lower
    bound_value = ${lower_value_threshold_concentration_enclosure_2}
  []
[]

[BCs]
  [concentration_H_enclosure_1_fixed]
    type = FunctionDirichletBC
    variable = concentration_H_enclosure_1
    boundary = left
    function = 'function_BC_concentration_H_enclosure_1'
  []
[]

[Functions]
  [function_BC_concentration_H_enclosure_1]
    type = ParsedFunction
    expression = 'exp(-${tau_constant_BC}/t)* ${initial_concentration_H_enclosure_1}'
  []
[]

[Kernels]
  # Diffusion equation for H in enclosure 1
  [H_time_derivative_enclosure_1]
    type = TimeDerivative
    variable = concentration_H_enclosure_1
    block = '1'
  []
  [H_diffusion_enclosure_1]
    type = MatDiffusion
    variable = concentration_H_enclosure_1
    diffusivity = diffusivity_gas
    block = '1'
  []
  # Diffusion equation for H in enclosure 2
  [H_time_derivative_enclosure_2]
    type = TimeDerivative
    variable = concentration_H_enclosure_2
    block = '2'
  []
  [H_diffusion_enclosure_2]
    type = MatDiffusion
    variable = concentration_H_enclosure_2
    diffusivity = diffusivity_YHx
    block = '2'
  []
[]

[AuxVariables]
  [temperature]
    initial_condition = '${temperature}'
  []
  [pressure_H2_enclosure_1]
    initial_condition = '${initial_pressure_H2_enclosure_1}'
  []
  [bounds_dummy_concentration_H_enclosure_1]
    order = FIRST
    family = LAGRANGE
  []
  [bounds_dummy_concentration_H_enclosure_2]
    order = FIRST
    family = LAGRANGE
  []
[]

[AuxKernels]
  [pressure_H2_enclosure_1]
    type = ParsedAux
    variable = pressure_H2_enclosure_1
    coupled_variables = 'concentration_H_enclosure_1 temperature'
    expression = '${R} * temperature * concentration_H_enclosure_1/2'
    block = 1
    execute_on = 'initial timestep_end'
  []
[]

[Materials]
  [diffusivity_YHx]
    type = DerivativeParsedMaterial
    property_name = diffusivity_YHx
    coupled_variables = 'temperature'
    expression = '${diffusivity_Do} * exp(-${fparse diffusivity_Ea*N_a/R}/temperature)'
    outputs = exodus
  []
  [diffusivity_gas]
    type = DerivativeParsedMaterial
    property_name = diffusivity_gas
    material_property_names = diffusivity_YHx
    expression = '${diffusivity_ratio_gas_YHx}*diffusivity_YHx'
    outputs = exodus
  []
  [reaction_rate_surface_YHx_1]
    type = ADDerivativeParsedMaterial
    property_name = reaction_rate_surface_YHx
    coupled_variables = 'temperature'
    expression = '${reaction_rate_0} * exp(-${fparse reaction_rate_Ea/boltzmann_constant}/temperature)' # 1/s
    block = '1'
  []
  [reaction_rate_surface_YHx_2]
    type = ADDerivativeParsedMaterial
    property_name = reaction_rate_surface_YHx
    coupled_variables = 'temperature'
    expression = '${reaction_rate_0} * exp(-${fparse reaction_rate_Ea/boltzmann_constant}/temperature)' # 1/s
    block = '2'
  []
[]

[InterfaceKernels]
  [interface_diffusion]
    type = InterfaceDiffusion
    variable = concentration_H_enclosure_2
    neighbor_var = concentration_H_enclosure_1
    boundary = interface_2
    D = diffusivity_YHx
    D_neighbor = diffusivity_gas
  []
  [interface_reaction_YHx_PCT]
    type = ADMatInterfaceReactionYHxPCT
    variable = concentration_H_enclosure_2
    neighbor_var = concentration_H_enclosure_1
    neighbor_temperature = temperature
    density = ${density_Y}
    boundary = interface_2
    forward_rate = 'reaction_rate_surface_YHx'
    backward_rate = 'reaction_rate_surface_YHx'
  []
[]

[Postprocessors]
  [temperature]
    type = ElementAverageValue
    variable = temperature
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_H2_enclosure_1]
    type = ElementAverageValue
    variable = pressure_H2_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [concentration_H_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = concentration_H_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_H_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface_2
    variable = concentration_H_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [atomic_fraction_H_enclosure_2_at_interface]
    type = ParsedPostprocessor
    pp_names = 'concentration_H_enclosure_2_at_interface'
    expression = 'concentration_H_enclosure_2_at_interface / ${density_Y}'
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_ratio_H2]
    type = ParsedPostprocessor
    expression = 'concentration_H_enclosure_1_at_interface / sqrt(concentration_H_enclosure_2_at_interface)'
    pp_names = 'concentration_H_enclosure_1_at_interface concentration_H_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'
  []
  [concentration_H_encl_1_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_H_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [concentration_H_encl_2_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_H_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []

  # postprocessors for mass conservation
  [mass_conservation_sum_encl1_encl2]
    type = LinearCombinationPostprocessor
    pp_names = 'concentration_H_encl_1_inventory concentration_H_encl_2_inventory'
    pp_coefs = '1                                1'
    execute_on = 'initial timestep_end'
  []
[]

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

[Executioner]
  type = Transient
  end_time = ${simulation_time}
  dtmax = ${dt_max}
  nl_max_its = 16
  l_max_its = 30
  nl_rel_tol = 1e-7
  nl_abs_tol = 4e-12
  scheme = 'bdf2'
  solve_type = 'Newton'
  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'
  line_search = 'none'
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = ${dt_init}
    optimal_iterations = 13
    iteration_window = 1
    growth_factor = 1.2
    cutback_factor = 0.9
    cutback_factor_at_failure = 0.9
  []
[]

[Outputs]
  file_base = ${output_file_base}
  csv = true
  exodus = true
[]

(test/tests/yttrium_hydrogen_system/YHx_PCT.i)

# Physical constants
R = '${units 8.31446261815324 J/mol/K}' # ideal gas constant from PhysicalConstants.h
N_a = '${units 6.02214076e23 1/mol}' # Avogadro's number from PhysicalConstants.h
boltzmann_constant = '${units 1.380649e-23 J/K}' # Boltzmann constant from PhysicalConstants.h

# Simulation conditions and material properties
temperature = '${units 1200 K}'
density_Y = '${units 48605 mol/m^3}'
initial_pressure_H2_enclosure_1 = '${units 1e4 Pa}'
initial_concentration_H_enclosure_1 = '${units ${fparse 2*initial_pressure_H2_enclosure_1 / (R*temperature)} mol/m^3}'
initial_atomic_fraction = 1.8 # (-)
initial_concentration_H_enclosure_2 = '${units ${fparse initial_atomic_fraction*density_Y} mol/m^3}'

# diffusivity from Majer et al., Journal of Alloys and Compounds 330-332 (2002) 438-442.
diffusivity_Do = '${units 1.e-8 m^2/s}'
diffusivity_Ea = '${units 0.38 eV -> J}'
diffusivity_ratio_gas_YHx = '${fparse initial_concentration_H_enclosure_2 / initial_concentration_H_enclosure_1 * 10}' # this ratio is large and helps InterfaceDiffusion due to the ratio of concentrations
# Surface reaction rate from P. W. Fisher, M. Tanase, Journal of Nuclear Materials 122-123 (1984) 1536-1540.
reaction_rate_0 = '${units 4.95e5 1/s}'
reaction_rate_Ea = '${units 1.52 eV -> J}'

# Domain size and mesh parameters
domain_length = '${units 1 m}'
num_nodes = 8

# time
simulation_time = '${units 1e9 s}'
dt_max = '${fparse simulation_time/100}'
dt_init = '${units 1e-3 s}'
tau_constant_BC = '${fparse dt_init*2e-2}' # the smaller, the faster the up-ramp for the pressure BC

# convergence parameters
lower_value_threshold_concentration_enclosure_1 = -1e-20
lower_value_threshold_concentration_enclosure_2 = -1e-20

# file base
output_file_base = 'YHx_PCT_out'

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

[Variables]
  [concentration_H_enclosure_1]
    block = 1
    initial_condition = '${initial_concentration_H_enclosure_1}'
  []
  [concentration_H_enclosure_2]
    block = 2
    initial_condition = '${initial_concentration_H_enclosure_2}'
  []
[]

[Bounds]
  # To prevent negative concentrations
  [concentration_H_enclosure_1_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_concentration_H_enclosure_1
    bounded_variable = concentration_H_enclosure_1
    bound_type = lower
    bound_value = ${lower_value_threshold_concentration_enclosure_1}
  []
  [concentration_H_enclosure_2_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_concentration_H_enclosure_2
    bounded_variable = concentration_H_enclosure_2
    bound_type = lower
    bound_value = ${lower_value_threshold_concentration_enclosure_2}
  []
[]

[BCs]
  [concentration_H_enclosure_1_fixed]
    type = FunctionDirichletBC
    variable = concentration_H_enclosure_1
    boundary = left
    function = 'function_BC_concentration_H_enclosure_1'
  []
[]

[Functions]
  [function_BC_concentration_H_enclosure_1]
    type = ParsedFunction
    expression = 'exp(-${tau_constant_BC}/t)* ${initial_concentration_H_enclosure_1}'
  []
[]

[Kernels]
  # Diffusion equation for H in enclosure 1
  [H_time_derivative_enclosure_1]
    type = TimeDerivative
    variable = concentration_H_enclosure_1
    block = '1'
  []
  [H_diffusion_enclosure_1]
    type = MatDiffusion
    variable = concentration_H_enclosure_1
    diffusivity = diffusivity_gas
    block = '1'
  []
  # Diffusion equation for H in enclosure 2
  [H_time_derivative_enclosure_2]
    type = TimeDerivative
    variable = concentration_H_enclosure_2
    block = '2'
  []
  [H_diffusion_enclosure_2]
    type = MatDiffusion
    variable = concentration_H_enclosure_2
    diffusivity = diffusivity_YHx
    block = '2'
  []
[]

[AuxVariables]
  [temperature]
    initial_condition = '${temperature}'
  []
  [pressure_H2_enclosure_1]
    initial_condition = '${initial_pressure_H2_enclosure_1}'
  []
  [bounds_dummy_concentration_H_enclosure_1]
    order = FIRST
    family = LAGRANGE
  []
  [bounds_dummy_concentration_H_enclosure_2]
    order = FIRST
    family = LAGRANGE
  []
[]

[AuxKernels]
  [pressure_H2_enclosure_1]
    type = ParsedAux
    variable = pressure_H2_enclosure_1
    coupled_variables = 'concentration_H_enclosure_1 temperature'
    expression = '${R} * temperature * concentration_H_enclosure_1/2'
    block = 1
    execute_on = 'initial timestep_end'
  []
[]

[Materials]
  [diffusivity_YHx]
    type = DerivativeParsedMaterial
    property_name = diffusivity_YHx
    coupled_variables = 'temperature'
    expression = '${diffusivity_Do} * exp(-${fparse diffusivity_Ea*N_a/R}/temperature)'
    outputs = exodus
  []
  [diffusivity_gas]
    type = DerivativeParsedMaterial
    property_name = diffusivity_gas
    material_property_names = diffusivity_YHx
    expression = '${diffusivity_ratio_gas_YHx}*diffusivity_YHx'
    outputs = exodus
  []
  [reaction_rate_surface_YHx_1]
    type = ADDerivativeParsedMaterial
    property_name = reaction_rate_surface_YHx
    coupled_variables = 'temperature'
    expression = '${reaction_rate_0} * exp(-${fparse reaction_rate_Ea/boltzmann_constant}/temperature)' # 1/s
    block = '1'
  []
  [reaction_rate_surface_YHx_2]
    type = ADDerivativeParsedMaterial
    property_name = reaction_rate_surface_YHx
    coupled_variables = 'temperature'
    expression = '${reaction_rate_0} * exp(-${fparse reaction_rate_Ea/boltzmann_constant}/temperature)' # 1/s
    block = '2'
  []
[]

[InterfaceKernels]
  [interface_diffusion]
    type = InterfaceDiffusion
    variable = concentration_H_enclosure_2
    neighbor_var = concentration_H_enclosure_1
    boundary = interface_2
    D = diffusivity_YHx
    D_neighbor = diffusivity_gas
  []
  [interface_reaction_YHx_PCT]
    type = ADMatInterfaceReactionYHxPCT
    variable = concentration_H_enclosure_2
    neighbor_var = concentration_H_enclosure_1
    neighbor_temperature = temperature
    density = ${density_Y}
    boundary = interface_2
    forward_rate = 'reaction_rate_surface_YHx'
    backward_rate = 'reaction_rate_surface_YHx'
  []
[]

[Postprocessors]
  [temperature]
    type = ElementAverageValue
    variable = temperature
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_H2_enclosure_1]
    type = ElementAverageValue
    variable = pressure_H2_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [concentration_H_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = concentration_H_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_H_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface_2
    variable = concentration_H_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [atomic_fraction_H_enclosure_2_at_interface]
    type = ParsedPostprocessor
    pp_names = 'concentration_H_enclosure_2_at_interface'
    expression = 'concentration_H_enclosure_2_at_interface / ${density_Y}'
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_ratio_H2]
    type = ParsedPostprocessor
    expression = 'concentration_H_enclosure_1_at_interface / sqrt(concentration_H_enclosure_2_at_interface)'
    pp_names = 'concentration_H_enclosure_1_at_interface concentration_H_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'
  []
  [concentration_H_encl_1_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_H_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [concentration_H_encl_2_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_H_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []

  # postprocessors for mass conservation
  [mass_conservation_sum_encl1_encl2]
    type = LinearCombinationPostprocessor
    pp_names = 'concentration_H_encl_1_inventory concentration_H_encl_2_inventory'
    pp_coefs = '1                                1'
    execute_on = 'initial timestep_end'
  []
[]

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

[Executioner]
  type = Transient
  end_time = ${simulation_time}
  dtmax = ${dt_max}
  nl_max_its = 16
  l_max_its = 30
  nl_rel_tol = 1e-7
  nl_abs_tol = 4e-12
  scheme = 'bdf2'
  solve_type = 'Newton'
  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'
  line_search = 'none'
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = ${dt_init}
    optimal_iterations = 13
    iteration_window = 1
    growth_factor = 1.2
    cutback_factor = 0.9
    cutback_factor_at_failure = 0.9
  []
[]

[Outputs]
  file_base = ${output_file_base}
  csv = true
  exodus = true
[]

(test/tests/yttrium_hydrogen_system/YHx_PCT.i)

# Physical constants
R = '${units 8.31446261815324 J/mol/K}' # ideal gas constant from PhysicalConstants.h
N_a = '${units 6.02214076e23 1/mol}' # Avogadro's number from PhysicalConstants.h
boltzmann_constant = '${units 1.380649e-23 J/K}' # Boltzmann constant from PhysicalConstants.h

# Simulation conditions and material properties
temperature = '${units 1200 K}'
density_Y = '${units 48605 mol/m^3}'
initial_pressure_H2_enclosure_1 = '${units 1e4 Pa}'
initial_concentration_H_enclosure_1 = '${units ${fparse 2*initial_pressure_H2_enclosure_1 / (R*temperature)} mol/m^3}'
initial_atomic_fraction = 1.8 # (-)
initial_concentration_H_enclosure_2 = '${units ${fparse initial_atomic_fraction*density_Y} mol/m^3}'

# diffusivity from Majer et al., Journal of Alloys and Compounds 330-332 (2002) 438-442.
diffusivity_Do = '${units 1.e-8 m^2/s}'
diffusivity_Ea = '${units 0.38 eV -> J}'
diffusivity_ratio_gas_YHx = '${fparse initial_concentration_H_enclosure_2 / initial_concentration_H_enclosure_1 * 10}' # this ratio is large and helps InterfaceDiffusion due to the ratio of concentrations
# Surface reaction rate from P. W. Fisher, M. Tanase, Journal of Nuclear Materials 122-123 (1984) 1536-1540.
reaction_rate_0 = '${units 4.95e5 1/s}'
reaction_rate_Ea = '${units 1.52 eV -> J}'

# Domain size and mesh parameters
domain_length = '${units 1 m}'
num_nodes = 8

# time
simulation_time = '${units 1e9 s}'
dt_max = '${fparse simulation_time/100}'
dt_init = '${units 1e-3 s}'
tau_constant_BC = '${fparse dt_init*2e-2}' # the smaller, the faster the up-ramp for the pressure BC

# convergence parameters
lower_value_threshold_concentration_enclosure_1 = -1e-20
lower_value_threshold_concentration_enclosure_2 = -1e-20

# file base
output_file_base = 'YHx_PCT_out'

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

[Variables]
  [concentration_H_enclosure_1]
    block = 1
    initial_condition = '${initial_concentration_H_enclosure_1}'
  []
  [concentration_H_enclosure_2]
    block = 2
    initial_condition = '${initial_concentration_H_enclosure_2}'
  []
[]

[Bounds]
  # To prevent negative concentrations
  [concentration_H_enclosure_1_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_concentration_H_enclosure_1
    bounded_variable = concentration_H_enclosure_1
    bound_type = lower
    bound_value = ${lower_value_threshold_concentration_enclosure_1}
  []
  [concentration_H_enclosure_2_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_concentration_H_enclosure_2
    bounded_variable = concentration_H_enclosure_2
    bound_type = lower
    bound_value = ${lower_value_threshold_concentration_enclosure_2}
  []
[]

[BCs]
  [concentration_H_enclosure_1_fixed]
    type = FunctionDirichletBC
    variable = concentration_H_enclosure_1
    boundary = left
    function = 'function_BC_concentration_H_enclosure_1'
  []
[]

[Functions]
  [function_BC_concentration_H_enclosure_1]
    type = ParsedFunction
    expression = 'exp(-${tau_constant_BC}/t)* ${initial_concentration_H_enclosure_1}'
  []
[]

[Kernels]
  # Diffusion equation for H in enclosure 1
  [H_time_derivative_enclosure_1]
    type = TimeDerivative
    variable = concentration_H_enclosure_1
    block = '1'
  []
  [H_diffusion_enclosure_1]
    type = MatDiffusion
    variable = concentration_H_enclosure_1
    diffusivity = diffusivity_gas
    block = '1'
  []
  # Diffusion equation for H in enclosure 2
  [H_time_derivative_enclosure_2]
    type = TimeDerivative
    variable = concentration_H_enclosure_2
    block = '2'
  []
  [H_diffusion_enclosure_2]
    type = MatDiffusion
    variable = concentration_H_enclosure_2
    diffusivity = diffusivity_YHx
    block = '2'
  []
[]

[AuxVariables]
  [temperature]
    initial_condition = '${temperature}'
  []
  [pressure_H2_enclosure_1]
    initial_condition = '${initial_pressure_H2_enclosure_1}'
  []
  [bounds_dummy_concentration_H_enclosure_1]
    order = FIRST
    family = LAGRANGE
  []
  [bounds_dummy_concentration_H_enclosure_2]
    order = FIRST
    family = LAGRANGE
  []
[]

[AuxKernels]
  [pressure_H2_enclosure_1]
    type = ParsedAux
    variable = pressure_H2_enclosure_1
    coupled_variables = 'concentration_H_enclosure_1 temperature'
    expression = '${R} * temperature * concentration_H_enclosure_1/2'
    block = 1
    execute_on = 'initial timestep_end'
  []
[]

[Materials]
  [diffusivity_YHx]
    type = DerivativeParsedMaterial
    property_name = diffusivity_YHx
    coupled_variables = 'temperature'
    expression = '${diffusivity_Do} * exp(-${fparse diffusivity_Ea*N_a/R}/temperature)'
    outputs = exodus
  []
  [diffusivity_gas]
    type = DerivativeParsedMaterial
    property_name = diffusivity_gas
    material_property_names = diffusivity_YHx
    expression = '${diffusivity_ratio_gas_YHx}*diffusivity_YHx'
    outputs = exodus
  []
  [reaction_rate_surface_YHx_1]
    type = ADDerivativeParsedMaterial
    property_name = reaction_rate_surface_YHx
    coupled_variables = 'temperature'
    expression = '${reaction_rate_0} * exp(-${fparse reaction_rate_Ea/boltzmann_constant}/temperature)' # 1/s
    block = '1'
  []
  [reaction_rate_surface_YHx_2]
    type = ADDerivativeParsedMaterial
    property_name = reaction_rate_surface_YHx
    coupled_variables = 'temperature'
    expression = '${reaction_rate_0} * exp(-${fparse reaction_rate_Ea/boltzmann_constant}/temperature)' # 1/s
    block = '2'
  []
[]

[InterfaceKernels]
  [interface_diffusion]
    type = InterfaceDiffusion
    variable = concentration_H_enclosure_2
    neighbor_var = concentration_H_enclosure_1
    boundary = interface_2
    D = diffusivity_YHx
    D_neighbor = diffusivity_gas
  []
  [interface_reaction_YHx_PCT]
    type = ADMatInterfaceReactionYHxPCT
    variable = concentration_H_enclosure_2
    neighbor_var = concentration_H_enclosure_1
    neighbor_temperature = temperature
    density = ${density_Y}
    boundary = interface_2
    forward_rate = 'reaction_rate_surface_YHx'
    backward_rate = 'reaction_rate_surface_YHx'
  []
[]

[Postprocessors]
  [temperature]
    type = ElementAverageValue
    variable = temperature
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_H2_enclosure_1]
    type = ElementAverageValue
    variable = pressure_H2_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [concentration_H_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = concentration_H_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_H_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface_2
    variable = concentration_H_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [atomic_fraction_H_enclosure_2_at_interface]
    type = ParsedPostprocessor
    pp_names = 'concentration_H_enclosure_2_at_interface'
    expression = 'concentration_H_enclosure_2_at_interface / ${density_Y}'
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_ratio_H2]
    type = ParsedPostprocessor
    expression = 'concentration_H_enclosure_1_at_interface / sqrt(concentration_H_enclosure_2_at_interface)'
    pp_names = 'concentration_H_enclosure_1_at_interface concentration_H_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'
  []
  [concentration_H_encl_1_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_H_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [concentration_H_encl_2_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_H_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []

  # postprocessors for mass conservation
  [mass_conservation_sum_encl1_encl2]
    type = LinearCombinationPostprocessor
    pp_names = 'concentration_H_encl_1_inventory concentration_H_encl_2_inventory'
    pp_coefs = '1                                1'
    execute_on = 'initial timestep_end'
  []
[]

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

[Executioner]
  type = Transient
  end_time = ${simulation_time}
  dtmax = ${dt_max}
  nl_max_its = 16
  l_max_its = 30
  nl_rel_tol = 1e-7
  nl_abs_tol = 4e-12
  scheme = 'bdf2'
  solve_type = 'Newton'
  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'
  line_search = 'none'
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = ${dt_init}
    optimal_iterations = 13
    iteration_window = 1
    growth_factor = 1.2
    cutback_factor = 0.9
    cutback_factor_at_failure = 0.9
  []
[]

[Outputs]
  file_base = ${output_file_base}
  csv = true
  exodus = true
[]

Description
Test
Example Input File Syntax
Input Parameters
Input Files
References