PKAFissionFragmentNeutronics

PKAFissionFragmentNeutronics generates primary knock-on atoms (PKA) originating from fission reactions. The difference to PKAFissionFragmentEmpirical is that it uses isotopic fission rates computed by neutronics calculations and ENDF data for sampling fission product species.

The partial fission rate for isotope $var element = document.getElementById("moose-equation-88343efb-61cf-4a13-8679-7e6726c63974");katex.render("i", element, {displayMode:false,throwOnError:false});$ is denoted by $var element = document.getElementById("moose-equation-949971b8-f1f5-4b28-b187-570a3d3c73e8");katex.render("F_i(\\vec{r}, \\vec{\\rho})", element, {displayMode:false,throwOnError:false});$ . It is computed by:

$var element = document.getElementById("moose-equation-9a79ab2f-1eb0-4785-9900-01ca5d55ffa6");katex.render("F_i(\\vec{r}, \\vec{\\rho}) = \\sum\\limits_{g=1}^G N_i(\\vec{r},\\vec{\\rho})\\sigma_{f,g,i} \\phi_g(\\vec{r}),", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-80c11971-419d-4bde-a7ae-9626033c754d");katex.render("\\vec{r}", element, {displayMode:false,throwOnError:false});$ is the location in the macroscopic, neutronics domain and $var element = document.getElementById("moose-equation-3f3d68c9-4c9d-4094-9cc2-38c92321c28d");katex.render("\\vec{\\rho}", element, {displayMode:false,throwOnError:false});$ is the location in the microscopic domain. These two locations are separated because their scale is significantly separated by 3-4 orders of magnitude. Changes with $var element = document.getElementById("moose-equation-60253a8a-22a9-46e7-891e-3af2fb5a5158");katex.render("\\vec{r}", element, {displayMode:false,throwOnError:false});$ are understood to be smooth changes of the average composition of the microscopic domain (orders of centimeters), while changes with $var element = document.getElementById("moose-equation-e15e6e33-820f-420d-9c60-dfa16613570d");katex.render("\\vec{\\rho}", element, {displayMode:false,throwOnError:false});$ captures compositional changes between grains on the microscopic domain, i.e. on the orders of micro-meters. The notation of separating the spatial dependence into two scales follows homogenization theory. $var element = document.getElementById("moose-equation-eb96adf9-3065-4ebf-b4a7-8fe1b75ae356");katex.render("N_i", element, {displayMode:false,throwOnError:false});$ is the number density of isotope $var element = document.getElementById("moose-equation-1cea6cec-e46e-4284-a923-99dc8c0a7c75");katex.render("i", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-988d8062-2d36-45e4-aeb3-b74c91656b60");katex.render("\\sigma_{f,g,i}", element, {displayMode:false,throwOnError:false});$ is the microscopic fission cross section of isotope $var element = document.getElementById("moose-equation-aba60326-bb47-46e2-90d0-7c48ba52af85");katex.render("i", element, {displayMode:false,throwOnError:false});$ in group $var element = document.getElementById("moose-equation-d773be46-9969-431d-a164-10e7be6d27d5");katex.render("g", element, {displayMode:false,throwOnError:false});$ , and $var element = document.getElementById("moose-equation-571bc1fc-8dd1-4e05-8d0b-3329ad87f946");katex.render("\\phi_g", element, {displayMode:false,throwOnError:false});$ is the scalar flux in group $var element = document.getElementById("moose-equation-f07602a6-5da2-4d47-8ca9-d37c6c922f3f");katex.render("g", element, {displayMode:false,throwOnError:false});$ .

Denoting the microscopic domain as $var element = document.getElementById("moose-equation-d5a33064-044e-4cc8-b3ad-525897f87bcc");katex.render("\\Theta_m", element, {displayMode:false,throwOnError:false});$ located at $var element = document.getElementById("moose-equation-fc21b5b7-c662-4682-a92d-9f8f825c0501");katex.render("\\vec{r}_m", element, {displayMode:false,throwOnError:false});$ , we can define a slowly varying average of the number densities:

$var element = document.getElementById("moose-equation-fddf708a-e233-4dfb-99b9-74a1e9b170c1");katex.render(" N_i(\\vec{r}) = \\frac{1}{\\Theta_m} \\int_{\\Theta_m} N_i(\\vec{r},\\vec{\\rho}) d\\vec{\\rho}.", element, {displayMode:true,throwOnError:false});$

$var element = document.getElementById("moose-equation-a5f717dc-6389-4c65-9076-41a98bf037a4");katex.render("N_i(\\vec{r})", element, {displayMode:false,throwOnError:false});$ is the number density provided to neutronics calculations. The nuclide fission rates computed by the neutronics calculation is consequently the slowly varying average given by:

$var element = document.getElementById("moose-equation-07c1d9e0-5ac4-4b52-a374-5db92cf932b6");katex.render("F_i (\\vec{r}) = \\sum\\limits_{g=1}^G N_i(\\vec{r})\\sigma_{f,g,i}.", element, {displayMode:true,throwOnError:false});$

The PKAFissionFragmentNeutronics accepts the values of $var element = document.getElementById("moose-equation-47c2d596-4e38-4c83-9e31-dd24b4ce1d14");katex.render("F_i (\\vec{r}_m)", element, {displayMode:false,throwOnError:false});$ as the partial_reaction_rates parameter. The fission rate density in the microscopic domain is computed by:

$var element = document.getElementById("moose-equation-27ea10d0-9a27-4f5d-850b-e11bb9a0a470");katex.render(" F_i (\\vec{r}_m, \\vec{\\rho}) = \\frac{N_i(\\vec{r}_m, \\vec{\\rho})}{N_i(\\vec{r}_m)} F_i(\\vec{r}_m),", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-0dcd627d-79c0-4f28-8fd1-bc3684ea220b");katex.render("N_i(\\vec{r}_m, \\vec{\\rho}) = \\gamma_i(\\vec{\\rho}) / \\Omega(\\vec{\\rho})", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-5a1ace9a-d820-4b86-a33f-3356544cc974");katex.render("\\Omega(\\vec{\\rho})", element, {displayMode:false,throwOnError:false});$ is the site volume that is provided to the MyTRIMRasterizer and $var element = document.getElementById("moose-equation-d76811ec-176c-4094-ab36-65d2803fc805");katex.render("\\gamma_i", element, {displayMode:false,throwOnError:false});$ is the $var element = document.getElementById("moose-equation-94034c00-d4e7-45ef-8ee5-2ecb388ee7b9");katex.render("i", element, {displayMode:false,throwOnError:false});$ -th variable provided to the MyTRIMRasterizer's var parameter. Hence, we have:

$var element = document.getElementById("moose-equation-637c2ace-32d5-4523-945e-922e69f58992");katex.render(" F_i (\\vec{r}_m, \\vec{\\rho}) = \\frac{\\gamma_i(\\vec{\\rho})}{\\Omega(\\vec{r}) N_i(\\vec{r}_m)} F_i(\\vec{r}_m),", element, {displayMode:true,throwOnError:false});$

PKAFissionFragmentNeutronics accepts the values of $var element = document.getElementById("moose-equation-610c795f-8e7b-417a-a368-d4da05a3ac7b");katex.render("N_i(\\vec{r}_m)", element, {displayMode:false,throwOnError:false});$ in the averaged_number_densities parameter.

The expected number of fissions in a mesh element with index $var element = document.getElementById("moose-equation-fff8d08f-1fa9-42ee-9860-b5d2723bf436");katex.render("j", element, {displayMode:false,throwOnError:false});$ and volume $var element = document.getElementById("moose-equation-2551f068-857f-489a-ab83-f6c6cff18cbd");katex.render("V_j", element, {displayMode:false,throwOnError:false});$ is given by:

$var element = document.getElementById("moose-equation-d6319202-5704-4b83-840b-da96bab3ddec");katex.render("C_i = \\Delta t \\int_{V_j} F_i (\\vec{r}_m, \\vec{\\rho}) d\\vec{\\rho}.", element, {displayMode:true,throwOnError:false});$

Non-integer results are rounded up with a probability of the $var element = document.getElementById("moose-equation-9d9564ad-f80a-4a42-a5f0-2584325757e8");katex.render("C_i - \\text{int}(C_i)", element, {displayMode:false,throwOnError:false});$ ; otherwise rounded down. Two PKAs are created from each fission event. The algorithm for sampling their type, energy, and direction of motion is described in Schunert and others (2017).

Input Parameters

averaged_number_densitiesThe number density of the species averaged over the domain.
C++ Type:std::vector<PostprocessorName>
Unit:(no unit assumed)
Controllable:No
Description:The number density of the species averaged over the domain.
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
partial_reaction_ratesPartial neutronic reaction rates per unit volume [sum_g xs_{r,g,i} * phi_g], r: reaction type, g: energy group, i: nuclide id.Provide number density as variable in rasterizer!
C++ Type:std::vector<PostprocessorName>
Unit:(no unit assumed)
Controllable:No
Description:Partial neutronic reaction rates per unit volume [sum_g xs_{r,g,i} * phi_g], r: reaction type, g: energy group, i: nuclide id.Provide number density as variable in rasterizer!

Optional Parameters

allow_duplicate_execution_on_initialFalseIn the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).
Default:False
C++ Type:bool
Controllable:No
Description:In the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).
execution_order_group0Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.
Default:0
C++ Type:int
Controllable:No
Description:Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.
force_postauxFalseForces the UserObject to be executed in POSTAUX
Default:False
C++ Type:bool
Controllable:No
Description:Forces the UserObject to be executed in POSTAUX
force_preauxFalseForces the UserObject to be executed in PREAUX
Default:False
C++ Type:bool
Controllable:No
Description:Forces the UserObject to be executed in PREAUX
force_preicFalseForces the UserObject to be executed in PREIC during initial setup
Default:False
C++ Type:bool
Controllable:No
Description:Forces the UserObject to be executed in PREIC during initial setup

Execution Scheduling Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
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

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

(tests/radiation_damage/coupled_fission_mockup/damage_sub.i)
(tests/userobjects/mytrim/UO2_pellet_initial.i)
(tests/radiation_damage/coupled_fission_mockup/damage_sub_3D.i)

References

S. Schunert and others. Heat Source Characterization In A TREAT Fuel Particle Using Coupled Neutronics Binary Collision Monte-Carlo Calculations. In M&C 2017 International Conference on Mathematics & Computational Methods Applied to Nuclear Science & Engineering. April 2017.

@inproceedings{SchunertTREATHeatSource,
    author = "Schunert, S. and others",
    title = "{Heat Source Characterization In A TREAT Fuel Particle Using Coupled Neutronics Binary Collision Monte-Carlo Calculations}",
    booktitle = "{M\\&C 2017 International Conference on Mathematics \\& Computational Methods Applied to Nuclear Science \\& Engineering}",
    month = "April",
    year = "2017"
}

(tests/radiation_damage/coupled_fission_mockup/damage_sub.i)

# dimensions are in Angstrom (A = 1e-10 m)
# typical treat fuel particle is rf = 20 mu-m = 2e5 A
# typical distance between fuel particles is D = 0.893 * rf / cube_root(1/2571)
# D = 244 mu-m = 2.44e6 A.
# the typical fuel particle distance is much larger than the mfp so we cut the domain
# at 4 * rf = 80 mu-m = 8e5 A
#
# _NOTE_ this calculation is actually for 3D but this is a 2D setup to save
# runtime
[Mesh]
  type = MyTRIMMesh
  dim = 2
  nx = 200
  ny = 200
  xmin = -8.0e5
  xmax = 8.0e5
  ymin = -8.0e5
  ymax = 8.0e5
  elem_type = QUAD4
  uniform_refine = 0
[]

[Problem]
  kernel_coverage_check = false
[]

[Variables]
  [c]
    [InitialCondition]
      type = ConstantIC
      value = 1
    []
  []
[]

[AuxVariables]
  [int_U]
    order = CONSTANT
    family = MONOMIAL
  []
  [vac_U]
    order = CONSTANT
    family = MONOMIAL
  []
  [int_O]
    order = CONSTANT
    family = MONOMIAL
  []
  [vac_O]
    order = CONSTANT
    family = MONOMIAL
  []
  [int_C]
    order = CONSTANT
    family = MONOMIAL
  []
  [vac_C]
    order = CONSTANT
    family = MONOMIAL
  []
  [c_U]
    order = CONSTANT
    family = MONOMIAL
    [InitialCondition]
      type = FunctionIC
      function = func_c_U
    []
  []
  [c_O]
    order = CONSTANT
    family = MONOMIAL
    [InitialCondition]
      type = FunctionIC
      function = func_c_O
    []
  []
  [c_C]
    order = CONSTANT
    family = MONOMIAL
    [InitialCondition]
      type = FunctionIC
      function = func_c_C
    []
  []
  [pka_count]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [int_U_aux]
    variable = int_U
    type = MyTRIMElementResultAux
    runner = runner
    ivar = 0
    defect = INT
  []
  [vac_U_aux]
    variable = vac_U
    type = MyTRIMElementResultAux
    runner = runner
    ivar = 0
    defect = VAC
  []
  [int_O_aux]
    variable = int_O
    type = MyTRIMElementResultAux
    runner = runner
    ivar = 1
    defect = INT
  []
  [vac_O_aux]
    variable = vac_O
    type = MyTRIMElementResultAux
    runner = runner
    ivar = 1
    defect = VAC
  []
  [int_C_aux]
    variable = int_C
    type = MyTRIMElementResultAux
    runner = runner
    ivar = 2
    defect = INT
  []
  [vac_C_aux]
    variable = vac_C
    type = MyTRIMElementResultAux
    runner = runner
    ivar = 2
    defect = VAC
  []
[]

[UserObjects]
  [neutronics_fission_generator]
    type = PKAFissionFragmentNeutronics
    partial_reaction_rates = '2.0e-11 0 0'
  []
  [rasterizer]
    type = MyTRIMRasterizer
    var = 'c_U  c_O  c_C'
    M = '235  16    12'
    Z = '92   8     6'
    site_volume = 0.0404 # nm^3 per UO2 unit
    pka_generator = neutronics_fission_generator
  []
  [runner]
    type = MyTRIMElementRun
    rasterizer = rasterizer
  []
[]

[Postprocessors]
  #[./pka_info]
  #  type = MyTRIMPKAInfo
  #  rasterizer = rasterizer
  #  value_type = TOTAL_ENERGY
  #[../]
  [int_U_pp]
    type = ElementIntegralVariablePostprocessor
    variable = int_U
    execute_on = timestep_end
  []
  [vac_U_pp]
    type = ElementIntegralVariablePostprocessor
    variable = vac_U
    execute_on = timestep_end
  []
  [int_O_pp]
    type = ElementIntegralVariablePostprocessor
    variable = int_O
    execute_on = timestep_end
  []
  [vac_O_pp]
    type = ElementIntegralVariablePostprocessor
    variable = vac_O
    execute_on = timestep_end
  []
  [int_C_pp]
    type = ElementIntegralVariablePostprocessor
    variable = int_C
    execute_on = timestep_end
  []
  [vac_C_pp]
    type = ElementIntegralVariablePostprocessor
    variable = vac_C
    execute_on = timestep_end
  []
  [n_pkas]
    type = ElementIntegralVariablePostprocessor
    variable = pka_count
    execute_on = timestep_end
  []
[]

[Functions]
  [func_c_U]
    type = ParsedFunction
    expression = 'r := sqrt(x*x+y*y); rho := if(r/rf<1.0e-10,1.0e-10,r/rf); th := tanh(a*(rho-1/rho)); 0.5 * (1+th) * ev + 0.5 * (1-th) * iv'
    symbol_names = 'rf    a ev iv'
    symbol_values = '2.0e5 8  0 0.3333333'
  []
  [func_c_O]
    type = ParsedFunction
    expression = 'r := sqrt(x*x+y*y); rho := if(r/rf<1.0e-10,1.0e-10,r/rf); th := tanh(a*(rho-1/rho)); 0.5 * (1+th) * ev + 0.5 * (1-th) * iv'
    symbol_names = 'rf    a ev iv'
    symbol_values = '2.0e5 8  0 0.66666666'
  []
  [func_c_C]
    # ev computed as follows:
    # rho_G = 1.72 g/cc (see Kun MO's paper)
    # rho_UO2 = 10.963 g/cc
    # MUO2 = 235 + 2*16 = 267
    # MG = 12
    # xG = rho_G / rho_UO2 * MUO2 / MG = 3.4908
    type = ParsedFunction
    expression = 'r := sqrt(x*x+y*y); rho := if(r/rf<1.0e-10,1.0e-10,r/rf); th := tanh(a*(rho-1/rho)); 0.5 * (1+th) * ev + 0.5 * (1-th) * iv'
    symbol_names = 'rf    a ev     iv'
    symbol_values = '2.0e5 8 3.4908 0'
  []
[]

[Executioner]
  type = Transient
  num_steps = 1
  nl_abs_tol = 1e-10
[]

[Outputs]
  exodus = true
  csv = true
[]

(tests/userobjects/mytrim/UO2_pellet_initial.i)

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 60
  ny = 60
  xmin = -100
  xmax = 100
  ymin = -100
  ymax = 100
  elem_type = QUAD4
  uniform_refine = 0
[]

[Variables]
  [c]
    initial_condition = 1.0
  []
[]

[AuxVariables]
  [int]
    order = CONSTANT
    family = MONOMIAL
  []
  [vac]
    order = CONSTANT
    family = MONOMIAL
  []
  [fuel_conc]
    order = CONSTANT
    family = MONOMIAL
    [InitialCondition]
      type = FunctionIC
      function = U02_conc
    []
  []
[]

[Kernels]
  [dt]
    type = TimeDerivative
    variable = c
  []
  [diff]
    type = Diffusion
    variable = c
  []
[]

[AuxKernels]
  [int]
    variable = int
    type = MyTRIMElementResultAux
    runner = runner
    ivar = 0
    defect = INT
  []
  [vac]
    variable = vac
    type = MyTRIMElementResultAux
    runner = runner
    ivar = 0
    defect = VAC
  []
[]

#[BCs]
#  [./Periodic]
#    [./all]
#      auto_direction = 'x y'
#    [../]
#  [../]
#[]

[UserObjects]
  [thermal_fission]
    type = PKAFissionFragmentNeutronics
  []
  [rasterizer]
    type = MyTRIMRasterizer
    var = c
    M = 40
    Z = 20
    site_volume = 0.0404 # nm^3 per UO2 unit
    pka_generator = thermal_fission
  []
  [runner]
    type = MyTRIMElementRun
    rasterizer = rasterizer
  []
[]

[Postprocessors]
  [int]
    type = ElementIntegralVariablePostprocessor
    variable = int
    execute_on = timestep_end
  []
  [vac]
    type = ElementIntegralVariablePostprocessor
    variable = vac
    execute_on = timestep_end
  []
[]

[Functions]
  [U02_conc]
    type = ParsedFunction
    expression = 'r:=sqrt(x*x+y*y); Rlow:=low*R; Rhigh:=high*R; Z:= (r-high*R)/(r-low*R); if(r<Rlow,1,
                                                                                     if(r>Rhigh,0,
                                                                                     -(exp(2*Z)-1)/(exp(2*Z)+1)))'
    #-(1/(high-low))*(r/R-1)+0.5))'
    symbol_names = 'R low high'
    symbol_values = '20 0.9 1.1'
  []
[]

[Executioner]
  type = Transient
  num_steps = 1
  nl_abs_tol = 1e-10
[]

[Outputs]
  exodus = true
  csv = true
[]

(tests/radiation_damage/coupled_fission_mockup/damage_sub_3D.i)

# dimensions are in Angstrom (A = 1e-10 m)
# typical treat fuel particle is rf = 20 mu-m = 2e5 A
# typical distance between fuel particles is D = 0.893 * rf / cube_root(1/2571)
# D = 244 mu-m = 2.44e6 A.
# the typical fuel particle distance is much larger than the mfp so we cut the domain
# at 4 * rf = 80 mu-m = 8e5 A
#
# runtime
[Mesh]
  type = MyTRIMMesh
  dim = 3
  nx = 100
  ny = 100
  nz = 100
  xmin = -4.0e5
  xmax = 4.0e5
  ymin = -4.0e5
  ymax = 4.0e5
  zmin = -4.0e5
  zmax = 4.0e5
  uniform_refine = 0
[]

[Problem]
  kernel_coverage_check = false
[]

[Variables]
  [c]
    [InitialCondition]
      type = ConstantIC
      value = 1
    []
  []
[]

[AuxVariables]
  [int_U]
    order = CONSTANT
    family = MONOMIAL
  []
  [vac_U]
    order = CONSTANT
    family = MONOMIAL
  []
  [int_O]
    order = CONSTANT
    family = MONOMIAL
  []
  [vac_O]
    order = CONSTANT
    family = MONOMIAL
  []
  [int_C]
    order = CONSTANT
    family = MONOMIAL
  []
  [vac_C]
    order = CONSTANT
    family = MONOMIAL
  []
  [c_U]
    order = CONSTANT
    family = MONOMIAL
    [InitialCondition]
      type = FunctionIC
      function = func_c_U
    []
  []
  [c_O]
    order = CONSTANT
    family = MONOMIAL
    [InitialCondition]
      type = FunctionIC
      function = func_c_O
    []
  []
  [c_C]
    order = CONSTANT
    family = MONOMIAL
    [InitialCondition]
      type = FunctionIC
      function = func_c_C
    []
  []
  [pka_count]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [int_U_aux]
    variable = int_U
    type = MyTRIMElementResultAux
    runner = runner
    ivar = 0
    defect = INT
  []
  [vac_U_aux]
    variable = vac_U
    type = MyTRIMElementResultAux
    runner = runner
    ivar = 0
    defect = VAC
  []
  [int_O_aux]
    variable = int_O
    type = MyTRIMElementResultAux
    runner = runner
    ivar = 1
    defect = INT
  []
  [vac_O_aux]
    variable = vac_O
    type = MyTRIMElementResultAux
    runner = runner
    ivar = 1
    defect = VAC
  []
  [int_C_aux]
    variable = int_C
    type = MyTRIMElementResultAux
    runner = runner
    ivar = 2
    defect = INT
  []
  [vac_C_aux]
    variable = vac_C
    type = MyTRIMElementResultAux
    runner = runner
    ivar = 2
    defect = VAC
  []
[]

[UserObjects]
  [neutronics_fission_generator]
    type = PKAFissionFragmentNeutronics
    relative_density = 1
    partial_reaction_rates = '1.0e-14 0 0'
  []
  [rasterizer]
    type = MyTRIMRasterizer
    var = 'c_U  c_O  c_C'
    M = '235  16    12'
    Z = '92   8     6'
    site_volume = 0.0404 # nm^3 per UO2 unit
    pka_generator = neutronics_fission_generator
    print_pka_statistics = true
  []
  [runner]
    type = MyTRIMElementRun
    rasterizer = rasterizer
  []
[]

[Postprocessors]
  [int_U_pp]
    type = ElementIntegralVariablePostprocessor
    variable = int_U
    execute_on = timestep_end
  []
  [vac_U_pp]
    type = ElementIntegralVariablePostprocessor
    variable = vac_U
    execute_on = timestep_end
  []
  [int_O_pp]
    type = ElementIntegralVariablePostprocessor
    variable = int_O
    execute_on = timestep_end
  []
  [vac_O_pp]
    type = ElementIntegralVariablePostprocessor
    variable = vac_O
    execute_on = timestep_end
  []
  [int_C_pp]
    type = ElementIntegralVariablePostprocessor
    variable = int_C
    execute_on = timestep_end
  []
  [vac_C_pp]
    type = ElementIntegralVariablePostprocessor
    variable = vac_C
    execute_on = timestep_end
  []
  [n_pkas]
    type = ElementIntegralVariablePostprocessor
    variable = pka_count
    execute_on = timestep_end
  []
[]

[Functions]
  [func_c_U]
    type = ParsedFunction
    expression = 'r := sqrt(x*x+y*y+z*z); rho := if(r/rf<1.0e-10,1.0e-10,r/rf); th := tanh(a*(rho-1/rho)); 0.5 * (1+th) * ev + 0.5 * (1-th) * iv'
    symbol_names = 'rf    a ev iv'
    symbol_values = '2.0e5 8  0 0.3333333'
  []
  [func_c_O]
    type = ParsedFunction
    expression = 'r := sqrt(x*x+y*y+z*z); rho := if(r/rf<1.0e-10,1.0e-10,r/rf); th := tanh(a*(rho-1/rho)); 0.5 * (1+th) * ev + 0.5 * (1-th) * iv'
    symbol_names = 'rf    a ev iv'
    symbol_values = '2.0e5 8  0 0.66666666'
  []
  [func_c_C]
    # ev computed as follows:
    # rho_G = 1.72 g/cc (see Kun MO's paper)
    # rho_UO2 = 10.963 g/cc
    # MUO2 = 235 + 2*16 = 267
    # MG = 12
    # xG = rho_G / rho_UO2 * MUO2 / MG = 3.4908
    type = ParsedFunction
    expression = 'r := sqrt(x*x+y*y+z*z); rho := if(r/rf<1.0e-10,1.0e-10,r/rf); th := tanh(a*(rho-1/rho)); 0.5 * (1+th) * ev + 0.5 * (1-th) * iv'
    symbol_names = 'rf    a ev     iv'
    symbol_values = '2.0e5 8 3.4908 0'
  []
[]

[Executioner]
  type = Transient
  num_steps = 1
  nl_abs_tol = 1e-10
[]

[Outputs]
  exodus = true
  csv = true
[]

Input Parameters
Input Files
References