MNCreepUpdate

Calculates thermal and irradiation creep rate for mixed mononitrides (UN and UPuN) fuel.

Description

MNCreepUpdate calculates the thermal and irradiation creep of MN due to applied stresses: $var element = document.getElementById("moose-equation-6efab1c5-3b21-44bf-936a-642ec9435755");katex.render(" \\dot{\\epsilon}_{total} = \\dot{\\epsilon}_{th} + \\dot{\\epsilon}_{irr}.", element, {displayMode:true,throwOnError:false});$

The thermal creep rate at theoretical density, $var element = document.getElementById("moose-equation-e148eae6-1005-446c-9281-f59e9998000b");katex.render("\\dot{\\epsilon}_{TD}", element, {displayMode:false,throwOnError:false});$ , is from Hayes et al. (1990): $var element = document.getElementById("moose-equation-4e6ab85a-4ae6-4f94-99fa-4bd89e8b477d");katex.render("\\dot{\\epsilon}_{TD} = 2.054\\times 10^{-3} \\sigma^{4.5} \\exp{\\left(-\\frac{39369.5}{T}\\right)},", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-b905217c-1f93-4fbd-8fc9-0ec8b84b81be");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the temperature (K) and $var element = document.getElementById("moose-equation-5ca2cf67-6489-4b0e-8862-07aa081bff2e");katex.render("\\sigma", element, {displayMode:false,throwOnError:false});$ is the effective stress (MPa). A porosity correction can be applied according to, $var element = document.getElementById("moose-equation-b5f7326a-17d1-4ce3-9f14-cd7bf2547451");katex.render(" f(\\phi) = \\frac{0.987}{(1-\\phi)^{27.6}} \\exp(-8.65\\phi),", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-6ca6926d-972c-405a-a135-6dcd16507ce6");katex.render("\\phi", element, {displayMode:false,throwOnError:false});$ is the porosity. The correction is multiplied, i.e., $var element = document.getElementById("moose-equation-d835c0c8-8046-45c6-9ef1-44ab9653c5b5");katex.render(" \\dot{\\epsilon}_{th} = \\dot{\\epsilon}_{TD} \\times f(\\phi).", element, {displayMode:true,throwOnError:false});$

The creep due to irradiation is from Brucklacher and Dienst (1972), as compiled in Blank (2006), based on a typical porosity of 4%, $var element = document.getElementById("moose-equation-df27c349-9fd4-4c4d-96fc-9ad97fba3acf");katex.render(" \\dot{\\epsilon}_{irr} = 2\\times 10^{-22} \\sigma F,", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-aa3fee18-5fff-4514-9ef7-02f5bf7ee293");katex.render("\\sigma", element, {displayMode:false,throwOnError:false});$ is the effective stress in MPa, $var element = document.getElementById("moose-equation-d4290f7c-6ceb-4784-961e-431d074cf7a6");katex.render("F", element, {displayMode:false,throwOnError:false});$ is the fission rate (fissions/cm $var element = document.getElementById("moose-equation-ad0f0edc-e93f-401a-af9b-314e5f878ba1");katex.render("^3", element, {displayMode:false,throwOnError:false});$ /s), and the creep $var element = document.getElementById("moose-equation-f07ab80b-5548-4220-abb9-8849b53b3083");katex.render("\\dot{\\epsilon}_{irr}", element, {displayMode:false,throwOnError:false});$ is given in 1/h. Note that the units will be converted in the code to use the fission rate in fissions/m $var element = document.getElementById("moose-equation-1b1ac0eb-b880-4c5a-8774-601d3827a35b");katex.render("^3", element, {displayMode:false,throwOnError:false});$ /s, the effective stress in Pa, and return the creep rate in 1/s.

Example Input Syntax

note:Need for inelastic strain return mapping stress

MNCreepUpdate must be run in conjunction with the inelastic strain return mapping stress calculator, as shown below.

[Materials<<<{"href": "../../../syntax/Materials/index.html"}>>>]
  [creep]
    type = MNCreepUpdate<<<{"description": "Calculates thermal and irradiation creep rate for mixed mononitrides (UN and UPuN) fuel.", "href": "MNCreepUpdate.html"}>>>
    temperature<<<{"description": "The coupled temperature (K)"}>>> = temp
    porosity<<<{"description": "Porosity material property name"}>>> = 0.1
    fission_rate<<<{"description": "Fission rate material property name"}>>> = 1.042e20
  []

  [radial_return_stress]
    type = ComputeMultipleInelasticStress<<<{"description": "Compute state (stress and internal parameters such as plastic strains and internal parameters) using an iterative process.  Combinations of creep models and plastic models may be used.", "href": "../ComputeMultipleInelasticStress.html"}>>>
    inelastic_models<<<{"description": "The material objects to use to calculate stress and inelastic strains. Note: specify creep models first and plasticity models second."}>>> = 'creep'
  []
[]

Input Parameters

fission_rateFission rate material property name
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Fission rate material property name
porosityPorosity material property name
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Porosity material property name
temperatureThe coupled temperature (K)
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The coupled temperature (K)

Required Parameters

absolute_tolerance1e-11Absolute convergence tolerance for Newton iteration
Default:1e-11
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Absolute convergence tolerance for Newton iteration
acceptable_multiplier10Factor applied to relative and absolute tolerance for acceptable convergence if iterations are no longer making progress
Default:10
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Factor applied to relative and absolute tolerance for acceptable convergence if iterations are no longer making progress
adaptive_substeppingFalseUse adaptive substepping, where the number of substeps is successively doubled until the return mapping model successfully converges or the maximum number of substeps is reached.
Default:False
C++ Type:bool
Controllable:No
Description:Use adaptive substepping, where the number of substeps is successively doubled until the return mapping model successfully converges or the maximum number of substeps is reached.
automatic_differentiation_return_mappingFalseWhether to use automatic differentiation to compute the derivative.
Default:False
C++ Type:bool
Controllable:No
Description:Whether to use automatic differentiation to compute the derivative.
base_nameOptional parameter that defines a prefix for all material properties related to this stress update model. This allows for multiple models of the same type to be used without naming conflicts.
C++ Type:std::string
Controllable:No
Description:Optional parameter that defines a prefix for all material properties related to this stress update model. This allows for multiple models of the same type to be used without naming conflicts.
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
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.
execute_onLINEARThe list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html.
Default:LINEAR
C++ Type:ExecFlagEnum
Options:XFEM_MARK, NONE, INITIAL, LINEAR, NONLINEAR_CONVERGENCE, NONLINEAR, POSTCHECK, TIMESTEP_END, TIMESTEP_BEGIN, MULTIAPP_FIXED_POINT_END, MULTIAPP_FIXED_POINT_BEGIN, FINAL, CUSTOM
Controllable:No
Description:The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html.
max_inelastic_increment0.0001The maximum inelastic strain increment allowed in a time step
Default:0.0001
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The maximum inelastic strain increment allowed in a time step
maximum_number_substeps25The maximum number of substeps allowed before cutting the time step.
Default:25
C++ Type:unsigned int
Controllable:No
Description:The maximum number of substeps allowed before cutting the time step.
relative_tolerance1e-08Relative convergence tolerance for Newton iteration
Default:1e-08
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Relative convergence tolerance for Newton iteration
use_substep_integration_errorFalseIf true, it establishes a substep size that will yield, at most,the creep numerical integration error given by substep_strain_tolerance.
Default:False
C++ Type:bool
Controllable:No
Description:If true, it establishes a substep size that will yield, at most,the creep numerical integration error given by substep_strain_tolerance.
use_substeppingNONEWhether and how to use substepping
Default:NONE
C++ Type:MooseEnum
Options:NONE, ERROR_BASED, INCREMENT_BASED
Controllable:No
Description:Whether and how to use substepping
value_range_behaviorEXCEPTIONWhat to do if input value is outside the range of applicability.
Default:EXCEPTION
C++ Type:MooseEnum
Options:ERROR, WARN, IGNORE, EXCEPTION
Controllable:No
Description:What to do if input value is outside the range of applicability.

Optional Parameters

apply_strainTrueFlag to apply strain. Used for testing.
Default:True
C++ Type:bool
Controllable:No
Description:Flag to apply strain. Used for testing.
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.
effective_inelastic_strain_nameeffective_creep_strainName of the material property that stores the effective inelastic strain
Default:effective_creep_strain
C++ Type:std::string
Controllable:No
Description:Name of the material property that stores the effective inelastic strain
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
substep_strain_tolerance0.1Maximum ratio of the initial elastic strain increment at start of the return mapping solve to the maximum inelastic strain allowable in a single substep. Reduce this value to increase the number of substeps
Default:0.1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Maximum ratio of the initial elastic strain increment at start of the return mapping solve to the maximum inelastic strain allowable in a single substep. Reduce this value to increase the number of substeps
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

internal_solve_full_iteration_historyFalseSet true to output full internal Newton iteration history at times determined by `internal_solve_output_on`. If false, only a summary is output.
Default:False
C++ Type:bool
Controllable:No
Description:Set true to output full internal Newton iteration history at times determined by `internal_solve_output_on`. If false, only a summary is output.
internal_solve_output_onon_errorWhen to output internal Newton solve information
Default:on_error
C++ Type:MooseEnum
Options:never, on_error, always
Controllable:No
Description:When to output internal Newton solve information

Debug 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

Input Files

(assessment/nitride/PW-CANEL/SNAP50/analysis/PW_SNAP50_UN_Pin_base.i)
(test/tests/solid_mechanics/mn_creep/exact.i)
(test/tests/solid_mechanics/mn_creep/ad_cycle.i)
(assessment/nitride/EBRII/SP1/analysis/SP1_base.i)
(assessment/nitride/JOYO/JOYO_L4C4/analysis/JOYO_Pin_base.i)
(test/tests/solid_mechanics/mn_creep/ad.i)
(assessment/nitride/MTR/SNAP50/analysis/SNAP50_base.i)
(assessment/nitride/EBRII/SP1/analysis/SP1_Fuel_focused_base.i)
(test/tests/solid_mechanics/mn_creep/nonad.i)
(assessment/nitride/EBRII/K4/analysis/base.i)

References

Hubert Blank. Nonoxide Ceramic Nuclear Fuels, chapter, pages. John Wiley & Sons, Ltd, 2006. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/9783527603978.mst0108, doi:10.1002/9783527603978.mst0108.

@inbook{blank:2006,
    author = "Blank, Hubert",
    publisher = "John Wiley \\& Sons, Ltd",
    isbn = "9783527603978",
    title = "Nonoxide Ceramic Nuclear Fuels",
    booktitle = "Materials Science and Technology",
    chapter = "",
    pages = "",
    doi = "10.1002/9783527603978.mst0108",
    url = "https://onlinelibrary.wiley.com/doi/abs/10.1002/9783527603978.mst0108",
    year = "2006"
}

D. Brucklacher and W. Dienst. Creep behavior of ceramic nuclear fuels under neutron irradiation. Journal of Nuclear Materials, 42(3):285–296, 1972. doi:10.1016/0022-3115(72)90079-7.

@article{brucklacher_et_al_1972,
    author = "Brucklacher, D. and Dienst, W.",
    title = "Creep behavior of ceramic nuclear fuels under neutron irradiation",
    journal = "Journal of Nuclear Materials",
    year = "1972",
    volume = "42",
    number = "3",
    pages = "285-296",
    doi = "10.1016/0022-3115(72)90079-7"
}

S.L. Hayes, J.K. Thomas, and K.L. Peddicord. Material property correlations for uranium mononitride: ii. mechanical properties. Journal of Nuclear Materials, 171(2):271 – 288, 1990. URL: http://www.sciencedirect.com/science/article/pii/002231159090375W, doi:https://doi.org/10.1016/0022-3115(90)90375-W.

@ARTICLE{hayes1990b,
    author = "Hayes, S.L. and Thomas, J.K. and Peddicord, K.L.",
    title = "Material property correlations for uranium mononitride: II. Mechanical properties",
    journal = "Journal of Nuclear Materials",
    volume = "171",
    number = "2",
    pages = "271 - 288",
    year = "1990",
    issn = "0022-3115",
    doi = "https://doi.org/10.1016/0022-3115(90)90375-W",
    url = "http://www.sciencedirect.com/science/article/pii/002231159090375W"
}

(test/tests/solid_mechanics/mn_creep/nonad.i)

# Tests the MNCreepUpdate models with other non-AD thermo-mechanical models

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

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

[Variables]
  [temp]
    initial_condition = 900
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [pellets]
    add_variables = true
    strain = FINITE
    incremental = true
    generate_output = 'stress_yy strain_xx strain_yy'
  []
[]

[Kernels]
  [heat_ie]
    type = TimeDerivative
    variable = temp
  []
[]

[BCs]
  [Pressure]
    [u_top_pull]
      boundary = top
      factor = -2e7
    []
  []
  [u_bottom_fix]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 7.5e10
    poissons_ratio = 0.3
  []
  [creep]
    type = MNCreepUpdate
    temperature = temp
    porosity = 0.1
    fission_rate = 1.042e20
  []
  [radial_return_stress]
    type = ComputeMultipleInelasticStress
    inelastic_models = 'creep'
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  num_steps = 10
  dt = 10000
[]

[Postprocessors]
   [average_stress_yy]
     type = ElementAverageValue
     variable = stress_yy
   []
   [average_strain_xx]
     type = ElementAverageValue
     variable = strain_xx
   []
   [average_strain_yy]
     type = ElementAverageValue
     variable = strain_yy
   []
   [disp_x]
     type = NodalVariableValue
     nodeid = 7
     variable = disp_x
  []
  [disp_y]
     type = NodalVariableValue
     nodeid = 7
     variable = disp_y
  []
[]

[Outputs]
  exodus = true
[]

(assessment/nitride/PW-CANEL/SNAP50/analysis/PW_SNAP50_UN_Pin_base.i)

[GlobalParams]
  order = FIRST
  energy_per_fission = ${energy_per_fission}
  displacements = 'disp_x disp_y'
  value_range_behavior = IGNORE
[]

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

[Mesh]
  coord_type = RZ
  [smeared_pellet_mesh]
    type = FuelPinMeshGenerator
    clad_thickness = ${cladding_liner_thickness}
    pellet_outer_radius = ${fuel_radius}
    pellet_height = ${fuel_height}
    clad_top_gap_height = ${plenum_height}
    clad_gap_width = ${cladding_gap_width}
    bottom_clad_height = ${cladding_bottom_top_plug_length}
    top_clad_height = ${cladding_bottom_top_plug_length}
    clad_bot_gap_height = ${gap_bottom_length}
    liner_thickness = ${liner_thickness}
    clad_mesh_density = customize
    pellet_mesh_density = customize
    nx_p = 6
    ny_p = 260
    nx_c = 4
    ny_c = 260
    ny_cu = 3
    ny_cl = 3
    nx_liner = 2
    pellet_quantity = 1
    elem_type = QUAD4
  []
  patch_size = 30
  patch_update_strategy = auto
  partitioner = centroid
  centroid_partitioner_direction = y
[]

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

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 ${time_end_ramp_up} ${time_start_ramp_down} ${time_end_ramp_down}'
    y = '0 ${avg_lin_power} ${avg_lin_power} 0'
  []
  [coolant_wall_temperature]
    type = ParsedFunction
    expression = 'if(t<${time_end_ramp_up}, ${initial_temperature} + t*(${a}*(y-${x_0})^2 + ${b} - ${initial_temperature})/(${time_end_ramp_up}),if(t<${time_start_ramp_down},${a}*(y-${x_0})^2 + ${b}, if(t<${time_end_ramp_down}, ${a}*(y-${x_0})^2 + ${b} + (t-${time_start_ramp_down})*(${initial_temperature}-(${a}*(y-${x_0})^2 + ${b}))/(${time_end_ramp_down}-${time_start_ramp_down}),${initial_temperature})))     '
  []
  [coolant_press_ramp]
    type = PiecewiseLinear
    x = '0       ${run_time}'
    y = '${coolant_pressure} ${coolant_pressure}' # unknown, Li coolant
  []
  [plenum_pressure]
    type = ConstantFunction
    value = ${initial_plenum_pressure}
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [fuel]
    block = pellet
    strain = FINITE
    add_variables = true
    generate_output = 'firstinv_strain stress_xx stress_yy stress_zz vonmises_stress hydrostatic_stress creep_strain_xx creep_strain_yy creep_strain_zz elastic_strain_xx elastic_strain_yy elastic_strain_zz strain_xx strain_yy strain_zz'
    extra_vector_tags = 'ref'
    eigenstrain_names = 'fuel_thermal_expansion solid_swelling_eigenstrain'
    temperature = temperature
  []
  [clad]
    block = clad
    strain = FINITE
    incremental = true
    add_variables = true
    generate_output = 'firstinv_strain stress_xx stress_yy stress_zz vonmises_stress hydrostatic_stress elastic_strain_xx elastic_strain_yy elastic_strain_zz strain_xx strain_yy strain_zz'
    extra_vector_tags = 'ref'
    eigenstrain_names = 'clad_thermal_expansion'
    temperature = temperature
  []
  [liner]
    block = liner
    strain = FINITE
    incremental = true
    add_variables = true
    generate_output = 'firstinv_strain stress_xx stress_yy stress_zz vonmises_stress hydrostatic_stress elastic_strain_xx elastic_strain_yy elastic_strain_zz strain_xx strain_yy strain_zz'
    extra_vector_tags = 'ref'
    eigenstrain_names = 'liner_thermal_expansion'
    temperature = temperature
  []
[]

[Kernels]
  [gravity]
    type = Gravity
    variable = disp_y
    value = -9.81
    extra_vector_tags = 'ref'
    block = 'pellet liner clad'
  []
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
    block = 'pellet liner clad'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
    block = 'pellet liner clad'
  []
  [heat_source]
    type = FissionRateHeatSource
    variable = temperature
    fission_rate = fission_rate
    extra_vector_tags = 'ref'
    energy_deposited_in_fuel = 0.95
    block = 'pellet'
  []
[]

[ThermalContactMortar]
  [thermal_contact]
    secondary_variable = temperature
    primary_boundary = '5' # clad_inside_right
    secondary_boundary = '10' # pellet_outer_radial_surface
    initial_moles = initial_moles # coupling to a postprocessor which supplies the initial plenum/gap gas mass
  []
[]

[Contact]
  [mechanical]
    model = frictionless
    formulation = mortar
    primary = 5 # clad_inside_right
    secondary = 10 # pellet_outer_radial_surface
    c_normal = 1e+11
    correct_edge_dropping = true
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = '12' # centerline
    value = 0.0
  []
  [no_y_clad]
    type = DirichletBC
    variable = disp_y
    boundary = '1' # clad_outside_bottom
    value = 0.0
  []
  [no_y_fuel]
    type = DirichletBC
    variable = disp_y
    boundary = '20' # bottom_of_bottom_pellet
    value = 0.0
  []
  [Pressure]
    [coolantPressure]
      boundary = '1 2 3' # clad_outside_bottom, clad_outside_right, clad_outside_top
      function = coolant_press_ramp
    []
  []
  [PlenumPressure] # apply plenum pressure on clad inner walls and pellet surfaces
    [plenumPressure]
      boundary = 9 # inside_surfaces
      initial_pressure = ${initial_plenum_pressure}
      startup_time = 0
      R = ${R}
      initial_temperature = ${initial_temperature}
      temperature = ave_temperature_interior # coupling to post processor to get gas temperature approximation
      output_initial_moles = initial_moles # coupling to post processor to get initial fill gas mass
      volume = gas_volume # coupling to post processor to get gas volume
      output = plenum_pressure # coupling to post processor to output plenum/gap pressure
    []
  []
  [clad_outer_temperature]
    type = FunctionDirichletBC
    boundary = '1 2 3' # clad_outside_bottom, clad_outside_right, clad_outside_top
    function = coolant_wall_temperature
    variable = temperature
  []
[]

[Materials]
  [porosity]
    type = GenericConstantMaterial
    block = pellet
    prop_names = porosity
    prop_values = ${initial_porosity}
    outputs = all
  []
  [fission_rate]
    type = FissionRate
    block = pellet
    rod_linear_power = power_history
    axial_power_profile = 1
    pellet_radius = ${fuel_radius}
    outputs = all
  []
  [fuel_thermal]
    type = MNThermal
    block = pellet
    temperature = temperature
    porosity = porosity
    outputs = all
  []
  [fuel_elasticity_tensor]
    block = pellet
    type = MNElasticityTensor
    temperature = temperature
    porosity = porosity
    output_properties = 'youngs_modulus poissons_ratio'
    outputs = all
  []
  [fuel_thermal_expansion]
    block = pellet
    type = MNThermalExpansionEigenstrain
    eigenstrain_name = fuel_thermal_expansion
    stress_free_temperature = ${stress_free_temperature_MN}
    temperature = temperature
  []
  [fuel_creep]
    block = pellet
    type = MNCreepUpdate
    max_inelastic_increment = 1e-4
    temperature = temperature
    porosity = porosity
    fission_rate = fission_rate
    outputs = all
  []
  [burnup]
    type = Burnup
    block = pellet
    atoms_heavy_metal_per_volume = ${atoms_heavy_metal_per_volume}
    outputs = all
  []
  [burnup_swelling]
    type = BurnupDependentEigenstrain
    block = pellet
    eigenstrain_name = 'solid_swelling_eigenstrain'
    burnup = burnup
    outputs = all
  []
  [fuel_radial_return_stress]
    block = pellet
    type = ComputeMultipleInelasticStress
    inelastic_models = 'fuel_creep'
  []
  [fuel_density]
    block = pellet
    type = StrainAdjustedDensity
    strain_free_density = ${density}
  []

  [clad_elasticity_tensor]
    block = clad
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 68.9e9
    poissons_ratio = 0.4
  []
  [clad_thermal_expansion]
    block = clad
    type = ComputeThermalExpansionEigenstrain
    eigenstrain_name = clad_thermal_expansion
    thermal_expansion_coeff = 7.54e-6
    stress_free_temperature = ${stress_free_temperature_NbZr}
    temperature = temperature
  []
  [clad_stress]
    block = clad
    type = ComputeFiniteStrainElasticStress
  []
  [clad_thermal]
    type = HeatConductionMaterial
    block = clad
    thermal_conductivity = 41.9
    specific_heat = 270
  []
  [clad_density]
    block = clad
    type = StrainAdjustedDensity
    strain_free_density = 8590
  []

  [liner_elasticity_tensor]
    block = liner
    type = TungstenElasticityTensor
    temperature = temperature
  []
  [liner_thermal_expansion]
    block = liner
    type = TungstenThermalExpansionEigenstrain
    eigenstrain_name = liner_thermal_expansion
    stress_free_temperature = ${stress_free_temperature_W}
    temperature = temperature
  []
  [liner_stress]
    block = liner
    type = ComputeFiniteStrainElasticStress
  []
  [liner_thermal]
    block = liner
    type = TungstenThermal
    temperature = temperature
  []
  [liner_density]
    block = liner
    type = StrainAdjustedDensity
    strain_free_density = 19300
  []
[]

[Dampers]
  [disp_x]
    type = MaxIncrement
    variable = disp_x
    max_increment = 1e-4
    min_damping = 1e-4
  []
  [disp_y]
    type = MaxIncrement
    variable = disp_y
    max_increment = 1e-3
    min_damping = 1e-4
  []
  [temperature]
    type = MaxIncrement
    variable = temperature
    max_increment = 50
    min_damping = 1e-4
  []
[]

[Executioner]
  type = Transient
  # With mortar contact
  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -mat_mffd_err -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu    superlu_dist   1e-5          NONZERO               1e-15'
  snesmf_reuse_base = false
  line_search = 'none'

  l_max_its = 60
  nl_max_its = 20
  nl_rel_tol = 5e-6
  nl_abs_tol = 5e-9

  end_time = ${run_time}
  dtmin = 1
  dtmax = 5e4

  [TimeStepper]
    type = IterationAdaptiveDT
    timestep_limiting_postprocessor = time_step_limit
    force_step_every_function_point = true
    timestep_limiting_function = power_history
    dt = 10
  []
[]

[Postprocessors]
  # elemental temperatures
  [temperature_fuel_avg]
    type = ElementAverageValue
    variable = temperature
    block = pellet
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_max]
    type = ElementExtremeValue
    variable = temperature
    block = pellet
  []
  [temperature_fuel_min]
    type = ElementExtremeValue
    variable = temperature
    block = pellet
    value_type = min
  []

  [temperature_clad_avg]
    type = ElementAverageValue
    variable = temperature
    block = clad
    execute_on = 'initial timestep_end'
  []
  [temperature_clad_max]
    type = ElementExtremeValue
    variable = temperature
    block = clad
  []
  [temperature_clad_min]
    type = ElementExtremeValue
    variable = temperature
    block = clad
    value_type = min
  []

  [temperature_liner_avg]
    type = ElementAverageValue
    variable = temperature
    block = liner
    execute_on = 'initial timestep_end'
  []
  [temperature_liner_max]
    type = ElementExtremeValue
    variable = temperature
    block = liner
  []
  [temperature_liner_min]
    type = ElementExtremeValue
    variable = temperature
    block = liner
    value_type = min
  []

  [ave_temperature_interior] # average temperature of the cladding interior and all pellet exteriors
    type = SideAverageValue
    boundary = '9' # inside_surfaces
    variable = temperature
    execute_on = 'initial linear'
  []
  [temperature_fuel_centerline_avg]
    type = AxisymmetricCenterlineAverageValue
    boundary = '12' # centerline
    variable = temperature
  []
  [temperature_fuel_centerline_max]
    type = NodalExtremeValue
    boundary = '12' # centerline
    variable = temperature
  []
  [temperature_fuel_centerline_min]
    type = NodalExtremeValue
    boundary = '12' # centerline
    variable = temperature
    value_type = min
  []
  [temperature_fuel_surface_avg]
    type = SideAverageValue
    boundary = '10' # pellet_outer_radial_surface
    variable = temperature
  []
  [temperature_fuel_surface_max]
    type = NodalExtremeValue
    boundary = '10' # pellet_outer_radial_surface
    variable = temperature
  []
  [temperature_fuel_surface_min]
    type = NodalExtremeValue
    boundary = '10' # pellet_outer_radial_surface
    variable = temperature
    value_type = min
  []

  # stresses
  [stress_vonmises_fuel_avg]
    type = ElementAverageValue
    variable = vonmises_stress
    block = pellet
  []
  [stress_vonmises_fuel_max]
    type = ElementExtremeValue
    variable = vonmises_stress
    block = pellet
  []
  [stress_vonmises_fuel_min]
    type = ElementExtremeValue
    variable = vonmises_stress
    value_type = min
    block = pellet
  []
  [stress_hydro_fuel_avg]
    type = ElementAverageValue
    variable = hydrostatic_stress
    block = pellet
  []
  [stress_hydro_fuel_max]
    type = ElementExtremeValue
    variable = hydrostatic_stress
    block = pellet
  []
  [stress_hydro_fuel_min]
    type = ElementExtremeValue
    variable = hydrostatic_stress
    value_type = min
    block = pellet
  []

  # strain information
  [strain_solid_swelling_fuel_avg]
    type = ElementAverageValue
    variable = solid_swelling
    block = pellet
  []
  [strain_volumetric_fuel_avg]
    type = ElementAverageValue
    variable = firstinv_strain
    block = pellet
  []
  [strain_axial_fuel_avg]
    type = ParsedPostprocessor
    pp_names = 'disp_y_fuel_top_surface_avg disp_y_fuel_bottom_surface_avg'
    expression = '(disp_y_fuel_top_surface_avg - disp_y_fuel_bottom_surface_avg) / ${fuel_height}'
  []
  [disp_y_fuel_top_surface_avg]
    type = SideAverageValue
    variable = disp_y
    boundary = '11' # top_of_top_pellet
  []
  [disp_y_fuel_top_surface_max]
    type = NodalExtremeValue
    variable = disp_y
    boundary = '11' # top_of_top_pellet
  []
  [disp_y_fuel_bottom_surface_avg]
    type = SideAverageValue
    variable = disp_y
    boundary = '20' # bottom_of_bottom_pellet
  []
  [disp_y_fuel_bottom_surface_max]
    type = NodalExtremeValue
    variable = disp_y
    boundary = '20' # bottom_of_bottom_pellet
  []
  [disp_x_fuel_radial_surface_max]
    type = NodalExtremeValue
    variable = disp_x
    boundary = '10' # pellet_outer_radial_surface
  []
  [disp_x_fuel_radial_surface_avg]
    type = SideAverageValue
    variable = disp_x
    boundary = '10' # pellet_outer_radial_surface
  []

  [burnup_avg]
    type = ElementAverageValue
    block = pellet
    variable = burnup
  []
  [burnup_max]
    type = ElementExtremeValue
    block = pellet
    variable = burnup
  []
  [fission_rate_avg]
    type = ElementAverageValue
    variable = fission_rate
    block = pellet
  []
  [porosity_fuel_avg]
    type = ElementAverageValue
    variable = porosity
    block = pellet
  []
  [porosity_fuel_max]
    type = ElementExtremeValue
    variable = porosity
    block = pellet
  []
  [porosity_fuel_min]
    type = ElementExtremeValue
    variable = porosity
    value_type = min
    block = pellet
  []

  [swelling_dia_percent]
    type = ParsedPostprocessor
    pp_names = 'disp_x_fuel_radial_surface_avg'
    expression = 'disp_x_fuel_radial_surface_avg / 3.14159 / ${fuel_radius} * 100'
  []
  [fuel_volume]
    type = VolumePostprocessor
    block = pellet
    execute_on = 'TIMESTEP_END INITIAL'
    use_displaced_mesh = true
  []
  [swelling_vol_percent]
    type = ParsedPostprocessor
    pp_names = 'fuel_volume'
    expression = '(fuel_volume - ${fuel_volume}) / ${fuel_volume} * 100'
  []
  [gas_volume]
    type = InternalVolume
    boundary = '9' # inside_surfaces
    execute_on = 'initial linear'
  []

  # extras
  [actual_time_step_limit]
    type = MaterialTimeStepPostprocessor
    outputs = none
    block = 'pellet'
  []
  [time_step_limit]
    type = ParsedPostprocessor
    expression = 'if(actual_time_step_limit > 1e6, 1e6, actual_time_step_limit)'
    pp_names = 'actual_time_step_limit'
  []
[]

[PerformanceMetricOutputs]
  outputs = 'performance_metrics performance_metrics_sync exodus console base_out'
[]

[Outputs]
  perf_graph = true
  csv = true
  sync_times = '1e2 1e3 5e3 1e4 5e4 1e5 5e5 1e6 2.5e6 5e6 1e7 2e7 3e7 ${fparse run_time -1e3} ${run_time}'
  file_base = '${group_name}_nominal_Pin'
  [base_out]
    type = CSV
    file_base = '${group_name}_Pin_base_out'
    show = 'burnup_max temperature_fuel_avg temperature_clad_avg temperature_liner_avg ave_temperature_interior swelling_vol_percent disp_x_fuel_radial_surface_max'
    sync_only = true
  []
  [checkpoint]
    type = Checkpoint
    time_step_interval =  10
  []
  [exodus]
    type = Exodus
  []
  [sync]
    type = CSV
    file_base = '${group_name}_nominal_Pin_sync'
    sync_only = true
  []
  [console]
    type = Console
    show = 'time_step_size temperature_fuel_avg temperature_fuel_centerline_max   stress_vonmises_fuel_max stress_hydro_fuel_max stress_hydro_fuel_min strain_axial_fuel_avg burnup_avg fission_rate_avg porosity_fuel_avg time_step_limit disp_x_fuel_radial_surface_max disp_x_fuel_radial_surface_avg swelling_dia_percent swelling_vol_percent fuel_volume'
  []
  [performance_metrics]
    type = CSV
    file_base = '${group_name}_nominal_performance_metrics_Pin'
    show = 'simulation_alive_time number_linear_iterations number_nonlinear_iterations time_step_size total_linear_iterations total_nonlinear_iterations physical_memory_use number_dofs number_nonlinear_variables residual_compute_time jacobian_compute_time'
  []
  [performance_metrics_sync]
    type = CSV
    sync_only = true
    file_base = '${group_name}_nominal_performance_metrics_Pin_sync'
    show = 'simulation_alive_time number_linear_iterations number_nonlinear_iterations time_step_size total_linear_iterations total_nonlinear_iterations physical_memory_use number_dofs number_nonlinear_variables residual_compute_time jacobian_compute_time'
  []
[]

[Debug]
  show_var_residual = 'disp_x disp_y temperature'
  show_var_residual_norms = true
[]

(test/tests/solid_mechanics/mn_creep/exact.i)

# Tests the MNCreepUpdate model against hand calculations

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  coord_type = RZ
  use_displaced_mesh = false
  [gen]
    type = ExamplePatchMeshGenerator
    dim = 2
  []
[]

[AuxVariables]
  [temp]
    initial_condition = 600
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [pellets]
    add_variables = true
    strain = SMALL
    incremental = true
    generate_output = 'stress_yy strain_xx strain_yy vonmises_stress'
  []
[]

[AuxKernels]
  [temp_aux]
    type = FunctionAux
    variable = temp
    function = '800 + .25 * t'
  []
[]

[BCs]
  [Pressure]
    [u_top_pull]
      boundary = top
      factor = -2e8
      function = 't / 1e4'
    []
  []
  [u_bottom_fix]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  []
[]

[Functions]
  [porosity_ramp]
    type = PiecewiseLinear
    x = '0 1000'
    y = '0 0.2'
  []
  [fission_rate_ramp]
    type = PiecewiseLinear
    x = '0    900'
    y = '1e23 0.1'
  []
[]

[Materials]
  [porosity]
    type = GenericFunctionMaterial
    prop_names = 'porosity fission_rate'
    prop_values = 'porosity_ramp fission_rate_ramp'
    outputs = all
  []
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 7.5e10
    poissons_ratio = 0.3
  []
  [creep]
    type = MNCreepUpdate
    temperature = temp
    porosity = porosity
    fission_rate = fission_rate
    outputs = all
  []
  [radial_return_stress]
    type = ComputeMultipleInelasticStress
    inelastic_models = 'creep'
  []
[]

[Executioner]
  type = Transient
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  num_steps = 10
  dt = 100
  dtmin = 100
[]

[Postprocessors]
  [vonmises_avg]
    type = ElementAverageValue
    variable = vonmises_stress
  []
  [strain_xx_avg]
    type = ElementAverageValue
    variable = strain_xx
  []
  [strain_yy_avg]
    type = ElementAverageValue
    variable = strain_yy
  []
  [porosity_avg]
    type = ElementAverageValue
    variable = porosity
  []
  [fission_rate_avg]
    type = ElementAverageValue
    variable = fission_rate
  []
  [creep_rate_avg]
    type = ElementAverageValue
    variable = creep_rate
  []
  [temp_avg]
    type = ElementAverageValue
    variable = temp
  []
  [dt]
    type = TimestepSize
    outputs = none
  []

  [thermal_strain_avg_old]
    type = ElementAverageValue
    variable = thermal_creep_strain
    execute_on = TIMESTEP_BEGIN
    outputs = none
  []
  [thermal_strain_avg]
    type = ElementAverageValue
    variable = thermal_creep_strain
  []
  [thermal_strain_inc]
    type = ParsedPostprocessor
    pp_names = 'temp_avg vonmises_avg porosity_avg dt'
    expression = 'thermal_creep_rate := 2.054e-3 * (vonmises_avg / 1.0e6)^4.5 * exp(-39369.5 / temp_avg);
                porosity_correction := 0.987 / (1 - porosity_avg)^27.6 * exp(-8.65 * porosity_avg);
                thermal_creep_rate * porosity_correction * dt'
    outputs = none
  []
  [thermal_strain_exact]
    type = ParsedPostprocessor
    pp_names = 'thermal_strain_inc thermal_strain_avg_old'
    expression = 'thermal_strain_inc + thermal_strain_avg_old'
    outputs = none
  []
  [thermal_strain_diff]
    type = ParsedPostprocessor
    pp_names = 'thermal_strain_exact thermal_strain_avg'
    expression = '(thermal_strain_exact - thermal_strain_avg) / thermal_strain_exact'
    outputs = none
  []
  [thermal_strain_diff_max]
    type = TimeExtremeValue
    postprocessor = 'thermal_strain_diff'
    value_type = 'abs_max'
  []

  [irradiation_strain_avg_old]
    type = ElementAverageValue
    variable = irradiation_creep_strain
    execute_on = TIMESTEP_BEGIN
    outputs = none
  []
  [irradiation_strain_avg]
    type = ElementAverageValue
    variable = irradiation_creep_strain
  []
  [irradiation_strain_inc]
    type = ParsedPostprocessor
    pp_names = 'vonmises_avg fission_rate_avg dt'
    expression = '(2.0e-22 / 3600.0 ) * (fission_rate_avg / 1.0e6) * (vonmises_avg / 1.0e6) * dt'
    outputs = none
  []
  [irradiation_strain_exact]
    type = ParsedPostprocessor
    pp_names = 'irradiation_strain_inc irradiation_strain_avg_old'
    expression = 'irradiation_strain_inc + irradiation_strain_avg_old'
    outputs = none
  []
  [irradiation_strain_diff]
    type = ParsedPostprocessor
    pp_names = 'irradiation_strain_exact irradiation_strain_avg'
    expression = '(irradiation_strain_exact - irradiation_strain_avg) / irradiation_strain_exact'
    outputs = none
  []
  [irradiation_strain_diff_max]
    type = TimeExtremeValue
    postprocessor = 'irradiation_strain_diff'
    value_type = 'abs_max'
  []

  [creep_rate_exact]
    type = ParsedPostprocessor
    pp_names = 'thermal_strain_inc irradiation_strain_inc dt'
    expression = '(thermal_strain_inc + irradiation_strain_inc) / dt'
    outputs = none
  []
  [creep_rate_diff]
    type = ParsedPostprocessor
    pp_names = 'creep_rate_exact creep_rate_avg'
    expression = '(creep_rate_exact - creep_rate_avg) / creep_rate_exact'
    outputs = none
  []
  [creep_rate_max_diff]
    type = TimeExtremeValue
    postprocessor = 'creep_rate_diff'
    value_type = 'abs_max'
  []
[]

[Outputs]
  csv = true
[]

(test/tests/solid_mechanics/mn_creep/ad_cycle.i)

# Tests the ADMNCreepUpdate models with other AD thermo-mechanical models for a variety of conditions for mixed mono-carbide fuel that can be cycled through ADComputeMultipleInelasticStress.

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

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

[Variables]
  [temp]
    initial_condition = 900
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [pellets]
    add_variables = true
    strain = FINITE
    incremental = true
    generate_output = 'stress_yy strain_xx strain_yy'
    use_automatic_differentiation = true
  []
[]

[Kernels]
  [heat_ie]
    type = ADTimeDerivative
    variable = temp
  []
[]

[BCs]
  [Pressure]
    [u_top_pull]
      boundary = top
      factor = -2e7
      use_automatic_differentiation = true
    []
  []
  [u_bottom_fix]
    type = ADDirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  []
[]

[Materials]
  [elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    youngs_modulus = 7.5e10
    poissons_ratio = 0.3
  []
  [creep_1]
    type = ADMNCreepUpdate
    base_name = 1
    temperature = temp
    porosity = 0.1
    fission_rate = 1.042e20
  []
  [creep_2]
    type = ADMNCreepUpdate
    base_name = 2
    temperature = temp
    porosity = 0.1
    fission_rate = 1.042e20
  []
  [radial_return_stress]
    type = ADComputeMultipleInelasticStress
    inelastic_models = 'creep_1 creep_2'
    cycle_models = true
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  num_steps = 10
  dt = 10000
[]

[Postprocessors]
   [average_stress_yy]
     type = ElementAverageValue
     variable = stress_yy
   []
   [average_strain_xx]
     type = ElementAverageValue
     variable = strain_xx
   []
   [average_strain_yy]
     type = ElementAverageValue
     variable = strain_yy
   []
   [disp_x]
     type = NodalVariableValue
     nodeid = 7
     variable = disp_x
  []
  [disp_y]
     type = NodalVariableValue
     nodeid = 7
     variable = disp_y
  []
[]

[Outputs]
  exodus = true
[]

(assessment/nitride/EBRII/SP1/analysis/SP1_base.i)

fuel_radius = '${fparse fuel_diameter / 2}' # m
fuel_volume = '${fparse pi * fuel_radius^2 * fuel_height}' # m
cladding_radial_gap = '${fparse diametral_gap / 2}'
pin_height = ${fparse fuel_height + plenum_height + 2 * cladding_bottom_top_plug_length} # m

R = 8.31446261815324
A_U = 0.238 # kg/mol
A_Pu = 0.239 # kg/mol
A_N = 0.014 # kg/mol
avo = 6.0221408e23 # atoms per mole
th_density = 14400 # kg/m3
initial_porosity = '${fparse 1.0 - fraction_th_density}'
density = '${fparse fraction_th_density * th_density}' # kg/m3
M_avg = '${fparse x_N * A_N + x_Pu * A_Pu + (1.0 - x_Pu) * A_U}' # kg / mol
atoms_heavy_metal_per_volume = '${fparse density / M_avg * avo}' # mol / m3

# Power history
avg_lin_power = '${fparse specific_power * density * pi * fuel_radius^2}' # W/m
time_end_ramp_up = '${fparse 5 * 3600}' # s, arbitrary 5 hour ramp
time_start_ramp_down = '${fparse time_end_ramp_up + run_time}' # s
time_end_ramp_down = '${fparse time_start_ramp_down + time_end_ramp_up}' # s
total_time = '${fparse time_end_ramp_down + 3600}' # s

[GlobalParams]
  order = FIRST
  energy_per_fission = 3.412e-11 # J/fission
  displacements = 'disp_x disp_y'
  value_range_behavior = IGNORE
  min_damping = 1e-4
  volumetric_locking_correction = true
  absolute_value_vector_tags = ref
[]

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

[Mesh]
  coord_type = RZ
  [smeared_pellet_mesh]
    type = FuelPinMeshGenerator
    clad_thickness = ${cladding_thickness}
    pellet_outer_radius = ${fuel_radius}
    pellet_height = ${fuel_height}
    clad_top_gap_height = ${plenum_height}
    clad_gap_width = ${cladding_radial_gap}
    top_clad_height = ${cladding_bottom_top_plug_length}
    bottom_clad_height = ${cladding_bottom_top_plug_length}
    clad_bot_gap_height = ${gap_bottom_length}
    liner_thickness = ${liner_thickness}
    clad_mesh_density = customize
    pellet_mesh_density = customize
    nx_p = 6
    ny_p = 150
    nx_c = 4
    ny_c = 150
    ny_cu = 3
    ny_cl = 3
    nx_liner = 2
    pellet_quantity = 1
    elem_type = QUAD4
  []
  patch_size = 30
  patch_update_strategy = always
  partitioner = centroid
  centroid_partitioner_direction = y
[]

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

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 ${time_end_ramp_up} ${time_start_ramp_down} ${time_end_ramp_down}'
    y = '0 ${avg_lin_power} ${avg_lin_power} 0'
  []
  [axial_power_function] # estimated, actual peaking unknown
    type = PowerPeakingFunction
    fit = EBRII_ROW_3
    pellet_length = '${fuel_height}'
    pellet_y_start = ${fparse cladding_bottom_top_plug_length + gap_bottom_length}
  []
  [coolant_wall_temperature]
    # This is fit from coolant wall temperature provided in Dutt, Cox, Millhollen, "Performance of Refractory Alloy-clad Fuel Pins" (1984)
    type = ParsedFunction
    expression = 'full_temp := 2.46e7*y^4 - 4.836e6*y^3 + 1.09e5*y^2 + 1.19e4*y + 1.13e3;
                  if(t < ${time_end_ramp_up}, ${initial_temperature} + t * (full_temp - ${initial_temperature}) / (${time_end_ramp_up}), if(t < ${time_start_ramp_down}, full_temp, if(t < ${time_end_ramp_down}, full_temp + (t - ${time_start_ramp_down}) * (${initial_temperature} - full_temp) / (${time_end_ramp_down} - ${time_start_ramp_down}), ${initial_temperature})))'
  []
  [gas_diffusivity_function]
    # x corresponds to temperature [K] and y corresponds to fission rate [fsn/m3/s]
    type = ParsedFunction
    expression = 'kBT := 1.380649e-23 / 1.602176634e-19 * x;
                  D1 := ${xe_D10} * exp(-1.0 * ${xe_Q1} / kBT);
                  D2 := (y / 1e19)^0.5 * ${xe_D20} * exp(${xe_D2Q1} / kBT + ${xe_D2Q2} / kBT / kBT);
                  D3 := 1.85e-39 * y;
                  D1 * ${D1_xe_scalar} + D2 * ${D2_xe_scalar} + D3 * ${D3_xe_scalar}'
  []
  [vacancy_diffusivity_function]
    # x corresponds to temperature [K] and y corresponds to fission rate [fsn/m3/s]
    type = ParsedFunction
    expression = 'kBT := 1.380649e-23 / 1.602176634e-19 * x;
                  D1 := ${D10} * exp(-1.0 * ${Q1} / kBT);
                  D2 := (y / 1e19)^0.5 * ${D20} * exp(${D2Q1} / kBT + ${D2Q2} / kBT / kBT);
                  D2b := (y / 1e19)^0.5 * ${D2b0} * exp(${D2bQ1} / kBT + ${D2bQ2} / kBT / kBT);
                  D1 * ${D1_scalar} + (D2 + D2b) * ${D2_scalar}'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [fuel]
    block = pellet
    strain = FINITE
    add_variables = true
    generate_output = 'firstinv_strain stress_xx stress_yy stress_zz vonmises_stress hydrostatic_stress creep_strain_xx creep_strain_yy creep_strain_zz elastic_strain_xx elastic_strain_yy elastic_strain_zz strain_xx strain_yy strain_zz'
    eigenstrain_names = 'fuel_thermal_expansion solid_swelling_eigenstrain gaseous_swelling_eigenstrain'
    temperature = temperature
    use_automatic_differentiation = true
  []
  [clad]
    block = clad
    strain = FINITE
    add_variables = true
    generate_output = 'firstinv_strain stress_xx stress_yy stress_zz vonmises_stress hydrostatic_stress elastic_strain_xx elastic_strain_yy elastic_strain_zz strain_xx strain_yy strain_zz'
    eigenstrain_names = 'clad_thermal_expansion'
    temperature = temperature
    use_automatic_differentiation = true
  []
  [liner]
    block = liner
    strain = FINITE
    add_variables = true
    generate_output = 'firstinv_strain stress_xx stress_yy stress_zz vonmises_stress hydrostatic_stress elastic_strain_xx elastic_strain_yy elastic_strain_zz strain_xx strain_yy strain_zz'
    eigenstrain_names = 'liner_thermal_expansion'
    temperature = temperature
    use_automatic_differentiation = true
  []
[]

[Kernels]
  [heat]
    type = ADHeatConduction
    variable = temperature
    block = 'pellet liner clad'
  []
  [heat_ie]
    type = ADHeatConductionTimeDerivative
    variable = temperature
    block = 'pellet liner clad'
  []
  [heat_source]
    type = ADFissionRateHeatSource
    variable = temperature
    fission_rate = fission_rate
    energy_deposited_in_fuel = 0.95
    block = 'pellet'
  []
[]

[ThermalContactMortar]
  [thermal_contact]
    secondary_variable = temperature
    primary_boundary = clad_inside_right
    secondary_boundary = pellet_outer_radial_surface
    initial_moles = initial_moles
    gas_released = fg_released
    jump_distance_model = LANNING
    plenum_pressure = plenum_pressure
    contact_pressure = mechanical_normal_lm
    use_automatic_differentiation = true
    # use mechanical contact subdomains
    primary_subdomain = mechanical_primary_subdomain
    secondary_subdomain = mechanical_secondary_subdomain
  []
[]

[Contact]
  [mechanical]
    model = frictionless
    formulation = mortar
    primary = clad_inside_right
    secondary = pellet_outer_radial_surface
    c_normal = 1e+11
    correct_edge_dropping = true
  []
[]

[BCs]
  [no_x_all]
    type = ADDirichletBC
    variable = disp_x
    boundary = centerline
    value = 0.0
  []
  [no_y_clad]
    type = ADDirichletBC
    variable = disp_y
    boundary = clad_outside_bottom
    value = 0.0
  []
  [no_y_fuel]
    type = ADDirichletBC
    variable = disp_y
    boundary = bottom_central_pellet_node
    value = 0.0
  []
  [Pressure]
    [coolantPressure]
      boundary = 'clad_outside_bottom clad_outside_right clad_outside_top'
      function = ${coolant_pressure}
      use_automatic_differentiation = true
    []
  []
  [PlenumPressure] # apply plenum pressure on clad inner walls and pellet surfaces
    [plenumPressure]
      boundary = inside_surfaces
      initial_pressure = ${initial_plenum_pressure}
      startup_time = 0
      R = ${R}
      initial_temperature = ${initial_temperature}
      temperature = ave_temperature_interior # coupling to post processor to get gas temperature approximation
      output_initial_moles = initial_moles # coupling to post processor to get initial fill gas mass
      volume = gas_volume # coupling to post processor to get gas volume
      output = plenum_pressure # coupling to post processor to output plenum/gap pressure
      use_automatic_differentiation = true
    []
  []
  [clad_outer_temperature]
    type = ADFunctionDirichletBC
    boundary = 'clad_outside_right'
    function = coolant_wall_temperature
    variable = temperature
  []
[]

[Materials]
  [fission_rate]
    type = ADFissionRate
    block = pellet
    rod_linear_power = power_history
    axial_power_profile = axial_power_function
    pellet_radius = ${fuel_radius}
    outputs = all
  []
  [burnup]
    type = ADBurnup
    block = pellet
    atoms_heavy_metal_per_volume = ${atoms_heavy_metal_per_volume}
    outputs = all
  []
  [fuel_density]
    block = pellet
    type = ADStrainAdjustedDensity
    strain_free_density = ${density}
  []
  [porosity_pp]
    type = ADParsedMaterial
    property_name = porosity_pp
    postprocessor_names = porosity_fuel_avg
    expression = porosity_fuel_avg
  []
  [fuel_thermal]
    type = ADMNThermal
    block = pellet
    temperature = temperature
    porosity = porosity_pp # This is a hack until thermal mortar can handle stateful materials
    outputs = all
  []
  [fuel_porosity]
    type = ADPorosityFromStrain
    block = pellet
    initial_porosity = ${initial_porosity}
    inelastic_strain = 'gaseous_swelling_eigenstrain'
    outputs = all
  []
  [fuel_elasticity_tensor]
    block = pellet
    type = ADMNElasticityTensor
    temperature = temperature
    use_old_porosity = true
    porosity = porosity
    output_properties = 'youngs_modulus poissons_ratio'
    outputs = all
  []
  [fuel_thermal_expansion]
    block = pellet
    type = ADMNThermalExpansionEigenstrain
    eigenstrain_name = fuel_thermal_expansion
    stress_free_temperature = ${initial_temperature}
    temperature = temperature
  []
  [fuel_radial_return_stress]
    block = pellet
    type = ADComputeMultipleInelasticStress
    inelastic_models = 'fuel_creep'
  []
  [fuel_creep]
    block = pellet
    type = ADMNCreepUpdate
    max_inelastic_increment = 1e-3
    temperature = temperature
    porosity = porosity
    fission_rate = fission_rate
    outputs = all
  []
  [burnup_swelling]
    type = ADBurnupDependentEigenstrain
    block = pellet
    eigenstrain_name = 'solid_swelling_eigenstrain'
    swelling_factor = 0.5 # 0.5% solid fission product swelling per % FIMA
    burnup = burnup
    outputs = all
  []
  [gaseous_swelling]
    type = ADParsedMaterial
    block = pellet
    material_property_names = 'deltav_v0_bubble_bulk deltav_v0_bd deltav_v0_bubble_intra_dislocation'
    property_name = 'gaseous_swelling'
    expression = 'deltav_v0_bubble_bulk + deltav_v0_bd + deltav_v0_bubble_intra_dislocation'
    outputs = all
  []
  [gaseous_swelling_eigenstrain]
    type = ADComputeVolumetricEigenstrain
    block = pellet
    volumetric_materials = 'deltav_v0_bubble_bulk deltav_v0_bd deltav_v0_bubble_intra_dislocation'
    eigenstrain_name = 'gaseous_swelling_eigenstrain'
  []
  [vacancy_GB_diffusion]
    type = ADParsedMaterial
    block = pellet
    property_name = vacancy_GB_diffusion
    coupled_variables = 'temperature'
    expression = 'kBT := 1.380649e-23 / 1.602176634e-19 * temperature;
                  ${D1_scalar} * ${D10} * exp(-1.0 * ${Q1} / kBT) * 1e6'
  []
  [fission_gas_behavior]
    type = ADUNSifgrs
    block = pellet
    temperature = temperature
    fission_rate_material = fission_rate
    ig_bubble_coarsening = WITH_COARSENING
    grain_radius_const = ${grain_radius}
    dislocation_density_material = dislocation_density
    vacancy_diffusivity_function = vacancy_diffusivity_function
    gas_diffusivity_function = gas_diffusivity_function
    outputs = all
    initial_porosity = ${initial_porosity}
    fract_yield = 0.475
    shear_modulus = shear_modulus
    dislocation_bubble_nucleation_factor = 5e5
    vacancy_GB_diffusivity = vacancy_GB_diffusion
  []
  [dislocation_density]
    type = ADParsedMaterial
    block = pellet
    property_name = dislocation_density
    expression = ${dislocation_density}
    outputs = all
  []

  [clad_elasticity_tensor]
    block = clad
    type = ADComputeIsotropicElasticityTensor
    youngs_modulus = 68.9e9
    poissons_ratio = 0.4
  []
  [clad_thermal_expansion]
    block = clad
    type = ADComputeThermalExpansionEigenstrain
    eigenstrain_name = clad_thermal_expansion
    thermal_expansion_coeff = 7.54e-6
    stress_free_temperature = ${initial_temperature}
    temperature = temperature
  []
  [clad_stress]
    block = clad
    type = ADComputeFiniteStrainElasticStress
  []
  [clad_thermal]
    type = ADHeatConductionMaterial
    block = clad
    thermal_conductivity = 41.9
    specific_heat = 270
  []
  [clad_density]
    block = clad
    type = ADStrainAdjustedDensity
    strain_free_density = 8590
  []

  [liner_elasticity_tensor]
    block = liner
    type = ADTungstenElasticityTensor
    temperature = temperature
  []
  [liner_thermal_expansion]
    block = liner
    type = ADTungstenThermalExpansionEigenstrain
    eigenstrain_name = liner_thermal_expansion
    stress_free_temperature = ${initial_temperature}
    temperature = temperature
  []
  [liner_stress]
    block = liner
    type = ADComputeFiniteStrainElasticStress
  []
  [liner_thermal]
    block = liner
    type = ADTungstenThermal
    temperature = temperature
  []
  [liner_density]
    block = liner
    type = ADStrainAdjustedDensity
    strain_free_density = 19300
  []
[]

[Dampers]
  [disp_x]
    type = MaxIncrement
    variable = disp_x
    max_increment = 1e-4
  []
  [disp_y]
    type = MaxIncrement
    variable = disp_y
    max_increment = 1e-3
  []
  [temperature]
    type = MaxIncrement
    variable = temperature
    max_increment = 50
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options = '-snes_converged_reason -ksp_converged_reason'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -mat_mffd_err -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu    superlu_dist   1e-5          NONZERO               1e-15'
  line_search = contact

  l_max_its = 60
  nl_max_its = 20
  nl_rel_tol = 5e-6
  nl_abs_tol = 5e-9
  nl_div_tol = -1
  nl_abs_div_tol = -1

  end_time = ${total_time}
  dtmin = 1
  dtmax = 5e5
  verbose = true

  automatic_scaling = true
  compute_scaling_once = false
  ignore_variables_for_autoscaling = 'thermal_contact_thermal_lm mechanical_normal_lm'

  [TimeStepper]
    type = IterationAdaptiveDT
    timestep_limiting_postprocessor = time_step_limit
    force_step_every_function_point = true
    timestep_limiting_function = power_history
    dt = 1e2
    post_function_sync_dt = 1e3
  []
[]

[Postprocessors]
  # elemental temperatures
  [temperature_fuel_avg]
    type = ElementAverageValue
    variable = temperature
    block = pellet
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_max]
    type = ElementExtremeValue
    variable = temperature
    block = pellet
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_min]
    type = ElementExtremeValue
    variable = temperature
    block = pellet
    value_type = min
    execute_on = 'initial timestep_end'
  []

  [temperature_cladding_avg]
    type = ElementAverageValue
    variable = temperature
    block = clad
    execute_on = 'initial timestep_end'
  []
  [temperature_cladding_max]
    type = ElementExtremeValue
    variable = temperature
    block = clad
    execute_on = 'initial timestep_end'
  []
  [temperature_cladding_min]
    type = ElementExtremeValue
    variable = temperature
    block = clad
    value_type = min
    execute_on = 'initial timestep_end'
  []

  [temperature_liner_avg]
    type = ElementAverageValue
    variable = temperature
    block = liner
    execute_on = 'initial timestep_end'
  []
  [temperature_liner_max]
    type = ElementExtremeValue
    variable = temperature
    block = liner
    execute_on = 'initial timestep_end'
  []
  [temperature_liner_min]
    type = ElementExtremeValue
    variable = temperature
    block = liner
    value_type = min
    execute_on = 'initial timestep_end'
  []

  [ave_temperature_interior] # average temperature of the cladding interior and all pellet exteriors
    type = SideAverageValue
    boundary = inside_surfaces
    variable = temperature
    execute_on = 'initial linear'
  []
  [temperature_fuel_centerline_avg]
    type = AxisymmetricCenterlineAverageValue
    boundary = centerline
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_centerline_max]
    type = NodalExtremeValue
    boundary = centerline
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_centerline_min]
    type = NodalExtremeValue
    boundary = centerline
    variable = temperature
    value_type = min
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_surface_avg]
    type = SideAverageValue
    boundary = pellet_outer_radial_surface
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_surface_max]
    type = NodalExtremeValue
    boundary = pellet_outer_radial_surface
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_surface_min]
    type = NodalExtremeValue
    boundary = pellet_outer_radial_surface
    variable = temperature
    value_type = min
    execute_on = 'initial timestep_end'
  []

  # stresses
  [stress_vonmises_fuel_avg]
    type = ElementAverageValue
    variable = vonmises_stress
    block = pellet
  []
  [stress_vonmises_fuel_max]
    type = ElementExtremeValue
    variable = vonmises_stress
    block = pellet
  []
  [stress_vonmises_fuel_min]
    type = ElementExtremeValue
    variable = vonmises_stress
    value_type = min
    block = pellet
  []
  [stress_hydro_fuel_avg]
    type = ElementAverageValue
    variable = hydrostatic_stress
    block = pellet
  []
  [stress_hydro_fuel_max]
    type = ElementExtremeValue
    variable = hydrostatic_stress
    block = pellet
  []
  [stress_hydro_fuel_min]
    type = ElementExtremeValue
    variable = hydrostatic_stress
    value_type = min
    block = pellet
  []

  # strain information
  [strain_solid_swelling_fuel_avg]
    type = ElementAverageValue
    variable = solid_swelling
    block = pellet
  []
  [strain_volumetric_fuel_avg]
    type = ElementAverageValue
    variable = firstinv_strain
    block = pellet
  []
  [strain_axial_fuel_avg]
    type = ParsedPostprocessor
    pp_names = 'disp_y_fuel_top_surface_avg disp_y_fuel_bottom_surface_avg'
    expression = '(disp_y_fuel_top_surface_avg - disp_y_fuel_bottom_surface_avg) / ${fuel_height}'
  []
  [disp_y_fuel_top_surface_avg]
    type = SideAverageValue
    variable = disp_y
    boundary = top_of_top_pellet
  []
  [disp_y_fuel_top_surface_max]
    type = NodalExtremeValue
    variable = disp_y
    boundary = top_of_top_pellet
  []
  [disp_y_fuel_bottom_surface_avg]
    type = SideAverageValue
    variable = disp_y
    boundary = bottom_of_bottom_pellet
  []
  [disp_y_fuel_bottom_surface_max]
    type = NodalExtremeValue
    variable = disp_y
    boundary = bottom_of_bottom_pellet
  []
  [disp_x_fuel_radial_surface_max]
    type = NodalExtremeValue
    variable = disp_x
    boundary = pellet_outer_radial_surface
  []
  [disp_x_fuel_radial_surface_avg]
    type = SideAverageValue
    variable = disp_x
    boundary = pellet_outer_radial_surface
  []
  [disp_x_cladding_radial_surface_max]
    type = NodalExtremeValue
    variable = disp_x
    boundary = clad_outside_right
  []
  [disp_x_cladding_radial_surface_avg]
    type = SideAverageValue
    variable = disp_x
    boundary = clad_outside_right
  []

  [burnup_avg]
    type = ElementAverageValue
    block = pellet
    variable = burnup
  []
  [burnup_max]
    type = ElementExtremeValue
    block = pellet
    variable = burnup
  []
  [fission_rate_avg]
    type = ElementAverageValue
    variable = fission_rate
    block = pellet
    execute_on = 'initial timestep_end'
  []
  [fission_rate_max]
    type = ElementExtremeValue
    variable = fission_rate
    block = pellet
    execute_on = 'initial timestep_end'
  []
  [dislocation_density_avg]
    type = ElementAverageValue
    variable = dislocation_density
    block = pellet
    execute_on = 'initial timestep_end'
  []
  [porosity_fuel_avg]
    type = ElementAverageValue
    variable = porosity
    block = pellet
    execute_on = 'initial timestep_end'
  []
  [porosity_fuel_max]
    type = ElementExtremeValue
    variable = porosity
    block = pellet
    execute_on = 'initial timestep_end'
  []
  [porosity_fuel_min]
    type = ElementExtremeValue
    variable = porosity
    value_type = min
    block = pellet
    execute_on = 'initial timestep_end'
  []

  [swelling_dia_percent]
    type = ParsedPostprocessor
    pp_names = 'disp_x_fuel_radial_surface_avg'
    expression = 'disp_x_fuel_radial_surface_avg / ${fparse pi * fuel_radius} * 100'
  []
  [fuel_volume]
    type = VolumePostprocessor
    block = pellet
    execute_on = 'TIMESTEP_END INITIAL'
    use_displaced_mesh = true
  []
  [swelling_vol_percent]
    type = ParsedPostprocessor
    pp_names = 'fuel_volume'
    expression = '(fuel_volume - ${fuel_volume}) / ${fuel_volume} * 100'
  []
  [gas_volume]
    type = InternalVolume
    boundary = inside_surfaces
    execute_on = 'initial linear'
  []

  # fission gas information
  [deltav_v0_bubble_bulk]
    type = ADElementAverageMaterialProperty
    mat_prop = deltav_v0_bubble_bulk
    block = pellet
  []
  [deltav_v0_bd]
    type = ADElementAverageMaterialProperty
    mat_prop = deltav_v0_bd
    block = pellet
  []
  [deltav_v0_bubble_intra_dislocation]
    type = ADElementAverageMaterialProperty
    mat_prop = deltav_v0_bubble_intra_dislocation
    block = pellet
  []
  [fg_produced]
    type = ADElementIntegralMaterialProperty
    mat_prop = gas_concentration_generated_total
    block = pellet
  []
  [fg_released]
    type = ADElementIntegralMaterialProperty
    mat_prop = gas_concentration_release_total
    block = pellet
  []
  [fgr_percent]
    type = ParsedPostprocessor
    pp_names = 'fg_released fg_produced'
    expression = 'fg_released / fg_produced * 100'
  []
  [fg_grain_boundary]
    type = ADElementIntegralMaterialProperty
    mat_prop = gas_concentration_GB_bubble_volume
    block = pellet
  []
  [fg_matrix]
    type = ADElementIntegralMaterialProperty
    mat_prop = gas_concentration_matrix_intra
    block = pellet
  []

  [fg_intra_bubble]
    type = ADElementIntegralMaterialProperty
    mat_prop = gas_concentration_bubble_intra
    block = pellet
  []
  [fg_dislocation]
    type = ADElementIntegralMaterialProperty
    mat_prop = gas_concentration_bubble_intra_dislocation
    block = pellet
  []

  [gas_conservation]
    type = ParsedPostprocessor
    pp_names = 'fg_produced fg_released fg_dislocation fg_intra_bubble fg_matrix fg_grain_boundary'
    expression = 'fg_produced - fg_released - fg_dislocation - fg_intra_bubble - fg_matrix - fg_grain_boundary'
  []
  [GBcoverage_max]
    type = ElementExtremeValue
    variable = GBCoverage
    block = pellet
  []
  [GBcoverage_min]
    type = ElementExtremeValue
    variable = GBCoverage
    value_type = min
    block = pellet
  []
  [GBcoverage_avg]
    type = ElementAverageValue
    variable = GBCoverage
    block = pellet
  []

  # extras
  [actual_time_step_limit]
    type = MaterialTimeStepPostprocessor
    outputs = none
    block = 'pellet'
  []
  [time_step_limit]
    type = ParsedPostprocessor
    expression = 'if(actual_time_step_limit > 1e6, 1e6, actual_time_step_limit)'
    pp_names = 'actual_time_step_limit'
  []
  [power_history]
    type = FunctionValuePostprocessor
    function = power_history
  []
  [coolant_temperature]
    type = FunctionValuePostprocessor
    function = coolant_wall_temperature
  []
[]

[VectorPostprocessors]
  [centerline]
    type = SideValueSampler
    boundary = centerline
    variable = 'temperature fission_rate'
    sort_by = y
    outputs = base
  []
  [cladding_surface]
    type = SideValueSampler
    boundary = clad_outside_right
    variable = 'temperature fission_rate disp_x'
    sort_by = y
    outputs = base
  []
  [midplane]
    type = LineValueSampler
    start_point = '0 ${fparse pin_height / 2} 0'
    end_point = '${fuel_radius} ${fparse pin_height / 2} 0'
    variable = 'temperature fission_rate'
    num_points = 11
    sort_by = x
    outputs = base
  []
[]

[PerformanceMetricOutputs]
  outputs = 'out console base'
[]

[Outputs]
  perf_graph = true
  sync_times = '1e2 1e3 5e3 1e4 5e4 1e5 5e5 1e6 2e6 3e6 4e6 5e6 6e6 7e6 8e6 9e6 1e7 ${time_end_ramp_up} ${time_start_ramp_down} ${time_end_ramp_down} ${total_time}'
  file_base = '${group_name}_out'
  [out]
    type = Exodus
    time_step_interval = 10
  []
  [sync]
    type = CSV
    file_base = '${group_name}_sync_out'
    show = 'burnup_avg temperature_fuel_avg temperature_fuel_max temperature_cladding_avg temperature_cladding_max temperature_liner_avg swelling_vol_percent fgr_percent fission_rate_avg fission_rate_max porosity_fuel_avg porosity_fuel_max'
    sync_only = true
  []
  [PIE_out]
    type = CSV
    file_base = '${group_name}_PIE_out'
    show = 'burnup_avg fgr_percent swelling_vol_percent'
    execute_on = 'FINAL'
  []
  [checkpoint]
    type = Checkpoint
    time_step_interval = 10
  []
  [base]
    type = CSV
    file_base = '${group_name}_csvs/${group_name}_base_out'
  []
  [console]
    type = Console
    show = 'temperature_cladding_avg temperature_cladding_max GBcoverage_max GBcoverage_avg fgr_percent gas_conservation time_step_limit time_step_size temperature_fuel_avg temperature_fuel_max stress_vonmises_fuel_max stress_hydro_fuel_max stress_hydro_fuel_min strain_axial_fuel_avg burnup_avg fission_rate_avg porosity_fuel_avg disp_x_fuel_radial_surface_max disp_x_cladding_radial_surface_max swelling_vol_percent'
  []
[]

[Debug]
  show_var_residual = 'disp_x disp_y temperature'
  show_var_residual_norms = true
[]

(assessment/nitride/JOYO/JOYO_L4C4/analysis/JOYO_Pin_base.i)

[GlobalParams]
  order = FIRST
  energy_per_fission = ${energy_per_fission}
  displacements = 'disp_x disp_y'
  value_range_behavior = IGNORE
[]

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

[Mesh]
  coord_type = RZ
  [smeared_pellet_mesh]
    type = FuelPinMeshGenerator
    clad_thickness = ${cladding_thickness}
    pellet_outer_radius = ${fuel_radius}
    pellet_height = ${fuel_height}
    clad_top_gap_height = ${plenum_height}
    clad_gap_width = ${cladding_gap_width}
    bottom_clad_height = ${cladding_bottom_top_plug_length}
    top_clad_height = ${cladding_bottom_top_plug_length}
    clad_bot_gap_height = ${gap_bottom_length}
    clad_mesh_density = customize
    pellet_mesh_density = customize
    nx_p = 6
    ny_p = 260
    nx_c = 4
    ny_c = 260
    ny_cu = 3
    ny_cl = 3
    pellet_quantity = 1
    elem_type = QUAD4
  []
  patch_size = 30
  patch_update_strategy = auto
  partitioner = centroid
  centroid_partitioner_direction = y
[]

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

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 ${time_end_ramp_up} ${time_start_ramp_down} ${time_end_ramp_down}'
    y = '0 ${avg_lin_power} ${avg_lin_power} 0'
  []
  [coolant_wall_temperature]
    type = ParsedFunction
    expression = 'if(t<${time_end_ramp_up}, ${initial_temperature} + t*(${a}*y^2 + ${b}*y + ${c} - ${initial_temperature})/(${time_end_ramp_up}),if(t<${time_start_ramp_down},${a}*y^2 + ${b}*y + ${c}, if(t<${time_end_ramp_down}, ${a}*y^2 + ${b}*y + ${c} + (t-${time_start_ramp_down})*(${initial_temperature}-(${a}*y^2 + ${b}*y + ${c}))/(${time_end_ramp_down}-${time_start_ramp_down}),${initial_temperature})))'
  []
  [coolant_press_ramp]
    type = PiecewiseLinear
    x = '0       ${run_time}'
    y = '${coolant_pressure} ${coolant_pressure}'
  []
  [plenum_pressure]
    type = ConstantFunction
    value = ${initial_plenum_pressure}
  []
  [fast_neutron_flux_function]
    type = PiecewiseLinear
    x = '0 ${time_end_ramp_up} ${run_time}'
    y = '0 ${fast_neutron_flux} ${fast_neutron_flux}'
  []
  [axial_peaking_factors]
    type = PiecewiseBilinear
    x = '${pos_1} ${pos_2} ${pos_3} ${pos_4} ${pos_5} ${pos_6} ${pos_7} ${pos_8} ${pos_9} ${pos_10} ${pos_11} ${pos_12} ${pos_13} ${pos_14} ${pos_15} ${pos_16}'
    y = '0 ${run_time}'
    z = '${pow_1} ${pow_2} ${pow_3} ${pow_4} ${pow_5} ${pow_6} ${pow_7} ${pow_8} ${pow_9} ${pow_10} ${pow_11} ${pow_12} ${pow_13} ${pow_14} ${pow_15} ${pow_16} ${pow_1} ${pow_2} ${pow_3} ${pow_4} ${pow_5} ${pow_6} ${pow_7} ${pow_8} ${pow_9} ${pow_10} ${pow_11} ${pow_12} ${pow_13} ${pow_14} ${pow_15} ${pow_16}'
    scale_factor = 1
    axis = 1
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [fuel]
    block = pellet
    strain = FINITE
    add_variables = true
    generate_output = 'firstinv_strain stress_xx stress_yy stress_zz vonmises_stress hydrostatic_stress creep_strain_xx creep_strain_yy creep_strain_zz elastic_strain_xx elastic_strain_yy elastic_strain_zz strain_xx strain_yy strain_zz'
    extra_vector_tags = 'ref'
    eigenstrain_names = 'fuel_thermal_expansion solid_swelling_eigenstrain'
    temperature = temperature
  []
  [clad]
    block = clad
    strain = FINITE
    incremental = true
    add_variables = true
    generate_output = 'firstinv_strain stress_xx stress_yy stress_zz vonmises_stress hydrostatic_stress elastic_strain_xx elastic_strain_yy elastic_strain_zz strain_xx strain_yy strain_zz'
    extra_vector_tags = 'ref'
    eigenstrain_names = 'clad_thermal_expansion'
    temperature = temperature
  []
[]

[Kernels]
  [gravity]
    type = Gravity
    variable = disp_y
    value = -9.81
    extra_vector_tags = 'ref'
    block = 'pellet clad'
  []
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
    block = 'pellet clad'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
    block = 'pellet clad'
  []
  [heat_source]
    type = FissionRateHeatSource
    variable = temperature
    fission_rate = fission_rate
    extra_vector_tags = 'ref'
    energy_deposited_in_fuel = 0.95
    block = 'pellet'
  []
[]

[ThermalContactMortar]
  [thermal_contact]
    secondary_variable = temperature
    primary_boundary = '5' # clad_inside_right
    secondary_boundary = '10' # pellet_outer_radial_surface
    initial_moles = initial_moles # coupling to a postprocessor which supplies the initial plenum/gap gas mass
  []
[]

[Contact]
  [mechanical]
    model = frictionless
    formulation = mortar
    primary = 5 # clad_inside_right
    secondary = 10 # pellet_outer_radial_surface
    c_normal = 1e+11
    correct_edge_dropping = true
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = '12' # centerline
    value = 0.0
  []
  [no_y_clad]
    type = DirichletBC
    variable = disp_y
    boundary = '1' # clad_outside_bottom
    value = 0.0
  []
  [no_y_fuel]
    type = DirichletBC
    variable = disp_y
    boundary = '20' # bottom_of_bottom_pellet
    value = 0.0
  []
  [Pressure]
    [coolantPressure]
      boundary = '1 2 3' # clad_outside_bottom, clad_outside_right, clad_outside_top
      function = coolant_press_ramp
    []
  []
  [PlenumPressure] # apply plenum pressure on clad inner walls and pellet surfaces
    [plenumPressure]
      boundary = 9 # inside_surfaces
      initial_pressure = ${initial_plenum_pressure}
      startup_time = 0
      R = ${R}
      initial_temperature = ${initial_temperature}
      temperature = ave_temperature_interior # coupling to post processor to get gas temperature approximation
      output_initial_moles = initial_moles # coupling to post processor to get initial fill gas mass
      volume = gas_volume # coupling to post processor to get gas volume
      output = plenum_pressure # coupling to post processor to output plenum/gap pressure
    []
  []
  [clad_outer_temperature]
    type = FunctionDirichletBC
    boundary = '1 2 3' # clad_outside_bottom, clad_outside_right, clad_outside_top
    function = coolant_wall_temperature
    variable = temperature
  []
[]

[Materials]
  [porosity]
    type = GenericConstantMaterial
    block = pellet
    prop_names = porosity
    prop_values = ${initial_porosity}
    outputs = all
  []
  [fission_rate]
    type = FissionRate
    block = pellet
    rod_linear_power = power_history
    axial_power_profile = axial_peaking_factors
    pellet_radius = ${fuel_radius}
    outputs = all
  []
  [fuel_thermal]
    type = MNThermal
    block = pellet
    temperature = temperature
    porosity = porosity
    outputs = all
  []
  [fuel_elasticity_tensor]
    block = pellet
    type = MNElasticityTensor
    temperature = temperature
    porosity = porosity
    output_properties = 'youngs_modulus poissons_ratio'
    outputs = all
  []
  [fuel_thermal_expansion]
    block = pellet
    type = MNThermalExpansionEigenstrain
    eigenstrain_name = fuel_thermal_expansion
    stress_free_temperature = ${initial_temperature}
    temperature = temperature
  []
  [fuel_creep]
    block = pellet
    type = MNCreepUpdate
    max_inelastic_increment = 1e-4
    temperature = temperature
    porosity = porosity
    fission_rate = fission_rate
    outputs = all
    value_range_behavior = EXCEPTION
  []
  [burnup]
    type = Burnup
    block = pellet
    atoms_heavy_metal_per_volume = ${atoms_heavy_metal_per_volume}
    outputs = all
  []
  [burnup_swelling]
    type = BurnupDependentEigenstrain
    block = pellet
    eigenstrain_name = 'solid_swelling_eigenstrain'
    burnup = burnup
    outputs = all
  []
  [fuel_radial_return_stress]
    block = pellet
    type = ComputeMultipleInelasticStress
    inelastic_models = 'fuel_creep'
  []
  [fuel_density]
    block = pellet
    type = StrainAdjustedDensity
    strain_free_density = ${density}
  []

  [clad_elasticity_tensor]
    block = clad
    type = SS316ElasticityTensor
    temperature = temperature
    elastic_constants_model = legacy_ifr
  []
  [clad_thermal_expansion]
    block = clad
    type = SS316ThermalExpansionEigenstrain
    eigenstrain_name = clad_thermal_expansion
    stress_free_temperature = ${initial_temperature}
    temperature = temperature
  []
  [fast_neutron_flux]
    type = FastNeutronFlux
    calculate_fluence = false
    block = clad
    flux_function = fast_neutron_flux_function
  []
  [clad_creep]
    block = clad
    type = SS316CreepUpdate
    fast_neutron_flux = fast_neutron_flux
    temperature = temperature
  []
  [clad_stress]
    block = clad
    type = ComputeMultipleInelasticStress
    inelastic_models = 'clad_creep'
  []
  [clad_thermal]
    type = SS316Thermal
    block = clad
    temperature = temperature
  []
  [clad_density]
    block = clad
    type = DerivativeParsedMaterial
    coupled_variables = 'temperature'
    expression = '-4.454e-5*temperature^2 - 0.4297*temperature + 8089.4'
    property_name = 'density'
  []
[]

[Dampers]
  [disp_x]
    type = MaxIncrement
    variable = disp_x
    max_increment = 1e-4
    min_damping = 1e-4
  []
  [disp_y]
    type = MaxIncrement
    variable = disp_y
    max_increment = 1e-3
    min_damping = 1e-4
  []
  [temperature]
    type = MaxIncrement
    variable = temperature
    max_increment = 50
    min_damping = 1e-4
  []
[]

[Executioner]
  type = Transient
  # With mortar contact
  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -mat_mffd_err -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu    superlu_dist   1e-5          NONZERO               1e-15'
  snesmf_reuse_base = false
  verbose = true
  line_search = 'none'

  l_max_its = 60
  nl_max_its = 20
  nl_rel_tol = 5e-6
  nl_abs_tol = 5e-9

  end_time = ${run_time}
  dtmin = 1
  dtmax = 5e4

  [TimeStepper]
    type = IterationAdaptiveDT
    timestep_limiting_postprocessor = time_step_limit
    force_step_every_function_point = true
    timestep_limiting_function = power_history
    dt = 10
  []
[]

[Postprocessors]
  # elemental temperatures
  [temperature_fuel_avg]
    type = ElementAverageValue
    variable = temperature
    block = pellet
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_max]
    type = ElementExtremeValue
    variable = temperature
    block = pellet
  []
  [temperature_fuel_min]
    type = ElementExtremeValue
    variable = temperature
    block = pellet
    value_type = min
  []

  [temperature_clad_avg]
    type = ElementAverageValue
    variable = temperature
    block = clad
    execute_on = 'initial timestep_end'
  []
  [temperature_clad_max]
    type = ElementExtremeValue
    variable = temperature
    block = clad
  []
  [temperature_clad_min]
    type = ElementExtremeValue
    variable = temperature
    block = clad
    value_type = min
  []

  [ave_temperature_interior] # average temperature of the cladding interior and all pellet exteriors
    type = SideAverageValue
    boundary = '9' # inside_surfaces
    variable = temperature
    execute_on = 'initial linear'
  []
  [temperature_fuel_centerline_avg]
    type = AxisymmetricCenterlineAverageValue
    boundary = '12' # centerline
    variable = temperature
  []
  [temperature_fuel_centerline_max]
    type = NodalExtremeValue
    boundary = '12' # centerline
    variable = temperature
  []
  [temperature_fuel_centerline_min]
    type = NodalExtremeValue
    boundary = '12' # centerline
    variable = temperature
    value_type = min
  []
  [temperature_fuel_surface_avg]
    type = SideAverageValue
    boundary = '10' # pellet_outer_radial_surface
    variable = temperature
  []
  [temperature_fuel_surface_max]
    type = NodalExtremeValue
    boundary = '10' # pellet_outer_radial_surface
    variable = temperature
  []
  [temperature_fuel_surface_min]
    type = NodalExtremeValue
    boundary = '10' # pellet_outer_radial_surface
    variable = temperature
    value_type = min
  []

  # stresses
  [stress_vonmises_fuel_avg]
    type = ElementAverageValue
    variable = vonmises_stress
    block = pellet
  []
  [stress_vonmises_fuel_max]
    type = ElementExtremeValue
    variable = vonmises_stress
    block = pellet
  []
  [stress_vonmises_fuel_min]
    type = ElementExtremeValue
    variable = vonmises_stress
    value_type = min
    block = pellet
  []
  [stress_hydro_fuel_avg]
    type = ElementAverageValue
    variable = hydrostatic_stress
    block = pellet
  []
  [stress_hydro_fuel_max]
    type = ElementExtremeValue
    variable = hydrostatic_stress
    block = pellet
  []
  [stress_hydro_fuel_min]
    type = ElementExtremeValue
    variable = hydrostatic_stress
    value_type = min
    block = pellet
  []

  # strain information
  [strain_solid_swelling_fuel_avg]
    type = ElementAverageValue
    variable = solid_swelling
    block = pellet
  []
  [strain_volumetric_fuel_avg]
    type = ElementAverageValue
    variable = firstinv_strain
    block = pellet
  []
  [strain_axial_fuel_avg]
    type = ParsedPostprocessor
    pp_names = 'disp_y_fuel_top_surface_avg disp_y_fuel_bottom_surface_avg'
    expression = '(disp_y_fuel_top_surface_avg - disp_y_fuel_bottom_surface_avg) / ${fuel_height}'
  []
  [disp_y_fuel_top_surface_avg]
    type = SideAverageValue
    variable = disp_y
    boundary = '11' # top_of_top_pellet
  []
  [disp_y_fuel_top_surface_max]
    type = NodalExtremeValue
    variable = disp_y
    boundary = '11' # top_of_top_pellet
  []
  [disp_y_fuel_bottom_surface_avg]
    type = SideAverageValue
    variable = disp_y
    boundary = '20' # bottom_of_bottom_pellet
  []
  [disp_y_fuel_bottom_surface_max]
    type = NodalExtremeValue
    variable = disp_y
    boundary = '20' # bottom_of_bottom_pellet
  []
  [disp_x_fuel_radial_surface_max]
    type = NodalExtremeValue
    variable = disp_x
    boundary = '10' # pellet_outer_radial_surface
  []
  [disp_x_fuel_radial_surface_avg]
    type = SideAverageValue
    variable = disp_x
    boundary = '10' # pellet_outer_radial_surface
  []

  [burnup_avg]
    type = ElementAverageValue
    block = pellet
    variable = burnup
  []
  [burnup_max]
    type = ElementExtremeValue
    block = pellet
    variable = burnup
  []
  [fission_rate_avg]
    type = ElementAverageValue
    variable = fission_rate
    block = pellet
  []
  [porosity_fuel_avg]
    type = ElementAverageValue
    variable = porosity
    block = pellet
  []
  [porosity_fuel_max]
    type = ElementExtremeValue
    variable = porosity
    block = pellet
  []
  [porosity_fuel_min]
    type = ElementExtremeValue
    variable = porosity
    value_type = min
    block = pellet
  []

  [swelling_dia_percent]
    type = ParsedPostprocessor
    pp_names = 'disp_x_fuel_radial_surface_avg'
    expression = 'disp_x_fuel_radial_surface_avg / 3.14159 / ${fuel_radius} * 100'
  []
  [fuel_volume]
    type = VolumePostprocessor
    block = pellet
    execute_on = 'TIMESTEP_END INITIAL'
    use_displaced_mesh = true
  []
  [swelling_vol_percent]
    type = ParsedPostprocessor
    pp_names = 'fuel_volume'
    expression = '(fuel_volume - ${fuel_volume}) / ${fuel_volume} * 100'
  []
  [gas_volume]
    type = InternalVolume
    boundary = '9' # inside_surfaces
    execute_on = 'initial linear'
  []

  # extras
  [actual_time_step_limit]
    type = MaterialTimeStepPostprocessor
    outputs = none
    block = 'pellet'
  []
  [time_step_limit]
    type = ParsedPostprocessor
    expression = 'if(actual_time_step_limit > 1e6, 1e6, actual_time_step_limit)'
    pp_names = 'actual_time_step_limit'
  []
[]

[PerformanceMetricOutputs]
  outputs = 'performance_metrics performance_metrics_sync exodus console base_out'
[]

[Outputs]
  perf_graph = true
  csv = true
  sync_times = '1e2 1e3 5e3 1e4 5e4 1e5 5e5 1e6 2.5e6 5e6 1e7 2e7 ${fparse run_time -1e3} ${run_time}'
  file_base = '${group_name}_nominal_Pin'
  [base_out]
    type = CSV
    file_base = '${group_name}_Pin_base_out'
    show = 'burnup_max temperature_fuel_avg temperature_clad_avg ave_temperature_interior swelling_vol_percent disp_x_fuel_radial_surface_max'
    sync_only = true
  []
  [checkpoint]
    type = Checkpoint
    time_step_interval =  10
  []
  [exodus]
    type = Exodus
  []
  [sync]
    type = CSV
    file_base = '${group_name}_nominal_Pin_sync'
    sync_only = true
  []
  [console]
    type = Console
    show = 'time_step_size temperature_fuel_avg temperature_fuel_centerline_max stress_vonmises_fuel_max stress_hydro_fuel_max stress_hydro_fuel_min strain_axial_fuel_avg burnup_avg fission_rate_avg porosity_fuel_avg time_step_limit disp_x_fuel_radial_surface_max disp_x_fuel_radial_surface_avg swelling_dia_percent swelling_vol_percent fuel_volume'
  []
  [performance_metrics]
    type = CSV
    file_base = '${group_name}_nominal_performance_metrics_Pin'
    show = 'simulation_alive_time number_linear_iterations number_nonlinear_iterations time_step_size total_linear_iterations total_nonlinear_iterations physical_memory_use number_dofs number_nonlinear_variables residual_compute_time jacobian_compute_time'
  []
  [performance_metrics_sync]
    type = CSV
    sync_only = true
    file_base = '${group_name}_nominal_performance_metrics_Pin_sync'
    show = 'simulation_alive_time number_linear_iterations number_nonlinear_iterations time_step_size total_linear_iterations total_nonlinear_iterations physical_memory_use number_dofs number_nonlinear_variables residual_compute_time jacobian_compute_time'
  []
[]

[Debug]
  show_var_residual = 'disp_x disp_y temperature'
  show_var_residual_norms = true
[]

(test/tests/solid_mechanics/mn_creep/ad.i)

# Tests the ADMNCreepUpdate models with other AD thermo-mechanical models

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

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

[Variables]
  [temp]
    initial_condition = 900
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [pellets]
    add_variables = true
    strain = FINITE
    incremental = true
    generate_output = 'stress_yy strain_xx strain_yy'
    use_automatic_differentiation = true
  []
[]

[Kernels]
  [heat_ie]
    type = ADTimeDerivative
    variable = temp
  []
[]

[BCs]
  [Pressure]
    [u_top_pull]
      boundary = top
      factor = -2e7
      use_automatic_differentiation = true
    []
  []
  [u_bottom_fix]
    type = ADDirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  []
[]

[Materials]
  [elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    youngs_modulus = 7.5e10
    poissons_ratio = 0.3
  []
  [creep]
    type = ADMNCreepUpdate
    temperature = temp
    porosity = 0.1
    fission_rate = 1.042e20
  []
  [radial_return_stress]
    type = ADComputeMultipleInelasticStress
    inelastic_models = 'creep'
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  num_steps = 10
  dt = 10000
[]

[Postprocessors]
   [average_stress_yy]
     type = ElementAverageValue
     variable = stress_yy
   []
   [average_strain_xx]
     type = ElementAverageValue
     variable = strain_xx
   []
   [average_strain_yy]
     type = ElementAverageValue
     variable = strain_yy
   []
   [disp_x]
     type = NodalVariableValue
     nodeid = 7
     variable = disp_x
  []
  [disp_y]
     type = NodalVariableValue
     nodeid = 7
     variable = disp_y
  []
[]

[Outputs]
  exodus = true
[]

(assessment/nitride/MTR/SNAP50/analysis/SNAP50_base.i)

fuel_radius = '${fparse fuel_diameter / 2}' # m
fuel_volume = '${fparse pi * fuel_radius^2 * fuel_height}' # m
spacer_height = '${fparse (pin_height - fuel_height) / 2 - cladding_bottom_top_plug_length}'
cladding_radial_gap = '${fparse cladding_outer_diameter / 2 - cladding_thickness - liner_thickness - fuel_radius}'
plenum_height = ${spacer_height}
gap_bottom_length = ${spacer_height}

R = 8.31446261815324
A_U = 0.238 # kg/mol
A_Pu = 0.239 # kg/mol
A_N = 0.014 # kg/mol
avo = 6.0221408e23 # atoms per mole
th_density = 14400 # kg/m3
initial_porosity = '${fparse 1.0 - fraction_th_density}'
density = '${fparse fraction_th_density * th_density}' # kg/m3
M_avg = '${fparse x_N * A_N + x_Pu * A_Pu + (1.0 - x_Pu) * A_U}' # kg / mol
atoms_heavy_metal_per_volume = '${fparse density / M_avg * avo}' # mol / m3

# Power history
avg_lin_power = '${fparse volumetric_power * pi * fuel_radius^2}' # W/m
time_end_ramp_up = '${fparse 5 * 3600}' # s, arbitrary 5 hour ramp
time_start_ramp_down = '${fparse time_end_ramp_up + run_time}' # s
time_end_ramp_down = '${fparse time_start_ramp_down + time_end_ramp_up}' # s
total_time = '${fparse time_end_ramp_down + 3600}' # s

[GlobalParams]
  order = FIRST
  energy_per_fission = 3.412e-11 # J/fission
  displacements = 'disp_x disp_y'
  value_range_behavior = IGNORE
  min_damping = 1e-4
  volumetric_locking_correction = true
  absolute_value_vector_tags = ref
[]

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

[Mesh]
  coord_type = RZ
  [smeared_pellet_mesh]
    type = FuelPinMeshGenerator
    clad_thickness = ${cladding_thickness}
    pellet_outer_radius = ${fuel_radius}
    pellet_height = ${fuel_height}
    clad_top_gap_height = ${plenum_height}
    clad_gap_width = ${cladding_radial_gap}
    top_clad_height = ${cladding_bottom_top_plug_length}
    bottom_clad_height = ${cladding_bottom_top_plug_length}
    clad_bot_gap_height = ${gap_bottom_length}
    liner_thickness = ${liner_thickness}
    clad_mesh_density = customize
    pellet_mesh_density = customize
    nx_p = 10
    ny_p = 100
    nx_c = 4
    ny_c = 100
    ny_cu = 3
    ny_cl = 3
    nx_liner = 4
    pellet_quantity = 1
    elem_type = QUAD4
  []
  patch_size = 30
  patch_update_strategy = always
  partitioner = centroid
  centroid_partitioner_direction = y
[]

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

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 ${time_end_ramp_up} ${time_start_ramp_down} ${time_end_ramp_down}'
    y = '0 ${avg_lin_power} ${avg_lin_power} 0'
  []
  [axial_power_function] # estimated, actual peaking unknown
    type = PowerPeakingFunction
    fit = EBRII_ROW_3
    pellet_length = '${fuel_height}'
    pellet_y_start = ${fparse cladding_bottom_top_plug_length + gap_bottom_length}
  []
  [coolant_wall_temperature]
    type = ParsedFunction
    expression = 'a := ${bottom_temperature} * (y - ${mid_midplane}) * (y - ${top_midplane}) / (${bottom_midplane} - ${mid_midplane}) / (${bottom_midplane} - ${top_midplane});
                  b := ${mid_temperature} * (y - ${bottom_midplane}) * (y - ${top_midplane}) / (${mid_midplane} - ${bottom_midplane}) / (${mid_midplane} - ${top_midplane});
                  c := ${top_temperature} * (y - ${bottom_midplane}) * (y - ${mid_midplane}) / (${top_midplane} - ${bottom_midplane}) / (${top_midplane} - ${mid_midplane});
                  full_temp := a + b + c;
                  if(t < ${time_end_ramp_up}, ${initial_temperature} + t * (full_temp - ${initial_temperature}) / (${time_end_ramp_up}), if(t < ${time_start_ramp_down}, full_temp, if(t < ${time_end_ramp_down}, full_temp + (t - ${time_start_ramp_down}) * (${initial_temperature} - full_temp) / (${time_end_ramp_down} - ${time_start_ramp_down}), ${initial_temperature})))'
  []
  [gas_diffusivity_function]
    # x corresponds to temperature [K] and y corresponds to fission rate [fsn/m3/s]
    type = ParsedFunction
    expression = 'kBT := 1.380649e-23 / 1.602176634e-19 * x;
                  D1 := ${xe_D10} * exp(-1.0 * ${xe_Q1} / kBT);
                  D2 := (y / 1e19)^0.5 * ${xe_D20} * exp(${xe_D2Q1} / kBT + ${xe_D2Q2} / kBT / kBT);
                  D3 := 1.85e-39 * y;
                  D1 * ${D1_xe_scalar} + D2 * ${D2_xe_scalar} + D3 * ${D3_xe_scalar}'
  []
  [vacancy_diffusivity_function]
    # x corresponds to temperature [K] and y corresponds to fission rate [fsn/m3/s]
    type = ParsedFunction
    expression = 'kBT := 1.380649e-23 / 1.602176634e-19 * x;
                  D1 := ${D10} * exp(-1.0 * ${Q1} / kBT);
                  D2 := (y / 1e19)^0.5 * ${D20} * exp(${D2Q1} / kBT + ${D2Q2} / kBT / kBT);
                  D2b := (y / 1e19)^0.5 * ${D2b0} * exp(${D2bQ1} / kBT + ${D2bQ2} / kBT / kBT);
                  D1 * ${D1_scalar} + (D2 + D2b) * ${D2_scalar}'
  []
  [radial_power_function]
    type = ParsedFunction
    symbol_names = 'a b c d'
    symbol_values = '0.88 -0.29 0.16 0.68'
    expression = 'a*(x/${fuel_radius})^3 + b*(x/${fuel_radius})^2 + c*(x/${fuel_radius}) + d'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [fuel]
    block = pellet
    strain = FINITE
    add_variables = true
    generate_output = 'firstinv_strain stress_xx stress_yy stress_zz vonmises_stress hydrostatic_stress creep_strain_xx creep_strain_yy creep_strain_zz elastic_strain_xx elastic_strain_yy elastic_strain_zz strain_xx strain_yy strain_zz'
    eigenstrain_names = 'fuel_thermal_expansion solid_swelling_eigenstrain gaseous_swelling_eigenstrain'
    temperature = temperature
    use_automatic_differentiation = true
  []
  [clad]
    block = clad
    strain = FINITE
    add_variables = true
    generate_output = 'firstinv_strain stress_xx stress_yy stress_zz vonmises_stress hydrostatic_stress elastic_strain_xx elastic_strain_yy elastic_strain_zz strain_xx strain_yy strain_zz'
    eigenstrain_names = 'clad_thermal_expansion'
    temperature = temperature
    use_automatic_differentiation = true
  []
  [liner]
    block = liner
    strain = FINITE
    add_variables = true
    generate_output = 'firstinv_strain stress_xx stress_yy stress_zz vonmises_stress hydrostatic_stress elastic_strain_xx elastic_strain_yy elastic_strain_zz strain_xx strain_yy strain_zz'
    eigenstrain_names = 'liner_thermal_expansion'
    temperature = temperature
    use_automatic_differentiation = true
  []
[]

[Kernels]
  [heat]
    type = ADHeatConduction
    variable = temperature
    block = 'pellet liner clad'
  []
  [heat_ie]
    type = ADHeatConductionTimeDerivative
    variable = temperature
    block = 'pellet liner clad'
  []
  [heat_source]
    type = ADFissionRateHeatSource
    variable = temperature
    fission_rate = fission_rate
    energy_deposited_in_fuel = 0.95
    block = 'pellet'
  []
  [gamma_heating]
    type = ADMatBodyForce
    variable = temperature
    block = 'pellet liner clad'
    material_property = gamma_heating_density
  []
[]

[ThermalContactMortar]
  [thermal_contact]
    secondary_variable = temperature
    primary_boundary = clad_inside_right
    secondary_boundary = pellet_outer_radial_surface
    initial_moles = initial_moles
    gas_released = fg_released
    jump_distance_model = LANNING
    plenum_pressure = plenum_pressure
    contact_pressure = mechanical_normal_lm
    use_automatic_differentiation = true
    # use mechanical contact subdomains
    primary_subdomain = mechanical_primary_subdomain
    secondary_subdomain = mechanical_secondary_subdomain
  []
[]

[Contact]
  [mechanical]
    model = frictionless
    formulation = mortar
    primary = clad_inside_right
    secondary = pellet_outer_radial_surface
    c_normal = 1e+11
    correct_edge_dropping = true
  []
[]

[BCs]
  [no_x_all]
    type = ADDirichletBC
    variable = disp_x
    boundary = centerline
    value = 0.0
  []
  [no_y_clad]
    type = ADDirichletBC
    variable = disp_y
    boundary = clad_outside_bottom
    value = 0.0
  []
  [no_y_fuel]
    type = ADDirichletBC
    variable = disp_y
    boundary = bottom_central_pellet_node
    value = 0.0
  []
  [Pressure]
    [coolantPressure]
      boundary = 'clad_outside_bottom clad_outside_right clad_outside_top'
      function = ${coolant_pressure}
      use_automatic_differentiation = true
    []
  []
  [PlenumPressure] # apply plenum pressure on clad inner walls and pellet surfaces
    [plenumPressure]
      boundary = inside_surfaces
      initial_pressure = ${initial_plenum_pressure}
      startup_time = 0
      R = ${R}
      initial_temperature = ${initial_temperature}
      temperature = ave_temperature_interior # coupling to post processor to get gas temperature approximation
      output_initial_moles = initial_moles # coupling to post processor to get initial fill gas mass
      volume = gas_volume # coupling to post processor to get gas volume
      additional_volumes = spacer_volume
      temperature_of_additional_volumes = ave_temperature_interior
      output = plenum_pressure # coupling to post processor to output plenum/gap pressure
      use_automatic_differentiation = true
    []
  []
  [clad_outer_temperature]
    type = ADFunctionDirichletBC
    boundary = 'clad_outside_right'
    function = coolant_wall_temperature
    variable = temperature
  []
[]

[Materials]
  [gamma_heating_density]
    type = ADParsedMaterial
    block = 'pellet clad liner'
    material_property_names = 'density'
    property_name = gamma_heating_density
    expression = 'density * 500' # kg/m3 * W/kg
  []

  [fission_rate]
    type = ADFissionRate
    block = pellet
    rod_linear_power = power_history
    axial_power_profile = axial_power_function
    radial_power_profile = radial_power_function
    pellet_radius = ${fuel_radius}
    outputs = all
  []
  [burnup]
    type = ADBurnup
    block = pellet
    atoms_heavy_metal_per_volume = ${atoms_heavy_metal_per_volume}
    outputs = all
  []
  [fuel_density]
    block = pellet
    type = ADStrainAdjustedDensity
    strain_free_density = ${density}
  []
  [porosity_pp]
    type = ADParsedMaterial
    property_name = porosity_pp
    postprocessor_names = porosity_fuel_avg
    expression = porosity_fuel_avg
  []
  [fuel_thermal]
    type = ADMNThermal
    block = pellet
    temperature = temperature
    porosity = porosity_pp # This is a hack until thermal mortar can handle stateful materials
    outputs = all
  []
  [fuel_porosity]
    type = ADPorosityFromStrain
    block = pellet
    initial_porosity = ${initial_porosity}
    inelastic_strain = 'gaseous_swelling_eigenstrain'
    outputs = all
  []
  [fuel_elasticity_tensor]
    block = pellet
    type = ADMNElasticityTensor
    temperature = temperature
    use_old_porosity = true
    porosity = porosity
    output_properties = 'youngs_modulus poissons_ratio'
    outputs = all
  []
  [fuel_thermal_expansion]
    block = pellet
    type = ADMNThermalExpansionEigenstrain
    eigenstrain_name = fuel_thermal_expansion
    stress_free_temperature = ${initial_temperature}
    temperature = temperature
  []
  [fuel_radial_return_stress]
    block = pellet
    type = ADComputeMultipleInelasticStress
    inelastic_models = 'fuel_creep'
  []
  [fuel_creep]
    block = pellet
    type = ADMNCreepUpdate
    max_inelastic_increment = 1e-3
    temperature = temperature
    porosity = porosity
    fission_rate = fission_rate
    outputs = all
  []
  [burnup_swelling]
    type = ADBurnupDependentEigenstrain
    block = pellet
    eigenstrain_name = 'solid_swelling_eigenstrain'
    swelling_factor = 0.5 # 0.5% solid fission product swelling per % FIMA
    burnup = burnup
    outputs = all
  []
  [gaseous_swelling]
    type = ADParsedMaterial
    block = pellet
    material_property_names = 'deltav_v0_bubble_bulk deltav_v0_bd deltav_v0_bubble_intra_dislocation'
    property_name = 'gaseous_swelling'
    expression = 'deltav_v0_bubble_bulk + deltav_v0_bd + deltav_v0_bubble_intra_dislocation'
    outputs = all
  []
  [gaseous_swelling_eigenstrain]
    type = ADComputeVolumetricEigenstrain
    block = pellet
    volumetric_materials = 'deltav_v0_bubble_bulk deltav_v0_bd deltav_v0_bubble_intra_dislocation'
    eigenstrain_name = 'gaseous_swelling_eigenstrain'
  []
  [vacancy_GB_diffusion]
    type = ADParsedMaterial
    block = pellet
    property_name = vacancy_GB_diffusion
    coupled_variables = 'temperature'
    expression = 'kBT := 1.380649e-23 / 1.602176634e-19 * temperature;
                  ${D1_scalar} * ${D10} * exp(-1.0 * ${Q1} / kBT) * 1e6'
  []
  [fission_gas_behavior]
    type = ADUNSifgrs
    block = pellet
    temperature = temperature
    fission_rate_material = fission_rate
    ig_bubble_coarsening = WITH_COARSENING
    grain_radius_const = ${grain_radius}
    dislocation_density_material = dislocation_density
    vacancy_diffusivity_function = vacancy_diffusivity_function
    gas_diffusivity_function = gas_diffusivity_function
    outputs = all
    initial_porosity = ${initial_porosity}
    fract_yield = 0.475
    shear_modulus = shear_modulus
    dislocation_bubble_nucleation_factor = 5e5
    vacancy_GB_diffusivity = vacancy_GB_diffusion
  []
  [dislocation_density]
    type = ADParsedMaterial
    block = pellet
    property_name = dislocation_density
    expression = ${dislocation_density}
    outputs = all
  []

  [clad_elasticity_tensor]
    block = clad
    type = ADComputeIsotropicElasticityTensor
    youngs_modulus = 68.9e9
    poissons_ratio = 0.4
  []
  [clad_thermal_expansion]
    block = clad
    type = ADComputeThermalExpansionEigenstrain
    eigenstrain_name = clad_thermal_expansion
    thermal_expansion_coeff = 7.54e-6
    stress_free_temperature = ${initial_temperature}
    temperature = temperature
  []
  [clad_stress]
    block = clad
    type = ADComputeFiniteStrainElasticStress
  []
  [clad_thermal]
    type = ADHeatConductionMaterial
    block = clad
    thermal_conductivity = 41.9
    specific_heat = 270
  []
  [clad_density]
    block = clad
    type = ADStrainAdjustedDensity
    strain_free_density = 8590
  []

  [liner_elasticity_tensor]
    block = liner
    type = ADTungstenElasticityTensor
    temperature = temperature
  []
  [liner_thermal_expansion]
    block = liner
    type = ADTungstenThermalExpansionEigenstrain
    eigenstrain_name = liner_thermal_expansion
    stress_free_temperature = ${initial_temperature}
    temperature = temperature
  []
  [liner_stress]
    block = liner
    type = ADComputeFiniteStrainElasticStress
  []
  [liner_thermal]
    block = liner
    type = ADTungstenThermal
    temperature = temperature
  []
  [liner_density]
    block = liner
    type = ADStrainAdjustedDensity
    strain_free_density = 19300
  []
[]

[Dampers]
  [disp_x]
    type = MaxIncrement
    variable = disp_x
    max_increment = 1e-4
  []
  [disp_y]
    type = MaxIncrement
    variable = disp_y
    max_increment = 1e-3
  []
  [temperature]
    type = MaxIncrement
    variable = temperature
    max_increment = 50
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options = '-snes_converged_reason -ksp_converged_reason'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -mat_mffd_err -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu    superlu_dist   1e-5          NONZERO               1e-15'
  line_search = contact

  l_max_its = 60
  nl_max_its = 20
  nl_rel_tol = 5e-6
  nl_abs_tol = 5e-9
  nl_div_tol = -1
  nl_abs_div_tol = -1

  end_time = ${total_time}
  dtmin = 1
  dtmax = 5e5
  verbose = true

  automatic_scaling = true
  compute_scaling_once = false
  ignore_variables_for_autoscaling = 'thermal_contact_thermal_lm mechanical_normal_lm'

  [TimeStepper]
    type = IterationAdaptiveDT
    timestep_limiting_postprocessor = time_step_limit
    force_step_every_function_point = true
    timestep_limiting_function = power_history
    dt = 1e2
    post_function_sync_dt = 1e3
  []
[]

[Postprocessors]
  # elemental temperatures
  [temperature_fuel_avg]
    type = ElementAverageValue
    variable = temperature
    block = pellet
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_max]
    type = ElementExtremeValue
    variable = temperature
    block = pellet
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_min]
    type = ElementExtremeValue
    variable = temperature
    block = pellet
    value_type = min
    execute_on = 'initial timestep_end'
  []

  [temperature_cladding_avg]
    type = ElementAverageValue
    variable = temperature
    block = clad
    execute_on = 'initial timestep_end'
  []
  [temperature_cladding_max]
    type = ElementExtremeValue
    variable = temperature
    block = clad
    execute_on = 'initial timestep_end'
  []
  [temperature_cladding_min]
    type = ElementExtremeValue
    variable = temperature
    block = clad
    value_type = min
    execute_on = 'initial timestep_end'
  []

  [temperature_liner_avg]
    type = ElementAverageValue
    variable = temperature
    block = liner
    execute_on = 'initial timestep_end'
  []
  [temperature_liner_max]
    type = ElementExtremeValue
    variable = temperature
    block = liner
    execute_on = 'initial timestep_end'
  []
  [temperature_liner_min]
    type = ElementExtremeValue
    variable = temperature
    block = liner
    value_type = min
    execute_on = 'initial timestep_end'
  []

  [ave_temperature_interior] # average temperature of the cladding interior and all pellet exteriors
    type = SideAverageValue
    boundary = inside_surfaces
    variable = temperature
    execute_on = 'initial linear'
  []
  [temperature_fuel_centerline_avg]
    type = AxisymmetricCenterlineAverageValue
    boundary = centerline
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_centerline_max]
    type = NodalExtremeValue
    boundary = centerline
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_centerline_min]
    type = NodalExtremeValue
    boundary = centerline
    variable = temperature
    value_type = min
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_surface_avg]
    type = SideAverageValue
    boundary = pellet_outer_radial_surface
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_surface_max]
    type = NodalExtremeValue
    boundary = pellet_outer_radial_surface
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_surface_min]
    type = NodalExtremeValue
    boundary = pellet_outer_radial_surface
    variable = temperature
    value_type = min
    execute_on = 'initial timestep_end'
  []

  # stresses
  [stress_vonmises_fuel_avg]
    type = ElementAverageValue
    variable = vonmises_stress
    block = pellet
  []
  [stress_vonmises_fuel_max]
    type = ElementExtremeValue
    variable = vonmises_stress
    block = pellet
  []
  [stress_vonmises_fuel_min]
    type = ElementExtremeValue
    variable = vonmises_stress
    value_type = min
    block = pellet
  []
  [stress_hydro_fuel_avg]
    type = ElementAverageValue
    variable = hydrostatic_stress
    block = pellet
  []
  [stress_hydro_fuel_max]
    type = ElementExtremeValue
    variable = hydrostatic_stress
    block = pellet
  []
  [stress_hydro_fuel_min]
    type = ElementExtremeValue
    variable = hydrostatic_stress
    value_type = min
    block = pellet
  []

  # strain information
  [strain_solid_swelling_fuel_avg]
    type = ElementAverageValue
    variable = solid_swelling
    block = pellet
  []
  [strain_volumetric_fuel_avg]
    type = ElementAverageValue
    variable = firstinv_strain
    block = pellet
  []
  [strain_axial_fuel_avg]
    type = ParsedPostprocessor
    pp_names = 'disp_y_fuel_top_surface_avg disp_y_fuel_bottom_surface_avg'
    expression = '(disp_y_fuel_top_surface_avg - disp_y_fuel_bottom_surface_avg) / ${fuel_height}'
  []
  [disp_y_fuel_top_surface_avg]
    type = SideAverageValue
    variable = disp_y
    boundary = top_of_top_pellet
  []
  [disp_y_fuel_top_surface_max]
    type = NodalExtremeValue
    variable = disp_y
    boundary = top_of_top_pellet
  []
  [disp_y_fuel_bottom_surface_avg]
    type = SideAverageValue
    variable = disp_y
    boundary = bottom_of_bottom_pellet
  []
  [disp_y_fuel_bottom_surface_max]
    type = NodalExtremeValue
    variable = disp_y
    boundary = bottom_of_bottom_pellet
  []
  [disp_x_fuel_radial_surface_max]
    type = NodalExtremeValue
    variable = disp_x
    boundary = pellet_outer_radial_surface
  []
  [disp_x_fuel_radial_surface_avg]
    type = SideAverageValue
    variable = disp_x
    boundary = pellet_outer_radial_surface
  []
  [disp_x_cladding_radial_surface_max]
    type = NodalExtremeValue
    variable = disp_x
    boundary = clad_outside_right
  []
  [disp_x_cladding_radial_surface_avg]
    type = SideAverageValue
    variable = disp_x
    boundary = clad_outside_right
  []

  [burnup_avg]
    type = ElementAverageValue
    block = pellet
    variable = burnup
  []
  [burnup_max]
    type = ElementExtremeValue
    block = pellet
    variable = burnup
  []
  [fission_rate_avg]
    type = ElementAverageValue
    variable = fission_rate
    block = pellet
    execute_on = 'initial timestep_end'
  []
  [fission_rate_max]
    type = ElementExtremeValue
    variable = fission_rate
    block = pellet
    execute_on = 'initial timestep_end'
  []
  [dislocation_density_avg]
    type = ElementAverageValue
    variable = dislocation_density
    block = pellet
    execute_on = 'initial timestep_end'
  []
  [porosity_fuel_avg]
    type = ElementAverageValue
    variable = porosity
    block = pellet
    execute_on = 'initial timestep_end'
  []
  [porosity_fuel_max]
    type = ElementExtremeValue
    variable = porosity
    block = pellet
    execute_on = 'initial timestep_end'
  []
  [porosity_fuel_min]
    type = ElementExtremeValue
    variable = porosity
    value_type = min
    block = pellet
    execute_on = 'initial timestep_end'
  []

  [swelling_dia_percent]
    type = ParsedPostprocessor
    pp_names = 'disp_x_fuel_radial_surface_avg'
    expression = 'disp_x_fuel_radial_surface_avg / ${fparse pi * fuel_radius} * 100'
  []
  [fuel_volume]
    type = VolumePostprocessor
    block = pellet
    execute_on = 'TIMESTEP_END INITIAL'
    use_displaced_mesh = true
  []
  [swelling_vol_percent]
    type = ParsedPostprocessor
    pp_names = 'fuel_volume'
    expression = '(fuel_volume - ${fuel_volume}) / ${fuel_volume} * 100'
  []
  [gas_volume]
    type = InternalVolume
    boundary = inside_surfaces
    execute_on = 'initial linear'
  []
  [spacer_volume]
    type = ParsedPostprocessor
    pp_names = ''
    expression = '${fparse - pi / 4 * (spacer_diameter^2 - (spacer_diameter - spacer_thickness * 2))^2 * spacer_height * 2}'
    outputs = none
  []

  # fission gas information
  [deltav_v0_bubble_bulk]
    type = ADElementAverageMaterialProperty
    mat_prop = deltav_v0_bubble_bulk
    block = pellet
  []
  [deltav_v0_bd]
    type = ADElementAverageMaterialProperty
    mat_prop = deltav_v0_bd
    block = pellet
  []
  [deltav_v0_bubble_intra_dislocation]
    type = ADElementAverageMaterialProperty
    mat_prop = deltav_v0_bubble_intra_dislocation
    block = pellet
  []
  [fg_produced]
    type = ADElementIntegralMaterialProperty
    mat_prop = gas_concentration_generated_total
    block = pellet
  []
  [fg_released]
    type = ADElementIntegralMaterialProperty
    mat_prop = gas_concentration_release_total
    block = pellet
  []
  [fgr_percent]
    type = ParsedPostprocessor
    pp_names = 'fg_released fg_produced'
    expression = 'fg_released / fg_produced * 100'
  []
  [fg_grain_boundary]
    type = ADElementIntegralMaterialProperty
    mat_prop = gas_concentration_GB_bubble_volume
    block = pellet
  []
  [fg_matrix]
    type = ADElementIntegralMaterialProperty
    mat_prop = gas_concentration_matrix_intra
    block = pellet
  []

  [fg_intra_bubble]
    type = ADElementIntegralMaterialProperty
    mat_prop = gas_concentration_bubble_intra
    block = pellet
  []
  [fg_dislocation]
    type = ADElementIntegralMaterialProperty
    mat_prop = gas_concentration_bubble_intra_dislocation
    block = pellet
  []

  [gas_conservation]
    type = ParsedPostprocessor
    pp_names = 'fg_produced fg_released fg_dislocation fg_intra_bubble fg_matrix fg_grain_boundary'
    expression = 'fg_produced - fg_released - fg_dislocation - fg_intra_bubble - fg_matrix - fg_grain_boundary'
  []
  [GBcoverage_max]
    type = ElementExtremeValue
    variable = GBCoverage
    block = pellet
  []
  [GBcoverage_min]
    type = ElementExtremeValue
    variable = GBCoverage
    value_type = min
    block = pellet
  []
  [GBcoverage_avg]
    type = ElementAverageValue
    variable = GBCoverage
    block = pellet
  []

  # extras
  [actual_time_step_limit]
    type = MaterialTimeStepPostprocessor
    outputs = none
    block = 'pellet'
  []
  [time_step_limit]
    type = ParsedPostprocessor
    expression = 'if(actual_time_step_limit > 1e6, 1e6, actual_time_step_limit)'
    pp_names = 'actual_time_step_limit'
  []
  [power_history]
    type = FunctionValuePostprocessor
    function = power_history
  []
  [coolant_temperature]
    type = FunctionValuePostprocessor
    function = coolant_wall_temperature
  []
[]

[VectorPostprocessors]
  [centerline]
    type = SideValueSampler
    boundary = centerline
    variable = 'temperature fission_rate'
    sort_by = y
    outputs = base
  []
  [cladding_surface]
    type = SideValueSampler
    boundary = clad_outside_right
    variable = 'temperature fission_rate disp_x'
    sort_by = y
    outputs = base
  []
  [midplane]
    type = LineValueSampler
    start_point = '0 ${fparse pin_height / 2} 0'
    end_point = '${fuel_radius} ${fparse pin_height / 2} 0'
    variable = 'temperature fission_rate'
    num_points = 11
    sort_by = x
    outputs = base
  []
[]

[PerformanceMetricOutputs]
  outputs = 'out console base'
[]

[Outputs]
  perf_graph = true
  sync_times = '1e2 1e3 5e3 1e4 5e4 1e5 5e5 1e6 2e6 3e6 4e6 5e6 6e6 7e6 8e6 9e6 1e7 ${time_end_ramp_up} ${time_start_ramp_down} ${time_end_ramp_down} ${total_time}'
  file_base = '${group_name}_out'
  [out]
    type = Exodus
    time_step_interval = 10
  []
  [sync]
    type = CSV
    file_base = '${group_name}_sync_out'
    show = 'burnup_avg temperature_fuel_avg temperature_fuel_max temperature_cladding_avg temperature_cladding_max temperature_liner_avg swelling_vol_percent fgr_percent fission_rate_avg fission_rate_max porosity_fuel_avg porosity_fuel_max'
    sync_only = true
  []
  [PIE_out]
    type = CSV
    file_base = '${group_name}_PIE_out'
    show = 'burnup_avg fgr_percent swelling_vol_percent'
    execute_on = 'FINAL'
  []
  [checkpoint]
    type = Checkpoint
    time_step_interval = 10
  []
  [base]
    type = CSV
    file_base = '${group_name}_csvs/${group_name}_base_out'
  []
  [console]
    type = Console
    show = 'temperature_cladding_avg temperature_cladding_max GBcoverage_max GBcoverage_avg fgr_percent gas_conservation time_step_limit time_step_size temperature_fuel_avg temperature_fuel_max stress_vonmises_fuel_max stress_hydro_fuel_max stress_hydro_fuel_min strain_axial_fuel_avg burnup_avg fission_rate_avg porosity_fuel_avg disp_x_fuel_radial_surface_max disp_x_cladding_radial_surface_max swelling_vol_percent'
  []
[]

[Debug]
  show_var_residual = 'disp_x disp_y temperature'
  show_var_residual_norms = true
[]

(assessment/nitride/EBRII/SP1/analysis/SP1_Fuel_focused_base.i)

[GlobalParams]
  order = FIRST
  energy_per_fission = 3.412e-11 # J/fission
  displacements = 'disp_x disp_y'
  value_range_behavior = WARN
[]

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

[Mesh]
  coord_type = RZ
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmax = ${fuel_radius}
    ymax = ${fuel_height}
    nx = 10
    ny = 10
  []
  [rename]
    type = RenameBlockGenerator
    input = gen
    old_block = 0
    new_block = fuel
  []
  [corner]
    type = ExtraNodesetGenerator
    input = rename
    nodes = '0 0 0'
    new_boundary = 'bottom_corner'
  []
[]

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

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 ${time_end_ramp_up} ${time_start_ramp_down} ${time_end_ramp_down}'
    y = '0 ${avg_lin_power} ${avg_lin_power} 0'
  []
  [fuel_surface_temperature_fcn]
    type = PiecewiseLinear
    x = '0 ${time_end_ramp_up} ${time_start_ramp_down} ${time_end_ramp_down}'
    y = '${initial_temperature} ${fuel_surface_temperature} ${fuel_surface_temperature} ${initial_temperature}'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [fuel]
    block = fuel
    strain = FINITE
    add_variables = true
    generate_output = 'firstinv_strain stress_xx stress_yy stress_zz vonmises_stress hydrostatic_stress creep_strain_xx creep_strain_yy creep_strain_zz elastic_strain_xx elastic_strain_yy elastic_strain_zz strain_xx strain_yy strain_zz'
    extra_vector_tags = 'ref'
    eigenstrain_names = 'fuel_thermal_expansion solid_swelling_eigenstrain'
    temperature = temperature
    use_automatic_differentiation = true
  []
[]

[Kernels]
  [heat]
    type = ADHeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_ie]
    type = ADHeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_source]
    type = ADFissionRateHeatSource
    variable = temperature
    fission_rate = fission_rate
    extra_vector_tags = 'ref'
    energy_deposited_in_fuel = 0.95
  []
[]

[BCs]
  [no_x_fuel]
    type = ADDirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0.0
  []
  [no_y_fuel]
    type = ADDirichletBC
    variable = disp_y
    boundary = 'bottom_corner'
    value = 0.0
  []
  [temperature]
    type = ADFunctionDirichletBC
    variable = temperature
    boundary = right
    function = fuel_surface_temperature_fcn
  []
[]

[Materials]
  [porosity]
    type = ADGenericConstantMaterial
    block = fuel
    prop_names = porosity
    prop_values = ${initial_porosity}
    outputs = all
  []
  [fission_rate]
    type = ADFissionRate
    block = fuel
    rod_linear_power = power_history
    axial_power_profile = 1
    pellet_radius = ${fuel_radius}
    outputs = all
  []
  [fuel_thermal]
    type = ADMNThermal
    block = fuel
    temperature = temperature
    porosity = porosity
    outputs = all
  []
  [fuel_elasticity_tensor]
    block = fuel
    type = ADMNElasticityTensor
    temperature = temperature
    porosity = porosity
    output_properties = 'youngs_modulus poissons_ratio'
    outputs = all
  []
  [fuel_thermal_expansion]
    block = fuel
    type = ADMNThermalExpansionEigenstrain
    eigenstrain_name = fuel_thermal_expansion
    stress_free_temperature = ${initial_temperature}
    temperature = temperature
  []
  [fuel_creep]
    block = fuel
    type = ADMNCreepUpdate
    max_inelastic_increment = 1e-4
    temperature = temperature
    porosity = porosity
    fission_rate = fission_rate
    outputs = all
  []
  [burnup]
    type = ADBurnup
    block = fuel
    atoms_heavy_metal_per_volume = ${atoms_heavy_metal_per_volume}
    outputs = all
  []
  [burnup_swelling]
    type = ADBurnupDependentEigenstrain
    block = fuel
    eigenstrain_name = 'solid_swelling_eigenstrain'
    burnup = burnup
    outputs = all
  []
  [fuel_radial_return_stress]
    block = fuel
    type = ADComputeMultipleInelasticStress
    inelastic_models = 'fuel_creep'
  []
  [fuel_density]
    block = fuel
    type = ADStrainAdjustedDensity
    strain_free_density = ${density}
  []
[]

[Dampers]
  [disp_x]
    type = MaxIncrement
    variable = disp_x
    max_increment = 1e-4
    min_damping = 1e-4
  []
  [disp_y]
    type = MaxIncrement
    variable = disp_y
    max_increment = 1e-3
    min_damping = 1e-4
  []
  [temperature]
    type = MaxIncrement
    variable = temperature
    max_increment = 50
    min_damping = 1e-4
  []
[]

[Executioner]
  type = Transient
  solve_type = 'NEWTON'
  petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -snes_force_interation -mat_mffd_err -pc_factor_shift_type'
  petsc_options_value = 'lu       superlu_dist                  1                       1e-5         NONZERO'
  line_search = 'none'
  snesmf_reuse_base = false
  verbose = true

  l_max_its = 60
  nl_max_its = 20
  nl_rel_tol = 5e-6
  nl_abs_tol = 5e-9

  end_time = ${run_time}
  dtmin = 1
  dtmax = 5e4

  [TimeStepper]
    type = IterationAdaptiveDT
    timestep_limiting_postprocessor = time_step_limit
    force_step_every_function_point = true
    timestep_limiting_function = power_history
    dt = 10
  []
[]

[Postprocessors]
  # elemental temperatures
  [temperature_fuel_avg]
    type = ElementAverageValue
    variable = temperature
    block = fuel
    execute_on = 'initial timestep_end'
  []
  [temperature_fuel_max]
    type = ElementExtremeValue
    variable = temperature
    block = fuel
  []
  [temperature_fuel_min]
    type = ElementExtremeValue
    variable = temperature
    block = fuel
    value_type = min
  []

  [temperature_fuel_centerline_avg]
    type = AxisymmetricCenterlineAverageValue
    boundary = 'left'
    variable = temperature
  []
  [temperature_fuel_centerline_max]
    type = NodalExtremeValue
    boundary = 'left'
    variable = temperature
  []
  [temperature_fuel_centerline_min]
    type = NodalExtremeValue
    boundary = 'left'
    variable = temperature
    value_type = min
  []
  [temperature_fuel_surface_avg]
    type = SideAverageValue
    boundary = 'right'
    variable = temperature
  []
  [temperature_fuel_surface_max]
    type = NodalExtremeValue
    boundary = 'right'
    variable = temperature
  []
  [temperature_fuel_surface_min]
    type = NodalExtremeValue
    boundary = 'right'
    variable = temperature
    value_type = min
  []

  # stresses
  [stress_vonmises_fuel_avg]
    type = ElementAverageValue
    variable = vonmises_stress
    block = fuel
  []
  [stress_vonmises_fuel_max]
    type = ElementExtremeValue
    variable = vonmises_stress
    block = fuel
  []
  [stress_vonmises_fuel_min]
    type = ElementExtremeValue
    variable = vonmises_stress
    value_type = min
    block = fuel
  []
  [stress_hydro_fuel_avg]
    type = ElementAverageValue
    variable = hydrostatic_stress
    block = fuel
  []
  [stress_hydro_fuel_max]
    type = ElementExtremeValue
    variable = hydrostatic_stress
    block = fuel
  []
  [stress_hydro_fuel_min]
    type = ElementExtremeValue
    variable = hydrostatic_stress
    value_type = min
    block = fuel
  []

  # strain information
  [strain_solid_swelling_fuel_avg]
    type = ElementAverageValue
    variable = solid_swelling
    block = fuel
  []
  [strain_volumetric_fuel_avg]
    type = ElementAverageValue
    variable = firstinv_strain
    block = fuel
  []
  [strain_axial_fuel_avg]
    type = ParsedPostprocessor
    pp_names = 'disp_y_fuel_top_surface_avg disp_y_fuel_bottom_surface_avg'
    expression = '(disp_y_fuel_top_surface_avg - disp_y_fuel_bottom_surface_avg) / ${fuel_height}'
  []
  [disp_y_fuel_top_surface_avg]
    type = SideAverageValue
    variable = disp_y
    boundary = 'top'
  []
  [disp_y_fuel_top_surface_max]
    type = NodalExtremeValue
    variable = disp_y
    boundary = 'top'
  []
  [disp_y_fuel_bottom_surface_avg]
    type = SideAverageValue
    variable = disp_y
    boundary = 'bottom'
  []
  [disp_y_fuel_bottom_surface_max]
    type = NodalExtremeValue
    variable = disp_y
    boundary = 'bottom'
  []
  [disp_x_fuel_radial_surface_max]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 'right'
  []
  [disp_x_fuel_radial_surface_avg]
    type = SideAverageValue
    variable = disp_x
    boundary = 'right'
  []

  [burnup_avg]
    type = ElementAverageValue
    block = fuel
    variable = burnup
  []
  [burnup_max]
    type = ElementExtremeValue
    block = fuel
    variable = burnup
  []
  [fission_rate_avg]
    type = ElementAverageValue
    variable = fission_rate
    block = fuel
  []
  [porosity_fuel_avg]
    type = ElementAverageValue
    variable = porosity
    block = fuel
  []
  [porosity_fuel_max]
    type = ElementExtremeValue
    variable = porosity
    block = fuel
  []
  [porosity_fuel_min]
    type = ElementExtremeValue
    variable = porosity
    value_type = min
    block = fuel
  []

  [swelling_dia_percent]
    type = ParsedPostprocessor
    pp_names = 'disp_x_fuel_radial_surface_avg'
    expression = 'disp_x_fuel_radial_surface_avg / 3.14159 / ${fuel_radius} * 100'
  []
  [fuel_volume]
    type = VolumePostprocessor
    block = fuel
    execute_on = 'TIMESTEP_END INITIAL'
    use_displaced_mesh = true
  []
  [swelling_vol_percent]
    type = ParsedPostprocessor
    pp_names = 'fuel_volume'
    expression = '(fuel_volume - ${fuel_volume}) / ${fuel_volume} * 100'
  []

  # extras
  [actual_time_step_limit]
    type = MaterialTimeStepPostprocessor
    outputs = none
  []
  [time_step_limit]
    type = ParsedPostprocessor
    expression = 'if(actual_time_step_limit > 1e6, 1e6, actual_time_step_limit)'
    pp_names = 'actual_time_step_limit'
  []
[]

[PerformanceMetricOutputs]
  outputs = 'performance_metrics performance_metrics_sync exodus console base_out'
[]

[Outputs]
  perf_graph = true
  csv = true
  sync_times = '1e2 1e3 5e3 1e4 5e4 1e5 5e5 1e6 2.5e6 5e6 1e7 2e7 3e7 4e7 ${fparse run_time -1e3} ${run_time}'
  file_base = '${group_name}_nominal_Fuel_focused'
  [base_out]
    type = CSV
    file_base = '${group_name}_Fuel_focused_base_out'
    show = 'burnup_max temperature_fuel_avg'
    sync_only = true
  []
  [checkpoint]
    type = Checkpoint
    time_step_interval =  10
  []
  [exodus]
    type = Exodus
  []
  [sync]
    type = CSV
    file_base = '${group_name}_nominal_Fuel_focused_sync'
    sync_only = true
  []
  [console]
    type = Console
    show = 'time_step_size temperature_fuel_avg temperature_fuel_centerline_max   stress_vonmises_fuel_max stress_hydro_fuel_max stress_hydro_fuel_min strain_axial_fuel_avg burnup_avg fission_rate_avg porosity_fuel_avg time_step_limit disp_x_fuel_radial_surface_max disp_x_fuel_radial_surface_avg swelling_dia_percent swelling_vol_percent fuel_volume'
  []
  [performance_metrics]
    type = CSV
    file_base = '${group_name}_nominal_performance_metrics_Fuel_focused'
    show = 'simulation_alive_time number_linear_iterations number_nonlinear_iterations time_step_size total_linear_iterations total_nonlinear_iterations physical_memory_use number_dofs number_nonlinear_variables residual_compute_time jacobian_compute_time'
  []
  [performance_metrics_sync]
    type = CSV
    sync_only = true
    file_base = '${group_name}_nominal_performance_metrics_Fuel_focused_sync'
    show = 'simulation_alive_time number_linear_iterations number_nonlinear_iterations time_step_size total_linear_iterations total_nonlinear_iterations physical_memory_use number_dofs number_nonlinear_variables residual_compute_time jacobian_compute_time'
  []
[]

[Debug]
  show_var_residual = 'disp_x disp_y temperature'
  show_var_residual_norms = true
[]

(test/tests/solid_mechanics/mn_creep/nonad.i)

# Tests the MNCreepUpdate models with other non-AD thermo-mechanical models

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

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

[Variables]
  [temp]
    initial_condition = 900
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [pellets]
    add_variables = true
    strain = FINITE
    incremental = true
    generate_output = 'stress_yy strain_xx strain_yy'
  []
[]

[Kernels]
  [heat_ie]
    type = TimeDerivative
    variable = temp
  []
[]

[BCs]
  [Pressure]
    [u_top_pull]
      boundary = top
      factor = -2e7
    []
  []
  [u_bottom_fix]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 7.5e10
    poissons_ratio = 0.3
  []
  [creep]
    type = MNCreepUpdate
    temperature = temp
    porosity = 0.1
    fission_rate = 1.042e20
  []
  [radial_return_stress]
    type = ComputeMultipleInelasticStress
    inelastic_models = 'creep'
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  num_steps = 10
  dt = 10000
[]

[Postprocessors]
   [average_stress_yy]
     type = ElementAverageValue
     variable = stress_yy
   []
   [average_strain_xx]
     type = ElementAverageValue
     variable = strain_xx
   []
   [average_strain_yy]
     type = ElementAverageValue
     variable = strain_yy
   []
   [disp_x]
     type = NodalVariableValue
     nodeid = 7
     variable = disp_x
  []
  [disp_y]
     type = NodalVariableValue
     nodeid = 7
     variable = disp_y
  []
[]

[Outputs]
  exodus = true
[]

(assessment/nitride/EBRII/K4/analysis/base.i)

density = 13630.0
A_U = 0.238
A_Pu = 0.239
A_N = 0.014
X_Pu = 0.1
X_N = 0.5
avo = 6.0221408e23
M_avg = '${fparse X_N * A_N + X_Pu * A_Pu + (1.0 - X_N - X_Pu) * A_U}'
atoms = '${fparse density / M_avg * (1.0 - X_N) * avo}'

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

[]

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

[Mesh]
  coord_type = RZ
  use_displaced_mesh = false
  [smeared_pellet_mesh]
    type = FuelPinMeshGenerator
    clad_thickness = 0.51e-3
    pellet_outer_radius = 3.1e-3
    pellet_height = 343e-3
    clad_top_gap_height = 43.5e-3
    clad_gap_width = 0.325e-3
    bottom_clad_height = 8e-3
    top_clad_height = 8e-3
    clad_bot_gap_height = 0.2e-3 # unknown
    clad_mesh_density = customize
    pellet_mesh_density = customize
    nx_p = 6
    ny_p = 260
    nx_c = 4
    ny_c = 260
    ny_cu = 3
    ny_cl = 3
    pellet_quantity = 1
    elem_type = QUAD8
  []
  patch_size = 30
  patch_update_strategy = auto
  partitioner = centroid
  centroid_partitioner_direction = y
[]

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

[AuxVariables]
  [pellet_wall_temp]
    order = CONSTANT
    family = MONOMIAL
  []
  [radial_heat_flux]
    order = CONSTANT
    family = MONOMIAL
  []
  [clad_inner_wall_temp]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 1e5   35e6 35100000' # unknown, final burnup was 9.6 %
    y = '0 1     1    0'
  []
  [axial_peaking_factors]
    type = PiecewiseLinear
    axis = y
    x = '0    46e-3 111e-3 200e-3 248e-3 313e-3 352e-3'
    y = '73e3 73e3  82e3   85e3   82e3   74e3   74e3'
  []
  [fission_func]
    type = ParsedFunction
    symbol_names = 'axial_peaking_factors power_history'
    symbol_values = 'axial_peaking_factors power_history'
    expression = 'linear_heating_rate := axial_peaking_factors * power_history;
             pellet_cross_sectional_area := 3.14159 * pow( 3.1e-3, 2 );
             volumetric_power := linear_heating_rate / pellet_cross_sectional_area;
             energy_per_fission := 3.2e-11;
             volumetric_power / energy_per_fission'
  []

  [coolant_press_ramp]
    type = PiecewiseLinear
    x = '0       35100000'
    y = '0.151e6 0.151e6' # unknown, Na coolant
  []
  [coolant_inlet_temp]
    type = PiecewiseLinear
    x = '0   1e5 35e6 35100000'
    y = '350 644 644  350'
  []
  [coolant_outlet_temp]
    type = PiecewiseLinear
    x = '0   1e5 35e6 35100000'
    y = '350 746 746  350'
  []
  [coolant_wall_temp]
    type = ParsedFunction
    symbol_values = 'coolant_inlet_temp coolant_outlet_temp'
    symbol_names = 'coolant_inlet_temp coolant_outlet_temp'
    expression = 'rod_length := 343e-3 + 43.5e-3 + 2 * 8e-3;
             delta_temp := coolant_outlet_temp - coolant_inlet_temp;
             interpolated_temp := coolant_inlet_temp + y * delta_temp / rod_length;
             max( min( interpolated_temp, coolant_outlet_temp ), coolant_inlet_temp )'
  []
  [radial_heat_flux]
    type = ParsedFunction
    symbol_values = 'coolant_wall_temp power_history axial_peaking_factors'
    symbol_names = 'coolant_wall_temp power_history axial_peaking_factors'
    expression = 'linear_heating_rate := axial_peaking_factors * power_history;
             pellet_radius := 3.1e-3;
             pellet_circumference := 3.14159 * pellet_radius * 2;
             linear_heating_rate / pellet_circumference'
  []
  [clad_inner_wall_temp]
    type = ParsedFunction
    symbol_values = 'radial_heat_flux coolant_wall_temp'
    symbol_names = 'radial_heat_flux coolant_wall_temp'
    expression = 'conductivity_316ss := 22;
             clad_thickness := 0.51e-3;
             coolant_wall_temp + radial_heat_flux / conductivity_316ss * clad_thickness'
  []
  [pellet_wall_temp]
    type = ParsedFunction
    symbol_values = 'radial_heat_flux coolant_wall_temp clad_inner_wall_temp'
    symbol_names = 'radial_heat_flux coolant_wall_temp clad_inner_wall_temp'
    expression = 'conductivity_gap := 100;
             gap_thickness := 0.325e-3;
             clad_inner_wall_temp + radial_heat_flux / conductivity_gap * gap_thickness'
  []

  [plenum_pressure]
    type = ConstantFunction
    value = 12.4e6 # initial fill
  []

  [porosity_ramp]
    type = PiecewiseLinear
    x = '0   1e5 35e6 35100000'
    y = '0.03 0.03 0.05  0.05' # Pellets start at 96.8% TD, but no data on final porosity.
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [fuel]
    block = pellet
    strain = SMALL
    incremental = true
    extra_vector_tags = 'ref'
    eigenstrain_names = 'fuel_thermal_expansion fuel_swelling'
    use_automatic_differentiation = true
  []
  [clad]
    block = clad
    strain = SMALL
    incremental = true
    extra_vector_tags = 'ref'
    use_automatic_differentiation = true
  []
[]

[Kernels]
  [gravity]
    type = ADGravity
    variable = disp_y
    value = -9.81
    extra_vector_tags = 'ref'
  []
  [heat]
    type = ADHeatConduction
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_ie]
    type = ADHeatConductionTimeDerivative
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_source]
    type = ADFissionRateHeatSource
    variable = temp
    block = pellet
    fission_rate = fission_rate
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [pellet_wall_temp]
    type = FunctionAux
    variable = pellet_wall_temp
    function = pellet_wall_temp
  []
  [radial_heat_flux]
    type = FunctionAux
    variable = radial_heat_flux
    function = radial_heat_flux
  []
  [clad_inner_wall_temp]
    type = FunctionAux
    variable = clad_inner_wall_temp
    function = clad_inner_wall_temp
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 12
    value = 0.0
  []
  [no_y_fuel]
    type = DirichletBC
    variable = disp_y
    boundary = 20
    value = 0.0
  []
  [no_y_clad]
    type = DirichletBC
    variable = disp_y
    boundary = 1
    value = 0.0
  []
  [Pressure]
    [coolantPressure]
      boundary = '1 2 3'
      function = coolant_press_ramp
    []
    [plenumPressure]
      boundary = 9
      function = plenum_pressure
    []
  []
  [clad_outer_temp]
    type = FunctionDirichletBC
    boundary = '1 2 3'
    function = coolant_wall_temp
    variable = temp
  []
  [pellet_to_clad_cooling]
    type = FunctionDirichletBC
    boundary = 10
    function = pellet_wall_temp
    variable = temp
  []
[]

[Materials]
  [fission]
    type = ADGenericFunctionMaterial
    prop_names = fission_rate
    prop_values = fission_func
    outputs = all
  []

  [fuel_thermal]
    block = pellet
    type = ADMNThermal
    temperature = temp
    porosity = porosity
    outputs = all
  []
  [fuel_porosity]
    block = pellet
    type = ADGenericFunctionMaterial
    prop_names = porosity
    prop_values = porosity_ramp
    outputs = all
  []
  [fuel_elasticity_tensor]
    block = pellet
    type = ADMNElasticityTensor
    temperature = temp
    porosity = porosity
    output_properties = 'youngs_modulus poissons_ratio'
    outputs = all
  []
  [fuel_thermal_expansion]
    block = pellet
    type = ADMNThermalExpansionEigenstrain
    eigenstrain_name = fuel_thermal_expansion
    stress_free_temperature = 350
    temperature = temp

  []
  [fuel_creep]
    block = pellet
    type = ADMNCreepUpdate
    temperature = temp
    porosity = porosity
    fission_rate = fission_rate
    outputs = all
  []
  [burnup]
    type = ADBurnup
    block = pellet
    atoms_heavy_metal_per_volume = ${atoms}
    outputs = all
  []
  [fuel_swelling]
    block = pellet
    type = ADMNVolumetricSwellingEigenstrain
    eigenstrain_name = fuel_swelling
    temperature = temp
    initial_porosity = 0.03
    burnup = burnup
    outputs = all
  []
  [fuel_radial_return_stress]
    block = pellet
    type = ADComputeMultipleInelasticStress
    inelastic_models = 'fuel_creep'
  []
  [fuel_density]
    block = pellet
    type = ADStrainAdjustedDensity
    strain_free_density = ${density}
  []

  [clad_elasticity_tensor]
    block = clad
    type = ADComputeIsotropicElasticityTensor
    youngs_modulus = 190e9
    poissons_ratio = 0.265
  []
  [clad_stress]
    block = clad
    type = ADComputeLinearElasticStress
  []
  [clad_thermal]
    block = clad
    type = ADSS316Thermal
    temperature = temp
  []
  [clad_density]
    block = clad
    type = ADStrainAdjustedDensity
    strain_free_density = 8000
  []
[]

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

[Dampers]
  [limitT]
    type = MaxIncrement
    max_increment = 100.0
    variable = temp
  []
  [limitX]
    type = MaxIncrement
    max_increment = 1e-5
    variable = disp_x
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'
  petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -ksp_gmres_restart'
  petsc_options_value = 'lu       superlu_dist                  51'
  line_search = 'none'

  l_max_its = 50
  l_tol = 8e-3
  nl_max_its = 15
  nl_rel_tol = 1e-4
  nl_abs_tol = 1e-10

  start_time = -200
  n_startup_steps = 1
  end_time = 35100000
  dtmin = 1
  dtmax = 2e6

  [Quadrature]
    order = fifth
    side_order = seventh
  []

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 2e2
    time_t = ' 0   1e5 35e6 35100000'
    time_dt = '1e2 1e2 1e2  1e2'
    optimal_iterations = 10
    growth_factor = 2
    cutback_factor = 0.5
    iteration_window = 1
    linear_iteration_ratio = 100
  []
[]

[Postprocessors]
  [peak_burnup]
    type = ElementExtremeValue
    block = pellet
    value_type = max
    variable = burnup
  []

  [_dt]
    type = TimestepSize
  []
  [num_lin_it]
    type = NumLinearIterations
  []
  [num_nonlin_it]
    type = NumNonlinearIterations
  []
  [alive_time]
    type = PerfGraphData
    section_name = Root
    data_type = TOTAL
  []
[]

[PerformanceMetricOutputs]
  outputs = performance_metrics
[]

[Outputs]
  color = true
  exodus = true
  perf_graph = true
  [console]
    type = Console
    max_rows = 25
    time_step_interval =  1
    output_linear = true
  []
  [out]
    type = CSV
    show = 'peak_burnup'
  []
  [performance_metrics]
    type = CSV
    show = 'peak_burnup num_lin_it num_nonlin_it alive_time'
  []
[]

Description
Example Input Syntax
Input Parameters
Input Files
References