IFBAHeProduction

Computes the helium gas production as a result of an IFBA layer applied to the surface of a fuel rod.

Description

IFBAHeProduction computes the Helium gas production as a result of an IFBA layer applied to the surface of a fuel rod.

An integral fuel burnable absorber (IFBA) is used for optimizing fuel assembly reactivity and power distribution in a core. The IFBA is usually applied as a thin layer of ZrB $var element = document.getElementById("moose-equation-69829a77-0ea0-4d36-8fc1-a7736454d561");katex.render("_2", element, {displayMode:false,throwOnError:false});$ over some length of a fuel rod. The boron-10 isotope in the IFBA material absorbs a neutron and results in a lithium and helium atom according to the following reaction: $(1)var element = document.getElementById("moose-equation-192a1ec8-477b-40a9-a7d9-a770b9c1d3d8");katex.render(" {_5^{10}}B + {_0^1}n \\rightarrow {_3^7}Li + {_2^4}He", element, {displayMode:true,throwOnError:false});$

In addition, the IFBA layer is depleted very quickly and is typically used up in the first $var element = document.getElementById("moose-equation-a33c5ceb-7b56-4ec0-a27d-9d248b0144c1");katex.render("\\sim", element, {displayMode:false,throwOnError:false});$ 3 percent of burnup or $var element = document.getElementById("moose-equation-e0769dad-152f-46cc-b18a-a3caaf7d0537");katex.render("\\sim", element, {displayMode:false,throwOnError:false});$ 18 months of exposure.

Helium Production

Since the IFBA layer is normally on the order of a few microns thick, the helium atoms generated are assumed to be released immediately into the plenum.

Two models for the helium gas production (i.e., boron-10 depletion) have been implemented in BISON. The first model uses an equation based on burnup (Lee et al., 2012) and the second is the model used in FRAPCON (K. J. Geelhood et al., 2015).

Burnup Based Model

This burnup based equation for boron-10 depletion was generated using a DeCART depletion calculation (Lee et al., 2012) for boron in an IFBA rod. The effects of boron concentration and U-235 enrichment were studied. In addition, an improved approximation for the time dependency of U-235 number density decrease and fissile plutonium production resulted in the following relationship for the number density of boron-10: $(2)var element = document.getElementById("moose-equation-2f2f9fe5-b678-4b69-ac07-0e12fa4006a1");katex.render(" N_{{^{10}}B}(t) = N_{{^{10}}B}(0) exp\\bigg({-Bu\\over {a w_{25} + (b w_{25} + c) Bu}}\\bigg)", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-29ffc3ba-2fb8-4189-8de1-4a2eb6e8fd25");katex.render("N_{{^{10}}B}(t)", element, {displayMode:false,throwOnError:false});$ is the boron-10 number density at time $var element = document.getElementById("moose-equation-b854911f-ee79-470c-aab5-e5769a958a8c");katex.render("t", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-4a07cd80-106e-435a-aedf-3ab7c0d4d932");katex.render("N_{{^{10}}B}(0)", element, {displayMode:false,throwOnError:false});$ is the initial boron-10 number density, $var element = document.getElementById("moose-equation-81f8b1b9-1aeb-4144-9e28-1fbad58025bc");katex.render("w_{25}", element, {displayMode:false,throwOnError:false});$ is the U-235 enrichment, and $var element = document.getElementById("moose-equation-a95fd92d-694b-43d4-b829-598aaca9046c");katex.render("Bu", element, {displayMode:false,throwOnError:false});$ is the burnup. The values of the parameters are given in Table 1.

Table 1: Parameters for the Evolution of Boron-10 Density with Burnup, Eq. (2) (Lee et al., 2012)

Parameter	$var element = document.getElementById("moose-equation-b412ad02-a8ce-4f84-9d5e-8dc2aa3b8c84");katex.render("a", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-55b58b61-c9d4-4581-ba7e-5a16f88d0e06");katex.render("b", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-b400a460-1e7a-4d2d-8406-ba1d0613ea5a");katex.render("c", element, {displayMode:false,throwOnError:false});$
Value	1.59389	-0.00773	0.01051

FRAPCON Model

In FRAPCON the helium production rate is derived using an MCNP calculation for boron depletion. The empirical fit to this calculation is $(3)var element = document.getElementById("moose-equation-b4040b0d-b24f-4c27-afd3-e18320e3ca28");katex.render(" \\dot{P}_{He} = \\left(A_1 f_{IFBA} + A_2\\right) C^2_{B_{10}} + \\left(B_1 f_{IFBA} + B_2\\right) C_{B_{10}}", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-726534fe-c47a-484c-9591-b304548a1777");katex.render("\\dot{P}_{He}", element, {displayMode:false,throwOnError:false});$ is the helium production rate (number of atoms/cm $var element = document.getElementById("moose-equation-df2cff74-d65b-4607-b7bc-34f141c708c0");katex.render("^3", element, {displayMode:false,throwOnError:false});$ -s), $var element = document.getElementById("moose-equation-5754c46a-6a72-47a0-9af9-b95d787108e3");katex.render("f_{IFBA}", element, {displayMode:false,throwOnError:false});$ is the percent of IFBA rods in a reactor core (this value is limited to the range of 10 percent to 50 percent), and $var element = document.getElementById("moose-equation-edb60b9a-e48c-4879-80f3-f9204ecd6168");katex.render("C_{B_{10}}", element, {displayMode:false,throwOnError:false});$ is the Boron-10 enrichment (percent) (restricted to 0 to 90 percent). The values of the parameters are given in Table 2.

Table 2: Parameters for the Empirical Fit to the MCNP Calculation for He Production Rate, Eq. (3) (K. J. Geelhood et al., 2015)

Parameter	$var element = document.getElementById("moose-equation-8d115404-60bc-4aef-a4e7-aa3af6685a35");katex.render("A_1", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-c08092a0-a630-46e3-9e27-6ada51a2129a");katex.render("A_2", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-9d1fbb1a-31a4-4347-a5be-792efd083db7");katex.render("B_1", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-5359e4f1-54af-4797-9055-d4736ec86eb2");katex.render("B_2", element, {displayMode:false,throwOnError:false});$
Value	$var element = document.getElementById("moose-equation-612bad2e-21d1-4037-b2b1-c9edc0c9c502");katex.render("6.23309\\times 10^{-9}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-157b9159-bee6-4c54-b1b7-5fb345c08bdc");katex.render("7.02006\\times 10^{-7}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-2835de84-55ce-4522-a3a8-26e3fb07dc9f");katex.render("-1.35676\\times 10^{-7}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-7f5f9ba0-38a5-4c52-979a-a63430a0c1c9");katex.render("3.1506\\times 10^{-5}", element, {displayMode:false,throwOnError:false});$

The implementation of this empirical relationship in the FRAPCON code (K. J. Geelhood et al., 2015) uses the following form to calculate the helium production rate $(4)var element = document.getElementById("moose-equation-3f7b7585-8c03-4c63-8ac9-398cffeab480");katex.render(" \\dot{P}_{He} = \\left[(A_1 f_{IFBA} + A_2) C^2_{B_{10}} + (B_1 f_{IFBA} + B_2) C_{B_{10}}\\right] {q^{\\prime}\\over 5.64} \\times {\\rho_{ZrB_2}\\over 4.53}~\\textrm{in (atoms/$\\mathrm {cm^3}$-s)}", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-143dda2e-7a18-450e-84c6-da99373ff718");katex.render("q^{\\prime}", element, {displayMode:false,throwOnError:false});$ is the linear power in kW/ft, $var element = document.getElementById("moose-equation-0b8b7eeb-9aac-4f44-b31d-a361f46bcf8b");katex.render("\\rho_{ZrB_2}", element, {displayMode:false,throwOnError:false});$ is the density of ZrB $var element = document.getElementById("moose-equation-a2e33593-3322-4358-8e87-4e97a3d03d1e");katex.render("_2", element, {displayMode:false,throwOnError:false});$ (g/cm $var element = document.getElementById("moose-equation-fa586441-74c5-46f3-859d-40d3f28334fc");katex.render("^3", element, {displayMode:false,throwOnError:false});$ ), and the coefficients are given in Table 3.

Table 3: Parameters for the Numerical Implementation of the He Production Rate, Eq. (4) (K. J. Geelhood et al., 2015)

Parameter	$var element = document.getElementById("moose-equation-7e6ef9f8-7696-4977-b1ee-3b667cdf7103");katex.render("A_1", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-e408ed06-a6ba-45c3-9a2a-2d9d68aa026c");katex.render("A_2", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-1b1fe721-dcf5-4ea3-9b0d-b9e907893b25");katex.render("B_1", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-2d58c4b1-7f95-40f2-a6d7-88d9366755a1");katex.render("B_2", element, {displayMode:false,throwOnError:false});$
Value	$var element = document.getElementById("moose-equation-b4d0769b-9e79-4f0e-bab5-d5fb3b89531f");katex.render("-9.66127\\times 10^8", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-45c24f01-58df-41b0-808f-50e0f6d872ec");katex.render("-1.088109\\times 10^{11}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-8192c13c-6dae-473c-9c57-bc3a6c4ef037");katex.render("-2.10296\\times 10^{10}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-2a4566c9-7419-44e4-9bbf-039e77ea3d17");katex.render("4.8843\\times 10^{13}", element, {displayMode:false,throwOnError:false});$

Then the boron depletion rate is a function of $var element = document.getElementById("moose-equation-eae68db5-b7c3-475e-9fb1-9e81f731aa16");katex.render("-\\dot{P_{He}}", element, {displayMode:false,throwOnError:false});$ , given as $(5)var element = document.getElementById("moose-equation-08503547-39d0-4feb-88b3-4b8ba8c17709");katex.render(" {C_{B_t}\\over 100} = {C_{B_{t-\\Delta t}}\\over 100} - \\dot{P}_{He} \\Delta t/ \\left(5.84\\times 10^{22} \\rho_{ZrB_2} / \\rho_{90\\%TD}\\right)", element, {displayMode:true,throwOnError:false});$

Input Parameter Conventions

The parameters listed in Table 4 are used to calculate the initial number density of boron-10 atoms and therefore the limiting value for the moles of He gas generated by the IFBA material. Note that the units used by BISON differ from the units generally used in industry; the user must convert the input arguments to the appropriate BISON units.

Table 4: Input Parameters for IFBA Postprocessor

Parameter	Typical Units	BISON Units
ZrB $var element = document.getElementById("moose-equation-57e75833-4c8c-4c12-af64-4db4ef36f5dc");katex.render("_{2}", element, {displayMode:false,throwOnError:false});$ Loading	mg/in	kg/m
B-10 Loading	mg/in	kg/m
IFBA Length	cm	m
B-10 Enrichment	%	fraction
ZrB $var element = document.getElementById("moose-equation-aeae7911-2b4e-4d37-8eb5-14476e8a09e7");katex.render("_{2}", element, {displayMode:false,throwOnError:false});$ Density	% TD	kg/m $var element = document.getElementById("moose-equation-f8588cef-a1d2-4f74-bd90-641814fac834");katex.render("^3", element, {displayMode:false,throwOnError:false});$
or
ZrB $var element = document.getElementById("moose-equation-e14a9b81-ed00-4e07-a685-2e8684e4b9c5");katex.render("_{2}", element, {displayMode:false,throwOnError:false});$ Thickness	$var element = document.getElementById("moose-equation-3e0bdeb6-d6e8-4783-9869-4be429cf492e");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ m	m
Fuel Outer Radius	cm	m

The model specific parameters for the two equations are listed in Table 5.

Table 5: Model Specific Input Parameters for IFBA

Model	Parameter	Typical Units	BISON Units
Burnup Based Model (Eq. (2))	U-235 Enrichment	%	fraction
FRAPCON Model (Eq. (3))	IFBA Rod %	%	fraction

In addition, the burnup model Eq. (2) requires a postprocessor to provide the average burnup in the fuel and the FRAPCON model Eq. (3) requires a postprocessor to provide the rod average linear power.

Example Input Syntax

[Postprocessors<<<{"href": "../../syntax/Postprocessors/index.html"}>>>]
  [he_prod]
    type = IFBAHeProduction<<<{"description": "Computes the helium gas production as a result of an IFBA layer applied to the surface of a fuel rod.", "href": "IFBAHeProduction.html"}>>>
    zrb2_load<<<{"description": "Loading of ZrB2 on pellets (kg/m)"}>>> = 1.181e-4
    ifba_len<<<{"description": "Axial length of IFBA layer on pellets (m)"}>>> = 1.0e-2
    b10_enrich<<<{"description": "B-10 enrichment (fraction)"}>>> = 0.50
    zrb2_rel_dens<<<{"description": "Relative density of ZrB2 material (fraction)"}>>> = 0.7
    model<<<{"description": "He generation model: Choices are burnup or frapcon. Default is burnup."}>>> = burnup
    u235_enrich<<<{"description": "Initial U-235 enrichment (fraction)"}>>> = 0.045
    burnup<<<{"description": "The name of the postprocessor containing the rod average burnup"}>>> = average_burnup
  []
[]

Because this input file is using the burnup model Eq. (2), the average burnup in the fuel post processor is also required. Note that the name of the average burnup postprocessor is used as the argument to the burnup input parameter in the IFBAHeProduction post processor.

[Postprocessors<<<{"href": "../../syntax/Postprocessors/index.html"}>>>]
  [average_burnup]
    type = ElementAverageValue<<<{"description": "Computes the volumetric average of a variable", "href": "ElementAverageValue.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = '3'
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = burnup
  []
[]

Input Parameters

b10_enrichB-10 enrichment (fraction)
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:B-10 enrichment (fraction)
ifba_lenAxial length of IFBA layer on pellets (m)
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Axial length of IFBA layer on pellets (m)

Required Parameters

b10_load-1Loading of B-10 on pellets (kg/m)
Default:-1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Loading of B-10 on pellets (kg/m)
burnupThe name of the postprocessor containing the rod average burnup
C++ Type:PostprocessorName
Unit:(no unit assumed)
Controllable:No
Description:The name of the postprocessor containing the rod average burnup
debug0Debugging output flag: 0 = no output (default), 1 = debugging output to console
Default:0
C++ Type:int
Controllable:No
Description:Debugging output flag: 0 = no output (default), 1 = debugging output to console
fuel_out_rad-1Fuel pellet outer radius (m)
Default:-1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Fuel pellet outer radius (m)
ifba_rod_pct-1Percentage of fuel rods containing IFBA liners (fraction)
Default:-1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Percentage of fuel rods containing IFBA liners (fraction)
modelburnupHe generation model: Choices are burnup or frapcon. Default is burnup.
Default:burnup
C++ Type:MooseEnum
Options:burnup, frapcon
Controllable:No
Description:He generation model: Choices are burnup or frapcon. Default is burnup.
rod_ave_lin_powThe power function
C++ Type:FunctionName
Unit:(no unit assumed)
Controllable:No
Description:The power function
u235_enrich-1Initial U-235 enrichment (fraction)
Default:-1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Initial U-235 enrichment (fraction)
zrb2_load-1Loading of ZrB2 on pellets (kg/m)
Default:-1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Loading of ZrB2 on pellets (kg/m)
zrb2_rel_dens-1Relative density of ZrB2 material (fraction)
Default:-1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Relative density of ZrB2 material (fraction)
zrb2_thick-1Thickness of ZrB2 layer on pellets (m)
Default:-1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Thickness of ZrB2 layer on pellets (m)

Optional Parameters

allow_duplicate_execution_on_initialFalseIn the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).
Default:False
C++ Type:bool
Controllable:No
Description:In the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).
execute_onTIMESTEP_ENDThe list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html.
Default:TIMESTEP_END
C++ Type:ExecFlagEnum
Options:XFEM_MARK, NONE, INITIAL, LINEAR, NONLINEAR_CONVERGENCE, NONLINEAR, POSTCHECK, TIMESTEP_END, TIMESTEP_BEGIN, MULTIAPP_FIXED_POINT_END, MULTIAPP_FIXED_POINT_BEGIN, FINAL, CUSTOM, TRANSFER
Controllable:No
Description:The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html.
execution_order_group0Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.
Default:0
C++ Type:int
Controllable:No
Description:Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.
force_postauxFalseForces the UserObject to be executed in POSTAUX
Default:False
C++ Type:bool
Controllable:No
Description:Forces the UserObject to be executed in POSTAUX
force_preauxFalseForces the UserObject to be executed in PREAUX
Default:False
C++ Type:bool
Controllable:No
Description:Forces the UserObject to be executed in PREAUX
force_preicFalseForces the UserObject to be executed in PREIC during initial setup
Default:False
C++ Type:bool
Controllable:No
Description:Forces the UserObject to be executed in PREIC during initial setup

Execution Scheduling 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.
outputsVector of output names where you would like to restrict the output of variables(s) associated with this object
C++ Type:std::vector<OutputName>
Controllable:No
Description:Vector of output names where you would like to restrict the output of variables(s) associated with this object
use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters

prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.

Material Property Retrieval Parameters

Input Files

(assessment/LWR/validation/LOCA_Studsvik/analysis/rod_196/Studsvik_196_part1.i)
(assessment/LWR/validation/LOCA_Studsvik/analysis/rod_196/Studsvik_196_part2_1p5d_fr_ffrd.i)
(assessment/LWR/validation/LOCA_Studsvik/analysis/rod_196/Studsvik_196_part2.i)
(assessment/LWR/validation/LOCA_Studsvik/analysis/rod_196/Studsvik_196_part1_1p5d_fr_ffrd.i)
(test/tests/ifba_he_production/fill_gas_xenon_w_ifba.i)
(test/tests/ifba_he_production/ifba_examp_template.i)
(test/tests/ifba_he_production/ifba_only_template.i)

References

W. G. Luscher K. J. Geelhood, P. A. Raynaud, and I. E. Porter. FRAPCON-4.0: A Computer Code for the Calculation of Steady-State, Thermal-Mechanical Behavior of Oxide Fuel Rods for High Burnup. Technical Report PNNL-19418 Vol.1 Rev. 2, Pacific Northwest National Laboratory, 2015.

@TECHREPORT{frap2015,
    author = "K. J. Geelhood, W. G. Luscher and Raynaud, P. A. and Porter, I. E.",
    title = "{FRAPCON}-4.0: {A Computer Code for the Calculation of Steady-State, Thermal-Mechanical Behavior of Oxide Fuel Rods for High Burnup}",
    year = "2015",
    number = "PNNL-19418 Vol.1 Rev. 2",
    institution = "Pacific Northwest National Laboratory"
}

C. Lee, K.H. Lee, J.Y. Cho, S.Y. Park, and Y.S. Yang. A study on the boron depletion in 16x16 boron mixed fuel rods. In Transactions of the Korean Nuclear Society Autumn Meeting. Gyeongju, Korea, 2012.

@inproceedings{lee2010,
    author = "Lee, C. and Lee, K.H. and Cho, J.Y. and Park, S.Y. and Yang, Y.S.",
    title = "A Study on the boron depletion in 16x16 boron mixed fuel rods",
    booktitle = "Transactions of the Korean Nuclear Society Autumn Meeting",
    year = "2012",
    address = "Gyeongju, Korea"
}

(test/tests/ifba_he_production/fill_gas_xenon_w_ifba.i)

#
# 2-D RZ One Pellet Test - IFBA using Xenon as fill gas
#
# This test is of a single pellet with cladding and a specified initial
# pressure of Xe fill gas. In addition, an IFBA layer is added which will
# generate He gas to be added to the plenum. The postprocessor interior_temp
# should be the same as the pure Xe test case initially and as the He gas is
# added to the plenum from the IFBA, the interior_temp value should approach
# the He fill gas test case (both in the doc subdirectory).
#
# This model demonstrates that the gas conductance for the plenum is being
# updated for the He gas generated by the IFBA layer.
#

initial_fuel_density = 10431.0 #95% TD (TD = 10980)

[GlobalParams]
  density = ${initial_fuel_density}
  order = SECOND
  family = LAGRANGE
  energy_per_fission = 3.2e-11 # J/fission (205 Mev)
  displacements = 'disp_x disp_y'
  temperature = temp
[]

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

[Mesh]
  coord_type = RZ
  [smeared_pellet_mesh]
    type = FuelPinMeshGenerator
    clad_mesh_density = customize
    pellet_mesh_density = customize
    ny_p = 1
    nx_p = 1
    nx_c = 1
    ny_cu = 1
    ny_c = 1
    ny_cl = 1
    clad_thickness = 5.6e-4
    pellet_outer_radius = 0.0041
    pellet_height = 0.01
    pellet_quantity = 1
    clad_bot_gap_height = 1e-3
    bottom_clad_height = 2.24e-3
    top_clad_height = 2.24e-3
    clad_gap_width = 8e-5
    plenum_fuel_ratio = 0.150
    elem_type = QUAD8
  []
  partitioner = centroid
  centroid_partitioner_direction = y
  patch_size = 5
[]

[Variables]
  [temp]
    initial_condition = 298
  []
[]

[AuxVariables]
  [fission_rate]
    block = '3'
  []
  [burnup]
    block = '3'
  []
  [grain_radius]
    block = '3'
    initial_condition = 5e-6
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 1.0e4 1.0e8'
    y = '0  1.0  1.0'
    scale_factor = 20e3 # 20 kW/m peak power.
  []
  [coolant_pressure_ramp]
    type = PiecewiseLinear
    x = '0 10000'
    y = '0 1'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [fuel]
    block = 3
    strain = FINITE
    incremental = true
    add_variables = true
    decomposition_method = EigenSolution
    extra_vector_tags = 'ref'
    eigenstrain_names = fuel_thermal_strain
  []
  [clad]
    block = 1
    strain = FINITE
    incremental = true
    add_variables = true
    decomposition_method = EigenSolution
    extra_vector_tags = 'ref'
    eigenstrain_names = clad_thermal_strain
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_source_fuel]
    type = NeutronHeatSource
    variable = temp
    block = '3'
    fission_rate = fission_rate
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [fissionrate]
    type = FissionRateGeneral
    fission_rate_formulation = GENERIC
    variable = fission_rate
    block = '3'
    value = 5.3548e+14
    fission_rate_function = power_history
  []
  [burnup]
    type = BurnupAux
    variable = burnup
    block = '3'
    fission_rate = fission_rate
    molecular_weight = 0.270
  []
  [grain_radius]
    type = GrainRadiusAux
    block = '3'
    variable = grain_radius
    temperature = temp
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    penalty = 1e+14 #1e7
    model = frictionless
    tangential_tolerance = 5e-4
    normal_smoothing_distance = 0.1
    normalize_penalty = true
  []
[]

[ThermalContact]
  [pellet_clad_thermal]
    type = GasGapHeatTransfer
    variable = temp
    primary = 5
    secondary = 10
    gas_released = 'fis_gas_released he_prod'
    initial_moles = initial_moles
    jump_distance_model = LANNING
    layer_thickness = layer_thickness
    plenum_pressure = plenum_pressure
    contact_pressure = contact_pressure
    initial_gas_types = Xe
    initial_fractions = 1
    released_gas_types = 'Kr Xe;
                          He'
    released_fractions = '0.153 0.847;
                          1'
    roughness_coef = 3.2
    roughness_secondary = 1e-6
    roughness_primary = 2e-6
    emissivity_primary = 0.8
    emissivity_secondary = 0.8
    quadrature = true
    normal_smoothing_distance = 0.1
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 1003
    value = 0.0
  []
  [no_y_clad_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 1
    value = 0.0
  []
  [no_y_fuel_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 20
    value = 0.0
  []
  [no_x_fuel]
    type = DirichletBC
    variable = disp_x
    boundary = 1005
    value = 0.0
  []
  [Clad_Temp]
    type = DirichletBC
    variable = temp
    boundary = '2'
    value = 580.0
  []
  [Pressure]
    [coolantPressure]
      boundary = '2'
      factor = 15.5e6
      function = coolant_pressure_ramp
    []
  []
  [PlenumPressure]
    [plenumPressure]
      boundary = 9
      initial_pressure = 0.50e6
      startup_time = 0.0
      material_input = 'fis_gas_released he_prod'
      output_initial_moles = initial_moles
      temperature = interior_temp
      volume = gas_volume
      output = plenum_pressure
    []
  []
[]

[Materials]
  [fuel_thermal]
    type = UO2Thermal
    block = '3'
    temperature = temp
    burnup = burnup
    thermal_conductivity_model = NFIR
  []
  [fuel_elasticity_tensor]
    type = UO2ElasticityTensor
    block = 3
  []
  [fuel_elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 3
  []
  [fuel_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = 3
    thermal_expansion_coeff = 10.0e-6
    stress_free_temperature = 298
    eigenstrain_name = fuel_thermal_strain
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = '3'
    temperature = temp
    fission_rate = fission_rate
    grain_radius = grain_radius
    gbs_model = true
    burnup = burnup
    diff_coeff_option = TURNBULL_D1_D2
  []
  [clad_thermal]
    type = HeatConductionMaterial
    block = 1
    thermal_conductivity = 16.0
    specific_heat = 330.0
  []
  [fclad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = 1
    youngs_modulus = 7.5e10
    poissons_ratio = 0.3
  []
  [clad_elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 1
  []
  [clad_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = 1
    thermal_expansion_coeff = 5.0e-6
    stress_free_temperature = 298
    eigenstrain_name = clad_thermal_strain
  []
  [clad_density]
    type = StrainAdjustedDensity
    block = 1
    strain_free_density = 6551.0
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = '3'
    strain_free_density = ${initial_fuel_density}
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

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

  line_search = 'none'

  l_max_its = 25
  nl_max_its = 40
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-8

  dtmax = 1.0e6
  dtmin = 1.0
  end_time = 5.3e7 # 1.7 years (~3% burnup)

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e3
    optimal_iterations = 30
    iteration_window = 4
    time_t = '0   1e4 1e8'
    time_dt = '1e4 1e5 1e6'
    timestep_limiting_function = power_history
    force_step_every_function_point = true
  []
  [Quadrature]
    order = fifth
    side_order = seventh
  []
  verbose = true
[]

[Postprocessors]
  [clad_inner_vol]
    type = InternalVolume
    boundary = 7
    execute_on = 'initial linear'
  []
  [pellet_volume]
    type = InternalVolume
    boundary = 8
    execute_on = 'initial linear'
  []
  [gas_volume]
    type = InternalVolume
    boundary = 9
    execute_on = 'initial linear'
  []
  [interior_temp]
    type = SideAverageValue
    boundary = 9 # cladding interior and pellet exterior
    variable = temp
    execute_on = 'initial linear'
  []
  [fis_gas_produced] # fission gas produced (moles)
    type = ElementIntegralFisGasGeneratedSifgrs
    block = '3'
  []
  [fis_gas_released]
    type = ElementIntegralFisGasReleasedSifgrs
    block = '3'
  []
  [fis_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    block = '3'
  []
  [fis_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    block = '3'
  []
  [power_history]
    type = FunctionValuePostprocessor
    function = power_history
  []
  [dt]
    type = TimestepSize
  []
  [residual]
    type = Residual
  []
  [nl_its]
    type = NumNonlinearIterations
  []
  [lin_its]
    type = NumLinearIterations
  []
  [average_burnup]
    type = ElementAverageValue
    block = '3'
    variable = burnup
  []
  [burnup]
    type = ElementAverageValue
    block = '3'
    variable = burnup
  []
  [average_fissionrate]
    type = ElementAverageValue
    block = '3'
    variable = fission_rate
  []
  [rod_total_power]
    type = ElementIntegralPower
    variable = temp
    fission_rate = fission_rate
    block = '3'
  []
  [rod_input_power]
    type = FunctionValuePostprocessor
    function = power_history
    scale_factor = 0.01 # change: length of fuel stack in meters (1 pellet height)
  []
  [he_prod]
    type = IFBAHeProduction
    zrb2_load = 1.181e-4
    ifba_len = 1.0e-2
    b10_enrich = 0.50
    zrb2_rel_dens = 0.7
    model = burnup
    u235_enrich = 0.045
    burnup = average_burnup
  []
[]

[Outputs]
  time_step_interval =  1
  exodus = false
  [console]
    type = Console
    solve_log = true
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    show = 'average_burnup burnup he_prod interior_temp plenum_pressure'
    file_base = fill_gas_xenon_w_ifba_check
  []
  [out]
    type = CSV
    delimiter = ' '
  []
[]

(test/tests/ifba_he_production/fill_gas_xenon_w_ifba.i)

#
# 2-D RZ One Pellet Test - IFBA using Xenon as fill gas
#
# This test is of a single pellet with cladding and a specified initial
# pressure of Xe fill gas. In addition, an IFBA layer is added which will
# generate He gas to be added to the plenum. The postprocessor interior_temp
# should be the same as the pure Xe test case initially and as the He gas is
# added to the plenum from the IFBA, the interior_temp value should approach
# the He fill gas test case (both in the doc subdirectory).
#
# This model demonstrates that the gas conductance for the plenum is being
# updated for the He gas generated by the IFBA layer.
#

initial_fuel_density = 10431.0 #95% TD (TD = 10980)

[GlobalParams]
  density = ${initial_fuel_density}
  order = SECOND
  family = LAGRANGE
  energy_per_fission = 3.2e-11 # J/fission (205 Mev)
  displacements = 'disp_x disp_y'
  temperature = temp
[]

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

[Mesh]
  coord_type = RZ
  [smeared_pellet_mesh]
    type = FuelPinMeshGenerator
    clad_mesh_density = customize
    pellet_mesh_density = customize
    ny_p = 1
    nx_p = 1
    nx_c = 1
    ny_cu = 1
    ny_c = 1
    ny_cl = 1
    clad_thickness = 5.6e-4
    pellet_outer_radius = 0.0041
    pellet_height = 0.01
    pellet_quantity = 1
    clad_bot_gap_height = 1e-3
    bottom_clad_height = 2.24e-3
    top_clad_height = 2.24e-3
    clad_gap_width = 8e-5
    plenum_fuel_ratio = 0.150
    elem_type = QUAD8
  []
  partitioner = centroid
  centroid_partitioner_direction = y
  patch_size = 5
[]

[Variables]
  [temp]
    initial_condition = 298
  []
[]

[AuxVariables]
  [fission_rate]
    block = '3'
  []
  [burnup]
    block = '3'
  []
  [grain_radius]
    block = '3'
    initial_condition = 5e-6
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 1.0e4 1.0e8'
    y = '0  1.0  1.0'
    scale_factor = 20e3 # 20 kW/m peak power.
  []
  [coolant_pressure_ramp]
    type = PiecewiseLinear
    x = '0 10000'
    y = '0 1'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [fuel]
    block = 3
    strain = FINITE
    incremental = true
    add_variables = true
    decomposition_method = EigenSolution
    extra_vector_tags = 'ref'
    eigenstrain_names = fuel_thermal_strain
  []
  [clad]
    block = 1
    strain = FINITE
    incremental = true
    add_variables = true
    decomposition_method = EigenSolution
    extra_vector_tags = 'ref'
    eigenstrain_names = clad_thermal_strain
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_source_fuel]
    type = NeutronHeatSource
    variable = temp
    block = '3'
    fission_rate = fission_rate
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [fissionrate]
    type = FissionRateGeneral
    fission_rate_formulation = GENERIC
    variable = fission_rate
    block = '3'
    value = 5.3548e+14
    fission_rate_function = power_history
  []
  [burnup]
    type = BurnupAux
    variable = burnup
    block = '3'
    fission_rate = fission_rate
    molecular_weight = 0.270
  []
  [grain_radius]
    type = GrainRadiusAux
    block = '3'
    variable = grain_radius
    temperature = temp
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    penalty = 1e+14 #1e7
    model = frictionless
    tangential_tolerance = 5e-4
    normal_smoothing_distance = 0.1
    normalize_penalty = true
  []
[]

[ThermalContact]
  [pellet_clad_thermal]
    type = GasGapHeatTransfer
    variable = temp
    primary = 5
    secondary = 10
    gas_released = 'fis_gas_released he_prod'
    initial_moles = initial_moles
    jump_distance_model = LANNING
    layer_thickness = layer_thickness
    plenum_pressure = plenum_pressure
    contact_pressure = contact_pressure
    initial_gas_types = Xe
    initial_fractions = 1
    released_gas_types = 'Kr Xe;
                          He'
    released_fractions = '0.153 0.847;
                          1'
    roughness_coef = 3.2
    roughness_secondary = 1e-6
    roughness_primary = 2e-6
    emissivity_primary = 0.8
    emissivity_secondary = 0.8
    quadrature = true
    normal_smoothing_distance = 0.1
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 1003
    value = 0.0
  []
  [no_y_clad_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 1
    value = 0.0
  []
  [no_y_fuel_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 20
    value = 0.0
  []
  [no_x_fuel]
    type = DirichletBC
    variable = disp_x
    boundary = 1005
    value = 0.0
  []
  [Clad_Temp]
    type = DirichletBC
    variable = temp
    boundary = '2'
    value = 580.0
  []
  [Pressure]
    [coolantPressure]
      boundary = '2'
      factor = 15.5e6
      function = coolant_pressure_ramp
    []
  []
  [PlenumPressure]
    [plenumPressure]
      boundary = 9
      initial_pressure = 0.50e6
      startup_time = 0.0
      material_input = 'fis_gas_released he_prod'
      output_initial_moles = initial_moles
      temperature = interior_temp
      volume = gas_volume
      output = plenum_pressure
    []
  []
[]

[Materials]
  [fuel_thermal]
    type = UO2Thermal
    block = '3'
    temperature = temp
    burnup = burnup
    thermal_conductivity_model = NFIR
  []
  [fuel_elasticity_tensor]
    type = UO2ElasticityTensor
    block = 3
  []
  [fuel_elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 3
  []
  [fuel_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = 3
    thermal_expansion_coeff = 10.0e-6
    stress_free_temperature = 298
    eigenstrain_name = fuel_thermal_strain
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = '3'
    temperature = temp
    fission_rate = fission_rate
    grain_radius = grain_radius
    gbs_model = true
    burnup = burnup
    diff_coeff_option = TURNBULL_D1_D2
  []
  [clad_thermal]
    type = HeatConductionMaterial
    block = 1
    thermal_conductivity = 16.0
    specific_heat = 330.0
  []
  [fclad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = 1
    youngs_modulus = 7.5e10
    poissons_ratio = 0.3
  []
  [clad_elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 1
  []
  [clad_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = 1
    thermal_expansion_coeff = 5.0e-6
    stress_free_temperature = 298
    eigenstrain_name = clad_thermal_strain
  []
  [clad_density]
    type = StrainAdjustedDensity
    block = 1
    strain_free_density = 6551.0
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = '3'
    strain_free_density = ${initial_fuel_density}
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

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

  line_search = 'none'

  l_max_its = 25
  nl_max_its = 40
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-8

  dtmax = 1.0e6
  dtmin = 1.0
  end_time = 5.3e7 # 1.7 years (~3% burnup)

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e3
    optimal_iterations = 30
    iteration_window = 4
    time_t = '0   1e4 1e8'
    time_dt = '1e4 1e5 1e6'
    timestep_limiting_function = power_history
    force_step_every_function_point = true
  []
  [Quadrature]
    order = fifth
    side_order = seventh
  []
  verbose = true
[]

[Postprocessors]
  [clad_inner_vol]
    type = InternalVolume
    boundary = 7
    execute_on = 'initial linear'
  []
  [pellet_volume]
    type = InternalVolume
    boundary = 8
    execute_on = 'initial linear'
  []
  [gas_volume]
    type = InternalVolume
    boundary = 9
    execute_on = 'initial linear'
  []
  [interior_temp]
    type = SideAverageValue
    boundary = 9 # cladding interior and pellet exterior
    variable = temp
    execute_on = 'initial linear'
  []
  [fis_gas_produced] # fission gas produced (moles)
    type = ElementIntegralFisGasGeneratedSifgrs
    block = '3'
  []
  [fis_gas_released]
    type = ElementIntegralFisGasReleasedSifgrs
    block = '3'
  []
  [fis_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    block = '3'
  []
  [fis_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    block = '3'
  []
  [power_history]
    type = FunctionValuePostprocessor
    function = power_history
  []
  [dt]
    type = TimestepSize
  []
  [residual]
    type = Residual
  []
  [nl_its]
    type = NumNonlinearIterations
  []
  [lin_its]
    type = NumLinearIterations
  []
  [average_burnup]
    type = ElementAverageValue
    block = '3'
    variable = burnup
  []
  [burnup]
    type = ElementAverageValue
    block = '3'
    variable = burnup
  []
  [average_fissionrate]
    type = ElementAverageValue
    block = '3'
    variable = fission_rate
  []
  [rod_total_power]
    type = ElementIntegralPower
    variable = temp
    fission_rate = fission_rate
    block = '3'
  []
  [rod_input_power]
    type = FunctionValuePostprocessor
    function = power_history
    scale_factor = 0.01 # change: length of fuel stack in meters (1 pellet height)
  []
  [he_prod]
    type = IFBAHeProduction
    zrb2_load = 1.181e-4
    ifba_len = 1.0e-2
    b10_enrich = 0.50
    zrb2_rel_dens = 0.7
    model = burnup
    u235_enrich = 0.045
    burnup = average_burnup
  []
[]

[Outputs]
  time_step_interval =  1
  exodus = false
  [console]
    type = Console
    solve_log = true
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    show = 'average_burnup burnup he_prod interior_temp plenum_pressure'
    file_base = fill_gas_xenon_w_ifba_check
  []
  [out]
    type = CSV
    delimiter = ' '
  []
[]

(assessment/LWR/validation/LOCA_Studsvik/analysis/rod_196/Studsvik_196_part1.i)


initial_fuel_density = 10431.0

[GlobalParams]
  density = ${initial_fuel_density}
  initial_porosity = 0.05
  order = SECOND
  family = LAGRANGE
  energy_per_fission = 3.2e-11 # J/fission
  volumetric_locking_correction = false
  displacements = 'disp_x disp_y'
[]

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

[Mesh]
  coord_type = RZ
  [smeared_mesh]
    type = FuelPinMeshGenerator
    clad_top_gap_height = 0.0248576
    pellet_height = 0.2606424
    pellet_quantity = 1
    clad_bot_gap_height = 0.0145
    pellet_outer_radius = 3.92e-3
    clad_gap_width = 80e-6
    clad_thickness = 0.57e-3
    clad_mesh_density = customize
    pellet_mesh_density = customize
    nx_c = 5
    ny_c = 50
    nx_p = 11
    ny_p = 60
    elem_type = QUAD8
  []
  patch_update_strategy = auto
  patch_size = 10 # For contact algorithm
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[Variables]
  # Define dependent variables and initial conditions
  [temperature]
    initial_condition = 295.0 # set initial temp to coolant inlet
  []
[]

[AuxVariables]
  # Define auxilary variables
  [fast_neutron_flux]
    block = clad
  []
  [fast_neutron_fluence]
    block = clad
  []
  [grain_radius]
    block = pellet
    initial_condition = 10e-6
  []
  [creep_strain_rate]
    order = CONSTANT
    family = MONOMIAL
  []
  [effective_creep_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [hoop_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase] # Fraction of beta phase in Zry
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness] # ZrO2 scale thickness (m)
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total] # Current oxigen weight fraction (oxide+metal) (/)
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total] # Gained oxygen weight fraction (oxide+metal) (/)
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress] # Hoop stress at cladding burst
    order = CONSTANT
    family = MONOMIAL
  []
  [burst] # Did cladding burst occur?
    order = CONSTANT
    family = MONOMIAL
  []

  [gap_cond]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    data_file = power_history.csv
    format = columns
    scale_factor = 1
  []
  [axial_peaking_factors]
    type = ParsedFunction
    expression = 1
  []
  [pressure_ramp] # reads and interpolates input data defining amplitude curve for fill gas pressure
    type = PiecewiseLinear
    x = '-200 0  86400  47386400  47472800  47559200  47645600  94945600  95032000'
    y = '0.0065371 1 1 1 1 1 1 1 0.0065371'
    scale_factor = 15.5e6
  []
  [forced_times]
    type = PiecewiseLinear
    data_file = timestep_limiting.csv
    scale_factor = 1
    format = columns
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [fuel]
    block = pellet
    add_variables = true
    strain = FINITE
    eigenstrain_names = 'fuel_thermal_eigenstrain fuel_relocation_eigenstrain fuel_volumetric_eigenstrain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
    decomposition_method = EigenSolution
    extra_vector_tags = 'ref'
    temperature = temperature
  []
  [clad]
    block = clad
    add_variables = true
    strain = FINITE
    eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_eigenstrain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz creep_strain_zz strain_zz'
    extra_vector_tags = 'ref'
    decomposition_method = EigenSolution
    temperature = temperature
  []
[]

[Kernels]
  [gravity]
    type = Gravity
    variable = disp_y
    value = -9.81
  []
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_source]
    type = NeutronHeatSource
    variable = temperature
    extra_vector_tags = 'ref'
    block = pellet
    burnup_function = burnup
  []
[]

[Burnup]
  [burnup]
    block = pellet
    rod_ave_lin_pow = power_history # using the power function defined above
    axial_power_profile = axial_peaking_factors # using the axial power profile function defined above
    num_radial = 80
    num_axial = 11
    fuel_pin_geometry = fuel_pin_geometry
    fuel_volume_ratio = 1.0 # for use with dished pellets (ratio of actual volume to cylinder volume)
    order = CONSTANT
    family = MONOMIAL
    RPF = RPF
    isotopes = 'U235 U238 Pu239 Pu240 Pu241 Pu242'
    isotope_fractions = '0.05 0.95 0 0 0 0'
  []
[]

[AuxKernels]
  # Define auxilliary kernels for each of the aux variables
  [fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    block = clad
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    factor = 3e13
    execute_on = timestep_begin
  []
  [fast_neutron_fluence]
    type = FastNeutronFluenceAux
    variable = fast_neutron_fluence
    block = clad
    fast_neutron_flux = fast_neutron_flux
    execute_on = timestep_begin
  []
  [grain_radius]
    type = GrainRadiusAux
    block = pellet
    variable = grain_radius
    temperature = temperature
    execute_on = linear
  []
  [creep_strain_rate]
    type = MaterialRealAux
    property = creep_rate
    variable = creep_strain_rate
    block = clad
  []
  [effective_creep_strain]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = effective_creep_strain
    block = clad
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    block = clad
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = MaterialRealAux
    boundary = 2
    variable = scale_thickness
    property = oxide_scale_thickness
  []
  [ofract_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
  []
  [ofgain_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    boundary = 2
    variable = burst_stress
    property = burst_stress
  []
  [hasburst]
    type = MaterialRealAux
    boundary = 2
    variable = burst
    property = failed
    execute_on = timestep_end
  []
  [conductance]
    type = MaterialRealAux
    property = gap_conductance
    variable = gap_cond
    boundary = 10
    execute_on = 'linear'
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    formulation = kinematic
    model = frictionless
    penalty = 1e7
  []
[]

[ThermalContact]
  [thermal_contact]
    type = GasGapHeatTransfer
    variable = temperature
    primary = 5
    secondary = 10
    initial_moles = initial_moles
    gas_released = 'fission_gas_released he_prod'
    released_gas_types = 'Kr Xe;
                          He'
    released_fractions = '0.153 0.847;
                          1'
    quadrature = true
    contact_pressure = contact_pressure
    refab_gas_types = He
    refab_fractions = 1
    refab_time = 95032000
    refab_type = 0
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 12
    value = 0.0
  []
  [no_y_clad_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = '1'
    value = 0.0
  []
  [no_y_fuel_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = '1020'
    value = 0.0
  []
  [Pressure]
    [coolantPressure]
      boundary = '1 2 3'
      function = pressure_ramp
    []
  []
  [PlenumPressure]
    [plenumPressure]
      boundary = 9
      initial_pressure = 3.44738e6
      startup_time = 0
      R = 8.3143
      output_initial_moles = initial_moles
      temperature = plenum_temp
      volume = plenum_volume
      material_input = 'fission_gas_released he_prod'
      output = plenum_pressure
      refab_time = 95032000
      refab_pressure = 8.2e6
      refab_temperature = 295.0
      refab_volume = 1.04e-05
      cladding_failure_status = burst
      equilibrium_pressure = equilibrium_pressure
      additional_volumes = additional_volume
      temperature_of_additional_volumes = addition_temperature
    []
  []
[]

[UserObjects]
  [fuel_pin_geometry]
    type = FuelPinGeometry
  []
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
[]

[PlenumTemperature]
  [plenum_temp]
    boundary = 5
    inner_surfaces = '5'
    outer_surfaces = '10'
    temperature = temperature
  []
[]

[CoolantChannel]
  [convective_clad_surface] # apply convective boundary to clad outer surface
    boundary = 2
    variable = temperature
    inlet_temperature = 580
    inlet_pressure = 15.5e6 # Pa
    inlet_massflux = 3800 # kg/m^2-sec
    rod_diameter = 0.00914 # m
    rod_pitch = 1.26e-2 # m
    compute_enthalpy = false
    linear_heat_rate = power_history
    axial_power_profile = axial_peaking_factors
    output_properties = 'coolant_channel_htype coolant_channel_hmode'
  []
[]

[Materials]
  [uo2_pulverization]
    type = UO2Pulverization
    block = pellet
    layered_average_contact_pressure = contact_pressure
    temperature = temperature
    burnup_function = burnup
    output_properties = pulverized
    outputs = all
  []

  # Define material behavior models and input material property data
  [fuel_thermal] # temperature and burnup dependent thermal properties of UO2 (BISON kernel)
    type = UO2Thermal
    block = pellet
    thermal_conductivity_model = NFIR
    temperature = temperature
    burnup_function = burnup
  []
  [fuel_elasticity_tensor]
    type = UO2IsotropicDamageElasticityTensor
    block = pellet
    fragmentation_model = BARANI
    temperature = temperature
    rod_ave_lin_pow = power_history
  []
  [fuel_elastic_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'fuel_creep'
    block = pellet
  []
  [fuel_creep]
    type = UO2CreepUpdate
    block = pellet
    temperature = temperature
    fission_rate = fission_rate
    initial_grain_radius = 10.0e-6
    oxygen_to_metal_ratio = 2.0
  []
  [fuel_relocation]
    type = UO2RelocationEigenstrain
    block = pellet
    burnup_function = burnup
    fuel_pin_geometry = fuel_pin_geometry
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    burnup_relocation_stop = 0.024
    relocation_activation1 = 5000
    relocation_model = ESCORE_modified
    eigenstrain_name = fuel_relocation_eigenstrain
  []
  [fuel_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = pellet
    thermal_expansion_coeff = 10.0e-6
    temperature = temperature
    stress_free_temperature = 295.0
    eigenstrain_name = fuel_thermal_eigenstrain
  []
  [fuel_volumetric_swelling]
    type = UO2VolumetricSwellingEigenstrain
    gas_swelling_model_type = SIFGRS
    block = pellet
    temperature = temperature
    burnup_function = burnup
    initial_fuel_density = 10431.0
    eigenstrain_name = fuel_volumetric_eigenstrain
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = pellet
    temperature = temperature
    burnup_function = burnup
    grain_radius = grain_radius
    gbs_model = true
  []

  [clad_density]
    type = StrainAdjustedDensity
    block = clad
    strain_free_density = 6550.
  []
  [clad_thermal]
    block = clad
    type = ZryThermal
    temperature = temperature
  []
  [clad_thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = clad
    temperature = temperature
    stress_free_temperature = 295.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [clad_elasticity_tensor]
    type = ZryElasticityTensor
    block = clad
    temperature = temperature
  []
  [zry_thermal_creep]
    type = ZryCreepLOCAUpdate
    block = clad
    temperature = temperature
    model_irradiation_creep = true
    model_primary_creep = true
    model_thermal_creep = true
    max_inelastic_increment = 5e-4
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = zirlo
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'zry_thermal_creep'
    block = clad
  []
  [clad_irradiation_growth]
    type = ZryIrradiationGrowthEigenstrain
    block = clad
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = zirlo
    eigenstrain_name = clad_irradiation_eigenstrain
  []
  [clad_phase]
    type = ZrPhase
    block = clad
    temperature = temperature
    numerical_method = 2
  []
  [clad_oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 4.18e-03
    clad_outer_radius = 4.75e-03
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = overstrain
    #    effective_strain_rate_creep = creep_strain_rate
    #    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = stress_zz
    hoop_creep_strain = creep_strain_zz
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfract_total
    temperature = temperature
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = pellet
    strain_free_density = ${initial_fuel_density}
  []
[]

[Dampers]
  [limitT]
    type = BoundingValueElementDamper
    min_value = 290.0
    max_value = 3000.0
    variable = temperature
  []
  [limitX]
    type = MaxIncrement
    max_increment = 1e-5
    variable = disp_x
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

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

  line_search = 'none'

  l_max_its = 50
  l_tol = 8e-3
  nl_max_its = 15
  nl_rel_tol = 1e-4
  nl_abs_tol = 1e-8
  start_time = -10
  n_startup_steps = 1
  end_time = 95032000
  dtmax = 1e6
  dtmin = 1e-6

  [TimeStepper]
    type = IterationAdaptiveDT
    timestep_limiting_postprocessor = material_timestep
    dt = 10
    optimal_iterations = 20
    iteration_window = 4
    linear_iteration_ratio = 100
    growth_factor = 2
    cutback_factor = .5
    timestep_limiting_function = forced_times
    force_step_every_function_point = true
  []
  [Quadrature]
    order = FIFTH
    side_order = SEVENTH
  []
[]

[Postprocessors]
  [ave_temp_interior]
    type = SideAverageValue
    boundary = 9
    variable = temperature
    execute_on = 'initial linear'
  []
  [clad_inner_vol]
    type = InternalVolume
    boundary = 7
    #outputs = exodus
    execute_on = 'initial timestep_end'
  []
  [fission_gas_produced] # fission gas produced (moles)
    type = ElementIntegralFisGasGeneratedSifgrs
    block = pellet
    execute_on = 'linear'
  []
  [fission_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    block = pellet
    outputs = exodus
    execute_on = 'linear'
  []
  [fission_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    block = pellet
    outputs = exodus
    execute_on = 'linear'
  []
  [flux_from_clad] # area integrated heat flux from the cladding
    type = SideDiffusiveFluxIntegral
    variable = temperature
    boundary = 5
    diffusivity = thermal_conductivity
  []
  [flux_from_fuel] # area integrated heat flux from the fuel
    type = SideDiffusiveFluxIntegral
    variable = temperature
    boundary = 10
    diffusivity = thermal_conductivity
  []

  [rod_total_power]
    type = ElementIntegralPower
    variable = temperature
    burnup_function = burnup
    block = pellet
  []
  [rod_input_power]
    type = FunctionValuePostprocessor
    function = power_history
    scale_factor = 0.1186 # rod height
  []
  [max_fuel_temp]
    type = NodalExtremeValue
    block = pellet
    value_type = max
    variable = temperature
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = clad
    value_type = max
    variable = temperature
  []
  [max_clad_hoop_strain]
    type = ElementExtremeValue
    block = clad
    value_type = max
    variable = strain_zz
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = clad
  []
  [burst]
    type = ElementExtremeValue
    value_type = max
    variable = burst
    block = clad
    execute_on = 'initial timestep_end'
  []
  [he_prod]
    type = IFBAHeProduction
    b10_load = 9.27165354e-5
    b10_enrich = 0.5
    burnup = average_burnup
    zrb2_thick = 10e-6
    fuel_out_rad = 9.32e-3
    ifba_len = 0.3
    u235_enrich = 0.05
  []
  [volume_pulverized]
    type = ElementIntegralMaterialProperty
    mat_prop = pulverized
    block = pellet
  []
  [max_fuel_temp_periphery]
    type = NodalExtremeValue
    value_type = max
    variable = temperature
    boundary = 10
  []
  [additional_volume]
    type = FunctionValuePostprocessor
    function = 8.5e-6
    execute_on = 'initial linear'
  []
  [addition_temperature]
    type = FunctionValuePostprocessor
    function = 300.0
    execute_on = 'initial linear'
  []
  [equilibrium_pressure]
    type = FunctionValuePostprocessor
    function = 101325.0
    execute_on = 'initial linear'
  []
[]

[PerformanceMetricOutputs]
[]

[StandardLWRFuelRodOutputs]
  temperature = temperature
  fuel_pellet_blocks = 3
[]

[Outputs]
  perf_graph = true
  exodus = true
  color = false
  csv = true
  [checkpoint]
    type = Checkpoint
    num_files = 2
  []
  [chkfile]
    type = CSV
    execute_on = FINAL
    show = 'volume_pulverized'
  []
[]

(assessment/LWR/validation/LOCA_Studsvik/analysis/rod_196/Studsvik_196_part2_1p5d_fr_ffrd.i)


initial_fuel_density = 10431.0

[GlobalParams]
  density = ${initial_fuel_density}
  initial_porosity = 0.05
  order = SECOND
  family = LAGRANGE
  energy_per_fission = 3.2e-11 # J/fission
  volumetric_locking_correction = false
  displacements = 'disp_x'
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
  restart_file_base = 'Studsvik_196_part1_1p5d_fr_ffrd_checkpoint_cp/LATEST'
[]

[Mesh]
  coord_type = RZ
  [layered1D_mesh]
    type = Layered1DMeshGenerator
    slices_per_block = 10
    clad_gap_width = 80e-6
    plenum_height = 0.0393576
    pellet_outer_radius = 3.92e-3
    clad_thickness = 0.57e-3
    fuel_height = 0.2606424
    # nx_c = 2
    # nx_p = 11
    elem_type = EDGE3
  []
  patch_update_strategy = auto
  patch_size = 10 # For contact algorithm
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[Variables]
  [temperature]
  []
[]

[AuxVariables]
  # Define auxilary variables
  [strain_yy_0]
    order = CONSTANT
    family = MONOMIAL
  []
  [fast_neutron_flux]
    block = clad
  []
  [fast_neutron_fluence]
    block = clad
  []
  [grain_radius]
    block = fuel
  []
  [creep_strain_rate]
    order = CONSTANT
    family = MONOMIAL
  []
  [effective_creep_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [hoop_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase] # Fraction of beta phase in Zry
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness] # ZrO2 scale thickness (m)
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total] # Current oxigen weight fraction (oxide+metal) (/)
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total] # Gained oxygen weight fraction (oxide+metal) (/)
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress] # Hoop stress at cladding burst
    order = CONSTANT
    family = MONOMIAL
  []
  [burst] # Did cladding burst occur?
    order = CONSTANT
    family = MONOMIAL
  []

  [gap_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [tangential_contact_pressure_aux]
    block = fuel
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    data_file = power_history.csv
    format = columns
    scale_factor = 1
  []
  [axial_peaking_factors]
    type = ParsedFunction
    expression = 1
  []
  [pressure_ramp] # reads and interpolates input data defining amplitude curve for fill gas pressure
    type = PiecewiseLinear
    x = '-200 0  86400  47386400  47472800  47559200  47645600  94945600  95032000'
    y = '0.0065371 1 1 1 1 1 1 1 0.0065371'
    scale_factor = 15.5e6
  []
  [clad_surface_temperature]
    type = PiecewiseBilinear
    axis = 1
    data_file = clad_temperature.csv
  []
  [forced_times]
    type = PiecewiseLinear
    data_file = timestep_limiting.csv
    scale_factor = 1
    format = columns
  []
  [clad_axial_pressure]
    type = CladdingAxialPressureFunction
    plenum_pressure = plenum_pressure
    coolant_pressure = pressure_ramp
    coolant_pressure_scaling_factor = 1.0
    fuel_pin_geometry = fuel_pin_geometry
  []
  [fuel_axial_pressure]
    type = ParsedFunction
    expression = plenum_pressure
    symbol_names = plenum_pressure
    symbol_values = plenum_pressure
  []
[]

[Physics]
  [SolidMechanics]
    [Layered1D]
      [fuel]
        block = fuel
        add_variables = true
        add_scalar_variables = true
        strain = FINITE
        out_of_plane_strain_name = strain_yy
        eigenstrain_names = 'fuel_thermal_eigenstrain fuel_volumetric_eigenstrain axial_relocation_eigenstrain'
        generate_output = 'stress_xx stress_yy stress_zz vonmises_stress strain_xx hoop_stress creep_strain_zz strain_zz'
        extra_vector_tags = 'ref'
        fuel_pin_geometry = fuel_pin_geometry
        group_scalar_vars_in_reference_residual = true
        mesh_generator = layered1D_mesh
        decomposition_method = EigenSolution
        temperature = temperature
        out_of_plane_pressure_function = fuel_axial_pressure
        layer_friction_user_object = 1DFriction_secondary
      []
      [clad]
        block = clad
        add_variables = true
        add_scalar_variables = true
        strain = FINITE
        out_of_plane_strain_name = strain_yy
        eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_eigenstrain'
        generate_output = 'stress_xx stress_yy stress_zz vonmises_stress strain_xx hoop_stress creep_strain_zz strain_zz'
        extra_vector_tags = 'ref'
        fuel_pin_geometry = fuel_pin_geometry
        group_scalar_vars_in_reference_residual = true
        mesh_generator = layered1D_mesh
        decomposition_method = EigenSolution
        temperature = temperature
        out_of_plane_pressure_function = clad_axial_pressure
        layer_friction_user_object = 1DFriction_primary
      []
    []
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_source]
    type = NeutronHeatSource
    variable = temperature
    extra_vector_tags = 'ref'
    block = fuel
    burnup_function = burnup
    axial_relocation_object = axial_relocation
  []
[]

[Burnup]
  [burnup]
    block = fuel
    rod_ave_lin_pow = power_history # using the power function defined above
    axial_power_profile = axial_peaking_factors # using the axial power profile function defined above
    num_radial = 80
    num_axial = 11
    fuel_pin_geometry = fuel_pin_geometry
    fuel_volume_ratio = 1.0 # for use with dished pellets (ratio of actual volume to cylinder volume)
    order = CONSTANT
    family = MONOMIAL
    RPF = RPF
    isotopes = 'U235 U238 Pu239 Pu240 Pu241 Pu242'
    isotope_fractions = '0.05 0.95 0 0 0 0'
  []
[]

[AuxKernels]
  # Define auxilliary kernels for each of the aux variables
  [tangential_contact_pressure_aux]
    type = SpatialUserObjectAux
    variable = tangential_contact_pressure_aux
    user_object = 1DFriction_secondary
    block = fuel
    execute_on = 'TIMESTEP_END'
  []
  [fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    block = clad
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    factor = 3e13
    execute_on = timestep_begin
  []
  [fast_neutron_fluence]
    type = FastNeutronFluenceAux
    variable = fast_neutron_fluence
    block = clad
    fast_neutron_flux = fast_neutron_flux
    execute_on = timestep_begin
  []
  [grain_radius]
    type = GrainRadiusAux
    block = fuel
    variable = grain_radius
    temperature = temperature
    execute_on = linear
  []
  [creep_strain_rate]
    type = MaterialRealAux
    property = creep_rate
    variable = creep_strain_rate
    block = clad
  []
  [effective_creep_strain]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = effective_creep_strain
    block = clad
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    block = clad
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = MaterialRealAux
    boundary = 2
    variable = scale_thickness
    property = oxide_scale_thickness
  []
  [ofract_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
  []
  [ofgain_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    boundary = 2
    variable = burst_stress
    property = burst_stress
  []
  [hasburst]
    type = MaterialRealAux
    boundary = 2
    variable = burst
    property = failed
    execute_on = timestep_end
  []
  [conductance]
    type = MaterialRealAux
    property = gap_conductance
    variable = gap_cond
    boundary = 10
    execute_on = 'linear'
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    formulation = kinematic
    model = frictionless
    penalty = 1e7
  []
[]

[ThermalContact]
  [thermal_contact]
    type = GasGapHeatTransfer
    variable = temperature
    primary = 5
    secondary = 10
    initial_moles = initial_moles
    gas_released = 'fission_gas_released he_prod'
    released_gas_types = 'Kr Xe;
                          He'
    released_fractions = '0.153 0.847;
                          1'
    quadrature = true
    contact_pressure = contact_pressure
    refab_gas_types = He
    refab_fractions = 1
    refab_time = 95032000
    refab_type = 0
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 12
    value = 0.0
  []
  [Pressure]
    [coolantPressure]
      boundary = '2'
      function = pressure_ramp
    []
  []
  [PlenumPressure]
    [plenumPressure]
      boundary = 9
      initial_pressure = 3.44738e6
      startup_time = 0
      R = 8.3143
      output_initial_moles = initial_moles
      temperature = plenum_temp
      volume = plenum_volume
      material_input = 'fission_gas_released he_prod'
      output = plenum_pressure
      refab_time = 95032000
      refab_pressure = 8.2e6
      refab_temperature = 295.0
      refab_volume = 1.04e-05
      cladding_failure_status = burst
      equilibrium_pressure = equilibrium_pressure
      additional_volumes = additional_volume
      temperature_of_additional_volumes = addition_temperature
    []
  []

  [clad_temp]
    type = FunctionDirichletBC
    function = clad_surface_temperature
    variable = temperature
    boundary = 2
  []
[]

[UserObjects]
  [layered_average_hoop_strain]
    type = LayeredAverage
    block = clad
    num_layers = 10
    direction = y
    variable = strain_zz
  []
  [cladding_strain_yy]
    type = LayeredAverage
    block = clad
    num_layers = 11
    direction = y
    variable = strain_yy
    execute_on = 'initial timestep_end'
  []
  [fuel_strain_yy]
    type = LayeredAverage
    block = fuel
    num_layers = 10
    direction = y
    variable = strain_yy
    execute_on = 'initial timestep_end'
  []
  [1DContactStressOOP_fuel]
    type = Layered1DContactInterfaceStress
    direction = y
    stress_name = stress
    num_layers = 10
    # If we do not provide the numbers below, it will look at the mesh, in all blocks to set the layer number. Then, it will
    # be wrong because the cladding has more height and won't be able to identify layers in the fuel.
    direction_min = 0.01306
    direction_max = 0.24761028

    block = fuel
    execute_on = 'LINEAR NONLINEAR'
  []
  [1DContactStressOOP_cladding]
    type = Layered1DContactInterfaceStress
    direction = y
    stress_name = stress
    num_layers = 10
    # If we do not provide the numbers below, it will look at the mesh, in all blocks to set the layer number. Then, it will
    # be wrong because the cladding has more height and won't be able to identify layers in the fuel.
    direction_min = 0.01306
    direction_max = 0.24761028

    block = clad
    execute_on = 'LINEAR NONLINEAR'
  []
  [1DFriction_secondary]
    type = Layered1DFrictionalForce
    force_postaux = true
    contact_pressure = contact_pressure
    direction = y
    boundary = pellet_outer_radial_surface
    num_layers = 10
    interface_oop_stress_provider_fuel = 1DContactStressOOP_fuel
    interface_oop_stress_provider_cladding = 1DContactStressOOP_cladding
    is_secondary_side = true
    tangential_pressure = tangential_contact_pressure_aux
    friction_coefficient = 0.2
    thickness = 0.02606424
    penalty_factor = 1.0e13
    # If we do not provide the numbers below, it will look at the mesh, in all blocks to set the layer number. Then, it will
    # be wrong because the cladding has more height and won't be able to identify layers in the fuel.
    direction_min = 0.01306
    direction_max = 0.24761028

    scalar_var_name_base_fuel = scalar_strain_yy_fuel
    scalar_num_variable_fuel = 10
    scalar_var_name_base_cladding = scalar_strain_yy_clad
    scalar_num_variable_cladding = 10

    execute_on = 'LINEAR NONLINEAR'
  []
  [1DFriction_primary]
    type = Layered1DFrictionalForce
    force_postaux = true
    contact_pressure = contact_pressure
    direction = y

    boundary = clad_inside_right
    num_layers = 10
    # If we do not provide the numbers below, it will look at the mesh, in all blocks to set the layer number. Then, it will
    # be wrong because the cladding has more height and won't be able to identify layers in the fuel.
    direction_min = 0.0165094
    direction_max = 0.24761028

    interface_oop_stress_provider_fuel = 1DContactStressOOP_fuel
    interface_oop_stress_provider_cladding = 1DContactStressOOP_cladding
    is_secondary_side = false
    secondary_side_frictional_user_object = 1DFriction_secondary
    friction_coefficient = 0.2
    thickness = 0.02606424
    penalty_factor = 1.0e13
    scalar_var_name_base_fuel = scalar_strain_yy_fuel
    scalar_num_variable_fuel = 10
    scalar_var_name_base_cladding = scalar_strain_yy_clad
    scalar_num_variable_cladding = 10

    execute_on = 'LINEAR NONLINEAR'
  []
  [terminator]
    type = Terminator
    expression = 'max_axial_relocation_strain > 0.25'
  []
[]

[PlenumTemperature]
  [plenum_temp]
    boundary = 5
    inner_surfaces = '5'
    outer_surfaces = '10'
    temperature = temperature
  []
[]

[CoolantChannel]
  [convective_clad_surface] # apply convective boundary to clad outer surface
    boundary = 2
    variable = temperature
    inlet_temperature = 580
    inlet_pressure = 15.5e6 # Pa
    inlet_massflux = 3800 # kg/m^2-sec
    rod_diameter = 0.00914 # m
    rod_pitch = 1.26e-2 # m
    compute_enthalpy = false
    linear_heat_rate = power_history
    axial_power_profile = axial_peaking_factors
    output_properties = 'coolant_channel_htype coolant_channel_hmode'
  []
[]

[Materials]
  [fuel_dispersal]
    type = UO2Dispersal
    block = fuel
    axial_relocation_object = axial_relocation
    layered_average_burnup = layered_average_burnup
    layered_average_hoop_strain = layered_average_hoop_strain
    dispersal_model = ONE_MM_TWO_PERCENT_STRAIN
  []

  # Define material behavior models and input material property data
  [fuel_thermal] # temperature and burnup dependent thermal properties of UO2 (BISON kernel)
    type = UO2Thermal
    block = fuel
    thermal_conductivity_model = NFIR
    temperature = temperature
    burnup_function = burnup
    axial_relocation_object = axial_relocation
    gap_thermal_conductivity = layered_average_gap_conductivity
  []
  [fuel_elasticity_tensor]
    type = UO2IsotropicDamageElasticityTensor
    block = fuel
    fragmentation_model = BARANI
    rod_ave_lin_pow = power_history
    temperature = temperature
    axial_relocation_object = axial_relocation
  []
  [fuel_elastic_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'fuel_creep'
    block = fuel
  []
  [fuel_creep]
    type = UO2CreepUpdate
    block = fuel
    temperature = temperature
    fission_rate = fission_rate
    initial_grain_radius = 10.0e-6
    oxygen_to_metal_ratio = 2.0
  []
  [fuel_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = fuel
    thermal_expansion_coeff = 10.0e-6
    temperature = temperature
    stress_free_temperature = 295.0
    eigenstrain_name = fuel_thermal_eigenstrain
  []
  [fuel_volumetric_swelling]
    type = UO2VolumetricSwellingEigenstrain
    gas_swelling_model_type = SIFGRS
    block = fuel
    temperature = temperature
    burnup_function = burnup
    initial_fuel_density = 10431.0
    eigenstrain_name = fuel_volumetric_eigenstrain
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = fuel
    temperature = temperature
    burnup_function = burnup
    grain_radius = grain_radius
    gbs_model = true
  []

  [clad_density]
    type = StrainAdjustedDensity
    block = clad
    strain_free_density = 6550.
  []
  [clad_thermal]
    block = clad
    type = ZryThermal
    temperature = temperature
  []
  [clad_thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = clad
    temperature = temperature
    stress_free_temperature = 295.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [clad_elasticity_tensor]
    type = ZryElasticityTensor
    block = clad
    temperature = temperature
  []
  [zry_thermal_creep]
    type = ZryCreepLOCAUpdate
    block = clad
    temperature = temperature
    model_irradiation_creep = true
    model_primary_creep = true
    model_thermal_creep = true
    max_inelastic_increment = 5e-4
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = zirlo
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'zry_thermal_creep'
    block = clad
  []
  [clad_irradiation_growth]
    type = ZryIrradiationGrowthEigenstrain
    block = clad
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = zirlo
    eigenstrain_name = clad_irradiation_eigenstrain
  []
  [clad_phase]
    type = ZrPhase
    block = clad
    temperature = temperature
    numerical_method = 2
  []
  [clad_oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 4.18e-03
    clad_outer_radius = 4.75e-03
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = overstrain
    hoop_stress = hoop_stress
    hoop_creep_strain = creep_strain_zz
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfract_total
    temperature = temperature
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = fuel
    strain_free_density = ${initial_fuel_density}
  []
[]

[VectorPostprocessors]
  [cladding_outer]
    type = NodalValueSampler
    boundary = 5
    variable = disp_x
    sort_by = y
  []
[]

[AxialRelocation]
  [relocation]
    rod_ave_lin_pow = power_history
    axial_direction = y
    fuel_blocks = fuel
    clad_blocks = clad
    contact_pressure_variable = contact_pressure
    out_of_plane_strain_variable = strain_yy_0
    penetration_variable = penetration
    clad_inner_volume_addition = 0
    burnup_variable = burnup
    temperature = temperature
    axial_relocation_output_options = MASS_FRACTION
    mesh_generator = layered1D_mesh
    # CHANGE
    gap_thickness_threshold = 0.000050
  []
[]

[Postprocessors]
  [volume_fuel_dispersed]
    type = LayeredElementIntegralMaterialProperty
    block = fuel
    mat_prop = dispersed
    fuel_pin_geometry = fuel_pin_geometry
    execute_on = 'initial timestep_end'
  []
  [mass_fuel_dispersed]
    type = ParsedPostprocessor
    pp_names = volume_fuel_dispersed
    expression = '10431 * volume_fuel_dispersed'
    execute_on = 'initial timestep_end'
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

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

  line_search = 'none'

  l_max_its = 50
  l_tol = 8e-3
  nl_max_its = 15
  nl_rel_tol = 1e-4
  nl_abs_tol = 1e-8
  n_startup_steps = 1
  end_time = 95033429.6
  dtmax = 20
  dtmin = 1e-6

  [TimeStepper]
    type = IterationAdaptiveDT
    timestep_limiting_postprocessor = material_timestep
    dt = 10
    optimal_iterations = 20
    iteration_window = 4
    linear_iteration_ratio = 100
    growth_factor = 2
    cutback_factor = .5
    timestep_limiting_function = forced_times
    force_step_every_function_point = true
  []
  [Quadrature]
    order = FIFTH
    side_order = SEVENTH
  []
[]

[Postprocessors]
  [ave_temp_interior]
    type = SideAverageValue
    boundary = 9
    variable = temperature
    execute_on = 'initial linear'
  []
  [clad_inner_vol]
    type = InternalVolume
    boundary = 7
    #outputs = exodus
    execute_on = 'initial timestep_end'
  []
  [fission_gas_produced] # fission gas produced (moles)
    type = ElementIntegralFisGasGeneratedSifgrs
    block = fuel
    execute_on = 'linear'
  []
  [fission_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    block = fuel
    outputs = exodus
    execute_on = 'linear'
  []
  [fission_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    block = fuel
    outputs = exodus
    execute_on = 'linear'
  []
  [flux_from_clad] # area integrated heat flux from the cladding
    type = SideDiffusiveFluxIntegral
    variable = temperature
    boundary = 5
    diffusivity = thermal_conductivity
  []
  [flux_from_fuel] # area integrated heat flux from the fuel
    type = SideDiffusiveFluxIntegral
    variable = temperature
    boundary = 10
    diffusivity = thermal_conductivity
  []

  [rod_total_power]
    type = ElementIntegralPower
    variable = temperature
    burnup_function = burnup
    block = fuel
  []
  [max_fuel_temp]
    type = NodalExtremeValue
    block = fuel
    value_type = max
    variable = temperature
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = clad
    value_type = max
    variable = temperature
  []
  [max_clad_hoop_strain]
    type = ElementExtremeValue
    block = clad
    value_type = max
    variable = strain_zz
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = clad
  []
  [max_axial_relocation_strain]
    type = ElementExtremeValue
    value_type = max
    variable = axial_relocation_strain
    block = fuel
    execute_on = 'initial timestep_end'
  []
  [he_prod]
    type = IFBAHeProduction
    b10_load = 9.27165354e-5
    b10_enrich = 0.5
    burnup = average_burnup
    zrb2_thick = 10e-6
    fuel_out_rad = 9.32e-3
    ifba_len = 0.3
    u235_enrich = 0.05
  []
  [burst]
    type = ElementExtremeValue
    value_type = max
    variable = burst
    block = clad
    execute_on = 'initial timestep_end'
  []
  [volume_pulverized]
    type = ElementIntegralMaterialProperty
    mat_prop = pulverized
    block = fuel
  []
  [max_fuel_temp_periphery]
    type = NodalExtremeValue
    value_type = max
    variable = temperature
    boundary = 10
  []
  [additional_volume]
    type = FunctionValuePostprocessor
    function = 8.5e-6
    execute_on = 'initial linear'
  []
  [addition_temperature]
    type = FunctionValuePostprocessor
    function = 300.0
    execute_on = 'initial linear'
  []
  [equilibrium_pressure]
    type = FunctionValuePostprocessor
    function = 101325.0
    execute_on = 'initial linear'
  []
[]

[PerformanceMetricOutputs]
[]

[StandardLWRFuelRodOutputs]
  temperature = temperature
  layered = true
  fuel_pellet_blocks = 'fuel'
  fuel_pin_geometry = fuel_pin_geometry
[]

[Outputs]
  perf_graph = true
  exodus = true
  color = false
  csv = true
  [chkfile]
    type = CSV
    execute_on = FINAL
    show = 'volume_pulverized'
  []
[]

(assessment/LWR/validation/LOCA_Studsvik/analysis/rod_196/Studsvik_196_part2.i)


initial_fuel_density = 10431.0

[GlobalParams]
  density = ${initial_fuel_density}
  initial_porosity = 0.05
  order = SECOND
  family = LAGRANGE
  energy_per_fission = 3.2e-11 # J/fission
  volumetric_locking_correction = false
  displacements = 'disp_x disp_y'
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
  restart_file_base = 'Studsvik_196_part1_checkpoint_cp/LATEST'
[]

[Mesh]
  coord_type = RZ
  [smeared_mesh]
    type = FuelPinMeshGenerator
    clad_top_gap_height = 0.0248576
    pellet_height = 0.2606424
    pellet_quantity = 1
    clad_bot_gap_height = 0.0145
    pellet_outer_radius = 3.92e-3
    clad_gap_width = 80e-6
    clad_thickness = 0.57e-3
    clad_mesh_density = customize
    pellet_mesh_density = customize
    nx_c = 5
    ny_c = 50
    nx_p = 11
    ny_p = 60
    elem_type = QUAD8
  []
  patch_update_strategy = auto
  patch_size = 10 # For contact algorithm
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[Variables]
  [temperature]
  []
[]

[AuxVariables]
  # Define auxilary variables
  [fast_neutron_flux]
    block = clad
  []
  [fast_neutron_fluence]
    block = clad
  []
  [grain_radius]
    block = pellet
  []
  [creep_strain_rate]
    order = CONSTANT
    family = MONOMIAL
  []
  [effective_creep_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [hoop_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase] # Fraction of beta phase in Zry
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness] # ZrO2 scale thickness (m)
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total] # Current oxigen weight fraction (oxide+metal) (/)
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total] # Gained oxygen weight fraction (oxide+metal) (/)
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress] # Hoop stress at cladding burst
    order = CONSTANT
    family = MONOMIAL
  []
  [burst] # Did cladding burst occur?
    order = CONSTANT
    family = MONOMIAL
  []

  [gap_cond]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    data_file = power_history.csv
    format = columns
    scale_factor = 1
  []
  [axial_peaking_factors]
    type = ParsedFunction
    expression = 1
  []
  [pressure_ramp] # reads and interpolates input data defining amplitude curve for fill gas pressure
    type = PiecewiseLinear
    x = '-200 0  86400  47386400  47472800  47559200  47645600  94945600  95032000'
    y = '0.0065371 1 1 1 1 1 1 1 0.0065371'
    scale_factor = 15.5e6
  []
  [clad_surface_temperature]
    type = PiecewiseBilinear
    axis = 1
    data_file = clad_temperature.csv
  []
  [forced_times]
    type = PiecewiseLinear
    data_file = timestep_limiting.csv
    scale_factor = 1
    format = columns
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [fuel]
    block = pellet
    add_variables = true
    strain = FINITE
    eigenstrain_names = 'fuel_thermal_eigenstrain fuel_relocation_eigenstrain fuel_volumetric_eigenstrain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
    decomposition_method = EigenSolution
    extra_vector_tags = 'ref'
    temperature = temperature
  []
  [clad]
    block = clad
    add_variables = true
    strain = FINITE
    eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_eigenstrain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz creep_strain_zz strain_zz'
    extra_vector_tags = 'ref'
    decomposition_method = EigenSolution
    temperature = temperature
  []
[]

[Kernels]
  [gravity]
    type = Gravity
    variable = disp_y
    value = -9.81
  []
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_source]
    type = NeutronHeatSource
    variable = temperature
    extra_vector_tags = 'ref'
    block = pellet
    burnup_function = burnup
  []
[]

[Burnup]
  [burnup]
    block = pellet
    rod_ave_lin_pow = power_history # using the power function defined above
    axial_power_profile = axial_peaking_factors # using the axial power profile function defined above
    num_radial = 80
    num_axial = 11
    fuel_pin_geometry = fuel_pin_geometry
    fuel_volume_ratio = 1.0 # for use with dished pellets (ratio of actual volume to cylinder volume)
    order = CONSTANT
    family = MONOMIAL
    RPF = RPF
    isotopes = 'U235 U238 Pu239 Pu240 Pu241 Pu242'
    isotope_fractions = '0.05 0.95 0 0 0 0'
  []
[]

[AuxKernels]
  # Define auxilliary kernels for each of the aux variables
  [fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    block = clad
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    factor = 3e13
    execute_on = timestep_begin
  []
  [fast_neutron_fluence]
    type = FastNeutronFluenceAux
    variable = fast_neutron_fluence
    block = clad
    fast_neutron_flux = fast_neutron_flux
    execute_on = timestep_begin
  []
  [grain_radius]
    type = GrainRadiusAux
    block = pellet
    variable = grain_radius
    temperature = temperature
    execute_on = linear
  []
  [creep_strain_rate]
    type = MaterialRealAux
    property = creep_rate
    variable = creep_strain_rate
    block = clad
  []
  [effective_creep_strain]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = effective_creep_strain
    block = clad
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    block = clad
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = MaterialRealAux
    boundary = 2
    variable = scale_thickness
    property = oxide_scale_thickness
  []
  [ofract_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
  []
  [ofgain_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    boundary = 2
    variable = burst_stress
    property = burst_stress
  []
  [hasburst]
    type = MaterialRealAux
    boundary = 2
    variable = burst
    property = failed
    execute_on = timestep_end
  []
  [conductance]
    type = MaterialRealAux
    property = gap_conductance
    variable = gap_cond
    boundary = 10
    execute_on = 'linear'
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    formulation = kinematic
    model = frictionless
    penalty = 1e7
  []
[]

[ThermalContact]
  [thermal_contact]
    type = GasGapHeatTransfer
    variable = temperature
    primary = 5
    secondary = 10
    initial_moles = initial_moles
    gas_released = 'fission_gas_released he_prod'
    released_gas_types = 'Kr Xe;
                          He'
    released_fractions = '0.153 0.847;
                          1'
    quadrature = true
    contact_pressure = contact_pressure
    refab_gas_types = He
    refab_fractions = 1
    refab_time = 95032000
    refab_type = 0
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 12
    value = 0.0
  []
  [no_y_clad_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = '1'
    value = 0.0
  []
  [no_y_fuel_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = '1020'
    value = 0.0
  []
  [Pressure]
    [coolantPressure]
      boundary = '1 2 3'
      function = pressure_ramp
    []
  []
  [PlenumPressure]
    [plenumPressure]
      boundary = 9
      initial_pressure = 3.44738e6
      startup_time = 0
      R = 8.3143
      output_initial_moles = initial_moles
      temperature = plenum_temp
      volume = plenum_volume
      material_input = 'fission_gas_released he_prod'
      output = plenum_pressure
      refab_time = 95032000
      refab_pressure = 8.2e6
      refab_temperature = 295.0
      refab_volume = 1.04e-05
      cladding_failure_status = burst
      equilibrium_pressure = equilibrium_pressure
      additional_volumes = additional_volume
      temperature_of_additional_volumes = addition_temperature
    []
  []

  [clad_temp]
    type = FunctionDirichletBC
    function = clad_surface_temperature
    variable = temperature
    boundary = 2
  []
[]

[UserObjects]
  [fuel_pin_geometry]
    type = FuelPinGeometry
  []
  # [terminator]
  #   type = Terminator
  #   expression = 'burst > 0'
  # []
[]

[PlenumTemperature]
  [plenum_temp]
    boundary = 5
    inner_surfaces = '5'
    outer_surfaces = '10'
    temperature = temperature
  []
[]

[CoolantChannel]
  [convective_clad_surface] # apply convective boundary to clad outer surface
    boundary = 2
    variable = temperature
    inlet_temperature = 580
    inlet_pressure = 15.5e6 # Pa
    inlet_massflux = 3800 # kg/m^2-sec
    rod_diameter = 0.00914 # m
    rod_pitch = 1.26e-2 # m
    compute_enthalpy = false
    linear_heat_rate = power_history
    axial_power_profile = axial_peaking_factors
    output_properties = 'coolant_channel_htype coolant_channel_hmode'
  []
[]

[Materials]
  [uo2_pulverization]
    type = UO2Pulverization
    block = pellet
    layered_average_contact_pressure = contact_pressure
    temperature = temperature
    burnup_function = burnup
    output_properties = pulverized
    outputs = all
  []

  # Define material behavior models and input material property data
  [fuel_thermal] # temperature and burnup dependent thermal properties of UO2 (BISON kernel)
    type = UO2Thermal
    block = pellet
    thermal_conductivity_model = NFIR
    temperature = temperature
    burnup_function = burnup
  []
  [fuel_elasticity_tensor]
    type = UO2IsotropicDamageElasticityTensor
    block = pellet
    fragmentation_model = BARANI
    rod_ave_lin_pow = power_history
    temperature = temperature
  []
  [fuel_elastic_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'fuel_creep'
    block = pellet
  []
  [fuel_creep]
    type = UO2CreepUpdate
    block = pellet
    temperature = temperature
    fission_rate = fission_rate
    initial_grain_radius = 10.0e-6
    oxygen_to_metal_ratio = 2.0
  []
  [fuel_relocation]
    type = UO2RelocationEigenstrain
    block = pellet
    burnup_function = burnup
    fuel_pin_geometry = fuel_pin_geometry
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    burnup_relocation_stop = 0.024
    relocation_activation1 = 5000
    relocation_model = ESCORE_modified
    eigenstrain_name = fuel_relocation_eigenstrain
  []
  [fuel_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = pellet
    thermal_expansion_coeff = 10.0e-6
    temperature = temperature
    stress_free_temperature = 295.0
    eigenstrain_name = fuel_thermal_eigenstrain
  []
  [fuel_volumetric_swelling]
    type = UO2VolumetricSwellingEigenstrain
    gas_swelling_model_type = SIFGRS
    block = pellet
    temperature = temperature
    burnup_function = burnup
    initial_fuel_density = 10431.0
    eigenstrain_name = fuel_volumetric_eigenstrain
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = pellet
    temperature = temperature
    burnup_function = burnup
    grain_radius = grain_radius
    gbs_model = true
  []

  [clad_density]
    type = StrainAdjustedDensity
    block = clad
    strain_free_density = 6550.
  []
  [clad_thermal]
    block = clad
    type = ZryThermal
    temperature = temperature
  []
  [clad_thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = clad
    temperature = temperature
    stress_free_temperature = 295.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [clad_elasticity_tensor]
    type = ZryElasticityTensor
    block = clad
    temperature = temperature
  []
  [zry_thermal_creep]
    type = ZryCreepLOCAUpdate
    block = clad
    temperature = temperature
    model_irradiation_creep = true
    model_primary_creep = true
    model_thermal_creep = true
    max_inelastic_increment = 5e-4
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = zirlo
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'zry_thermal_creep'
    block = clad
  []
  [clad_irradiation_growth]
    type = ZryIrradiationGrowthEigenstrain
    block = clad
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = zirlo
    eigenstrain_name = clad_irradiation_eigenstrain
  []
  [clad_phase]
    type = ZrPhase
    block = clad
    temperature = temperature
    numerical_method = 2
  []
  [clad_oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 4.18e-03
    clad_outer_radius = 4.75e-03
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = overstrain
    #    effective_strain_rate_creep = creep_strain_rate
    #    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = stress_zz
    hoop_creep_strain = creep_strain_zz
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfract_total
    temperature = temperature
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = pellet
    strain_free_density = ${initial_fuel_density}
  []

[]

[Dampers]
  [limitT]
    type = BoundingValueElementDamper
    min_value = 290.0
    max_value = 3000.0
    variable = temperature
  []
  [limitX]
    type = MaxIncrement
    max_increment = 1e-5
    variable = disp_x
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

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

  line_search = 'none'

  l_max_its = 50
  l_tol = 8e-3
  nl_max_its = 15
  nl_rel_tol = 1e-4
  nl_abs_tol = 1e-8
  # n_startup_steps = 1
  end_time = 95033429.6
  dtmax = 20
  dtmin = 1e-6

  [TimeStepper]
    type = IterationAdaptiveDT
    timestep_limiting_postprocessor = material_timestep
    dt = 10
    optimal_iterations = 20
    iteration_window = 4
    linear_iteration_ratio = 100
    growth_factor = 2
    cutback_factor = .5
    timestep_limiting_function = forced_times
    force_step_every_function_point = true
  []
  [Quadrature]
    order = FIFTH
    side_order = SEVENTH
  []
[]

[Postprocessors]
  [ave_temp_interior]
    type = SideAverageValue
    boundary = 9
    variable = temperature
    execute_on = 'initial linear'
  []
  [clad_inner_vol]
    type = InternalVolume
    boundary = 7
    #outputs = exodus
    execute_on = 'initial timestep_end'
  []
  [fission_gas_produced] # fission gas produced (moles)
    type = ElementIntegralFisGasGeneratedSifgrs
    block = pellet
    execute_on = 'linear'
  []
  [fission_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    block = pellet
    outputs = exodus
    execute_on = 'linear'
  []
  [fission_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    block = pellet
    outputs = exodus
    execute_on = 'linear'
  []
  [flux_from_clad] # area integrated heat flux from the cladding
    type = SideDiffusiveFluxIntegral
    variable = temperature
    boundary = 5
    diffusivity = thermal_conductivity
  []
  [flux_from_fuel] # area integrated heat flux from the fuel
    type = SideDiffusiveFluxIntegral
    variable = temperature
    boundary = 10
    diffusivity = thermal_conductivity
  []

  [rod_total_power]
    type = ElementIntegralPower
    variable = temperature
    burnup_function = burnup
    block = pellet
  []
  [rod_input_power]
    type = FunctionValuePostprocessor
    function = power_history
    scale_factor = 0.1186 # rod height
  []
  [max_fuel_temp]
    type = NodalExtremeValue
    block = pellet
    value_type = max
    variable = temperature
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = clad
    value_type = max
    variable = temperature
  []
  [max_clad_hoop_strain]
    type = ElementExtremeValue
    block = clad
    value_type = max
    variable = strain_zz
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = clad
  []
  [burst]
    type = ElementExtremeValue
    value_type = max
    variable = burst
    block = clad
    execute_on = 'initial timestep_end'
  []
  [he_prod]
    type = IFBAHeProduction
    b10_load = 9.27165354e-5
    b10_enrich = 0.5
    burnup = average_burnup
    zrb2_thick = 10e-6
    fuel_out_rad = 9.32e-3
    ifba_len = 0.3
    u235_enrich = 0.05
  []
  [volume_pulverized]
    type = ElementIntegralMaterialProperty
    mat_prop = pulverized
    block = pellet
  []
  [max_fuel_temp_periphery]
    type = NodalExtremeValue
    value_type = max
    variable = temperature
    boundary = 10
  []
  [additional_volume]
    type = FunctionValuePostprocessor
    function = 8.5e-6
    execute_on = 'initial linear'
  []
  [addition_temperature]
    type = FunctionValuePostprocessor
    function = 300.0
    execute_on = 'initial linear'
  []
  [equilibrium_pressure]
    type = FunctionValuePostprocessor
    function = 101325.0
    execute_on = 'initial linear'
  []
[]

[PerformanceMetricOutputs]
[]

[StandardLWRFuelRodOutputs]
  temperature = temperature
  fuel_pellet_blocks = 3
[]

[Outputs]
  perf_graph = true
  exodus = true
  color = false
  csv = true
  [chkfile]
    type = CSV
    execute_on = FINAL
    show = 'volume_pulverized'
  []
[]

(assessment/LWR/validation/LOCA_Studsvik/analysis/rod_196/Studsvik_196_part1_1p5d_fr_ffrd.i)


initial_fuel_density = 10431.0

[GlobalParams]
  density = ${initial_fuel_density}
  initial_porosity = 0.05
  order = SECOND
  family = LAGRANGE
  energy_per_fission = 3.2e-11 # J/fission
  volumetric_locking_correction = false
  displacements = 'disp_x'
[]

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

[Mesh]
  coord_type = RZ
  [layered1D_mesh]
    type = Layered1DMeshGenerator
    slices_per_block = 10
    clad_gap_width = 80e-6
    plenum_height = 0.0393576
    pellet_outer_radius = 3.92e-3
    clad_thickness = 0.57e-3
    fuel_height = 0.2606424
    # nx_c = 2
    # nx_p = 11
    elem_type = EDGE3
  []
  patch_update_strategy = auto
  patch_size = 10 # For contact algorithm
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[Variables]
  # Define dependent variables and initial conditions
  [temperature]
    initial_condition = 295.0 # set initial temp to coolant inlet
  []
[]

[AuxVariables]
  # Define auxilary variables
  [strain_yy_0]
    order = CONSTANT
    family = MONOMIAL
  []
  [tangential_contact_pressure_aux]
    block = fuel
  []
  [fast_neutron_flux]
    block = clad
  []
  [fast_neutron_fluence]
    block = clad
  []
  [grain_radius]
    block = fuel
    initial_condition = 10e-6
  []
  [creep_strain_rate]
    order = CONSTANT
    family = MONOMIAL
  []
  [effective_creep_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [hoop_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase] # Fraction of beta phase in Zry
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness] # ZrO2 scale thickness (m)
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total] # Current oxigen weight fraction (oxide+metal) (/)
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total] # Gained oxygen weight fraction (oxide+metal) (/)
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress] # Hoop stress at cladding burst
    order = CONSTANT
    family = MONOMIAL
  []
  [burst] # Did cladding burst occur?
    order = CONSTANT
    family = MONOMIAL
  []

  [gap_cond]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    data_file = power_history.csv
    format = columns
    scale_factor = 1
  []
  [axial_peaking_factors]
    type = ParsedFunction
    expression = 1
  []
  [pressure_ramp] # reads and interpolates input data defining amplitude curve for fill gas pressure
    type = PiecewiseLinear
    x = '-200 0  86400  47386400  47472800  47559200  47645600  94945600  95032000'
    y = '0.0065371 1 1 1 1 1 1 1 0.0065371'
    scale_factor = 15.5e6
  []
  [forced_times]
    type = PiecewiseLinear
    data_file = timestep_limiting.csv
    scale_factor = 1
    format = columns
  []
  [clad_axial_pressure]
    type = CladdingAxialPressureFunction
    plenum_pressure = plenum_pressure
    coolant_pressure = pressure_ramp
    coolant_pressure_scaling_factor = 1.0
    fuel_pin_geometry = fuel_pin_geometry
  []
  [fuel_axial_pressure]
    type = ParsedFunction
    expression = plenum_pressure
    symbol_names = plenum_pressure
    symbol_values = plenum_pressure
  []
[]

[Physics]
  [SolidMechanics]
    [Layered1D]
      [fuel]
        block = fuel
        add_variables = true
        add_scalar_variables = true
        strain = FINITE
        out_of_plane_strain_name = strain_yy
        eigenstrain_names = 'fuel_thermal_eigenstrain fuel_volumetric_eigenstrain axial_relocation_eigenstrain'
        generate_output = 'stress_xx stress_yy stress_zz vonmises_stress strain_xx hoop_stress creep_strain_zz strain_zz'
        extra_vector_tags = 'ref'
        fuel_pin_geometry = fuel_pin_geometry
        group_scalar_vars_in_reference_residual = true
        mesh_generator = layered1D_mesh
        decomposition_method = EigenSolution
        temperature = temperature
        out_of_plane_pressure_function = fuel_axial_pressure
        layer_friction_user_object = 1DFriction_secondary
      []
      [clad]
        block = clad
        add_variables = true
        add_scalar_variables = true
        strain = FINITE
        out_of_plane_strain_name = strain_yy
        eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_eigenstrain'
        generate_output = 'stress_xx stress_yy stress_zz vonmises_stress strain_xx hoop_stress creep_strain_zz strain_zz'
        extra_vector_tags = 'ref'
        fuel_pin_geometry = fuel_pin_geometry
        group_scalar_vars_in_reference_residual = true
        mesh_generator = layered1D_mesh
        decomposition_method = EigenSolution
        temperature = temperature
        out_of_plane_pressure_function = clad_axial_pressure
        layer_friction_user_object = 1DFriction_primary
      []
    []
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_source]
    type = NeutronHeatSource
    variable = temperature
    extra_vector_tags = 'ref'
    block = fuel
    burnup_function = burnup
    axial_relocation_object = axial_relocation
  []
[]

[Burnup]
  [burnup]
    block = fuel
    rod_ave_lin_pow = power_history # using the power function defined above
    axial_power_profile = axial_peaking_factors # using the axial power profile function defined above
    num_radial = 80
    num_axial = 11
    fuel_pin_geometry = fuel_pin_geometry
    fuel_volume_ratio = 1.0 # for use with dished fuels (ratio of actual volume to cylinder volume)
    order = CONSTANT
    family = MONOMIAL
    RPF = RPF
    isotopes = 'U235 U238 Pu239 Pu240 Pu241 Pu242'
    isotope_fractions = '0.05 0.95 0 0 0 0'
  []
[]

[AuxKernels]
  # Define auxilliary kernels for each of the aux variables
  [fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    block = clad
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    factor = 3e13
    execute_on = timestep_begin
  []
  [fast_neutron_fluence]
    type = FastNeutronFluenceAux
    variable = fast_neutron_fluence
    block = clad
    fast_neutron_flux = fast_neutron_flux
    execute_on = timestep_begin
  []
  [grain_radius]
    type = GrainRadiusAux
    block = fuel
    variable = grain_radius
    temperature = temperature
    execute_on = linear
  []
  [creep_strain_rate]
    type = MaterialRealAux
    property = creep_rate
    variable = creep_strain_rate
    block = clad
  []
  [effective_creep_strain]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = effective_creep_strain
    block = clad
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    block = clad
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = MaterialRealAux
    boundary = 2
    variable = scale_thickness
    property = oxide_scale_thickness
  []
  [ofract_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
  []
  [ofgain_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    boundary = 2
    variable = burst_stress
    property = burst_stress
  []
  [hasburst]
    type = MaterialRealAux
    boundary = 2
    variable = burst
    property = failed
    execute_on = timestep_end
  []
  [conductance]
    type = MaterialRealAux
    property = gap_conductance
    variable = gap_cond
    boundary = 10
    execute_on = 'linear'
  []
  [tangential_contact_pressure_aux]
    type = SpatialUserObjectAux
    variable = tangential_contact_pressure_aux
    user_object = 1DFriction_secondary
    block = fuel
    execute_on = 'TIMESTEP_END'
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    formulation = kinematic
    model = frictionless
    penalty = 1e7
  []
[]

[ThermalContact]
  [thermal_contact]
    type = GasGapHeatTransfer
    variable = temperature
    primary = 5
    secondary = 10
    initial_moles = initial_moles
    gas_released = 'fission_gas_released he_prod'
    released_gas_types = 'Kr Xe;
                          He'
    released_fractions = '0.153 0.847;
                          1'
    quadrature = true
    contact_pressure = contact_pressure
    refab_gas_types = He
    refab_fractions = 1
    refab_time = 95032000
    refab_type = 0
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 12
    value = 0.0
  []
  [Pressure]
    [coolantPressure]
      boundary = '2'
      function = pressure_ramp
    []
  []
  [PlenumPressure]
    [plenumPressure]
      boundary = 9
      initial_pressure = 3.44738e6
      startup_time = 0
      R = 8.3143
      output_initial_moles = initial_moles
      temperature = plenum_temp
      volume = plenum_volume
      material_input = 'fission_gas_released he_prod'
      output = plenum_pressure
      refab_time = 95032000
      refab_pressure = 8.2e6
      refab_temperature = 295.0
      refab_volume = 1.04e-05
      cladding_failure_status = burst
      equilibrium_pressure = equilibrium_pressure
      additional_volumes = additional_volume
      temperature_of_additional_volumes = addition_temperature
    []
  []
[]

[UserObjects]
  [layered_average_hoop_strain]
    type = LayeredAverage
    block = clad
    num_layers = 10
    direction = y
    variable = strain_zz
  []
  [cladding_strain_yy]
    type = LayeredAverage
    block = clad
    num_layers = 11
    direction = y
    variable = strain_yy
    execute_on = 'initial timestep_end'
  []
  [fuel_strain_yy]
    type = LayeredAverage
    block = fuel
    num_layers = 10
    direction = y
    variable = strain_yy
    execute_on = 'initial timestep_end'
  []
  # [fuel_pin_geometry]
  #   type = Layered1DFuelPinGeometry
  #   mesh_generator = layered1D_mesh
  # []
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  # We could have two element UOs to obtain interface stress
  [1DContactStressOOP_fuel]
    type = Layered1DContactInterfaceStress
    direction = y
    stress_name = stress
    num_layers = 10
    # If we do not provide the numbers below, it will look at the mesh, in all blocks to set the layer number. Then, it will
    # be wrong because the cladding has more height and won't be able to identify layers in the fuel.
    direction_min = 0.01306
    direction_max = 0.24761028

    block = fuel
    execute_on = 'LINEAR NONLINEAR'
  []
  [1DContactStressOOP_cladding]
    type = Layered1DContactInterfaceStress
    direction = y
    stress_name = stress
    num_layers = 10
    # If we do not provide the numbers below, it will look at the mesh, in all blocks to set the layer number. Then, it will
    # be wrong because the cladding has more height and won't be able to identify layers in the fuel.
    direction_min = 0.01306
    direction_max = 0.24761028

    block = clad
    execute_on = 'LINEAR NONLINEAR'
  []
  [1DFriction_secondary]
    type = Layered1DFrictionalForce
    force_postaux = true
    contact_pressure = contact_pressure
    direction = y
    boundary = pellet_outer_radial_surface
    num_layers = 10
    interface_oop_stress_provider_fuel = 1DContactStressOOP_fuel
    interface_oop_stress_provider_cladding = 1DContactStressOOP_cladding
    is_secondary_side = true
    tangential_pressure = tangential_contact_pressure_aux
    friction_coefficient = 0.2
    thickness = 0.02606424
    penalty_factor = 1.0e13
    # If we do not provide the numbers below, it will look at the mesh, in all blocks to set the layer number. Then, it will
    # be wrong because the cladding has more height and won't be able to identify layers in the fuel.
    direction_min = 0.01306
    direction_max = 0.24761028

    scalar_var_name_base_fuel = scalar_strain_yy_fuel
    scalar_num_variable_fuel = 10
    scalar_var_name_base_cladding = scalar_strain_yy_clad
    scalar_num_variable_cladding = 10

    execute_on = 'LINEAR NONLINEAR'
  []
  [1DFriction_primary]
    type = Layered1DFrictionalForce
    force_postaux = true
    contact_pressure = contact_pressure
    direction = y

    boundary = clad_inside_right
    num_layers = 10
    # If we do not provide the numbers below, it will look at the mesh, in all blocks to set the layer number. Then, it will
    # be wrong because the cladding has more height and won't be able to identify layers in the fuel.
    direction_min = 0.0165094
    direction_max = 0.24761028

    interface_oop_stress_provider_fuel = 1DContactStressOOP_fuel
    interface_oop_stress_provider_cladding = 1DContactStressOOP_cladding
    is_secondary_side = false
    secondary_side_frictional_user_object = 1DFriction_secondary
    friction_coefficient = 0.2
    thickness = 0.02606424
    penalty_factor = 1.0e13
    scalar_var_name_base_fuel = scalar_strain_yy_fuel
    scalar_num_variable_fuel = 10
    scalar_var_name_base_cladding = scalar_strain_yy_clad
    scalar_num_variable_cladding = 10

    execute_on = 'LINEAR NONLINEAR'
  []
[]

[PlenumTemperature]
  [plenum_temp]
    boundary = 5
    inner_surfaces = '5'
    outer_surfaces = '10'
    temperature = temperature
  []
[]

[CoolantChannel]
  [convective_clad_surface] # apply convective boundary to clad outer surface
    boundary = 2
    variable = temperature
    inlet_temperature = 580
    inlet_pressure = 15.5e6 # Pa
    inlet_massflux = 3800 # kg/m^2-sec
    rod_diameter = 0.00914 # m
    rod_pitch = 1.26e-2 # m
    compute_enthalpy = false
    linear_heat_rate = power_history
    axial_power_profile = axial_peaking_factors
    output_properties = 'coolant_channel_htype coolant_channel_hmode'
  []
[]

[Materials]
  # [uo2_pulverization]
  #   type = UO2Pulverization
  #   block = fuel
  #   layered_average_contact_pressure = contact_pressure
  #   temperature = temperature
  #   burnup_function = burnup
  #   output_properties = pulverized
  #   outputs = all
  # []
  [fuel_dispersal]
    type = UO2Dispersal
    block = fuel
    axial_relocation_object = axial_relocation
    layered_average_burnup = layered_average_burnup
    layered_average_hoop_strain = layered_average_hoop_strain
    dispersal_model = ONE_MM_TWO_PERCENT_STRAIN
  []

  # Define material behavior models and input material property data
  [fuel_thermal] # temperature and burnup dependent thermal properties of UO2 (BISON kernel)
    type = UO2Thermal
    block = fuel
    thermal_conductivity_model = NFIR
    temperature = temperature
    burnup_function = burnup
    axial_relocation_object = axial_relocation
    gap_thermal_conductivity = layered_average_gap_conductivity
  []
  [fuel_elasticity_tensor]
    type = UO2IsotropicDamageElasticityTensor
    block = fuel
    fragmentation_model = BARANI
    temperature = temperature
    rod_ave_lin_pow = power_history
    axial_relocation_object = axial_relocation
  []
  [fuel_elastic_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'fuel_creep'
    block = fuel
  []
  [fuel_creep]
    type = UO2CreepUpdate
    block = fuel
    temperature = temperature
    fission_rate = fission_rate
    initial_grain_radius = 10.0e-6
    oxygen_to_metal_ratio = 2.0
  []
  # [fuel_relocation]
  #   type = UO2RelocationEigenstrain
  #   block = fuel
  #   burnup_function = burnup
  #   fuel_pin_geometry = fuel_pin_geometry
  #   rod_ave_lin_pow = power_history
  #   axial_power_profile = axial_peaking_factors
  #   burnup_relocation_stop = 0.024
  #   relocation_activation1 = 5000
  #   relocation_model = ESCORE_modified
  #   eigenstrain_name = fuel_relocation_eigenstrain
  # []
  [fuel_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = fuel
    thermal_expansion_coeff = 10.0e-6
    temperature = temperature
    stress_free_temperature = 295.0
    eigenstrain_name = fuel_thermal_eigenstrain
  []
  [fuel_volumetric_swelling]
    type = UO2VolumetricSwellingEigenstrain
    gas_swelling_model_type = SIFGRS
    block = fuel
    temperature = temperature
    burnup_function = burnup
    initial_fuel_density = 10431.0
    eigenstrain_name = fuel_volumetric_eigenstrain
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = fuel
    temperature = temperature
    burnup_function = burnup
    grain_radius = grain_radius
    gbs_model = true
  []

  [clad_density]
    type = StrainAdjustedDensity
    block = clad
    strain_free_density = 6550.
  []
  [clad_thermal]
    block = clad
    type = ZryThermal
    temperature = temperature
  []
  [clad_thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = clad
    temperature = temperature
    stress_free_temperature = 295.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [clad_elasticity_tensor]
    type = ZryElasticityTensor
    block = clad
    temperature = temperature
  []
  [zry_thermal_creep]
    type = ZryCreepLOCAUpdate
    block = clad
    temperature = temperature
    model_irradiation_creep = true
    model_primary_creep = true
    model_thermal_creep = true
    max_inelastic_increment = 5e-4
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = zirlo
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'zry_thermal_creep'
    block = clad
  []
  [clad_irradiation_growth]
    type = ZryIrradiationGrowthEigenstrain
    block = clad
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = zirlo
    eigenstrain_name = clad_irradiation_eigenstrain
  []
  [clad_phase]
    type = ZrPhase
    block = clad
    temperature = temperature
    numerical_method = 2
  []
  [clad_oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 4.18e-03
    clad_outer_radius = 4.75e-03
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = overstrain
    #    effective_strain_rate_creep = creep_strain_rate
    #    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    hoop_creep_strain = creep_strain_zz
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfract_total
    temperature = temperature
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = fuel
    strain_free_density = ${initial_fuel_density}
  []
[]

[VectorPostprocessors]
  [cladding_outer]
    type = NodalValueSampler
    boundary = 5
    variable = disp_x
    sort_by = y
  []
[]

[AxialRelocation]
  [relocation]
    rod_ave_lin_pow = power_history
    axial_direction = y
    fuel_blocks = fuel
    clad_blocks = clad
    contact_pressure_variable = contact_pressure
    out_of_plane_strain_variable = strain_yy_0
    penetration_variable = penetration
    clad_inner_volume_addition = 0
    burnup_variable = burnup
    temperature = temperature
    axial_relocation_output_options = MASS_FRACTION
    mesh_generator = layered1D_mesh
    # CHANGE
    gap_thickness_threshold = 0.000050
  []
[]

[Postprocessors]
  [volume_fuel_dispersed]
    type = LayeredElementIntegralMaterialProperty
    block = fuel
    mat_prop = dispersed
    fuel_pin_geometry = fuel_pin_geometry
    execute_on = 'initial timestep_end'
  []
  [mass_fuel_dispersed]
    type = ParsedPostprocessor
    pp_names = volume_fuel_dispersed
    expression = '10431 * volume_fuel_dispersed'
    execute_on = 'initial timestep_end'
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

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

  line_search = 'none'

  l_max_its = 50
  l_tol = 8e-3
  nl_max_its = 15
  nl_rel_tol = 1e-4
  nl_abs_tol = 1e-8
  start_time = -10
  n_startup_steps = 1
  end_time = 95032000
  dtmax = 1e6
  dtmin = 1e-6

  [TimeStepper]
    type = IterationAdaptiveDT
    timestep_limiting_postprocessor = material_timestep
    dt = 10
    optimal_iterations = 20
    iteration_window = 4
    linear_iteration_ratio = 100
    growth_factor = 2
    cutback_factor = .5
    timestep_limiting_function = forced_times
    force_step_every_function_point = true
  []
  [Quadrature]
    order = FIFTH
    side_order = SEVENTH
  []
[]

[Postprocessors]
  [ave_temp_interior]
    type = SideAverageValue
    boundary = 9
    variable = temperature
    execute_on = 'initial linear'
  []
  [clad_inner_vol]
    type = InternalVolume
    boundary = 7
    #outputs = exodus
    execute_on = 'initial timestep_end'
  []
  [fission_gas_produced] # fission gas produced (moles)
    type = ElementIntegralFisGasGeneratedSifgrs
    block = fuel
    execute_on = 'linear'
  []
  [fission_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    block = fuel
    outputs = exodus
    execute_on = 'linear'
  []
  [fission_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    block = fuel
    outputs = exodus
    execute_on = 'linear'
  []
  [flux_from_clad] # area integrated heat flux from the cladding
    type = SideDiffusiveFluxIntegral
    variable = temperature
    boundary = 5
    diffusivity = thermal_conductivity
  []
  [flux_from_fuel] # area integrated heat flux from the fuel
    type = SideDiffusiveFluxIntegral
    variable = temperature
    boundary = 10
    diffusivity = thermal_conductivity
  []

  [rod_total_power]
    type = ElementIntegralPower
    variable = temperature
    burnup_function = burnup
    block = fuel
  []
  [max_fuel_temp]
    type = NodalExtremeValue
    block = fuel
    value_type = max
    variable = temperature
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = clad
    value_type = max
    variable = temperature
  []
  [max_clad_hoop_strain]
    type = ElementExtremeValue
    block = clad
    value_type = max
    variable = strain_zz
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = clad
  []
  [burst]
    type = ElementExtremeValue
    value_type = max
    variable = burst
    block = clad
    execute_on = 'initial timestep_end'
  []
  [he_prod]
    type = IFBAHeProduction
    b10_load = 9.27165354e-5
    b10_enrich = 0.5
    burnup = average_burnup
    zrb2_thick = 10e-6
    fuel_out_rad = 9.32e-3
    ifba_len = 0.3
    u235_enrich = 0.05
  []
  [volume_pulverized]
    type = ElementIntegralMaterialProperty
    mat_prop = pulverized
    block = fuel
  []
  [max_fuel_temp_periphery]
    type = NodalExtremeValue
    value_type = max
    variable = temperature
    boundary = 10
  []
  [additional_volume]
    type = FunctionValuePostprocessor
    function = 8.5e-6
    execute_on = 'initial linear'
  []
  [addition_temperature]
    type = FunctionValuePostprocessor
    function = 300.0
    execute_on = 'initial linear'
  []
  [equilibrium_pressure]
    type = FunctionValuePostprocessor
    function = 101325.0
    execute_on = 'initial linear'
  []
[]

[PerformanceMetricOutputs]
[]

[StandardLWRFuelRodOutputs]
  layered = true
  fuel_pin_geometry = fuel_pin_geometry
  fuel_pellet_blocks = 'fuel'
[]

[Outputs]
  perf_graph = true
  exodus = true
  color = false
  csv = true
  [checkpoint]
    type = Checkpoint
    num_files = 2
  []
  [chkfile]
    type = CSV
    execute_on = FINAL
    show = 'volume_pulverized'
  []
[]

(test/tests/ifba_he_production/fill_gas_xenon_w_ifba.i)

#
# 2-D RZ One Pellet Test - IFBA using Xenon as fill gas
#
# This test is of a single pellet with cladding and a specified initial
# pressure of Xe fill gas. In addition, an IFBA layer is added which will
# generate He gas to be added to the plenum. The postprocessor interior_temp
# should be the same as the pure Xe test case initially and as the He gas is
# added to the plenum from the IFBA, the interior_temp value should approach
# the He fill gas test case (both in the doc subdirectory).
#
# This model demonstrates that the gas conductance for the plenum is being
# updated for the He gas generated by the IFBA layer.
#

initial_fuel_density = 10431.0 #95% TD (TD = 10980)

[GlobalParams]
  density = ${initial_fuel_density}
  order = SECOND
  family = LAGRANGE
  energy_per_fission = 3.2e-11 # J/fission (205 Mev)
  displacements = 'disp_x disp_y'
  temperature = temp
[]

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

[Mesh]
  coord_type = RZ
  [smeared_pellet_mesh]
    type = FuelPinMeshGenerator
    clad_mesh_density = customize
    pellet_mesh_density = customize
    ny_p = 1
    nx_p = 1
    nx_c = 1
    ny_cu = 1
    ny_c = 1
    ny_cl = 1
    clad_thickness = 5.6e-4
    pellet_outer_radius = 0.0041
    pellet_height = 0.01
    pellet_quantity = 1
    clad_bot_gap_height = 1e-3
    bottom_clad_height = 2.24e-3
    top_clad_height = 2.24e-3
    clad_gap_width = 8e-5
    plenum_fuel_ratio = 0.150
    elem_type = QUAD8
  []
  partitioner = centroid
  centroid_partitioner_direction = y
  patch_size = 5
[]

[Variables]
  [temp]
    initial_condition = 298
  []
[]

[AuxVariables]
  [fission_rate]
    block = '3'
  []
  [burnup]
    block = '3'
  []
  [grain_radius]
    block = '3'
    initial_condition = 5e-6
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 1.0e4 1.0e8'
    y = '0  1.0  1.0'
    scale_factor = 20e3 # 20 kW/m peak power.
  []
  [coolant_pressure_ramp]
    type = PiecewiseLinear
    x = '0 10000'
    y = '0 1'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [fuel]
    block = 3
    strain = FINITE
    incremental = true
    add_variables = true
    decomposition_method = EigenSolution
    extra_vector_tags = 'ref'
    eigenstrain_names = fuel_thermal_strain
  []
  [clad]
    block = 1
    strain = FINITE
    incremental = true
    add_variables = true
    decomposition_method = EigenSolution
    extra_vector_tags = 'ref'
    eigenstrain_names = clad_thermal_strain
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_source_fuel]
    type = NeutronHeatSource
    variable = temp
    block = '3'
    fission_rate = fission_rate
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [fissionrate]
    type = FissionRateGeneral
    fission_rate_formulation = GENERIC
    variable = fission_rate
    block = '3'
    value = 5.3548e+14
    fission_rate_function = power_history
  []
  [burnup]
    type = BurnupAux
    variable = burnup
    block = '3'
    fission_rate = fission_rate
    molecular_weight = 0.270
  []
  [grain_radius]
    type = GrainRadiusAux
    block = '3'
    variable = grain_radius
    temperature = temp
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    penalty = 1e+14 #1e7
    model = frictionless
    tangential_tolerance = 5e-4
    normal_smoothing_distance = 0.1
    normalize_penalty = true
  []
[]

[ThermalContact]
  [pellet_clad_thermal]
    type = GasGapHeatTransfer
    variable = temp
    primary = 5
    secondary = 10
    gas_released = 'fis_gas_released he_prod'
    initial_moles = initial_moles
    jump_distance_model = LANNING
    layer_thickness = layer_thickness
    plenum_pressure = plenum_pressure
    contact_pressure = contact_pressure
    initial_gas_types = Xe
    initial_fractions = 1
    released_gas_types = 'Kr Xe;
                          He'
    released_fractions = '0.153 0.847;
                          1'
    roughness_coef = 3.2
    roughness_secondary = 1e-6
    roughness_primary = 2e-6
    emissivity_primary = 0.8
    emissivity_secondary = 0.8
    quadrature = true
    normal_smoothing_distance = 0.1
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 1003
    value = 0.0
  []
  [no_y_clad_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 1
    value = 0.0
  []
  [no_y_fuel_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 20
    value = 0.0
  []
  [no_x_fuel]
    type = DirichletBC
    variable = disp_x
    boundary = 1005
    value = 0.0
  []
  [Clad_Temp]
    type = DirichletBC
    variable = temp
    boundary = '2'
    value = 580.0
  []
  [Pressure]
    [coolantPressure]
      boundary = '2'
      factor = 15.5e6
      function = coolant_pressure_ramp
    []
  []
  [PlenumPressure]
    [plenumPressure]
      boundary = 9
      initial_pressure = 0.50e6
      startup_time = 0.0
      material_input = 'fis_gas_released he_prod'
      output_initial_moles = initial_moles
      temperature = interior_temp
      volume = gas_volume
      output = plenum_pressure
    []
  []
[]

[Materials]
  [fuel_thermal]
    type = UO2Thermal
    block = '3'
    temperature = temp
    burnup = burnup
    thermal_conductivity_model = NFIR
  []
  [fuel_elasticity_tensor]
    type = UO2ElasticityTensor
    block = 3
  []
  [fuel_elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 3
  []
  [fuel_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = 3
    thermal_expansion_coeff = 10.0e-6
    stress_free_temperature = 298
    eigenstrain_name = fuel_thermal_strain
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = '3'
    temperature = temp
    fission_rate = fission_rate
    grain_radius = grain_radius
    gbs_model = true
    burnup = burnup
    diff_coeff_option = TURNBULL_D1_D2
  []
  [clad_thermal]
    type = HeatConductionMaterial
    block = 1
    thermal_conductivity = 16.0
    specific_heat = 330.0
  []
  [fclad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = 1
    youngs_modulus = 7.5e10
    poissons_ratio = 0.3
  []
  [clad_elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 1
  []
  [clad_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = 1
    thermal_expansion_coeff = 5.0e-6
    stress_free_temperature = 298
    eigenstrain_name = clad_thermal_strain
  []
  [clad_density]
    type = StrainAdjustedDensity
    block = 1
    strain_free_density = 6551.0
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = '3'
    strain_free_density = ${initial_fuel_density}
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

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

  line_search = 'none'

  l_max_its = 25
  nl_max_its = 40
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-8

  dtmax = 1.0e6
  dtmin = 1.0
  end_time = 5.3e7 # 1.7 years (~3% burnup)

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e3
    optimal_iterations = 30
    iteration_window = 4
    time_t = '0   1e4 1e8'
    time_dt = '1e4 1e5 1e6'
    timestep_limiting_function = power_history
    force_step_every_function_point = true
  []
  [Quadrature]
    order = fifth
    side_order = seventh
  []
  verbose = true
[]

[Postprocessors]
  [clad_inner_vol]
    type = InternalVolume
    boundary = 7
    execute_on = 'initial linear'
  []
  [pellet_volume]
    type = InternalVolume
    boundary = 8
    execute_on = 'initial linear'
  []
  [gas_volume]
    type = InternalVolume
    boundary = 9
    execute_on = 'initial linear'
  []
  [interior_temp]
    type = SideAverageValue
    boundary = 9 # cladding interior and pellet exterior
    variable = temp
    execute_on = 'initial linear'
  []
  [fis_gas_produced] # fission gas produced (moles)
    type = ElementIntegralFisGasGeneratedSifgrs
    block = '3'
  []
  [fis_gas_released]
    type = ElementIntegralFisGasReleasedSifgrs
    block = '3'
  []
  [fis_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    block = '3'
  []
  [fis_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    block = '3'
  []
  [power_history]
    type = FunctionValuePostprocessor
    function = power_history
  []
  [dt]
    type = TimestepSize
  []
  [residual]
    type = Residual
  []
  [nl_its]
    type = NumNonlinearIterations
  []
  [lin_its]
    type = NumLinearIterations
  []
  [average_burnup]
    type = ElementAverageValue
    block = '3'
    variable = burnup
  []
  [burnup]
    type = ElementAverageValue
    block = '3'
    variable = burnup
  []
  [average_fissionrate]
    type = ElementAverageValue
    block = '3'
    variable = fission_rate
  []
  [rod_total_power]
    type = ElementIntegralPower
    variable = temp
    fission_rate = fission_rate
    block = '3'
  []
  [rod_input_power]
    type = FunctionValuePostprocessor
    function = power_history
    scale_factor = 0.01 # change: length of fuel stack in meters (1 pellet height)
  []
  [he_prod]
    type = IFBAHeProduction
    zrb2_load = 1.181e-4
    ifba_len = 1.0e-2
    b10_enrich = 0.50
    zrb2_rel_dens = 0.7
    model = burnup
    u235_enrich = 0.045
    burnup = average_burnup
  []
[]

[Outputs]
  time_step_interval =  1
  exodus = false
  [console]
    type = Console
    solve_log = true
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    show = 'average_burnup burnup he_prod interior_temp plenum_pressure'
    file_base = fill_gas_xenon_w_ifba_check
  []
  [out]
    type = CSV
    delimiter = ' '
  []
[]

(test/tests/ifba_he_production/ifba_examp_template.i)

#
# 2-D RZ One Pellet Test - Coarse mesh example of IFBA layer
#
# This is an input template for a fast running example using the IFBA
# postprocessor. All of the possible ways to specify the IFBA layer are run
# using this template in a regression test format.
#
# The expected ouputs for each test depends on the model equation being used
# to calculate the He produced. For the burnup based equation, the He moles
# released at the end of the calculation is 1.4897e-6. A hand calculation is
# reproduced in the Excel spreadsheet IFBA_He_Calc included in the test
# directory. The burnup equation result computed for the same inputs is
# 1.4902e-6.
#
# Using the FRAPCON equation calculates a rate of He production, so comparing
# the first couple of time steps of the simulation to the hand calculation is
# more straightforward. Comparing the BISON results to the hand calculation is
#
#   Time(s)       He Prod (BISON)       He Prod (Excel)
#   1000            1.01465e-10           1.01465e-10
#   3000            7.10250e-10           7.18769e-10
#

initial_fuel_density = 10431.0 #95% TD (TD = 10980)

[GlobalParams]
  density = ${initial_fuel_density}
  displacements = 'disp_x disp_y'
  order = SECOND
  energy_per_fission = 3.2e-11 # J/fission (205 Mev)
  temperature = temp
  volumetric_locking_correction = false
[]

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

[Mesh]
  coord_type = RZ
  [smeared_pellet_mesh]
    type = FuelPinMeshGenerator
    clad_mesh_density = customize
    pellet_mesh_density = customize
    ny_p = 1
    nx_p = 1
    nx_c = 1
    ny_cu = 1
    ny_c = 1
    ny_cl = 1
    clad_thickness = 5.6e-4
    pellet_outer_radius = 0.0041
    pellet_height = 0.01
    pellet_quantity = 1
    clad_bot_gap_height = 1e-3
    bottom_clad_height = 2.24e-3
    top_clad_height = 2.24e-3
    clad_gap_width = 8e-5
    plenum_fuel_ratio = 0.150
    elem_type = QUAD8
  []
  partitioner = centroid
  centroid_partitioner_direction = y
  patch_size = 5
[]

[Variables]
  [temp]
    initial_condition = 298
  []
[]

[AuxVariables]
  [fission_rate]
    block = '3'
  []
  [burnup]
    block = '3'
  []
  [grain_radius]
    block = '3'
    initial_condition = 5e-6
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 1.0e4 1.0e8'
    y = '0  1.0  1.0'
    scale_factor = 20e3 # 20 kW/m peak power.
  []
  [coolant_pressure_ramp]
    type = PiecewiseLinear
    x = '0 10000'
    y = '0 1'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [fuel]
    block = 3
    strain = FINITE
    incremental = true
    add_variables = true
    decomposition_method = EigenSolution
    extra_vector_tags = 'ref'
    eigenstrain_names = 'fuel_thermal_strain'
  []
  [clad]
    block = 1
    strain = FINITE
    incremental = true
    add_variables = true
    decomposition_method = EigenSolution
    extra_vector_tags = 'ref'
    eigenstrain_names = 'clad_thermal_strain'
  []

[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_source_fuel]
    type = NeutronHeatSource
    variable = temp
    block = '3'
    fission_rate = fission_rate
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [fissionrate]
    type = FissionRateGeneral
    fission_rate_formulation = GENERIC
    variable = fission_rate
    block = '3'
    value = 5.3548e+14
    fission_rate_function = power_history
  []
  [burnup]
    type = BurnupAux
    variable = burnup
    block = '3'
    fission_rate = fission_rate
    molecular_weight = 0.270
  []
  [grain_radius]
    type = GrainRadiusAux
    block = '3'
    variable = grain_radius
    temperature = temp
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    penalty = 1e+14 #1e7
    model = frictionless
    tangential_tolerance = 5e-4
    normal_smoothing_distance = 0.1
    normalize_penalty = true
  []
[]

[ThermalContact]
  [pellet_clad_thermal]
    type = GasGapHeatTransfer
    variable = temp
    primary = 5
    secondary = 10
    gas_released = 'fis_gas_released he_prod'
    initial_moles = initial_moles
    jump_distance_model = LANNING
    layer_thickness = layer_thickness
    plenum_pressure = plenum_pressure
    contact_pressure = contact_pressure
    released_gas_types = 'Kr Xe;
                          He'
    released_fractions = '0.153 0.847;
                          1'
    roughness_coef = 3.2
    roughness_secondary = 1e-6
    roughness_primary = 2e-6
    emissivity_primary = 0.8
    emissivity_secondary = 0.8
    quadrature = true
    normal_smoothing_distance = 0.1
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 1003
    value = 0.0
  []
  [no_y_clad_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 1
    value = 0.0
  []
  [no_y_fuel_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 20
    value = 0.0
  []
  [no_x_fuel]
    type = DirichletBC
    variable = disp_x
    boundary = 1005
    value = 0.0
  []
  [Clad_Temp]
    type = DirichletBC
    variable = temp
    boundary = '2'
    value = 580.0
  []
  [Pressure]
    [coolantPressure]
      boundary = '2'
      factor = 15.5e6
      function = coolant_pressure_ramp
    []
  []
  [PlenumPressure]
    [plenumPressure]
      boundary = 9
      initial_pressure = 0.50e6
      refab_pressure = 0.50e6
      startup_time = 0.0
      material_input = 'fis_gas_released he_prod'
      output_initial_moles = initial_moles
      temperature = interior_temp
      volume = gas_volume
      output = plenum_pressure
      displacements = 'disp_x disp_y'
    []
  []
[]

[Materials]
  [fuel_thermal]
    type = UO2Thermal
    block = '3'
    temperature = temp
    burnup = burnup
    thermal_conductivity_model = NFIR
  []
  [fuel_elasticity_tensor]
    type = UO2ElasticityTensor
    block = 3
  []
  [fuel_elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 3
  []
  [fuel_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = 3
    thermal_expansion_coeff = 10.0e-6
    stress_free_temperature = 298
    eigenstrain_name = 'fuel_thermal_strain'
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = '3'
    temperature = temp
    fission_rate = fission_rate
    grain_radius = grain_radius
    gbs_model = true
    burnup = burnup
    diff_coeff_option = TURNBULL_D1_D2
  []
  [clad_thermal]
    type = HeatConductionMaterial
    block = 1
    thermal_conductivity = 16.0
    specific_heat = 330.0
  []
  [fclad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = 1
    youngs_modulus = 7.5e10
    poissons_ratio = 0.3
  []
  [clad_elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 1
  []
  [clad_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = 1
    thermal_expansion_coeff = 5.0e-6
    stress_free_temperature = 298
    eigenstrain_name = 'clad_thermal_strain'
  []
  [clad_density]
    type = StrainAdjustedDensity
    block = 1
    strain_free_density = 6551.0
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = '3'
    strain_free_density = ${initial_fuel_density}
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

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

  line_search = 'none'

  l_max_its = 25
  nl_max_its = 40
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-8

  dtmax = 1.0e6
  dtmin = 1.0
  end_time = 2.5e6

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e3
    optimal_iterations = 30
    iteration_window = 4
    time_t = '0   1e4 1e8'
    time_dt = '1e4 1e5 1e6'
    timestep_limiting_function = power_history
    force_step_every_function_point = true
  []
  [Quadrature]
    order = fifth
    side_order = seventh
  []
  verbose = true
[]

[Postprocessors]
  [clad_inner_vol]
    type = InternalVolume
    boundary = 7
    execute_on = 'initial linear'
  []
  [pellet_volume]
    type = InternalVolume
    boundary = 8
    execute_on = 'initial linear'
  []
  [gas_volume]
    type = InternalVolume
    boundary = 9
    execute_on = 'initial linear'
  []
  [interior_temp]
    type = SideAverageValue
    boundary = 9 # cladding interior and pellet exterior
    variable = temp
    execute_on = 'initial linear'
  []
  [fis_gas_produced] # fission gas produced (moles)
    type = ElementIntegralFisGasGeneratedSifgrs
    block = '3'
  []
  [fis_gas_released]
    type = ElementIntegralFisGasReleasedSifgrs
    block = '3'
  []
  [fis_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    block = '3'
  []
  [fis_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    block = '3'
  []
  [power_history]
    type = FunctionValuePostprocessor
    function = power_history
  []
  [dt]
    type = TimestepSize
  []
  [residual]
    type = Residual
  []
  [nl_its]
    type = NumNonlinearIterations
  []
  [lin_its]
    type = NumLinearIterations
  []
  [average_burnup]
    type = ElementAverageValue
    block = '3'
    variable = burnup
  []
  [burnup]
    type = ElementAverageValue
    block = '3'
    variable = burnup
  []
  [average_fissionrate]
    type = ElementAverageValue
    block = '3'
    variable = fission_rate
  []
  [rod_total_power]
    type = ElementIntegralPower
    variable = temp
    fission_rate = fission_rate
    block = '3'
  []
  [rod_input_power]
    type = FunctionValuePostprocessor
    function = power_history
    scale_factor = 0.01 #BWR change: length of fuel stack in meters (5*pellet height)
  []
  [he_prod]
    type = IFBAHeProduction
  []
[]

[Outputs]
  time_step_interval =  1
  exodus = false
  [console]
    type = Console
    solve_log = true
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    show = 'average_burnup burnup he_prod interior_temp plenum_pressure'
  []
  [outfile]
    type = CSV
    delimiter = ' '
  []
[]

(test/tests/ifba_he_production/ifba_only_template.i)

#
# 2-D RZ One Pellet IFBA Test - IFBA He generation only
#
# This test is of a single pellet with cladding and a specified initial
# pressure of He fill gas and IFBA layer .
# The initial loading of B-10 is converted to He gas and adds to the
# plenum pressure. No fission gas production is included in this model. This
# allows the effect of the IFBA layer to be seen clearly.
#
# The power is ramped up and held constant to heat the fill gas and establish
# an initial "hot" pressure. Since there is no fission gas production or
# release in this model, the pressure at temperature should be able to be
# calculated and compared to the BISON result.
#
# This case builds on the baseline case. The amount of He added due to IFBA
# can be calculated and the BISON result checked.
#
# This input template is used for a set of tests exercising the main input
# options for the IFBA postprocessor.
#

initial_fuel_density = 10431.0 #95% TD (TD = 10980)

[GlobalParams]
  density = ${initial_fuel_density}
  displacements = 'disp_x disp_y'
  order = SECOND
  energy_per_fission = 3.2e-11 # J/fission (205 Mev)
  temperature = temp
  volumetric_locking_correction = false
[]

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

[Mesh]
  coord_type = RZ
  [smeared_pellet_mesh]
    type = FuelPinMeshGenerator
    clad_mesh_density = customize
    pellet_mesh_density = customize
    ny_p = 1
    nx_p = 1
    nx_c = 1
    ny_cu = 1
    ny_c = 1
    ny_cl = 1
    clad_thickness = 5.6e-4
    pellet_outer_radius = 0.0041
    pellet_height = 0.01
    pellet_quantity = 1
    clad_bot_gap_height = 1e-3
    bottom_clad_height = 2.24e-3
    top_clad_height = 2.24e-3
    clad_gap_width = 8e-5
    plenum_fuel_ratio = 0.150
    elem_type = QUAD8
  []
  partitioner = centroid
  centroid_partitioner_direction = y
  patch_size = 5
[]

[Variables]
  [temp]
    initial_condition = 298
  []
[]

[AuxVariables]
  [fission_rate]
    block = '3'
  []
  [burnup]
    block = '3'
  []
  [grain_radius]
    block = '3'
    initial_condition = 5e-6
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 1.0e4 1.0e8'
    y = '0  1.0  1.0'
    scale_factor = 25e3 # 25 kW/m peak power.
  []
  [coolant_pressure_ramp]
    type = PiecewiseLinear
    x = '0 10000'
    y = '0 1'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [fuel]
    block = 3
    strain = FINITE
    incremental = true
    add_variables = true
    decomposition_method = EigenSolution
    extra_vector_tags = 'ref'
    eigenstrain_names = 'fuel_thermal_strain'
  []
  [clad]
    block = 1
    strain = FINITE
    incremental = true
    add_variables = true
    decomposition_method = EigenSolution
    extra_vector_tags = 'ref'
    eigenstrain_names = 'clad_thermal_strain'
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_source_fuel]
    type = NeutronHeatSource
    variable = temp
    block = '3'
    fission_rate = fission_rate
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [fissionrate]
    type = FissionRateGeneral
    fission_rate_formulation = GENERIC
    variable = fission_rate
    block = '3'
    value = 5.3548e+14
    fission_rate_function = power_history
  []
  [burnup]
    type = BurnupAux
    variable = burnup
    block = '3'
    fission_rate = fission_rate
    molecular_weight = 0.270
  []
  [grain_radius]
    type = GrainRadiusAux
    block = '3'
    variable = grain_radius
    temperature = temp
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    penalty = 1e+14 #1e7
    model = frictionless
    tangential_tolerance = 5e-4
    normal_smoothing_distance = 0.1
    normalize_penalty = true
  []
[]

[ThermalContact]
  [pellet_clad_thermal]
    type = GasGapHeatTransfer
    variable = temp
    primary = 5
    secondary = 10
    initial_moles = initial_moles
    jump_distance_model = LANNING
    layer_thickness = layer_thickness
    plenum_pressure = plenum_pressure
    contact_pressure = contact_pressure
    roughness_coef = 3.2
    roughness_secondary = 1e-6
    roughness_primary = 2e-6
    emissivity_primary = 0.8
    emissivity_secondary = 0.8
    quadrature = true
    normal_smoothing_distance = 0.1
    gas_released = null
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 1003
    value = 0.0
  []
  [no_y_clad_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 1
    value = 0.0
  []
  [no_y_fuel_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 20
    value = 0.0
  []
  [no_x_fuel]
    type = DirichletBC
    variable = disp_x
    boundary = 1005
    value = 0.0
  []
  [Clad_Temp]
    type = DirichletBC
    variable = temp
    boundary = '2'
    value = 580.0
  []
  [Pressure]
    [coolantPressure]
      boundary = '2'
      factor = 15.5e6
      function = coolant_pressure_ramp
    []
  []
  [PlenumPressure]
    [plenumPressure]
      boundary = 9
      initial_pressure = 0.50e6
      startup_time = 0.0
      material_input = he_prod
      output_initial_moles = initial_moles
      temperature = interior_temp
      volume = gas_volume
      output = plenum_pressure
      displacements = 'disp_x disp_y'
    []
  []
[]

[Materials]
  [fuel_thermal]
    type = UO2Thermal
    block = '3'
    temperature = temp
    burnup = burnup
    thermal_conductivity_model = NFIR
  []
  [fuel_elasticity_tensor]
    type = UO2ElasticityTensor
    block = 3
  []
  [fuel_elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 3
  []
  [fuel_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = 3
    thermal_expansion_coeff = 10.0e-6
    stress_free_temperature = 298
    eigenstrain_name = 'fuel_thermal_strain'
  []
  [clad_thermal]
    type = HeatConductionMaterial
    block = 1
    thermal_conductivity = 16.0
    specific_heat = 330.0
  []
  [fclad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = 1
    youngs_modulus = 7.5e10
    poissons_ratio = 0.3
  []
  [clad_elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 1
  []
  [clad_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = 1
    thermal_expansion_coeff = 5.0e-6
    stress_free_temperature = 298
    eigenstrain_name = 'clad_thermal_strain'
  []
  [clad_density]
    type = StrainAdjustedDensity
    block = 1
    strain_free_density = 6551.0
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = '3'
    strain_free_density = ${initial_fuel_density}
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

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

  line_search = 'none'

  l_max_its = 25
  nl_max_its = 40
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-8

  dtmax = 1.0e6
  dtmin = 1.0
  end_time = 5.3e7 # 1.7 years (~3% burnup)

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e3
    optimal_iterations = 30
    iteration_window = 4
    time_t = '0   1e4 1e8'
    time_dt = '1e4 1e6 1e6'
    timestep_limiting_function = power_history
    force_step_every_function_point = true
  []
  [Quadrature]
    order = fifth
    side_order = seventh
  []
  verbose = true
[]

[Postprocessors]
  [clad_inner_vol]
    type = InternalVolume
    boundary = 7
    execute_on = 'initial linear'
  []
  [pellet_volume]
    type = InternalVolume
    boundary = 8
    execute_on = 'initial linear'
  []
  [gas_volume]
    type = InternalVolume
    boundary = 9
    execute_on = 'initial linear'
  []
  [interior_temp]
    type = SideAverageValue
    boundary = 9 # cladding interior and pellet exterior
    variable = temp
    execute_on = 'initial linear'
  []
  [power_history]
    type = FunctionValuePostprocessor
    function = power_history
  []
  [dt]
    type = TimestepSize
  []
  [residual]
    type = Residual
  []
  [nl_its]
    type = NumNonlinearIterations
  []
  [lin_its]
    type = NumLinearIterations
  []
  [average_burnup]
    type = ElementAverageValue
    block = '3'
    variable = burnup
  []
  [burnup]
    type = ElementAverageValue
    block = '3'
    variable = burnup
  []
  [average_fissionrate]
    type = ElementAverageValue
    block = '3'
    variable = fission_rate
  []
  [rod_total_power]
    type = ElementIntegralPower
    variable = temp
    fission_rate = fission_rate
    block = '3'
  []
  [rod_input_power]
    type = FunctionValuePostprocessor
    function = power_history
    scale_factor = 0.05175 #BWR change: length of fuel stack in meters (5*pellet height)
  []
  [he_prod]
    type = IFBAHeProduction
  []
  [null]
    type = FunctionValuePostprocessor
    function = 0
  []
[]

[Outputs]
  time_step_interval =  1
  exodus = false
  [console]
    type = Console
    solve_log = true
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    show = 'average_burnup burnup he_prod interior_temp plenum_pressure'
  []
  [out]
    type = CSV
    delimiter = ' '
  []
[]

Description
Helium Production
Input Parameter Conventions
Example Input Syntax
Input Parameters
Input Files
References