ZrDiffusivityUPuZr

Computes Fickian and Soret diffusivity.

Description

Zirconium redistribution in U-Zr and U-Pu-Zr based fuels impacts both fuel mechanical and thermal performance. For the detailed theoretical information of the calculation of the UPuZr phase diffusivity and impact on redistribution, refer to the Zirconium Redistribution page.

Numerical Implementation in BISON

Using the phase fractions calculated in PhaseUPuZr, ZrDiffusivityUPuZr calculates the Fickian and Soret diffusion coefficients using $var element = document.getElementById("moose-equation-1a92c99f-c837-4d35-9601-a86fc359d8d5");katex.render(" \\begin{aligned} D(x,T) \\quad = & \\quad{{D}_{\\pi}}(x,T) \\\\ S(x,T) \\quad = & \\quad{{D}_{\\pi}}(x,T)\\left( \\frac{x(1-{{x}_{\\text{Pu}}}-x)}{1-{{x}_{\\text{Pu}}}} \\right)\\frac{Q_{\\pi}^{\\star}}{R{{T}^{2}}} \\\\ S(x,T) \\quad = & \\quad{{f}_{\\beta}}{{D}_{\\beta }}({{X}_{\\beta}}(T),T)\\left( \\frac{{{X}_{\\beta}}(T)(1-{{x}_{\\text{Pu}}}-{{X}_{\\beta}}(T))}{1-{{x}_{\\text{Pu}}}} \\right)\\left( \\frac{\\Delta {{H}_{\\beta}}+Q_{\\beta}^{*}}{R{{T}^{2}}} \\right) \\\\ & +(1-{{f}_{\\beta }}){{D}_{\\gamma}}({{X}_{\\gamma}}(T),T)\\left( \\frac{{{X}_{\\gamma }}(T)(1-{{x}_{\\text{Pu}}}-{{X}_{\\gamma }}(T))}{1-{{x}_{\\text{Pu}}}} \\right)\\left( \\frac{\\Delta {{H}_{\\gamma }}+Q_{\\gamma}^{*}}{R{{T}^{2}}} \\right) \\end{aligned}", element, {displayMode:true,throwOnError:false});$ The artificial two-phase diffusivity parameter $var element = document.getElementById("moose-equation-85581100-e689-4c11-b688-e39789cd174d");katex.render("p", element, {displayMode:false,throwOnError:false});$ is used to set the diffusivity in the two-phase regions, and is typically set to $var element = document.getElementById("moose-equation-1c85cf72-9276-415b-a283-f264134b9d59");katex.render("p=0.2", element, {displayMode:false,throwOnError:false});$ . A scalar factor $var element = document.getElementById("moose-equation-7d9c50c3-b165-48b3-936d-88507fe5956f");katex.render("D_{0,scalar}", element, {displayMode:false,throwOnError:false});$ is also included to tune the diffusivities within each phase. It should be noted that $var element = document.getElementById("moose-equation-d79c6d87-3a34-444d-969c-a70045cd0a25");katex.render("D_{0,scalar}", element, {displayMode:false,throwOnError:false});$ is multiplied against all diffusivity values, including the Soret calculation, except for the artificial two-phase regions. In these regions, the only scalar utilized is the $var element = document.getElementById("moose-equation-6028cd9c-beec-4bd4-9c88-42741a6fc44c");katex.render("p", element, {displayMode:false,throwOnError:false});$ values.

Example Input Syntax

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [zr_diff]
    type = ZrDiffusivityUPuZr<<<{"description": "Computes Fickian and Soret diffusivity.", "href": "ZrDiffusivityUPuZr.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 0
    temperature<<<{"description": "Coupled temperature."}>>> = temp
    X_Pu<<<{"description": "Coupled plutonium atom fraction."}>>> = 0.16
    X_Zr<<<{"description": "Coupled zirconium atom fraction."}>>> = X_Zr
    outputs<<<{"description": "Vector of output names where you would like to restrict the output of variables(s) associated with this object"}>>> = 'all'
  []
[]

Input Parameters

X_PuCoupled plutonium atom fraction.
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Coupled plutonium atom fraction.
X_ZrCoupled zirconium atom fraction.
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Coupled zirconium atom fraction.
temperatureCoupled temperature.
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Coupled temperature.

Required Parameters

D0_alpha3e-06Alpha phase diffusivity prefactor
Default:3e-06
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Alpha phase diffusivity prefactor
D0_beta1.14e-05Beta phase diffusivity prefactor
Default:1.14e-05
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Beta phase diffusivity prefactor
D0_delta3e-06Delta phase diffusivity prefactor
Default:3e-06
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Delta phase diffusivity prefactor
D0_scale_alpha1Alpha diffusion coefficient multiplier
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Alpha diffusion coefficient multiplier
D0_scale_beta1Delta diffusion coefficient multiplier
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Delta diffusion coefficient multiplier
D0_scale_delta1Gamma diffusion coefficient multiplier
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Gamma diffusion coefficient multiplier
D0_scale_gamma1Beta diffusion coefficient multiplier
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Beta diffusion coefficient multiplier
Q_alpha170000Alpha phase activation energy
Default:170000
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Alpha phase activation energy
Q_beta180000Beta phase activation energy
Default:180000
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Beta phase activation energy
Q_delta150000Delta phase activation energy
Default:150000
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Delta phase activation energy
Q_scale_alpha1Alpha diffusion activation energy multiplier
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Alpha diffusion activation energy multiplier
Q_scale_beta1Delta diffusion activation energy multiplier
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Delta diffusion activation energy multiplier
Q_scale_delta1Gamma diffusion activation energy multiplier
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Gamma diffusion activation energy multiplier
Q_scale_gamma1Beta diffusion activation energy multiplier
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Beta diffusion activation energy multiplier
Qstar_alpha200000Alpha phase heat of transport
Default:200000
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Alpha phase heat of transport
Qstar_beta450000Beta phase heat of transport
Default:450000
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Beta phase heat of transport
Qstar_delta160000Delta phase heat of transport
Default:160000
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Delta phase heat of transport
Qstar_gamma-200000Gamma phase heat of transport
Default:-200000
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Gamma phase heat of transport
S_scale_alpha1Alpha soret diffusion coefficient multiplier
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Alpha soret diffusion coefficient multiplier
S_scale_beta1Delta soret diffusion coefficient multiplier
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Delta soret diffusion coefficient multiplier
S_scale_delta1Gamma soret diffusion coefficient multiplier
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Gamma soret diffusion coefficient multiplier
S_scale_gamma1Beta soret diffusion coefficient multiplier
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Beta soret diffusion coefficient multiplier
blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
computeTrueWhen false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
Default:True
C++ Type:bool
Controllable:No
Description:When false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
constant_onNONEWhen ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
Default:NONE
C++ Type:MooseEnum
Options:NONE, ELEMENT, SUBDOMAIN
Controllable:No
Description:When ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
declare_suffixAn optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
limit_twophase_soretFalseFlag to smoothly transition Soret diffusivity to zero in the two phase regions as fractional two phase concentration goes from 0.5 to 0.33
Default:False
C++ Type:bool
Controllable:No
Description:Flag to smoothly transition Soret diffusivity to zero in the two phase regions as fractional two phase concentration goes from 0.5 to 0.33
p_alpha0.2Artificial diffusion coefficient for alpha phase
Default:0.2
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Artificial diffusion coefficient for alpha phase
p_beta0.2Artificial diffusion coefficient for beta phase
Default:0.2
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Artificial diffusion coefficient for beta phase
p_delta0.2Artificial diffusion coefficient for delta phase
Default:0.2
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Artificial diffusion coefficient for delta phase
p_gamma0.2Artificial diffusion coefficient for gamma phase
Default:0.2
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Artificial diffusion coefficient for gamma phase
verboseFalsePrint diagnostic information
Default:False
C++ Type:bool
Controllable:No
Description:Print diagnostic information

Optional Parameters

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

Advanced Parameters

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

Outputs Parameters

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

Material Property Retrieval Parameters

Input Files

(test/tests/zirconium_diffusion/ad_exact.i)
(examples/metal_fuel/uzr/pin.i)
(test/tests/zrdiffusivity_upuzr/ad_patch.i)
(test/tests/zirconium_diffusion/exact.i)
(test/tests/zirconium_diffusion/ad_test.i)
(test/tests/zrdiffusivity_upuzr/exact.i)
(test/tests/zirconium_diffusion/test.i)
(test/tests/zrdiffusivity_upuzr/patch.i)

(test/tests/zrdiffusivity_upuzr/patch.i)

# Calculates Fickian and Soret Diffusivty by changing temperature and zirconium atom fraction

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

[Variables]
  [temp]
    initial_condition = 773
  []
  [X_Zr]
    initial_condition = 0.01
  []
[]

[Kernels]
  [temp_diff]
    type = Diffusion
    variable = temp
  []
  [Zr_diff]
    type = Diffusion
    variable = X_Zr
  []
[]

[AuxVariables]
  [alpha_solvus]
    order = CONSTANT
    family = MONOMIAL
  []
  [beta_solvus]
    order = CONSTANT
    family = MONOMIAL
  []
  [gamma_solvus]
    order = CONSTANT
    family = MONOMIAL
  []
  [delta_solvus]
    order = CONSTANT
    family = MONOMIAL
  []
  [alpha_2phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [beta_2phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [gamma_2phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [delta_2phase]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [alpha_solvus]
    type = MaterialStdVectorAux
    variable = alpha_solvus
    property = phase_solvus
    index = 0
  []
  [beta_solvus]
    type = MaterialStdVectorAux
    variable = beta_solvus
    property = phase_solvus
    index = 1
  []
  [gamma_solvus]
    type = MaterialStdVectorAux
    variable = gamma_solvus
    property = phase_solvus
    index = 2
  []
  [delta_solvus]
    type = MaterialStdVectorAux
    variable = delta_solvus
    property = phase_solvus
    index = 3
  []
  [alpha_2phase]
    type = MaterialStdVectorAux
    variable = alpha_2phase
    property = twophase_fractions
    index = 0
  []
  [beta_2phase]
    type = MaterialStdVectorAux
    variable = beta_2phase
    property = twophase_fractions
    index = 1
  []
  [gamma_2phase]
    type = MaterialStdVectorAux
    variable = gamma_2phase
    property = twophase_fractions
    index = 2
  []
  [delta_2phase]
    type = MaterialStdVectorAux
    variable = delta_2phase
    property = twophase_fractions
    index = 3
  []
[]

[BCs]
  [temp_BC]
    type = FunctionDirichletBC
    boundary = 'top bottom left right'
    variable = temp
    function = temp_ramp
  []
  [conc_BC]
    type = FunctionDirichletBC
    boundary = 'top bottom left right'
    variable = X_Zr
    function = conc_ramp
  []
[]

[Functions]
  [temp_ramp]
    type = PiecewiseLinear
    x = '0   50  100 150 200'
    y = '950 850 850 950 850'
  []
  [conc_ramp]
    type = PiecewiseLinear
    x = '0    50   100  150  200'
    y = '0.01 0.01 0.1  0.55 0.8'
  []
[]

[Materials]
  [phase]
    type = PhaseUPuZr
    block = 0
    temperature = temp
    X_Pu = 0.16
    X_Zr = X_Zr
    outputs = 'all'
    output_properties = 'alpha beta gamma delta phase phase_region'
    lever_weight_coeffs = '.1 .9' # different for testing, usually '0 1'
    calc_H = false
  []
  [zr_diff]
    type = ZrDiffusivityUPuZr
    block = 0
    temperature = temp
    X_Pu = 0.16
    X_Zr = X_Zr
    outputs = 'all'
  []
[]

[Executioner]
  type = Transient
  end_time = 200
[]

[Postprocessors]
  [temp]
    type = ElementAverageValue
    variable = temp
    execute_on = 'initial timestep_end'
  []
  [Zr]
    type = ElementAverageValue
    variable = X_Zr
  []

  [phase_region]
    type = ElementAverageValue
    variable = phase_region
  []

  [alpha]
    type = ElementAverageValue
    variable = alpha
  []
  [beta]
    type = ElementAverageValue
    variable = beta
  []
  [gamma]
    type = ElementAverageValue
    variable = gamma
  []
  [delta]
    type = ElementAverageValue
    variable = delta
  []

  [alpha_2phase]
    type = ElementAverageValue
    variable = alpha_2phase
  []
  [beta_2phase]
    type = ElementAverageValue
    variable = beta_2phase
  []
  [gamma_2phase]
    type = ElementAverageValue
    variable = gamma_2phase
  []
  [delta_2phase]
    type = ElementAverageValue
    variable = delta_2phase
  []

  [alpha_solvus]
    type = ElementAverageValue
    variable = alpha_solvus
  []
  [beta_solvus]
    type = ElementAverageValue
    variable = beta_solvus
  []
  [gamma_solvus]
    type = ElementAverageValue
    variable = gamma_solvus
  []
  [delta_solvus]
    type = ElementAverageValue
    variable = delta_solvus
  []
  [zr_diff]
    type = ElementAverageValue
    variable = Zr_diffusivity
  []
  [zr_soret]
    type = ElementAverageValue
    variable = Zr_Soret
  []
[]

[Outputs]
  csv = true
[]

(test/tests/zirconium_diffusion/ad_exact.i)

# This input file mirrors the method in moose/modules/misc/tests/kernels/thermo_diffusion.i

# Steady-state test for the ZirconiumDiffusion kernel.
#
# This test applies a constant temperature gradient to drive thermo-diffusion
# in the variable u. At steady state, the thermo-diffusion is balanced by
# diffusion due to Fick's Law, so the total flux is
#
#   J = -D grad(u) - S grad(T) = 0
#
# If there are no fluxes at the boundaries, then there is no background flux and
# these two terms must balance each other everywhere:
#
#   grad(u) = -(S / D) grad(T)
#
# Using S = D * u(1-u) * Qstar/R/T^2, and setting Qstar=R the dx can be eliminated to give
#
#   du = -u(1-u) dT/T^2
#
# This can be solved to give the profile for u as a function of temperature:
#
#   u = exp(T^-1 + C) / ( 1 + exp(T^-1 + C ) )
#
# Here, we are using simple heat conduction with Dirichlet boundaries on 0 <= x <= 1
# to give a linear profile for temperature: T = x + 1:
#
#   u = exp(1/(x+1) + C) / ( 1 + exp(1/(x+1) + C ) )
#
# Using the BC u(x=0) = 0.5, then C = -1.
# This analytical result is tracked by the aux variable "correct_u".
#
#
# Postprocessor Values:
# +----------------+----------------+----------------+----------------+
# | time           | correct_u      | temp           | u              |
# +----------------+----------------+----------------+----------------+
# |   0.000000e+00 |   4.243783e-01 |   1.000000e+00 |   0.000000e+00 |
# |   1.000000e+00 |   4.243783e-01 |   1.500000e+00 |   4.243658e-01 |
# +----------------+----------------+----------------+----------------+

[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = 1
    ymin = 0
    ymax = 1
    nx = 10
    ny = 1
  []
[]

[Variables]
  [u]
    initial_condition = 0.5
  []
  [temp]
    initial_condition = 1
  []
[]

[Kernels]
  [zr_diff]
    type = ADZirconiumDiffusion
    variable = u
    temperature = temp
  []

  # Heat diffusion gives a linear temperature profile to drive the Soret diffusion.
  [diffT]
    type = ADDiffusion
    variable = temp
  []
[]

[BCs]
  [left]
    type = ADDirichletBC
    variable = u
    boundary = left
    value = 0.5
  []

  [leftt]
    type = ADDirichletBC
    variable = temp
    boundary = left
    value = 1
  []
  [rightt]
    type = ADDirichletBC
    variable = temp
    boundary = right
    value = 2
  []
[]

[Materials]
  [zr_diff]
    type = ADZrDiffusivityUPuZr
    block = 0
    temperature = temp
    X_Pu = 0
    X_Zr = u
    Qstar_alpha = 8.31446
    Q_scale_alpha = 0
  []
  [phase]
    type = ADPhaseUPuZr
    block = 0
    temperature = 800
    X_Pu = 0
    X_Zr = 0.01
  []
[]

[Functions]
  # The correct (analytical) result for the solution of u
  [concentration_profile]
    type = ParsedFunction
    expression = 'exp(1.0/(x+1.0)-1) / ( 1.0 + exp(1.0/(x+1.0)-1) )'
  []
[]

[AuxVariables]
  [correct_u]
  []
[]

[AuxKernels]
  # This kernel is just being used to plot the correct result for u.
  [copy_from_function]
    type = FunctionAux
    variable = correct_u
    function = concentration_profile
    execute_on = 'initial timestep_begin'
  []
[]

[Executioner]
  type = Steady

  nl_abs_tol = 1e-50
  nl_rel_tol = 1e-12
  l_tol = 1e-10
[]

[Postprocessors]
  [u]
    type = ElementAverageValue
    variable = u
    execute_on = 'initial timestep_end'
  []
  [temp]
    type = ElementAverageValue
    variable = temp
    execute_on = 'initial timestep_end'
  []
  [correct_u]
    type = ElementAverageValue
    variable = correct_u
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  exodus = true
[]

(examples/metal_fuel/uzr/pin.i)

initial_fuel_density = 15800

[GlobalParams]
  energy_per_fission = 3.2e-11  # J/fission
  temperature = temp
[]

[Mesh]
  coord_type = RZ
  [mesh]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 300
    xmax = 0.002
  []
[]

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

[Variables]
  [temp]
    initial_condition = 298
  []
  [X_Zr]
    scaling = 1e10
    initial_condition = 0.225
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 100000 3901600 4333600 6148000 7948800  '
    y = '0 40016 40344 42312 41984 41984  '
  []
  [temp_history]
    type = PiecewiseLinear
    x = '0 100000 3901600 4333600 6148000 7948800  '
    y = '298.0 805.5 803.8 811.1 808.5 807.4'
  []
[]

[Kernels]
  [heat]
    type = ConstitutiveHeatConduction
    variable = temp
    thermal_conductivity = 'thermal_conductivity'
    thermal_conductivity_args = 'temp X_Zr'
    thermal_conductivity_derivs = 'thermal_conductivity_dT thermal_conductivity_dZr'
    extra_vector_tags = 'ref'
  []
  [heat_ie]
    type = ConstitutiveHeatConductionTimeDerivative
    variable = temp
    specific_heat = 'specific_heat'
    specific_heat_args = 'temp'
    specific_heat_derivs = 'specific_heat_dT'
    extra_vector_tags = 'ref'
  []
  [heat_source]
    type = FissionRateHeatSource
    variable = temp
    fission_rate_args = 'X_Zr'
    fission_rate_derivs = 'fission_rate_dZr'
    extra_vector_tags = 'ref'
  []
  [ZrDiffusion]
    type = ZirconiumDiffusion
    variable = X_Zr
    temperature = temp
    extra_vector_tags = 'ref'
  []
  [diffusion_ie]
    type = TimeDerivative
    variable = X_Zr
    extra_vector_tags = 'ref'
  []
[]

[BCs]
  [temp_rightside]
    type = FunctionDirichletBC
    variable = temp
    boundary = right
    function = temp_history
  []
[]

[Materials]
  [fission_rate]
    type = UPuZrFissionRate
    rod_linear_power = power_history
    axial_power_profile = 1.0
    pellet_radius = 0.002
    initial_X_Zr = 0.225
    X_Zr = X_Zr
    X_Pu_function = 0.163
  []
  [burnup]
    type = UPuZrBurnup
    initial_X_Zr = 0.225
    initial_X_Pu = 0.163
    density = ${initial_fuel_density}
  []
  [metal_fuel_thermal]
    type = UPuZrThermal
    spheat_model = savage
    thcond_model = lanl
    porosity = porosity
    X_Zr = X_Zr
    X_Pu = 0.163
  []
  [swelling]
    type = UPuZrVolumetricSwellingEigenstrain
    burnup = 0
    hydrostatic_stress = 0
    eigenstrain_name = 'swelling'
    output_properties = 'porosity gaseous_porosity'
  []
  [fuel_density]
    type = GenericConstantMaterial
    prop_names = density
    prop_values = ${initial_fuel_density}
  []
  [phase]
    type = PhaseUPuZr
    X_Zr = X_Zr
    X_Pu = 0.163
    AB_temp = 965.15
    CD_temp = 995.15
  []
  [zr_diff]
    type = ZrDiffusivityUPuZr
    X_Zr = X_Zr
    X_Pu = 0.163
    p_alpha = 0.20
    p_delta = 0.20
    p_beta  = 0.20
    p_gamma = 0.20
    D0_scale_alpha = 15.0
    D0_scale_delta = 1.0
    D0_scale_beta = 2.0
    D0_scale_gamma = 7.0
    outputs = all
  []
[]

[Dampers]
  [bounded_zr]
    type = BoundingValueNodalDamper
    variable = X_Zr
    max_value = 0.83
    min_value = 0
  []
[]

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

[Executioner]
  type = Transient

  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  line_search = 'none'

  end_time = 3e6
  dtmin = 1e-2
  dtmax = 1e6
  nl_max_its = 30
  l_max_its = 100
  nl_abs_tol = 1e-9
  nl_rel_tol = 1e-9
  l_tol = 1e-05

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e3
    optimal_iterations = 8
    iteration_window = 2
    growth_factor = 2
    cutback_factor = .5
    timestep_limiting_function = power_history
    force_step_every_function_point = true
  []
[]

[Postprocessors]
  [dt]
    type = TimestepSize
  []
  [num_lin]
    type = NumLinearIterations
  []
  [num_nonlin]
    type = NumNonlinearIterations
  []
  [ave_temp]
    type = ElementAverageValue
    variable = temp
    block = 0
  []
  [X_Zr]
    type = NodalExtremeValue
    variable = X_Zr
    block = 0
    value_type = max
  []
[]

[Outputs]
  perf_graph = true
  #exodus = true

  [csv_output]
    type = CSV
    precision = 8
    file_base = pin_out
    hide = 'dt num_lin num_nonlin'
    execute_on = final
  []
[]

(test/tests/zrdiffusivity_upuzr/ad_patch.i)

# Calculates Fickian and Soret Diffusivty by changing temperature and zirconium atom fraction

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

[Variables]
  [temp]
    initial_condition = 773
  []
  [X_Zr]
    initial_condition = 0.01
  []
[]

[Kernels]
  [temp_diff]
    type = ADDiffusion
    variable = temp
  []
  [Zr_diff]
    type = ADDiffusion
    variable = X_Zr
  []
[]

[AuxVariables]
  [alpha_solvus]
    order = CONSTANT
    family = MONOMIAL
  []
  [beta_solvus]
    order = CONSTANT
    family = MONOMIAL
  []
  [gamma_solvus]
    order = CONSTANT
    family = MONOMIAL
  []
  [delta_solvus]
    order = CONSTANT
    family = MONOMIAL
  []
  [alpha_2phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [beta_2phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [gamma_2phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [delta_2phase]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [alpha_solvus]
    type = ADMaterialStdVectorAux
    variable = alpha_solvus
    property = phase_solvus
    index = 0
  []
  [beta_solvus]
    type = ADMaterialStdVectorAux
    variable = beta_solvus
    property = phase_solvus
    index = 1
  []
  [gamma_solvus]
    type = ADMaterialStdVectorAux
    variable = gamma_solvus
    property = phase_solvus
    index = 2
  []
  [delta_solvus]
    type = ADMaterialStdVectorAux
    variable = delta_solvus
    property = phase_solvus
    index = 3
  []
  [alpha_2phase]
    type = ADMaterialStdVectorAux
    variable = alpha_2phase
    property = twophase_fractions
    index = 0
  []
  [beta_2phase]
    type = ADMaterialStdVectorAux
    variable = beta_2phase
    property = twophase_fractions
    index = 1
  []
  [gamma_2phase]
    type = ADMaterialStdVectorAux
    variable = gamma_2phase
    property = twophase_fractions
    index = 2
  []
  [delta_2phase]
    type = ADMaterialStdVectorAux
    variable = delta_2phase
    property = twophase_fractions
    index = 3
  []
[]

[BCs]
  [temp_BC]
    type = FunctionDirichletBC
    boundary = 'top bottom left right'
    variable = temp
    function = temp_ramp
  []
  [conc_BC]
    type = FunctionDirichletBC
    boundary = 'top bottom left right'
    variable = X_Zr
    function = conc_ramp
  []
[]

[Functions]
  [temp_ramp]
    type = PiecewiseLinear
    x = '0   50  100 150 200'
    y = '950 850 850 950 850'
  []
  [conc_ramp]
    type = PiecewiseLinear
    x = '0    50   100  150  200'
    y = '0.01 0.01 0.1  0.55 0.8'
  []
[]

[Materials]
  [phase]
    type = ADPhaseUPuZr
    block = 0
    temperature = temp
    X_Pu = 0.16
    X_Zr = X_Zr
    outputs = 'all'
    output_properties = 'alpha beta gamma delta phase phase_region'
    lever_weight_coeffs = '.1 .9' # different for testing, usually '0 1'
    calc_H = false
  []
  [zr_diff]
    type = ADZrDiffusivityUPuZr
    block = 0
    temperature = temp
    X_Pu = 0.16
    X_Zr = X_Zr
    outputs = 'all'
  []
[]

[Executioner]
  type = Transient
  end_time = 200
[]

[Postprocessors]
  [temp]
    type = ElementAverageValue
    variable = temp
    execute_on = 'initial timestep_end'
  []
  [Zr]
    type = ElementAverageValue
    variable = X_Zr
  []

  [phase_region]
    type = ElementAverageValue
    variable = phase_region
  []

  [alpha]
    type = ElementAverageValue
    variable = alpha
  []
  [beta]
    type = ElementAverageValue
    variable = beta
  []
  [gamma]
    type = ElementAverageValue
    variable = gamma
  []
  [delta]
    type = ElementAverageValue
    variable = delta
  []

  [alpha_2phase]
    type = ElementAverageValue
    variable = alpha_2phase
  []
  [beta_2phase]
    type = ElementAverageValue
    variable = beta_2phase
  []
  [gamma_2phase]
    type = ElementAverageValue
    variable = gamma_2phase
  []
  [delta_2phase]
    type = ElementAverageValue
    variable = delta_2phase
  []

  [alpha_solvus]
    type = ElementAverageValue
    variable = alpha_solvus
  []
  [beta_solvus]
    type = ElementAverageValue
    variable = beta_solvus
  []
  [gamma_solvus]
    type = ElementAverageValue
    variable = gamma_solvus
  []
  [delta_solvus]
    type = ElementAverageValue
    variable = delta_solvus
  []
  [zr_diff]
    type = ElementAverageValue
    variable = Zr_diffusivity
  []
  [zr_soret]
    type = ElementAverageValue
    variable = Zr_Soret
  []
[]

[Outputs]
  csv = true
[]

(test/tests/zirconium_diffusion/exact.i)

# This input file mirrors the method in moose/modules/misc/tests/kernels/thermo_diffusion.i

# Steady-state test for the ZirconiumDiffusion kernel.
#
# This test applies a constant temperature gradient to drive thermo-diffusion
# in the variable u. At steady state, the thermo-diffusion is balanced by
# diffusion due to Fick's Law, so the total flux is
#
#   J = -D grad(u) - S grad(T) = 0
#
# If there are no fluxes at the boundaries, then there is no background flux and
# these two terms must balance each other everywhere:
#
#   grad(u) = -(S / D) grad(T)
#
# Using S = D * u(1-u) * Qstar/R/T^2, and setting Qstar=R the dx can be eliminated to give
#
#   du = -u(1-u) dT/T^2
#
# This can be solved to give the profile for u as a function of temperature:
#
#   u = exp(T^-1 + C) / ( 1 + exp(T^-1 + C ) )
#
# Here, we are using simple heat conduction with Dirichlet boundaries on 0 <= x <= 1
# to give a linear profile for temperature: T = x + 1:
#
#   u = exp(1/(x+1) + C) / ( 1 + exp(1/(x+1) + C ) )
#
# Using the BC u(x=0) = 0.5, then C = -1.
# This analytical result is tracked by the aux variable "correct_u".
#
#
# Postprocessor Values:
# +----------------+----------------+----------------+----------------+
# | time           | correct_u      | temp           | u              |
# +----------------+----------------+----------------+----------------+
# |   0.000000e+00 |   4.243783e-01 |   1.000000e+00 |   0.000000e+00 |
# |   1.000000e+00 |   4.243783e-01 |   1.500000e+00 |   4.243658e-01 |
# +----------------+----------------+----------------+----------------+

[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = 1
    ymin = 0
    ymax = 1
    nx = 10
    ny = 1
  []
[]

[Variables]
  [u]
    initial_condition = 0.5
  []
  [temp]
    initial_condition = 1
  []
[]

[Kernels]
  [zr_diff]
    type = ZirconiumDiffusion
    variable = u
    temperature = temp
  []

  # Heat diffusion gives a linear temperature profile to drive the Soret diffusion.
  [diffT]
    type = Diffusion
    variable = temp
  []
[]

[BCs]
  [left]
    type = DirichletBC
    variable = u
    boundary = left
    value = 0.5
  []

  [leftt]
    type = DirichletBC
    variable = temp
    boundary = left
    value = 1
  []
  [rightt]
    type = DirichletBC
    variable = temp
    boundary = right
    value = 2
  []
[]

[Materials]
  [zr_diff]
    type = ZrDiffusivityUPuZr
    block = 0
    temperature = temp
    X_Pu = 0
    X_Zr = u
    Qstar_alpha = 8.31446
    Q_scale_alpha = 0
  []
  [phase]
    type = PhaseUPuZr
    block = 0
    temperature = 800
    X_Pu = 0
    X_Zr = 0.01
  []
[]

[Functions]
  # The correct (analytical) result for the solution of u
  [concentration_profile]
    type = ParsedFunction
    expression = 'exp(1.0/(x+1.0)-1) / ( 1.0 + exp(1.0/(x+1.0)-1) )'
  []
[]

[AuxVariables]
  [correct_u]
  []
[]

[AuxKernels]
  # This kernel is just being used to plot the correct result for u.
  [copy_from_function]
    type = FunctionAux
    variable = correct_u
    function = concentration_profile
    execute_on = 'initial timestep_begin'
  []
[]

[Executioner]
  type = Steady

  nl_abs_tol = 1e-50
  nl_rel_tol = 1e-12
  l_tol = 1e-10
[]

[Postprocessors]
  [u]
    type = ElementAverageValue
    variable = u
    execute_on = 'initial timestep_end'
  []
  [temp]
    type = ElementAverageValue
    variable = temp
    execute_on = 'initial timestep_end'
  []
  [correct_u]
    type = ElementAverageValue
    variable = correct_u
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/zirconium_diffusion/ad_test.i)

# This is an example of an axial slice from the T179 EBR-II X419 experiment,
# which had 61U-16Pu-23Zr atom%, 71U-19Pu-10Zr wt%.
#
# Most importantly, this test shows that the new calculation method for zirconium
# redistribution matches the old method, to a maximum difference of 3e-3.

[Mesh]
  coord_type = RZ
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 104
    ny = 1
    xmax = 0.00216
    ymax = 2.16e-05
  []
[]

[Variables]
  [T]
  []
  [zirconium]
    initial_condition = 0.22 #initial mole fraction of Zr in the system
    scaling = 1e17
  []
[]

[ICs]
  [T]
    type = FunctionIC
    variable = T
    function = outer2
  []
[]

[Kernels]
  [heat]
    type = ADHeatConduction
    variable = T
  []
  [ie]
    type = ADHeatConductionTimeDerivative
    variable = T
  []
  [ZrDiffusion]
    type = ADZirconiumDiffusion
    variable = zirconium
    temperature = T
  []
  [diffusion_ie]
    type = ADTimeDerivative
    variable = zirconium
  []
  [heat_source]
    type = ADCoupledForce
    variable = T
    v = power
  []
[]

[AuxVariables]
  [fission_rate]
  []
  [power]
  []
[]

[AuxKernels]
  [fission_rate_aux]
    type = ConstantAux
    variable = fission_rate
    value = 8.74117559524e+19
  []
  [power_aux]
    type = ParsedAux
    coupled_variables = 'zirconium fission_rate'
    constant_names = 'A B epf'
    constant_expressions = '-1.12884007 1.24834482 3.2e-11'
    expression = 'fission_rate * epf * (A*zirconium + B)'
    variable = power
  []
[]

[Functions]
  [outer2]
    type = ParsedFunction
    expression = '918.15-3.2e7*(x*x)' #r^2 initial temperature profile
  []
[]

[BCs]
  [temp_rightside]
    type = ADDirichletBC
    variable = T
    boundary = right
    value = 818.15
  []
[]

[Materials]
  [porosity]
    type = ADGenericConstantMaterial
    block = 0
    prop_names = porosity
    prop_values = 0.0
  []
  [phase]
    type = ADPhaseUPuZr
    block = 0
    temperature = T
    X_Pu = 0.16
    X_Zr = zirconium
    AB_temp = 936.15
    CD_temp = 968.15
  []
  [zr_diff]
    type = ADZrDiffusivityUPuZr
    block = 0
    temperature = T
    X_Pu = 0.16
    X_Zr = zirconium
    p_alpha = ${fparse 0.35 / 100}
    p_delta = ${fparse 0.35 / 2}
    p_beta = ${fparse 0.35 / 55}
    p_gamma = ${fparse 0.35 / 2}
    D0_scale_alpha = 100.0
    D0_scale_delta = 2.0
    D0_scale_beta = 55.0
    D0_scale_gamma = 2.0
  []
  [thermal_upuzr]
    type = ADUPuZrThermal
    block = 0
    temperature = T
    X_Pu = 0.16
    X_Zr = zirconium
    thcond_model = galloway
    spheat_model = savage
    A_Pu = .244
    porosity = porosity
  []
  [fuel_density]
    type = ADParsedMaterial
    block = 0
    property_name = density
    expression = 10431.0
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-ksp_gmres_restart -pc_type    -sub_pc_type -pc_hypre_type'
  petsc_options_value = ' 1001               hypre         hypre        boomeramg'
  line_search = 'none'

  l_tol = 8e-3
  nl_max_its = 15
  nl_rel_tol = 1e-9
  nl_abs_tol = 1e-18

  l_max_its = 1000
  dtmax = 1e5
  dtmin = 1e-2
  end_time = 1e5

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1.0e3
    optimal_iterations = 25
    iteration_window = 3
    linear_iteration_ratio = 35
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/zrdiffusivity_upuzr/exact.i)

# Calculates Fickian and Soret Diffusivty by changing temperature and zirconium atom fraction
# Check associated excel sheet for exact solution

D0_alpha = 1e-6
D0_beta = 1e-5
D0_delta = 2e-06
Q_alpha = 160000
Q_beta = 170000
Q_delta = 140000
Qstar_alpha = 190000
Qstar_beta = 400000
Qstar_delta = 150000
Qstar_gamma = -150000
D0_scale_alpha = 2
D0_scale_beta = 3
D0_scale_delta = 4
D0_scale_gamma = 5
Q_scale_alpha = 0.8
Q_scale_beta = 0.9
Q_scale_delta = 1.1
Q_scale_gamma = 1.2
S_scale_alpha = 6
S_scale_beta = 7
S_scale_delta = 8
S_scale_gamma = 9
p_beta = 0.05
p_alpha = 0.1
p_delta = 0.15
p_gamma = 0.25
R = 8.31446
X_Pu = 0.16
H_alpha = 150000
H_beta = -55000
H_gamma = -150000
H_delta = -3500

[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 1
  []
[]

[Problem]
  solve = false
[]

[AuxVariables]
  [temp]
    initial_condition = 773
  []
  [X_Zr]
    initial_condition = 0.01
  []
[]

[AuxKernels]
  [temp]
    type = FunctionAux
    variable = temp
    function = temp_ramp
  []
  [zr]
    type = FunctionAux
    variable = X_Zr
    function = conc_ramp
  []
[]

[Functions]
  [temp_ramp]
    type = PiecewiseLinear
    x = '0   50  100 150 200'
    y = '950 850 850 950 850'
  []
  [conc_ramp]
    type = PiecewiseLinear
    x = '0    50   100  150  200'
    y = '0.01 0.01 0.1  0.55 0.8'
  []
[]

[AuxVariables]
  [alpha_solvus]
    order = CONSTANT
    family = MONOMIAL
  []
  [beta_solvus]
    order = CONSTANT
    family = MONOMIAL
  []
  [gamma_solvus]
    order = CONSTANT
    family = MONOMIAL
  []
  [delta_solvus]
    order = CONSTANT
    family = MONOMIAL
  []
  [alpha_2phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [beta_2phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [gamma_2phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [delta_2phase]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [alpha_solvus]
    type = MaterialStdVectorAux
    variable = alpha_solvus
    property = phase_solvus
    index = 0
  []
  [beta_solvus]
    type = MaterialStdVectorAux
    variable = beta_solvus
    property = phase_solvus
    index = 1
  []
  [gamma_solvus]
    type = MaterialStdVectorAux
    variable = gamma_solvus
    property = phase_solvus
    index = 2
  []
  [delta_solvus]
    type = MaterialStdVectorAux
    variable = delta_solvus
    property = phase_solvus
    index = 3
  []
  [alpha_2phase]
    type = MaterialStdVectorAux
    variable = alpha_2phase
    property = twophase_fractions
    index = 0
  []
  [beta_2phase]
    type = MaterialStdVectorAux
    variable = beta_2phase
    property = twophase_fractions
    index = 1
  []
  [gamma_2phase]
    type = MaterialStdVectorAux
    variable = gamma_2phase
    property = twophase_fractions
    index = 2
  []
  [delta_2phase]
    type = MaterialStdVectorAux
    variable = delta_2phase
    property = twophase_fractions
    index = 3
  []
[]

[Materials]
  [phase]
    type = PhaseUPuZr
    temperature = temp
    X_Pu = ${X_Pu}
    X_Zr = X_Zr
    outputs = 'all'
    output_properties = 'alpha beta gamma delta phase phase_region'
    lever_weight_coeffs = '.1 .9' # different for testing, usually '0 1'
    calc_H = false
    H_values = '${H_alpha} ${H_beta} ${H_gamma} ${H_delta}'
  []
  [zr_diff]
    type = ZrDiffusivityUPuZr
    temperature = temp
    X_Pu = ${X_Pu}
    X_Zr = X_Zr
    D0_alpha = ${D0_alpha}
    D0_beta = ${D0_beta}
    D0_delta = ${D0_delta}
    Q_alpha = ${Q_alpha}
    Q_beta = ${Q_beta}
    Q_delta = ${Q_delta}
    Qstar_alpha = ${Qstar_alpha}
    Qstar_beta = ${Qstar_beta}
    Qstar_delta = ${Qstar_delta}
    Qstar_gamma = ${Qstar_gamma}
    D0_scale_alpha = ${D0_scale_alpha}
    D0_scale_beta = ${D0_scale_beta}
    D0_scale_delta = ${D0_scale_delta}
    D0_scale_gamma = ${D0_scale_gamma}
    Q_scale_alpha = ${Q_scale_alpha}
    Q_scale_beta = ${Q_scale_beta}
    Q_scale_delta = ${Q_scale_delta}
    Q_scale_gamma = ${Q_scale_gamma}
    S_scale_alpha = ${S_scale_alpha}
    S_scale_beta = ${S_scale_beta}
    S_scale_delta = ${S_scale_delta}
    S_scale_gamma = ${S_scale_gamma}
    p_beta = ${p_beta}
    p_alpha = ${p_alpha}
    p_delta = ${p_delta}
    p_gamma = ${p_gamma}
    limit_twophase_soret = true
    outputs = 'all'
  []
[]

[Executioner]
  type = Transient
  end_time = 200
[]

[Postprocessors]
  [temp]
    type = ElementAverageValue
    variable = temp
    execute_on = 'initial timestep_end'
  []
  [X_Zr]
    type = ElementAverageValue
    variable = X_Zr
  []

  [phase_region]
    type = ElementAverageValue
    variable = phase_region
  []

  [alpha]
    type = ElementAverageValue
    variable = alpha
  []
  [beta]
    type = ElementAverageValue
    variable = beta
  []
  [gamma]
    type = ElementAverageValue
    variable = gamma
  []
  [delta]
    type = ElementAverageValue
    variable = delta
  []

  [alpha_2phase]
    type = ElementAverageValue
    variable = alpha_2phase
  []
  [beta_2phase]
    type = ElementAverageValue
    variable = beta_2phase
  []
  [gamma_2phase]
    type = ElementAverageValue
    variable = gamma_2phase
  []
  [delta_2phase]
    type = ElementAverageValue
    variable = delta_2phase
  []

  [alpha_solvus]
    type = ElementAverageValue
    variable = alpha_solvus
  []
  [beta_solvus]
    type = ElementAverageValue
    variable = beta_solvus
  []
  [gamma_solvus]
    type = ElementAverageValue
    variable = gamma_solvus
  []
  [delta_solvus]
    type = ElementAverageValue
    variable = delta_solvus
  []
  [D_calc]
    type = ElementAverageValue
    variable = Zr_diffusivity
  []
  [S_calc]
    type = ElementAverageValue
    variable = Zr_Soret
  []

  [D_alpha]
    type = ParsedPostprocessor
    pp_names = temp
    expression = '${D0_scale_alpha} * ${D0_alpha} * exp(-${Q_alpha} * ${Q_scale_alpha} / ${R} / temp)'
    outputs = none
  []
  [D_beta]
    type = ParsedPostprocessor
    pp_names = temp
    expression = '${D0_scale_beta} * ${D0_beta} * exp(-${Q_beta} * ${Q_scale_beta} / ${R} / temp)'
    outputs = none
  []
  [D_delta]
    type = ParsedPostprocessor
    pp_names = temp
    expression = '${D0_scale_delta} * ${D0_delta} * exp(-${Q_delta} * ${Q_scale_delta} / ${R} / temp)'
    outputs = none
  []
  [D_gamma]
    type = ParsedPostprocessor
    pp_names = 'temp X_Zr'
    expression = '${D0_scale_gamma} * exp(-11.7432 - 18.5358 * X_Zr + 21.0226 * X_Zr^2) * exp(-(1.28e5 - 1.07e5 * X_Zr + 1.74e5 * X_Zr^2) * ${Q_scale_gamma} / ${R} / temp)'
    outputs = none
  []
  [D_gamma_solvus]
    type = ParsedPostprocessor
    pp_names = 'temp gamma_solvus'
    expression = '${D0_scale_gamma} * exp(-11.7432 - 18.5358 * gamma_solvus + 21.0226 * gamma_solvus^2) * exp(-(1.28e5 - 1.07e5 * gamma_solvus + 1.74e5 * gamma_solvus^2) * ${Q_scale_gamma} / ${R} / temp)'
    outputs = none
  []

  [D_exact]
    type = ParsedPostprocessor
    pp_names = 'D_alpha D_beta D_delta D_gamma D_gamma_solvus alpha beta gamma delta alpha_2phase beta_2phase gamma_2phase delta_2phase'
    expression = 'D_alpha * ((alpha - alpha_2phase) + alpha_2phase * ${p_alpha}) + D_beta * ((beta - beta_2phase) + beta_2phase * ${p_beta}) + D_delta * ((delta - delta_2phase) + delta_2phase * ${p_delta}) + D_gamma * (gamma - gamma_2phase) + D_gamma_solvus * gamma_2phase * ${p_gamma}'
  []
  [D_diff]
    type = ParsedPostprocessor
    pp_names = 'D_exact D_calc'
    expression = '(D_exact - D_calc) / D_exact'
    outputs = none
  []
  [D_max_diff]
    type = TimeExtremeValue
    postprocessor = 'D_diff'
    value_type = abs_max
  []

  [S_twophase_mixing_alpha]
    type = ParsedPostprocessor
    pp_names = 'alpha_2phase'
    expression = 'x := (alpha_2phase - 0.33) / (0.5 - 0.33);
                relx := if(x < 0, 0, if(x > 1, 1, x));
                relx^3 * (relx * (relx * 6 - 15) + 10)'
    outputs = none
  []
  [S_twophase_mixing_beta]
    type = ParsedPostprocessor
    pp_names = 'beta_2phase'
    expression = 'x := (beta_2phase - 0.33) / (0.5 - 0.33);
                relx := if(x < 0, 0, if(x > 1, 1, x));
                relx^3 * (relx * (relx * 6 - 15) + 10)'
    outputs = none
  []
  [S_twophase_mixing_delta]
    type = ParsedPostprocessor
    pp_names = 'delta_2phase'
    expression = 'x := (delta_2phase - 0.33) / (0.5 - 0.33);
                relx := if(x < 0, 0, if(x > 1, 1, x));
                relx^3 * (relx * (relx * 6 - 15) + 10)'
    outputs = none
  []
  [S_twophase_mixing_gamma]
    type = ParsedPostprocessor
    pp_names = 'gamma_2phase'
    expression = 'x := (gamma_2phase - 0.33) / (0.5 - 0.33);
                relx := if(x < 0, 0, if(x > 1, 1, x));
                relx^3 * (relx * (relx * 6 - 15) + 10)'
    outputs = none
  []

  [S_exact]
    type = ParsedPostprocessor
    pp_names = 'temp X_Zr D_alpha D_beta D_delta D_gamma D_gamma_solvus alpha beta gamma delta alpha_2phase beta_2phase gamma_2phase delta_2phase alpha_solvus beta_solvus gamma_solvus delta_solvus S_twophase_mixing_alpha S_twophase_mixing_beta S_twophase_mixing_delta S_twophase_mixing_gamma'
    expression = 'D_Qstar := D_alpha * (alpha - alpha_2phase) * ${Qstar_alpha} + D_beta * (beta - beta_2phase) * ${Qstar_beta} + D_delta * (delta - delta_2phase) * ${Qstar_delta} + D_gamma * (gamma - gamma_2phase) * ${Qstar_gamma};
                S_single_phase := D_Qstar * X_Zr * (1.0 - X_Zr - ${X_Pu}) / (1 - ${X_Pu}) / ${R} / temp / temp;
                S_two_phase_alpha := alpha_2phase * D_alpha * alpha_solvus * (1 - ${X_Pu} - alpha_solvus) * ${S_scale_alpha} * (${H_alpha} + ${Qstar_alpha}) / (1 - ${X_Pu}) / ${R} / temp / temp;
                S_two_phase_beta := beta_2phase * D_beta * beta_solvus * (1 - ${X_Pu} - beta_solvus) * ${S_scale_beta} * (${H_beta} + ${Qstar_beta}) / (1 - ${X_Pu}) / ${R} / temp / temp;
                S_two_phase_delta := delta_2phase * D_delta * delta_solvus * (1 - ${X_Pu} - delta_solvus) * ${S_scale_delta} * (${H_delta} + ${Qstar_delta}) / (1 - ${X_Pu}) / ${R} / temp / temp;
                S_two_phase_gamma := gamma_2phase * D_gamma_solvus * gamma_solvus * (1 - ${X_Pu} - gamma_solvus) * ${S_scale_gamma} * (${H_gamma} + ${Qstar_gamma}) / (1 - ${X_Pu}) / ${R} / temp / temp;
                S_single_phase + S_two_phase_alpha * S_twophase_mixing_alpha + S_two_phase_beta * S_twophase_mixing_beta + S_two_phase_delta * S_twophase_mixing_delta + S_two_phase_gamma * S_twophase_mixing_gamma'
  []
  [S_diff]
    type = ParsedPostprocessor
    pp_names = 'S_exact S_calc'
    expression = '(S_exact - S_calc) / S_exact'
    outputs = none
  []
  [S_max_diff]
    type = TimeExtremeValue
    postprocessor = 'S_diff'
    value_type = abs_max
  []

[]

[Outputs]
  csv = true
[]

(test/tests/zirconium_diffusion/test.i)

# This is an example of an axial slice from the T179 EBR-II X419 experiment,
# which had 61U-16Pu-23Zr atom%, 71U-19Pu-10Zr wt%.
#
# Most importantly, this test shows that the new calculation method for zirconium
# redistribution matches the old method, to a maximum difference of 3e-3.

[Mesh]
  coord_type = RZ
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 104
    ny = 1
    xmax = 0.00216
    ymax = 2.16e-05
  []
[]

[Variables]
  [T]
  []
  [zirconium]
    initial_condition = 0.22 #initial mole fraction of Zr in the system
    scaling = 1e17
  []
[]

[ICs]
  [T]
    type = FunctionIC
    variable = T
    function = outer2
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = T
  []
  [ie]
    type = HeatConductionTimeDerivative
    variable = T
  []
  [ZrDiffusion]
    type = ZirconiumDiffusion
    variable = zirconium
    temperature = T
  []
  [diffusion_ie]
    type = TimeDerivative
    variable = zirconium
  []
  [heat_source]
    type = CoupledForce
    variable = T
    v = power
  []
[]

[AuxVariables]
  [fission_rate]
  []
  [power]
  []
[]

[AuxKernels]
  [fission_rate_aux]
    type = ConstantAux
    variable = fission_rate
    value = 8.74117559524e+19
  []
  [power_aux]
    type = ParsedAux
    coupled_variables = 'zirconium fission_rate'
    constant_names = 'A B epf'
    constant_expressions = '-1.12884007 1.24834482 3.2e-11'
    expression = 'fission_rate * epf * (A*zirconium + B)'
    variable = power
  []
[]

[Functions]
  [outer2]
    type = ParsedFunction
    expression = '918.15-3.2e7*(x*x)' #r^2 initial temperature profile
  []
[]

[BCs]
  [temp_rightside]
    type = DirichletBC
    variable = T
    boundary = right
    value = 818.15
  []
[]

[Materials]
  [porosity]
    type = GenericConstantMaterial
    block = 0
    prop_names = porosity
    prop_values = 0.0
  []
  [phase]
    type = PhaseUPuZr
    block = 0
    temperature = T
    X_Pu = 0.16
    X_Zr = zirconium
    AB_temp = 936.15
    CD_temp = 968.15
  []
  [zr_diff]
    type = ZrDiffusivityUPuZr
    block = 0
    temperature = T
    X_Pu = 0.16
    X_Zr = zirconium
    p_alpha = ${fparse 0.35 / 100}
    p_delta = ${fparse 0.35 / 2}
    p_beta = ${fparse 0.35 / 55}
    p_gamma = ${fparse 0.35 / 2}
    D0_scale_alpha = 100.0
    D0_scale_delta = 2.0
    D0_scale_beta = 55.0
    D0_scale_gamma = 2.0
  []
  [thermal_upuzr]
    type = UPuZrThermal
    block = 0
    temperature = T
    X_Pu = 0.16
    X_Zr = zirconium
    thcond_model = galloway
    spheat_model = savage
    A_Pu = .244
    porosity = porosity
  []
  [fuel_density]
    type = ParsedMaterial
    block = 0
    property_name = density
    expression = 10431.0
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-ksp_gmres_restart -pc_type    -sub_pc_type -pc_hypre_type'
  petsc_options_value = ' 1001               hypre         hypre        boomeramg'
  line_search = 'none'

  l_tol = 8e-3
  nl_max_its = 15
  nl_rel_tol = 1e-9
  nl_abs_tol = 1e-18

  l_max_its = 1000
  dtmax = 1e5
  dtmin = 1e-2
  end_time = 1e5

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1.0e3
    optimal_iterations = 25
    iteration_window = 3
    linear_iteration_ratio = 35
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/zrdiffusivity_upuzr/patch.i)

# Calculates Fickian and Soret Diffusivty by changing temperature and zirconium atom fraction

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

[Variables]
  [temp]
    initial_condition = 773
  []
  [X_Zr]
    initial_condition = 0.01
  []
[]

[Kernels]
  [temp_diff]
    type = Diffusion
    variable = temp
  []
  [Zr_diff]
    type = Diffusion
    variable = X_Zr
  []
[]

[AuxVariables]
  [alpha_solvus]
    order = CONSTANT
    family = MONOMIAL
  []
  [beta_solvus]
    order = CONSTANT
    family = MONOMIAL
  []
  [gamma_solvus]
    order = CONSTANT
    family = MONOMIAL
  []
  [delta_solvus]
    order = CONSTANT
    family = MONOMIAL
  []
  [alpha_2phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [beta_2phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [gamma_2phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [delta_2phase]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [alpha_solvus]
    type = MaterialStdVectorAux
    variable = alpha_solvus
    property = phase_solvus
    index = 0
  []
  [beta_solvus]
    type = MaterialStdVectorAux
    variable = beta_solvus
    property = phase_solvus
    index = 1
  []
  [gamma_solvus]
    type = MaterialStdVectorAux
    variable = gamma_solvus
    property = phase_solvus
    index = 2
  []
  [delta_solvus]
    type = MaterialStdVectorAux
    variable = delta_solvus
    property = phase_solvus
    index = 3
  []
  [alpha_2phase]
    type = MaterialStdVectorAux
    variable = alpha_2phase
    property = twophase_fractions
    index = 0
  []
  [beta_2phase]
    type = MaterialStdVectorAux
    variable = beta_2phase
    property = twophase_fractions
    index = 1
  []
  [gamma_2phase]
    type = MaterialStdVectorAux
    variable = gamma_2phase
    property = twophase_fractions
    index = 2
  []
  [delta_2phase]
    type = MaterialStdVectorAux
    variable = delta_2phase
    property = twophase_fractions
    index = 3
  []
[]

[BCs]
  [temp_BC]
    type = FunctionDirichletBC
    boundary = 'top bottom left right'
    variable = temp
    function = temp_ramp
  []
  [conc_BC]
    type = FunctionDirichletBC
    boundary = 'top bottom left right'
    variable = X_Zr
    function = conc_ramp
  []
[]

[Functions]
  [temp_ramp]
    type = PiecewiseLinear
    x = '0   50  100 150 200'
    y = '950 850 850 950 850'
  []
  [conc_ramp]
    type = PiecewiseLinear
    x = '0    50   100  150  200'
    y = '0.01 0.01 0.1  0.55 0.8'
  []
[]

[Materials]
  [phase]
    type = PhaseUPuZr
    block = 0
    temperature = temp
    X_Pu = 0.16
    X_Zr = X_Zr
    outputs = 'all'
    output_properties = 'alpha beta gamma delta phase phase_region'
    lever_weight_coeffs = '.1 .9' # different for testing, usually '0 1'
    calc_H = false
  []
  [zr_diff]
    type = ZrDiffusivityUPuZr
    block = 0
    temperature = temp
    X_Pu = 0.16
    X_Zr = X_Zr
    outputs = 'all'
  []
[]

[Executioner]
  type = Transient
  end_time = 200
[]

[Postprocessors]
  [temp]
    type = ElementAverageValue
    variable = temp
    execute_on = 'initial timestep_end'
  []
  [Zr]
    type = ElementAverageValue
    variable = X_Zr
  []

  [phase_region]
    type = ElementAverageValue
    variable = phase_region
  []

  [alpha]
    type = ElementAverageValue
    variable = alpha
  []
  [beta]
    type = ElementAverageValue
    variable = beta
  []
  [gamma]
    type = ElementAverageValue
    variable = gamma
  []
  [delta]
    type = ElementAverageValue
    variable = delta
  []

  [alpha_2phase]
    type = ElementAverageValue
    variable = alpha_2phase
  []
  [beta_2phase]
    type = ElementAverageValue
    variable = beta_2phase
  []
  [gamma_2phase]
    type = ElementAverageValue
    variable = gamma_2phase
  []
  [delta_2phase]
    type = ElementAverageValue
    variable = delta_2phase
  []

  [alpha_solvus]
    type = ElementAverageValue
    variable = alpha_solvus
  []
  [beta_solvus]
    type = ElementAverageValue
    variable = beta_solvus
  []
  [gamma_solvus]
    type = ElementAverageValue
    variable = gamma_solvus
  []
  [delta_solvus]
    type = ElementAverageValue
    variable = delta_solvus
  []
  [zr_diff]
    type = ElementAverageValue
    variable = Zr_diffusivity
  []
  [zr_soret]
    type = ElementAverageValue
    variable = Zr_Soret
  []
[]

[Outputs]
  csv = true
[]

Description
Numerical Implementation in BISON
Example Input Syntax
Input Parameters
Input Files