ADMatInterfaceReactionZr2FeHxPCT

Implements a reaction to establish ReactionRate=k_f*u-k_b*v to compute the surface H concentration in Zr2FeHx 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 zirconium-iron hydride (Zr $var element = document.getElementById("moose-equation-31755ae3-3ebd-494b-93ec-72bb27099079");katex.render("_{2}", element, {displayMode:false,throwOnError:false});$ FeH $var element = document.getElementById("moose-equation-a7385e98-6664-4506-80be-5ddd7bb82cb6");katex.render("_{x}", element, {displayMode:false,throwOnError:false});$ ) is described by pressure-composition-temperature (PCT) data. In TMAP8, the PCT curve can be imposed by an interface kernel which then dictates the material composition. Wherein, the ADMatInterfaceReactionZr2FeHxPCT interface kernel imposes the surface concentration of H in Zr $var element = document.getElementById("moose-equation-03de4087-1379-432a-82ab-3d130638a9d9");katex.render("_{2}", element, {displayMode:false,throwOnError:false});$ FeH $var element = document.getElementById("moose-equation-8956dd74-0ec5-45f0-aab7-85e6c00a5323");katex.render("_{x}", element, {displayMode:false,throwOnError:false});$ based on the input pressure (Pa) and temperature (K). ADMatInterfaceReactionZr2FeHxPCT is related to ADMatInterfaceReaction, but incorporates the Zr $var element = document.getElementById("moose-equation-96f318dc-8b84-4135-a712-69733200cd9b");katex.render("_{2}", element, {displayMode:false,throwOnError:false});$ FeH $var element = document.getElementById("moose-equation-90c6a641-b4dc-4b17-80f0-220fe836a612");katex.render("_{x}", element, {displayMode:false,throwOnError:false});$ PCT curves. At the interface between a solid (main) and a gas (neighbor), it imposes: $(1)var element = document.getElementById("moose-equation-e63258d4-2bec-40c2-b6d2-46e3e2f003a6");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-3c4ab50c-def2-4e6a-bcb4-c8a0181409d9");katex.render("C_s", element, {displayMode:false,throwOnError:false});$ is the surface H concentration in mol/m $var element = document.getElementById("moose-equation-82a8802e-333f-4b81-a03b-76576c6b6e72");katex.render("^3", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-9f473ebb-2fc3-4778-b1c1-82a1e2761c43");katex.render("K_b", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-dfbef657-e62d-4f3c-b088-8a7bea83050f");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-9a49b428-d2ab-4233-a959-08e05383539c");katex.render("f_{at}(T,P)", element, {displayMode:false,throwOnError:false});$ is the composition in atomic fraction of H in Zr $var element = document.getElementById("moose-equation-200c111c-1563-4cbe-ae2c-5cddf5a0c105");katex.render("_{2}", element, {displayMode:false,throwOnError:false});$ FeH $var element = document.getElementById("moose-equation-095ee421-3861-4179-90b1-7626839892c8");katex.render("_{x}", element, {displayMode:false,throwOnError:false});$ given a gas temperature $var element = document.getElementById("moose-equation-54584469-e688-4fa3-9e93-a3a1d4dbdaad");katex.render("T", element, {displayMode:false,throwOnError:false});$ and gas pressure $var element = document.getElementById("moose-equation-5f2f8a85-4dce-40c2-a2af-d8640509bfda");katex.render("P", element, {displayMode:false,throwOnError:false});$ , and $var element = document.getElementById("moose-equation-d5d27705-0950-4471-9b74-59c20dd150b2");katex.render("\\rho", element, {displayMode:false,throwOnError:false});$ is the zirconium-iron molar density in mol/m $var element = document.getElementById("moose-equation-841799fc-ccbd-4e73-bf88-3e17ef52d817");katex.render("^3", element, {displayMode:false,throwOnError:false});$ . Note that the neighbor pressure is given as a molecular concentration (i.e., not an atomic concentration) and is converted to pressure within the interface kernel.

Figure 1 shows the experimental data from Yang et al. (2025) along with the TMAP8 fitts and simulation results.

PCT data for Zr$_{2}$FeH$_{x}$ from [!cite](yang2025potential) with fits implemented in TMAP8 and test results.

Figure 1: PCT data for Zr $var element = document.getElementById("moose-equation-bcc8f4f6-1cac-4f66-8855-9c8f0bcd4a2b");katex.render("_{2}", element, {displayMode:false,throwOnError:false});$ FeH $var element = document.getElementById("moose-equation-ee37c460-5952-403b-a5b4-f07a68668996");katex.render("_{x}", element, {displayMode:false,throwOnError:false});$ from Yang et al. (2025) with fits implemented in TMAP8 and test results.

The pressure isotherm is captured as: $(2)var element = document.getElementById("moose-equation-705480b4-e8bd-461c-a842-ee565e66519a");katex.render(" f_{at}(T,P) = 5.0 - 8.32 \\times 10^{-3}\\left[1 \\times 10^{-3} +\\exp(-2.49 - 7.62\\times 10^{-3} T + (5.63 \\times 10^{-2} + 1.72 \\times 10^{-4} T) (\\log\\left(P - P_{lim}\\right)))\\right]^{-1},", element, {displayMode:true,throwOnError:false});$

with the pressure limit $var element = document.getElementById("moose-equation-8e85136f-c963-4095-9ba1-b77bb9eb40c1");katex.render("P_{lim}", element, {displayMode:false,throwOnError:false});$ set to be a constant value written as $(3)var element = document.getElementById("moose-equation-4d2fc124-fde6-479f-b114-a3bdd4748ddd");katex.render(" P_{lim} = 5\\ \\text{Pa}.", element, {displayMode:true,throwOnError:false});$

The validity of this present fit is between: $(4)var element = document.getElementById("moose-equation-69c121c0-b3df-415a-b121-85e2663afd90");katex.render(" 7 < P\\,\\text{[Pa]} < 5 \\times 10^{5}.", element, {displayMode:true,throwOnError:false});$

Test

(test/tests/Zr2Fe_hydrogen_system/Zr2FeHx_PCT.i) tests the implantation of the Zr $var element = document.getElementById("moose-equation-5edc7700-de46-4caa-ba07-35cbd08772e0");katex.render("_{2}", element, {displayMode:false,throwOnError:false});$ FeH $var element = document.getElementById("moose-equation-16293339-d918-4761-9253-a87591e69533");katex.render("_{x}", element, {displayMode:false,throwOnError:false});$ PCT curves in TMAP8 using constant testing conditions, i.e pressure and temperature. The domain contains two blocks: gas (left) and Zr $var element = document.getElementById("moose-equation-176c7217-e13c-4811-992c-e3d7b8b360e3");katex.render("_{2}", element, {displayMode:false,throwOnError:false});$ FeH $var element = document.getElementById("moose-equation-a0bf2ca2-ce85-4f9f-acb8-529befaed8a1");katex.render("_{x}", element, {displayMode:false,throwOnError:false});$ (right) with an interface between the two blocks. The diffusion is given by Yu et al. (2024) and the surface reaction rate $var element = document.getElementById("moose-equation-3e8fc1f5-95d9-47de-874c-708a7a7736fc");katex.render("K", element, {displayMode:false,throwOnError:false});$ is taken from Yang et al. (2025) ( $var element = document.getElementById("moose-equation-03abfe1e-c23c-4335-9bc8-7b9af55000d4");katex.render("K_f=K_b=K", element, {displayMode:false,throwOnError:false});$ ). Note that the diffusion is for ZrH $var element = document.getElementById("moose-equation-0ba1ac4c-791f-4c7a-a9b5-35af032ebd89");katex.render("_{1.8}", element, {displayMode:false,throwOnError:false});$ because no diffusion data exist for Zr $var element = document.getElementById("moose-equation-399f7048-7d77-4f38-80fe-d77eb9a4ec1d");katex.render("_{2}", element, {displayMode:false,throwOnError:false});$ FeH $var element = document.getElementById("moose-equation-1342d753-7e2a-4101-a38d-76c6c969f8dd");katex.render("_{x}", 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 ADMatInterfaceReactionZr2FeHxPCT to model the steady-state condition for the hydrogen concentration at the surface $var element = document.getElementById("moose-equation-121621be-bc46-40c2-a942-8b3ec6e04463");katex.render("C_s", element, {displayMode:false,throwOnError:false});$ defined by: $(5)var element = document.getElementById("moose-equation-db705bbd-d699-4477-8c2f-d7116e98ac7b");katex.render(" \\frac{d C_s}{dt} = 0 = K (f_{at}(T,P) \\rho - C_s).", element, {displayMode:true,throwOnError:false});$

The results of the high pressure test for ( $var element = document.getElementById("moose-equation-1e159268-30f3-4510-9456-19ceec74c598");katex.render("T", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-99595910-334f-4a44-855c-9dce15c39751");katex.render("P", element, {displayMode:false,throwOnError:false});$ ) =(598.15 K, $var element = document.getElementById("moose-equation-0aa8953c-1378-4ecc-985d-beda45073e8c");katex.render("1 \\times 10^{3}", element, {displayMode:false,throwOnError:false});$ Pa), (623.15 K, $var element = document.getElementById("moose-equation-a2122f6d-bca4-4546-a978-19ff726378f1");katex.render("1 \\times 10^{4}", element, {displayMode:false,throwOnError:false});$ Pa), (648.15 K, $var element = document.getElementById("moose-equation-67800e22-9cec-4b76-8835-a31ca8095264");katex.render("1 \\times 10^{2}", element, {displayMode:false,throwOnError:false});$ Pa), and (648.15 K, $var element = document.getElementById("moose-equation-88f6e5e4-ade3-471e-a69f-c756e76f2f51");katex.render("1 \\times 10^{5}", element, {displayMode:false,throwOnError:false});$ Pa) are shown in Figure 1 and are a good fit.

Overall PCT Testing

(test/tests/Zr2Fe_hydrogen_system/Zr2FeHx_PCT.i) can also test the entire PCT modelling capabilities in TMAP8. The model follows the same structure as the previoius section, but has an arbitrarily high diffusion value to acheive steady-state quickly. The model also includes a linear pressure increase to cover the entire PCT curve.

$var element = document.getElementById("moose-equation-6af9be46-61f5-45a9-baf6-81e5ca712343");katex.render("P = P_{\\text{initial}} + t \\frac{P_{\\text{max}} - P_{\\text{initial}}}{t_{\\text{end}}}.", element, {displayMode:true,throwOnError:false});$

The testing conditions include ( $var element = document.getElementById("moose-equation-c4064bbe-ceea-4f7b-b106-26ce98f1c6ce");katex.render("T", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-eacd180b-f47f-4d80-8708-57d02246bd84");katex.render("P_{initial}", element, {displayMode:false,throwOnError:false});$ ) = (598.15 K, $var element = document.getElementById("moose-equation-88d9a657-9ddd-494a-a2f9-5edbbda0622f");katex.render("7", element, {displayMode:false,throwOnError:false});$ Pa), (623.15 K, $var element = document.getElementById("moose-equation-85aebe3e-d7c7-423c-941f-5d840a7d69f4");katex.render("7", element, {displayMode:false,throwOnError:false});$ Pa), and (648.15 K, $var element = document.getElementById("moose-equation-e718ec80-bc04-4653-ba9b-d1d54760fbd0");katex.render("7", element, {displayMode:false,throwOnError:false});$ Pa). The Figure 1 shows the PCT fit against the experimental data. Evidently, the fit exhibits moderate error, but the deviations remain within an acceptable range for modeling purposes.

Example Input File Syntax

[InterfaceKernels<<<{"href": "../../syntax/InterfaceKernels/index.html"}>>>]
  [interface_reaction_Zr2FeHx_PCT]
    type = ADMatInterfaceReactionZr2FeHxPCT<<<{"description": "Implements a reaction to establish ReactionRate=k_f*u-k_b*v to compute the surface H concentration in Zr2FeHx 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": "ADMatInterfaceReactionZr2FeHxPCT.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_Zr2Fe}
    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_Zr2FeHx'
    backward_rate<<<{"description": "Backward reaction rate coefficient (1/s)."}>>> = 'reaction_rate_surface_Zr2FeHx'
  []
[]

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).

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_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)
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/Zr2Fe_hydrogen_system/Zr2FeHx_PCT.i)

References

Zhiyi Yang, Yuxiao Jia, Yang Liu, Xuezhang Xiao, Tiao Ying, Xingwen Feng, Yan Shi, Changan Chen, Wenhua Luo, and Lixin Chen. A potential hydrogen isotope storage material Zr2Fe: deep exploration on phase transition behaviors and disproportionation mechanism. Energy Materials, 5(1):500011, Jan 2025. doi:10.20517/energymater.2024.83.

BibTeX

@article{yang2025potential,
    author = "Yang, Zhiyi and Jia, Yuxiao and Liu, Yang and Xiao, Xuezhang and Ying, Tiao and Feng, Xingwen and Shi, Yan and Chen, Changan and Luo, Wenhua and Chen, Lixin",
    title = "A potential hydrogen isotope storage material {Zr2Fe}: deep exploration on phase transition behaviors and disproportionation mechanism",
    journal = "Energy Materials",
    volume = "5",
    number = "1",
    pages = "500011",
    year = "2025",
    month = "Jan",
    doi = "10.20517/energymater.2024.83",
    publisher = "OAE Publishing Inc."
}

RIS

TY  - JOUR
AU  - Yang, Zhiyi
AU  - Jia, Yuxiao
AU  - Liu, Yang
AU  - Xiao, Xuezhang
AU  - Ying, Tiao
AU  - Feng, Xingwen
AU  - Shi, Yan
AU  - Chen, Changan
AU  - Luo, Wenhua
AU  - Chen, Lixin
TI  - A potential hydrogen isotope storage material Zr2Fe: deep exploration on phase transition behaviors and disproportionation mechanism
JO  - Energy Materials
PY  - 2025
VL  - 5
IS  - 1
PB  - OAE Publishing Inc.
DO  - 10.20517/energymater.2024.83
SP  - 500011
ER  -

Plain Text

Zhiyi Yang, Yuxiao Jia, Yang Liu, Xuezhang Xiao, Tiao Ying, Xingwen Feng, Yan Shi, Changan Chen, Wenhua Luo, and Lixin Chen. A potential hydrogen isotope storage material Zr2Fe: deep exploration on phase transition behaviors and disproportionation mechanism. Energy Materials, 5(1):500011, Jan 2025. doi:10.20517/energymater.2024.83.

Feifei Yu, Xia Xiang, Xiaotao Zu, and Shuanglin Hu. Hydrogen diffusion in zirconium hydrides from on-the-fly machine learning molecular dynamics. International Journal of Hydrogen Energy, 56:1057–1066, Feb 2024. doi:10.1016/j.ijhydene.2023.12.241.

BibTeX

@article{yu2024hydrogen,
    author = "Yu, Feifei and Xiang, Xia and Zu, Xiaotao and Hu, Shuanglin",
    title = "Hydrogen Diffusion in Zirconium Hydrides from On-the-Fly Machine Learning Molecular Dynamics",
    journal = "International Journal of Hydrogen Energy",
    volume = "56",
    month = "Feb",
    year = "2024",
    pages = "1057--1066",
    doi = "10.1016/j.ijhydene.2023.12.241",
    publisher = "Elsevier"
}

RIS

TY  - JOUR
AU  - Yu, Feifei
AU  - Xiang, Xia
AU  - Zu, Xiaotao
AU  - Hu, Shuanglin
TI  - Hydrogen Diffusion in Zirconium Hydrides from On-the-Fly Machine Learning Molecular Dynamics
JO  - International Journal of Hydrogen Energy
PY  - 2024
VL  - 56
PB  - Elsevier
DO  - 10.1016/j.ijhydene.2023.12.241
SP  - 1057
EP  - 1066
ER  -

Plain Text

Feifei Yu, Xia Xiang, Xiaotao Zu, and Shuanglin Hu. Hydrogen diffusion in zirconium hydrides from on-the-fly machine learning molecular dynamics. International Journal of Hydrogen Energy, 56:1057–1066, Feb 2024. doi:10.1016/j.ijhydene.2023.12.241.

(test/tests/Zr2Fe_hydrogen_system/Zr2FeHx_PCT.i)

###########################################
# Testing PCT Zr2FeHx modelling capability and performance
# This input file works in two independent ways:
# (A) : Evaluate the equilibrium of hydrogen and Zr2Fe at a set Pressure boundary and Temperature
# (B) : Evaluate the equilibrium of hydrogen and Zr2Fe at an increasing Pressure boundary and set Temperature
###########################################

# 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 598.15 K}'
density_Zr2Fe = '${units 29375.60 mol/m^3}'


initial_pressure_H2_enclosure_1 = '${units 7 Pa}'
initial_concentration_H_enclosure_1 = '${units ${fparse 2*initial_pressure_H2_enclosure_1 / (R*temperature)} mol/m^3}'
initial_atomic_fraction= '${fparse 5.0 - 8.320e-03  / ( 1e-03 + exp(-2.4851 - 7.6091e-03 * temperature + (5.6264e-02 + 1.7197e-04 * temperature) * log(max(initial_pressure_H2_enclosure_1 - 5, 1e-10))))}'
initial_concentration_H_enclosure_2 = '${units ${fparse initial_atomic_fraction*density_Zr2Fe} mol/m^3}'

# diffusivity from:
#Zr:Yu, Feifei, et al. “Hydrogen Diffusion in Zirconium Hydrides from On-the-Fly Machine Learning Molecular Dynamics.” International Journal of Hydrogen Energy, vol. 56, Feb. 2024, pp. 1057–66. DOI.org (Crossref), https://doi.org/10.1016/j.ijhydene.2023.12.241

diffusivity_Do = '${units 1.53e-07 m^2/s}'
diffusivity_Ea = '${units 0.61 eV -> J}'

diffusivity_ratio_gas_Zr2FeHx = '${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:
#Yang, Zhiyi, et al. “A Potential Hydrogen Isotope Storage Material Zr2Fe: Deep Exploration on Phase Transition Behaviors and Disproportionation Mechanism.” Energy Materials, vol. 5, no. 1, Jan. 2025. DOI.org (Crossref), https://doi.org/10.20517/energymater.2024.83.
reaction_rate_0 = '${units 4.28e02 1/s}'
reaction_rate_Ea = '${units 0.9108 eV -> J}'

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

# time
simulation_time = '${units 1e9 s}'
dt_max = '${fparse simulation_time/100}'
dt_init = '${units 1e-3 s}'

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

# file base
output_file_base ='Zr2FeHx_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(-${dt_init}*2e-2/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_Zr2FeHx
    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_Zr2FeHx]
    type = DerivativeParsedMaterial
    property_name = diffusivity_Zr2FeHx
    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_Zr2FeHx
    expression = '${diffusivity_ratio_gas_Zr2FeHx}*diffusivity_Zr2FeHx'
    outputs = exodus
  []
  [reaction_rate_surface_Zr2FeHx_1]
    type = ADDerivativeParsedMaterial
    property_name = reaction_rate_surface_Zr2FeHx
    coupled_variables = 'temperature'
    expression = '${reaction_rate_0} * exp(-${fparse reaction_rate_Ea/boltzmann_constant}/temperature)' # 1/s
    block = '1'
  []
  [reaction_rate_surface_Zr2FeHx_2]
    type = ADDerivativeParsedMaterial
    property_name = reaction_rate_surface_Zr2FeHx
    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_Zr2FeHx
    D_neighbor = diffusivity_gas
  []
  [interface_reaction_Zr2FeHx_PCT]
    type = ADMatInterfaceReactionZr2FeHxPCT
    variable = concentration_H_enclosure_2
    neighbor_var = concentration_H_enclosure_1
    neighbor_temperature = temperature
    density = ${density_Zr2Fe}
    boundary = interface_2
    forward_rate = 'reaction_rate_surface_Zr2FeHx'
    backward_rate = 'reaction_rate_surface_Zr2FeHx'
  []
[]

[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_Zr2Fe}'
    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-2
  nl_abs_tol = 4e-15
  scheme = 'bdf2'
  solve_type = 'Newton'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_type -snes_type'
  petsc_options_value = 'lu       mumps                      vinewtonrsls' # This petsc option helps prevent negative concentrations with bounds'
  line_search = 'none'
  automatic_scaling = true
  compute_scaling_once = true
  [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]
    type = Exodus
    execute_on = 'timestep_end'
  []
[]

(test/tests/Zr2Fe_hydrogen_system/Zr2FeHx_PCT.i)

###########################################
# Testing PCT Zr2FeHx modelling capability and performance
# This input file works in two independent ways:
# (A) : Evaluate the equilibrium of hydrogen and Zr2Fe at a set Pressure boundary and Temperature
# (B) : Evaluate the equilibrium of hydrogen and Zr2Fe at an increasing Pressure boundary and set Temperature
###########################################

# 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 598.15 K}'
density_Zr2Fe = '${units 29375.60 mol/m^3}'


initial_pressure_H2_enclosure_1 = '${units 7 Pa}'
initial_concentration_H_enclosure_1 = '${units ${fparse 2*initial_pressure_H2_enclosure_1 / (R*temperature)} mol/m^3}'
initial_atomic_fraction= '${fparse 5.0 - 8.320e-03  / ( 1e-03 + exp(-2.4851 - 7.6091e-03 * temperature + (5.6264e-02 + 1.7197e-04 * temperature) * log(max(initial_pressure_H2_enclosure_1 - 5, 1e-10))))}'
initial_concentration_H_enclosure_2 = '${units ${fparse initial_atomic_fraction*density_Zr2Fe} mol/m^3}'

# diffusivity from:
#Zr:Yu, Feifei, et al. “Hydrogen Diffusion in Zirconium Hydrides from On-the-Fly Machine Learning Molecular Dynamics.” International Journal of Hydrogen Energy, vol. 56, Feb. 2024, pp. 1057–66. DOI.org (Crossref), https://doi.org/10.1016/j.ijhydene.2023.12.241

diffusivity_Do = '${units 1.53e-07 m^2/s}'
diffusivity_Ea = '${units 0.61 eV -> J}'

diffusivity_ratio_gas_Zr2FeHx = '${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:
#Yang, Zhiyi, et al. “A Potential Hydrogen Isotope Storage Material Zr2Fe: Deep Exploration on Phase Transition Behaviors and Disproportionation Mechanism.” Energy Materials, vol. 5, no. 1, Jan. 2025. DOI.org (Crossref), https://doi.org/10.20517/energymater.2024.83.
reaction_rate_0 = '${units 4.28e02 1/s}'
reaction_rate_Ea = '${units 0.9108 eV -> J}'

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

# time
simulation_time = '${units 1e9 s}'
dt_max = '${fparse simulation_time/100}'
dt_init = '${units 1e-3 s}'

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

# file base
output_file_base ='Zr2FeHx_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(-${dt_init}*2e-2/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_Zr2FeHx
    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_Zr2FeHx]
    type = DerivativeParsedMaterial
    property_name = diffusivity_Zr2FeHx
    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_Zr2FeHx
    expression = '${diffusivity_ratio_gas_Zr2FeHx}*diffusivity_Zr2FeHx'
    outputs = exodus
  []
  [reaction_rate_surface_Zr2FeHx_1]
    type = ADDerivativeParsedMaterial
    property_name = reaction_rate_surface_Zr2FeHx
    coupled_variables = 'temperature'
    expression = '${reaction_rate_0} * exp(-${fparse reaction_rate_Ea/boltzmann_constant}/temperature)' # 1/s
    block = '1'
  []
  [reaction_rate_surface_Zr2FeHx_2]
    type = ADDerivativeParsedMaterial
    property_name = reaction_rate_surface_Zr2FeHx
    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_Zr2FeHx
    D_neighbor = diffusivity_gas
  []
  [interface_reaction_Zr2FeHx_PCT]
    type = ADMatInterfaceReactionZr2FeHxPCT
    variable = concentration_H_enclosure_2
    neighbor_var = concentration_H_enclosure_1
    neighbor_temperature = temperature
    density = ${density_Zr2Fe}
    boundary = interface_2
    forward_rate = 'reaction_rate_surface_Zr2FeHx'
    backward_rate = 'reaction_rate_surface_Zr2FeHx'
  []
[]

[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_Zr2Fe}'
    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-2
  nl_abs_tol = 4e-15
  scheme = 'bdf2'
  solve_type = 'Newton'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_type -snes_type'
  petsc_options_value = 'lu       mumps                      vinewtonrsls' # This petsc option helps prevent negative concentrations with bounds'
  line_search = 'none'
  automatic_scaling = true
  compute_scaling_once = true
  [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]
    type = Exodus
    execute_on = 'timestep_end'
  []
[]

(test/tests/Zr2Fe_hydrogen_system/Zr2FeHx_PCT.i)

###########################################
# Testing PCT Zr2FeHx modelling capability and performance
# This input file works in two independent ways:
# (A) : Evaluate the equilibrium of hydrogen and Zr2Fe at a set Pressure boundary and Temperature
# (B) : Evaluate the equilibrium of hydrogen and Zr2Fe at an increasing Pressure boundary and set Temperature
###########################################

# 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 598.15 K}'
density_Zr2Fe = '${units 29375.60 mol/m^3}'


initial_pressure_H2_enclosure_1 = '${units 7 Pa}'
initial_concentration_H_enclosure_1 = '${units ${fparse 2*initial_pressure_H2_enclosure_1 / (R*temperature)} mol/m^3}'
initial_atomic_fraction= '${fparse 5.0 - 8.320e-03  / ( 1e-03 + exp(-2.4851 - 7.6091e-03 * temperature + (5.6264e-02 + 1.7197e-04 * temperature) * log(max(initial_pressure_H2_enclosure_1 - 5, 1e-10))))}'
initial_concentration_H_enclosure_2 = '${units ${fparse initial_atomic_fraction*density_Zr2Fe} mol/m^3}'

# diffusivity from:
#Zr:Yu, Feifei, et al. “Hydrogen Diffusion in Zirconium Hydrides from On-the-Fly Machine Learning Molecular Dynamics.” International Journal of Hydrogen Energy, vol. 56, Feb. 2024, pp. 1057–66. DOI.org (Crossref), https://doi.org/10.1016/j.ijhydene.2023.12.241

diffusivity_Do = '${units 1.53e-07 m^2/s}'
diffusivity_Ea = '${units 0.61 eV -> J}'

diffusivity_ratio_gas_Zr2FeHx = '${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:
#Yang, Zhiyi, et al. “A Potential Hydrogen Isotope Storage Material Zr2Fe: Deep Exploration on Phase Transition Behaviors and Disproportionation Mechanism.” Energy Materials, vol. 5, no. 1, Jan. 2025. DOI.org (Crossref), https://doi.org/10.20517/energymater.2024.83.
reaction_rate_0 = '${units 4.28e02 1/s}'
reaction_rate_Ea = '${units 0.9108 eV -> J}'

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

# time
simulation_time = '${units 1e9 s}'
dt_max = '${fparse simulation_time/100}'
dt_init = '${units 1e-3 s}'

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

# file base
output_file_base ='Zr2FeHx_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(-${dt_init}*2e-2/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_Zr2FeHx
    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_Zr2FeHx]
    type = DerivativeParsedMaterial
    property_name = diffusivity_Zr2FeHx
    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_Zr2FeHx
    expression = '${diffusivity_ratio_gas_Zr2FeHx}*diffusivity_Zr2FeHx'
    outputs = exodus
  []
  [reaction_rate_surface_Zr2FeHx_1]
    type = ADDerivativeParsedMaterial
    property_name = reaction_rate_surface_Zr2FeHx
    coupled_variables = 'temperature'
    expression = '${reaction_rate_0} * exp(-${fparse reaction_rate_Ea/boltzmann_constant}/temperature)' # 1/s
    block = '1'
  []
  [reaction_rate_surface_Zr2FeHx_2]
    type = ADDerivativeParsedMaterial
    property_name = reaction_rate_surface_Zr2FeHx
    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_Zr2FeHx
    D_neighbor = diffusivity_gas
  []
  [interface_reaction_Zr2FeHx_PCT]
    type = ADMatInterfaceReactionZr2FeHxPCT
    variable = concentration_H_enclosure_2
    neighbor_var = concentration_H_enclosure_1
    neighbor_temperature = temperature
    density = ${density_Zr2Fe}
    boundary = interface_2
    forward_rate = 'reaction_rate_surface_Zr2FeHx'
    backward_rate = 'reaction_rate_surface_Zr2FeHx'
  []
[]

[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_Zr2Fe}'
    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-2
  nl_abs_tol = 4e-15
  scheme = 'bdf2'
  solve_type = 'Newton'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_type -snes_type'
  petsc_options_value = 'lu       mumps                      vinewtonrsls' # This petsc option helps prevent negative concentrations with bounds'
  line_search = 'none'
  automatic_scaling = true
  compute_scaling_once = true
  [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]
    type = Exodus
    execute_on = 'timestep_end'
  []
[]

(test/tests/Zr2Fe_hydrogen_system/Zr2FeHx_PCT.i)

###########################################
# Testing PCT Zr2FeHx modelling capability and performance
# This input file works in two independent ways:
# (A) : Evaluate the equilibrium of hydrogen and Zr2Fe at a set Pressure boundary and Temperature
# (B) : Evaluate the equilibrium of hydrogen and Zr2Fe at an increasing Pressure boundary and set Temperature
###########################################

# 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 598.15 K}'
density_Zr2Fe = '${units 29375.60 mol/m^3}'


initial_pressure_H2_enclosure_1 = '${units 7 Pa}'
initial_concentration_H_enclosure_1 = '${units ${fparse 2*initial_pressure_H2_enclosure_1 / (R*temperature)} mol/m^3}'
initial_atomic_fraction= '${fparse 5.0 - 8.320e-03  / ( 1e-03 + exp(-2.4851 - 7.6091e-03 * temperature + (5.6264e-02 + 1.7197e-04 * temperature) * log(max(initial_pressure_H2_enclosure_1 - 5, 1e-10))))}'
initial_concentration_H_enclosure_2 = '${units ${fparse initial_atomic_fraction*density_Zr2Fe} mol/m^3}'

# diffusivity from:
#Zr:Yu, Feifei, et al. “Hydrogen Diffusion in Zirconium Hydrides from On-the-Fly Machine Learning Molecular Dynamics.” International Journal of Hydrogen Energy, vol. 56, Feb. 2024, pp. 1057–66. DOI.org (Crossref), https://doi.org/10.1016/j.ijhydene.2023.12.241

diffusivity_Do = '${units 1.53e-07 m^2/s}'
diffusivity_Ea = '${units 0.61 eV -> J}'

diffusivity_ratio_gas_Zr2FeHx = '${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:
#Yang, Zhiyi, et al. “A Potential Hydrogen Isotope Storage Material Zr2Fe: Deep Exploration on Phase Transition Behaviors and Disproportionation Mechanism.” Energy Materials, vol. 5, no. 1, Jan. 2025. DOI.org (Crossref), https://doi.org/10.20517/energymater.2024.83.
reaction_rate_0 = '${units 4.28e02 1/s}'
reaction_rate_Ea = '${units 0.9108 eV -> J}'

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

# time
simulation_time = '${units 1e9 s}'
dt_max = '${fparse simulation_time/100}'
dt_init = '${units 1e-3 s}'

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

# file base
output_file_base ='Zr2FeHx_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(-${dt_init}*2e-2/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_Zr2FeHx
    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_Zr2FeHx]
    type = DerivativeParsedMaterial
    property_name = diffusivity_Zr2FeHx
    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_Zr2FeHx
    expression = '${diffusivity_ratio_gas_Zr2FeHx}*diffusivity_Zr2FeHx'
    outputs = exodus
  []
  [reaction_rate_surface_Zr2FeHx_1]
    type = ADDerivativeParsedMaterial
    property_name = reaction_rate_surface_Zr2FeHx
    coupled_variables = 'temperature'
    expression = '${reaction_rate_0} * exp(-${fparse reaction_rate_Ea/boltzmann_constant}/temperature)' # 1/s
    block = '1'
  []
  [reaction_rate_surface_Zr2FeHx_2]
    type = ADDerivativeParsedMaterial
    property_name = reaction_rate_surface_Zr2FeHx
    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_Zr2FeHx
    D_neighbor = diffusivity_gas
  []
  [interface_reaction_Zr2FeHx_PCT]
    type = ADMatInterfaceReactionZr2FeHxPCT
    variable = concentration_H_enclosure_2
    neighbor_var = concentration_H_enclosure_1
    neighbor_temperature = temperature
    density = ${density_Zr2Fe}
    boundary = interface_2
    forward_rate = 'reaction_rate_surface_Zr2FeHx'
    backward_rate = 'reaction_rate_surface_Zr2FeHx'
  []
[]

[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_Zr2Fe}'
    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-2
  nl_abs_tol = 4e-15
  scheme = 'bdf2'
  solve_type = 'Newton'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_type -snes_type'
  petsc_options_value = 'lu       mumps                      vinewtonrsls' # This petsc option helps prevent negative concentrations with bounds'
  line_search = 'none'
  automatic_scaling = true
  compute_scaling_once = true
  [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]
    type = Exodus
    execute_on = 'timestep_end'
  []
[]

Description
Test
Example Input File Syntax
Input Parameters
Input Files
References