UO2ThermalExpansionMATPROEigenstrain

Computes eigenstrain due to thermal expansion in UO2 fuel using MATPRO correlations.

Description

This class, UO2ThermalExpansionMATPROEigenstrain, calculates the strain due to thermal expansion using correlations from MATPRO. The correlation in MATPRO is a function of Plutonium content in the fuel, the fraction of the fuel that is in a molten state, and the temperature of the fuel. These correlations are based on a reference temperature of 300K, but this class permits the use of arbitrary values for the stress-free temperature. In BISON we assume a Pu content of 0.0 and molten fuel fraction of 0.0 in computing the thermal expansion of the UO $var element = document.getElementById("moose-equation-64bf42ec-4d55-42d2-a774-99ec8d81d78f");katex.render("_2", element, {displayMode:false,throwOnError:false});$ fuel.

This material uses the MATPRO fthexp function to determine the thermal expansion strains based on temperature (Allison et al., 1993).

An alternative model for thermal expansion in UO $var element = document.getElementById("moose-equation-796f5db1-2453-4fd2-adbf-2ec4cefe92e7");katex.render("_2", element, {displayMode:false,throwOnError:false});$ fuel proposed by Martin (1988) is also available in BISON, and is provided in UO2ThermalExpansionMartinEigenstrain. The MATPRO and Martin correlations produce very similar results at lower temperatures, but the Martin model predicts higher thermal strains at higher temperatures (above about 2200 K). Figure 1 shows a comparison of the free thermal expansion predicted by the two models.

Figure 1: Comparison of linear free expansion of UO $var element = document.getElementById("moose-equation-cd38c247-688f-40bd-8f02-08bfdac1ce7a");katex.render("_2", element, {displayMode:false,throwOnError:false});$ as a function of temperature as represented by the MATPRO and Martin correlations. This is consistent with a similar comparison shown in Figure 6 of Martin (1988).

Example Input Syntax

[Materials<<<{"href": "../../../syntax/Materials/index.html"}>>>]
  [thermal_strain]
    type = UO2ThermalExpansionMATPROEigenstrain<<<{"description": "Computes eigenstrain due to thermal expansion in UO2 fuel using MATPRO correlations.", "href": "UO2ThermalExpansionMATPROEigenstrain.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 0
    temperature<<<{"description": "Coupled temperature"}>>> = temperature
    stress_free_temperature<<<{"description": "Reference temperature at which there is no thermal expansion for thermal eigenstrain calculation"}>>> = 300
    eigenstrain_name<<<{"description": "Material property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator."}>>> = thermal_eigenstrain
  []
[]

The eigenstrain name must also be passed to the strain calculator, and an example of this is shown below:

[Physics<<<{"href": "../../../syntax/Physics/index.html"}>>>]
  [SolidMechanics<<<{"href": "../../../syntax/Physics/SolidMechanics/index.html"}>>>]
    [QuasiStatic<<<{"href": "../../../syntax/Physics/SolidMechanics/QuasiStatic/index.html"}>>>]
      [all]
        strain<<<{"description": "Strain formulation"}>>> = FINITE
        eigenstrain_names<<<{"description": "List of eigenstrains to be applied in this strain calculation"}>>> = thermal_eigenstrain
      []
    []
  []
[]

Input Parameters

eigenstrain_nameMaterial property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator.
C++ Type:std::string
Controllable:No
Description:Material property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator.
stress_free_temperatureReference temperature at which there is no thermal expansion for thermal eigenstrain calculation
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Reference temperature at which there is no thermal expansion for thermal eigenstrain calculation

Required Parameters

base_nameOptional parameter that allows the user to define multiple mechanics material systems on the same block, i.e. for multiple phases
C++ Type:std::string
Controllable:No
Description:Optional parameter that allows the user to define multiple mechanics material systems on the same block, i.e. for multiple phases
blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
computeTrueWhen false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
Default:True
C++ Type:bool
Controllable:No
Description:When false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
constant_onNONEWhen ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
Default:NONE
C++ Type:MooseEnum
Options:NONE, ELEMENT, SUBDOMAIN
Controllable:No
Description:When ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
declare_suffixAn optional suffix parameter that can be appended to any declared 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 declared properties. The suffix will be prepended with a '_' character.
mean_thermal_expansion_coefficient_nameName of the mean coefficient of thermal expansion.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Name of the mean coefficient of thermal expansion.
temperatureCoupled temperature
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Coupled temperature
use_old_temperatureFalseFlag to optionally use the temperature value from the previous timestep.
Default:False
C++ Type:bool
Controllable:No
Description:Flag to optionally use the temperature value from the previous timestep.

Optional Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters

output_propertiesList of material properties, from this material, to output (outputs must also be defined to an output type)
C++ Type:std::vector<std::string>
Controllable:No
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
outputsnone Vector of output names where you would like to restrict the output of variables(s) associated with this object
Default:none
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

Outputs 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

thermal_expansion_scale_factor1The scaling factor on the thermal expansion strain.
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The scaling factor on the thermal expansion strain.

Advanced: Scaling Factor Parameters

Input Files

(test/tests/solid_mechanics/uo2_eigenstrains/uo2_matpro_cte.i)
(assessment/LWR/validation/RIA_CABRI_REP_Na4/analysis/REP_Na_4/REP_Na_4.i)
(assessment/LWR/validation/IFA_677/analysis/IFA_677_Base.i)
(assessment/LWR/validation/RIA_CABRI_REP_Na4/analysis/REP_Na_4/RIA/REP_Na_4_RIA.i)
(examples/accident_tolerant_fuel/uo2_fecral/uo2_fecral.i)
(test/tests/solid_mechanics/uo2_eigenstrains/ad_uo2_matpro_cte.i)
(examples/accident_tolerant_fuel/uo2_coated_zircaloy/uo2_coated_zircaloy.i)
(assessment/LWR/validation/Super_Ramp/analysis/Super_Ramp_Base.i)
(assessment/LWR/validation/Super_Ramp/analysis/PK62/PK62_weighted_gap_VCP.i)

References

C. M. Allison, G. A. Berna, R. Chambers, E. W. Coryell, K. L. Davis, D. L. Hagrman, D. T. Hagrman, N. L. Hampton, J. K. Hohorst, R. E. Mason, M. L. McComas, K. A. McNeil, R. L. Miller, C. S. Olsen, G. A. Reymann, and L. J. Siefken. SCDAP/RELAP5/MOD3.1 code manual, volume IV: MATPRO-A library of materials properties for light-water-reactor accident analysis. Technical Report NUREG/CR-6150, EGG-2720, Idaho National Engineering Laboratory, 1993.

@TECHREPORT{hagrman93,
    author = "Allison, C. M. and Berna, G. A. and Chambers, R. and Coryell, E. W. and Davis, K. L. and Hagrman, D. L. and Hagrman, D. T. and Hampton, N. L. and Hohorst, J. K. and Mason, R. E. and McComas, M. L. and McNeil, K. A. and Miller, R. L. and Olsen, C. S. and Reymann, G. A. and Siefken, L. J.",
    editor = "Hagrman, D. T.",
    title = "{SCDAP/RELAP5/MOD3.1} Code Manual, Volume {IV}: {MATPRO}-{A} Library of Materials Properties for Light-Water-Reactor Accident Analysis",
    institution = "Idaho National Engineering Laboratory",
    year = "1993",
    number = "NUREG/CR-6150, EGG-2720"
}

D. G. Martin. The thermal expansion of solid UO$_2$ and (U, Pu) mixed oxides - A review and recommendations. Journal of Nuclear Materials, 152(2):94–101, 1988.

@ARTICLE{martin88,
    author = "Martin, D. G.",
    title = "The Thermal Expansion of Solid {UO}$\_2$ and ({U}, {Pu}) Mixed Oxides - {A} Review and Recommendations",
    journal = "Journal of Nuclear Materials",
    year = "1988",
    volume = "152",
    pages = "94-101",
    number = "2"
}

(test/tests/solid_mechanics/uo2_eigenstrains/uo2_matpro_cte.i)

#--------------------------------------------------------------------------------
#
# This test case is designed to test thermal expansion model of UO2 fuel in a
# single element 1x1 square.
#
# - Temperature ramps from 300 K to 1500 K over 2 seconds with a nominal timestep of
#    0.1 seconds.  The equation for thermal expansion of uo2, with no Pu content and
#    no molten fraction is:
#     thermal_strain = 1.0e-5 * temperature - 3.0e-3 + 4.0e-2 * exp (-6.9e-20/(bk * temperature)),
#    where bk is the Boltzmann constant and is equal to 1.38e-23
#
#  The results for this nominal temperature change are given below.
#
#  Temperature (K)    Thermal Strain from MATPRO (m/m)       BISON-Computed Thermal Strain (m/m)
#       240              -6.000e-04                             -6.000023e-04
#       300               ~0.0 (reference temperature)           6.840728e-19
#       360              6.0004E-04                              6.000349e-04
#       420              1.2003E-03                              1.200268e-03
#       480              1.8012E-03                              1.801195e-03
#       540              2.4038E-03                              2.403807e-03
#       600              3.0096E-03                              3.009612e-03
#       660              3.6205E-03                              3.620507e-03
#       720              4.2386E-03                              4.238557e-03
#       780              4.8658E-03                              4.865782e-03
#       840              5.5040E-03                              5.503983e-03
#       900              6.1546E-03                              6.154634e-03
#       960              6.8188E-03                              6.818829e-03
#      1020              7.4973E-03                              7.497278e-03
#      1080              8.1903E-03                              8.190333e-03
#      1140              8.8980E-03                              8.898032e-03
#      1200              9.6202E-03                              9.620152e-03
#      1260              1.0356E-02                              1.035625e-02
#      1320              1.1106E-02                              1.110574e-02
#      1380              1.1868E-02                              1.186789e-02
#      1440              1.2642E-02                              1.264192e-02
#      1500              1.3427E-02                              1.342696e-02
#--------------------------------------------------------------------------------

[GlobalParams]
  displacements = 'disp_x disp_y' # disp_z'
  volumetric_locking_correction = true
[]

[Mesh]
  coord_type = RZ
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
  []
[]

[Variables]
  [disp_x]
    order = FIRST
    family = LAGRANGE
  []
  [disp_y]
    order = FIRST
    family = LAGRANGE
  []
[]

[AuxVariables]
  [thermal_eigenstrain_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [thermal_eigenstrain_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [temperature]
    order = FIRST
    family = LAGRANGE
    initial_condition = 180.0
  []
[]

[Functions]
  [temperaturefcn]
    type = ParsedFunction
    expression = 180+6.0e2*t
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = FINITE
    eigenstrain_names = thermal_eigenstrain
  []
[]

[AuxKernels]
  [fuel_thermal_strain_xx]
    type = RankTwoAux
    rank_two_tensor = thermal_eigenstrain
    variable = thermal_eigenstrain_xx
    index_j = 0
    index_i = 0
    execute_on = timestep_end
  []
  [fuel_thermal_strain_yy]
    type = RankTwoAux
    rank_two_tensor = thermal_eigenstrain
    variable = thermal_eigenstrain_yy
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  []
  [total_strain_xx]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_xx
    index_i = 0
    index_j = 0
    execute_on = timestep_end
  []
  [temperature]
    type = FunctionAux
    variable = temperature
    function = temperaturefcn
    execute_on = timestep_begin
  []
[]

[BCs]
  [u_bottom_fix]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  []
  [u_left_fix]
    type = DirichletBC
    variable = disp_x
    boundary = left
    value = 0.0
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = 0
    youngs_modulus = 2e11
    poissons_ratio = 0.3
  []

  [thermal_strain]
    type = UO2ThermalExpansionMATPROEigenstrain
    block = 0
    temperature = temperature
    stress_free_temperature = 300
    eigenstrain_name = thermal_eigenstrain
  []

  [elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 0
  []
[]

[Postprocessors]
  [temperature]
    type = AverageNodalVariableValue
    variable = temperature
  []
  [thermal_eigenstrain_xx]
    type = ElementAverageValue
    variable = thermal_eigenstrain_xx
  []
  [thermal_eigenstrain_yy]
    type = ElementAverageValue
    variable = thermal_eigenstrain_yy
  []

  [total_strain_xx]
    type = ElementAverageValue
    variable = total_strain_xx
  []
[]

[Executioner]
  type = Transient

  l_max_its = 100
  nl_max_its = 100
  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-10
  l_tol = 1e-6
  start_time = 0.0
  end_time = 2.2
  dt = 0.1
[]

[Outputs]
  csv = true
[]

(test/tests/solid_mechanics/uo2_eigenstrains/uo2_martin_cte.i)

#--------------------------------------------------------------------------------
#
# This test case is designed to test thermal expansion of UO2 fuel using the
# UO2ThermalExpansionMartinEigenstrain model, which is based of the model of
# Martin (1988).
#
# The Martin paper provides mean and instantaneous thermal expansion functions.
# The BISON implementation of this model uses the instantaneous version because
# there are errros in the mean equation for the higher temperature regime. This
# makes it difficult to compare the model against an analytic solution. The
# generated plot of the function qualitatively matches the function plotted in the
# original paper.

[GlobalParams]
  displacements = 'disp_x disp_y' # disp_z'
  volumetric_locking_correction = true
[]

[Mesh]
  coord_type = RZ
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
  []
[]

[Variables]
  [disp_x]
    order = FIRST
    family = LAGRANGE
  []
  [disp_y]
    order = FIRST
    family = LAGRANGE
  []
[]

[AuxVariables]
  [thermal_eigenstrain_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [thermal_eigenstrain_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [temperature]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300.0
  []
[]

[Functions]
  [temperaturefcn]
    type = ParsedFunction
    expression = 300+6.0e2*t
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = FINITE
    eigenstrain_names = thermal_eigenstrain
  []
[]

[AuxKernels]
  [fuel_thermal_strain_xx]
    type = RankTwoAux
    rank_two_tensor = thermal_eigenstrain
    variable = thermal_eigenstrain_xx
    index_j = 0
    index_i = 0
    execute_on = timestep_end
  []
  [fuel_thermal_strain_yy]
    type = RankTwoAux
    rank_two_tensor = thermal_eigenstrain
    variable = thermal_eigenstrain_yy
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  []
  [total_strain_xx]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_xx
    index_i = 0
    index_j = 0
    execute_on = timestep_end
  []
  [temperature]
    type = FunctionAux
    variable = temperature
    function = temperaturefcn
    execute_on = timestep_begin
  []
[]

[BCs]
  [u_bottom_fix]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  []
  [u_left_fix]
    type = DirichletBC
    variable = disp_x
    boundary = left
    value = 0.0
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = 0
    youngs_modulus = 2e11
    poissons_ratio = 0.3
  []

  [thermal_strain]
    type = UO2ThermalExpansionMartinEigenstrain
    block = 0
    temperature = temperature
    stress_free_temperature = 300
    eigenstrain_name = thermal_eigenstrain
  []

  [elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 0
  []
[]

[Postprocessors]
  [temperature]
    type = AverageNodalVariableValue
    variable = temperature
  []
  [total_strain_xx]
    type = ElementAverageValue
    variable = total_strain_xx
  []
[]

[Executioner]
  type = Transient

  l_max_its = 100
  nl_max_its = 100
  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-10
  l_tol = 1e-6
  start_time = 0.0
  end_time = 4.6
  dt = 0.1
[]

[Outputs]
  csv = true
  execute_on = timestep_end
[]

(test/tests/solid_mechanics/uo2_eigenstrains/uo2_matpro_cte.i)

#--------------------------------------------------------------------------------
#
# This test case is designed to test thermal expansion model of UO2 fuel in a
# single element 1x1 square.
#
# - Temperature ramps from 300 K to 1500 K over 2 seconds with a nominal timestep of
#    0.1 seconds.  The equation for thermal expansion of uo2, with no Pu content and
#    no molten fraction is:
#     thermal_strain = 1.0e-5 * temperature - 3.0e-3 + 4.0e-2 * exp (-6.9e-20/(bk * temperature)),
#    where bk is the Boltzmann constant and is equal to 1.38e-23
#
#  The results for this nominal temperature change are given below.
#
#  Temperature (K)    Thermal Strain from MATPRO (m/m)       BISON-Computed Thermal Strain (m/m)
#       240              -6.000e-04                             -6.000023e-04
#       300               ~0.0 (reference temperature)           6.840728e-19
#       360              6.0004E-04                              6.000349e-04
#       420              1.2003E-03                              1.200268e-03
#       480              1.8012E-03                              1.801195e-03
#       540              2.4038E-03                              2.403807e-03
#       600              3.0096E-03                              3.009612e-03
#       660              3.6205E-03                              3.620507e-03
#       720              4.2386E-03                              4.238557e-03
#       780              4.8658E-03                              4.865782e-03
#       840              5.5040E-03                              5.503983e-03
#       900              6.1546E-03                              6.154634e-03
#       960              6.8188E-03                              6.818829e-03
#      1020              7.4973E-03                              7.497278e-03
#      1080              8.1903E-03                              8.190333e-03
#      1140              8.8980E-03                              8.898032e-03
#      1200              9.6202E-03                              9.620152e-03
#      1260              1.0356E-02                              1.035625e-02
#      1320              1.1106E-02                              1.110574e-02
#      1380              1.1868E-02                              1.186789e-02
#      1440              1.2642E-02                              1.264192e-02
#      1500              1.3427E-02                              1.342696e-02
#--------------------------------------------------------------------------------

[GlobalParams]
  displacements = 'disp_x disp_y' # disp_z'
  volumetric_locking_correction = true
[]

[Mesh]
  coord_type = RZ
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
  []
[]

[Variables]
  [disp_x]
    order = FIRST
    family = LAGRANGE
  []
  [disp_y]
    order = FIRST
    family = LAGRANGE
  []
[]

[AuxVariables]
  [thermal_eigenstrain_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [thermal_eigenstrain_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [temperature]
    order = FIRST
    family = LAGRANGE
    initial_condition = 180.0
  []
[]

[Functions]
  [temperaturefcn]
    type = ParsedFunction
    expression = 180+6.0e2*t
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = FINITE
    eigenstrain_names = thermal_eigenstrain
  []
[]

[AuxKernels]
  [fuel_thermal_strain_xx]
    type = RankTwoAux
    rank_two_tensor = thermal_eigenstrain
    variable = thermal_eigenstrain_xx
    index_j = 0
    index_i = 0
    execute_on = timestep_end
  []
  [fuel_thermal_strain_yy]
    type = RankTwoAux
    rank_two_tensor = thermal_eigenstrain
    variable = thermal_eigenstrain_yy
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  []
  [total_strain_xx]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_xx
    index_i = 0
    index_j = 0
    execute_on = timestep_end
  []
  [temperature]
    type = FunctionAux
    variable = temperature
    function = temperaturefcn
    execute_on = timestep_begin
  []
[]

[BCs]
  [u_bottom_fix]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  []
  [u_left_fix]
    type = DirichletBC
    variable = disp_x
    boundary = left
    value = 0.0
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = 0
    youngs_modulus = 2e11
    poissons_ratio = 0.3
  []

  [thermal_strain]
    type = UO2ThermalExpansionMATPROEigenstrain
    block = 0
    temperature = temperature
    stress_free_temperature = 300
    eigenstrain_name = thermal_eigenstrain
  []

  [elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 0
  []
[]

[Postprocessors]
  [temperature]
    type = AverageNodalVariableValue
    variable = temperature
  []
  [thermal_eigenstrain_xx]
    type = ElementAverageValue
    variable = thermal_eigenstrain_xx
  []
  [thermal_eigenstrain_yy]
    type = ElementAverageValue
    variable = thermal_eigenstrain_yy
  []

  [total_strain_xx]
    type = ElementAverageValue
    variable = total_strain_xx
  []
[]

[Executioner]
  type = Transient

  l_max_its = 100
  nl_max_its = 100
  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-10
  l_tol = 1e-6
  start_time = 0.0
  end_time = 2.2
  dt = 0.1
[]

[Outputs]
  csv = true
[]

(assessment/LWR/validation/RIA_CABRI_REP_Na4/analysis/REP_Na_4/REP_Na_4.i)

# REP Na 4 Rodlet Base Irradiation


initial_fuel_density = 10476.35

[GlobalParams]
  density = ${initial_fuel_density} # assumed TD = 10970
  order = SECOND
  family = LAGRANGE
  energy_per_fission = 3.2e-11 # J/fission
  displacements = 'disp_x disp_y'
[]

[Problem]
  type = AugmentedLagrangianContactProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
  maximum_lagrangian_update_iterations = 200
  acceptable_iterations = 30
  acceptable_multiplier = 5
[]

[Mesh]
  coord_type = RZ
  patch_size = 40
  #patch_update_strategy = auto
  #partitioner = centroid
  #centroid_partitioner_direction = y
  [mesh]
    type = FileMeshGenerator
    file = REP_Na4.e
  []
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [temp]
    initial_condition = 293.15
  []
[]

[UserObjects]
  [pin_geometry]
    type = FuelPinGeometry
    clad_inner_wall = 5
    clad_outer_wall = 2
    clad_top = 3
    clad_bottom = 1
    pellet_exteriors = 8
  []
[]

[AuxVariables]
  [fast_neutron_flux]
    block = 'clad'
  []
  [fast_neutron_fluence]
    block = 'clad'
  []
  [grain_radius]
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    initial_condition = 5.0e-6 # Assume Grain size 10 microns
  []
  [gap_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [fuel_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [swelling_strain]
    order = CONSTANT
    family = MONOMIAL
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
  []
  [hoop_creep_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [hoop_plastic_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [axial_creep_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [hoop_elastic_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [axial_plastic_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_hoop_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_axial_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [axial_elastic_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [clad_coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_temp]
    order = CONSTANT
    family = MONOMIAL
  []
  [clad_coolant_flux]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_channel_hmode]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_channel_htype]
    order = CONSTANT
    family = MONOMIAL
  []
  [critical_heat_flux]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxide_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [plastic_strain_mag]
    order = CONSTANT
    family = MONOMIAL
    block = 'clad'
  []
  [SED]
    order = CONSTANT
    family = MONOMIAL
    block = 'clad'
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [power_profile]
    type = PiecewiseLinear # reads and interpolates an input file containing rod average linear power vs time
    data_file = REPNa4_power_history_Full.csv
    format = columns
    scale_factor = 1
  []
  [axial_peaking_factors] # reads and interpolates an input file containing the axial power profile vs time
    type = PiecewiseBilinear
    data_file = REPNa4_axial_peaking_Full.csv
    scale_factor = 1
    axis = 1
  []
  [pressure_ramp] # inlet coolant pressure evolution
    type = PiecewiseLinear
    scale_factor = 1.0
    xy_data = '0               101325
               8640            15499970
               124675200       15499970
               124718400       101325
               125193600       101325
               125193610       101325
               125193650       500008
               125193700       500008
               125193900       500008
               125194000       101325
               125194100       101325'
  []

  [temp_ramp] # inlet coolant temp evolution
    type = PiecewiseLinear
    scale_factor = 0.985
    xy_data = '0               293.15
               8640            591
               20476800        591
               21859200        600
               47692800        600
               51840000        593
               72144000        593
               73440000        586
               96940800        586
               99360000        583
               124675200       583
               124761600       293.15
               125193600       293.15
               125193650       553.15
               125193900.0     553.150
               125194000.0     293.150
               125194100.0     293.150'
  []

  [burnup_GWd]
    type = ParsedFunction
    expression = bu*950
    symbol_names = 'bu'
    symbol_values = 'average_burnup'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [pellets]
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    strain = FINITE
    eigenstrain_names = 'fuel_thermal_strain fuel_relocation_strain
      fuel_volumetric_strain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz axial_stress'
    extra_vector_tags = 'ref'
  []
  [clad]
    block = 'clad'
    strain = FINITE
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_strain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz plastic_strain_xx plastic_strain_yy plastic_strain_zz
      creep_strain_xx creep_strain_yy creep_strain_zz elastic_strain_xx
      elastic_strain_yy elastic_strain_zz hoop_stress axial_stress'
    extra_vector_tags = 'ref'
  []
[]

[Kernels]
  [heat] # gradient term in heat conduction equation
    type = HeatConduction
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_ie] # time term in heat conduction equation
    type = HeatConductionTimeDerivative
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_source_fuel] # source term in heat conduction equation
    type = NeutronHeatSource
    variable = temp
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    burnup_function = burnup
    extra_vector_tags = 'ref'
  []
[]

[Burnup]
  [burnup]
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    rod_ave_lin_pow = power_profile
    axial_power_profile = axial_peaking_factors
    num_radial = 80
    num_axial = 11
    a_upper = 0.5678974
    a_lower = 0.0045
    fuel_inner_radius = 0
    fuel_outer_radius = 0.0040959
    fuel_volume_ratio = 1
    isotopes = 'U235 U238 Pu239 Pu240 Pu241 Pu242'
    isotope_fractions = '0.0449 0.9551 0 0 0 0'
    RPF = RPF
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    axial_power_profile = axial_peaking_factors
    rod_ave_lin_pow = power_profile
    factor = 3e13 #n/m2-s
    block = 'clad'
    execute_on = timestep_begin
  []
  [fast_neutron_fluence]
    type = FastNeutronFluenceAux
    variable = fast_neutron_fluence
    fast_neutron_flux = fast_neutron_flux
    block = 'clad'
    execute_on = timestep_begin
  []
  [grain_radius]
    type = GrainRadiusAux
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    variable = grain_radius
    temperature = temp
    execute_on = linear
  []
  [gap_conductance]
    type = MaterialRealAux
    property = gap_conductance
    variable = gap_cond
    boundary = 10
  []
  [fuel_conductance]
    type = MaterialRealAux
    property = thermal_conductivity
    variable = fuel_cond
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
  []
  [swelling_strain]
    type = MaterialRealAux
    property = volumetric_swelling_strain
    variable = swelling_strain
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
  []
  [hoop_creep_strain]
    type = RankTwoAux
    rank_two_tensor = creep_strain
    variable = hoop_creep_strain
    index_i = 2
    index_j = 2
    execute_on = timestep_end
    block = 'clad'
  []
  [axial_creep_strain]
    type = RankTwoAux
    rank_two_tensor = creep_strain
    variable = axial_creep_strain
    index_i = 1
    index_j = 1
    execute_on = timestep_end
    block = 'clad'
  []
  [hoop_plastic_strain]
    type = RankTwoAux
    rank_two_tensor = plastic_strain
    variable = hoop_plastic_strain
    index_i = 2
    index_j = 2
    execute_on = timestep_end
    block = 'clad'
  []
  [axial_plastic_strain]
    type = RankTwoAux
    rank_two_tensor = plastic_strain
    variable = axial_plastic_strain
    index_i = 1
    index_j = 1
    execute_on = timestep_end
    block = 'clad'
  []
  [total_hoop_strain]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_hoop_strain
    index_i = 2
    index_j = 2
    execute_on = timestep_end
    block = 'clad'
  []
  [total_axial_strain]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_axial_strain
    index_i = 1
    index_j = 1
    execute_on = timestep_end
    block = 'clad'
  []
  [hoop_elastic_strain]
    type = RankTwoAux
    rank_two_tensor = elastic_strain
    variable = hoop_elastic_strain
    index_i = 2
    index_j = 2
    execute_on = timestep_end
    block = 'clad'
  []
  [axial_elastic_strain]
    type = RankTwoAux
    rank_two_tensor = elastic_strain
    variable = axial_elastic_strain
    index_i = 1
    index_j = 1
    execute_on = timestep_end
    block = 'clad'
  []
  [creep_strain_mag]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    block = clad
    execute_on = timestep_end
  []
  [plastic_strain_mag]
    type = MaterialRealAux
    property = effective_plastic_strain
    variable = plastic_strain_mag
    block = clad
    execute_on = timestep_end
  []
  [clad_coolant_htc]
    type = MaterialRealAux
    property = coolant_channel_htc
    variable = clad_coolant_htc
    boundary = 2
  []
  [coolant_temp]
    type = MaterialRealAux
    property = coolant_temperature
    variable = coolant_temp
    boundary = 2
  []
  [clad_coolant_flux]
    type = MaterialRealAux
    property = output_heat_flux
    variable = clad_coolant_flux
    boundary = 2
  []
  [coolant_channel_hmode]
    type = MaterialRealAux
    property = coolant_channel_hmode
    variable = coolant_channel_hmode
    boundary = 2
  []
  [coolant_channel_htype]
    type = MaterialRealAux
    property = coolant_channel_htype
    variable = coolant_channel_htype
    boundary = 2
  []
  [critical_heat_flux]
    type = MaterialRealAux
    property = critical_heat_flux
    variable = critical_heat_flux
    boundary = 2
  []
  [oxide]
    type = MaterialRealAux
    property = oxide_scale_thickness
    variable = oxide_thickness
    boundary = 2
  []
  [SED]
    type = MaterialRealAux
    variable = SED
    property = strain_energy_density
    execute_on = timestep_end
  []
  [ofract_total]
    type = MaterialRealAux
    property = current_oxygen_weight_frac_total
    variable = oxywtfract_total
    execute_on = timestep_end
    boundary = 2
  []
  [ofgain_total]
    type = MaterialRealAux
    property = oxygen_weight_frac_gained_total
    variable = oxywtfgain_total
    execute_on = timestep_end
    boundary = 2
  []
  [fract_bphase]
    type = MaterialRealAux
    property = fract_beta_phase
    variable = fract_beta_phase
    block = 'clad'
  []
[]

# Define mechanical contact between the fuel (sideset=10) and the clad (sideset=5)
[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 210 #10
    penalty = 1e9
    model = coulomb
    formulation = augmented_lagrange
    friction_coefficient = 10.0
    tangential_tolerance = 1e-3
    normal_smoothing_distance = 0.1
    al_penetration_tolerance = 1e-6
    al_incremental_slip_tolerance = 1e-6
    al_frictional_force_tolerance = 5e-2
  []

  [pellet_clad_mechanical_2]
    primary = 5
    secondary = 410
    penalty = 1e9
    model = coulomb
    formulation = augmented_lagrange
    friction_coefficient = 0.0
    tangential_tolerance = 1e-3
    normal_smoothing_distance = 0.1
    al_penetration_tolerance = 1e-6
    al_incremental_slip_tolerance = 1e-6
    al_frictional_force_tolerance = 5e-2
  []
[]

# Define thermal contact between the fuel (sideset=10) and the clad (sideset=5)
[ThermalContact]
  [thermal_contact]
    type = GasGapHeatTransfer
    variable = temp
    primary = 5
    secondary = 10
    initial_moles = initial_moles # coupling to a postprocessor which supplies the initial plenum/gap gas mass
    gas_released = fission_gas_released # coupling to a postprocessor which supplies the fission gas addition
    jump_distance_model = LANNING
    plenum_pressure = plenum_pressure
    contact_pressure = contact_pressure
    roughness_primary = 2.0e-6
    roughness_secondary = 0.5e-6
    roughness_coef = 3.2
    normal_smoothing_distance = 0.1
    quadrature = true
    emissivity_primary = 0.800 #Emissivity for fuel
    emissivity_secondary = 0.325 #Emissivity for clad
    refab_time = 125107200
    refab_gas_types = He
    refab_fractions = 1
    contact_coef = 20 #10 default
  []
[]

[BCs]
  # pin pellets and clad along axis of symmetry (y)
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 12
    value = 0.0
  []
  # pin clad bottom in the axial direction (y)
  [no_y_clad_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 1
    value = 0.0
  []
  # pin fuel bottom in the axial direction (y)
  [no_y_fuel_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 20
    value = 0.0
  []
  [Pressure]
    #  apply coolant pressure on clad outer walls
    [coolantPressure]
      boundary = '1 2 3'
      factor = 1
      function = pressure_ramp # use the pressure_ramp function defined above
    []
  []
  [PlenumPressure]
    #  apply plenum pressure on clad inner walls and pellet surfaces
    [plenumPressure]
      boundary = 9
      initial_temperature = 293.15
      initial_pressure = 2.60e6
      startup_time = 0
      R = 8.3143
      output_initial_moles = initial_moles # coupling to post processor to get initial fill gas mass
      temperature = ave_temp_interior #plenumTemp
      volume = gas_volume # coupling to post processor to get gas volume
      material_input = fission_gas_released # coupling to post processor to get fission gas added
      output = plenum_pressure # coupling to post processor to output plenum/gap pressure
      displacements = 'disp_x disp_y'
      #      extra_vector_tags = 'ref'
      refab_time = 125107200
      refab_pressure = 0.301e6
      refab_temperature = 293.15
      refab_volume = 2.0e-6
    []
  []
[]

[CoolantChannel]
  [convective_clad_surface] # apply convective boundary to clad outer surface
    boundary = '1 2 3'
    variable = temp
    inlet_temperature = temp_ramp # K
    inlet_pressure = pressure_ramp # Pa
    inlet_massflux = 3244.044104 # kg/m^2-sec
    rod_diameter = 0.00951 # m
    rod_pitch = 1.26e-2 # m
    coolant_material = 'water'
    compute_enthalpy = true
    oxide_thickness = oxide_thickness # coupled oxide_thickness
    number_axial_zone = 50
  []
  #
  #    [convective_clad_surface_sodium] # apply convective boundary to clad outer surface
  #      boundary = '1 2 3'
  #      variable = temp
  #      inlet_temperature = temp_ramp      # K
  #      inlet_pressure    = pressure_ramp   # Pa
  #      inlet_massflux    = 3533     # kg/m^2-sec Based on flow rate provided and flow area and estimated density of 885.1 kg/m^3
  #      flow_area = 8.74855e-5         #m^2
  #      heated_diameter = 1.172526e-2   #m
  #      hydraulic_diameter = 4.7e-3     #m
  #      heated_perimeter = 2.984513e-2  #m
  #      coolant_material = 'sodium'
  #      compute_enthalpy = true
  #      heat_transfer_mode = 0
  #      oxide_thickness  = oxide_thickness # coupled oxide_thickness
  #      number_axial_zone = 50
  #      rod_diameter = 0.0095   # m
  #      htc_correlation_type = 2
  #    []
[]

[Materials]
  [fuel_density]
    type = StrainAdjustedDensity
    strain_free_density = ${initial_fuel_density}
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
  []
  [fuel_thermal]
    type = UO2Thermal
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    temperature = temp
    burnup_function = burnup
    thermal_conductivity_model = NFIR
    initial_porosity = 0.045
  []
  [fuel_elasticity_tensor]
    type = UO2ElasticityTensor
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    matpro_poissons_ratio = 1
    matpro_youngs_modulus = 1
    temperature = temp
  []
  [fuel_elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
  []
  [fuel_thermal_expansion]
    type = UO2ThermalExpansionMATPROEigenstrain
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    temperature = temp
    stress_free_temperature = 293.15
    eigenstrain_name = fuel_thermal_strain
  []
  [fuel_relocation]
    type = UO2RelocationEigenstrain
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    burnup_function = burnup
    rod_ave_lin_pow = power_profile
    axial_power_profile = axial_peaking_factors
    relocation_activation1 = 5000
    burnup_relocation_stop = 0.0208
    relocation_model = ESCORE_modified
    eigenstrain_name = fuel_relocation_strain
    fuel_pin_geometry = 'pin_geometry'
  []
  [fuel_volumetric_swelling]
    type = UO2VolumetricSwellingEigenstrain
    gas_swelling_model_type = SIFGRS
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    temperature = temp
    burnup_function = burnup
    initial_fuel_density = 10476.35
    total_densification = 0.00675
    initial_porosity = 0.045
    eigenstrain_name = fuel_volumetric_strain
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    temperature = temp
    fission_rate = fission_rate
    initial_porosity = 0.045
    grain_radius = grain_radius
    gbs_model = true
    transient_option = NO_TRANSIENT
  []

  [clad_density]
    type = StrainAdjustedDensity
    block = 'clad'
    strain_free_density = 6550
  []
  [clad_thermal]
    type = ZryThermal
    block = 'clad'
    temperature = temp
  []
  [clad_elasticity_tensor]
    type = ZryElasticityTensor
    temperature = temp
    matpro_poissons_ratio = true
    matpro_youngs_modulus = true
    cold_work_factor = 0.5
    fast_neutron_fluence = fast_neutron_fluence
    block = 'clad'
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep clad_zryplasticity'
    block = 'clad'
  []
  [clad_zrycreep]
    type = ZryCreepLimbackHoppeUpdate
    block = 'clad'
    temperature = temp
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    creeprate_scale_factor = 1
    model_irradiation_creep = 1
    model_primary_creep = 1
    model_thermal_creep = 1
    max_inelastic_increment = 0.0002
  []
  [clad_zryplasticity]
    type = ZryPlasticityUpdate
    block = 'clad'
    temperature = temp
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    cold_work_factor = 0.5
    plasticity_model_type = MATPRO
    zircaloy_alloy_type = 4
    max_inelastic_increment = 0.0002
  []
  [clad_thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = 'clad'
    temperature = temp
    stress_free_temperature = 293.15
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [clad_irradiation_swelling]
    type = ZryIrradiationGrowthEigenstrain
    block = 'clad'
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = stress_relief_annealed
    eigenstrain_name = clad_irradiation_strain
  []
  [clad_oxidation]
    type = ZryOxidation
    boundary = 2
    clad_inner_radius = 0.00417789
    clad_outer_radius = 0.00475615
    use_coolant_channel = true
    temperature = temp
    fast_neutron_flux = fast_neutron_flux
    oxygen_weight_fraction_initial = 0.0012
  []
  [phase]
    type = ZrPhase
    block = 'clad'
    numerical_method = 2
    temperature = temp
  []
  [StrainEnergyDensity]
    type = StrainEnergyDensity
    block = 'clad'
    incremental = 1
  []
[]

[Dampers]
  [BoundingValueNodalDamper]
    type = BoundingValueNodalDamper
    max_value = 3200 # The maximum permissible iterative value for the variable.
    min_value = 200 # The minimum permissible iterative value for the variable.
    variable = temp # The name of the variable that this damper operates on
  []
  [contact_slip]
    type = ContactSlipDamper
    primary = 5
    secondary = 10
    min_damping_factor = 0.05
  []
[]

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

[Executioner]
  type = Transient

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

  line_search = 'contact'
  contact_line_search_allowed_lambda_cuts = 0
  contact_line_search_ltol = 0.5

  verbose = true
  l_max_its = 100
  l_tol = 1e-3

  nl_max_its = 40
  nl_rel_tol = 1e-4
  nl_abs_tol = 1e-10

  start_time = 0
  end_time = 125193600 #125194100 is the end time for the RIA
  dtmax = 1e6
  dtmin = 1e-7

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1.0e3
    optimal_iterations = 100
    iteration_window = 2
    linear_iteration_ratio = 100
    timestep_limiting_function = power_profile
    max_function_change = 1e6
    force_step_every_function_point = true
    timestep_limiting_postprocessor = material_timestep
    time_t = '125193610 125193620 125193630 125193640 125193650 125193660 125193670 125193680'
    time_dt = '10 10 10 10 10 10 10 10'
    growth_factor = 5
  []
  [Quadrature]
    order = FIFTH #SEVENTH
    side_order = SEVENTH #Comment out if order = SEVENTH
  []
[]

[Postprocessors]
  [ave_temp_interior] # average temperature of the cladding interior and all pellet exteriors
    type = SideAverageValue
    boundary = 9 #For RIA the node number is ##***8479***##
    variable = temp
    execute_on = 'initial linear'
  []
  [clad_inner_vol] # volume inside of cladding
    type = InternalVolume
    boundary = 7
    execute_on = 'initial linear'
  []
  [avg_clad_temp] # average temperature of cladding interior
    type = SideAverageValue
    boundary = 7
    variable = temp
    execute_on = 'initial timestep_end'
  []
  [max_fuel_temp]
    type = NodalExtremeValue
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    value_type = max
    variable = temp
    execute_on = 'initial timestep_end'
  []
  [min_fuel_temp]
    type = NodalExtremeValue
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    value_type = min
    variable = temp
    execute_on = 'initial timestep_end'
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = 'clad'
    value_type = max
    variable = temp
    execute_on = 'initial timestep_end'
  []
  [min_clad_temp]
    type = NodalExtremeValue
    block = 'clad'
    value_type = min
    variable = temp
    execute_on = 'initial timestep_end'
  []
  [fis_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    outputs = exodus
  []
  [fis_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    outputs = exodus
  []
  [gas_volume] # gas volume
    type = InternalVolume
    boundary = 9
    execute_on = 'initial linear'
  []
  [flux_from_clad] # area integrated heat flux from the cladding
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 5
    diffusivity = thermal_conductivity
  []
  [flux_from_fuel] # area integrated heat flux from the fuel
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 10
    diffusivity = thermal_conductivity
  []
  [1_rod_input_power]
    type = FunctionValuePostprocessor
    function = power_profile
  []
  [3_burnup_GWd]
    type = FunctionValuePostprocessor
    function = burnup_GWd
  []
  [vonmises_stress_fuel]
    type = ElementAverageValue
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    variable = vonmises_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = 'clad'
    variable = vonmises_stress
  []
  [z_average_RPF]
    type = ElementAverageValue
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    variable = RPF
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = 'clad'
  []

  ## Nodal values
  [FCT] #fuel centerline temperature
    type = NodalVariableValue
    variable = temp
    nodeid = 3866 #(0, 0.303375, 0)
  []
  [FST] #fuel surface temperature
    type = NodalVariableValue
    variable = temp
    nodeid = 3823 #(0.0040959, 0.303375, 0)
  []
  [CIST] #clad inner surface temperature
    type = NodalVariableValue
    variable = temp
    nodeid = 9557 #(0.0041779, 0.305106, 0)
  []
  [COST] #clad outer surface temperature
    type = NodalVariableValue
    variable = temp
    nodeid = 9547 #(0.00475615, 0.305106, 0)
  []
  [gap]
    type = NodalVariableValue
    variable = penetration
    nodeid = 3823 #(0.0040959, 0.303375, 0)
    use_displaced_mesh = true
  []

  #######################################
  [qpoint_penetration] #FOCE
    type = ElementalVariableValue
    variable = qpoint_penetration
    elementid = 1200
    use_displaced_mesh = 1
  []

  [penetration] #FOCN
    type = NodalVariableValue
    variable = penetration
    nodeid = 3823
    use_displaced_mesh = 1
  []

  [contact_pressure] #FOCN
    type = NodalVariableValue
    variable = contact_pressure
    nodeid = 3823
    use_displaced_mesh = 1
  []

  [gap_cond] #FOCE
    type = ElementalVariableValue
    variable = gap_cond
    elementid = 1200
    use_displaced_mesh = 1
  []

  [creep_hoop_strain] #COCE
    type = ElementalVariableValue
    variable = hoop_creep_strain
    elementid = 2981
    use_displaced_mesh = 1
  []

  [elastic_hoop_strain] #COCE
    type = ElementalVariableValue
    variable = hoop_elastic_strain
    elementid = 2981
    use_displaced_mesh = 1
  []

  [plastic_hoop_strain] #COCE
    type = ElementalVariableValue
    variable = hoop_plastic_strain
    elementid = 2981
    use_displaced_mesh = 1
  []

  [total_hoop_strain] #COCE
    type = ElementalVariableValue
    variable = total_hoop_strain
    elementid = 2981
    use_displaced_mesh = 1
  []

  [clad_hoop_stress] #COCE
    type = ElementalVariableValue
    variable = hoop_stress
    elementid = 2981
    use_displaced_mesh = 1
  []

  [clad_axial_elongation] #COTN
    type = NodalVariableValue
    variable = disp_y
    nodeid = 10755
    use_displaced_mesh = 1
  []

  [clad_oxide_thickness] #COCE
    type = ElementalVariableValue
    variable = oxide_thickness
    elementid = 2981
    use_displaced_mesh = 1
  []

  [clad_coolant_htc] #COCE
    type = ElementalVariableValue
    variable = clad_coolant_htc
    elementid = 2981
    use_displaced_mesh = 1
  []

  [coolant_temp] #COCE
    type = ElementalVariableValue
    variable = coolant_temp
    elementid = 2981
    use_displaced_mesh = 1
  []

  [clad_coolant_flux] #COCE
    type = ElementalVariableValue
    variable = clad_coolant_flux
    elementid = 2981
    use_displaced_mesh = 1
  []

  [coolant_channel_hmode] #COCE
    type = ElementalVariableValue
    variable = coolant_channel_hmode
    elementid = 2981
    use_displaced_mesh = 1
  []

  [coolant_channel_htype] #COCE
    type = ElementalVariableValue
    variable = coolant_channel_htype
    elementid = 2981
    use_displaced_mesh = 1
  []

  [critical_heat_flux] #COCE
    type = ElementalVariableValue
    variable = critical_heat_flux
    elementid = 2981
    use_displaced_mesh = 1
  []

  [fuel_centerline_temp] #FICN
    type = NodalVariableValue
    variable = temp
    nodeid = 3866
  []

  [fuel_surface_temp] #FOCN
    type = NodalVariableValue
    variable = temp
    nodeid = 3823
  []

  [clad_inner_surface_temp] #CICN
    type = NodalVariableValue
    variable = temp
    nodeid = 9557
  []

  [clad_outer_surface_temp] #COCN
    type = NodalVariableValue
    variable = temp
    nodeid = 9547
  []

  [fuel_axial_elongation] #FOTN
    type = NodalVariableValue
    variable = disp_y
    nodeid = 7739
  []

  [clad_radial_elongation] #COCN
    type = NodalVariableValue
    variable = disp_x
    nodeid = 9547
  []

  [fuel_radial_elongation] #FOCN
    type = NodalVariableValue
    variable = disp_x
    nodeid = 3823
  []

  [SED_PPN_O] #COCE
    type = ElementalVariableValue
    variable = SED
    elementid = 2981
    use_displaced_mesh = 1
  []

  [SED_PPN_I] #CICE
    type = ElementalVariableValue
    variable = SED
    elementid = 2984
    use_displaced_mesh = 1
  []

  [zz_OFract_PPN_O] #COCE
    type = ElementalVariableValue
    variable = oxywtfract_total
    elementid = 2981
    use_displaced_mesh = 1
  []

  [zz_OGain_PPN_O] #COCE
    type = ElementalVariableValue
    variable = oxywtfgain_total
    elementid = 2981
    use_displaced_mesh = 1
  []
  #######################################

  [max_clad_SED]
    type = ElementExtremeValue
    block = 'clad'
    variable = SED
    value_type = max
  []

  #Post processor to calculate radial average enthalpy. This postprocessor isnt available yet in BISON
  [z_RAE]
    type = RadialAverageEnthalpy
    vector_postprocessor = rad_temp
    radial_direction = x
    axial_direction = y
    axial_position = 0.3
    temperature_name = temp
  []
  [peak_RAE]
    type = TimeExtremeValue
    postprocessor = z_RAE
  []
[]

[VectorPostprocessors]
  [clad_dia]
    type = NodalValueSampler
    variable = disp_x
    boundary = 2
    sort_by = y
    outputs = 'outfile_clad_diameter'
  []
  [pellet_dia]
    type = NodalValueSampler
    variable = disp_x
    boundary = 10
    sort_by = y
    outputs = 'outfile_pellet_diameter'
  []

  #Location of peak power node at appoximately 0.3 m in mesh
  [rad_temp]
    type = NodalValueSampler
    block = 3
    sort_by = y
    variable = temp
    execute_on = timestep_end
    outputs = 'outfile_radial_temp'
  []
[]

[StandardLWRFuelRodOutputs]
  fuel_pellet_blocks = 'pellet_type_1 pellet_type_2 pellet_type_3'
  temperature = temp
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  time_step_interval =  1
  csv = true
  #exodus = true
  color = false
  [outfile_clad_diameter]
    type = CSV
    sync_times = '125194100'
    sync_only = true
  []
  [outfile_pellet_diameter]
    type = CSV
    sync_times = '125194100'
    sync_only = true
  []
  [outfile_radial_temp]
    type = CSV
    end_time = -100000
  []

  [console]
    type = Console
    output_linear = true
    max_rows = 10
  []
  [chkfile]
    type = CSV
    show = 'ave_temp_interior fission_gas_released FCT average_burnup peak_RAE'
    execute_on = 'FINAL'
  []
  [exodus]
    type = Exodus
    time_step_interval =  4
    end_time = 125193700
  []
  [exodus_RIA]
    type = Exodus
    time_step_interval =  3
    start_time = 125193695
  []
  [checkpoint_RIA]
    type = Checkpoint
    file_base = recover_files_RIA
    sync_times = '124761600 125107200 125193600 125193650 125193700 125193700.06 125193700.07 125193700.08 125193700.09 125193700.10 125193700.20 125193700.30 125193700.40 125193700.50 125193700.60 125193700.70 125193700.80 125193700.90 125193701.00 125193701.25 125193701.50 125193701.75 125193702.00 125193702.25 125193702.50 125193702.75 125193703.00 125193704.00 125193705.00'
  []
[]

[Debug]
  show_var_residual = 'disp_x disp_y temp'
  show_var_residual_norms = true
  show_material_props = true
[]

(assessment/LWR/validation/IFA_677/analysis/IFA_677_Base.i)

# This is a partial input file that contains characteristics common to the entire assessment case
# NOTE: This file is not meant to be run on its own as it requires information contained in the .params files

# Fuel material properties
total_densification = 0.09e-2 # (-)
oxygen_to_metal_ratio = 2.002 # (-)

# Cladding material properties
cladding_density = 6550.0 # kg/m^3

# Cladding geometry
clad_inner_radius = 4.65e-3 # m
clad_outer_radius = 5.375e-3 # m

# Rod geometry
a_lower = 0.001025 # m
fuel_inner_radius = 0.0 # m
fuel_outer_radius = 0.4565e-02 # m
fuel_volume_ratio = 0.9940 # (-)
fuel_diameter = 9.13e-3 # m
diametral_gap = 170.e-6 # m


# Neutronics, power, and isotope fractions
energy_per_fission = 3.28451e-11 # J/fission
isotope_fraction_Pu239 = 0.0
isotope_fraction_Pu240 = 0.0
isotope_fraction_Pu241 = 0.0
isotope_fraction_Pu242 = 0.0

# Temperature conditions
initial_temperature = 293.0 # K
stress_free_temperature = 293.0 # K

# Coolant pressure ramp parameters
pressure_ramp_x = '-200 0'
pressure_ramp_y = '0.0298 1'
pressure_ramp_factor = 3.4e6 # (-)

# Coolant Channel parameters
inlet_pressure = 3.4e+06 # Pa
inlet_massflux = 450. # kg/m^2-s
rod_diameter = 10.75e-03 # m
rod_pitch = 46.e-03 # m
htc_correlation_type = 2

# Contact
contact_penalty = 1e14 # (-)
roughness_primary = 5.0e-6
roughness_secondary = 0.32e-6
roughness_coef = 3.2
normal_smoothing_distance = 0.1 # m

# Relocation
relocation_activation1 = 5000 # W/m

# Plenum parameters
initial_plenum_pressure = 1.35e6 # Pa
startup_time = 0 # s

# Physical constants
gravitational_acceleration_constant = -9.81 # m/s^2
ideal_gas_constant = 8.3143 # J/mol-K

# Numerical options
l_max_its = 50
l_tol = 8e-3
nl_max_its = 20
nl_rel_tol = 1e-4
start_time = -200 # s
n_startup_steps = 1
dtmax = 5e5 # s
dtmin = 0.1 # s

TimeStepper_dt = 2e2 # s
TimeStepper_max_function_change = 1e20

# Postprocessor parameters
burnup_scaling_factor = 950.0 # (-)
time_days_scale_factor = 1.157407407e-5 # (-)
midplane_oxide_thickness_elementid = 629
outer_midplane_clad_temperature_nodeid = 763
fuel_central_temperature_nodeid = 337
clad_elongation_nodeid = 826
fuel_elongation_nodeid = 504
upper_TC_temperature_nodeid = 633
lower_TC_temperature_nodeid = 31

[GlobalParams]
  displacements = 'disp_x disp_y'
  temperature = temperature
  order = FIRST
  family = LAGRANGE
  density = ${initial_fuel_density}
  initial_porosity = ${initial_fuel_porosity}
  energy_per_fission = ${energy_per_fission}
  volumetric_locking_correction = false
[]

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

[Mesh]
  coord_type = RZ
  patch_update_strategy = auto
  patch_size = 10
  partitioner = centroid
  centroid_partitioner_direction = y
  [mesh]
    type = FileMeshGenerator
    file = ${rod_mesh_file}
  []
[]

[Variables]
  [temperature]
    initial_condition = ${initial_temperature}
  []
[]

[AuxVariables]
  [fast_neutron_flux]
  []
  [fast_neutron_fluence]
  []
  [grain_radius]
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    initial_condition = ${initial_grain_radius}
  []
  [gap_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [effective_creep_strain]
    block = clad
    order = CONSTANT
    family = MONOMIAL
  []
  [oxide_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    data_file = ${power_history_data_file}
    scale_factor = 1e3
    format = columns
  []
  [corrected_power_history]
    type = ParsedFunction
    symbol_names = 'rod_average_burnup power_history'
    symbol_values = 'burnup_MWdkgU power_history'
    expression = 'if(rod_average_burnup < 14.5, power_history, 0.98 * power_history)'
  []
  [axial_peaking_factors]
    type = PiecewiseBilinear
    data_file = ${axial_peaking_data_file}
    scale_factor = 1
    axis = 1
  []
  [power_correction]
    type = PiecewiseBilinear
    data_file = ${power_correction_data_file}
    scale_factor = 1
    axis = 1
  []
  [corrected_axial_peaking_factors]
    type = CompositeFunction
    functions = 'axial_peaking_factors power_correction'
  []
  [pressure_ramp]
    type = PiecewiseLinear
    x = ${pressure_ramp_x}
    y = ${pressure_ramp_y}
  []
  [coolant_inlet_temperature]
    type = PiecewiseLinear
    data_file = ${coolant_inlet_temperature_data_file}
    scale_factor = 1
    format = columns
  []
  [fast_flux]
    type = PiecewiseLinear
    data_file = ${fast_neutron_flux_data_file}
    scale_factor = 1e4
    format = columns
  []
[]

[Physics]

  [SolidMechanics]
    [QuasiStatic]
      [pellets]
        block = 'pellet_type_1 pellet_type_2 pellet_type_3'
        add_variables = true
        strain = FINITE
        eigenstrain_names = 'fuel_relocation_strain fuel_thermal_eigenstrain fuel_volumetric_strain'
        generate_output = 'vonmises_stress hydrostatic_stress stress_xx stress_yy stress_zz'
        extra_vector_tags = 'ref'
      []
      [clad]
        block = clad
        add_variables = true
        strain = FINITE
        eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_strain'
        generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx strain_yy strain_zz'
        extra_vector_tags = 'ref'
      []
    []
  []
[]

[Kernels]
  [gravity]
    type = Gravity
    variable = disp_y
    value = ${gravitational_acceleration_constant}
  []
  [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_type_1 pellet_type_2 pellet_type_3'
    burnup_function = burnup
  []
[]

[Burnup]
  [burnup]
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    rod_ave_lin_pow = corrected_power_history
    axial_power_profile = corrected_axial_peaking_factors
    num_radial = 80
    num_axial = 11
    a_lower = ${a_lower}
    a_upper = ${a_upper}
    fuel_inner_radius = ${fuel_inner_radius}
    fuel_outer_radius = ${fuel_outer_radius}
    fuel_volume_ratio = ${fuel_volume_ratio}
    order = CONSTANT
    family = MONOMIAL
    RPF = RPF
    isotopes = 'U235 U238 Pu239 Pu240 Pu241 Pu242'
    isotope_fractions = '${isotope_fraction_U235} ${isotope_fraction_U238} ${isotope_fraction_Pu239} ${isotope_fraction_Pu240} ${isotope_fraction_Pu241} ${isotope_fraction_Pu242}'
  []
[]

[AuxKernels]
  [fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    block = clad
    function = fast_flux
    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_type_1 pellet_type_2 pellet_type_3'
    variable = grain_radius
    temperature = temperature
    execute_on = linear
  []
  [conductance]
    type = MaterialRealAux
    property = gap_conductance
    variable = gap_cond
    boundary = 10
    execute_on = linear
  []
  [effective_creep_strain]
    type = MaterialRealAux
    block = 1
    property = effective_creep_strain
    variable = effective_creep_strain
    execute_on = 'timestep_end'
  []
  [oxide]
    type = MaterialRealAux
    property = oxide_scale_thickness
    variable = oxide_thickness
    boundary = 2
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    formulation = penalty
    model = frictionless
    normalize_penalty = true
    penalty = ${contact_penalty}
    normal_smoothing_distance = ${normal_smoothing_distance}
  []
[]

[ThermalContact]
  [thermal_contact]
    type = GasGapHeatTransfer
    variable = temperature
    primary = 5
    secondary = 10
    initial_moles = initial_moles
    gas_released = fission_gas_released
    contact_pressure = contact_pressure
    quadrature = true
    roughness_primary = ${roughness_primary}
    roughness_secondary = ${roughness_secondary}
    roughness_coef = ${roughness_coef}
  []
[]

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

[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 = '20'
    value = 0.0
  []
  [Pressure]
    [coolantPressure]
      boundary = '1 2 3'
      factor = ${pressure_ramp_factor}
      function = pressure_ramp
    []
  []
  [PlenumPressure]
    [plenumPressure]
      boundary = '9'
      initial_pressure = ${initial_plenum_pressure}
      startup_time = ${startup_time}
      R = ${ideal_gas_constant}
      output_initial_moles = initial_moles
      temperature = plenum_temp
      volume = gas_volume
      material_input = fission_gas_released
      output = plenum_pressure
    []
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = '1 2 3'
    variable = temperature
    inlet_temperature = coolant_inlet_temperature
    inlet_pressure = ${inlet_pressure}
    inlet_massflux = ${inlet_massflux}
    rod_diameter = ${rod_diameter}
    rod_pitch = ${rod_pitch}
    htc_correlation_type = ${htc_correlation_type}
    compute_enthalpy = true
    linear_heat_rate = corrected_power_history
    axial_power_profile = corrected_axial_peaking_factors
    oxide_thickness = oxide_thickness
  []
[]

[Materials]
  [fuel_thermal]
    type = UO2Thermal
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    thermal_conductivity_model = HALDEN #NFIR
    temperature = temperature
    burnup_function = burnup
  []
  [fuel_volumetric_swelling]
    type = UO2VolumetricSwellingEigenstrain
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    burnup_function = burnup
    total_densification = ${total_densification}
    initial_fuel_density = ${initial_fuel_density}
    gas_swelling_model_type = SIFGRS
    eigenstrain_name = fuel_volumetric_strain
  []
  [fuel_relocation]
    type = UO2RelocationEigenstrain
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    relocation_model = ESCORE_modified
    burnup_function = burnup
    diameter = ${fuel_diameter}
    diametral_gap = ${diametral_gap}
    rod_ave_lin_pow = corrected_power_history
    axial_power_profile = corrected_axial_peaking_factors
    relocation_activation1 = ${relocation_activation1}
    eigenstrain_name = fuel_relocation_strain
  []
  [fuel_thermal_expansion]
    type = UO2ThermalExpansionMATPROEigenstrain
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    stress_free_temperature = ${stress_free_temperature}
    eigenstrain_name = fuel_thermal_eigenstrain
  []
  [fuel_elasticity_tensor]
    type = UO2IsotropicDamageElasticityTensor
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    rod_ave_lin_pow = corrected_power_history
    fragmentation_model = BARANI
  []
  [fuel_stress]
    type = ComputeMultipleInelasticStress
    inelastic_models = 'fuel_creep'
    tangent_operator = elastic
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
  []
  [fuel_creep]
    type = UO2CreepUpdate
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    temperature = temperature
    fission_rate = fission_rate
    initial_grain_radius = ${initial_grain_radius}
    oxygen_to_metal_ratio = ${oxygen_to_metal_ratio}
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    strain_free_density = ${initial_fuel_density}
  []
  [clad_thermal]
    type = ZryThermal
    block = clad
  []
  [clad_elasticity_tensor]
    type = ZryElasticityTensor
    block = clad
    fast_neutron_fluence = fast_neutron_fluence
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep'
    block = clad
  []
  [clad_zrycreep]
    type = ZryCreepLimbackHoppeUpdate
    block = clad
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    model_primary_creep = true
    model_thermal_creep = true
    model_irradiation_creep = true
  []
  [clad_thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = clad
    stress_free_temperature = ${stress_free_temperature}
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [clad_irradiation_growth]
    type = ZryIrradiationGrowthEigenstrain
    block = clad
    fast_neutron_fluence = fast_neutron_fluence
    eigenstrain_name = clad_irradiation_strain
  []
  [clad_density]
    type = StrainAdjustedDensity
    block = clad
    strain_free_density = ${cladding_density}
  []
  [clad_oxidation]
    type = ZryOxidation
    boundary = 2
    clad_inner_radius = ${clad_inner_radius}
    clad_outer_radius = ${clad_outer_radius}
    use_coolant_channel = true
    fast_neutron_flux = fast_neutron_flux
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    temperature = temperature
    burnup_function = burnup
    grain_radius = grain_radius
    gbs_model = true
    transient_option = MICROCRACKING_BURNUP
    ig_bubble_model = NUCLEATION_RESOLUTION
    diff_coeff_option = TURNBULL_D1_4D2_4D3
    doping_type = CR2O3_DOPED
    cr_doped_option = BEST_ESTIMATE_1773
  []
[]

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

[Dampers]
  [BoundingValueNodalDamper]
    type = MaxIncrement
    max_increment = 20
    variable = temperature
  []
[]

[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 = ${l_max_its}
  l_tol = ${l_tol}
  nl_max_its = ${nl_max_its}
  nl_rel_tol = ${nl_rel_tol}
  nl_abs_tol = ${nl_abs_tol}
  start_time = ${start_time}
  n_startup_steps = ${n_startup_steps}
  end_time = ${end_time}
  dtmax = ${dtmax}
  dtmin = ${dtmin}

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = ${TimeStepper_dt}
    timestep_limiting_function = power_history
    max_function_change = ${TimeStepper_max_function_change}
    force_step_every_function_point = true
  []

  [Quadrature]
    order = FIFTH
    side_order = SEVENTH
  []
[]

[Postprocessors]
  [alhr_input]
    type = FunctionValuePostprocessor
    function = corrected_power_history
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_max]
    type = NodalExtremeValue
    variable = temperature
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    execute_on = 'initial timestep_end'
  []
  [burnup]
    type = RodAverageBurnup
    burnup_function = burnup
  []
  [burnup_MWdkgU]
    type = ScalePostprocessor
    value = burnup
    scaling_factor = ${burnup_scaling_factor}
  []
  [time_days]
    type = FunctionValuePostprocessor
    function = t
    scale_factor = ${time_days_scale_factor}
  []
  [temperature_clad_max]
    type = NodalExtremeValue
    block = clad
    value_type = max
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [oxide_thickness_midplane]
    type = ElementalVariableValue
    elementid = ${midplane_oxide_thickness_elementid}
    variable = oxide_thickness
    execute_on = 'initial timestep_end'
  []
  [strain_clad_hoop_max]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = clad
  []
  [temperature_clad_outer_midplane]
    type = NodalVariableValue
    nodeid = ${outer_midplane_clad_temperature_nodeid}
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_central]
    type = NodalVariableValue
    variable = temperature
    nodeid = ${fuel_central_temperature_nodeid}
    execute_on = 'initial timestep_end'
  []
  [gas_volume]
    type = InternalVolume
    boundary = '9'
    execute_on = 'initial linear'
  []
  [fission_gas_produced]
    type = ElementIntegralFisGasGeneratedSifgrs
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    execute_on = linear
  []
  [fission_gas_released]
    type = ElementIntegralFisGasReleasedSifgrs
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    execute_on = linear
  []
  [fgr_percent]
    type = FGRPercent
    fission_gas_released = fission_gas_released
    fission_gas_generated = fission_gas_produced
    execute_on = linear
  []
  [clad_elongation]
    type = NodalVariableValue
    variable = disp_y
    nodeid = ${clad_elongation_nodeid}
  []
  [fuel_elongation]
    type = NodalVariableValue
    variable = disp_y
    nodeid = ${fuel_elongation_nodeid}
  []
  [upper_TC_temperature]
    type = NodalVariableValue
    variable = temperature
    nodeid = ${upper_TC_temperature_nodeid}
  []
  [lower_TC_temperature]
    type = NodalVariableValue
    variable = temperature
    nodeid = ${lower_TC_temperature_nodeid}
  []
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  color = false
  [console]
    type = Console
    max_rows = 15
  []
  [chkfile]
    type = CSV
    file_base = '${id}_chkfile'
    show = 'lower_TC_temperature upper_TC_temperature fgr_percent plenum_pressure strain_clad_hoop_max'
    execute_on = 'FINAL'
  []
  [csv]
    type = CSV
    file_base = '${id}_csv'
  []
  [exodus]
    type = Exodus
    file_base = '${id}_exodus'
  []
[]

(assessment/LWR/validation/RIA_CABRI_REP_Na4/analysis/REP_Na_4/RIA/REP_Na_4_RIA.i)

# REP Na 4 RIA


initial_fuel_density = 10476.35

[GlobalParams]
  density = ${initial_fuel_density} # assumed TD = 10970
  order = SECOND
  family = LAGRANGE
  energy_per_fission = 3.2e-11 # J/fission
  displacements = 'disp_x disp_y'
[]

[Problem]
  coord_type = RZ
  type = AugmentedLagrangianContactProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
  maximum_lagrangian_update_iterations = 200
  acceptable_iterations = 30
  acceptable_multiplier = 10
[]

[Mesh]
  patch_size = 40
  #patch_update_strategy = auto
  #partitioner = centroid
  #centroid_partitioner_direction = y
  [mesh]
    type = FileMeshGenerator
    file = REP_Na4.e
  []
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [temp]
  []
[]

[UserObjects]
  [pin_geometry]
    type = FuelPinGeometry
    clad_inner_wall = 5
    clad_outer_wall = 2
    clad_top = 3
    clad_bottom = 1
    pellet_exteriors = 8
  []
[]

[AuxVariables]
  [fast_neutron_flux]
    block = 'clad'
  []
  [fast_neutron_fluence]
    block = 'clad'
  []
  [grain_radius]
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
  []
  [gap_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [fuel_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [swelling_strain]
    order = CONSTANT
    family = MONOMIAL
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
  []
  [hoop_creep_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [axial_creep_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [hoop_plastic_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [axial_plastic_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_hoop_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_axial_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [hoop_elastic_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [axial_elastic_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [clad_coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_temp]
    order = CONSTANT
    family = MONOMIAL
  []
  [clad_coolant_flux]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_channel_hmode]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_channel_htype]
    order = CONSTANT
    family = MONOMIAL
  []
  [critical_heat_flux]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxide_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [plastic_strain_mag]
    order = CONSTANT
    family = MONOMIAL
    block = 'clad'
  []
  [SED]
    order = CONSTANT
    family = MONOMIAL
    block = 'clad'
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [power_profile]
    type = PiecewiseLinear # reads and interpolates an input file containing rod average linear power vs time
    data_file = REPNa4_power_history_Full.csv
    format = columns
    scale_factor = 1
  []
  [axial_peaking_factors] # reads and interpolates an input file containing the axial power profile vs time
    type = PiecewiseBilinear
    data_file = REPNa4_axial_peaking_Full.csv
    scale_factor = 1
    axis = 1
  []
  [pressure_ramp] # inlet coolant pressure evolution
    type = PiecewiseLinear
    format = columns
    scale_factor = 1.0
    xy_data = '0               101325
               8640            15499970
               124675200       15499970
               124718400       101325
               125193600       101325
               125193610       101325
               125193650       500008
               125193700       500008
               125193900       500008
               125194000       101325
               125194100       101325'
  []

  [temp_ramp] # inlet coolant temp evolution
    type = PiecewiseLinear
    format = columns
    scale_factor = 1.0
    xy_data = '0               293.15
               8640            591
               20476800        591
               21859200        600
               47692800        600
               51840000        593
               72144000        593
               73440000        586
               96940800        586
               99360000        583
               124675200.0     583
               124761600.0     293.150
               125193600.0     293.150
               125193650.0     553.150
               125193900.0     553.150
               125194000.0     293.150
               125194100.0     293.150'
  []

  [burnup_GWd]
    type = ParsedFunction
    expression = bu*950
    symbol_names = 'bu'
    symbol_values = 'average_burnup'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [pellets]
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    strain = FINITE
    eigenstrain_names = 'fuel_thermal_strain fuel_relocation_strain
      fuel_volumetric_strain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz axial_stress'
    extra_vector_tags = 'ref'
  []
  [clad]
    block = 'clad'
    strain = FINITE
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_strain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz plastic_strain_xx plastic_strain_yy plastic_strain_zz
      creep_strain_xx creep_strain_yy creep_strain_zz elastic_strain_xx
      elastic_strain_yy elastic_strain_zz hoop_stress axial_stress'
    extra_vector_tags = 'ref'
  []
[]

[Kernels]
  [heat] # gradient term in heat conduction equation
    type = HeatConduction
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_ie] # time term in heat conduction equation
    type = HeatConductionTimeDerivative
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_source_fuel] # source term in heat conduction equation
    type = NeutronHeatSource
    variable = temp
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    burnup_function = burnup
    extra_vector_tags = 'ref'
  []
[]

[Burnup]
  [burnup]
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    rod_ave_lin_pow = power_profile
    axial_power_profile = axial_peaking_factors
    num_radial = 80
    num_axial = 11
    a_upper = 0.5678974
    a_lower = 0.0045
    fuel_inner_radius = 0
    fuel_outer_radius = 0.0040959
    fuel_volume_ratio = 1
    isotopes = 'U235 U238 Pu239 Pu240 Pu241 Pu242'
    isotope_fractions = '0.0449 0.9551 0 0 0 0'
    RPF = RPF
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    axial_power_profile = axial_peaking_factors
    rod_ave_lin_pow = power_profile
    factor = 3e13 #n/m2-s
    block = 'clad'
    execute_on = timestep_begin
  []
  [fast_neutron_fluence]
    type = FastNeutronFluenceAux
    variable = fast_neutron_fluence
    fast_neutron_flux = fast_neutron_flux
    block = 'clad'
    execute_on = timestep_begin
  []
  [grain_radius]
    type = GrainRadiusAux
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    variable = grain_radius
    temperature = temp
    execute_on = linear
  []
  [gap_conductance]
    type = MaterialRealAux
    property = gap_conductance
    variable = gap_cond
    boundary = 10
  []
  [fuel_conductance]
    type = MaterialRealAux
    property = thermal_conductivity
    variable = fuel_cond
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
  []
  [swelling_strain]
    type = MaterialRealAux
    property = volumetric_swelling_strain
    variable = swelling_strain
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
  []
  [vonmises_stress]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = vonmises_stress
    scalar_type = VonMisesStress
    execute_on = timestep_end
  []
  [hoop_creep_strain]
    type = RankTwoAux
    rank_two_tensor = creep_strain
    variable = hoop_creep_strain
    index_i = 2
    index_j = 2
    execute_on = timestep_end
    block = 'clad'
  []
  [axial_creep_strain]
    type = RankTwoAux
    rank_two_tensor = creep_strain
    variable = axial_creep_strain
    index_i = 1
    index_j = 1
    execute_on = timestep_end
    block = 'clad'
  []
  [hoop_plastic_strain]
    type = RankTwoAux
    rank_two_tensor = plastic_strain
    variable = hoop_plastic_strain
    index_i = 2
    index_j = 2
    execute_on = timestep_end
    block = 'clad'
  []
  [axial_plastic_strain]
    type = RankTwoAux
    rank_two_tensor = plastic_strain
    variable = axial_plastic_strain
    index_i = 1
    index_j = 1
    execute_on = timestep_end
    block = 'clad'
  []
  [total_hoop_strain]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_hoop_strain
    index_i = 2
    index_j = 2
    execute_on = timestep_end
    block = 'clad'
  []
  [total_axial_strain]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_axial_strain
    index_i = 1
    index_j = 1
    execute_on = timestep_end
    block = 'clad'
  []
  [hoop_elastic_strain]
    type = RankTwoAux
    rank_two_tensor = elastic_strain
    variable = hoop_elastic_strain
    index_i = 2
    index_j = 2
    execute_on = timestep_end
    block = 'clad'
  []
  [axial_elastic_strain]
    type = RankTwoAux
    rank_two_tensor = elastic_strain
    variable = axial_elastic_strain
    index_i = 1
    index_j = 1
    execute_on = timestep_end
    block = 'clad'
  []
  [creep_strain_mag]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    block = clad
    execute_on = timestep_end
  []
  [plastic_strain_mag]
    type = MaterialRealAux
    property = effective_plastic_strain
    variable = plastic_strain_mag
    block = clad
    execute_on = timestep_end
  []
  [clad_coolant_htc]
    type = MaterialRealAux
    property = coolant_channel_htc
    variable = clad_coolant_htc
    boundary = 2
  []
  [coolant_temp]
    type = MaterialRealAux
    property = coolant_temperature
    variable = coolant_temp
    boundary = 2
  []
  [clad_coolant_flux]
    type = MaterialRealAux
    property = output_heat_flux
    variable = clad_coolant_flux
    boundary = 2
  []
  [coolant_channel_hmode]
    type = MaterialRealAux
    property = coolant_channel_hmode
    variable = coolant_channel_hmode
    boundary = 2
  []
  [coolant_channel_htype]
    type = MaterialRealAux
    property = coolant_channel_htype
    variable = coolant_channel_htype
    boundary = 2
  []
  [critical_heat_flux]
    type = MaterialRealAux
    property = critical_heat_flux
    variable = critical_heat_flux
    boundary = 2
  []
  [oxide]
    type = MaterialRealAux
    property = oxide_scale_thickness
    variable = oxide_thickness
    boundary = 2
  []
  [SED]
    type = MaterialRealAux
    variable = SED
    property = strain_energy_density
    execute_on = timestep_end
  []
  [ofract_total]
    type = MaterialRealAux
    property = current_oxygen_weight_frac_total
    variable = oxywtfract_total
    execute_on = timestep_end
    boundary = 2
  []
  [ofgain_total]
    type = MaterialRealAux
    property = oxygen_weight_frac_gained_total
    variable = oxywtfgain_total
    execute_on = timestep_end
    boundary = 2
  []
  [fract_bphase]
    type = MaterialRealAux
    property = fract_beta_phase
    variable = fract_beta_phase
    block = 'clad'
  []
[]

# Define mechanical contact between the fuel (sideset=10) and the clad (sideset=5)
[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 210 #10
    penalty = 1e9
    model = coulomb
    formulation = augmented_lagrange
    friction_coefficient = 0.3
    tangential_tolerance = 1e-3
    normal_smoothing_distance = 0.1
    al_penetration_tolerance = 1e-6
    al_incremental_slip_tolerance = 1e-6
    al_frictional_force_tolerance = 5e-2
  []

  [pellet_clad_mechanical_2]
    primary = 5
    secondary = 410
    penalty = 1e9
    model = coulomb
    formulation = augmented_lagrange
    friction_coefficient = 0.0
    tangential_tolerance = 1e-3
    normal_smoothing_distance = 0.1
    al_penetration_tolerance = 1e-6
    al_incremental_slip_tolerance = 1e-6
    al_frictional_force_tolerance = 5e-2
  []
[]

# Define thermal contact between the fuel (sideset=10) and the clad (sideset=5)
[ThermalContact]
  [thermal_contact]
    type = GasGapHeatTransfer
    variable = temp
    primary = 5
    secondary = 10
    initial_moles = initial_moles # coupling to a postprocessor which supplies the initial plenum/gap gas mass
    gas_released = fis_gas_released # coupling to a postprocessor which supplies the fission gas addition
    jump_distance_model = LANNING
    plenum_pressure = plenum_pressure
    contact_pressure = contact_pressure
    roughness_primary = 0.1e-6 #2.0e-6
    roughness_secondary = 0.1e-6 #0.5e-6
    roughness_coef = 3.2
    normal_smoothing_distance = 0.1
    quadrature = true
    emissivity_primary = 0.800 #Emissivity for fuel
    emissivity_secondary = 0.325 #Emissivity for clad
    refab_time = 125107200
    refab_gas_types = He
    refab_fractions = 1
    contact_coef = 20 #10 default
  []
[]

[BCs]
  # pin pellets and clad along axis of symmetry (y)
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 12
    value = 0.0
  []
  # pin clad bottom in the axial direction (y)
  [no_y_clad_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 1
    value = 0.0
  []
  # pin fuel bottom in the axial direction (y)
  [no_y_fuel_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 20
    value = 0.0
  []
  [Pressure]
    #  apply coolant pressure on clad outer walls
    [coolantPressure]
      boundary = '1 2 3'
      factor = 1
      function = pressure_ramp # use the pressure_ramp function defined above
    []
  []
  [PlenumPressure]
    #  apply plenum pressure on clad inner walls and pellet surfaces
    [plenumPressure]
      boundary = 9
      initial_temperature = 293.15
      initial_pressure = 2.60e6
      startup_time = 0
      R = 8.3143
      output_initial_moles = initial_moles # coupling to post processor to get initial fill gas mass
      temperature = ave_temp_interior #plenumTemp
      volume = gas_volume # coupling to post processor to get gas volume
      material_input = fis_gas_released # coupling to post processor to get fission gas added
      output = plenum_pressure # coupling to post processor to output plenum/gap pressure
      displacements = 'disp_x disp_y'
      #      extra_vector_tags = 'ref'
      refab_time = 125107200
      refab_pressure = 0.301e6
      refab_temperature = 293.15
      refab_volume = 2.0e-6
    []
  []
[]

[CoolantChannel]
  #    [convective_clad_surface_water] # apply convective boundary to clad outer surface
  #      boundary = '1 2 3'
  #      variable = temp
  #      inlet_temperature = temp_ramp      # K
  #      inlet_pressure    = pressure_ramp   # Pa
  #      inlet_massflux    = 3244.044104     # kg/m^2-sec
  #      rod_diameter      = 0.00951   # m
  #      rod_pitch         = 1.26e-2  # m
  #      coolant_material = 'water'
  #      compute_enthalpy = true
  #      oxide_thickness  = oxide_thickness # coupled oxide_thickness
  #      number_axial_zone = 50
  #    []
  #
  [convective_clad_surface] # apply convective boundary to clad outer surface
    boundary = '1 2 3'
    variable = temp
    inlet_temperature = temp_ramp # K
    inlet_pressure = pressure_ramp # Pa
    inlet_massflux = 3533 # kg/m^2-sec Based on flow rate provided and flow area and estimated density of 885.1 kg/m^3
    flow_area = 8.74855e-5 #m^2
    heated_diameter = 1.172526e-2 #m
    hydraulic_diameter = 4.7e-3 #m
    heated_perimeter = 2.984513e-2 #m
    coolant_material = 'sodium'
    compute_enthalpy = true
    heat_transfer_mode = 0
    oxide_thickness = oxide_thickness # coupled oxide_thickness
    number_axial_zone = 50
    rod_diameter = 0.0095 # m
    htc_correlation_type = 2
  []
[]

[Materials]
  [fuel_density]
    type = StrainAdjustedDensity
    #density = 10476.35
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
  []
  [fuel_thermal]
    type = UO2Thermal
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    temperature = temp
    burnup_function = burnup
    thermal_conductivity_model = NFIR
    initial_porosity = 0.045
  []
  [fuel_elasticity_tensor]
    type = UO2ElasticityTensor
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    matpro_poissons_ratio = 1
    matpro_youngs_modulus = 1
    temperature = temp
  []
  [fuel_elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
  []
  [fuel_thermal_expansion]
    type = UO2ThermalExpansionMATPROEigenstrain
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    temperature = temp
    stress_free_temperature = 293.15
    eigenstrain_name = fuel_thermal_strain
  []
  [fuel_relocation]
    type = UO2RelocationEigenstrain
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    burnup_function = burnup
    rod_ave_lin_pow = power_profile
    axial_power_profile = axial_peaking_factors
    relocation_activation1 = 5000
    burnup_relocation_stop = 0.0 #0.0208
    relocation_model = ESCORE_modified
    eigenstrain_name = fuel_relocation_strain
    fuel_pin_geometry = 'pin_geometry'
  []
  [fuel_volumetric_swelling]
    type = UO2VolumetricSwellingEigenstrain
    gas_swelling_model_type = SIFGRS
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    temperature = temp
    burnup_function = burnup
    initial_fuel_density = 10476.35
    total_densification = 0.00675
    initial_porosity = 0.045
    eigenstrain_name = fuel_volumetric_strain
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    temperature = temp
    fission_rate = fission_rate
    burnup_function = burnup
    initial_porosity = 0.045
    grain_radius_const = 5.0e-6
    gbs_model = false
    transient_option = MICROCRACKING_BURNUP
  []

  [clad_density]
    type = StrainAdjustedDensity
    block = 'clad'
    strain_free_density = 6550
  []
  [clad_thermal]
    type = ZryThermal
    block = 'clad'
    temperature = temp
  []
  [clad_elasticity_tensor]
    type = ZryElasticityTensor
    temperature = temp
    matpro_poissons_ratio = true
    matpro_youngs_modulus = true
    cold_work_factor = 0.5
    fast_neutron_fluence = fast_neutron_fluence
    block = 'clad'
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    #inelastic_models = 'clad_zrycreep clad_zryplasticity'
    inelastic_models = 'clad_zryplasticity'
    block = 'clad'
  []
  [clad_zrycreep]
    type = ZryCreepLimbackHoppeUpdate
    block = 'clad'
    temperature = temp
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    #creeprate_scale_factor = 1
    model_irradiation_creep = 1
    model_primary_creep = 1
    model_thermal_creep = 1
    max_inelastic_increment = 0.0001

    creeprate_scale_factor = 0
    enable = 0
  []
  [clad_zryplasticity]
    type = ZryPlasticityUpdate
    block = 'clad'
    temperature = temp
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    cold_work_factor = 0.5
    plasticity_model_type = MATPRO
    zircaloy_alloy_type = 4
    max_inelastic_increment = 0.0001
  []
  [clad_thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = 'clad'
    temperature = temp
    stress_free_temperature = 293.15
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [clad_irradiation_swelling]
    type = ZryIrradiationGrowthEigenstrain
    block = 'clad'
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = stress_relief_annealed
    eigenstrain_name = clad_irradiation_strain
    compute = 0
  []
  [clad_oxidation]
    type = ZryOxidation
    boundary = 2
    clad_inner_radius = 0.00417789
    clad_outer_radius = 0.00475615
    use_coolant_channel = true
    temperature = temp
    fast_neutron_flux = fast_neutron_flux
    oxygen_weight_fraction_initial = 0.0012
  []
  [phase]
    type = ZrPhase
    block = 'clad'
    numerical_method = 2
    temperature = temp
  []
  [StrainEnergyDensity]
    type = StrainEnergyDensity
    block = 'clad'
    incremental = 1
  []
[]

[Dampers]
  [BoundingValueNodalDamper]
    type = BoundingValueNodalDamper
    max_value = 3200 # The maximum permissible iterative value for the variable.
    min_value = 200 # The minimum permissible iterative value for the variable.
    variable = temp # The name of the variable that this damper operates on
  []
  [contact_slip]
    type = ContactSlipDamper
    primary = 5
    secondary = 10
    min_damping_factor = 0.05
  []
[]

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

[Executioner]
  type = Transient

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

  line_search = 'contact'
  contact_line_search_allowed_lambda_cuts = 0
  contact_line_search_ltol = 0.5

  verbose = true
  l_max_its = 100
  l_tol = 1e-3

  nl_max_its = 40
  nl_rel_tol = 1e-4
  nl_abs_tol = 1e-10

  start_time = 0
  end_time = 125194100 #125193600    #125194100 is the end time for the RIA
  dtmax = 10
  dtmin = 1e-7

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 10
    optimal_iterations = 20
    iteration_window = 4
    linear_iteration_ratio = 100
    timestep_limiting_function = power_profile
    max_function_change = 5e5
    force_step_every_function_point = true
    timestep_limiting_postprocessor = material_timestep
    time_t = '125193610 125193620 125193630 125193640 125193650 125193660 125193670 125193680'
    time_dt = '10 10 10 10 10 10 10 10'
  []
  [Quadrature]
    order = FIFTH #SEVENTH
    side_order = SEVENTH #Comment out if order = SEVENTH
  []
[]

[Postprocessors]
  #  [ave_temp_interior]            # average temperature of the cladding interior and all pellet exteriors
  #    type = SideAverageValue
  #    boundary = 9                   #For RIA the node number is ##***8479***##
  #    variable = temp
  #    execute_on = 'initial linear'
  #  []
  [ave_temp_interior]
    type = NodalVariableValue
    variable = temp
    nodeid = 8479
    execute_on = 'initial linear'
  []
  [clad_inner_vol] # volume inside of cladding
    type = InternalVolume
    boundary = 7
    execute_on = 'initial linear'
  []
  [avg_clad_temp] # average temperature of cladding interior
    type = SideAverageValue
    boundary = 7
    variable = temp
    execute_on = 'initial timestep_end'
  []
  [max_fuel_temp]
    type = NodalExtremeValue
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    value_type = max
    variable = temp
    execute_on = 'initial timestep_end'
  []
  [min_fuel_temp]
    type = NodalExtremeValue
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    value_type = min
    variable = temp
    execute_on = 'initial timestep_end'
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = 'clad'
    value_type = max
    variable = temp
    execute_on = 'initial timestep_end'
  []
  [min_clad_temp]
    type = NodalExtremeValue
    block = 'clad'
    value_type = min
    variable = temp
    execute_on = 'initial timestep_end'
  []
  [fis_gas_generated]
    type = ElementIntegralFisGasGeneratedSifgrs
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
  []
  [fis_gas_released]
    type = ElementIntegralFisGasReleasedSifgrs
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
  []
  [fis_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    outputs = exodus
  []
  [fis_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    outputs = exodus
  []
  [gas_volume] # gas volume
    type = InternalVolume
    boundary = 9
    execute_on = 'initial linear'
  []
  [flux_from_clad] # area integrated heat flux from the cladding
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 5
    diffusivity = thermal_conductivity
  []
  [flux_from_fuel] # area integrated heat flux from the fuel
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 10
    diffusivity = thermal_conductivity
  []
  [1_rod_input_power]
    type = FunctionValuePostprocessor
    function = power_profile
  []
  [3_burnup_GWd]
    type = FunctionValuePostprocessor
    function = burnup_GWd
  []
  [fis_gas_percent]
    type = FGRPercent
    fission_gas_released = fis_gas_released
    fission_gas_generated = fis_gas_generated
  []
  [vonmises_stress_fuel]
    type = ElementAverageValue
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    variable = vonmises_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = 'clad'
    variable = vonmises_stress
  []
  [z_average_RPF]
    type = ElementAverageValue
    block = 'pellet_type_1 pellet_type_2 pellet_type_3'
    variable = RPF
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = 'clad'
  []
  ## Nodal values
  [FCT] #fuel centerline temperature
    type = NodalVariableValue
    variable = temp
    nodeid = 3866 #(0, 0.303375, 0)
  []
  [FST] #fuel surface temperature
    type = NodalVariableValue
    variable = temp
    nodeid = 3823 #(0.0040959, 0.303375, 0)
  []
  [CIST] #clad inner surface temperature
    type = NodalVariableValue
    variable = temp
    nodeid = 9557 #(0.0041779, 0.305106, 0)
  []
  [COST] #clad outer surface temperature
    type = NodalVariableValue
    variable = temp
    nodeid = 9547 #(0.00475615, 0.305106, 0)
  []
  [gap]
    type = NodalVariableValue
    variable = penetration
    nodeid = 3823 #(0.0040959, 0.303375, 0)
    use_displaced_mesh = true
  []

  #######################################
  [qpoint_penetration] #FOCE
    type = ElementalVariableValue
    variable = qpoint_penetration
    elementid = 1200
    use_displaced_mesh = 1
  []

  [penetration] #FOCN
    type = NodalVariableValue
    variable = penetration
    nodeid = 3823
    use_displaced_mesh = 1
  []

  [contact_pressure] #FOCN
    type = NodalVariableValue
    variable = contact_pressure
    nodeid = 3823
    use_displaced_mesh = 1
  []

  [gap_cond] #FOCE
    type = ElementalVariableValue
    variable = gap_cond
    elementid = 1200
    use_displaced_mesh = 1
  []

  [creep_hoop_strain] #COCE
    type = ElementalVariableValue
    variable = hoop_creep_strain
    elementid = 2981
    use_displaced_mesh = 1
  []

  [elastic_hoop_strain] #COCE
    type = ElementalVariableValue
    variable = hoop_elastic_strain
    elementid = 2981
    use_displaced_mesh = 1
  []

  [plastic_hoop_strain] #COCE
    type = ElementalVariableValue
    variable = hoop_plastic_strain
    elementid = 2981
    use_displaced_mesh = 1
  []

  [total_hoop_strain] #COCE
    type = ElementalVariableValue
    variable = total_hoop_strain
    elementid = 2981
    use_displaced_mesh = 1
  []

  [clad_hoop_stress] #COCE
    type = ElementalVariableValue
    variable = hoop_stress
    elementid = 2981
    use_displaced_mesh = 1
  []

  [clad_axial_elongation] #COTN
    type = NodalVariableValue
    variable = disp_y
    nodeid = 10755
    use_displaced_mesh = 1
  []

  [clad_oxide_thickness] #COCE
    type = ElementalVariableValue
    variable = oxide_thickness
    elementid = 2981
    use_displaced_mesh = 1
  []

  [clad_coolant_htc] #COCE
    type = ElementalVariableValue
    variable = clad_coolant_htc
    elementid = 2981
    use_displaced_mesh = 1
  []

  [coolant_temp] #COCE
    type = ElementalVariableValue
    variable = coolant_temp
    elementid = 2981
    use_displaced_mesh = 1
  []

  [clad_coolant_flux] #COCE
    type = ElementalVariableValue
    variable = clad_coolant_flux
    elementid = 2981
    use_displaced_mesh = 1
  []

  [coolant_channel_hmode] #COCE
    type = ElementalVariableValue
    variable = coolant_channel_hmode
    elementid = 2981
    use_displaced_mesh = 1
  []

  [coolant_channel_htype] #COCE
    type = ElementalVariableValue
    variable = coolant_channel_htype
    elementid = 2981
    use_displaced_mesh = 1
  []

  [critical_heat_flux] #COCE
    type = ElementalVariableValue
    variable = critical_heat_flux
    elementid = 2981
    use_displaced_mesh = 1
  []

  [fuel_centerline_temp] #FICN
    type = NodalVariableValue
    variable = temp
    nodeid = 3866
  []

  [fuel_surface_temp] #FOCN
    type = NodalVariableValue
    variable = temp
    nodeid = 3823
  []

  [clad_inner_surface_temp] #CICN
    type = NodalVariableValue
    variable = temp
    nodeid = 9557
  []

  [clad_outer_surface_temp] #COCN
    type = NodalVariableValue
    variable = temp
    nodeid = 9547
  []

  [fuel_axial_elongation] #FOTN
    type = NodalVariableValue
    variable = disp_y
    nodeid = 7739
  []

  [clad_radial_elongation] #COCN
    type = NodalVariableValue
    variable = disp_x
    nodeid = 9547
  []

  [fuel_radial_elongation] #FOCN
    type = NodalVariableValue
    variable = disp_x
    nodeid = 3823
  []

  [SED_PPN_O] #COCE
    type = ElementalVariableValue
    variable = SED
    elementid = 2981
    use_displaced_mesh = 1
  []

  [SED_PPN_I] #CICE
    type = ElementalVariableValue
    variable = SED
    elementid = 2984
    use_displaced_mesh = 1
  []

  [zz_OFract_PPN_O] #COCE
    type = ElementalVariableValue
    variable = oxywtfract_total
    elementid = 2981
    use_displaced_mesh = 1
  []

  [zz_OGain_PPN_O] #COCE
    type = ElementalVariableValue
    variable = oxywtfgain_total
    elementid = 2981
    use_displaced_mesh = 1
  []
  #######################################

  [max_clad_SED]
    type = ElementExtremeValue
    block = 'clad'
    variable = SED
    value_type = max
  []

  #Post processor to calculate radial average enthalpy. This postprocessor isnt available yet in BISON
  [z_RAE]
    type = RadialAverageEnthalpy
    vector_postprocessor = rad_temp
    radial_direction = x
    axial_direction = y
    axial_position = 0.3
    temperature_name = temp
  []
  [peak_RAE]
    type = TimeExtremeValue
    postprocessor = z_RAE
  []
[]

[VectorPostprocessors]
  [clad_dia]
    type = NodalValueSampler
    variable = disp_x
    boundary = 2
    sort_by = y
    outputs = 'outfile_clad_diameter'
  []
  [pellet_dia]
    type = NodalValueSampler
    variable = disp_x
    boundary = 10
    sort_by = y
    outputs = 'outfile_pellet_diameter'
  []

  #Location of peak power node at appoximately 0.3 m in mesh
  [rad_temp]
    type = NodalValueSampler
    block = 3
    sort_by = y
    variable = temp
    execute_on = timestep_end
    outputs = 'outfile_radial_temp'
  []
[]

[StandardLWRFuelRodOutputs]
  fuel_pellet_blocks = 'pellet_type_1 pellet_type_2 pellet_type_3'
  temperature = temp
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  time_step_interval =  1
  csv = true
  #exodus = true
  color = false
  [outfile_clad_diameter]
    type = CSV
    sync_times = '125194100'
    sync_only = true
  []
  [outfile_pellet_diameter]
    type = CSV
    sync_times = '125194100'
    sync_only = true
  []
  [outfile_radial_temp]
    type = CSV
    end_time = -100000
  []

  [console]
    type = Console
    output_linear = true
    max_rows = 10
  []
  [checkpoint]
    type = Checkpoint
    num_files = 2
    file_base = recover_files
  []
  [chkfile]
    type = CSV
    show = 'ave_temp_interior fis_gas_released FCT average_burnup peak_RAE'
    execute_on = 'FINAL'
  []
  [exodus]
    type = Exodus
    time_step_interval =  4
    end_time = 125193700
  []
  [exodus_RIA]
    type = Exodus
    time_step_interval =  3
    start_time = 125193695
  []
  [checkpoint_RIA]
    type = Checkpoint
    file_base = recover_files_RIA
    sync_times = '124761600 125107200 125193600 125193650 125193700 125193700.06 125193700.07 125193700.08 125193700.09 125193700.10 125193700.20 125193700.30 125193700.40 125193700.50 125193700.60 125193700.70 125193700.80 125193700.90 125193701.00 125193701.25 125193701.50 125193701.75 125193702.00 125193702.25 125193702.50 125193702.75 125193703.00 125193704.00 125193705.00'
    sync_only = true
  []
[]

[Debug]
  show_var_residual = 'disp_x disp_y temp'
  show_var_residual_norms = true
  show_material_props = true
[]

(examples/accident_tolerant_fuel/uo2_fecral/uo2_fecral.i)


initial_fuel_density = 10431.0

[GlobalParams]
  # Set initial fuel density, other global parameters
  density = ${initial_fuel_density}
  order = SECOND
  family = LAGRANGE
  energy_per_fission = 3.2e-11 # J/fission
  displacements = 'disp_x disp_y'
[]

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

[Mesh]
  coord_type = RZ
  displacements = 'disp_x disp_y'
  patch_size = 10 # For contact algorithm
  patch_update_strategy = auto
  partitioner = centroid
  centroid_partitioner_direction = y
  [mesh]
    type = FileMeshGenerator
    file = uo2_fecral_smeared.e
  []
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [temp]
    initial_condition = 293.0
  []
[]

[UserObjects]
  [pin_geometry]
    type = FuelPinGeometry
    clad_inner_wall = 5
    clad_outer_wall = 2
    clad_top = 3
    clad_bottom = 1
    pellet_exteriors = 8
  []
[]

[AuxVariables]
  [fast_neutron_flux]
    block = clad
  []
  [fast_neutron_fluence]
    block = clad
  []
  [grain_radius]
    block = pellet_type_1
    initial_condition = 10e-6
  []
  [gap_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxide_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_hoop_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_hoop]
    order = CONSTANT
    family = MONOMIAL
  []
  [hoop_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_rate]
    order = CONSTANT
    family = MONOMIAL
  []
  [mass_gain]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 1e4 1e8'
    y = '0 2.5e4 2.5e4'
    scale_factor = 1
  []

  [axial_peaking_factors]
    type = ParsedFunction
    expression = 1
  []

  [pressure_ramp]
    type = PiecewiseLinear
    x = '-200 0 1e8'
    y = '6.537e-3 1 1'
    scale_factor = 15.5e6
  []

  [mass_flux_func]
    type = PiecewiseLinear
    x = '-200 0  1e8'
    y = '3800. 3800. 3800.'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [pellets]
    block = pellet_type_1
    strain = FINITE
    eigenstrain_names = 'fuel_relocation_strain fuel_thermal_strain fuel_volumetric_strain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
    extra_vector_tags = 'ref'
    temperature = temp
  []
  [clad]
    block = clad
    strain = FINITE
    eigenstrain_names = 'clad_thermal_eigenstrain clad_volumetric_strain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
    extra_vector_tags = 'ref'
    temperature = temp
  []
[]

[Kernels]
  [gravity]
    type = Gravity
    variable = disp_y
    value = -9.81
    extra_vector_tags = 'ref'
  []

  [heat]
    type = HeatConduction
    variable = temp
    extra_vector_tags = 'ref'
  []

  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temp
    extra_vector_tags = 'ref'
  []

  [heat_source]
    type = NeutronHeatSource
    variable = temp
    block = pellet_type_1
    burnup_function = burnup
    extra_vector_tags = 'ref'
  []
[]

[Burnup]
  [burnup]
    block = pellet_type_1
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    num_radial = 81
    num_axial = 11
    fuel_pin_geometry = pin_geometry
    fuel_volume_ratio = 1.0
    RPF = RPF
  []
[]

[AuxKernels]
  [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_type_1
    variable = grain_radius
    temperature = temp
    execute_on = linear
  []
  [hoop_stress]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = hoop_stress
    scalar_type = HoopStress
    execute_on = timestep_end
  []
  [total_hoop_strain]
    type = RankTwoScalarAux
    rank_two_tensor = total_strain
    variable = total_hoop_strain
    scalar_type = HoopStress
    execute_on = timestep_end
  []
  [creep_strain_hoop]
    type = RankTwoScalarAux
    rank_two_tensor = creep_strain
    variable = creep_strain_hoop
    scalar_type = HoopStress
    block = clad
    execute_on = timestep_end
  []
  [conductance]
    type = MaterialRealAux
    property = gap_conductance
    variable = gap_cond
    boundary = 10
  []
  [coolant_htc]
    type = MaterialRealAux
    property = coolant_channel_htc
    variable = coolant_htc
    boundary = 2
  []
  [creep_rate]
    type = MaterialRealAux
    variable = creep_rate
    property = creep_rate
    execute_on = timestep_end
    block = clad
  []
  [oxide]
    type = MaterialRealAux
    variable = oxide_thickness
    property = scale_thickness
    boundary = 2
  []
  [mass_gain]
    type = MaterialRealAux
    variable = mass_gain
    property = oxide_mass_gain
    boundary = 2
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    normal_smoothing_distance = 0.1
    penalty = 1e7
  []
[]

[ThermalContact]
  [thermal_contact]
    type = GasGapHeatTransfer
    variable = temp
    primary = 5
    secondary = 10
    initial_moles = initial_moles
    gas_released = fis_gas_released
    contact_pressure = contact_pressure
    quadrature = true
  []
[]

[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 = 2.0e6
      startup_time = 0
      R = 8.3143
      output_initial_moles = initial_moles
      temperature = ave_temp_interior
      volume = gas_volume
      material_input = fis_gas_released
      output = plenum_pressure
    []
  []

[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = '1 2 3'
    variable = temp
    inlet_temperature = 580 # K
    inlet_pressure = pressure_ramp # Pa
    inlet_massflux = mass_flux_func # kg/m^2-sec
    rod_diameter = 9.5e-3 # m
    rod_pitch = 1.26e-2 # m
    linear_heat_rate = power_history
    axial_power_profile = axial_peaking_factors
    oxide_thickness = oxide_thickness
  []
[]

[Materials]
  [fuel_thermal]
    type = UO2Thermal
    block = pellet_type_1
    thermal_conductivity_model = NFIR
    temperature = temp
    burnup_function = burnup
  []
  [fuel_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = pellet_type_1
    youngs_modulus = 2.0e11
    poissons_ratio = 0.345
  []
  [elastic_stress]
    type = ComputeSmearedCrackingStress
    block = pellet_type_1
    cracking_stress = 1.68e8
    inelastic_models = 'fuel_creep'
    softening_models = exponential_softening
    shear_retention_factor = 0.1
    max_stress_correction = 0
    cracked_elasticity_type = DIAGONAL
    output_properties = crack_damage
    outputs = exodus
  []
  [exponential_softening]
    type = ExponentialSoftening
  []
  [fuel_creep]
    type = UO2CreepUpdate
    block = pellet_type_1
    burnup_function = burnup
    temperature = temp
  []
  [fuel_relocation]
    type = UO2RelocationEigenstrain
    block = pellet_type_1
    burnup_function = burnup
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    fuel_pin_geometry = 'pin_geometry'
    relocation_activation1 = 5000
    relocation_model = ESCORE_modified
    eigenstrain_name = fuel_relocation_strain
  []
  [fuel_thermal_expansion]
    type = UO2ThermalExpansionMATPROEigenstrain
    block = pellet_type_1
    temperature = temp
    stress_free_temperature = 295.0
    eigenstrain_name = fuel_thermal_strain
  []
  [fuel_volumetric_swelling]
    type = UO2VolumetricSwellingEigenstrain
    gas_swelling_model_type = SIFGRS
    block = pellet_type_1
    temperature = temp
    burnup_function = burnup
    initial_fuel_density = 10431.0
    eigenstrain_name = fuel_volumetric_strain
  []
  [clad_thermal]
    type = FeCrAlThermal
    material = C35M
    block = clad
    temperature = temp
  []
  [clad_elasticity_tensor] # isotropic elasticity tensor for Zry cladding
    type = FeCrAlElasticityTensor
    temperature = temp
    fecral_material_type = C35M
    block = clad
  []
  [clad_stress] # stress update class to govern the return mapping algorithm for creep
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_creep clad_plasticity'
    block = clad
  []
  [clad_creep]
    type = FeCrAlCreepUpdate
    block = clad
    temperature = temp
    fecral_material_type = C35M
    fast_neutron_flux = fast_neutron_flux
    model_irradiation_creep = true
    model_thermal_creep = true
    max_inelastic_increment = 1e-4
  []
  [thermal_expansion]
    type = FeCrAlThermalExpansionEigenstrain
    block = clad
    temperature = temp
    fecral_material_type = C35M
    stress_free_temperature = 295.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [irradiation_swelling]
    type = FeCrAlVolumetricSwellingEigenstrain
    block = clad
    temperature = temp
    fast_neutron_fluence = fast_neutron_fluence
    eigenstrain_name = clad_volumetric_strain
  []
  [clad_plasticity]
    type = FeCrAlPlasticityUpdate
    block = clad
    hardening_constant = 2.5e9
    temperature = temp
    yield_stress = 500.0
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = pellet_type_1
    temperature = temp
    burnup_function = burnup
    grain_radius = grain_radius
    gbs_model = true
  []
  [clad_density]
    type = StrainAdjustedDensity
    block = clad
    strain_free_density = 7250.0
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = pellet_type_1
    strain_free_density = ${initial_fuel_density}
  []

  [failure_criterion]
    type = FeCrAlCladdingFailure
    boundary = '2 5'
    hoop_stress = hoop_stress
    failure_criterion = UTS
    temperature = temp
  []
  [oxidation]
    type = FeCrAlOxidation
    reactor_type = PWR
    boundary = 2
  []
[]

[Dampers]
  [limitT]
    type = BoundingValueNodalDamper
    max_value = 3200.0
    min_value = 293.0
    variable = temp
  []
  [limitX]
    type = MaxIncrement
    max_increment = 1e-5
    variable = disp_x
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

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

  line_search = 'none'

  l_max_its = 100
  l_tol = 8e-3

  nl_max_its = 25
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-10

  start_time = -200
  n_startup_steps = 1
  end_time = 1e8

  dtmax = 1e5
  dtmin = 1

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 2.0e2
    force_step_every_function_point = true
    timestep_limiting_function = power_history
    max_function_change = 5e5
    optimal_iterations = 10
    iteration_window = 2
    linear_iteration_ratio = 100
    growth_factor = 2.0
    timestep_limiting_postprocessor = material_timestep
  []

  [Quadrature]
    order = FIFTH
    side_order = SEVENTH
  []
[]

[Postprocessors]
  [ave_temp_interior]
    type = SideAverageValue
    boundary = 9
    variable = temp
    execute_on = 'initial linear'
  []

  [avg_clad_temp]
    type = SideAverageValue
    boundary = 7
    variable = temp
  []

  [fis_gas_produced]
    type = ElementIntegralFisGasGeneratedSifgrs
    block = pellet_type_1
  []

  [fis_gas_released]
    type = ElementIntegralFisGasReleasedSifgrs
    block = pellet_type_1
  []
  [fis_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    block = pellet_type_1
  []
  [fis_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    block = pellet_type_1
  []

  [gas_volume]
    type = InternalVolume
    boundary = 9
    execute_on = 'initial linear'
  []

  [_dt]
    type = TimestepSize
  []

  [num_lin_it]
    type = NumLinearIterations
  []
  [num_nonlin_it]
    type = NumNonlinearIterations
  []
  [tot_lin_it]
    type = CumulativeValuePostprocessor
    postprocessor = num_lin_it
  []
  [tot_nonlin_it]
    type = CumulativeValuePostprocessor
    postprocessor = num_nonlin_it
  []
  [alive_time]
    type = PerfGraphData
    section_name = Root
    data_type = TOTAL
  []

  [rod_total_power]
    type = ElementIntegralPower
    variable = temp
    burnup_function = burnup
    block = pellet_type_1
  []

  [alhr_input]
    type = FunctionValuePostprocessor
    function = power_history
  []
  [average_burnup]
    type = ElementAverageValue
    block = pellet_type_1
    variable = burnup
  []
  [oxide_thickness]
    type = ElementExtremeValue
    block = clad
    variable = oxide_thickness
  []
  [mass_gain]
    type = ElementExtremeValue
    block = clad
    variable = mass_gain
  []
  [fis_gas_percent]
    type = FGRPercent
    fission_gas_released = fis_gas_released
    fission_gas_generated = fis_gas_produced
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = clad
  []
[]

[Outputs]
  perf_graph = true
  time_step_interval =  1
  exodus = true
  csv = true
  print_linear_residuals = true
  color = false

  [console]
    type = Console
    max_rows = 25
  []
[]

(test/tests/solid_mechanics/uo2_eigenstrains/ad_uo2_matpro_cte.i)

#--------------------------------------------------------------------------------
#
# This test case is designed to test thermal expansion model of UO2 fuel in a
# single element 1x1 square.
#
# - Temperature ramps from 300 K to 1500 K over 2 seconds with a nominal timestep of
#    0.1 seconds.  The equation for thermal expansion of uo2, with no Pu content and
#    no molten fraction is:
#     thermal_strain = 1.0e-5 * temperature - 3.0e-3 + 4.0e-2 * exp (-6.9e-20/(bk * temperature)),
#    where bk is the Boltzmann constant and is equal to 1.38e-23
#
#  The results for this nominal temperature change are given below.
#
#  Temperature (K)    Thermal Strain from MATPRO (m/m)       BISON-Computed Thermal Strain (m/m)
#       240              -6.000e-04                             -6.000023e-04
#       300               ~0.0 (reference temperature)           6.840728e-19
#       360              6.0004E-04                              6.000349e-04
#       420              1.2003E-03                              1.200268e-03
#       480              1.8012E-03                              1.801195e-03
#       540              2.4038E-03                              2.403807e-03
#       600              3.0096E-03                              3.009612e-03
#       660              3.6205E-03                              3.620507e-03
#       720              4.2386E-03                              4.238557e-03
#       780              4.8658E-03                              4.865782e-03
#       840              5.5040E-03                              5.503983e-03
#       900              6.1546E-03                              6.154634e-03
#       960              6.8188E-03                              6.818829e-03
#      1020              7.4973E-03                              7.497278e-03
#      1080              8.1903E-03                              8.190333e-03
#      1140              8.8980E-03                              8.898032e-03
#      1200              9.6202E-03                              9.620152e-03
#      1260              1.0356E-02                              1.035625e-02
#      1320              1.1106E-02                              1.110574e-02
#      1380              1.1868E-02                              1.186789e-02
#      1440              1.2642E-02                              1.264192e-02
#      1500              1.3427E-02                              1.342696e-02
#--------------------------------------------------------------------------------

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  coord_type = RZ
  use_displaced_mesh = false
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
  []
[]

[AuxVariables]
  [temperature]
    order = FIRST
    family = LAGRANGE
    initial_condition = 180.0
  []
[]

[Functions]
  [temperaturefcn]
    type = ParsedFunction
    expression = 180+6.0e2*t
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    add_variables = true
    strain = SMALL
    eigenstrain_names = thermal_eigenstrain
    use_automatic_differentiation = true
  []
[]

[AuxKernels]
  [temperature]
    type = FunctionAux
    variable = temperature
    function = temperaturefcn
    execute_on = timestep_begin
  []
[]

[BCs]
  [u_bottom_fix]
    type = ADDirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  []
  [u_left_fix]
    type = ADDirichletBC
    variable = disp_x
    boundary = left
    value = 0.0
  []
[]

[Materials]
  [elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    block = 0
    youngs_modulus = 2e11
    poissons_ratio = 0.3
  []
  [thermal_strain]
    type = ADUO2ThermalExpansionMATPROEigenstrain
    block = 0
    temperature = temperature
    stress_free_temperature = 300
    eigenstrain_name = thermal_eigenstrain
  []
  [elastic_stress]
    type = ADComputeLinearElasticStress
    block = 0
  []
[]

[Executioner]
  type = Transient
  l_max_its = 100
  nl_max_its = 100
  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-10
  l_tol = 1e-6
  start_time = 0.0
  end_time = 2.2
  dt = 0.1
[]

(examples/accident_tolerant_fuel/uo2_coated_zircaloy/uo2_coated_zircaloy.i)


initial_fuel_density = 10431.0

[GlobalParams]
  # Set initial fuel density, other global parameters
  density = ${initial_fuel_density}
  order = SECOND
  family = LAGRANGE
  energy_per_fission = 3.2e-11 # J/fission
  displacements = 'disp_x disp_y'
[]

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

[Mesh]
  coord_type = RZ
  [smeared_pellet_mesh]
    type = FuelPinMeshGenerator
    clad_top_gap_height = 0.026
    pellet_height = 0.1186
    pellet_quantity = 1
    clad_bot_gap_height = 0
    pellet_outer_radius = 4.1e-3
    clad_gap_width = 80e-6
    clad_thickness = 0.57e-3
    coating_thickness = 40e-6
    clad_mesh_density = customize
    pellet_mesh_density = customize
    nx_c = 3
    ny_c = 40
    nx_p = 11
    ny_p = 40
    nx_coating = 2
    elem_type = QUAD8
  []
  patch_update_strategy = auto
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [temperature]
    initial_condition = 293.0
  []
[]

[UserObjects]
  [pin_geometry]
    type = FuelPinGeometry
    clad_inner_wall = 5
    clad_outer_wall = 2
    clad_top = 3
    clad_bottom = 1
    pellet_exteriors = 8
  []
[]

[AuxVariables]
  [fast_neutron_flux]
    block = clad
  []
  [fast_neutron_fluence]
    block = clad
  []
  [grain_radius]
    block = pellet
    initial_condition = 10e-6
  []
  [gap_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxide_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_hoop_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_hoop]
    order = CONSTANT
    family = MONOMIAL
  []
  [hoop_stress]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 1e4 1e8'
    y = '0 2.5e4 2.5e4'
    scale_factor = 1
  []

  [axial_peaking_factors]
    type = ParsedFunction
    expression = 1
  []

  [pressure_ramp]
    type = PiecewiseLinear
    x = '-200 0 1e8'
    y = '6.537e-3 1 1'
    scale_factor = 15.5e6
  []

  [mass_flux_func]
    type = PiecewiseLinear
    x = '-200 0  1e8'
    y = '3800. 3800. 3800.'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [pellets]
    block = pellet
    strain = FINITE
    eigenstrain_names = 'fuel_relocation_strain fuel_thermal_strain fuel_volumetric_strain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
    extra_vector_tags = 'ref'
  []
  [clad]
    block = clad
    strain = FINITE
    eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_strain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
    extra_vector_tags = 'ref'
  []
  [coating]
    block = coating
    strain = FINITE
    eigenstrain_names = 'coating_thermal_eigenstrain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
    extra_vector_tags = 'ref'
  []
[]

[Kernels]
  [gravity]
    type = Gravity
    variable = disp_y
    value = -9.81
    extra_vector_tags = 'ref'
  []

  [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
    block = pellet
    burnup_function = burnup
    extra_vector_tags = 'ref'
  []
[]

[Burnup]
  [burnup]
    block = pellet
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    num_radial = 81
    num_axial = 11
    fuel_pin_geometry = pin_geometry
    fuel_volume_ratio = 1.0
    RPF = RPF
  []
[]

[AuxKernels]
  [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
  []
  [hoop_stress]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = hoop_stress
    scalar_type = HoopStress
    execute_on = timestep_end
  []
  [total_hoop_strain]
    type = RankTwoScalarAux
    rank_two_tensor = total_strain
    variable = total_hoop_strain
    scalar_type = HoopStress
    execute_on = timestep_end
  []
  [creep_strain_hoop]
    type = RankTwoScalarAux
    rank_two_tensor = creep_strain
    variable = creep_strain_hoop
    scalar_type = HoopStress
    block = clad
    execute_on = timestep_end
  []
  [conductance]
    type = MaterialRealAux
    property = gap_conductance
    variable = gap_cond
    boundary = 10
  []
  [coolant_htc]
    type = MaterialRealAux
    property = coolant_channel_htc
    variable = coolant_htc
    boundary = 2
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    normal_smoothing_distance = 0.1
    penalty = 1e7
  []
[]

[ThermalContact]
  [thermal_contact]
    type = GasGapHeatTransfer
    variable = temperature
    primary = 5
    secondary = 10
    initial_moles = initial_moles
    gas_released = fis_gas_released
    contact_pressure = contact_pressure
    quadrature = true
  []
[]

[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 = 2.0e6
      startup_time = 0
      R = 8.3143
      output_initial_moles = initial_moles
      temperature = ave_temp_interior
      volume = gas_volume
      material_input = fis_gas_released
      output = plenum_pressure
    []
  []

[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = '1 2 3'
    variable = temperature
    inlet_temperature = 580
    inlet_pressure = pressure_ramp
    inlet_massflux = mass_flux_func
    rod_diameter = 9.54e-3
    rod_pitch = 1.26e-2
    linear_heat_rate = power_history
    axial_power_profile = axial_peaking_factors
    oxide_thickness = oxide_thickness
  []
[]

[Materials]
  # Fuel
  [fuel_thermal]
    type = UO2Thermal
    block = pellet
    thermal_conductivity_model = NFIR
    temperature = temperature
    burnup_function = burnup
  []
  [fuel_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = pellet
    youngs_modulus = 2.0e11
    poissons_ratio = 0.345
  []
  [elastic_stress]
    type = ComputeSmearedCrackingStress
    block = pellet
    cracking_stress = 1.68e8
    inelastic_models = 'fuel_creep'
    softening_models = exponential_softening
    shear_retention_factor = 0.1
    max_stress_correction = 0
    cracked_elasticity_type = DIAGONAL
    output_properties = crack_damage
    outputs = exodus
  []
  [exponential_softening]
    type = ExponentialSoftening
  []
  [fuel_creep]
    type = UO2CreepUpdate
    block = pellet
    burnup_function = burnup
    temperature = temperature
  []
  [fuel_relocation]
    type = UO2RelocationEigenstrain
    block = pellet
    burnup_function = burnup
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    fuel_pin_geometry = 'pin_geometry'
    relocation_activation1 = 5000
    relocation_model = ESCORE_modified
    eigenstrain_name = fuel_relocation_strain
  []
  [fuel_thermal_expansion]
    type = UO2ThermalExpansionMATPROEigenstrain
    block = pellet
    temperature = temperature
    stress_free_temperature = 295.0
    eigenstrain_name = fuel_thermal_strain
  []
  [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_strain
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = pellet
    temperature = temperature
    burnup_function = burnup
    grain_radius = grain_radius
    gbs_model = true
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = pellet
    strain_free_density = ${initial_fuel_density}
  []

  # Clad
  [clad_thermal]
    type = ZryThermal
    temperature = temperature
    block = clad
  []
  [clad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 7.5e10
    poissons_ratio = 0.3
    block = clad
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep clad_plasticity'
    block = clad
  []
  [clad_zrycreep]
    type = ZryCreepLimbackHoppeUpdate
    block = clad
    temperature = temperature
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    model_irradiation_creep = true
    model_primary_creep = true
    model_thermal_creep = true
    max_inelastic_increment = 1e-4
    zircaloy_material_type = stress_relief_annealed
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = clad
    temperature = temperature
    stress_free_temperature = 295.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [irradiation_swelling]
    type = ZryIrradiationGrowthEigenstrain
    block = clad
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = stress_relief_annealed
    eigenstrain_name = clad_irradiation_strain
  []
  [clad_plasticity]
    type = ZryPlasticityUpdate
    block = clad
    temperature = temperature
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    cold_work_factor = 0.5
    plasticity_model_type = MATPRO
    zircaloy_alloy_type = 4
  []
  [clad_density]
    type = StrainAdjustedDensity
    block = clad
    strain_free_density = 6511.0
  []

  # Coating
  [coat_thermal]
    type = ChromiumThermal
    block = coating
    temperature = temperature
  []
  [coating_elasticity_tensor]
    type = ChromiumElasticityTensor
    temperature = temperature
    block = coating
  []
  [coat_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'coat_creep coat_plasticity'
    block = coating
  []
  [coat_creep]
    type = ChromiumCreepUpdate
    temperature = temperature
    block = coating
  []
  [coat_plasticity]
    type = ChromiumPlasticityUpdate
    temperature = temperature
    block = coating
    fast_neutron_fluence = 0.0
    hardening_constant = 2e9
  []
  [coat_thermal_expansion]
    type = ChromiumThermalExpansionEigenstrain
    block = coating
    temperature = temperature
    stress_free_temperature = 293.0
    eigenstrain_name = coating_thermal_eigenstrain
  []
  [density_coat]
    type = StrainAdjustedDensity
    block = coating
    strain_free_density = 7190.0
  []
[]

[Dampers]
  [limitT]
    type = BoundingValueNodalDamper
    max_value = 3200.0
    min_value = 293.0
    variable = temperature
  []
  [limitX]
    type = MaxIncrement
    max_increment = 1e-5
    variable = disp_x
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

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

  line_search = 'none'

  l_max_its = 100
  l_tol = 8e-3

  nl_max_its = 25
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-10

  start_time = -200
  n_startup_steps = 1
  end_time = 5e7

  dtmax = 1e6
  dtmin = 1

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 2.0e2
    force_step_every_function_point = true
    timestep_limiting_function = power_history
    max_function_change = 5e5
    optimal_iterations = 10
    iteration_window = 2
    linear_iteration_ratio = 100
    growth_factor = 2.0
    timestep_limiting_postprocessor = material_timestep
  []

  [Quadrature]
    order = FIFTH
    side_order = SEVENTH
  []
[]

[Postprocessors]
  [ave_temp_interior]
    type = SideAverageValue
    boundary = 9
    variable = temperature
    execute_on = 'initial linear'
  []

  [avg_clad_temp]
    type = SideAverageValue
    boundary = 7
    variable = temperature
  []

  [fis_gas_produced]
    type = ElementIntegralFisGasGeneratedSifgrs
    block = pellet
  []

  [fis_gas_released]
    type = ElementIntegralFisGasReleasedSifgrs
    block = pellet
  []
  [fis_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    block = pellet
  []
  [fis_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    block = pellet
  []

  [gas_volume]
    type = InternalVolume
    boundary = 9
    execute_on = 'initial linear'
  []

  [_dt]
    type = TimestepSize
  []

  [num_lin_it]
    type = NumLinearIterations
  []
  [num_nonlin_it]
    type = NumNonlinearIterations
  []
  [tot_lin_it]
    type = CumulativeValuePostprocessor
    postprocessor = num_lin_it
  []
  [tot_nonlin_it]
    type = CumulativeValuePostprocessor
    postprocessor = num_nonlin_it
  []
  [alive_time]
    type = PerfGraphData
    section_name = Root
    data_type = TOTAL
  []

  [rod_total_power]
    type = ElementIntegralPower
    variable = temperature
    burnup_function = burnup
    block = pellet
  []

  [alhr_input]
    type = FunctionValuePostprocessor
    function = power_history
  []
  [average_burnup]
    type = ElementAverageValue
    block = pellet
    variable = burnup
  []
  [oxide_thickness]
    type = ElementExtremeValue
    block = clad
    variable = oxide_thickness
  []
  [fis_gas_percent]
    type = FGRPercent
    fission_gas_released = fis_gas_released
    fission_gas_generated = fis_gas_produced
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = clad
  []
[]

[Outputs]
  perf_graph = true
  time_step_interval =  1
  exodus = true
  csv = true
  print_linear_residuals = true
  color = false

  [console]
    type = Console
    max_rows = 25
  []
[]

(assessment/LWR/validation/Super_Ramp/analysis/Super_Ramp_Base.i)

#This is a partial input file base with information/features common to several experiments within this assessment
#NOTE: This file will NOT run on its own, it requires a PK#X.params file and a PK##.params file to run

# physical constants
R = 8.3143 # (J/K*mol) -- THIS SHOULD BE EDITED TO USE PHYSICALCONSTANTS' VALUE

# fuel isotope fractions and fission energy
energy_per_fission = 3.28451e-11 # (J/fission)
isotope_fraction_Pu239 = 0.0
isotope_fraction_Pu240 = 0.0
isotope_fraction_Pu241 = 0.0
isotope_fraction_Pu242 = 0.0

# rod geometry
clad_bot_gap_height = 1.0e-3 # (m)

# variable and kernel initial values
initial_temperature = 293.15 # (K)
gravity_constant = -9.81 # (m/s^2)

# fuel/cladding contact
friction_coefficient = 0.4
c_normal = 1e+12
c_tangential = 1e+24
tangential_lm_scaling = 1.0e-16
normal_lm_scaling = 1.0e-10
roughness_secondary = 1.0e-6 # (m)
roughness_primary = 2.0e-6 # (m)
roughness_coef = 3.2
jump_distance_model = LANNING
contact_pressure = mechanical_normal_lm
thermal_lm_scaling = 1.0e-2

# plenum parameters
initial_pressure = 2.25e6 # (Pa)
startup_time = 0 # (s)

# fuel/clad material properties
fuel_cracking_stress = 1.68e8
stress_free_temperature = 293.15 # (K)
clad_density = 6550.0 # (kg/m^3)

# numerical options
l_max_its = 100
l_tol = 1e-4
nl_max_its = 30
start_time = 0.0 # s
dtmin = 1.0 # s


[GlobalParams]
  density = ${fuel_density}
  displacements = 'disp_x disp_y'
  order = SECOND
  family = LAGRANGE
  energy_per_fission = ${energy_per_fission}
  volumetric_locking_correction = false
  initial_porosity = ${initial_fuel_porosity}
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
  converge_on = 'disp_x disp_y temperature'
[]

[Mesh]
  coord_type = RZ
  [smeared_pellet_mesh]
    type = FuelPinMeshGenerator
    pellet_quantity = ${number_pellets}
    pellet_height = ${pellet_height}
    pellet_outer_radius = ${pellet_outer_radius}
    pellet_mesh_density = customize
    nx_p = 11
    ny_p = 84
    clad_mesh_density = customize
    nx_c = 4
    ny_c = 84
    clad_gap_width = ${clad_gap_width}
    clad_thickness = ${clad_thickness}
    clad_bot_gap_height = ${clad_bot_gap_height}
    bottom_clad_height = ${bottom_clad_height}
    top_clad_height = ${top_clad_height}
    clad_top_gap_height = ${clad_top_gap_height}
    ny_cu = 3
    ny_cl = 3
    elem_type = QUAD8
  []
  patch_size = 10
  patch_update_strategy = iteration
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[Variables]
  [temperature]
    initial_condition = ${initial_temperature}
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    data_file = ${power_history_data_file}
    format = columns
  []
  [axial_power_factors]
    type = PiecewiseBilinear
    data_file = ${axial_power_factors_data_file}
    axis = 1
  []
  [clad_out_temp]
    type = PiecewiseLinear
    data_file = ${clad_out_temp_data_file}
    format = columns
  []
  [axial_temp_factors]
    type = PiecewiseBilinear
    data_file = ${axial_temp_factors_data_file}
    axis = 1
  []
  [clad_temp_bc]
    type = CompositeFunction
    functions = 'clad_out_temp axial_temp_factors'
  []
  [coolant_pressure]
    type = PiecewiseLinear
    data_file = ${coolant_pressure_data_file}
    format = columns
  []
  [fast_flux]
    type = PiecewiseLinear
    data_file = ${fast_neutron_flux_data_file}
    format = columns
  []
[]

[AuxVariables]
  [gap_conductance]
    block = 'mechanical_secondary_subdomain'
  []
  [fast_neutron_flux]
    block = clad
  []
  [fast_neutron_fluence]
    block = clad
  []
  [grain_radius]
    block = pellet
    initial_condition = ${initial_grain_radius}
  []
  [porosity]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
    initial_condition = ${initial_fuel_porosity}
  []
  [pellet_id]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [oxide_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
    block = clad
  []
  [gas_gen_3]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [gas_grn_3]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [gas_bdr_3]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [gas_rel_3]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [bbl_bdr_2]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [prs_bbl_bdr]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [prseq_bbl_bdr]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [rad_bbl_bdr]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [GBCoverage]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [sat_coverage]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [eff_diff_coeff]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [deltav_v0_bd]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [gaseous_porosity]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [thermal_conductivity]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [layered_average_contact_pressure]
    order = CONSTANT
    family = MONOMIAL
  []
  [volumetric_swelling_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [gap_cond]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [pellets]
    block = pellet
    temperature = temperature
    add_variables = true
    strain = FINITE
    eigenstrain_names = 'fuel_relocation_strain fuel_thermal_strain
      fuel_volumetric_strain'
    generate_output = 'vonmises_stress hydrostatic_stress stress_xx stress_yy
      stress_zz strain_xx strain_yy strain_zz'
    extra_vector_tags = 'ref'
  []
  [clad]
    block = clad
    temperature = temperature
    add_variables = true
    strain = FINITE
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_strain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz
      creep_strain_xx creep_strain_yy creep_strain_xy creep_strain_zz strain_xx
      strain_yy strain_zz hoop_stress'
    extra_vector_tags = 'ref'
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
    block = '1 3'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
    block = '1 3'
  []
  [heat_source]
    type = NeutronHeatSource
    variable = temperature
    block = pellet
    burnup_function = burnup
    extra_vector_tags = 'ref'
  []
  [gravity]
    type = Gravity
    variable = disp_y
    value = ${gravity_constant}
    block = '1 3'
  []
[]

[AuxKernels]
  [gap_conductance]
    type = GapConductanceMortar
    primary_boundary = 5
    secondary_boundary = 10
    primary_subdomain = 'mechanical_primary_subdomain'
    secondary_subdomain = 'mechanical_secondary_subdomain'
    heat_flux = thermal_contact_thermal_lm
    temperature = temperature
    variable = gap_conductance
  []
  [fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    block = clad
    function = fast_flux
    execute_on = timestep_begin
  []
  [fast_neutron_fluence]
    type = FastNeutronFluenceAux
    variable = fast_neutron_fluence
    fast_neutron_flux = fast_neutron_flux
    execute_on = timestep_begin
  []
  [grain_radius]
    type = GrainRadiusAux
    block = pellet
    variable = grain_radius
    temperature = temperature
    execute_on = linear
  []
  [porosity]
    type = PorosityAuxUO2
    block = pellet
    variable = porosity
    execute_on = linear
  []
  [pelletid]
    type = PelletIdAux
    block = pellet
    variable = pellet_id
    fuel_pin_geometry = pin_geometry
    number_pellets = ${number_pellets}
    execute_on = initial
  []
  [oxi_thickness]
    type = MaterialRealAux
    variable = oxide_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [creep_strain_mag]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fggen]
    type = MaterialRealAux
    variable = gas_gen_3
    property = gas_concentration_generated_total
  []
  [fggrn]
    type = MaterialRealAux
    variable = gas_grn_3
    property = gas_concentration_intra_total
  []
  [fgbdr]
    type = MaterialRealAux
    variable = gas_bdr_3
    property = gas_concentration_GB_bubble_volume
  []
  [fgrel]
    type = MaterialRealAux
    variable = gas_rel_3
    property = gas_concentration_release_total
  []
  [nbbl2]
    type = MaterialRealAux
    variable = bbl_bdr_2
    property = bubble_GB_surface_density
  []
  [prsbbl]
    type = MaterialRealAux
    variable = prs_bbl_bdr
    property = bubble_GB_pressure
  []
  [prseqbbl]
    type = MaterialRealAux
    variable = prseq_bbl_bdr
    property = bubble_GB_pressure_equilibrium
  []
  [radbbl]
    type = MaterialRealAux
    variable = rad_bbl_bdr
    property = bubble_radius_GB
  []
  [frcvrg]
    type = MaterialRealAux
    variable = GBCoverage
    property = GBCoverage
  []
  [stcvrg]
    type = MaterialRealAux
    variable = sat_coverage
    property = sat_coverage
  []
  [diffc]
    type = MaterialRealAux
    variable = eff_diff_coeff
    property = eff_diff_coeff
  []
  [dvv0bd]
    type = MaterialRealAux
    variable = deltav_v0_bd
    property = deltav_v0_bubble_GB
  []
  [gaspor]
    type = MaterialRealAux
    variable = gaseous_porosity
    property = gaseous_porosity
  []
  [fuel_conductivity]
    type = MaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
  []
  [layered_average_contact_pressure]
    type = SpatialUserObjectAux
    block = pellet
    variable = layered_average_contact_pressure
    execute_on = nonlinear
    user_object = layered_average_contact_pressure
  []
  [volumetric_swelling_strain]
    type = MaterialRealAux
    variable = volumetric_swelling_strain
    property = volumetric_swelling_strain
    block = pellet
    execute_on = nonlinear
  []
[]

[Burnup]
  [burnup]
    block = pellet
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_power_factors
    num_radial = 80
    num_axial = 20
    isotopes = 'U235 U238 Pu239 Pu240 Pu241 Pu242'
    isotope_fractions = '${isotope_fraction_U235} ${isotope_fraction_U238} ${isotope_fraction_Pu239} ${isotope_fraction_Pu240} ${isotope_fraction_Pu241} ${isotope_fraction_Pu242}'
    RPF = RPF
    fuel_pin_geometry = pin_geometry
    fuel_volume_ratio = 1.0
  []
[]

[Contact]
  [mechanical]
    model = coulomb
    formulation = mortar
    primary = 5
    secondary = 10
    friction_coefficient = ${friction_coefficient}
    c_normal = ${c_normal}
    c_tangential = ${c_tangential}
    tangential_lm_scaling = ${tangential_lm_scaling}
    normal_lm_scaling = ${normal_lm_scaling}
  []
[]

[ThermalContactMortar]
  [thermal_contact]
    secondary_variable = temperature
    primary_boundary = '5'
    secondary_boundary = '10'
    initial_moles = initial_moles
    gas_released = fission_gas_released
    roughness_secondary = ${roughness_secondary}
    roughness_primary = ${roughness_primary}
    roughness_coef = ${roughness_coef}
    plenum_pressure = plenum_pressure
    jump_distance_model = ${jump_distance_model}
    contact_pressure = ${contact_pressure}
    thermal_lm_scaling = ${thermal_lm_scaling}
  []
[]

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

[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
  []
  [temp]
    type = FunctionDirichletBC
    boundary = '1 2 3'
    variable = temperature
    function = clad_temp_bc
  []
  [Pressure]
    [coolantPressure]
      boundary = '1 2 3'
      function = coolant_pressure
    []
  []
  [PlenumPressure]
    [plenumPressure]
      boundary = 9
      initial_pressure = ${initial_pressure}
      startup_time = ${startup_time}
      R = ${R}
      output_initial_moles = initial_moles
      temperature = plenum_temp
      volume = plenum_volume
      material_input = fission_gas_released
      output = plenum_pressure
    []
  []
[]

[Materials]
  [fuel_thermal]
    type = UO2Thermal
    block = pellet
    thermal_conductivity_model = NFIR
    temperature = temperature
    burnup_function = burnup
  []
  [fuel_volumetric_swelling]
    type = UO2VolumetricSwellingEigenstrain
    gas_swelling_model_type = SIFGRS
    block = pellet
    temperature = temperature
    burnup_function = burnup
    initial_fuel_density = ${fuel_density}
    eigenstrain_name = fuel_volumetric_strain
  []
  [fuel_elasticity_tensor]
    type = UO2ElasticityTensor
    block = pellet
    density = ${fuel_density}
    temperature = temperature
    matpro_youngs_modulus = true
    matpro_poissons_ratio = true
  []
  [fuel_creep]
    type = UO2CreepUpdate
    block = pellet
    temperature = temperature
    burnup_function = burnup
    initial_grain_radius = ${initial_grain_radius}
  []
  [fuel_stress]
    type = ComputeSmearedCrackingStress
    block = pellet
    cracking_stress = ${fuel_cracking_stress}
    inelastic_models = 'fuel_creep'
    softening_models = exponential_softening
    shear_retention_factor = 0.1
    max_stress_correction = 0
    cracked_elasticity_type = DIAGONAL
    output_properties = crack_damage
    outputs = exodus
  []
  [exponential_softening]
    type = ExponentialSoftening
  []
  [fuel_relocation]
    type = UO2RelocationEigenstrain
    block = pellet
    burnup_function = burnup
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_power_factors
    relocation_activation1 = 5000
    relocation_model = ESCORE_modified
    eigenstrain_name = fuel_relocation_strain
    fuel_pin_geometry = pin_geometry
    model_relocation_recovery = true
    max_relocation_recovery_fraction = 0.5
    relocation_scaling_factor = 1
    volumetric_swelling_increment = vol_swell_increment
    layered_average_contact_pressure = layered_average_contact_pressure
    outputs = all
    output_properties = 'relocation_strain recovered_relocation_strain'
  []
  [fuel_thermal_expansion]
    type = UO2ThermalExpansionMATPROEigenstrain
    block = pellet
    temperature = temperature
    stress_free_temperature = ${stress_free_temperature}
    eigenstrain_name = fuel_thermal_strain
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = pellet
    hydrostatic_stress = hydrostatic_stress
    diff_coeff_option = TURNBULL_D1_4D2_D3
    transient_option = MICROCRACKING_BURNUP
    res_param_option = HETEROGENEOUS_WHITE
    ig_bubble_model = NUCLEATION_RESOLUTION
    ig_diff_algorithm = polypole2
    temperature = temperature
    burnup_function = burnup
    grain_radius = grain_radius
    gbs_model = true
    pellet_id = pellet_id
    pellet_brittle_zone = pbz
    ath_model = true
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_power_factors
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = pellet
    strain_free_density = ${fuel_density}
  []

  [clad_thermal]
    type = ZryThermal
    block = clad
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ZryElasticityTensor
    block = clad
    fast_neutron_fluence = fast_neutron_fluence
    temperature = temperature
    matpro_youngs_modulus = true
    matpro_poissons_ratio = true
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = nonlinear
    inelastic_models = 'clad_zrycreep'
    block = clad
  []
  [clad_zrycreep]
    type = ZryCreepLimbackHoppeUpdate
    block = clad
    temperature = temperature
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    model_irradiation_creep = true
    model_primary_creep = true
    model_thermal_creep = true
    zircaloy_material_type = stress_relief_annealed
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = clad
    temperature = temperature
    stress_free_temperature = ${stress_free_temperature}
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [irradiation_swelling]
    type = ZryIrradiationGrowthEigenstrain
    block = clad
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = stress_relief_annealed
    eigenstrain_name = clad_irradiation_strain
  []
  [clad_density]
    type = StrainAdjustedDensity
    block = clad
    strain_free_density = ${clad_density}
  []
  [clad_oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    fuel_pin_geometry = pin_geometry
  []
[]

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

[Executioner]
  type = Transient

  petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -mat_mffd_err -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = ${petsc_options_value}
  line_search = 'none'
  verbose = true

  l_max_its = ${l_max_its}
  l_tol = ${l_tol}

  nl_max_its = ${nl_max_its}
  nl_rel_tol = ${nl_rel_tol}
  nl_abs_tol = ${nl_abs_tol}

  start_time = ${start_time}
  n_startup_steps = 1
  end_time = ${end_time}
  dtmax = ${dtmax}
  dtmin = ${dtmin}

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 100
    optimal_iterations = ${optimal_iterations}
    iteration_window = ${iteration_window}
    timestep_limiting_function = power_history
    force_step_every_function_point = true
  []
[]

[Postprocessors]
  [avg_gap_conductance]
    type = ElementAverageValue
    block = 'mechanical_secondary_subdomain'
    variable = gap_conductance
    execute_on = 'initial timestep_end'
  []
  [central_fuel_temp]
    type = NodalVariableValue
    variable = temperature
    nodeid = 3110 #Mesh dependent!
    execute_on = 'initial timestep_end'
  []
  [max_fuel_temp]
    type = NodalExtremeValue
    block = pellet
    value_type = max
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = clad
    value_type = max
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [midplane_hoop_strain_inner_clad]
    type = ElementalVariableValue
    elementid = 209 # Mesh dependent
    variable = strain_zz
    execute_on = 'initial timestep_end'
  []
  [midplane_hoop_stress_inner_clad]
    type = ElementalVariableValue
    elementid = 209 # Mesh dependent
    variable = hoop_stress
    execute_on = 'initial timestep_end'
  []
  [midplane_contact_pressure]
    type = NodalVariableValue
    variable = contact_pressure
    nodeid = 3141 # Mesh dependent
    execute_on = 'initial timestep_end'
  []
  [total_rod_integral_power]
    type = ElementIntegralPower
    variable = temperature
    burnup_function = burnup
    block = pellet
  []
  [total_rod_input_power]
    type = FunctionValuePostprocessor
    function = power_history
    scale_factor = ${scale_factor}
  []
  [vol_swell_increment]
    type = SideAverageIncrementTensorComponent
    boundary = 10
    variable = volumetric_swelling_strain
    execute_on = nonlinear
  []
  [midplane_clad_outer_temp]
    type = NodalVariableValue
    nodeid = 757 # Mesh dependent
    variable = temperature
  []
  [midplane_clad_inner_temp]
    type = NodalVariableValue
    nodeid = 747 # Mesh dependent
    variable = temperature
  []
[]

[StandardLWRFuelRodOutputs]
  fuel_pellet_blocks = 3
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  color = false
  [csv]
    type = CSV
    file_base = '${id}_csv'
  []
  [exodus]
    type = Exodus
    file_base = '${id}_exodus'
  []
  [console]
    type = Console
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = '${id}_chkfile'
    show = 'average_burnup fission_gas_released_percentage central_fuel_temp midplane_contact_pressure'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [pin_geometry]
    type = FuelPinGeometry
  []
  [pbz]
    type = PelletBrittleZone
    block = pellet
    pellet_id = pellet_id
    temperature = temperature
    fuel_pin_geometry = pin_geometry
    number_pellets = ${number_pellets}
    execute_on = 'initial linear'
  []
  [layered_average_contact_pressure]
    type = LayeredSideAverage
    variable = contact_pressure
    direction = y
    num_layers = 1
    execute_on = timestep_end
    boundary = 10
  []
[]

(assessment/LWR/validation/Super_Ramp/analysis/PK62/PK62_weighted_gap_VCP.i)

# This assessment case requires increasing the default automatic differentiation
# vector (AD size). At the time of writing, the default size is 50. This test
# requires, at least, 56. The minimum size required for a test to run can be
# specified in the 'assesment' file via "min_ad_size". To configure MOOSE
# to increase the AD vector size, one can use the command below in the MOOSE
# folder to expand this size, e.g., to 200:
#
# ./configure --with-derivative-type=sparse --with-ad-indexing-type=global --with-derivative-size=200

# physical constants
R = 8.3143 # J/mol*K -- THIS SHOULD BE EDITED TO USE PHYSICALCONSTANTS' VALUE

# fuel isotope fractions and fission energy
energy_per_fission = 3.28451e-11 # J/fission
isotope_fraction_U235 = 0.02985
isotope_fraction_U238 = 0.97015
isotope_fraction_Pu239 = 0.0
isotope_fraction_Pu240 = 0.0
isotope_fraction_Pu241 = 0.0
isotope_fraction_Pu242 = 0.0

# rod geometry
pellet_quantity = 29 # (-)
pellet_height = 0.010862 # m
pellet_outer_radius = 4.57e-3 # m
clad_gap_width = 75.0e-6 # m
clad_thickness = 0.725e-3 # m
clad_bot_gap_height = 1.0e-3 # m
bottom_clad_height = 15.5e-3 # m
top_clad_height = 15.5e-3 # m
clad_top_gap_height = 32.5e-3 # m
fuel_volume_ratio = 1.0 # (-)
rod_input_power_scale_factor = 0.376004 # m (rod height)

# variable and kernel initial values
initial_temperature = 293.15 # K
gravity_constant = -9.81 # m/s^2

# fuel/cladding contact
c_normal = 1e+05 # (-)
jump_distance_model = 'LANNING'
roughness_primary = 2e-6 # (-)
roughness_secondary = 1e-6 # (-)
roughness_coef = 3.2 # (-)
relocation_activation1 = 5000 # W/m
max_relocation_recovery_fraction = 0.5 # (-)
relocation_scaling_factor = 1 # (-)

# plenum parameters
initial_plenum_pressure = 2.25e6 # Pa
startup_time = 0 # s

# fuel/clad material properties
initial_fuel_density = 10420.0 # kg/m^3
initial_fuel_porosity = 0.049 # (-)
initial_grain_radius = 17.16e-6 # m
fuel_cracking_stress = 1.68e8 # Pa
fuel_shear_retention_factor = 0.1 # (-)
fuel_max_stress_correction = 0 # (-)
stress_free_temperature = 293.15 # K
cladding_density = 6550.0 # kg/m^3

# numerical options
damper_max_temperature_increment = 150.0 # K (was 30 K)
l_max_its = 35
l_tol = 1e-5
nl_max_its = 30
nl_rel_tol = 2e-6
nl_abs_tol = 1e-5
start_time = 0.0 # s
n_startup_steps = 1
dtmax = 1.0e6 # s
dtmin = 1.0 # s
Timestepper_dt = 100 # s
Timestepper_optimal_iterations = 15
Timestepper_iteration_window = 3

# irradiation history
end_time = 76630068.0 # s

# data files
power_history_data_file = 'alhr_history.csv'
axial_power_factors_data_file = 'axial_lhr_factors.csv'
clad_out_temp_data_file = 'temp_outer_clad_history.csv'
axial_temp_factors_data_file = 'axial_temp_factors.csv'
coolant_pressure_data_file = 'pressure_coolant.csv'
fast_neutron_flux_data_file = 'fast_neutron_flux.csv'

[GlobalParams]
  density = ${initial_fuel_density}
  displacements = 'disp_x disp_y'
  order = FIRST
  family = LAGRANGE
  energy_per_fission = ${energy_per_fission}
  volumetric_locking_correction = true
  initial_porosity = ${initial_fuel_porosity}
[]

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

[Mesh]
  coord_type = RZ
  patch_size = 80
  patch_update_strategy = iteration
  [smeared_pellet_mesh]
    type = FuelPinMeshGenerator
    pellet_quantity = ${pellet_quantity}
    pellet_height = ${pellet_height}
    pellet_outer_radius = ${pellet_outer_radius}
    pellet_mesh_density = customize
    nx_p = 11
    ny_p = 168
    clad_mesh_density = customize
    nx_c = 4
    ny_c = 84
    clad_gap_width = ${clad_gap_width}
    clad_thickness = ${clad_thickness}
    clad_bot_gap_height = ${clad_bot_gap_height}
    bottom_clad_height = ${bottom_clad_height}
    top_clad_height = ${top_clad_height}
    clad_top_gap_height = ${clad_top_gap_height}
    ny_cu = 3
    ny_cl = 3
    elem_type = QUAD4
  []
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[Variables]
  [temperature]
    initial_condition = ${initial_temperature}
    block = '1 3'
  []
  [disp_x]
    order = FIRST
    family = LAGRANGE
    block = '1 3'
    scaling = 1.0e3
  []

  [disp_y]
    order = FIRST
    family = LAGRANGE
    block = '1 3'
    scaling = 1.0e3
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    data_file = ${power_history_data_file}
    format = columns
  []
  [axial_power_factors]
    type = PiecewiseBilinear
    data_file = ${axial_power_factors_data_file}
    axis = 1
  []
  [clad_out_temp]
    type = PiecewiseLinear
    data_file = ${clad_out_temp_data_file}
    format = columns
  []
  [axial_temp_factors]
    type = PiecewiseBilinear
    data_file = ${axial_temp_factors_data_file}
    axis = 1
  []
  [clad_temp_bc]
    type = CompositeFunction
    functions = 'clad_out_temp axial_temp_factors'
  []
  [coolant_pressure]
    type = PiecewiseLinear
    data_file = ${coolant_pressure_data_file}
    format = columns
  []
  [fast_flux]
    type = PiecewiseLinear
    data_file = ${fast_neutron_flux_data_file}
    format = columns
  []
  [axial_power_constant]
    type = ConstantFunction
    value = 1
  []
[]

[AuxVariables]
  [fast_neutron_flux]
    block = clad
  []
  [fast_neutron_fluence]
    block = clad
  []
  [grain_radius]
    block = pellet
    initial_condition = ${initial_grain_radius}
  []
  [porosity]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
    initial_condition = ${initial_fuel_porosity}
  []
  [pellet_id]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [oxide_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
    block = clad
  []
  [gas_gen_3]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [gas_grn_3]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [gas_bdr_3]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [gas_rel_3]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [bbl_bdr_2]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [prs_bbl_bdr]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [prseq_bbl_bdr]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [rad_bbl_bdr]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [GBCoverage]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [sat_coverage]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [eff_diff_coeff]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [deltav_v0_bd]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [gaseous_porosity]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [thermal_conductivity]
    order = CONSTANT
    family = MONOMIAL
    block = pellet
  []
  [gap_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [layered_average_contact_pressure]
    order = CONSTANT
    family = MONOMIAL
  []
  [volumetric_swelling_strain]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [pellets]
    block = pellet
    temperature = temperature
    strain = FINITE
    eigenstrain_names = 'fuel_relocation_strain fuel_thermal_strain
      fuel_volumetric_strain'
    generate_output = 'vonmises_stress hydrostatic_stress stress_xx stress_yy
      stress_zz strain_xx strain_yy strain_zz'
    extra_vector_tags = 'ref'
  []
  [clad]
    block = clad
    temperature = temperature
    strain = FINITE
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_strain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz
      creep_strain_xx creep_strain_yy creep_strain_xy creep_strain_zz strain_xx
      strain_yy strain_zz hoop_stress'
    extra_vector_tags = 'ref'
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
    block = '1 3'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
    block = '1 3'
  []
  [heat_source]
    type = NeutronHeatSource
    variable = temperature
    block = pellet
    burnup_function = burnup
    extra_vector_tags = 'ref'
  []
  [gravity]
    type = Gravity
    variable = disp_y
    value = ${gravity_constant}
  []
[]

[AuxKernels]
  [fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    block = clad
    function = fast_flux
    execute_on = timestep_begin
  []
  [fast_neutron_fluence]
    block = clad
    type = FastNeutronFluenceAux
    variable = fast_neutron_fluence
    fast_neutron_flux = fast_neutron_flux
    execute_on = timestep_begin
  []
  [grain_radius]
    type = GrainRadiusAux
    block = pellet
    variable = grain_radius
    temperature = temperature
    execute_on = linear
  []
  [porosity]
    type = PorosityAuxUO2
    block = pellet
    variable = porosity
    execute_on = linear
  []
  [pelletid]
    type = PelletIdAux
    block = pellet
    variable = pellet_id
    fuel_pin_geometry = pin_geometry
    number_pellets = ${pellet_quantity}
    execute_on = initial
  []
  [oxi_thickness]
    type = MaterialRealAux
    variable = oxide_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [creep_strain_mag]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fggen]
    type = MaterialRealAux
    variable = gas_gen_3
    property = gas_concentration_generated_total
  []
  [fggrn]
    type = MaterialRealAux
    variable = gas_grn_3
    property = gas_concentration_intra_total
  []
  [fgbdr]
    type = MaterialRealAux
    variable = gas_bdr_3
    property = gas_concentration_GB_bubble_volume
  []
  [fgrel]
    type = MaterialRealAux
    variable = gas_rel_3
    property = gas_concentration_release_total
  []
  [nbbl2]
    type = MaterialRealAux
    variable = bbl_bdr_2
    property = bubble_GB_surface_density
  []
  [prsbbl]
    type = MaterialRealAux
    variable = prs_bbl_bdr
    property = bubble_GB_pressure
  []
  [prseqbbl]
    type = MaterialRealAux
    variable = prseq_bbl_bdr
    property = bubble_GB_pressure_equilibrium
  []
  [radbbl]
    type = MaterialRealAux
    variable = rad_bbl_bdr
    property = bubble_radius_GB
  []
  [frcvrg]
    type = MaterialRealAux
    variable = GBCoverage
    property = GBCoverage
  []
  [stcvrg]
    type = MaterialRealAux
    variable = sat_coverage
    property = sat_coverage
  []
  [diffc]
    type = MaterialRealAux
    variable = eff_diff_coeff
    property = eff_diff_coeff
  []
  [dvv0bd]
    type = MaterialRealAux
    variable = deltav_v0_bd
    property = deltav_v0_bubble_GB
  []
  [gaspor]
    type = MaterialRealAux
    variable = gaseous_porosity
    property = gaseous_porosity
  []
  [fuel_conductivity]
    type = MaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
  []
  [layered_average_contact_pressure]
    type = SpatialUserObjectAux
    block = pellet
    variable = layered_average_contact_pressure
    execute_on = nonlinear
    user_object = layered_average_contact_pressure
  []
  [volumetric_swelling_strain]
    type = MaterialRealAux
    variable = volumetric_swelling_strain
    property = volumetric_swelling_strain
    block = pellet
    execute_on = nonlinear
  []
[]

[Burnup]
  [burnup]
    block = pellet
    fuel_volume_ratio = ${fuel_volume_ratio}
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_power_factors
    num_radial = 80
    num_axial = 20
    isotopes = 'U235 U238 Pu239 Pu240 Pu241 Pu242'
    isotope_fractions = '${isotope_fraction_U235} ${isotope_fraction_U238} ${isotope_fraction_Pu239} ${isotope_fraction_Pu240} ${isotope_fraction_Pu241} ${isotope_fraction_Pu242}'
    RPF = RPF
    fuel_pin_geometry = pin_geometry
  []
[]

[Contact]
  [mechanical]
    model = frictionless
    formulation = mortar
    primary = 5
    secondary = 10
    c_normal = ${c_normal}
  []
[]

[ThermalContactMortar]
  [thermal]
    secondary_variable = temperature
    primary_boundary = 5
    secondary_boundary = 10
    initial_moles = initial_moles
    gas_released = fission_gas_released
    jump_distance_model = ${jump_distance_model}
    plenum_pressure = plenum_pressure
    roughness_primary = ${roughness_primary}
    roughness_secondary = ${roughness_secondary}
    roughness_coef = ${roughness_coef}
    contact_pressure = mechanical_normal_lm
    layer_thickness = layer_thickness_action
  []
[]

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

[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
  []
  [temp]
    type = FunctionDirichletBC
    boundary = '1 2 3'
    variable = temperature
    function = clad_temp_bc
  []
  [Pressure]
    [coolantPressure]
      boundary = '1 2 3'
      function = coolant_pressure
    []
  []
  [PlenumPressure]
    [plenumPressure]
      boundary = 9
      initial_pressure = ${initial_plenum_pressure}
      startup_time = ${startup_time}
      R = ${R}
      output_initial_moles = initial_moles
      temperature = plenum_temp
      volume = plenum_volume
      material_input = fission_gas_released
      output = plenum_pressure
    []
  []
[]

[Materials]
  [fuel_thermal]
    type = UO2Thermal
    block = pellet
    thermal_conductivity_model = NFIR
    temperature = temperature
    burnup_function = burnup
  []
  [fuel_volumetric_swelling]
    type = UO2VolumetricSwellingEigenstrain
    gas_swelling_model_type = SIFGRS
    block = pellet
    temperature = temperature
    burnup_function = burnup
    initial_fuel_density = ${initial_fuel_density}
    eigenstrain_name = fuel_volumetric_strain
  []
  [fuel_elasticity_tensor]
    type = UO2ElasticityTensor
    block = pellet
    density = ${initial_fuel_density}
    temperature = temperature
    matpro_youngs_modulus = true
    matpro_poissons_ratio = true
  []
  [fuel_creep]
    type = UO2CreepUpdate
    block = pellet
    temperature = temperature
    burnup_function = burnup
    initial_grain_radius = ${initial_grain_radius}
  []
  [fuel_stress]
    type = ComputeSmearedCrackingStress
    block = pellet
    cracking_stress = ${fuel_cracking_stress}
    inelastic_models = 'fuel_creep'
    softening_models = exponential_softening
    shear_retention_factor = ${fuel_shear_retention_factor}
    max_stress_correction = ${fuel_max_stress_correction}
    cracked_elasticity_type = DIAGONAL
    output_properties = crack_damage
    outputs = exodus
  []
  [exponential_softening]
    type = ExponentialSoftening
  []
  [fuel_relocation]
    type = UO2RelocationEigenstrain
    block = pellet
    burnup_function = burnup
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_power_factors
    relocation_activation1 = ${relocation_activation1}
    relocation_model = ESCORE_modified
    eigenstrain_name = fuel_relocation_strain
    fuel_pin_geometry = pin_geometry
    model_relocation_recovery = true
    max_relocation_recovery_fraction = ${max_relocation_recovery_fraction}
    relocation_scaling_factor = ${relocation_scaling_factor}
    volumetric_swelling_increment = vol_swell_increment
    layered_average_contact_pressure = layered_average_contact_pressure
    outputs = all
    output_properties = 'relocation_strain recovered_relocation_strain'
  []
  [fuel_thermal_expansion]
    type = UO2ThermalExpansionMATPROEigenstrain
    block = pellet
    temperature = temperature
    stress_free_temperature = ${stress_free_temperature}
    eigenstrain_name = fuel_thermal_strain
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = pellet
    hydrostatic_stress = hydrostatic_stress
    diff_coeff_option = TURNBULL_D1_4D2_D3
    transient_option = MICROCRACKING_BURNUP
    res_param_option = HETEROGENEOUS_WHITE
    ig_bubble_model = NUCLEATION_RESOLUTION
    ig_diff_algorithm = polypole2
    temperature = temperature
    burnup_function = burnup
    grain_radius = grain_radius
    gbs_model = true
    pellet_id = pellet_id
    pellet_brittle_zone = pbz
    ath_model = true
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_power_factors
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = pellet
    strain_free_density = ${initial_fuel_density}
  []

  [clad_thermal]
    type = ZryThermal
    block = clad
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ZryElasticityTensor
    block = clad
    fast_neutron_fluence = fast_neutron_fluence
    temperature = temperature
    matpro_youngs_modulus = true
    matpro_poissons_ratio = true
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = nonlinear
    inelastic_models = 'clad_zrycreep'
    block = clad
  []
  [clad_zrycreep]
    type = ZryCreepLimbackHoppeUpdate
    block = clad
    temperature = temperature
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    model_irradiation_creep = true
    model_primary_creep = true
    model_thermal_creep = true
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = clad
    temperature = temperature
    stress_free_temperature = ${stress_free_temperature}
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [irradiation_swelling]
    type = ZryIrradiationGrowthEigenstrain
    block = clad
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = stress_relief_annealed
    eigenstrain_name = clad_irradiation_strain
  []
  [clad_density]
    type = StrainAdjustedDensity
    block = clad
    strain_free_density = ${cladding_density}
  []
  [clad_oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    fuel_pin_geometry = pin_geometry
  []
[]

#### This is the part where VCP is set up

[Preconditioning]
  [vcp]
    type = VCP
    full = true
    # condense out the normal LM from mechanical contact and the LM from thermal contact
    lm_variable = 'mechanical_normal_lm thermal_thermal_lm'
    primary_variable = 'disp_x temperature'
    # we use LU to solve the system after condensation
    # AMG and other solver types have shown convergence for simpler mechanical/thermo-mechanical contact problems but has not luck with this one
    preconditioner = 'LU'
    is_lm_coupling_diagonal = true
    adaptive_condensation = true
  []
[]

[Dampers]
  [maxincrement]
    type = MaxIncrement
    max_increment = ${damper_max_temperature_increment}
    variable = temperature
  []
[]

[Debug]
  show_var_residual_norms = true
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'
  petsc_options = '-snes_converged_reason -ksp_converged_reason'
  # IMPORTANT: we should not specify the pc_type below, otherwise the VCP setting will be overriden
  petsc_options_iname = ' -mat_mffd_err -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = '   1e-5          NONZERO               1e-15'
  snesmf_reuse_base = false

  line_search = 'none'

  verbose = true
  l_max_its = ${l_max_its}
  l_tol = ${l_tol}

  nl_max_its = ${nl_max_its}
  nl_rel_tol = ${nl_rel_tol}
  nl_abs_tol = ${nl_abs_tol}

  start_time = ${start_time}
  n_startup_steps = ${n_startup_steps}
  end_time = ${end_time}
  dtmax = ${dtmax}
  dtmin = ${dtmin}

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = ${Timestepper_dt}
    optimal_iterations = ${Timestepper_optimal_iterations}
    iteration_window = ${Timestepper_iteration_window}
    timestep_limiting_function = power_history
    force_step_every_function_point = true
  []
[]

[Postprocessors]
  [avg_gap_conductance]
    type = SideAverageValue
    boundary = 10
    variable = gap_cond
    execute_on = 'initial timestep_end'
  []
  [max_fuel_temp]
    type = NodalExtremeValue
    block = pellet
    value_type = max
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = clad
    value_type = max
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [total_rod_integral_power]
    type = ElementIntegralPower
    variable = temperature
    burnup_function = burnup
    block = pellet
  []
  [total_rod_input_power]
    type = FunctionValuePostprocessor
    function = power_history
    scale_factor = ${rod_input_power_scale_factor}
  []
  [vol_swell_increment]
    type = SideAverageIncrementTensorComponent
    boundary = 10
    variable = volumetric_swelling_strain
    execute_on = nonlinear
  []
  [contact_dof]
    type = ContactDOFSetSize
    variable = 'mechanical_normal_lm'
    subdomain = 'secondary_lower'
    execute_on = 'nonlinear timestep_end'
  []
  [nl_its]
    type = NumNonlinearIterations
  []
  [total_nl_its]
    type = CumulativeValuePostprocessor
    postprocessor = nl_its
  []
  [l_its]
    type = NumLinearIterations
  []
  [total_l_its]
    type = CumulativeValuePostprocessor
    postprocessor = l_its
  []
[]

[StandardLWRFuelRodOutputs]
  fuel_pellet_blocks = 3
[]

[Outputs]
  perf_graph = true
  csv = true
  exodus = true
  color = false
  [console]
    type = Console
    max_rows = 25
  []
  [chkfile]
    type = CSV
    show = 'average_burnup fission_gas_released_percentage'
    execute_on = 'FINAL'
  []
[]

[VectorPostprocessors]
  [temperature_post]
    type = NodalValueSampler
    variable = temperature
    boundary = '10'
    sort_by = y
  []
  [contact_post]
    type = NodalValueSampler
    variable = mechanical_normal_lm
    boundary = '10'
    sort_by = y
  []
  [disp_x]
    type = NodalValueSampler
    variable = disp_x
    boundary = '10'
    sort_by = y
  []
  [disp_y]
    type = NodalValueSampler
    variable = disp_y
    boundary = '10'
    sort_by = y
  []
[]

[UserObjects]
  [pin_geometry]
    type = FuelPinGeometry
  []
  [pbz]
    type = PelletBrittleZone
    block = pellet
    pellet_id = pellet_id
    temperature = temperature
    fuel_pin_geometry = pin_geometry
    number_pellets = ${pellet_quantity}
    execute_on = 'initial linear'
  []
  [layered_average_contact_pressure]
    type = LayeredSideAverage
    variable = mechanical_normal_lm
    direction = y
    num_layers = 1
    execute_on = timestep_end
    boundary = 10
  []
[]

Description
Example Input Syntax
Input Parameters
Input Files
References