Griffin steady state input

We will first cover the input for the main application, Griffin, which tackles the neutronics problem. While the Griffin manual is ultimately the most complete reference on this input, we will try to provide enough details here for a complete comprehension of this input.

We first define in the header:

• initial conditions for the temperatures. These may be specified directly in the AuxVariables block, and will be overriden by the values provided by the thermal hydraulics solve.

• the core power.

• the percent of power generated by decay heat. There is no decay heat generation in this input, though Griffin has the capability to handle decay heat generation.

The next block is a GlobalParams block. These parameters will be inserted in every block in the input, and are used to reduce the size of the input file. Here we specify the group cross section library. We will give additional details about the group cross sections in the Materials block.

[GlobalParams]
  library_file = "cross_sections/2-D_8Gt_multiregions_transient.xml"
  library_name = 2-D_8Gt
[]

The [TransportSystems] block is used to specify the solver parameters. We chose a diffusion solver as accuracy is generally satisfactory with graphite-moderated reactors, as we confirmed by benchmarking with Serpent (Giudicelli et al., 2021). We select the eigenvalue equation type as the steady state coupling is an eigenvalue calculation for neutron transport. The steady state is the eigenpair of the and the steady flux distribution. We also specify the boundary conditions in this block. This is a 2D RZ model, so the center of the geometry is an axis of symmetry. A vacuum boundary condition is placed on the other boundaries. This approximation is appropriate when the boundaries are sufficiently far from the active region.

[TransportSystems]
  particle = neutron
  equation_type = eigenvalue

  G = 8

  ReflectingBoundary = 'reflector_surface'
  VacuumBoundary = 'brick_surface top bottom'

  [diff]
    scheme = CFEM-Diffusion
    n_delay_groups = 6
    family = LAGRANGE
    order = FIRST
    # verbose = 2
  []
[]

The next part of the input specifies the geometry. The geometry in MOOSE is described by an unstructured mesh. MOOSE also contains basic mesh generation capabilities which can be leveraged to generated rectilinear meshes. However, for realistic reactor models, it is more common to use an external meshing software. We used Cubit to generate meshes for the Mk1-FHR. The mesh may be shared between applications, but it is generally recommended to tailor the mesh to the physics equations being solved. We specify the external mesh using a FileMeshGenerator.

[Mesh]
  [mesh_reader]
    type = FileMeshGenerator
    file = 'meshes/core_pronghorn.e'
  []
[]

This mesh is then modified using MOOSE mesh generators. We remove the brick insulation from the system using a BlockDeletionGenerator as its influence on the neutronics is negligible. Before that, we added a sideset, a collection of element sides, to the mesh. This new sideset is used to define the boundary condition for the outer boundary of the core.

[Mesh]
  [new_boundary]
    type = SideSetsBetweenSubdomainsGenerator
    input = mesh_reader
    primary_block = '9'
    paired_block = '10'
    new_boundary = 'brick_surface'
  []
[]

Variables and kernels need not be defined when using Griffin, as the TransportSystems block is an Action in MOOSE vocabulary that takes care of defining those. We will explain those blocks more in details for the Pronghorn input, which does not use an Action.

AuxVariables are used to define fields/variables that are not solved for. They may be used for coupling, for computing material properties, for outputting quantities of interests, and many other uses. In this input file, Tfuel and Tsalt are used to couple to the thermal hydraulics simulation. These fields are populated by Transfers from the Pronghorn app, further in the input file.

[AuxVariables]
  [Tfuel]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${initial_fuel_temperature}
    block = '1 2 3 4 5 6 7 8 9' #FIXME
  []

  [Tsalt]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${initial_salt_temperature}
    block = '1 2 3 4 5 6 7 8 9' #FIXME
  []
[]

The CR_insertion auxiliary variable is the insertion of the control rod in the system.

[AuxVariables]
  [CR_insertion]
    family = MONOMIAL
    order = CONSTANT
    block = '1 2'
  []
[]

The next three AuxVariables are the fission, decay heat and total power densities. They are computed from inst_power_density which is computed by Griffin. The decay heat is important to track for loss of flow transients, the total power density will be required as the heat source in the thermal hydraulics simulation.

[AuxVariables]
  [fission_power_density]
    # Energy fraction released Instantaneously from fission (W).
    family = MONOMIAL
    order = CONSTANT
  []

  [decay_heat_power_density]
    # Energy fraction released by fission products (W).
    family = MONOMIAL
    order = CONSTANT

    # Energy fraction released by fission products (W).
    family = MONOMIAL
    order = CONSTANT
  []
[]

The fluxes output by Griffin are normalized with regards to the solver fission source. To output the physical fluxes normalized to the real core power, which may be used to compute material damage in the fuel for example, we define the scaled_sflux_gi AuxVariables.

[AuxVariables]
  [scaled_sflux_g0]
    family = LAGRANGE
    order = FIRST
  []
[]

AuxKernels are used to operate on AuxVariables. They may scale, multiply, add and perform many other operations. They may be block restricted, as some AuxVariables are not defined over the entire domain, or they may inherit the same block restriction as the variable.

[AuxKernels]
  [fission_power_density]
    type = ScaleAux
    multiplier = ${fis_fract}
    variable = fission_power_density
    source_variable = inst_power_density
  []
  [decay_heat_power_density]
    type = ScaleAux
    multiplier = ${dh_fract}
    variable = decay_heat_power_density
    source_variable = inst_power_density
  []
  [total_power_density]
    type = ParsedAux
    expression = 'fission_power_density + decay_heat_power_density'
    coupled_variables = 'fission_power_density decay_heat_power_density'
    variable = total_power_density
    execute_on = 'INITIAL timestep_end'
  []
  [scale_flux_g0]
    type = ScaleAux
    multiplying_pp = power_scaling
    source_variable = sflux_g0
    variable = scaled_sflux_g0
  []
  [scale_flux_g1]
    type = ScaleAux
    multiplying_pp = power_scaling
    source_variable = sflux_g1
    variable = scaled_sflux_g1
  []
  [scale_flux_g2]
    type = ScaleAux
    multiplying_pp = power_scaling
    source_variable = sflux_g2
    variable = scaled_sflux_g2
  []
  [scale_flux_g3]
    type = ScaleAux
    multiplying_pp = power_scaling
    source_variable = sflux_g3
    variable = scaled_sflux_g3
  []
  [scale_flux_g4]
    type = ScaleAux
    multiplying_pp = power_scaling
    source_variable = sflux_g4
    variable = scaled_sflux_g4
  []
  [scale_flux_g5]
    type = ScaleAux
    multiplying_pp = power_scaling
    source_variable = sflux_g5
    variable = scaled_sflux_g5
  []
  [scale_flux_g6]
    type = ScaleAux
    multiplying_pp = power_scaling
    source_variable = sflux_g6
    variable = scaled_sflux_g6
  []
  [scale_flux_g7]
    type = ScaleAux
    multiplying_pp = power_scaling
    source_variable = sflux_g7
    variable = scaled_sflux_g7
  []
  [rod_insertion]
    type = FunctionAux
    function = control_rod_position
    variable = 'CR_insertion'
  []
[]

Functions may be defined in MOOSE in a Functions block. The control rod position is set in this input using a function. Since it's a steady state calculation, we could also have simply used an initial condition. The FunctionAux is using the function to set the CR_insertion auxiliary variable.

[Functions]
  [control_rod_position]
    type = PiecewiseLinear
    x = '0 10'
    y = '2 2'
  []
[]

Grffin uses the PowerDensity as a shorthand to define variables and postprocessors related to the power of the core. The following fields are used:

• power: to define the numerical value of the core power

• power_density_variable: a variable for the power density distribution

• integrated_power_postprocessor: a postprocessor to compute the total core power

• power_scaling_postprocessor: a postprocessor used to scale the fluxes to the real core power

The CONSTANT MONOMIAL basis is used to represent the power distribution. This will be able to match the discontinuities in the fission group cross sections

[PowerDensity]
  power = '${fparse total_power * fis_fract}'
  power_density_variable = inst_power_density
  integrated_power_postprocessor = total_power
  power_scaling_postprocessor = power_scaling
  family = MONOMIAL
  order = CONSTANT
[]

The group cross secions are distributed through the geometry using the Materials block. Each block in this section is a Material object defining the group cross section in a block. We will go over the pebble_bed material, the group cross section of the salt-pebble homogenized mixture in the bed.

The CoupledFeedbackNeutronicsMaterial is used to take into account the effect of the fuel and salt temperature on the group cross section. The cross section library name and the file containing it are specified using the GlobalParams block. The temperature dependence is captured using a tabulation. The names used in the tabulation are matched to the variables in the simulation using grid_names and grid_variables. The correct entry in the library is selected using the material_id. Since macroscopic group cross sections are used, the isotopes aren't specified and a pseudo-isotope is used.

[Materials]
  [inner_reflector]
    type = CoupledFeedbackNeutronicsMaterial #FIXME
    grid_names = 'Tfuel Tsalt CR'
    grid_variables = 'Tfuel Tsalt CR_insertion'
    plus = true
    diffusion_coefficient_scheme = local
    isotopes = 'pseudo'
    densities = '1.0'
    is_meter = true
    material_id = 100000
    block = 1
  []
[]

The boundary conditions were already specified by the TransportSystems action so no additional boundary condition is specified.

The Executioner block specifies how to solve the equation system. We choose an eigenvalue executioner, PicardEigen, as a criticality calculation is an eigenvalue problem.

The number of non-linear iterations and the non-linear relative and absolute convergence criteria may be respectively reduced and loosened to reduce the computational cost of the solution at the expense of its convergence.

[Executioner]
  type = Eigenvalue

  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart'
  petsc_options_value = 'hypre boomeramg 100'
  line_search = none
  snesmf_reuse_base = false

  # Linear/nonlinear iterations
  l_max_its = 100
  l_tol = 1e-3

  nl_max_its = 50
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-5

  # Power iterations, before the non-linear solve
  free_power_iterations = 2

  # Fixed point iterations
  fixed_point_abs_tol = 1e-6
  fixed_point_rel_tol = 1e-6
  fixed_point_max_its = 10
[]

The Pronghorn sub-application is created by the MultiApps block. Since we are seeking a steady state solution, we want the thermal hydraulics problem to be fully solved at every multiphysics iteration. This is accomplished by a FullSolveMultiApp. The execute_on field is set to timestep_end to have the thermal hydraulics solve be performed after the neutronics solve. This means that for the first multiphysics iteration, the temperature fields in the neutronics solve will be set to an initial guess, while the power distribution in the thermal hydraulics solve will be updated once already. 0 0 0 is specified for the position since the two meshes are aligned.

[MultiApps]
  [thermo]
    type = FullSolveMultiApp
    input_files = 'ss1_combined.i'
    execute_on = 'timestep_end' # once power distribution is computed
    positions = '0.0 0.0 0.0'
  []
[]

The coupling to the thermal hydraulics simulation is done by using a MultiAppProjectionTransfer of the power density. This transfer is conservative, in that the total power in both applications is preserved. This is important for the multiphysics coupling, as the power will not fluctuate in the thermal hydraulics simulation based on the power density profile and the difference between the neutronics and thermal hydraulics meshes. MultiAppInterpolationTransfer are used to transfer the solid and fluid phase temperature fields. Conservation is not as important in this direction, as discretization errors on the temperature do not propagate to large errors in the flux distribution.

[Transfers]
  [power_tosub]
    type = MultiAppProjectionTransfer
    to_multi_app = thermo
    source_variable = total_power_density
    variable = power_distribution
    to_postprocessors_to_be_preserved = power
    from_postprocessors_to_be_preserved = total_power
    execute_on = 'timestep_end' # overwrite initial condition
  []
  [Tcoolant_fromsub]
    type = MultiAppGeometricInterpolationTransfer
    from_multi_app = thermo
    source_variable = T_fluid
    variable = Tsalt
    execute_on = 'timestep_end'
  []
  [Tfuel_fromsub]
    type = MultiAppGeometricInterpolationTransfer
    from_multi_app = thermo
    source_variable = T_solid
    variable = Tfuel
    execute_on = 'timestep_end'
  []
  [num_fixed_point]
    type = MultiAppPostprocessorTransfer
    to_multi_app = thermo
    from_postprocessor = num_fixed_point
    to_postprocessor = num_fixed_point
  []
[]

The [RestartVariables] block is commented out, but it may be used to quickly recover the solution from a previous run instead of calculating it again. It is important that both main and auxiliary variables are recovered. As the multiapp is executed at the end of the timestep, the restarted solver would not have access to the coupled calculation temperature fields when it restarts.

Finally, the Outputs block indicates which types of outputs the simulation should return. We specify here three types:

• exodus. Exodus files contain the mesh, the (aux)variables and global quantities such as postprocessors. We can use Paraview to view the spatial dependence of each field. They may be used to restart simulations.

• comma-separated values, or csv. This file reports the values of the postprocessors at each timestep. They are useful for easily importing then plotting the simulation results using Python or Matlab.

• checkpoint files. Checkpoint files are used solely for restarting simulations. Unlike Exodus files, they also contain all the subapp information, so only the main app checkpoint files are necessary. Checkpoint files are output at a regular interval, specified by the num_files field.

[Outputs]
  exodus = true
  csv = true
  [Checkpoint]
    type = Checkpoint
    execute_on = 'FINAL'
  []
  # Reduce base output
  print_linear_converged_reason = false
  print_linear_residuals = false
  print_nonlinear_converged_reason = false
[]