Doppler Slab Analytic Benchmark

Problem Description

This test case models the 1D Doppler slab benchmark from Remley and Griesheimer (2019), which couples neutron transport and heat conduction in a semi-infinite slab with microscopic absorption cross section $var element = document.getElementById("moose-equation-d10624b6-4cb9-4847-b56c-bdba7d61f339");katex.render("\\sigma_{0}", element, {displayMode:false,throwOnError:false});$ , number density $var element = document.getElementById("moose-equation-397fadcf-f1a0-492c-8e3a-0721518d9d99");katex.render("N_{0}", element, {displayMode:false,throwOnError:false});$ , and thermal conductivity $var element = document.getElementById("moose-equation-a03b7fe5-2a83-4a72-a87f-873a5aeac514");katex.render("k", element, {displayMode:false,throwOnError:false});$ . A beam of neutrons (intensity $var element = document.getElementById("moose-equation-8abd9440-b6dd-4ba2-a235-369aec1c8d24");katex.render("\\psi_{0}", element, {displayMode:false,throwOnError:false});$ ) is incident on the left face ( $var element = document.getElementById("moose-equation-cc384b20-a266-4af9-aa34-16b6fd4f546c");katex.render("z = 0", element, {displayMode:false,throwOnError:false});$ ) of the slab, which is held at a fixed temperature of $var element = document.getElementById("moose-equation-649abbc5-09b1-4b36-9194-63962f77944a");katex.render("T_{0}", element, {displayMode:false,throwOnError:false});$ . The geometry of the benchmark problem can be found in Figure 1.

Doppler slab geometry, taken from [!cite](doppler_slab_benchmark).

Figure 1: Doppler slab geometry, taken from Remley and Griesheimer (2019).

The macroscopic absorption cross section in the slab takes one of two functional forms. The first is a linear dependence on temperature:

(1)var element = document.getElementById("moose-equation-ab6c27fb-dd88-4440-8b0c-efefcda510aa");katex.render("\\Sigma(z) = \\underbrace{N_{0}\\sigma_{0}}_{\\Sigma_{0}} + \\alpha\\left(T(z) - T_{0}\\right)\\text{.}", element, {displayMode:true,throwOnError:false});

The second is an inverse square-root dependence on temperature:

(2)var element = document.getElementById("moose-equation-f5fba032-27c0-4be0-9137-f32bc62bdd62");katex.render("\\Sigma(z) = \\underbrace{N_{0}\\sigma_{0}}_{\\Sigma_{0}}\\sqrt{\\frac{T_{0}}{T(z)}}\\text{.}", element, {displayMode:true,throwOnError:false});

Absorption reactions heat the slab with an energy released per reaction of $var element = document.getElementById("moose-equation-94b0efe2-fbd8-4123-8e25-12bc6ead7245");katex.render("q", element, {displayMode:false,throwOnError:false});$ . The analytic solution for the temperature (Eq. (3)) and flux (Eq. (4)) in the slab when the cross section is represented by Eq. (1) can be obtained:

(3)var element = document.getElementById("moose-equation-8e04131b-e9ba-4b0d-8e2b-b824e25c2c14");katex.render("T(z) = T_{0} + \\frac{2 \\phi_{0} T_{0}\\left[1 - \\exp\\left(-\\sqrt{A}\\Sigma_{0} z\\right)\\right]} {1 - \\exp\\left(-\\sqrt{A}\\Sigma_{0}z\\right) + \\sqrt{A}\\left[1 + \\exp\\left(-\\sqrt{A}\\Sigma_{0}z\\right)\\right]}\\text{,}", element, {displayMode:true,throwOnError:false});

(4)var element = document.getElementById("moose-equation-8961fc17-c44a-4d0d-a636-213e336f7235");katex.render("\\psi(z) = \\frac{4 \\psi_{0} A\\left[1 - \\exp\\left(-\\sqrt{A}\\Sigma_{0} z\\right)\\right]} {\\left[1 + \\sqrt{A} - \\left(1 - \\sqrt{A}\\right)\\exp\\left(-\\sqrt{A}\\Sigma_{0}z\\right)\\right]^{2}}\\text{.}", element, {displayMode:true,throwOnError:false});

$var element = document.getElementById("moose-equation-56bc7ef8-61f6-4e87-925c-e8b27dc0e660");katex.render("\\phi_{0}", element, {displayMode:false,throwOnError:false});$ is defined as:

(5)var element = document.getElementById("moose-equation-14f9d90e-9860-46de-803d-c23342d05197");katex.render("\\phi_{0} = \\frac{\\psi_{0} q}{k T_{0} \\Sigma_{0}}\\text{,}", element, {displayMode:true,throwOnError:false});

and $var element = document.getElementById("moose-equation-be67cf2a-7109-45fe-bc86-0ded01dc3fd8");katex.render("A", element, {displayMode:false,throwOnError:false});$ is:

(6)var element = document.getElementById("moose-equation-7d5dd48c-aadb-435d-95f0-d77eaa8f59d0");katex.render("A = 1 + \\frac{2\\alpha\\psi_{0}q}{k\\Sigma_{0}^2}\\text{.}", element, {displayMode:true,throwOnError:false});

For the case where Eq. (2) is used for the cross section, the following solutions for the temperature (Eq. (7)) and flux (Eq. (8)) can be obtained:

(7)var element = document.getElementById("moose-equation-db8d03ed-91e9-469a-912f-17b211c3297f");katex.render("T(z) = T_{0}B^2\\left[W\\left(\\frac{\\phi_{0}}{2B}\\exp\\left(\\frac{\\Sigma_{0}z + C}{B}\\right)\\right) + 1\\right]^2\\text{,}", element, {displayMode:true,throwOnError:false});

and

(8)var element = document.getElementById("moose-equation-0a23567f-355e-465e-815b-7101fe37e099");katex.render("\\psi(z) = \\frac{2k T_{0}\\Sigma_{0}B}{q}\\left[W\\left(\\frac{\\phi_{0}}{2B}\\exp\\left(\\frac{\\Sigma_{0}z + C}{B}\\right)\\right)\\right]\\text{.}", element, {displayMode:true,throwOnError:false});

$var element = document.getElementById("moose-equation-1bd8534b-6598-472d-b6a0-928478343346");katex.render("B", element, {displayMode:false,throwOnError:false});$ is defined as:

(9)var element = document.getElementById("moose-equation-0d831766-2716-4ff5-9ac5-d3d17c517099");katex.render("B = -\\left[1 + \\frac{\\phi_{0}}{2}\\right]\\text{,}", element, {displayMode:true,throwOnError:false});

and $var element = document.getElementById("moose-equation-ed3cdb2f-47fe-4195-acc6-c1d5fce4f56e");katex.render("C", element, {displayMode:false,throwOnError:false});$ is:

(10)var element = document.getElementById("moose-equation-a8c269ee-c248-4fdb-b4ae-a8363ebb8038");katex.render("C = \\frac{\\phi_{0}}{2}\\text{.}", element, {displayMode:true,throwOnError:false});

$var element = document.getElementById("moose-equation-42d1a792-edb0-451c-835b-eec01b05da5b");katex.render("W(\\mathcal{Z})", element, {displayMode:false,throwOnError:false});$ is the Lambert W function (also known as the product logarithm). When using linear temperature feedback, the benchmark parameters must be selected such that $var element = document.getElementById("moose-equation-5a04894e-d8ae-4663-9061-c55b37d8ed62");katex.render("A > 0", element, {displayMode:false,throwOnError:false});$ to prevent non- physical solutions (negative cross sections). There are no such constraints on the inverse-root feedback. The cannonical benchmark parameters can be found in Table 1.

Table 1: Doppler slab benchmark parameters proposed in Remley and Griesheimer (2019).

Parameter	Value
$var element = document.getElementById("moose-equation-1f7873dc-6b28-4961-a6fe-9de0f4d2bda5");katex.render("\\psi_{0}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-6b0b2c2b-dbd4-40b1-ac81-790a13cb200d");katex.render("6.65\\times 10^{11}", element, {displayMode:false,throwOnError:false});$ (n cm^-2 s^-1)
$var element = document.getElementById("moose-equation-7e689f7d-ae56-4b60-abed-08b8aee15ac4");katex.render("N_{0}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-e021796e-a0cb-4234-bd25-431d1f9dde89");katex.render("0.025", element, {displayMode:false,throwOnError:false});$ (atom barn^-1 cm^-1)
$var element = document.getElementById("moose-equation-e53659fa-6eea-4eda-9258-6a395ca284d2");katex.render("\\sigma_{0}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-781cd053-4e76-4c33-9ba5-decde5fc9352");katex.render("4", element, {displayMode:false,throwOnError:false});$ (barn)
$var element = document.getElementById("moose-equation-4e9ae05c-07c7-4378-b3c1-02c7f9f6fd04");katex.render("q", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-3af72bab-1006-43d2-bf4f-85812f28d27e");katex.render("10^{6}", element, {displayMode:false,throwOnError:false});$ (eV)
$var element = document.getElementById("moose-equation-d24bc84b-b718-4f05-a59f-e6699b27d8f9");katex.render("T_{0}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-92b6f8a5-d0df-4a62-8c66-ffd013c565c6");katex.render("293", element, {displayMode:false,throwOnError:false});$ (K)
$var element = document.getElementById("moose-equation-7a0b105a-36d1-4f67-912c-c27ce80ba62a");katex.render("k", element, {displayMode:false,throwOnError:false});$	0.6 W (m^-1 K^-1)
$var element = document.getElementById("moose-equation-778a9b87-778e-4ace-ad7b-31e92aaca6b5");katex.render("\\alpha", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-ad87c0f8-75b1-41c6-bfff-7e189f268377");katex.render("-0.0001", element, {displayMode:false,throwOnError:false});$ (cm^-1 K^-1)

Cardinal Multi-Group Monte Carlo Solutions

Cardinal's coupling between MOOSE's heat conduction solver, OpenMC, and NekRS running as a solid heat conduction solver was verified in Hegazy and Novak (2023) with the Doppler slab problem. 200 mesh elements with linear Lagrange basis functions were used in the heat conduction solution. The OpenMC simulation used 200 cells and a mesh tally with 200 bins to ensure temperature feedback and flux tallies matched the heat conduction mesh. As the benchmark problem is semi-infinite, a length must be chosen to truncate the problem at in order to apply a symmetry boundary condition (heat conduction) and vacuum condition (neutronics). This was selected to be $var element = document.getElementById("moose-equation-1a373af8-7093-45fe-8b6f-0a2964bf4dc6");katex.render("200", element, {displayMode:false,throwOnError:false});$ cm. Results for Eq. (1) can be found in Figure 2, and results for Eq. (2) can be found in Figure 3.

Cardinal simulations of the Doppler slab problem with linear temperature feedback, taken from [!cite](doppler_slab_mc_cardinal).

Figure 2: Cardinal simulations of the Doppler slab problem with linear temperature feedback, taken from Hegazy and Novak (2023).

Cardinal simulations of the Doppler slab problem with inverse root temperature feedback, taken from [!cite](doppler_slab_mc_cardinal).

Figure 3: Cardinal simulations of the Doppler slab problem with inverse root temperature feedback, taken from Hegazy and Novak (2023).

The maximum relative error for linear feedback is $var element = document.getElementById("moose-equation-5d75dbae-4dac-4e96-bce6-d342bc6839ef");katex.render("0.0767\\%", element, {displayMode:false,throwOnError:false});$ (temperature) and $var element = document.getElementById("moose-equation-a2b6d0bc-39f9-445a-8d9a-27b58c7f9008");katex.render("4.0100\\%", element, {displayMode:false,throwOnError:false});$ (flux). The maximim relative error for inverse root feedback is $var element = document.getElementById("moose-equation-27b72800-979b-45f7-bc05-c67b36028179");katex.render("0.1500\\%", element, {displayMode:false,throwOnError:false});$ (temperature) and $var element = document.getElementById("moose-equation-aa9d7354-4149-47a6-b9ec-b490486384bf");katex.render("0.7800\\%", element, {displayMode:false,throwOnError:false});$ (flux). Error in Cardinal results is dominated by (i) statistical noise near the right boundary due to neutron attenuation, and (ii) boundary condition effects at the right boundary caused by the finite (as opposed to semi-infinite) slab. More information regarding the problem setup and results can be found in Hegazy and Novak (2023).

Input files for the linear temperature feedback case are included in Cardinal to verify the heat conduction module, Cardinal's OpenMC coupling, and different temperature interpolation schemes present in OpenMC. These input files can be found in /test/tests/neutronics/mg/mc_doppler_slab_lin. The OpenMC and heat conduction input files for this problem can be found below. Please note that these input decks run with fewer particles and Picard iterations to make them suitable for inclusion in Cardinal's test suite. Actual verification results presented in the literature using these models use higher numbers of batches and particles.

joule_per_ev = 1.60218e-19

q      = ${fparse 1.0e6}       # Energy release per absorption (eV)
Sigma0 = ${fparse 4.0 * 0.025} # Initial macroscopic XS, 1/cm
T0     = 293.6                 # Surface temperature (K)
Y0     = 6.65e11               # Source intensity (n / cm2-s)
alpha  = -0.0001               # Linear doppler Coefficient (1 / K)

!include mesh.i

[AuxVariables]
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []

  # Variables added to convert into proper heat source term for the power in the solid
  [heat_source] # W / cm^3
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []

  # Change the unit of the 'flux' (neutrons / m^2 / s) into a volumetric power. This is
  # based on Eq. (5b) in the paper, which shows the volumetric power is q * Sigma * flux
  [compute_power]
    type = ParsedAux
    variable = heat_source
    coupled_variables = 'flux temp'
    expression = 'flux * ${q} * ${joule_per_ev} * ${Sigma0} * (1+(${alpha}/${Sigma0})*(temp-${T0}))'
    execute_on = 'timestep_end'
  []
[]

[ICs]
  [temp]
    type = ConstantIC
    variable = temp
    value = ${T0}
  []
[]

[Problem]
  type = OpenMCCellAverageProblem
  verbose = true
  source_strength = ${fparse Y0 * 2.0 * 2.0} # multiply by the left area to get units of neutrons / s

  cell_level = 0
  temperature_blocks = '0'
  initial_properties = xml

  particles = 1000
  batches = 100
  inactive_batches = 0

  relaxation = 'constant'

  [Tallies]
    [heating]
      type = CellTally
      score = 'flux'
      output = 'unrelaxed_tally_std_dev'
      check_tally_sum = false
    []
  []
[]

[Executioner]
  type = Transient
  dt = 1.0
  steady_state_detection = true
  check_aux = true
  num_steps = 3
[]

[MultiApps]
  [solid]
    type = TransientMultiApp
    app_type = CardinalApp
    input_files = 'solid.i'
    execute_on = timestep_end
    sub_cycling = true
  []
[]

[Transfers]
  [solid_temp]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = T
    variable = temp
    from_multi_app = solid
  []
  [source_to_solid]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    to_multi_app = solid
  []
[]

[Postprocessors]
  [max_heat_source]
    type = ElementExtremeValue
    variable = heat_source
    execute_on = 'timestep_begin'
  []
  [max_flux]
    type = ElementExtremeValue
    variable = flux
  []
  [max_temp]
    type = ElementExtremeValue
    variable = temp
  []
[]

[Outputs]
  csv = true
  execute_on = 'TIMESTEP_END'
[]

!include mesh.i

[Variables]
  [T]
    initial_condition = 293.6
  []
[]

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
    block = '0'
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = '0'
  []
[]

[BCs]
  [left_boundary]
    type = DirichletBC
    variable = T
    boundary = 'left'
    value = 293.6
  []
[]

[Materials]
  [absorbing]
    type = HeatConductionMaterial
    thermal_conductivity = ${fparse 0.6 / 100.0} # W/cm-K
    block = '0'
  []
[]

[Postprocessors]
  [T_solid]
    type = ElementExtremeValue
    variable = T
  []
[]

[AuxVariables]
  [power]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[Executioner]
  type = Transient
  nl_abs_tol = 1e-5
  nl_rel_tol = 1e-16
  petsc_options_value = 'hypre boomeramg'
  petsc_options_iname = '-pc_type -pc_hypre_type'
[]

Cardinal Random Ray Solutions

In addition to solutions with OpenMC's Monte Carlo solver, Cardinal also includes input files for the Doppler slab problem solved with The Random Ray Method (TRRM). These solutions use the same spatial and temperature discretization parameters as the Monte Carlo solutions with a few modifications to support TRRM. This includes the use of a cell where rays are sampled, a cell where they accumulate the specified source (as TRRM only supports volume sources), and a void cell at the end of the geometry right before the vacuum boundary condition. The inactive ray distance is set to zero (as there is no initialization bias), and the active ray length is set to 209 cm (such that they are guaranteed to terminate in the void region). Inactive batches are not required as there is no scattering or fission source to converge. The results for Eq. (1) can be found in Figure 4, and results for Eq. (2) can be found in Figure 5. Both of these solutions use 400 elements, 400 OpenMC cells, and linear temperature interpolation for cross sections.

Cardinal [!ac](TRRM) simulations of the Doppler slab problem with linear temperature feedback.

Figure 4: Cardinal TRRM simulations of the Doppler slab problem with linear temperature feedback.

Cardinal [!ac](TRRM) simulations of the Doppler slab problem with inverse root temperature feedback.

Figure 5: Cardinal TRRM simulations of the Doppler slab problem with inverse root temperature feedback.

In addition to verification via direct temperature feedback, density feedback in Cardinal was verified using the analytic expression for the inverse square root cross section. The results for Eq. (2) using density feedback can be found in Figure 6.

Cardinal [!ac](TRRM) simulations of the Doppler slab problem with inverse root temperature feedback using density feedback.

Figure 6: Cardinal TRRM simulations of the Doppler slab problem with inverse root temperature feedback using density feedback.

Finally, the mesh resolution was varied between 10 and 400 elements/cells to determine the rate of convergence of the TRRM solutions. The results can be found in Figure 7, which show the expected first order rate of convergence for both the linear and inverse square root cases. In the case where temperature feedback was used, it can be seen that the $var element = document.getElementById("moose-equation-96271926-e9ba-4bb9-b2b9-e8b7afbd86be");katex.render("L_{2}", element, {displayMode:false,throwOnError:false});$ error in the solution plateaus after 100 elements due to error in the temperature interpolation scheme and the setting of cell averaged temperature values. Direct density feedback does not show the same behavior as the density update is computed analytically without any need for interpolation.

Solution error in the [!ac](TRRM) solutions as a function of the average mesh element size.

Figure 7: Solution error in the TRRM solutions as a function of the average mesh element size.

The input files for the linear temperature feedback case can be found in /test/tests/neutronics/mg/rr_doppler_slab_lin, where they are used to regression test temperature feedback with TRRM. These can be found below.

joule_per_ev = 1.60218e-19

q      = ${fparse 1.0e6}                 # energy release per absorption (eV)
Sigma0 = ${fparse 4.0 * 0.025}           # Initial macroscopic XS, 1/cm
T0     = 293.6                           # surface temperature (K)
Y0     = 6.65e11                         # source intensity (n / cm2-s)
alpha  = -0.0001                         # Linear doppler Coefficient

!include mesh.i

[AuxVariables]
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []

  # Variables added to convert into proper heat source term for the power in the solid
  [heat_source] # W / cm^3
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []

  # Change the unit of the 'flux' (neutrons / m^2 / s) into a volumetric power. This is
  # based on Eq. (5b) in the paper, which shows the volumetric power is q * Sigma * flux
  [changing_units]
    type = ParsedAux
    variable = heat_source
    coupled_variables = 'flux temp'
    expression = 'flux * ${q} * ${joule_per_ev} * ${Sigma0} * (1+(${alpha}/${Sigma0})*(temp-${T0}))'
    execute_on = 'timestep_end'
  []
[]

[ICs]
  [temp]
    type = ConstantIC
    variable = temp
    value = ${T0}
  []
[]

[Problem]
  type = OpenMCCellAverageProblem
  verbose = false
  source_strength = ${fparse Y0 * 2.0 * 2.0} # multiply by the left area to get units of neutrons / s

  cell_level = 0
  temperature_blocks = '0'
  initial_properties = xml

  particles = 1000
  batches = 100
  inactive_batches = 0

  inactive_distance = 0.0
  active_distance = 209.0

  relaxation = 'constant'

  [Tallies]
    [heating]
      type = CellTally
      score = 'flux'
      output = 'unrelaxed_tally_std_dev'
      check_tally_sum = false
    []
  []
[]

[Executioner]
  type = Transient
  dt = 1.0
  steady_state_detection = true
  check_aux = true
  num_steps = 3
[]

[MultiApps]
  [solid]
    type = TransientMultiApp
    app_type = CardinalApp
    input_files = 'solid.i'
    execute_on = timestep_end
    sub_cycling = true
  []
[]

[Transfers]
  [solid_temp]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = T
    variable = temp
    from_multi_app = solid
  []
  [source_to_solid]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    to_multi_app = solid
  []
[]

[Postprocessors]
  [max_heat_source]
    type = ElementExtremeValue
    variable = heat_source
    execute_on = 'timestep_begin'
  []
  [max_flux]
    type = ElementExtremeValue
    variable = flux
  []
  [max_temp]
    type = ElementExtremeValue
    variable = temp
  []
[]

[Outputs]
  csv = true
  execute_on = 'FINAL'
[]

!include mesh.i

[Variables]
  [T]
    initial_condition = 293.6
  []
[]

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
    block = '0'
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = '0'
  []
[]

[BCs]
  [left_boundary]
    type = DirichletBC
    variable = T
    boundary = 'left'
    value = 293.6
  []
[]

[Materials]
  [absorbing]
    type = HeatConductionMaterial
    thermal_conductivity = ${fparse 0.6 / 100.0} # W/cm-K
    block = '0'
  []
[]

[Postprocessors]
  [T_solid]
    type = ElementExtremeValue
    variable = T
  []
[]

[AuxVariables]
  [power]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[Executioner]
  type = Transient
  nl_abs_tol = 1e-5
  nl_rel_tol = 1e-16
  petsc_options_value = 'hypre boomeramg'
  petsc_options_iname = '-pc_type -pc_hypre_type'
[]

The input files for the inverse root temperature feedback case can be found in /test/tests/neutronics/mg/rr_doppler_slab_root, where they are used to regression test density feedback with TRRM.

joule_per_ev = 1.60218e-19

q      = ${fparse 1.0e6}                 # energy release per absorption (eV)
Sigma0 = ${fparse 4.0 * 0.025}           # Initial macroscopic XS, 1/cm
T0     = 293.6                           # surface temperature (K)
Y0     = 6.65e11                         # source intensity (n / cm2-s)

!include mesh.i

[AuxVariables]
  [temp]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_density]
    family = MONOMIAL
    order = CONSTANT
  []

  # Variables added to convert into proper heat source term for the power in the solid
  [heat_source] # W / cm^3
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []

  # Change the unit of the 'flux' (neutrons / m^2 / s) into a volumetric power. This is
  # based on Eq. (5b) in the paper, which shows the volumetric power is q * Sigma * flux
  [changing_units]
    type = ParsedAux
    variable = heat_source
    coupled_variables = 'flux temp'
    expression = 'flux * ${q} * ${joule_per_ev} * ${Sigma0} * sqrt(${T0} / temp)'
    execute_on = 'timestep_end'
  []

  # Compute the "density" with inverse square root Doppler boadening
  [compute_density]
    type = ParsedAux
    variable = density
    coupled_variables = 'temp'
    expression = 'sqrt(${T0} / temp)'
    execute_on = 'timestep_end'
  []
[]

[ICs]
  [density]
    type = ConstantIC
    variable = density
    value = 1.0
  []
  [temp]
    type = ConstantIC
    variable = temp
    value = ${T0}
  []
[]

[Problem]
  type = OpenMCCellAverageProblem
  verbose = false
  source_strength = ${fparse Y0 * 2.0 * 2.0} # multiply by the left area to get units of neutrons / s

  cell_level = 0
  density_blocks = '0'
  mgxs_reference_densities_by_block = '1.0'
  initial_properties = xml

  particles = 1000
  batches = 100
  inactive_batches = 0

  inactive_distance = 0.0
  active_distance = 209.0

  relaxation = 'constant'

  [Tallies]
    [heating]
      type = CellTally
      score = 'flux'
      output = 'unrelaxed_tally_std_dev'
      check_tally_sum = false
    []
  []
[]

[Executioner]
  type = Transient
  dt = 1.0
  steady_state_detection = true
  check_aux = true
  num_steps = 3
[]

[MultiApps]
  [solid]
    type = TransientMultiApp
    app_type = CardinalApp
    input_files = 'solid.i'
    execute_on = timestep_end
    sub_cycling = true
  []
[]

[Transfers]
  [solid_temp]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = T
    variable = temp
    from_multi_app = solid
    error_on_miss = false
  []
  [source_to_solid]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    to_multi_app = solid
  []
[]

[Postprocessors]
  [max_heat_source]
    type = ElementExtremeValue
    variable = heat_source
    execute_on = 'timestep_begin'
  []
  [max_flux]
    type = ElementExtremeValue
    variable = flux
  []
[]

[Outputs]
  csv = true
  execute_on = 'FINAL'
[]

!include mesh.i

[Variables]
  [T]
    initial_condition = 293.6
  []
[]

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
    block = '0'
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = '0'
  []
[]

[BCs]
  [left_boundary]
    type = DirichletBC
    variable = T
    boundary = 'left'
    value = 293.6
  []
[]

[Materials]
  [absorbing]
    type = HeatConductionMaterial
    thermal_conductivity = ${fparse 0.6 / 100.0} # W/cm-K
    block = '0'
  []
[]

[Postprocessors]
  [T_solid]
    type = ElementExtremeValue
    variable = T
  []
[]

[AuxVariables]
  [power]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[Executioner]
  type = Transient
  nl_abs_tol = 1e-5
  nl_rel_tol = 1e-16
  petsc_options_value = 'hypre boomeramg'
  petsc_options_iname = '-pc_type -pc_hypre_type'
[]

References

Aya H. Hegazy and April J. Novak. Verification of the Cardinal Multiphysics Solver with the Doppler Slab Benchmark. In Proceedings of M&C. 2023.

BibTeX

@InProceedings{doppler_slab_mc_cardinal,
    author = "Hegazy, Aya H. and Novak, April J.",
    title = "{Verification of the Cardinal Multiphysics Solver with the Doppler Slab Benchmark}",
    booktitle = "{Proceedings of M\\&C}",
    year = "2023"
}

RIS

TY  - CPAPER
AU  - Hegazy, Aya H.
AU  - Novak, April J.
TI  - Verification of the Cardinal Multiphysics Solver with the Doppler Slab Benchmark
T2  - Proceedings of M&C
PY  - 2023
ER  -

Plain Text

Aya H. Hegazy and April J. Novak. Verification of the Cardinal Multiphysics Solver with the Doppler Slab Benchmark. In Proceedings of M&C. 2023.

Kyle E. Remley and David P. Griesheimer. A Fully Analytic Coupled Thermal-Neutronics Benchmark and its Application to Monte Carlo Simulation. In Proceedings of M&C. 2019.

BibTeX

@InProceedings{doppler_slab_benchmark,
    author = "Remley, Kyle E. and Griesheimer, David P.",
    title = "{A Fully Analytic Coupled Thermal-Neutronics Benchmark and its Application to Monte Carlo Simulation}",
    booktitle = "{Proceedings of M\\&C}",
    year = "2019"
}

RIS

TY  - CPAPER
AU  - Remley, Kyle E.
AU  - Griesheimer, David P.
TI  - A Fully Analytic Coupled Thermal-Neutronics Benchmark and its Application to Monte Carlo Simulation
T2  - Proceedings of M&C
PY  - 2019
ER  -

Plain Text

Kyle E. Remley and David P. Griesheimer. A Fully Analytic Coupled Thermal-Neutronics Benchmark and its Application to Monte Carlo Simulation. In Proceedings of M&C. 2019.

(test/tests/neutronics/mg/mc_doppler_slab_lin/openmc_base.i)

joule_per_ev = 1.60218e-19

q      = ${fparse 1.0e6}       # Energy release per absorption (eV)
Sigma0 = ${fparse 4.0 * 0.025} # Initial macroscopic XS, 1/cm
T0     = 293.6                 # Surface temperature (K)
Y0     = 6.65e11               # Source intensity (n / cm2-s)
alpha  = -0.0001               # Linear doppler Coefficient (1 / K)

!include mesh.i

[AuxVariables]
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []

  # Variables added to convert into proper heat source term for the power in the solid
  [heat_source] # W / cm^3
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []

  # Change the unit of the 'flux' (neutrons / m^2 / s) into a volumetric power. This is
  # based on Eq. (5b) in the paper, which shows the volumetric power is q * Sigma * flux
  [compute_power]
    type = ParsedAux
    variable = heat_source
    coupled_variables = 'flux temp'
    expression = 'flux * ${q} * ${joule_per_ev} * ${Sigma0} * (1+(${alpha}/${Sigma0})*(temp-${T0}))'
    execute_on = 'timestep_end'
  []
[]

[ICs]
  [temp]
    type = ConstantIC
    variable = temp
    value = ${T0}
  []
[]

[Problem]
  type = OpenMCCellAverageProblem
  verbose = true
  source_strength = ${fparse Y0 * 2.0 * 2.0} # multiply by the left area to get units of neutrons / s

  cell_level = 0
  temperature_blocks = '0'
  initial_properties = xml

  particles = 1000
  batches = 100
  inactive_batches = 0

  relaxation = 'constant'

  [Tallies]
    [heating]
      type = CellTally
      score = 'flux'
      output = 'unrelaxed_tally_std_dev'
      check_tally_sum = false
    []
  []
[]

[Executioner]
  type = Transient
  dt = 1.0
  steady_state_detection = true
  check_aux = true
  num_steps = 3
[]

[MultiApps]
  [solid]
    type = TransientMultiApp
    app_type = CardinalApp
    input_files = 'solid.i'
    execute_on = timestep_end
    sub_cycling = true
  []
[]

[Transfers]
  [solid_temp]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = T
    variable = temp
    from_multi_app = solid
  []
  [source_to_solid]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    to_multi_app = solid
  []
[]

[Postprocessors]
  [max_heat_source]
    type = ElementExtremeValue
    variable = heat_source
    execute_on = 'timestep_begin'
  []
  [max_flux]
    type = ElementExtremeValue
    variable = flux
  []
  [max_temp]
    type = ElementExtremeValue
    variable = temp
  []
[]

[Outputs]
  csv = true
  execute_on = 'TIMESTEP_END'
[]

(test/tests/neutronics/mg/mc_doppler_slab_lin/solid.i)

!include mesh.i

[Variables]
  [T]
    initial_condition = 293.6
  []
[]

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
    block = '0'
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = '0'
  []
[]

[BCs]
  [left_boundary]
    type = DirichletBC
    variable = T
    boundary = 'left'
    value = 293.6
  []
[]


[Materials]
  [absorbing]
    type = HeatConductionMaterial
    thermal_conductivity = ${fparse 0.6 / 100.0} # W/cm-K
    block = '0'
  []
[]

[Postprocessors]
  [T_solid]
    type = ElementExtremeValue
    variable = T
  []
[]

[AuxVariables]
  [power]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[Executioner]
  type = Transient
  nl_abs_tol = 1e-5
  nl_rel_tol = 1e-16
  petsc_options_value = 'hypre boomeramg'
  petsc_options_iname = '-pc_type -pc_hypre_type'
[]

(test/tests/neutronics/mg/rr_doppler_slab_lin/openmc_base.i)

joule_per_ev = 1.60218e-19

q      = ${fparse 1.0e6}                 # energy release per absorption (eV)
Sigma0 = ${fparse 4.0 * 0.025}           # Initial macroscopic XS, 1/cm
T0     = 293.6                           # surface temperature (K)
Y0     = 6.65e11                         # source intensity (n / cm2-s)
alpha  = -0.0001                         # Linear doppler Coefficient

!include mesh.i

[AuxVariables]
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []

  # Variables added to convert into proper heat source term for the power in the solid
  [heat_source] # W / cm^3
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []

  # Change the unit of the 'flux' (neutrons / m^2 / s) into a volumetric power. This is
  # based on Eq. (5b) in the paper, which shows the volumetric power is q * Sigma * flux
  [changing_units]
    type = ParsedAux
    variable = heat_source
    coupled_variables = 'flux temp'
    expression = 'flux * ${q} * ${joule_per_ev} * ${Sigma0} * (1+(${alpha}/${Sigma0})*(temp-${T0}))'
    execute_on = 'timestep_end'
  []
[]


[ICs]
  [temp]
    type = ConstantIC
    variable = temp
    value = ${T0}
  []
[]

[Problem]
  type = OpenMCCellAverageProblem
  verbose = false
  source_strength = ${fparse Y0 * 2.0 * 2.0} # multiply by the left area to get units of neutrons / s

  cell_level = 0
  temperature_blocks = '0'
  initial_properties = xml

  particles = 1000
  batches = 100
  inactive_batches = 0

  inactive_distance = 0.0
  active_distance = 209.0

  relaxation = 'constant'

  [Tallies]
    [heating]
      type = CellTally
      score = 'flux'
      output = 'unrelaxed_tally_std_dev'
      check_tally_sum = false
    []
  []
[]

[Executioner]
  type = Transient
  dt = 1.0
  steady_state_detection = true
  check_aux = true
  num_steps = 3
[]

[MultiApps]
  [solid]
    type = TransientMultiApp
    app_type = CardinalApp
    input_files = 'solid.i'
    execute_on = timestep_end
    sub_cycling = true
  []
[]

[Transfers]
  [solid_temp]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = T
    variable = temp
    from_multi_app = solid
  []
  [source_to_solid]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    to_multi_app = solid
  []
[]

[Postprocessors]
  [max_heat_source]
    type = ElementExtremeValue
    variable = heat_source
    execute_on = 'timestep_begin'
  []
  [max_flux]
    type = ElementExtremeValue
    variable = flux
  []
  [max_temp]
    type = ElementExtremeValue
    variable = temp
  []
[]

[Outputs]
  csv = true
  execute_on = 'FINAL'
[]

(test/tests/neutronics/mg/rr_doppler_slab_lin/solid.i)

!include mesh.i

[Variables]
  [T]
    initial_condition = 293.6
  []
[]

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
    block = '0'
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = '0'
  []
[]

[BCs]
  [left_boundary]
    type = DirichletBC
    variable = T
    boundary = 'left'
    value = 293.6
  []
[]

[Materials]
  [absorbing]
    type = HeatConductionMaterial
    thermal_conductivity = ${fparse 0.6 / 100.0} # W/cm-K
    block = '0'
  []
[]

[Postprocessors]
  [T_solid]
    type = ElementExtremeValue
    variable = T
  []
[]

[AuxVariables]
  [power]
    family = MONOMIAL
    order = CONSTANT
  []
[]


[Executioner]
  type = Transient
  nl_abs_tol = 1e-5
  nl_rel_tol = 1e-16
  petsc_options_value = 'hypre boomeramg'
  petsc_options_iname = '-pc_type -pc_hypre_type'
[]

(test/tests/neutronics/mg/rr_doppler_slab_root/openmc.i)

joule_per_ev = 1.60218e-19

q      = ${fparse 1.0e6}                 # energy release per absorption (eV)
Sigma0 = ${fparse 4.0 * 0.025}           # Initial macroscopic XS, 1/cm
T0     = 293.6                           # surface temperature (K)
Y0     = 6.65e11                         # source intensity (n / cm2-s)

!include mesh.i

[AuxVariables]
  [temp]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_density]
    family = MONOMIAL
    order = CONSTANT
  []

  # Variables added to convert into proper heat source term for the power in the solid
  [heat_source] # W / cm^3
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []

  # Change the unit of the 'flux' (neutrons / m^2 / s) into a volumetric power. This is
  # based on Eq. (5b) in the paper, which shows the volumetric power is q * Sigma * flux
  [changing_units]
    type = ParsedAux
    variable = heat_source
    coupled_variables = 'flux temp'
    expression = 'flux * ${q} * ${joule_per_ev} * ${Sigma0} * sqrt(${T0} / temp)'
    execute_on = 'timestep_end'
  []

  # Compute the "density" with inverse square root Doppler boadening
  [compute_density]
    type = ParsedAux
    variable = density
    coupled_variables = 'temp'
    expression = 'sqrt(${T0} / temp)'
    execute_on = 'timestep_end'
  []
[]


[ICs]
  [density]
    type = ConstantIC
    variable = density
    value = 1.0
  []
  [temp]
    type = ConstantIC
    variable = temp
    value = ${T0}
  []
[]

[Problem]
  type = OpenMCCellAverageProblem
  verbose = false
  source_strength = ${fparse Y0 * 2.0 * 2.0} # multiply by the left area to get units of neutrons / s

  cell_level = 0
  density_blocks = '0'
  mgxs_reference_densities_by_block = '1.0'
  initial_properties = xml

  particles = 1000
  batches = 100
  inactive_batches = 0

  inactive_distance = 0.0
  active_distance = 209.0

  relaxation = 'constant'

  [Tallies]
    [heating]
      type = CellTally
      score = 'flux'
      output = 'unrelaxed_tally_std_dev'
      check_tally_sum = false
    []
  []
[]

[Executioner]
  type = Transient
  dt = 1.0
  steady_state_detection = true
  check_aux = true
  num_steps = 3
[]

[MultiApps]
  [solid]
    type = TransientMultiApp
    app_type = CardinalApp
    input_files = 'solid.i'
    execute_on = timestep_end
    sub_cycling = true
  []
[]

[Transfers]
  [solid_temp]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = T
    variable = temp
    from_multi_app = solid
    error_on_miss = false
  []
  [source_to_solid]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    to_multi_app = solid
  []
[]

[Postprocessors]
  [max_heat_source]
    type = ElementExtremeValue
    variable = heat_source
    execute_on = 'timestep_begin'
  []
  [max_flux]
    type = ElementExtremeValue
    variable = flux
  []
[]

[Outputs]
  csv = true
  execute_on = 'FINAL'
[]

(test/tests/neutronics/mg/rr_doppler_slab_root/solid.i)

!include mesh.i

[Variables]
  [T]
    initial_condition = 293.6
  []
[]

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
    block = '0'
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = '0'
  []
[]

[BCs]
  [left_boundary]
    type = DirichletBC
    variable = T
    boundary = 'left'
    value = 293.6
  []
[]


[Materials]
  [absorbing]
    type = HeatConductionMaterial
    thermal_conductivity = ${fparse 0.6 / 100.0} # W/cm-K
    block = '0'
  []
[]

[Postprocessors]
  [T_solid]
    type = ElementExtremeValue
    variable = T
  []
[]

[AuxVariables]
  [power]
    family = MONOMIAL
    order = CONSTANT
  []
[]


[Executioner]
  type = Transient
  nl_abs_tol = 1e-5
  nl_rel_tol = 1e-16
  petsc_options_value = 'hypre boomeramg'
  petsc_options_iname = '-pc_type -pc_hypre_type'
[]

Problem Description
Cardinal Multi - Group Monte Carlo Solutions
Cardinal Random Ray Solutions
References