Multiphysics Model

Contact: Nicolas Stauff (nstauff.at.anl.gov), Yinbin Miao (ymiao.at.anl.gov)

note:Acknowledgement

This HP-MR model was built upon earlier work performed under ARPA-E MEITNER project and reported in the journal paper (Matthews et al., 2021), and some parts of the inputs are coming from these original models.

Mesh File

As a series of powerful reactor-focused mesh generators have been implemented into MOOSE as the Reactor Module (Shemon et al., 2022), the meshing tool used for the HP-MR simulation has been moved from Cubit to MOOSE's intrinsic mesh generator set.

While Sockeye uses its own simple implemented 2D axisymetric heat pipe meshing capabilities, both Griffin and BISON applications adopt meshes generated by MOOSE's intrinsic mesh generators in this HP-MR model. The input file used to generate the Griffin mesh is listed as follows.

################################################################
# This is the MOOSE input file to generate the mesh that can be
# used for the 1/6 core Heat-Pipe Microreactor Multiphysics
# simulations (Griffin neutronics simulation).
# Running this input requires MOOSE Reactor Module Objects
# Users should use
# --mesh-only HPMR_OneSixth_Core_meshgenerator_tri_rotate_bdry.e
# command line argument to generate
# exodus file for further Multiphysics simulations.
################################################################
[GlobalParams]
  create_outward_interface_boundaries = false
[]

[Mesh]
  [BatchMeshGeneratorAction]
    # We generate each control drum in 2 steps (concentric circles generation and azimuthal absorber/reflector splitting),
    # this syntax allows us to input the control drum (and reflector) generation in a very compact manner.
    [control_drums_ccg]
      mesh_generator_type = 'HexagonConcentricCircleAdaptiveBoundaryMeshGenerator'
      mesh_name_prefix = 'cd0'
      multi_batch_params_method = 'corresponding'
      # The main different parameters between the control drums are how the meshes are adapted to the reactor core
      # The batch parameters here are used to account for these differences
      batch_vector_input_param_names = 'meshes_to_adapt_to sides_to_adapt'
      batch_vector_input_param_types = 'MGNAME UINT'
      batch_vector_input_param_values = 'Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned|
                                         2 3 4;2 3;1 2 3;1 2;0 1 2;0 1;
                                         0 1 5;0 5;0 4 5;4 5;3 4 5;3 4'
      # The rest of the control drum parameters are the same for all control drums
      fixed_vector_input_param_names = 'num_sectors_per_side background_block_ids ring_radii ring_intervals ring_block_ids'
      fixed_vector_input_param_types = 'UINT USHORT REAL UINT USHORT'
      fixed_vector_input_param_values = '4 4 4 4 4 4;
                                         504;
                                         12.25 13.25;
                                         2 1;
                                         500 501 502'
      fixed_scalar_input_param_names = 'hexagon_size background_intervals preserve_volumes is_control_drum create_outward_interface_boundaries'
      fixed_scalar_input_param_types = 'REAL UINT BOOL BOOL BOOL'
      fixed_scalar_input_param_values = '13.376 2 true true false'
    []
    [control_drums_split]
      mesh_generator_type = 'AzimuthalBlockSplitGenerator'
      mesh_name_prefix = 'cd'
      multi_batch_params_method = 'corresponding'
      # Based on the locations of the control drums, we need to split the blocks in different ways
      batch_scalar_input_param_names = 'input start_angle'
      batch_scalar_input_param_types = 'MGNAME REAL'
      batch_scalar_input_param_values = 'cd0_0 cd0_1 cd0_2 cd0_3 cd0_4 cd0_5 cd0_6 cd0_7 cd0_8 cd0_9 cd0_10 cd0_11;
                                         15 345 315 285 255 225 195 165 135 105 75 45'
      fixed_vector_input_param_names = 'old_blocks new_block_ids'
      fixed_vector_input_param_types = 'SDNAME USHORT'
      fixed_vector_input_param_values = '502;503'
      fixed_scalar_input_param_names = 'angle_range'
      fixed_scalar_input_param_types = 'REAL'
      fixed_scalar_input_param_values = '90'
    []
    [reflector_blocks]
      mesh_generator_type = 'HexagonConcentricCircleAdaptiveBoundaryMeshGenerator'
      mesh_name_prefix = 'ref'
      multi_batch_params_method = 'corresponding'
      # Similar to the control drums, we need to adapt the reflector meshes to the reactor core based on the locations
      batch_vector_input_param_names = 'sides_to_adapt'
      batch_vector_input_param_types = 'UINT'
      batch_vector_input_param_values = '0;1;2;3;4;5'
      fixed_vector_input_param_names = 'meshes_to_adapt_to num_sectors_per_side background_block_ids'
      fixed_vector_input_param_types = 'MGNAME UINT USHORT'
      fixed_vector_input_param_values = 'Patterned;
                                         4 4 4 4 4 4;
                                         400 401'
      fixed_scalar_input_param_names = 'hexagon_size background_intervals create_outward_interface_boundaries'
      fixed_scalar_input_param_types = 'REAL UINT BOOL'
      fixed_scalar_input_param_values = '13.376 2 false'
    []
  []
  # Moderator pins
  [Mod_hex]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2 '
    background_intervals = 1
    background_block_ids = '10'
    polygon_size = 1.15
    polygon_size_style ='apothem'
    ring_radii = '0.825 0.92'
    ring_intervals = '2 1'
    ring_block_ids = '103 100 101' # 103 is tri mesh
    preserve_volumes = on
    quad_center_elements = false
  []
  # Heat Pipe pins
  [HP_hex]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2 '
    background_intervals = 1
    background_block_ids = '10'
    polygon_size = 1.15
    polygon_size_style ='apothem'
    ring_radii = '0.97 1.07'
    ring_intervals = '2 1'
    ring_block_ids = '203 200 201' # 203 is tri mesh
    preserve_volumes = on
    quad_center_elements = false
  []
  # Fuel pins
  [Fuel_hex]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2 '
    background_intervals = 1
    background_block_ids = '10'
    polygon_size = 1.15
    polygon_size_style ='apothem'
    ring_radii = '1'
    ring_intervals = '2'
    ring_block_ids = '303 301'  # 303 is tri mesh
    preserve_volumes = on
    quad_center_elements = false
  []
  # Assembly mesh with all the three types of pins
  [Patterned]
    type = PatternedHexMeshGenerator
    inputs = 'Fuel_hex HP_hex Mod_hex'
    hexagon_size = 13.376
    background_block_id = 10
    background_intervals = 1
    pattern = '1 0 1 0 1 0 1;
              0 2 0 2 0 2 0 0;
             1 0 1 0 1 0 1 2 1;
            0 2 0 2 0 2 0 0 0 0;
           1 0 1 0 1 0 1 2 1 2 1;
          0 2 0 2 0 2 0 0 0 0 0 0;
         1 0 1 0 1 0 1 2 1 2 1 2 1;
          0 2 0 2 0 2 0 0 0 0 0 0;
           1 0 1 0 1 0 1 2 1 2 1;
            0 2 0 2 0 2 0 0 0 0;
             1 0 1 0 1 0 1 2 1;
              0 2 0 2 0 2 0 0;
               1 0 1 0 1 0 1'
  []
  # Central void region
  [air_center]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    num_sectors_per_side= '4 4 4 4 4 4'
    meshes_to_adapt_to = 'Patterned Patterned Patterned Patterned Patterned Patterned'
    sides_to_adapt = '0 1 2 3 4 5'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = '600 601'
  []
  # Dummy assembly to help form a regular hexagonal core pattern
  # This will be deleted later
  [dummy]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    num_sectors_per_side= '4 4 4 4 4 4'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = '700 701'
       # external_boundary_id = 9998
  []
  # Stitching assemblies together to form the core mesh
  [core]
    type = PatternedHexMeshGenerator
    inputs = 'Patterned cd_0 cd_1 cd_2 cd_3 cd_4 cd_5 cd_6 cd_7 cd_8 cd_9 cd_10 cd_11 ref_0 ref_1 ref_2 ref_3 ref_4 ref_5 dummy air_center'
    # Pattern ID  0       1    2    3    4    5    6    7    8    9    10    11    12    13    14    15    16    17    18   19       20
    pattern_boundary = none
    generate_core_metadata = true
    pattern = '19 17 12 16 19;
             17 11  0  0  1 16;
           10  0  0  0  0  0  2;
          18 0  0   0  0  0  0 15;
        19  9  0  0  20  0  0  3 19;
          18 0  0   0  0  0  0 15;
            8  0  0  0  0  0  4;
             13  7  0  0  5 14;
               19 13  6 14 19'
    rotate_angle = 60
    assign_control_drum_id = true
  []
  # Delete the dummy assemblies
  [del_dummy]
    type = BlockDeletionGenerator
    block = '700 701'
    input = core
    new_boundary = 10000
  []
  # Add peripheral circular region
  [add_outer_shield]
    type = PeripheralRingMeshGenerator
    input = del_dummy
    peripheral_layer_num = 1
    peripheral_ring_radius = 115.0
    input_mesh_external_boundary = 10000
    peripheral_ring_block_id = 250
    peripheral_ring_block_name = outer_shield
 []
 # trim the full core to 1/6 core
 [del_1]
    type = PlaneDeletionGenerator
    point = '0 0 0'
    normal = '10 17.32 0'
    input = add_outer_shield
    new_boundary = 147
  []
  [del_2]
    type = PlaneDeletionGenerator
    point = '0 0 0'
    normal = '10 -17.32 0'
    input = del_1
    new_boundary = 147
  []
  # extrusion
  [extrude]
     type = AdvancedExtruderGenerator
     input = del_2
     heights = '20 160 20'
     # Use `num_layers = '6 16 6'` for BISON mesh
     num_layers = '1 8 1'
     # biased upper and lower reflector mesh, uncomment for BISON mesh
     # biases = '1.6 1.0 0.625'
     direction = '0 0 1'
     top_boundary = 2000
     bottom_boundary = 3000
  []
  # Define some special reflector blocks
  [reflector_bottom_quad]
    type = ParsedSubdomainMeshGenerator
    input = extrude
    combinatorial_geometry = 'z<=20'
    block_id = 1000
    excluded_subdomains = '103 203 303 250 400 401 500 501 502 503 504 600 601'
  []
  [reflector_bottom_tri]
    type = ParsedSubdomainMeshGenerator
    input = reflector_bottom_quad
    combinatorial_geometry = 'z<=20'
    block_id = 1003
    excluded_subdomains = '1000 250 400 401 500 501 502 503 504 600 601'
  []
  [reflector_top_quad]
    type = ParsedSubdomainMeshGenerator
    input = reflector_bottom_tri
    combinatorial_geometry = 'z>=180'
    block_id = 1000
    excluded_subdomains = '103 203 303 200 201 250 400 401 500 501 502 503 504 600 601'
  []
  [reflector_top_tri]
    type = ParsedSubdomainMeshGenerator
    input = reflector_top_quad
    combinatorial_geometry = 'z>=180'
    block_id = 1003
    excluded_subdomains = '1000 200 201 203 250 400 401 500 501 502 503 504 600 601'
  []
  # Assgin block names
  [rename_blocks]
    type = RenameBlockGenerator
    old_block = '     101       100         103          201        200           203            301       303        400            401         10      503      600         601        504          500                501          502         1000           1003'
    new_block = ' mod_ss moderator_quad moderator_tri hp_ss  heat_pipes_quad heat_pipes_tri fuel_quad fuel_tri reflector_tri reflector_quad monolith B4C air_gap_tri air_gap_quad air_gap_quad reflector_tri  reflector_quad reflector_quad reflector_quad reflector_tri'
    input = reflector_top_tri
  []
  # Assign boundary names
  [rename_boundaries]
    type = RenameBoundaryGenerator
    input = rename_blocks
    old_boundary = ' 10000 2000 3000 147'
    new_boundary = ' side top bottom side_mirror'
  []
  # Scale to m unit for general MOOSE application
  # Note that is_meter needs to be true in Griffin
  [scale]
    type = TransformGenerator
    input = rename_boundaries
    transform = SCALE
    vector_value = '1e-2 1e-2 1e-2'
  []
  # Rotate to the desired orientation for simulation
  [rotate_end]
    type = TransformGenerator
    input = scale
    transform = ROTATE
    vector_value = '0 0 -120'
  []
[]

Figure 1: A cartoon showing the step-by-step procedures to generate the HP-MR mesh.

A 1/6 HP-MR core was generated using the input file shown above with steps illustrated in Figure 1. This mesh contains stainless steel envelopes for moderators and heat pipes, while the helium gaps are not meshed. The mesh density in radial direction is high, as multiple small features (fuel rods, moderators, heat pipes and control drums) are involved.

Compared with the mesh used in the legacy model, which was created from Cubit:

An outer surface zone surrounding the reflectors was added to provide the appropriate convex boundary for neutronic calculation.
The mesh was adjusted to preserve the actual physical volume of each material region.

The same mesh used by Griffin could also be adopted by BISON. However, the relatively coarse Griffin mesh would lead to power balance issues. To be specific, near the external boundaries (i.e., top/bottom/peripheral surfaces) as well as the interfaces between the reactor matrix and heat pipes, the mesh is not fine enough to ensure accurate calculation of heat flux, which is proportional to the temperature gradient. This issue is especially important to the interfaces between the reactor matrix and heat pipe outer surfaces, because they are the major heat sink of the system. Therefore, for the BISON mesh, the boundary layer and biased meshing features available in the Reactor Module are used to create a finer mesh with a focus on these external boundaries and interfaces. Only minor modifications on the mesh generation input file are needed to achieve these features, as indicated in the comments in the Griffin mesh generation input file shown above. The differences between the Griffin and the BISON meshes are illustrated in Figure 2.

Figure 2: The difference between the meshes used by the Griffin model (left) and the BISON model (right).

A more detailed comparison can be found in Table 1, which shows that using a finer mesh for BISON significantly reduced the power imbalance issue.

Table 1: Percentages of power removed from different surfaces of the HPMR showing i the improvements made by using a finer mesh for BISON modeling.

	Total Energy Removed	Heat Pipes	Top/Bottom Surfaces	Symmetric Boundary	Peripheral Surface
Coarse Mesh for BISON, standalone BISON with flat power profile	93.31%	90.04%	2.84%	0.40%	0.02%
Fine Mesh for BISON, standalone BISON with flat power profile	98.66%	96.05%	2.19%	0.40%	0.03%
Fine Mesh for BISON, steady-state multiphysics simulation	99.30%	97.50%	1.37%	0.41%	0.03%

The input file used to generate the fine BISON mesh is listed as follows. Note the boundary layers and biasing setup in the input file. Additionally, the heat pipe related blocks are deleted from the BISON mesh and the heat pipe surfaces are defined as boundaries.

#####################################################################
# This is the MOOSE input file to generate the mesh that can be
# used for the 1/6 core Heat-Pipe Microreactor Multiphysics
# simulations (BISON thermal simulation).
# Running this input requires MOOSE Reactor Module Objects
# Users should use
# --mesh-only HPMR_OneSixth_Core_meshgenerator_tri_rotate_bdry_fine.e
# command line argument to generate
# exodus file for further Multiphysics simulations.
#####################################################################
[GlobalParams]
  # we use ':=' to override the values already defined in the !include file
  create_outward_interface_boundaries := false
[]

[Mesh]
  [HP_hex]
    # inner background boundary layer, only for BISON mesh
    background_inner_boundary_layer_bias = 1.5
    background_inner_boundary_layer_intervals = 3
    background_inner_boundary_layer_width = 0.03
  []
  [add_outer_shield]
    peripheral_layer_num := 2
    # Peripheral boundary layer, only for BISON mesh
    peripheral_outer_boundary_layer_bias = 0.625
    peripheral_outer_boundary_layer_intervals = 3
    peripheral_outer_boundary_layer_width = 2
  []
  [extrude]
    num_layers := '6 16 6'
    # biased upper and lower reflector mesh, only for BISON mesh
    biases = '1.6 1.0 0.625'
  []
  [split_hp_ss]
    type = SubdomainBoundingBoxGenerator
    input = rotate_end
    block_id = 1201
    block_name = 'hp_ss_up'
    restricted_subdomains = 'hp_ss'
    bottom_left = '-100 -100 1.8'
    top_right = '100 100 3.0'
  []
  [add_exterior_ht_low]
    type = SideSetsBetweenSubdomainsGenerator
    input = split_hp_ss
    paired_block = 'hp_ss'
    primary_block = 'monolith'
    new_boundary = 'heat_pipe_ht_surf_low'
  []
  [add_exterior_ht_up]
    type = SideSetsBetweenSubdomainsGenerator
    input = add_exterior_ht_low
    paired_block = 'hp_ss_up'
    primary_block = 'reflector_quad'
    new_boundary = 'heat_pipe_ht_surf_up'
  []
  [add_exterior_bot]
    type = SideSetsBetweenSubdomainsGenerator
    input = add_exterior_ht_up
    paired_block = 'heat_pipes_quad heat_pipes_tri hp_ss'
    primary_block = 'reflector_quad reflector_tri'
    new_boundary = 'heat_pipe_ht_surf_bot'
  []
  [merge_hp_surf]
    type = RenameBoundaryGenerator
    input = add_exterior_bot
    old_boundary = 'heat_pipe_ht_surf_low heat_pipe_ht_surf_up'
    new_boundary = 'heat_pipe_ht_surf heat_pipe_ht_surf'
  []
  [remove_hp]
    type = BlockDeletionGenerator
    input = merge_hp_surf
    block = 'hp_ss hp_ss_up heat_pipes_quad heat_pipes_tri'
  []

  # remove extra nodesets to limit the size of the mesh
  [clean_up]
    type = BoundaryDeletionGenerator
    input = remove_hp
    boundary_names = '1 3'
  []
[]

# This input defines the mesh generation for the Griffin mesh
!include HPMR_OneSixth_Core_meshgenerator_tri.i

Multiphysics Model Setup

Cross-Section Generation Using Serpent Code

The first step is to generate homogenized multi-group cross-sections using Serpent-2. The generated cross-sections are then converted into an XML-format file for compatibility with Griffin.

%----------------------------------------------------------------------------------------------------
% ANL HPMR Serpent (Version 2.1.32) Model Input
% If using or referring to this model, please cite as explained in
% https://mooseframework.inl.gov/virtual_test_bed/citing.html
%Contacts: Ahmed Abdelhameed (aabdelhameed.at.anl.gov), Nicolas Stauff (nstauff.at.anl.gov)
% -----------------------------------------------------------------------------
set title "ANL MicroReactor Core - pf - PF = 40"

pbed 801 33 "TRISO_U900_PF40_R100.inp"

% --- Compact --------------------------------------------------------------------------------------
surf  1 hexyprism 0.0 0.0 13.376      0 200 %fueled core
surf 91 hexyprism 0.0 0.0 130.0       0 200 %fueled core

surf 95l pz  20
surf 95u pz  180
surf 96l pz  18.0
surf 96u pz  182.0

surf  2 cyl 0.0 0.0 1.00   0.0 200
surf  3 cyl 0.0 0.0 0.825  0.0 200
surf 35 cyl 0.0 0.0 0.875  0.0 200
surf 36 cyl 0.0 0.0 0.900  0.0 200
surf 3g cyl 0.0 0.0 0.920  0.0 200
surf 5g cyl 0.0 0.0 1.07   0.0 200
surf  5 cyl 0.0 0.0 1.05   0.0 200
surf 51 cyl 0.0 0.0 0.97   0.0 200
surf 52 cyl 0.0 0.0 0.90   0.0 200
surf 53 cyl 0.0 0.0 0.80   0.0 200

surf inf inf

particle 900
fuel      2.1250e-02
buffer    3.1250e-02
PyC1      3.5250e-02
SiC       3.8750e-02
PyC2      4.2750e-02
matrix_pin

% infinite cells defining material universes
cell 51 802 moderator  -inf
cell 52 33  matrix_pin -inf
cell 53 803    matrix  -inf
cell 55 804       air  -inf
cell 56 805  beryllium -inf
cell 57 806  shell_mod -inf
cell 64 813  shell_hp  -inf
cell 58 820  shell_air_center -inf

% mixed wick and coo_vap and liqu
cell 65 815 hp_vp_liq_wick -inf

% mixed air and ss in moderator
cell 68 816 shell_air_mod -inf

% mixed air and ss in heatpipe
cell 72 817 shell_air_hp   -inf

%  control drums
cell 61 810 beryllium_drum -inf
cell 62 811   B4C_drum -inf
cell 63 812   B4C_central -inf

% fuel
cell 11   3 fill 801        -2  95l -95u % fuel compact
cell 12   3 fill 803         2  95l -95u
cell 125  3 fill 805                 95u
cell 126  3 fill 805                -95l

% Yan simplified model removed caps for the HP and moderator pins extended into the Be reflectors
% moderator
cell 13   2 fill 802    -3      95l -95u % active length
cell 14   2 fill 816     3 -3g  95l -95u
cell 15   2 fill 803     3g     95l -95u
cell 151  2 fill 805       -36  96l -95l % 0.05 cm shell
cell 152  2 fill 805       -36  95u -96u % 0.05 cm shell
cell 151a 2 fill 805   -3g  36  96l -95l % 0.05 cm shell
cell 152a 2 fill 805   -3g  36  95u -96u % 0.05 cm shell
cell 153  2 fill 805        3g  96l -95l % 0.05 cm shell
cell 154  2 fill 805        3g  95u -96u % 0.05 cm shell
cell 155  2 fill 805                 96u
cell 156  2 fill 805                -96l

% heat pipe cell
cell 16i  1 fill 815      -51 95l
cell 16o  1 fill 817   51 -5g 95l
cell 17   1 fill 803       5g 95l -95u
cell 175  1 fill 805       5g 95u
cell 176  1 fill 805         -96l
cell 177  1 fill 805      -5  96l -95l
cell 178  1 fill 805       5g 96l -95l
cell 178a 1 fill 805   -5g 5  96l -95l

% monolith filled cell
% cell 18  8 fill 804     -inf
cell 18  8 fill 820     -inf

% monolith filled cell
cell 19  9 fill 803     -inf

% monolith filled cell
cell 20 10 fill 805           -inf

% monolith filled cell
% cell 21 11 fill 804           -inf
cell 21 11 fill 820           -inf

% --- Assembly lattice

lat 20  2  0.0 0.0 19 19 2.3
% New Assembly Configuration - 2 time less moderator - more fuel
 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
  9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
   9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
    9 9 9 9 9 9 9 9 9 1 3 1 3 1 3 1 9 9 9
     9 9 9 9 9 9 9 9 3 3 2 3 2 3 2 3 9 9 9
      9 9 9 9 9 9 9 1 2 1 3 1 3 1 3 1 9 9 9
       9 9 9 9 9 9 3 3 3 3 2 3 2 3 2 3 9 9 9
        9 9 9 9 9 1 2 1 2 1 3 1 3 1 3 1 9 9 9
         9 9 9 9 3 3 3 3 3 3 2 3 2 3 2 3 9 9 9
          9 9 9 1 2 1 2 1 2 1 3 1 3 1 3 1 9 9 9
           9 9 9 3 3 3 3 3 3 2 3 2 3 2 3 9 9 9 9
            9 9 9 1 2 1 2 1 3 1 3 1 3 1 9 9 9 9 9
             9 9 9 3 3 3 3 2 3 2 3 2 3 9 9 9 9 9 9
              9 9 9 1 2 1 3 1 3 1 3 1 9 9 9 9 9 9 9
               9 9 9 3 3 2 3 2 3 2 3 9 9 9 9 9 9 9 9
                9 9 9 1 3 1 3 1 3 1 9 9 9 9 9 9 9 9 9
                 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
                  9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
                   9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9

trans 20 0.0 0.0 0.0 0.0 0.0 30.

lat 80  2  0.0 0.0 1 1 26.752
 10  % Be

lat 90  2  0.0 0.0 1 1 26.752
  8  % air

lat 91  2  0.0 0.0 1 1 26.752
  11  % air o  B4C

lat 50  2  0.0 0.0 13 13 26.752
 90 90 90 90 90 90 90 90 90 90 90 90 90
  90 90 90 90 90 90 90 90 90 90 90 90 90
   90 90 90 90 90 90 90 80 90 80 90 90 90
    90 90 90 90 90 80 90 20 20 90 80 90 90
     90 90 90 90 90 20 20 20 20 20 90 90 90
      90 90 90 80 20 20 20 20 20 20 80 90 90
       90 90 90 90 20 20 91 20 20 90 90 90 90
        90 90 80 20 20 20 20 20 20 80 90 90 90
         90 90 90 20 20 20 20 20 90 90 90 90 90
          90 90 80 90 20 20 90 80 90 90 90 90 90
           90 90 90 80 90 80 90 90 90 90 90 90 90
            90 90 90 90 90 90 90 90 90 90 90 90 90
             90 90 90 90 90 90 90 90 90 90 90 90 90

% -----------DRUMS DEFINITTION---------------------------------------------------------------------------------------

cell 102 0 fill 50 -91 731 732 733 734 735 736 737 738 739 740 741 742
cell 180 0 fill       731 -731
cell 181 0 fill       732 -732
cell 182 0 fill       733 -733
cell 183 0 fill       734 -734
cell 184 0 fill       735 -735
cell 185 0 fill       736 -736
cell 186 0 fill       737 -737
cell 187 0 fill       738 -738
cell 188 0 fill       739 -739
cell 189 0 fill       740 -740
cell 190 0 fill       741 -741
cell 191 0 fill       742 -742

surf 731 cyl 80.2560   0.00 13.2500   0.00 200.00
surf 831 cyl 80.2560   0.00 12.2500   0.00 200.00
surf 771 plane 1.0000   0.00 0.0  88.92
surf 871 plane 1.0000  -0.00 0.0  71.59
cell 801  731 fill 810   -731   -771
cell 802  731 fill 811    831  -731  771
cell 803  731 fill 810   -831   771
surf 732 cyl 40.1280  69.50 13.2500   0.00 200.00
surf 832 cyl 40.1280  69.50 12.2500   0.00 200.00
surf 772 plane 1.0000   1.73 0.0 177.84
surf 872 plane 1.0000   1.73 0.0 143.19
cell 806  732 fill 810   -732   -772
cell 807  732 fill 811    832  -732  772
cell 808  732 fill 810   -832   772
surf 733 cyl -40.1280  69.50 13.2500   0.00 200.00
surf 833 cyl -40.1280  69.50 12.2500   0.00 200.00
surf 773 plane 1.0000  -1.73 0.0 -177.84
surf 873 plane 1.0000  -1.73 0.0 -143.19
cell 811  733 fill 810   -733   773
cell 812  733 fill 811    833  -733  -773
cell 813  733 fill 810   -833   -773
surf 734 cyl -80.2560   0.00 13.2500   0.00 200.00
surf 834 cyl -80.2560   0.00 12.2500   0.00 200.00
surf 774 plane 1.0000  -0.00 0.0 -88.92
surf 874 plane 1.0000   0.00 0.0 -71.59
cell 816  734 fill 810   -734   774
cell 817  734 fill 811    834  -734  -774
cell 818  734 fill 810   -834   -774
surf 735 cyl -40.1280 -69.50 13.2500   0.00 200.00
surf 835 cyl -40.1280 -69.50 12.2500   0.00 200.00
surf 775 plane 1.0000   1.73 0.0 -177.84
surf 875 plane 1.0000   1.73 0.0 -143.19
cell 821  735 fill 810   -735   775
cell 822  735 fill 811    835  -735  -775
cell 823  735 fill 810   -835   -775
surf 736 cyl 40.1280 -69.50 13.2500   0.00 200.00
surf 836 cyl 40.1280 -69.50 12.2500   0.00 200.00
surf 776 plane 1.0000  -1.73 0.0 177.84
surf 876 plane 1.0000  -1.73 0.0 143.19
cell 826  736 fill 810   -736   -776
cell 827  736 fill 811    836  -736  776
cell 828  736 fill 810   -836   776
surf 737 cyl 80.2563  46.34 13.2500   0.00 200.00
surf 837 cyl 80.2563  46.34 12.2500   0.00 200.00
surf 777 plane 1.0000   0.58 0.0 117.01
surf 877 plane 1.0000   0.58 0.0  97.01
cell 831  737 fill 810   -737   -777
cell 832  737 fill 811    837  -737  777
cell 833  737 fill 810   -837   777
surf 738 cyl 0.0000  92.67 13.2500   0.00 200.00
surf 838 cyl 0.0000  92.67 12.2500   0.00 200.00
surf 778 plane 0.0000   1.00 0.0 101.33
surf 878 plane 0.0000   1.00 0.0  84.01
cell 836  738 fill 810   -738   -778
cell 837  738 fill 811    838  -738  778
cell 838  738 fill 810   -838   778
surf 739 cyl -80.2563  46.34 13.2500   0.00 200.00
surf 839 cyl -80.2563  46.34 12.2500   0.00 200.00
surf 779 plane 1.0000  -0.58 0.0 -117.01
surf 879 plane 1.0000  -0.58 0.0 -97.01
cell 841  739 fill 810   -739   779
cell 842  739 fill 811    839  -739  -779
cell 843  739 fill 810   -839   -779
surf 740 cyl -80.2563 -46.34 13.2500   0.00 200.00
surf 840 cyl -80.2563 -46.34 12.2500   0.00 200.00
surf 780 plane 1.0000   0.58 0.0 -117.01
surf 880 plane 1.0000   0.58 0.0 -97.01
cell 846  740 fill 810   -740   780
cell 847  740 fill 811    840  -740  -780
cell 848  740 fill 810   -840   -780
surf 741 cyl -0.0000 -92.67 13.2500   0.00 200.00
surf 841 cyl -0.0000 -92.67 12.2500   0.00 200.00
surf 781 plane 0.0000   1.00 0.0 -101.33
surf 881 plane 1.0000 1219076296694385.50 0.0 -102414528894101376.00
cell 851  741 fill 810   -741   781
cell 852  741 fill 811    841  -741  -781
cell 853  741 fill 810   -841   -781
surf 742 cyl 80.2563 -46.34 13.2500   0.00 200.00
surf 842 cyl 80.2563 -46.34 12.2500   0.00 200.00
surf 782 plane 1.0000  -0.58 0.0 117.01
surf 882 plane 1.0000  -0.58 0.0  97.01
cell 856  742 fill 810   -742   -782
cell 857  742 fill 811    842  -742  782
cell 858  742 fill 810   -842   782

% --------------------------------------------------------------------------------------------------
cell 104 0 outside  91

set ures 1

% --------------------------------------------------------------------------------------------------
% Material data:

% materials

mat air -0.18e-3   rgb 255 255 255   tmp  700
 2004.01c  1

mat shell_mod -7.90   rgb 133 133 133    tmp  700
 6012.01c  1.9010E-03
 14028.01c  9.2693E-03
 14029.01c  4.7251E-04
 14030.01c  3.1166E-04
 15031.01c  4.1322E-04
 16032.01c  2.4710E-04
 16033.01c  1.9511E-06
 16034.01c  1.1056E-05
 16036.01c  3.0016E-08
 24050.01c  7.9116E-03
 24052.01c  1.5257E-01
 24053.01c  1.7300E-02
 24054.01c  4.3063E-03
 25055.01c  1.0280E-02
 26054.01c  3.9029E-02
 26056.01c  6.1213E-01
 26057.01c  1.4144E-02
 26058.01c  1.3343E-03
 28058.01c  7.7516E-02
 28060.01c  2.9859E-02
 28061.01c  1.2981E-03
 28062.01c  4.1389E-03
 28064.01c  1.0544E-03
 42092.01c  2.1259E-03
 42094.01c  1.3299E-03
 42095.01c  2.3030E-03
 42096.01c  2.4191E-03
 42097.01c  1.3903E-03
 42098.01c  3.5249E-03
42100.01c  1.4135E-03

mat shell_hp -7.90   rgb 133 133 133    tmp  700
 6012.01c  1.9010E-03
14028.01c  9.2693E-03
14029.01c  4.7251E-04
14030.01c  3.1166E-04
15031.01c  4.1322E-04
16032.01c  2.4710E-04
16033.01c  1.9511E-06
16034.01c  1.1056E-05
16036.01c  3.0016E-08
24050.01c  7.9116E-03
24052.01c  1.5257E-01
24053.01c  1.7300E-02
24054.01c  4.3063E-03
25055.01c  1.0280E-02
26054.01c  3.9029E-02
26056.01c  6.1213E-01
26057.01c  1.4144E-02
26058.01c  1.3343E-03
28058.01c  7.7516E-02
28060.01c  2.9859E-02
28061.01c  1.2981E-03
28062.01c  4.1389E-03
28064.01c  1.0544E-03
42092.01c  2.1259E-03
42094.01c  1.3299E-03
42095.01c  2.3030E-03
42096.01c  2.4191E-03
42097.01c  1.3903E-03
42098.01c  3.5249E-03
42100.01c  1.4135E-03

mat shell_air_mod 2.290E-02 rgb 133 133 133    tmp  700
   2004.01c   1.983E-05
   6012.01c   4.349E-05
   14028.01c   2.121E-04
   14029.01c   1.081E-05
   14030.01c   7.130E-06
   15031.01c   9.454E-06
   16032.01c   5.653E-06
   16033.01c   4.464E-08
   16034.01c   2.530E-07
   16036.01c   6.867E-10
   24050.01c   1.810E-04
   24052.01c   3.491E-03
   24053.01c   3.958E-04
   24054.01c   9.852E-05
   25055.01c   2.352E-04
   26054.01c   8.929E-04
   26056.01c   1.400E-02
   26057.01c   3.236E-04
   26058.01c   3.053E-05
   28058.01c   1.773E-03
   28060.01c   6.831E-04
   28061.01c   2.970E-05
   28062.01c   9.469E-05
   28064.01c   2.412E-05
   42092.01c   4.864E-05
   42094.01c   3.043E-05
   42095.01c   5.269E-05
   42096.01c   5.535E-05
   42097.01c   3.181E-05
   42098.01c   8.065E-05
   42100.01c   3.234E-05

mat shell_air_hp 6.771E-02 rgb 133 133 133    tmp  700
2004.01c  5.629E-06
6012.01c  1.287E-04
14028.01c  6.276E-04
14029.01c  3.199E-05
14030.01c  2.110E-05
15031.01c  2.798E-05
16032.01c  1.673E-05
16033.01c  1.321E-07
16034.01c  7.486E-07
16036.01c  2.032E-09
24050.01c  5.357E-04
24052.01c  1.033E-02
24053.01c  1.171E-03
24054.01c  2.916E-04
25055.01c  6.960E-04
26054.01c  2.643E-03
26056.01c  4.145E-02
26057.01c  9.576E-04
26058.01c  9.034E-05
28058.01c  5.248E-03
28060.01c  2.022E-03
28061.01c  8.789E-05
28062.01c  2.802E-04
28064.01c  7.139E-05
42092.01c  1.439E-04
42094.01c  9.004E-05
42095.01c  1.559E-04
42096.01c  1.638E-04
42097.01c  9.413E-05
42098.01c  2.387E-04
42100.01c  9.570E-05

mat shell_air_center 8.815E-04 rgb 133 133 133  tmp  700
2004.01c   2.68E-05
6012.01c   1.62E-06
14028.01c   7.92E-06
14029.01c   4.04E-07
14030.01c   2.66E-07
15031.01c   3.53E-07
16032.01c   2.11E-07
16033.01c   1.67E-09
16034.01c   9.45E-09
16036.01c   2.57E-11
24050.01c   6.76E-06
24052.01c   1.30E-04
24053.01c   1.48E-05
24054.01c   3.68E-06
25055.01c   8.79E-06
26054.01c   3.34E-05
26056.01c   5.23E-04
26057.01c   1.21E-05
26058.01c   1.14E-06
28058.01c   6.63E-05
28060.01c   2.55E-05
28061.01c   1.11E-06
28062.01c   3.54E-06
28064.01c   9.01E-07
42092.01c   1.82E-06
42094.01c   1.14E-06
42095.01c   1.97E-06
42096.01c   2.07E-06
42097.01c   1.19E-06
42098.01c   3.01E-06
42100.01c   1.21E-06

mat coo_vap  -1.11e-4 rgb 250 50 250    tmp  700
19039.01c   0.93258
19040.01c   0.00012
19041.01c   0.06730

mat coo_liq  -0.705 rgb 50 250 50    tmp  700
19039.01c   0.93258
19040.01c   0.00012
19041.01c   0.06730

mat wick  -2.753 rgb 250 250 50     tmp  700
 6012.01c  5.589E-04
14028.01c  2.725E-03
14029.01c  1.389E-04
14030.01c  9.163E-05
15031.01c  1.215E-04
16032.01c  7.265E-05
16033.01c  5.736E-07
16034.01c  3.250E-06
16036.01c  8.825E-09
24050.01c  2.326E-03
24052.01c  4.485E-02
24053.01c  5.086E-03
24054.01c  1.266E-03
25055.01c  3.022E-03
26054.01c  1.147E-02
26056.01c  1.800E-01
26057.01c  4.158E-03
26058.01c  3.923E-04
28058.01c  2.279E-02
28060.01c  8.778E-03
28061.01c  3.816E-04
28062.01c  1.217E-03
28064.01c  3.100E-04
42092.01c  6.250E-04
42094.01c  3.910E-04
42095.01c  6.771E-04
42096.01c  7.112E-04
42097.01c  4.087E-04
42098.01c  1.036E-03
42100.01c  4.156E-04
19039.01c  6.584E-01
19040.01c  8.472E-05
19041.01c  4.751E-02

mat hp_vp_liq_wick  8.324E-03  rgb 250 250 50     tmp  700
   6012.01c   3.808E-06
   14028.01c   1.856E-05
   14029.01c   9.463E-07
   14030.01c   6.242E-07
   15031.01c   8.277E-07
   16032.01c   4.949E-07
   16033.01c   3.908E-09
   16034.01c   2.214E-08
   16036.01c   6.012E-11
   24050.01c   1.585E-05
   24052.01c   3.055E-04
   24053.01c   3.465E-05
   24054.01c   8.625E-06
   25055.01c   2.059E-05
   26054.01c   7.814E-05
   26056.01c   1.226E-03
   26057.01c   2.833E-05
   26058.01c   2.673E-06
   28058.01c   1.553E-04
   28060.01c   5.980E-05
   28061.01c   2.600E-06
   28062.01c   8.291E-06
   28064.01c   2.112E-06
   42092.01c   4.258E-06
   42094.01c   2.664E-06
   42095.01c   4.613E-06
   42096.01c   4.845E-06
   42097.01c   2.784E-06
   42098.01c   7.058E-06
   42100.01c   2.831E-06
   19039.01c   5.895E-03
   19040.01c   7.586E-07
   19041.01c   4.254E-04

therm h-yh2 h-yh2.84t
therm y-yh2 y-yh2.84t
mat moderator   -4.0850  moder h-yh2 1001 moder y-yh2 39089  tmp  700.0   rgb 73 64 171
39089.01c  0.357142857
1001.01c   0.642857143

% -- Fuel

div fuel sep 0
mat fuel   -10.7440 burn 1    rgb 255 0 0  tmp  700
  92235.01c    6.8794e-02
  92238.01c    2.7604e-01
  6012.01c    1.3793e-01
  8016.01c    5.1724e-01

% --- Carbon buffer layer:

mat buffer   -1.0400   rgb 0 200 200  moder grph_f 6012    tmp  700
 6012.01c    1.0

% --- Pyrolytic carbon layer:

mat PyC1     -1.8820    rgb 0 200 200  moder grph_f 6012    tmp  700
 6012.01c    1.0

% --- Pyrolytic carbon layer:

mat PyC2     -1.8820  rgb 0 200 200    moder grph_f 6012     tmp  700
 6012.01c    1.0

% --- Silicon carbide layer:

mat SiC     -3.1710  rgb 0 100 0     tmp  700
14028.01c    0.4611
14029.01c    0.0234
14030.01c    0.0154
 6012.01c    0.5

% --- Graphite matrix:
therm grph_f grph.84t
mat matrix_pin   -1.8060  rgb 220 220 220    moder grph_f 6012    tmp  700
 6012.01c  0.9999997
 5010.01c  0.0000003

therm grph_m grph.84t
mat matrix     -1.8060 rgb 200 200 200    moder grph_m 6012    tmp  700
 6012.01c  0.9999997
 5010.01c  0.0000003

therm bemet be-met.84t
mat beryllium  -1.848   rgb 0 255 0 moder bemet 4009  tmp  700
 4009.01c  1.00

mat beryllium_drum  -1.848   rgb 50 220 50 moder bemet 4009  tmp  700
 4009.01c  1.00

%
mat B4C_drum  -2.510  rgb 250 0 0   tmp  700
 5010.01c 0.76
 5011.01c 0.04
 6012.01c 0.2
%
mat B4C_central  -1.25  rgb 150 20 20   tmp  700
 5010.01c 0.76
 5011.01c 0.04
 6012.01c 0.2
% -----------------------------------------------------------------------------

% Calculation parameters:

% -----------------------------------------------------------------------------

% --- Geometry plot:
plot 3 10000 10000   21.0
plot 2 10000 10000    5.0
mesh 3 5000 5000
mesh 2 5000 5000

set mcvol 100000000
set xenon 1

% --- Libraries:
% Add your own libraries here
set acelib "sss_endfb8_serpent_yan.xsdir"  % ACE library for cross-sections
set declib "sss_endfb7.dec" % Decay library
set nfylib "sss_endfb7.nfy" % Neutron fission yield library

% --- Boundary conditions:

set bc 1 1 1

% --- Power density

set power 2e6  % W

% --- Histories:
set pop 100000 500 30

set gcu 801 802 803 805 815 816 817 810 811 820

set nfg  11  8.00000e-08 1.80000e-07 6.25000e-07 1.30000e-06 4.00000e-06
             1.48730e-04 9.11800e-03 1.83000e-01 5.00000e-01 1.35300e+00
% -----------------------------------------------------------------------------

Note on LFS

Large File Storage (LFS) is used for the following file: TRISO_U900_PF40_R100

This file defines the distribution, and radius of the TRISO particles to achieve 40% packing fraction. Make sure to install git lfs, then fetch and pull these files to be able to run the Serpent inputs.

The Monte Carlo Serpent-2 model has been set up to simulate the double heterogeneity of the HP-MR full core explicitly with the most recent ENDF/B-8.0 cross section data sets. This Monte Carlo model is used not only to provide reference solutions to the Griffin neutronics model for multiphysics simulations, but also to generate condensed multigroup cross sections for Griffin.

In particular, multigroup cross sections were generated for each material zone by tallying the average reaction rates and group fluxes within that material region in Serpent-2 criticality calculations. Different sets of cross sections were prepared with fuel and moderator temperatures ( $var element = document.getElementById("moose-equation-8d18655a-e041-45a9-b658-5f3d43e6b64f");katex.render("T_f", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-571d3d3b-b417-472c-8759-077ba0aefd1c");katex.render("T_m", element, {displayMode:false,throwOnError:false});$ ) assumed at 600 K, 700 K, 800 K, 1000 K, and 1200 K respectively.

These cross-section sets were then read by the ISOXML utility module in Griffin to create a bi-dimensional cross-section table, which was used in the multiphysics coupled transient simulations.

Griffin Model

Griffin is used to govern the neutronics of the HP-MR as the parent application of the multiphysics simulation.

The DFEM-SN(1,3) neutronics solver with CMFD acceleration in Griffin is used in this model.

################################################################################
## NEAMS Micro-Reactor Application Driver                                     ##
## Heat Pipe Microreactor Steady State                                        ##
## Griffin Main Application input file                                        ##
## DFEM-SN (1, 3) with CMFD acceleration                                      ##
################################################################################

[Mesh]
  [fmg]
    type = FileMeshGenerator
    file = '../mesh/gold/HPMR_OneSixth_Core_meshgenerator_tri_rotate_bdry.e'
  []
  [fmg_id]
    type = SubdomainExtraElementIDGenerator
    input = fmg
    subdomains = '200 203 100 103 301 303 10 503 600 601 201 101 400 401 250'
    extra_element_id_names = 'material_id equivalence_id'
    extra_element_ids = '815 815 802 802 801 801 803 811 820 820 817 816 805 805 820;
                         815 815 802 802 801 801 803 811 820 820 817 816 805 805 820'
  []
  [coarse_mesh]
    type = GeneratedMeshGenerator
    dim = 3
    nx = 10
    ny = 10
    nz = 10
    xmin = -0.1
    xmax = 1.1
    ymin = -0.1
    ymax = 1.2
    zmin = -0.0
    zmax = 2.1
  []
  [assign_coarse_id]
    type = CoarseMeshExtraElementIDGenerator
    input = fmg_id
    coarse_mesh = coarse_mesh
    extra_element_id_name = coarse_element_id
  []
  uniform_refine = 0
[]

[Executioner]
  type = SweepUpdate

  richardson_abs_tol = 1e-6 #-8
  richardson_rel_tol = 1e-8
  richardson_value = eigenvalue

  richardson_max_its = 1000
  inner_solve_type = GMRes
  max_inner_its = 20

  fixed_point_max_its = 1
  custom_pp = eigenvalue
  custom_rel_tol = 1e-6
  force_fixed_point_solve = true

  cmfd_acceleration = true
  coarse_element_id = coarse_element_id
  prolongation_type = multiplicative
  max_diffusion_coefficient = 1
[]

[TransportSystems]
  particle = neutron
  equation_type = eigenvalue

  G = 11
  VacuumBoundary = '10000 2000 3000'
  ReflectingBoundary = '147'

  [sn]
    scheme = DFEM-SN
    family = MONOMIAL
    order = first
    AQtype = Gauss-Chebyshev
    NPolar = 1
    NAzmthl = 3
    NA = 2
    sweep_type = asynchronous_parallel_sweeper
    using_array_variable = true
    collapse_scattering = true
    n_delay_groups = 6
  []
[]

[AuxVariables]
  [Tf]
    initial_condition = 873.15
    order = CONSTANT
    family = MONOMIAL
  []
  [Tm]
    initial_condition = 873.15
    order = CONSTANT
    family = MONOMIAL
  []
[]

[GlobalParams]
  library_file = '../isoxml/fullcore_xml_G11_endfb8_ss_tr.xml'
  library_name = fullcore_xml_G11_endfb8_ss_tr
  isotopes = 'pseudo'
  densities = 1.0
  is_meter = true
  plus = true
  dbgmat = false
  grid_names = 'Tfuel Tmod'
  grid_variables = 'Tf Tm'
[]

[Materials]
  [mod]
    type = CoupledFeedbackMatIDNeutronicsMaterial
    block = '200 203 100 103 301 303 10 503 600 601 201 101 400 401 250'
  []
[]

[PowerDensity]
  power = 345.6e3
  power_density_variable = power_density
  integrated_power_postprocessor = integrated_power
[]

[MultiApps]
  [bison]
    type = FullSolveMultiApp
    app_type = BisonApp
    positions = '0 0 0'
    input_files = HPMR_thermo_ss.i
    execute_on = 'timestep_end'
    keep_solution_during_restore = true
  []
[]

[Transfers]
  [to_sub_power_density]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    to_multi_app = bison
    source_variable = power_density
    variable = power_density
    from_postprocessors_to_be_preserved = integrated_power
    to_postprocessors_to_be_preserved = power

    search_value_conflicts = false
  []
  [from_sub_temp_fuel]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = bison
    variable = Tf
    source_variable = Tfuel
    execute_on = 'initial timestep_end'
    displaced_source_mesh = false
    displaced_target_mesh = false
    use_displaced_mesh = false

    search_value_conflicts = false
  []
  [from_sub_temp_mod]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = bison
    variable = Tm
    source_variable = Tmod
    execute_on = 'initial timestep_end'
    displaced_source_mesh = false
    displaced_target_mesh = false
    use_displaced_mesh = false

    search_value_conflicts = false
  []
[]

[UserObjects]
  [ss]
    type = TransportSolutionVectorFile
    transport_system = sn
    writing = true
    execute_on = final
  []
[]

[Postprocessors]
  [scaled_power_avg]
    type = ElementAverageValue
    block = 'fuel_quad fuel_tri'
    variable = power_density
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_avg]
    type = ElementAverageValue
    variable = Tf
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tf
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_min]
    type = ElementExtremeValue
    value_type = min
    variable = Tf
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [mod_temp_avg]
    type = ElementAverageValue
    variable = Tm
    block = 'moderator_quad moderator_tri'
    execute_on = 'initial timestep_end'
  []
  [mod_temp_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tm
    block = 'moderator_quad moderator_tri'
    execute_on = 'initial timestep_end'
  []
  [mod_temp_min]
    type = ElementExtremeValue
    value_type = min
    variable = Tm
    block = 'moderator_quad moderator_tri'
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  csv = true
  exodus = true
  perf_graph = true
[]

BISON Model

BISON is used to simulate thermal physics of the 1/6 HP-MR core as a child application, as it is equipped with a comprehensive set of thermophysical properties of relevant nuclear materials, especially the TRISO-loaded graphite and 316SS.

As the heat pipe behavior is simulated by Sockeye, the heat pipe blocks are removed from the mesh used by BISON. The thermal behavior of the rest of the HP-MR core is dominated by heat conduction. Convective boundary conditions were defined at the top and bottom of the model with an external temperature of 800 K and a heat transfer coefficient (HTC) of 100 W/m²•K.

################################################################################
## NEAMS Micro-Reactor Application Driver                                     ##
## Heat Pipe Microreactor Steady State                                        ##
## BISON Child Application input file                                         ##
## Thermal Only Physics                                                       ##
################################################################################

R_clad_o = 0.0105 # heat pipe outer radius
R_hp_hole = 0.0107 # heat pipe + gap
num_sides = 12 # number of sides of heat pipe as a result of mesh polygonization
alpha = '${fparse 2 * pi / num_sides}'
perimeter_correction = '${fparse 0.5 * alpha / sin(0.5 * alpha)}' # polygonization correction factor for perimeter
area_correction = '${fparse sqrt(alpha / sin(alpha))}' # polygonization correction factor for area
corr_factor = '${fparse R_hp_hole / R_clad_o * area_correction / perimeter_correction}'

[GlobalParams]
  flux_conversion_factor = 1
[]

[Mesh]
  [fmg]
    type = FileMeshGenerator
    file = '../mesh/gold/HPMR_OneSixth_Core_meshgenerator_tri_rotate_bdry_fine.e'
  []
[]

[Variables]
  [temp]
    initial_condition = 800
  []
[]

[Kernels]
  # The time derivative kernel has been removed
  # We are obtaining a steady-state solution and it's preferrable not to march in time
  # to steady state as marching is expensive.
  [heat_conduction]
    type = HeatConduction
    variable = temp
  []
  [heat_source_fuel]
    type = CoupledForce
    variable = temp
    block = 'fuel_quad fuel_tri'
    v = power_density
  []
[]

[AuxVariables]
  [power_density]
    block = 'fuel_quad fuel_tri'
    family = L2_LAGRANGE
    order = FIRST
    initial_condition = 3.4e6
  []
  [Tfuel]
    block = 'fuel_quad fuel_tri'
  []
  [Tmod]
    block = 'moderator_quad moderator_tri'
  []
  [fuel_thermal_conductivity]
    block = 'fuel_quad fuel_tri'
    order = CONSTANT
    family = MONOMIAL
  []
  [fuel_specific_heat]
    block = 'fuel_quad fuel_tri'
    order = CONSTANT
    family = MONOMIAL
  []
  [monolith_thermal_conductivity]
    block = 'monolith '
    order = CONSTANT
    family = MONOMIAL
  []
  [monolith_specific_heat]
    block = 'monolith '
    order = CONSTANT
    family = MONOMIAL
  []
  [flux_uo] #auxvariable to hold heat pipe surface flux from UserObject
    block = 'reflector_quad monolith'
  []
  [flux_uo_corr] #auxvariable to hold corrected flux_uo
    block = 'reflector_quad monolith'
  []
  [hp_temp_aux]
    block = 'reflector_quad monolith'
  []
[]

[AuxKernels]
  [assign_tfuel]
    type = NormalizationAux
    variable = Tfuel
    source_variable = temp
    execute_on = 'timestep_end'
  []
  [assign_tmod]
    type = NormalizationAux
    variable = Tmod
    source_variable = temp
    execute_on = 'timestep_end'
  []
  [fuel_thermal_conductivity]
    type = MaterialRealAux
    variable = fuel_thermal_conductivity
    property = thermal_conductivity
    execute_on = timestep_end
  []
  [fuel_specific_heat]
    type = MaterialRealAux
    variable = fuel_specific_heat
    property = specific_heat
    execute_on = timestep_end
  []
  [monolith_thermal_conductivity]
    type = MaterialRealAux
    variable = monolith_thermal_conductivity
    property = thermal_conductivity
    execute_on = timestep_end
  []
  [monolith_specific_heat]
    type = MaterialRealAux
    variable = monolith_specific_heat
    property = specific_heat
    execute_on = timestep_end
  []
  [flux_uo]
    type = SpatialUserObjectAux
    variable = flux_uo
    user_object = flux_uo
  []
  [flux_uo_corr]
    type = NormalizationAux
    variable = flux_uo_corr
    source_variable = flux_uo
    normal_factor = '${fparse corr_factor}'
  []
[]

[BCs]
  [outside_bc]
    type = ConvectiveFluxFunction # (Robin BC)
    variable = temp
    boundary = 'bottom top side'
    coefficient = 1e3 # W/K/m^2
    T_infinity = 800 # K air temperature at the top of the core
  []
  [hp_temp]
    type = CoupledConvectiveHeatFluxBC
    boundary = 'heat_pipe_ht_surf'
    variable = temp
    T_infinity = hp_temp_aux
    htc = 750
  []
[]

[Materials]
  [fuel_matrix_thermal]
    type = GraphiteMatrixThermal
    block = 'fuel_quad fuel_tri'
    graphite_grade = IG_110
    packing_fraction = 0.4
    specific_heat_scale_factor = 1.0
    thermal_conductivity_scale_factor = 1.0
    fast_neutron_fluence = 0 #6.75E+24 # this value is nuetron fluence over 0.1MeV
    temperature = temp
  []
  [monolith_matrix_thermal]
    type = GraphiteMatrixThermal
    block = 'monolith '
    graphite_grade = IG_110
    packing_fraction = 0
    specific_heat_scale_factor = 1.0
    thermal_conductivity_scale_factor = 1.0
    fast_neutron_fluence = 0 #6.75E+24 # this value is nuetron fluence over 0.1MeV
    temperature = temp
  []
  [moderator_thermal]
    type = HeatConductionMaterial
    block = 'moderator_quad moderator_tri'
    temp = temp
    thermal_conductivity = 20 # W/m/K
    specific_heat = 500 # random value
  []
  [airgap_thermal]
    type = HeatConductionMaterial
    block = 'air_gap_tri air_gap_quad outer_shield' # Helium gap
    temp = temp
    thermal_conductivity = 0.15 # W/m/K
    specific_heat = 5197 # random value
  []
  [axial_reflector_thermal]
    type = HeatConductionMaterial
    block = 'reflector_tri reflector_quad'
    temp = temp
    thermal_conductivity = 199 # W/m/K
    specific_heat = 1867 # random value
  []
  [B4C_thermal]
    type = HeatConductionMaterial
    block = 'B4C'
    temp = temp
    thermal_conductivity = 92 # W/m/K
    specific_heat = 960 # random value
  []
  [SS_thermal]
    type = SS316Thermal
    temperature = temp
    block = mod_ss
  []
  [fuel_density]
    type = Density
    block = 'fuel_quad fuel_tri'
    density = 2276.5
  []
  [moderator_density]
    type = Density
    block = 'moderator_quad moderator_tri'
    density = 4.3e3
  []
  [monolith_density]
    type = Density
    block = 'monolith'
    density = 1806
  []
  [airgap_density]
    type = Density
    block = 'air_gap_tri air_gap_quad outer_shield' #helium
    density = 180
  []
  [axial_reflector_density]
    type = Density
    block = 'reflector_tri reflector_quad'
    density = 1848
  []
  [B4C_density]
    type = Density
    block = B4C
    density = 2510
  []
  [SS_density]
    type = Density
    density = 7990
    block = mod_ss
  []
[]

[MultiApps]
  [sockeye]
    type = TransientMultiApp
    app_type = SockeyeApp
    positions_file = 'hp_centers.txt'
    input_files = 'HPMR_sockeye_ss.i'
    execute_on = 'timestep_begin' # execute on timestep begin because hard to have a good initial guess on heat flux
    max_procs_per_app = 1
    output_in_position = true
    sub_cycling = true
  []
[]

[Transfers]
  [from_sockeye_temp]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = sockeye
    source_variable = hp_temp_aux
    variable = hp_temp_aux
    execute_on = 'timestep_begin'
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [to_sockeye_flux]
    type = MultiAppGeneralFieldUserObjectTransfer
    variable = master_flux
    to_multi_app = sockeye
    execute_on = 'timestep_begin'
    source_user_object = flux_uo
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
[]

[UserObjects]
  [flux_uo]
    type = NearestPointLayeredSideDiffusiveFluxAverage
    direction = z
    num_layers = 100
    points_file = 'hp_centers.txt'
    variable = temp
    diffusivity = thermal_conductivity
    execute_on = linear
    boundary = 'heat_pipe_ht_surf'
  []
[]

[Executioner]
  type = Transient

  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart '
  petsc_options_value = 'hypre boomeramg 100'

  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-8

  start_time = -5e4 # negative start time so we can start running from t = 0
  end_time = 0
  dtmin = 1
  # Multiple steps are needed to allow iteration between the thermo and heat pipe applications
  num_steps = 5
[]

[Postprocessors]
  [hp_heat_integral]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 'heat_pipe_ht_surf'
    diffusivity = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [hp_end_heat_integral]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 'heat_pipe_ht_surf_bot'
    diffusivity = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [ext_side_integral]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 'side'
    diffusivity = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [mirror_side_integral]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 'side_mirror'
    diffusivity = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [tb_integral]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 'top bottom'
    diffusivity = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [total_sink]
    type = SumPostprocessor
    values = 'hp_heat_integral hp_end_heat_integral ext_side_integral mirror_side_integral tb_integral'
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_avg]
    type = ElementAverageValue
    variable = temp
    block = 'fuel_quad fuel_tri'
  []
  [fuel_temp_max]
    type = ElementExtremeValue
    variable = temp
    block = 'fuel_quad fuel_tri'
  []
  [fuel_temp_min]
    type = ElementExtremeValue
    variable = temp
    block = 'fuel_quad fuel_tri'
    value_type = min
  []
  [mod_temp_avg]
    type = ElementAverageValue
    variable = temp
    block = 'moderator_quad moderator_tri'
  []
  [mod_temp_max]
    type = ElementExtremeValue
    variable = temp
    block = 'moderator_quad moderator_tri'
  []
  [mod_temp_min]
    type = ElementExtremeValue
    variable = temp
    block = 'moderator_quad moderator_tri'
    value_type = min
  []
  [monolith_temp_avg]
    type = ElementAverageValue
    variable = temp
    block = monolith
  []
  [monolith_temp_max]
    type = ElementExtremeValue
    variable = temp
    block = monolith
  []
  [monolith_temp_min]
    type = ElementExtremeValue
    variable = temp
    block = monolith
    value_type = min
  []
  [heatpipe_surface_temp_avg]
    type = SideAverageValue
    variable = temp
    boundary = heat_pipe_ht_surf
  []
  [power]
    type = ElementIntegralVariablePostprocessor
    block = 'fuel_quad fuel_tri'
    variable = power_density
    execute_on = 'initial timestep_end transfer'
  []
  [fuel_cP]
    type = ElementExtremeValue
    value_type = 'max'
    variable = fuel_specific_heat
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [fuel_k]
    type = ElementExtremeValue
    value_type = 'max'
    variable = fuel_thermal_conductivity
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [monolith_cP]
    type = ElementExtremeValue
    value_type = 'max'
    variable = monolith_specific_heat
    block = 'monolith '
    execute_on = 'initial timestep_end'
  []
  [monolith_k]
    type = ElementExtremeValue
    value_type = 'max'
    variable = monolith_thermal_conductivity
    block = 'monolith '
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  perf_graph = true
  exodus = true
  color = true
  csv = true
  checkpoint = true
[]

Sockeye Model

For each heat pipe in the 1/6 HP-MR core, a Sockeye grandchild application is used to calculate its thermal performance.

The effective thermal conductivity model, i.e., a 2D axisymmetric conduction model with a very high thermal conductivity of 2×10⁵ W/m•K is applied to the vapor core. A heat flux boundary condition is applied to the exterior of the casing in the evaporator section, which is provided by the bulk conduction model. A convective boundary condition is applied to the exterior of the envelope in the condenser section, with an external temperature of 800 K and a heat transfer coefficient (HTC) of 10⁶ W/m²•K (the value is intentionally set high to achieve an effective Dirichlet boundary condition).

################################################################################
## NEAMS Micro-Reactor Application Driver                                     ##
## Heat Pipe Microreactor Steady State                                        ##
## Sockeye Grandchild Application input file                                  ##
## Effective Heat Conduction Model with Operation Limits                      ##
################################################################################

# Total heat removed/added to heat pipe
Q_hp = 1800.

# Wick characteristics
R_pore = 15.0e-6
D_h_pore = '${fparse 2.0 * R_pore}'
permeability = 2e-9
porosity = 0.70

# Envelope ("env")
# SS316. Incropera & DeWitt, 3rd ed, Table A.1 @ 900K (627C)
# Density (kg/m3)
rho_env = 8238.
# Thermal conductivity (W/m-K)
k_env = 23.05
# Specific heat capacity (J/kg-K)
cp_env = 589.

# Potassium vapor
# From Appendix C of text book "Heat Pipe Design and Technology"
# Sat. Vapor Potassium at T = 650C (923 K)
# Density (kg/m3)
rho_vapor = 0.193
# Effective "super" conductivity (W/m-K)
#k_vapor = 1.0e5 # Reaches 600 seconds
k_vapor = 1.0e6
# Specific heat capacity (J/kg-K)
cp_vapor = 5320.

# Potassium liquid
# From Appendix C of text book "Heat Pipe Design and Technology"
# Sat. Liquid Potassium at T = 650C (923 K)
# Density (kg/m3)
rho_liquid = 695.4
# Thermal conductivity (W/m-K)
k_liquid = 40.08
# Specific heat capacity (J/kg-K)
# From Table 1.1, no temperature data given
cp_liquid = 810.
# Melting point, Table 3.1, lists 62C, rounding up
# T_melting = 340.

# Wick, homogenize envelope and fluid
# Density (kg/m3)
rho_wick = '${fparse porosity * rho_liquid + (1.0 - porosity) * rho_env}'
# Thermal conductivity (W/m-K)
k_wick = '${fparse porosity * k_liquid + (1.0 - porosity) * k_env}'
# Specific heat capacity (J/kg-K)
# From Table 1.1, no temperature data given
cp_wick = '${fparse porosity * cp_liquid + (1.0 - porosity) * cp_env}'

# Elevations and lengths
# Note: For blackbox model -- manually update "length" input
length_evap = 180.0e-2
length_adia = 30.0e-2
length_cond = 90.0e-2

# Mesh density
# The dimensions are nicely divisible by 3 cm mesh.
nelem_base_evap = 50
nelem_base_adia = 10
nelem_base_cond = 30
mesh_density = 3
nelem_evap = '${fparse mesh_density*nelem_base_evap}'
nelem_adia = '${fparse mesh_density*nelem_base_adia}'
nelem_cond = '${fparse mesh_density*nelem_base_cond}'

# Envelope thickness
t_env = 0.08e-2
# Liquid annulus thickness
t_ann = 0.07e-2
# Wick thickness
t_wick = 0.1e-2

# Radial geometry
# Envelope outer
R_hp_o = 1.05e-2
D_hp_o = '${fparse 2.0 * R_hp_o}'
# Inner Envelope/outer annulus
R_hp_i = '${fparse R_hp_o - t_env}'
D_hp_i = '${fparse 2.0 * R_hp_i}'
# Inner annulus/wick outer
R_wick_o = '${fparse R_hp_i - t_ann}'
D_wick_o = '${fparse 2.0 * R_wick_o}'
# Inner wick/vapor core outer
R_wick_i = '${fparse R_wick_o - t_wick}'
D_wick_i = '${fparse 2.0 * R_wick_i}'

# BCs for condenser
T_ext_cond = 800.
htc_ext_cond = 1.0e6

# Evaporator parameters
S_evap = '${fparse pi * D_hp_o * length_evap}'
q_evap = '${fparse Q_hp / S_evap}'

[GlobalParams]
  scaling_factor_temperature = 1e-2
  fp_2phase = fp_2phase
[]

[FluidProperties]
  [fp_2phase]
    type = PotassiumTwoPhaseFluidProperties
    emit_on_nan = none
  []
[]

[Components]
  [hp]
    type = HeatPipeConduction

    # Common to both HeatPipe2Phase and HeatPipeBlackbox
    position = '0 0 0'
    orientation = '0 0 1'
    length = '${length_evap} ${length_adia} ${length_cond}'
    n_elems = '${nelem_evap} ${nelem_adia} ${nelem_cond}'
    gravity_vector = '0 0 -9.8'
    D_wick_i = ${D_wick_i}
    D_wick_o = ${D_wick_o}
    R_pore = ${R_pore}
    porosity = ${porosity}
    permeability = ${permeability}

    # HeatPipeConduction only
    # Axial dimensions (for heat transfer & analytic limits)
    axial_region_names = 'evap adia cond'
    L_evap = ${length_evap}
    L_adia = ${length_adia}
    L_cond = ${length_cond}
    # Radial dimensions, mesh, and materials for heat transfer problem
    D_clad_o = ${D_hp_o}
    D_clad_i = ${D_hp_i}
    D_h_pore = ${D_h_pore}
    # Mesh
    n_elems_clad = 4
    n_elems_wick = 8
    n_elems_core = 10

    # Thermal properties
    sp_vapor = sp_vapor
    sp_liquid = sp_wick_ann
    sp_wick = sp_wick_ann
    sp_clad = sp_clad
    k_core = ${k_vapor}
    k_eff = ${k_wick}

    fp_2phase = fp_2phase
    evaporator_at_start_end = true
    # Initial temperature of block
    initial_T = ${T_ext_cond}
    T_ref = T_inner_avg

    # To evaluate the constant properties
    T_ref_density = 1000

    make_pressure_corrections = true
  []

  [condenser_boundary]
    type = HSBoundaryAmbientConvection
    boundary = 'hp:cond:outer'
    hs = hp
    T_ambient = ${T_ext_cond}
    htc_ambient = ${htc_ext_cond} #large value to approach an effective DirichletBC
    scale = 1
  []
  [evaporator_boundary]
    type = HSBoundaryExternalAppConvection
    boundary = 'hp:evap:outer'
    hs = hp
    T_ext = virtual_Text
    htc_ext = virtual_htc
  []
[]

[SolidProperties]
  [sp_vapor]
    type = ThermalFunctionSolidProperties
    rho = ${rho_vapor}
    cp = ${cp_vapor}
    k = ${k_vapor}
  []
  [sp_wick_ann]
    type = ThermalFunctionSolidProperties
    rho = ${rho_wick}
    cp = ${cp_wick}
    k = ${k_wick}
  []
  [sp_clad]
    type = ThermalFunctionSolidProperties
    rho = ${rho_env}
    cp = ${cp_env}
    k = ${k_env}
  []
[]

[UserObjects]
  [surf_T]
    type = LayeredSideAverage
    direction = z
    num_layers = 100
    variable = T_solid
    boundary = 'hp:evap:outer'
  []
[]

[AuxKernels]
  [hp_var]
    type = SpatialUserObjectAux
    variable = hp_temp_aux
    user_object = surf_T
  []
  [virtual_Text]
    type = ParsedAux
    variable = virtual_Text
    coupled_variables = 'T_solid master_flux virtual_htc'
    expression = 'master_flux/virtual_htc + T_solid'
  []
[]

[Functions]
  [scale_fcn]
    type = ParsedFunction
    symbol_names = 'catastrophic_pp recoverable_pp operational_pp'
    symbol_values = 'catastrophic_pp recoverable_pp operational_pp'
    expression = 'catastrophic_pp*recoverable_pp*operational_pp'
  []
[]

[AuxVariables]
  [T_wall_var]
    initial_condition = ${T_ext_cond}
  []
  [operational_aux]
    initial_condition = 1
  []
  [master_flux]
    initial_condition = ${q_evap}
  []
  [hp_temp_aux]
    initial_condition = ${T_ext_cond}
  []
  [virtual_Text]
    initial_condition = ${T_ext_cond}
  []
  [virtual_htc]
    initial_condition = 1.0
  []
[]

[Postprocessors]
  [Integral_BC_Total]
    type = SumPostprocessor
    values = 'condenser_boundary_integral evaporator_boundary_integral'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [ZeroPP]
    type = EmptyPostprocessor
  []
  [Integral_BC_Cond]
    type = DifferencePostprocessor
    value1 = ZeroPP
    value2 = condenser_boundary_integral
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Integral_BC_RelErr]
    type = RelativeDifferencePostprocessor
    value1 = evaporator_boundary_integral
    value2 = Integral_BC_Cond
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [operational_pp]
    type = ElementAverageValue
    variable = operational_aux
    execute_on = 'initial timestep_begin TIMESTEP_END'
  []
  # set `catastrophic_pp` as if it is recoverable for a solving to steady-state simulation
  # this MUST be changed back to catestrophic in transient simulations
  [catastrophic_pp]
    type = HeatRemovalRateLimitScale
    heat_addition_pps = 'evaporator_boundary_integral'
    limit_condenser_side = false
    recoverable_heat_removal_limit_pps = 'hp_boiling_limit hp_capillary_limit hp_entrainment_limit'
    catastrophic_heat_removal_limit_pps = ''
    T = T_inner_avg
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [recoverable_pp]
    type = HeatRemovalRateLimitScale
    heat_addition_pps = 'evaporator_boundary_integral'
    limit_condenser_side = false
    catastrophic_heat_removal_limit_pps = ''
    recoverable_heat_removal_limit_pps = 'hp_sonic_limit hp_viscous_limit'
    T = T_inner_avg
    execute_on = 'INITIAL linear nonlinear TIMESTEP_END'
  []
  [T_evap_inner]
    type = SideAverageValue
    boundary = hp:evap:inner
    variable = T_solid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_cond_inner]
    type = SideAverageValue
    boundary = hp:cond:inner
    variable = T_solid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_evap_outer]
    type = SideAverageValue
    boundary = hp:evap:outer
    variable = T_solid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_cond_outer]
    type = SideAverageValue
    boundary = hp:cond:outer
    variable = T_solid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_wall_var_avg]
    type = ElementAverageValue
    variable = T_wall_var
    execute_on = 'Initial timestep_end'
  []
  [T_inner_avg]
    type = SideAverageValue
    variable = T_solid
    boundary = hp:inner
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_inner_max]
    type = NodalExtremeValue
    variable = T_solid
    boundary = hp:inner
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_inner_min]
    type = NodalExtremeValue
    variable = T_solid
    boundary = hp:inner
    execute_on = 'INITIAL TIMESTEP_END'
    value_type = min
  []
  [DT_outer]
    type = DifferencePostprocessor
    value1 = T_evap_outer
    value2 = T_cond_outer
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [DT_inner]
    type = DifferencePostprocessor
    value1 = T_evap_inner
    value2 = T_cond_inner
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [A_int_master_flux]
    type = SideIntegralVariablePostprocessor
    variable = master_flux
    boundary = 'hp:evap:inner'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [A_int_T_ext]
    type = SideIntegralVariablePostprocessor
    variable = virtual_Text
    boundary = 'hp:evap:inner'
    execute_on = 'INITIAL LINEAR'
  []
  [A_avg_T_aux]
    type = AverageNodalVariableValue
    variable = hp_temp_aux
    boundary = 'hp:evap:inner'
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[VectorPostprocessors]
  [env_vpp]
    type = NodalValueSampler
    variable = T_solid
    block = 'hp:clad'
    sort_by = x
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

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

[Executioner]
  type = Transient

  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  scheme = bdf2
  line_search = none

  # ensure nl_abs_tol >> nl_rel_tol
  nl_abs_tol = 1e-6
  nl_rel_tol = 1e-8
  nl_max_its = 100

  l_tol = 1e-3
  l_max_its = 100

  start_time = -5e4 # negative start time so we can start running from t = 0
  end_time = 0
  dtmin = 1
  dt = 5000
[]

[Outputs]
  [console]
    type = Console
    max_rows = 5
    execute_postprocessors_on = 'INITIAL TIMESTEP_END FINAL FAILED'
  []
  [csv]
    type = CSV
    execute_on = 'INITIAL TIMESTEP_END FINAL FAILED'
    execute_vector_postprocessors_on = 'FINAL FAILED'
  []
[]

MultiApps Hierarchy

The three MOOSE applications that govern different physics respectively are coupled together leveraging the MOOSE MultiApps system. Each code has its own mesh and corresponding space and time scales. The MOOSE MultiApp system is used to tightly couple the different codes via fixed point iterations as illustrated below in Figure 3.

At the top level, Griffin calculates the power density based on its neutronics results given the fuel and moderator temperatures taken from the BISON child application. The calculated power density is also provided to BISON to update the temperature calculation. On the other hand, for the surface of each heat pipe, BISON calculates the layered average heat flux along the normal direction of the surface and transfers the 1D field values to the respective Sockeye sub-application so that Sockeye can use them as its boundary condition. Meanwhile, BISON takes the heat pipe surface temperature calculated by Sockeye back for its own boundary condition on the heat pipe surface.

Figure 3: MultiApps hierarchy of the HP-MR multiphysics simulations.

Simulation Cases

Three different scenarios are simulated in this HP-MR model. The most fundamental case is the steady state operation of the HP-MR, where the thermal power of the full reactor core is preset as 2.07 MW (i.e., 345.6 kW per 1/6 core). Using the steady state simulation results as the reference case, two different transient scenarios are also simulated: a load following transient and a transient event initiated by the failure of a single heat pipe.

Steady State Operation

For steady state simulation, the Griffin parent application performs an Eigenvalue calculation with a fixed preset integrated power (345.6 kW) to compute $var element = document.getElementById("moose-equation-0812b113-ab95-43e0-b202-6ef1b0fc6f77");katex.render("k_{eff}", element, {displayMode:false,throwOnError:false});$ as well as the shape of power distribution. In the meantime, the BISON child application works with multiple Sockeye grandchild sub-applications to perform transient thermal calculation using the given power density distribution until a thermal equilibrium is achieved. Once the fixed point iteration between different levels of MultiApps is converged, the steady state operation condition of the HP-MR at the given power level is obtained.

Optionally, the fixed preset integrated power level can be raised to force the HP-MR to operate at an overpower level. This is useful to demonstrate the capabilities of the NEAMS codes to predict some types of power transients that only occurs in the case of overpowered operation, such as the heat pipe cascade failure event to be demonstrated in this model.

Load Following Transient

The load following transient is initiated by introducing a drop in the heat removal capability of the secondary coolant loop that transfer the heat from the condensor regions of the heat pipes to either a heat exchanger or directly to the energy coversion component. During steady state operation, a high heat transfer coefficient (HTC) of 10⁶ W/m²•K is used on the outer surface of the heat pipe condensor region. To initiate the load following transient event, the HTC value is reduced by 99.99% to 10² W/m²•K to mimic a load reduction event. The same transient solving approach is used for BISON and Sockeye applications. However, instead of an Eigenvalue calculation, Griffin here also performs transient calculation. Based on the changes in $var element = document.getElementById("moose-equation-cfa092b8-61c4-4c0c-8f8f-06fb054d51ab");katex.render("k_{eff}", element, {displayMode:false,throwOnError:false});$ induced by temperature evolution calculated by BISON/Sockeye, Griffin predicts the corresponding kinetics of neutronics and consequent power distribution.

Transient Initiated by a Single Heat Pipe Failure

Another type of transient to be simulated in this model is more localized. Instead of reducing the condensor region HTC of all heat pipes, the condensor region HTC of a single heat pipe is nearly eliminated (i.e. reduced from 10⁶ W/m²•K to 10^-2 W/m²•K) to mimic a single heat pipe failure event, while the other heat pipes are kept unchanged (note that actually there are six heat pipes that are compromised due to the 1/6 symmetric mesh used here in this model). The failed heat pipe is intentionally selected to be in the hottest region of the HP-MR (i.e., in the center region of the center assembly).

Localized Single Heat Pipe Failure

The favored consequence of a single heat pipe failure in an HP-MR is that the extraneous heat that was removed by the failed heat pipe can be removed by the surrounding functioning heat pipes. In that case, the effects of the failed heat pipe are localized and the average core temperature only experiences a very minor increase. As a result, the HP-MR can still operate normally with a minor decrease in power. Ideally, the designed HP-MR power need to be set at an appropriate level that a single heat failure can be localized. Therefore, this is the expected simulation result for the HP-MR operating at the designed power level.

Heat Pipe Cascade Failure

On the other hand, when the HP-MR is operating at an overpower level (e.g., 720 kW in this model), the heat pipes in the hottest region are already operaing at a heat removal rate very close to their operating limits. Therefore, if a single heat pipe fails, its surrounding heat pipes will exceed their operating limits and then fail subsequently. This kind of failure events will then spread to a wider range and lead to heat pipe cascade failure. As a significant fraction of the heat pipes in an HP-MR fail, the average core temperature increases prominently and leads to a major power drop. Consequently, the HP-MR can no longer maintain its normal operation. A reliable HP-MR model should be able to predict such an unfavored transient event. In this model, the HP-MR is also intentionally set to run at an overpower of 720 kW before a single heat pipe failure is introduced to demonstrate such modeling capability.

Guidance of Running Different Simulation Cases

Fist of all, the user needs to run the mesh generation input file to generate the needed EXODUS mesh file that is required by the Griffin simulation. The mesh generation input file only uses objects available in the MOOSE Framework and the Reactor Module, so the input file can be run by any MOOSE applications that are compiled with the Reactor Module, such as Griffin. The example command to generate the mesh EXODUS file is listed as follows:

cd /mrad/mesh
griffin-opt -i HPMR_OneSixth_Core_meshgenerator_tri.i --mesh-only HPMR_OneSixth_Core_meshgenerator_tri_rotate_bdry.e

The mesh used by BISON is finer than the Griffin mesh, so some modifications need to be made to file. Also, the heat pipe surface boundaries need to be defined, and the heat pipe blocks themselves need to be removed. Thus, a separate input file is provided for the BISON mesh generation. Users can run the following command to generate the BISON mesh as an EXODUS file:

cd /mrad/mesh
griffin-opt -i HPMR_OneSixth_Core_meshgenerator_tri_fine.i --mesh-only HPMR_OneSixth_Core_meshgenerator_tri_rotate_bdry_fine.e

Then, the steady state simulation can be run. As the output files of the steady state simulation are already required as the initial status of the transient simulations, it is also required to be run before any transient simulations. The multiphysics steady state simulation uses Griffin, BISON and Sockeye objects. Therefore, users should use a super application that covers these three applications, such as DireWolf or BlueCRAB (compiled with Sockeye included). The example command to run the steady state simulation is shown as folllows:

cd /mrad/steady
mpirun -n 40 dire_wolf-opt -i HPMR_dfem_griffin_ss.i

With the EXODUS and checkpoint output file of the steady state simulation, the transient simulatons can be run. Before running any real cases, it is recommended to verify that the steady state obtained is correct. To do this, a null transient (i.e., a time dependent simulation with exactly the same conditions as the steady state simulation). If the steady state is correct, the null transient should yield constant values for all time-dependent quantities. The application needed to run transient simulations should be the same as the one used to obtain the steady state. The example command to run the null transient simulation is shown as folllows:

cd /mrad/transient_null
mpirun -n 40 dire_wolf-opt -i HPMR_dfem_griffin_trN.i

Once the steady state obtained is verified, the real transient cases can be run using the similar approach.

cd /mrad/load_following
mpirun -n 40 dire_wolf-opt -i HPMR_dfem_griffin_tr.i

cd /mrad/heat_pipe_failure
mpirun -n 40 dire_wolf-opt -i HPMR_dfem_griffin_tr.i

Overpower Simulation

The default power values in the input files in this repository are all set at the nominal power (i.e., 345.6 kW for the 1/6 HP-MR core). To run the overpower case to simulate the heat pipe cascade failure, the following values need to be changed to the 720 kW overpower level. The steady state configuration must be re-generated before the overpower transient may be run.

[PowerDensity]
  power = 345.6e3
  power_density_variable = power_density
  integrated_power_postprocessor = integrated_power
[]

[PowerDensity]
  power = 345.6e3
  power_density_variable = power_density
  integrated_power_postprocessor = integrated_power
[]

[PowerDensity]
  power = 345.6e3
  power_density_variable = power_density
  integrated_power_postprocessor = integrated_power
[]

References

Christopher Matthews, Vincent Laboure, Mark DeHart, Joshua Hansel, David Andrs, Yaqi Wang, Javier Ortensi, and Richard C Martineau. Coupled multiphysics simulations of heat pipe microreactors using direwolf. Nuclear Technology, 207(7):1142–1162, 2021.

@article{matthews2021coupled,
    author = "Matthews, Christopher and Laboure, Vincent and DeHart, Mark and Hansel, Joshua and Andrs, David and Wang, Yaqi and Ortensi, Javier and Martineau, Richard C",
    title = "Coupled multiphysics simulations of heat pipe microreactors using DireWolf",
    journal = "Nuclear Technology",
    volume = "207",
    number = "7",
    pages = "1142--1162",
    year = "2021",
    publisher = "Taylor \\& Francis"
}

E Shemon, Y Miao, S Kumar, K Mo, YS Jung, A Oaks, Guillaume Giudicelli, Logan Harbour, and Roy Stogner. Moose reactor module meshing enhancements to support reactor analysis. Technical Report ANL/NSE-22/65, Argonne National Lab.(ANL), Argonne, IL (United States), 2022.

@techreport{shemon2022moose,
    author = "Shemon, E and Miao, Y and Kumar, S and Mo, K and Jung, YS and Oaks, A and Giudicelli, Guillaume and Harbour, Logan and Stogner, Roy",
    title = "MOOSE Reactor Module Meshing Enhancements to Support Reactor Analysis",
    number = "ANL/NSE-22/65",
    year = "2022",
    institution = "Argonne National Lab.(ANL), Argonne, IL (United States)"
}

(microreactors/mrad/mesh/HPMR_OneSixth_Core_meshgenerator_tri.i)

################################################################
# This is the MOOSE input file to generate the mesh that can be
# used for the 1/6 core Heat-Pipe Microreactor Multiphysics
# simulations (Griffin neutronics simulation).
# Running this input requires MOOSE Reactor Module Objects
# Users should use
# --mesh-only HPMR_OneSixth_Core_meshgenerator_tri_rotate_bdry.e
# command line argument to generate
# exodus file for further Multiphysics simulations.
################################################################
[GlobalParams]
  create_outward_interface_boundaries = false
[]

[Mesh]
  [BatchMeshGeneratorAction]
    # We generate each control drum in 2 steps (concentric circles generation and azimuthal absorber/reflector splitting),
    # this syntax allows us to input the control drum (and reflector) generation in a very compact manner.
    [control_drums_ccg]
      mesh_generator_type = 'HexagonConcentricCircleAdaptiveBoundaryMeshGenerator'
      mesh_name_prefix = 'cd0'
      multi_batch_params_method = 'corresponding'
      # The main different parameters between the control drums are how the meshes are adapted to the reactor core
      # The batch parameters here are used to account for these differences
      batch_vector_input_param_names = 'meshes_to_adapt_to sides_to_adapt'
      batch_vector_input_param_types = 'MGNAME UINT'
      batch_vector_input_param_values = 'Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned|
                                         2 3 4;2 3;1 2 3;1 2;0 1 2;0 1;
                                         0 1 5;0 5;0 4 5;4 5;3 4 5;3 4'
      # The rest of the control drum parameters are the same for all control drums
      fixed_vector_input_param_names = 'num_sectors_per_side background_block_ids ring_radii ring_intervals ring_block_ids'
      fixed_vector_input_param_types = 'UINT USHORT REAL UINT USHORT'
      fixed_vector_input_param_values = '4 4 4 4 4 4;
                                         504;
                                         12.25 13.25;
                                         2 1;
                                         500 501 502'
      fixed_scalar_input_param_names = 'hexagon_size background_intervals preserve_volumes is_control_drum create_outward_interface_boundaries'
      fixed_scalar_input_param_types = 'REAL UINT BOOL BOOL BOOL'
      fixed_scalar_input_param_values = '13.376 2 true true false'
    []
    [control_drums_split]
      mesh_generator_type = 'AzimuthalBlockSplitGenerator'
      mesh_name_prefix = 'cd'
      multi_batch_params_method = 'corresponding'
      # Based on the locations of the control drums, we need to split the blocks in different ways
      batch_scalar_input_param_names = 'input start_angle'
      batch_scalar_input_param_types = 'MGNAME REAL'
      batch_scalar_input_param_values = 'cd0_0 cd0_1 cd0_2 cd0_3 cd0_4 cd0_5 cd0_6 cd0_7 cd0_8 cd0_9 cd0_10 cd0_11;
                                         15 345 315 285 255 225 195 165 135 105 75 45'
      fixed_vector_input_param_names = 'old_blocks new_block_ids'
      fixed_vector_input_param_types = 'SDNAME USHORT'
      fixed_vector_input_param_values = '502;503'
      fixed_scalar_input_param_names = 'angle_range'
      fixed_scalar_input_param_types = 'REAL'
      fixed_scalar_input_param_values = '90'
    []
    [reflector_blocks]
      mesh_generator_type = 'HexagonConcentricCircleAdaptiveBoundaryMeshGenerator'
      mesh_name_prefix = 'ref'
      multi_batch_params_method = 'corresponding'
      # Similar to the control drums, we need to adapt the reflector meshes to the reactor core based on the locations
      batch_vector_input_param_names = 'sides_to_adapt'
      batch_vector_input_param_types = 'UINT'
      batch_vector_input_param_values = '0;1;2;3;4;5'
      fixed_vector_input_param_names = 'meshes_to_adapt_to num_sectors_per_side background_block_ids'
      fixed_vector_input_param_types = 'MGNAME UINT USHORT'
      fixed_vector_input_param_values = 'Patterned;
                                         4 4 4 4 4 4;
                                         400 401'
      fixed_scalar_input_param_names = 'hexagon_size background_intervals create_outward_interface_boundaries'
      fixed_scalar_input_param_types = 'REAL UINT BOOL'
      fixed_scalar_input_param_values = '13.376 2 false'
    []
  []
  # Moderator pins
  [Mod_hex]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2 '
    background_intervals = 1
    background_block_ids = '10'
    polygon_size = 1.15
    polygon_size_style ='apothem'
    ring_radii = '0.825 0.92'
    ring_intervals = '2 1'
    ring_block_ids = '103 100 101' # 103 is tri mesh
    preserve_volumes = on
    quad_center_elements = false
  []
  # Heat Pipe pins
  [HP_hex]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2 '
    background_intervals = 1
    background_block_ids = '10'
    polygon_size = 1.15
    polygon_size_style ='apothem'
    ring_radii = '0.97 1.07'
    ring_intervals = '2 1'
    ring_block_ids = '203 200 201' # 203 is tri mesh
    preserve_volumes = on
    quad_center_elements = false
  []
  # Fuel pins
  [Fuel_hex]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2 '
    background_intervals = 1
    background_block_ids = '10'
    polygon_size = 1.15
    polygon_size_style ='apothem'
    ring_radii = '1'
    ring_intervals = '2'
    ring_block_ids = '303 301'  # 303 is tri mesh
    preserve_volumes = on
    quad_center_elements = false
  []
  # Assembly mesh with all the three types of pins
  [Patterned]
    type = PatternedHexMeshGenerator
    inputs = 'Fuel_hex HP_hex Mod_hex'
    hexagon_size = 13.376
    background_block_id = 10
    background_intervals = 1
    pattern = '1 0 1 0 1 0 1;
              0 2 0 2 0 2 0 0;
             1 0 1 0 1 0 1 2 1;
            0 2 0 2 0 2 0 0 0 0;
           1 0 1 0 1 0 1 2 1 2 1;
          0 2 0 2 0 2 0 0 0 0 0 0;
         1 0 1 0 1 0 1 2 1 2 1 2 1;
          0 2 0 2 0 2 0 0 0 0 0 0;
           1 0 1 0 1 0 1 2 1 2 1;
            0 2 0 2 0 2 0 0 0 0;
             1 0 1 0 1 0 1 2 1;
              0 2 0 2 0 2 0 0;
               1 0 1 0 1 0 1'
  []
  # Central void region
  [air_center]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    num_sectors_per_side= '4 4 4 4 4 4'
    meshes_to_adapt_to = 'Patterned Patterned Patterned Patterned Patterned Patterned'
    sides_to_adapt = '0 1 2 3 4 5'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = '600 601'
  []
  # Dummy assembly to help form a regular hexagonal core pattern
  # This will be deleted later
  [dummy]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    num_sectors_per_side= '4 4 4 4 4 4'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = '700 701'
       # external_boundary_id = 9998
  []
  # Stitching assemblies together to form the core mesh
  [core]
    type = PatternedHexMeshGenerator
    inputs = 'Patterned cd_0 cd_1 cd_2 cd_3 cd_4 cd_5 cd_6 cd_7 cd_8 cd_9 cd_10 cd_11 ref_0 ref_1 ref_2 ref_3 ref_4 ref_5 dummy air_center'
    # Pattern ID  0       1    2    3    4    5    6    7    8    9    10    11    12    13    14    15    16    17    18   19       20
    pattern_boundary = none
    generate_core_metadata = true
    pattern = '19 17 12 16 19;
             17 11  0  0  1 16;
           10  0  0  0  0  0  2;
          18 0  0   0  0  0  0 15;
        19  9  0  0  20  0  0  3 19;
          18 0  0   0  0  0  0 15;
            8  0  0  0  0  0  4;
             13  7  0  0  5 14;
               19 13  6 14 19'
    rotate_angle = 60
    assign_control_drum_id = true
  []
  # Delete the dummy assemblies
  [del_dummy]
    type = BlockDeletionGenerator
    block = '700 701'
    input = core
    new_boundary = 10000
  []
  # Add peripheral circular region
  [add_outer_shield]
    type = PeripheralRingMeshGenerator
    input = del_dummy
    peripheral_layer_num = 1
    peripheral_ring_radius = 115.0
    input_mesh_external_boundary = 10000
    peripheral_ring_block_id = 250
    peripheral_ring_block_name = outer_shield
 []
 # trim the full core to 1/6 core
 [del_1]
    type = PlaneDeletionGenerator
    point = '0 0 0'
    normal = '10 17.32 0'
    input = add_outer_shield
    new_boundary = 147
  []
  [del_2]
    type = PlaneDeletionGenerator
    point = '0 0 0'
    normal = '10 -17.32 0'
    input = del_1
    new_boundary = 147
  []
  # extrusion
  [extrude]
     type = AdvancedExtruderGenerator
     input = del_2
     heights = '20 160 20'
     # Use `num_layers = '6 16 6'` for BISON mesh
     num_layers = '1 8 1'
     # biased upper and lower reflector mesh, uncomment for BISON mesh
     # biases = '1.6 1.0 0.625'
     direction = '0 0 1'
     top_boundary = 2000
     bottom_boundary = 3000
  []
  # Define some special reflector blocks
  [reflector_bottom_quad]
    type = ParsedSubdomainMeshGenerator
    input = extrude
    combinatorial_geometry = 'z<=20'
    block_id = 1000
    excluded_subdomains = '103 203 303 250 400 401 500 501 502 503 504 600 601'
  []
  [reflector_bottom_tri]
    type = ParsedSubdomainMeshGenerator
    input = reflector_bottom_quad
    combinatorial_geometry = 'z<=20'
    block_id = 1003
    excluded_subdomains = '1000 250 400 401 500 501 502 503 504 600 601'
  []
  [reflector_top_quad]
    type = ParsedSubdomainMeshGenerator
    input = reflector_bottom_tri
    combinatorial_geometry = 'z>=180'
    block_id = 1000
    excluded_subdomains = '103 203 303 200 201 250 400 401 500 501 502 503 504 600 601'
  []
  [reflector_top_tri]
    type = ParsedSubdomainMeshGenerator
    input = reflector_top_quad
    combinatorial_geometry = 'z>=180'
    block_id = 1003
    excluded_subdomains = '1000 200 201 203 250 400 401 500 501 502 503 504 600 601'
  []
  # Assgin block names
  [rename_blocks]
    type = RenameBlockGenerator
    old_block = '     101       100         103          201        200           203            301       303        400            401         10      503      600         601        504          500                501          502         1000           1003'
    new_block = ' mod_ss moderator_quad moderator_tri hp_ss  heat_pipes_quad heat_pipes_tri fuel_quad fuel_tri reflector_tri reflector_quad monolith B4C air_gap_tri air_gap_quad air_gap_quad reflector_tri  reflector_quad reflector_quad reflector_quad reflector_tri'
    input = reflector_top_tri
  []
  # Assign boundary names
  [rename_boundaries]
    type = RenameBoundaryGenerator
    input = rename_blocks
    old_boundary = ' 10000 2000 3000 147'
    new_boundary = ' side top bottom side_mirror'
  []
  # Scale to m unit for general MOOSE application
  # Note that is_meter needs to be true in Griffin
  [scale]
    type = TransformGenerator
    input = rename_boundaries
    transform = SCALE
    vector_value = '1e-2 1e-2 1e-2'
  []
  # Rotate to the desired orientation for simulation
  [rotate_end]
    type = TransformGenerator
    input = scale
    transform = ROTATE
    vector_value = '0 0 -120'
  []
[]

(microreactors/mrad/mesh/HPMR_OneSixth_Core_meshgenerator_tri_fine.i)

#####################################################################
# This is the MOOSE input file to generate the mesh that can be
# used for the 1/6 core Heat-Pipe Microreactor Multiphysics
# simulations (BISON thermal simulation).
# Running this input requires MOOSE Reactor Module Objects
# Users should use
# --mesh-only HPMR_OneSixth_Core_meshgenerator_tri_rotate_bdry_fine.e
# command line argument to generate
# exodus file for further Multiphysics simulations.
#####################################################################
[GlobalParams]
  # we use ':=' to override the values already defined in the !include file
  create_outward_interface_boundaries := false
[]

[Mesh]
  [HP_hex]
    # inner background boundary layer, only for BISON mesh
    background_inner_boundary_layer_bias = 1.5
    background_inner_boundary_layer_intervals = 3
    background_inner_boundary_layer_width = 0.03
  []
  [add_outer_shield]
    peripheral_layer_num := 2
    # Peripheral boundary layer, only for BISON mesh
    peripheral_outer_boundary_layer_bias = 0.625
    peripheral_outer_boundary_layer_intervals = 3
    peripheral_outer_boundary_layer_width = 2
  []
  [extrude]
    num_layers := '6 16 6'
    # biased upper and lower reflector mesh, only for BISON mesh
    biases = '1.6 1.0 0.625'
  []
  [split_hp_ss]
    type = SubdomainBoundingBoxGenerator
    input = rotate_end
    block_id = 1201
    block_name = 'hp_ss_up'
    restricted_subdomains = 'hp_ss'
    bottom_left = '-100 -100 1.8'
    top_right = '100 100 3.0'
  []
  [add_exterior_ht_low]
    type = SideSetsBetweenSubdomainsGenerator
    input = split_hp_ss
    paired_block = 'hp_ss'
    primary_block = 'monolith'
    new_boundary = 'heat_pipe_ht_surf_low'
  []
  [add_exterior_ht_up]
    type = SideSetsBetweenSubdomainsGenerator
    input = add_exterior_ht_low
    paired_block = 'hp_ss_up'
    primary_block = 'reflector_quad'
    new_boundary = 'heat_pipe_ht_surf_up'
  []
  [add_exterior_bot]
    type = SideSetsBetweenSubdomainsGenerator
    input = add_exterior_ht_up
    paired_block = 'heat_pipes_quad heat_pipes_tri hp_ss'
    primary_block = 'reflector_quad reflector_tri'
    new_boundary = 'heat_pipe_ht_surf_bot'
  []
  [merge_hp_surf]
    type = RenameBoundaryGenerator
    input = add_exterior_bot
    old_boundary = 'heat_pipe_ht_surf_low heat_pipe_ht_surf_up'
    new_boundary = 'heat_pipe_ht_surf heat_pipe_ht_surf'
  []
  [remove_hp]
    type = BlockDeletionGenerator
    input = merge_hp_surf
    block = 'hp_ss hp_ss_up heat_pipes_quad heat_pipes_tri'
  []

  # remove extra nodesets to limit the size of the mesh
  [clean_up]
    type = BoundaryDeletionGenerator
    input = remove_hp
    boundary_names = '1 3'
  []
[]

# This input defines the mesh generation for the Griffin mesh
!include HPMR_OneSixth_Core_meshgenerator_tri.i

(microreactors/mrad/Serpent_Model/serpent_input.i)

%----------------------------------------------------------------------------------------------------
% ANL HPMR Serpent (Version 2.1.32) Model Input
% If using or referring to this model, please cite as explained in
% https://mooseframework.inl.gov/virtual_test_bed/citing.html
%Contacts: Ahmed Abdelhameed (aabdelhameed.at.anl.gov), Nicolas Stauff (nstauff.at.anl.gov)
% -----------------------------------------------------------------------------
set title "ANL MicroReactor Core - pf - PF = 40"

pbed 801 33 "TRISO_U900_PF40_R100.inp"

% --- Compact --------------------------------------------------------------------------------------
surf  1 hexyprism 0.0 0.0 13.376      0 200 %fueled core
surf 91 hexyprism 0.0 0.0 130.0       0 200 %fueled core

surf 95l pz  20
surf 95u pz  180
surf 96l pz  18.0
surf 96u pz  182.0

surf  2 cyl 0.0 0.0 1.00   0.0 200
surf  3 cyl 0.0 0.0 0.825  0.0 200
surf 35 cyl 0.0 0.0 0.875  0.0 200
surf 36 cyl 0.0 0.0 0.900  0.0 200
surf 3g cyl 0.0 0.0 0.920  0.0 200
surf 5g cyl 0.0 0.0 1.07   0.0 200
surf  5 cyl 0.0 0.0 1.05   0.0 200
surf 51 cyl 0.0 0.0 0.97   0.0 200
surf 52 cyl 0.0 0.0 0.90   0.0 200
surf 53 cyl 0.0 0.0 0.80   0.0 200

surf inf inf

particle 900
fuel      2.1250e-02
buffer    3.1250e-02
PyC1      3.5250e-02
SiC       3.8750e-02
PyC2      4.2750e-02
matrix_pin


% infinite cells defining material universes
cell 51 802 moderator  -inf
cell 52 33  matrix_pin -inf
cell 53 803    matrix  -inf
cell 55 804       air  -inf
cell 56 805  beryllium -inf
cell 57 806  shell_mod -inf
cell 64 813  shell_hp  -inf
cell 58 820  shell_air_center -inf

% mixed wick and coo_vap and liqu
cell 65 815 hp_vp_liq_wick -inf

% mixed air and ss in moderator
cell 68 816 shell_air_mod -inf

% mixed air and ss in heatpipe
cell 72 817 shell_air_hp   -inf

%  control drums
cell 61 810 beryllium_drum -inf
cell 62 811   B4C_drum -inf
cell 63 812   B4C_central -inf

% fuel
cell 11   3 fill 801        -2  95l -95u % fuel compact
cell 12   3 fill 803         2  95l -95u
cell 125  3 fill 805                 95u
cell 126  3 fill 805                -95l

% Yan simplified model removed caps for the HP and moderator pins extended into the Be reflectors
% moderator
cell 13   2 fill 802    -3      95l -95u % active length
cell 14   2 fill 816     3 -3g  95l -95u
cell 15   2 fill 803     3g     95l -95u
cell 151  2 fill 805       -36  96l -95l % 0.05 cm shell
cell 152  2 fill 805       -36  95u -96u % 0.05 cm shell
cell 151a 2 fill 805   -3g  36  96l -95l % 0.05 cm shell
cell 152a 2 fill 805   -3g  36  95u -96u % 0.05 cm shell
cell 153  2 fill 805        3g  96l -95l % 0.05 cm shell
cell 154  2 fill 805        3g  95u -96u % 0.05 cm shell
cell 155  2 fill 805                 96u
cell 156  2 fill 805                -96l

% heat pipe cell
cell 16i  1 fill 815      -51 95l
cell 16o  1 fill 817   51 -5g 95l
cell 17   1 fill 803       5g 95l -95u
cell 175  1 fill 805       5g 95u
cell 176  1 fill 805         -96l
cell 177  1 fill 805      -5  96l -95l
cell 178  1 fill 805       5g 96l -95l
cell 178a 1 fill 805   -5g 5  96l -95l

% monolith filled cell
% cell 18  8 fill 804     -inf
cell 18  8 fill 820     -inf

% monolith filled cell
cell 19  9 fill 803     -inf

% monolith filled cell
cell 20 10 fill 805           -inf

% monolith filled cell
% cell 21 11 fill 804           -inf
cell 21 11 fill 820           -inf

% --- Assembly lattice

lat 20  2  0.0 0.0 19 19 2.3
% New Assembly Configuration - 2 time less moderator - more fuel
 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
  9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
   9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
    9 9 9 9 9 9 9 9 9 1 3 1 3 1 3 1 9 9 9
     9 9 9 9 9 9 9 9 3 3 2 3 2 3 2 3 9 9 9
      9 9 9 9 9 9 9 1 2 1 3 1 3 1 3 1 9 9 9
       9 9 9 9 9 9 3 3 3 3 2 3 2 3 2 3 9 9 9
        9 9 9 9 9 1 2 1 2 1 3 1 3 1 3 1 9 9 9
         9 9 9 9 3 3 3 3 3 3 2 3 2 3 2 3 9 9 9
          9 9 9 1 2 1 2 1 2 1 3 1 3 1 3 1 9 9 9
           9 9 9 3 3 3 3 3 3 2 3 2 3 2 3 9 9 9 9
            9 9 9 1 2 1 2 1 3 1 3 1 3 1 9 9 9 9 9
             9 9 9 3 3 3 3 2 3 2 3 2 3 9 9 9 9 9 9
              9 9 9 1 2 1 3 1 3 1 3 1 9 9 9 9 9 9 9
               9 9 9 3 3 2 3 2 3 2 3 9 9 9 9 9 9 9 9
                9 9 9 1 3 1 3 1 3 1 9 9 9 9 9 9 9 9 9
                 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
                  9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
                   9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9

trans 20 0.0 0.0 0.0 0.0 0.0 30.

lat 80  2  0.0 0.0 1 1 26.752
 10  % Be

lat 90  2  0.0 0.0 1 1 26.752
  8  % air

lat 91  2  0.0 0.0 1 1 26.752
  11  % air o  B4C

lat 50  2  0.0 0.0 13 13 26.752
 90 90 90 90 90 90 90 90 90 90 90 90 90
  90 90 90 90 90 90 90 90 90 90 90 90 90
   90 90 90 90 90 90 90 80 90 80 90 90 90
    90 90 90 90 90 80 90 20 20 90 80 90 90
     90 90 90 90 90 20 20 20 20 20 90 90 90
      90 90 90 80 20 20 20 20 20 20 80 90 90
       90 90 90 90 20 20 91 20 20 90 90 90 90
        90 90 80 20 20 20 20 20 20 80 90 90 90
         90 90 90 20 20 20 20 20 90 90 90 90 90
          90 90 80 90 20 20 90 80 90 90 90 90 90
           90 90 90 80 90 80 90 90 90 90 90 90 90
            90 90 90 90 90 90 90 90 90 90 90 90 90
             90 90 90 90 90 90 90 90 90 90 90 90 90

% -----------DRUMS DEFINITTION---------------------------------------------------------------------------------------


cell 102 0 fill 50 -91 731 732 733 734 735 736 737 738 739 740 741 742
cell 180 0 fill       731 -731
cell 181 0 fill       732 -732
cell 182 0 fill       733 -733
cell 183 0 fill       734 -734
cell 184 0 fill       735 -735
cell 185 0 fill       736 -736
cell 186 0 fill       737 -737
cell 187 0 fill       738 -738
cell 188 0 fill       739 -739
cell 189 0 fill       740 -740
cell 190 0 fill       741 -741
cell 191 0 fill       742 -742


surf 731 cyl 80.2560   0.00 13.2500   0.00 200.00
surf 831 cyl 80.2560   0.00 12.2500   0.00 200.00
surf 771 plane 1.0000   0.00 0.0  88.92
surf 871 plane 1.0000  -0.00 0.0  71.59
cell 801  731 fill 810   -731   -771
cell 802  731 fill 811    831  -731  771
cell 803  731 fill 810   -831   771
surf 732 cyl 40.1280  69.50 13.2500   0.00 200.00
surf 832 cyl 40.1280  69.50 12.2500   0.00 200.00
surf 772 plane 1.0000   1.73 0.0 177.84
surf 872 plane 1.0000   1.73 0.0 143.19
cell 806  732 fill 810   -732   -772
cell 807  732 fill 811    832  -732  772
cell 808  732 fill 810   -832   772
surf 733 cyl -40.1280  69.50 13.2500   0.00 200.00
surf 833 cyl -40.1280  69.50 12.2500   0.00 200.00
surf 773 plane 1.0000  -1.73 0.0 -177.84
surf 873 plane 1.0000  -1.73 0.0 -143.19
cell 811  733 fill 810   -733   773
cell 812  733 fill 811    833  -733  -773
cell 813  733 fill 810   -833   -773
surf 734 cyl -80.2560   0.00 13.2500   0.00 200.00
surf 834 cyl -80.2560   0.00 12.2500   0.00 200.00
surf 774 plane 1.0000  -0.00 0.0 -88.92
surf 874 plane 1.0000   0.00 0.0 -71.59
cell 816  734 fill 810   -734   774
cell 817  734 fill 811    834  -734  -774
cell 818  734 fill 810   -834   -774
surf 735 cyl -40.1280 -69.50 13.2500   0.00 200.00
surf 835 cyl -40.1280 -69.50 12.2500   0.00 200.00
surf 775 plane 1.0000   1.73 0.0 -177.84
surf 875 plane 1.0000   1.73 0.0 -143.19
cell 821  735 fill 810   -735   775
cell 822  735 fill 811    835  -735  -775
cell 823  735 fill 810   -835   -775
surf 736 cyl 40.1280 -69.50 13.2500   0.00 200.00
surf 836 cyl 40.1280 -69.50 12.2500   0.00 200.00
surf 776 plane 1.0000  -1.73 0.0 177.84
surf 876 plane 1.0000  -1.73 0.0 143.19
cell 826  736 fill 810   -736   -776
cell 827  736 fill 811    836  -736  776
cell 828  736 fill 810   -836   776
surf 737 cyl 80.2563  46.34 13.2500   0.00 200.00
surf 837 cyl 80.2563  46.34 12.2500   0.00 200.00
surf 777 plane 1.0000   0.58 0.0 117.01
surf 877 plane 1.0000   0.58 0.0  97.01
cell 831  737 fill 810   -737   -777
cell 832  737 fill 811    837  -737  777
cell 833  737 fill 810   -837   777
surf 738 cyl 0.0000  92.67 13.2500   0.00 200.00
surf 838 cyl 0.0000  92.67 12.2500   0.00 200.00
surf 778 plane 0.0000   1.00 0.0 101.33
surf 878 plane 0.0000   1.00 0.0  84.01
cell 836  738 fill 810   -738   -778
cell 837  738 fill 811    838  -738  778
cell 838  738 fill 810   -838   778
surf 739 cyl -80.2563  46.34 13.2500   0.00 200.00
surf 839 cyl -80.2563  46.34 12.2500   0.00 200.00
surf 779 plane 1.0000  -0.58 0.0 -117.01
surf 879 plane 1.0000  -0.58 0.0 -97.01
cell 841  739 fill 810   -739   779
cell 842  739 fill 811    839  -739  -779
cell 843  739 fill 810   -839   -779
surf 740 cyl -80.2563 -46.34 13.2500   0.00 200.00
surf 840 cyl -80.2563 -46.34 12.2500   0.00 200.00
surf 780 plane 1.0000   0.58 0.0 -117.01
surf 880 plane 1.0000   0.58 0.0 -97.01
cell 846  740 fill 810   -740   780
cell 847  740 fill 811    840  -740  -780
cell 848  740 fill 810   -840   -780
surf 741 cyl -0.0000 -92.67 13.2500   0.00 200.00
surf 841 cyl -0.0000 -92.67 12.2500   0.00 200.00
surf 781 plane 0.0000   1.00 0.0 -101.33
surf 881 plane 1.0000 1219076296694385.50 0.0 -102414528894101376.00
cell 851  741 fill 810   -741   781
cell 852  741 fill 811    841  -741  -781
cell 853  741 fill 810   -841   -781
surf 742 cyl 80.2563 -46.34 13.2500   0.00 200.00
surf 842 cyl 80.2563 -46.34 12.2500   0.00 200.00
surf 782 plane 1.0000  -0.58 0.0 117.01
surf 882 plane 1.0000  -0.58 0.0  97.01
cell 856  742 fill 810   -742   -782
cell 857  742 fill 811    842  -742  782
cell 858  742 fill 810   -842   782

% --------------------------------------------------------------------------------------------------
cell 104 0 outside  91

set ures 1


% --------------------------------------------------------------------------------------------------
% Material data:

% materials

mat air -0.18e-3   rgb 255 255 255   tmp  700
 2004.01c  1

mat shell_mod -7.90   rgb 133 133 133    tmp  700
 6012.01c  1.9010E-03
 14028.01c  9.2693E-03
 14029.01c  4.7251E-04
 14030.01c  3.1166E-04
 15031.01c  4.1322E-04
 16032.01c  2.4710E-04
 16033.01c  1.9511E-06
 16034.01c  1.1056E-05
 16036.01c  3.0016E-08
 24050.01c  7.9116E-03
 24052.01c  1.5257E-01
 24053.01c  1.7300E-02
 24054.01c  4.3063E-03
 25055.01c  1.0280E-02
 26054.01c  3.9029E-02
 26056.01c  6.1213E-01
 26057.01c  1.4144E-02
 26058.01c  1.3343E-03
 28058.01c  7.7516E-02
 28060.01c  2.9859E-02
 28061.01c  1.2981E-03
 28062.01c  4.1389E-03
 28064.01c  1.0544E-03
 42092.01c  2.1259E-03
 42094.01c  1.3299E-03
 42095.01c  2.3030E-03
 42096.01c  2.4191E-03
 42097.01c  1.3903E-03
 42098.01c  3.5249E-03
42100.01c  1.4135E-03

mat shell_hp -7.90   rgb 133 133 133    tmp  700
 6012.01c  1.9010E-03
14028.01c  9.2693E-03
14029.01c  4.7251E-04
14030.01c  3.1166E-04
15031.01c  4.1322E-04
16032.01c  2.4710E-04
16033.01c  1.9511E-06
16034.01c  1.1056E-05
16036.01c  3.0016E-08
24050.01c  7.9116E-03
24052.01c  1.5257E-01
24053.01c  1.7300E-02
24054.01c  4.3063E-03
25055.01c  1.0280E-02
26054.01c  3.9029E-02
26056.01c  6.1213E-01
26057.01c  1.4144E-02
26058.01c  1.3343E-03
28058.01c  7.7516E-02
28060.01c  2.9859E-02
28061.01c  1.2981E-03
28062.01c  4.1389E-03
28064.01c  1.0544E-03
42092.01c  2.1259E-03
42094.01c  1.3299E-03
42095.01c  2.3030E-03
42096.01c  2.4191E-03
42097.01c  1.3903E-03
42098.01c  3.5249E-03
42100.01c  1.4135E-03

mat shell_air_mod 2.290E-02 rgb 133 133 133    tmp  700
   2004.01c   1.983E-05
   6012.01c   4.349E-05
   14028.01c   2.121E-04
   14029.01c   1.081E-05
   14030.01c   7.130E-06
   15031.01c   9.454E-06
   16032.01c   5.653E-06
   16033.01c   4.464E-08
   16034.01c   2.530E-07
   16036.01c   6.867E-10
   24050.01c   1.810E-04
   24052.01c   3.491E-03
   24053.01c   3.958E-04
   24054.01c   9.852E-05
   25055.01c   2.352E-04
   26054.01c   8.929E-04
   26056.01c   1.400E-02
   26057.01c   3.236E-04
   26058.01c   3.053E-05
   28058.01c   1.773E-03
   28060.01c   6.831E-04
   28061.01c   2.970E-05
   28062.01c   9.469E-05
   28064.01c   2.412E-05
   42092.01c   4.864E-05
   42094.01c   3.043E-05
   42095.01c   5.269E-05
   42096.01c   5.535E-05
   42097.01c   3.181E-05
   42098.01c   8.065E-05
   42100.01c   3.234E-05

mat shell_air_hp 6.771E-02 rgb 133 133 133    tmp  700
2004.01c  5.629E-06
6012.01c  1.287E-04
14028.01c  6.276E-04
14029.01c  3.199E-05
14030.01c  2.110E-05
15031.01c  2.798E-05
16032.01c  1.673E-05
16033.01c  1.321E-07
16034.01c  7.486E-07
16036.01c  2.032E-09
24050.01c  5.357E-04
24052.01c  1.033E-02
24053.01c  1.171E-03
24054.01c  2.916E-04
25055.01c  6.960E-04
26054.01c  2.643E-03
26056.01c  4.145E-02
26057.01c  9.576E-04
26058.01c  9.034E-05
28058.01c  5.248E-03
28060.01c  2.022E-03
28061.01c  8.789E-05
28062.01c  2.802E-04
28064.01c  7.139E-05
42092.01c  1.439E-04
42094.01c  9.004E-05
42095.01c  1.559E-04
42096.01c  1.638E-04
42097.01c  9.413E-05
42098.01c  2.387E-04
42100.01c  9.570E-05


mat shell_air_center 8.815E-04 rgb 133 133 133  tmp  700
2004.01c   2.68E-05
6012.01c   1.62E-06
14028.01c   7.92E-06
14029.01c   4.04E-07
14030.01c   2.66E-07
15031.01c   3.53E-07
16032.01c   2.11E-07
16033.01c   1.67E-09
16034.01c   9.45E-09
16036.01c   2.57E-11
24050.01c   6.76E-06
24052.01c   1.30E-04
24053.01c   1.48E-05
24054.01c   3.68E-06
25055.01c   8.79E-06
26054.01c   3.34E-05
26056.01c   5.23E-04
26057.01c   1.21E-05
26058.01c   1.14E-06
28058.01c   6.63E-05
28060.01c   2.55E-05
28061.01c   1.11E-06
28062.01c   3.54E-06
28064.01c   9.01E-07
42092.01c   1.82E-06
42094.01c   1.14E-06
42095.01c   1.97E-06
42096.01c   2.07E-06
42097.01c   1.19E-06
42098.01c   3.01E-06
42100.01c   1.21E-06

mat coo_vap  -1.11e-4 rgb 250 50 250    tmp  700
19039.01c   0.93258
19040.01c   0.00012
19041.01c   0.06730

mat coo_liq  -0.705 rgb 50 250 50    tmp  700
19039.01c   0.93258
19040.01c   0.00012
19041.01c   0.06730

mat wick  -2.753 rgb 250 250 50     tmp  700
 6012.01c  5.589E-04
14028.01c  2.725E-03
14029.01c  1.389E-04
14030.01c  9.163E-05
15031.01c  1.215E-04
16032.01c  7.265E-05
16033.01c  5.736E-07
16034.01c  3.250E-06
16036.01c  8.825E-09
24050.01c  2.326E-03
24052.01c  4.485E-02
24053.01c  5.086E-03
24054.01c  1.266E-03
25055.01c  3.022E-03
26054.01c  1.147E-02
26056.01c  1.800E-01
26057.01c  4.158E-03
26058.01c  3.923E-04
28058.01c  2.279E-02
28060.01c  8.778E-03
28061.01c  3.816E-04
28062.01c  1.217E-03
28064.01c  3.100E-04
42092.01c  6.250E-04
42094.01c  3.910E-04
42095.01c  6.771E-04
42096.01c  7.112E-04
42097.01c  4.087E-04
42098.01c  1.036E-03
42100.01c  4.156E-04
19039.01c  6.584E-01
19040.01c  8.472E-05
19041.01c  4.751E-02

mat hp_vp_liq_wick  8.324E-03  rgb 250 250 50     tmp  700
   6012.01c   3.808E-06
   14028.01c   1.856E-05
   14029.01c   9.463E-07
   14030.01c   6.242E-07
   15031.01c   8.277E-07
   16032.01c   4.949E-07
   16033.01c   3.908E-09
   16034.01c   2.214E-08
   16036.01c   6.012E-11
   24050.01c   1.585E-05
   24052.01c   3.055E-04
   24053.01c   3.465E-05
   24054.01c   8.625E-06
   25055.01c   2.059E-05
   26054.01c   7.814E-05
   26056.01c   1.226E-03
   26057.01c   2.833E-05
   26058.01c   2.673E-06
   28058.01c   1.553E-04
   28060.01c   5.980E-05
   28061.01c   2.600E-06
   28062.01c   8.291E-06
   28064.01c   2.112E-06
   42092.01c   4.258E-06
   42094.01c   2.664E-06
   42095.01c   4.613E-06
   42096.01c   4.845E-06
   42097.01c   2.784E-06
   42098.01c   7.058E-06
   42100.01c   2.831E-06
   19039.01c   5.895E-03
   19040.01c   7.586E-07
   19041.01c   4.254E-04

therm h-yh2 h-yh2.84t
therm y-yh2 y-yh2.84t
mat moderator   -4.0850  moder h-yh2 1001 moder y-yh2 39089  tmp  700.0   rgb 73 64 171
39089.01c  0.357142857
1001.01c   0.642857143

% -- Fuel

div fuel sep 0
mat fuel   -10.7440 burn 1    rgb 255 0 0  tmp  700
  92235.01c    6.8794e-02
  92238.01c    2.7604e-01
  6012.01c    1.3793e-01
  8016.01c    5.1724e-01

% --- Carbon buffer layer:

mat buffer   -1.0400   rgb 0 200 200  moder grph_f 6012    tmp  700
 6012.01c    1.0

% --- Pyrolytic carbon layer:

mat PyC1     -1.8820    rgb 0 200 200  moder grph_f 6012    tmp  700
 6012.01c    1.0

% --- Pyrolytic carbon layer:

mat PyC2     -1.8820  rgb 0 200 200    moder grph_f 6012     tmp  700
 6012.01c    1.0

% --- Silicon carbide layer:

mat SiC     -3.1710  rgb 0 100 0     tmp  700
14028.01c    0.4611
14029.01c    0.0234
14030.01c    0.0154
 6012.01c    0.5

% --- Graphite matrix:
therm grph_f grph.84t
mat matrix_pin   -1.8060  rgb 220 220 220    moder grph_f 6012    tmp  700
 6012.01c  0.9999997
 5010.01c  0.0000003

therm grph_m grph.84t
mat matrix     -1.8060 rgb 200 200 200    moder grph_m 6012    tmp  700
 6012.01c  0.9999997
 5010.01c  0.0000003

therm bemet be-met.84t
mat beryllium  -1.848   rgb 0 255 0 moder bemet 4009  tmp  700
 4009.01c  1.00

mat beryllium_drum  -1.848   rgb 50 220 50 moder bemet 4009  tmp  700
 4009.01c  1.00

%
mat B4C_drum  -2.510  rgb 250 0 0   tmp  700
 5010.01c 0.76
 5011.01c 0.04
 6012.01c 0.2
%
mat B4C_central  -1.25  rgb 150 20 20   tmp  700
 5010.01c 0.76
 5011.01c 0.04
 6012.01c 0.2
% -----------------------------------------------------------------------------

% Calculation parameters:

% -----------------------------------------------------------------------------

% --- Geometry plot:
plot 3 10000 10000   21.0
plot 2 10000 10000    5.0
mesh 3 5000 5000
mesh 2 5000 5000

set mcvol 100000000
set xenon 1

% --- Libraries:
% Add your own libraries here
set acelib "sss_endfb8_serpent_yan.xsdir"  % ACE library for cross-sections
set declib "sss_endfb7.dec" % Decay library
set nfylib "sss_endfb7.nfy" % Neutron fission yield library

% --- Boundary conditions:

set bc 1 1 1

% --- Power density

set power 2e6  % W

% --- Histories:
set pop 100000 500 30

set gcu 801 802 803 805 815 816 817 810 811 820

set nfg  11  8.00000e-08 1.80000e-07 6.25000e-07 1.30000e-06 4.00000e-06
             1.48730e-04 9.11800e-03 1.83000e-01 5.00000e-01 1.35300e+00
% -----------------------------------------------------------------------------

(microreactors/mrad/steady/HPMR_dfem_griffin_ss.i)

################################################################################
## NEAMS Micro-Reactor Application Driver                                     ##
## Heat Pipe Microreactor Steady State                                        ##
## Griffin Main Application input file                                        ##
## DFEM-SN (1, 3) with CMFD acceleration                                      ##
################################################################################

[Mesh]
  [fmg]
    type = FileMeshGenerator
    file = '../mesh/gold/HPMR_OneSixth_Core_meshgenerator_tri_rotate_bdry.e'
  []
  [fmg_id]
    type = SubdomainExtraElementIDGenerator
    input = fmg
    subdomains = '200 203 100 103 301 303 10 503 600 601 201 101 400 401 250'
    extra_element_id_names = 'material_id equivalence_id'
    extra_element_ids = '815 815 802 802 801 801 803 811 820 820 817 816 805 805 820;
                         815 815 802 802 801 801 803 811 820 820 817 816 805 805 820'
  []
  [coarse_mesh]
    type = GeneratedMeshGenerator
    dim = 3
    nx = 10
    ny = 10
    nz = 10
    xmin = -0.1
    xmax = 1.1
    ymin = -0.1
    ymax = 1.2
    zmin = -0.0
    zmax = 2.1
  []
  [assign_coarse_id]
    type = CoarseMeshExtraElementIDGenerator
    input = fmg_id
    coarse_mesh = coarse_mesh
    extra_element_id_name = coarse_element_id
  []
  uniform_refine = 0
[]

[Executioner]
  type = SweepUpdate

  richardson_abs_tol = 1e-6 #-8
  richardson_rel_tol = 1e-8
  richardson_value = eigenvalue

  richardson_max_its = 1000
  inner_solve_type = GMRes
  max_inner_its = 20

  fixed_point_max_its = 1
  custom_pp = eigenvalue
  custom_rel_tol = 1e-6
  force_fixed_point_solve = true

  cmfd_acceleration = true
  coarse_element_id = coarse_element_id
  prolongation_type = multiplicative
  max_diffusion_coefficient = 1
[]

[TransportSystems]
  particle = neutron
  equation_type = eigenvalue

  G = 11
  VacuumBoundary = '10000 2000 3000'
  ReflectingBoundary = '147'

  [sn]
    scheme = DFEM-SN
    family = MONOMIAL
    order = first
    AQtype = Gauss-Chebyshev
    NPolar = 1
    NAzmthl = 3
    NA = 2
    sweep_type = asynchronous_parallel_sweeper
    using_array_variable = true
    collapse_scattering = true
    n_delay_groups = 6
  []
[]

[AuxVariables]
  [Tf]
    initial_condition = 873.15
    order = CONSTANT
    family = MONOMIAL
  []
  [Tm]
    initial_condition = 873.15
    order = CONSTANT
    family = MONOMIAL
  []
[]

[GlobalParams]
  library_file = '../isoxml/fullcore_xml_G11_endfb8_ss_tr.xml'
  library_name = fullcore_xml_G11_endfb8_ss_tr
  isotopes = 'pseudo'
  densities = 1.0
  is_meter = true
  plus = true
  dbgmat = false
  grid_names = 'Tfuel Tmod'
  grid_variables = 'Tf Tm'
[]

[Materials]
  [mod]
    type = CoupledFeedbackMatIDNeutronicsMaterial
    block = '200 203 100 103 301 303 10 503 600 601 201 101 400 401 250'
  []
[]

[PowerDensity]
  power = 345.6e3
  power_density_variable = power_density
  integrated_power_postprocessor = integrated_power
[]

[MultiApps]
  [bison]
    type = FullSolveMultiApp
    app_type = BisonApp
    positions = '0 0 0'
    input_files = HPMR_thermo_ss.i
    execute_on = 'timestep_end'
    keep_solution_during_restore = true
  []
[]

[Transfers]
  [to_sub_power_density]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    to_multi_app = bison
    source_variable = power_density
    variable = power_density
    from_postprocessors_to_be_preserved = integrated_power
    to_postprocessors_to_be_preserved = power

    search_value_conflicts = false
  []
  [from_sub_temp_fuel]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = bison
    variable = Tf
    source_variable = Tfuel
    execute_on = 'initial timestep_end'
    displaced_source_mesh = false
    displaced_target_mesh = false
    use_displaced_mesh = false

    search_value_conflicts = false
  []
  [from_sub_temp_mod]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = bison
    variable = Tm
    source_variable = Tmod
    execute_on = 'initial timestep_end'
    displaced_source_mesh = false
    displaced_target_mesh = false
    use_displaced_mesh = false

    search_value_conflicts = false
  []
[]

[UserObjects]
  [ss]
    type = TransportSolutionVectorFile
    transport_system = sn
    writing = true
    execute_on = final
  []
[]

[Postprocessors]
  [scaled_power_avg]
    type = ElementAverageValue
    block = 'fuel_quad fuel_tri'
    variable = power_density
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_avg]
    type = ElementAverageValue
    variable = Tf
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tf
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_min]
    type = ElementExtremeValue
    value_type = min
    variable = Tf
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [mod_temp_avg]
    type = ElementAverageValue
    variable = Tm
    block = 'moderator_quad moderator_tri'
    execute_on = 'initial timestep_end'
  []
  [mod_temp_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tm
    block = 'moderator_quad moderator_tri'
    execute_on = 'initial timestep_end'
  []
  [mod_temp_min]
    type = ElementExtremeValue
    value_type = min
    variable = Tm
    block = 'moderator_quad moderator_tri'
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  csv = true
  exodus = true
  perf_graph = true
[]

(microreactors/mrad/steady/HPMR_thermo_ss.i)

################################################################################
## NEAMS Micro-Reactor Application Driver                                     ##
## Heat Pipe Microreactor Steady State                                        ##
## BISON Child Application input file                                         ##
## Thermal Only Physics                                                       ##
################################################################################

R_clad_o = 0.0105 # heat pipe outer radius
R_hp_hole = 0.0107 # heat pipe + gap
num_sides = 12 # number of sides of heat pipe as a result of mesh polygonization
alpha = '${fparse 2 * pi / num_sides}'
perimeter_correction = '${fparse 0.5 * alpha / sin(0.5 * alpha)}' # polygonization correction factor for perimeter
area_correction = '${fparse sqrt(alpha / sin(alpha))}' # polygonization correction factor for area
corr_factor = '${fparse R_hp_hole / R_clad_o * area_correction / perimeter_correction}'

[GlobalParams]
  flux_conversion_factor = 1
[]

[Mesh]
  [fmg]
    type = FileMeshGenerator
    file = '../mesh/gold/HPMR_OneSixth_Core_meshgenerator_tri_rotate_bdry_fine.e'
  []
[]

[Variables]
  [temp]
    initial_condition = 800
  []
[]

[Kernels]
  # The time derivative kernel has been removed
  # We are obtaining a steady-state solution and it's preferrable not to march in time
  # to steady state as marching is expensive.
  [heat_conduction]
    type = HeatConduction
    variable = temp
  []
  [heat_source_fuel]
    type = CoupledForce
    variable = temp
    block = 'fuel_quad fuel_tri'
    v = power_density
  []
[]

[AuxVariables]
  [power_density]
    block = 'fuel_quad fuel_tri'
    family = L2_LAGRANGE
    order = FIRST
    initial_condition = 3.4e6
  []
  [Tfuel]
    block = 'fuel_quad fuel_tri'
  []
  [Tmod]
    block = 'moderator_quad moderator_tri'
  []
  [fuel_thermal_conductivity]
    block = 'fuel_quad fuel_tri'
    order = CONSTANT
    family = MONOMIAL
  []
  [fuel_specific_heat]
    block = 'fuel_quad fuel_tri'
    order = CONSTANT
    family = MONOMIAL
  []
  [monolith_thermal_conductivity]
    block = 'monolith '
    order = CONSTANT
    family = MONOMIAL
  []
  [monolith_specific_heat]
    block = 'monolith '
    order = CONSTANT
    family = MONOMIAL
  []
  [flux_uo] #auxvariable to hold heat pipe surface flux from UserObject
    block = 'reflector_quad monolith'
  []
  [flux_uo_corr] #auxvariable to hold corrected flux_uo
    block = 'reflector_quad monolith'
  []
  [hp_temp_aux]
    block = 'reflector_quad monolith'
  []
[]

[AuxKernels]
  [assign_tfuel]
    type = NormalizationAux
    variable = Tfuel
    source_variable = temp
    execute_on = 'timestep_end'
  []
  [assign_tmod]
    type = NormalizationAux
    variable = Tmod
    source_variable = temp
    execute_on = 'timestep_end'
  []
  [fuel_thermal_conductivity]
    type = MaterialRealAux
    variable = fuel_thermal_conductivity
    property = thermal_conductivity
    execute_on = timestep_end
  []
  [fuel_specific_heat]
    type = MaterialRealAux
    variable = fuel_specific_heat
    property = specific_heat
    execute_on = timestep_end
  []
  [monolith_thermal_conductivity]
    type = MaterialRealAux
    variable = monolith_thermal_conductivity
    property = thermal_conductivity
    execute_on = timestep_end
  []
  [monolith_specific_heat]
    type = MaterialRealAux
    variable = monolith_specific_heat
    property = specific_heat
    execute_on = timestep_end
  []
  [flux_uo]
    type = SpatialUserObjectAux
    variable = flux_uo
    user_object = flux_uo
  []
  [flux_uo_corr]
    type = NormalizationAux
    variable = flux_uo_corr
    source_variable = flux_uo
    normal_factor = '${fparse corr_factor}'
  []
[]

[BCs]
  [outside_bc]
    type = ConvectiveFluxFunction # (Robin BC)
    variable = temp
    boundary = 'bottom top side'
    coefficient = 1e3 # W/K/m^2
    T_infinity = 800 # K air temperature at the top of the core
  []
  [hp_temp]
    type = CoupledConvectiveHeatFluxBC
    boundary = 'heat_pipe_ht_surf'
    variable = temp
    T_infinity = hp_temp_aux
    htc = 750
  []
[]

[Materials]
  [fuel_matrix_thermal]
    type = GraphiteMatrixThermal
    block = 'fuel_quad fuel_tri'
    graphite_grade = IG_110
    packing_fraction = 0.4
    specific_heat_scale_factor = 1.0
    thermal_conductivity_scale_factor = 1.0
    fast_neutron_fluence = 0 #6.75E+24 # this value is nuetron fluence over 0.1MeV
    temperature = temp
  []
  [monolith_matrix_thermal]
    type = GraphiteMatrixThermal
    block = 'monolith '
    graphite_grade = IG_110
    packing_fraction = 0
    specific_heat_scale_factor = 1.0
    thermal_conductivity_scale_factor = 1.0
    fast_neutron_fluence = 0 #6.75E+24 # this value is nuetron fluence over 0.1MeV
    temperature = temp
  []
  [moderator_thermal]
    type = HeatConductionMaterial
    block = 'moderator_quad moderator_tri'
    temp = temp
    thermal_conductivity = 20 # W/m/K
    specific_heat = 500 # random value
  []
  [airgap_thermal]
    type = HeatConductionMaterial
    block = 'air_gap_tri air_gap_quad outer_shield' # Helium gap
    temp = temp
    thermal_conductivity = 0.15 # W/m/K
    specific_heat = 5197 # random value
  []
  [axial_reflector_thermal]
    type = HeatConductionMaterial
    block = 'reflector_tri reflector_quad'
    temp = temp
    thermal_conductivity = 199 # W/m/K
    specific_heat = 1867 # random value
  []
  [B4C_thermal]
    type = HeatConductionMaterial
    block = 'B4C'
    temp = temp
    thermal_conductivity = 92 # W/m/K
    specific_heat = 960 # random value
  []
  [SS_thermal]
    type = SS316Thermal
    temperature = temp
    block = mod_ss
  []
  [fuel_density]
    type = Density
    block = 'fuel_quad fuel_tri'
    density = 2276.5
  []
  [moderator_density]
    type = Density
    block = 'moderator_quad moderator_tri'
    density = 4.3e3
  []
  [monolith_density]
    type = Density
    block = 'monolith'
    density = 1806
  []
  [airgap_density]
    type = Density
    block = 'air_gap_tri air_gap_quad outer_shield' #helium
    density = 180
  []
  [axial_reflector_density]
    type = Density
    block = 'reflector_tri reflector_quad'
    density = 1848
  []
  [B4C_density]
    type = Density
    block = B4C
    density = 2510
  []
  [SS_density]
    type = Density
    density = 7990
    block = mod_ss
  []
[]

[MultiApps]
  [sockeye]
    type = TransientMultiApp
    app_type = SockeyeApp
    positions_file = 'hp_centers.txt'
    input_files = 'HPMR_sockeye_ss.i'
    execute_on = 'timestep_begin' # execute on timestep begin because hard to have a good initial guess on heat flux
    max_procs_per_app = 1
    output_in_position = true
    sub_cycling = true
  []
[]

[Transfers]
  [from_sockeye_temp]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = sockeye
    source_variable = hp_temp_aux
    variable = hp_temp_aux
    execute_on = 'timestep_begin'
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [to_sockeye_flux]
    type = MultiAppGeneralFieldUserObjectTransfer
    variable = master_flux
    to_multi_app = sockeye
    execute_on = 'timestep_begin'
    source_user_object = flux_uo
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
[]

[UserObjects]
  [flux_uo]
    type = NearestPointLayeredSideDiffusiveFluxAverage
    direction = z
    num_layers = 100
    points_file = 'hp_centers.txt'
    variable = temp
    diffusivity = thermal_conductivity
    execute_on = linear
    boundary = 'heat_pipe_ht_surf'
  []
[]

[Executioner]
  type = Transient

  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart '
  petsc_options_value = 'hypre boomeramg 100'

  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-8

  start_time = -5e4 # negative start time so we can start running from t = 0
  end_time = 0
  dtmin = 1
  # Multiple steps are needed to allow iteration between the thermo and heat pipe applications
  num_steps = 5
[]

[Postprocessors]
  [hp_heat_integral]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 'heat_pipe_ht_surf'
    diffusivity = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [hp_end_heat_integral]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 'heat_pipe_ht_surf_bot'
    diffusivity = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [ext_side_integral]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 'side'
    diffusivity = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [mirror_side_integral]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 'side_mirror'
    diffusivity = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [tb_integral]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 'top bottom'
    diffusivity = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [total_sink]
    type = SumPostprocessor
    values = 'hp_heat_integral hp_end_heat_integral ext_side_integral mirror_side_integral tb_integral'
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_avg]
    type = ElementAverageValue
    variable = temp
    block = 'fuel_quad fuel_tri'
  []
  [fuel_temp_max]
    type = ElementExtremeValue
    variable = temp
    block = 'fuel_quad fuel_tri'
  []
  [fuel_temp_min]
    type = ElementExtremeValue
    variable = temp
    block = 'fuel_quad fuel_tri'
    value_type = min
  []
  [mod_temp_avg]
    type = ElementAverageValue
    variable = temp
    block = 'moderator_quad moderator_tri'
  []
  [mod_temp_max]
    type = ElementExtremeValue
    variable = temp
    block = 'moderator_quad moderator_tri'
  []
  [mod_temp_min]
    type = ElementExtremeValue
    variable = temp
    block = 'moderator_quad moderator_tri'
    value_type = min
  []
  [monolith_temp_avg]
    type = ElementAverageValue
    variable = temp
    block = monolith
  []
  [monolith_temp_max]
    type = ElementExtremeValue
    variable = temp
    block = monolith
  []
  [monolith_temp_min]
    type = ElementExtremeValue
    variable = temp
    block = monolith
    value_type = min
  []
  [heatpipe_surface_temp_avg]
    type = SideAverageValue
    variable = temp
    boundary = heat_pipe_ht_surf
  []
  [power]
    type = ElementIntegralVariablePostprocessor
    block = 'fuel_quad fuel_tri'
    variable = power_density
    execute_on = 'initial timestep_end transfer'
  []
  [fuel_cP]
    type = ElementExtremeValue
    value_type = 'max'
    variable = fuel_specific_heat
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [fuel_k]
    type = ElementExtremeValue
    value_type = 'max'
    variable = fuel_thermal_conductivity
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [monolith_cP]
    type = ElementExtremeValue
    value_type = 'max'
    variable = monolith_specific_heat
    block = 'monolith '
    execute_on = 'initial timestep_end'
  []
  [monolith_k]
    type = ElementExtremeValue
    value_type = 'max'
    variable = monolith_thermal_conductivity
    block = 'monolith '
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  perf_graph = true
  exodus = true
  color = true
  csv = true
  checkpoint = true
[]

(microreactors/mrad/steady/HPMR_sockeye_ss.i)

################################################################################
## NEAMS Micro-Reactor Application Driver                                     ##
## Heat Pipe Microreactor Steady State                                        ##
## Sockeye Grandchild Application input file                                  ##
## Effective Heat Conduction Model with Operation Limits                      ##
################################################################################

# Total heat removed/added to heat pipe
Q_hp = 1800.

# Wick characteristics
R_pore = 15.0e-6
D_h_pore = '${fparse 2.0 * R_pore}'
permeability = 2e-9
porosity = 0.70

# Envelope ("env")
# SS316. Incropera & DeWitt, 3rd ed, Table A.1 @ 900K (627C)
# Density (kg/m3)
rho_env = 8238.
# Thermal conductivity (W/m-K)
k_env = 23.05
# Specific heat capacity (J/kg-K)
cp_env = 589.

# Potassium vapor
# From Appendix C of text book "Heat Pipe Design and Technology"
# Sat. Vapor Potassium at T = 650C (923 K)
# Density (kg/m3)
rho_vapor = 0.193
# Effective "super" conductivity (W/m-K)
#k_vapor = 1.0e5 # Reaches 600 seconds
k_vapor = 1.0e6
# Specific heat capacity (J/kg-K)
cp_vapor = 5320.

# Potassium liquid
# From Appendix C of text book "Heat Pipe Design and Technology"
# Sat. Liquid Potassium at T = 650C (923 K)
# Density (kg/m3)
rho_liquid = 695.4
# Thermal conductivity (W/m-K)
k_liquid = 40.08
# Specific heat capacity (J/kg-K)
# From Table 1.1, no temperature data given
cp_liquid = 810.
# Melting point, Table 3.1, lists 62C, rounding up
# T_melting = 340.

# Wick, homogenize envelope and fluid
# Density (kg/m3)
rho_wick = '${fparse porosity * rho_liquid + (1.0 - porosity) * rho_env}'
# Thermal conductivity (W/m-K)
k_wick = '${fparse porosity * k_liquid + (1.0 - porosity) * k_env}'
# Specific heat capacity (J/kg-K)
# From Table 1.1, no temperature data given
cp_wick = '${fparse porosity * cp_liquid + (1.0 - porosity) * cp_env}'

# Elevations and lengths
# Note: For blackbox model -- manually update "length" input
length_evap = 180.0e-2
length_adia = 30.0e-2
length_cond = 90.0e-2

# Mesh density
# The dimensions are nicely divisible by 3 cm mesh.
nelem_base_evap = 50
nelem_base_adia = 10
nelem_base_cond = 30
mesh_density = 3
nelem_evap = '${fparse mesh_density*nelem_base_evap}'
nelem_adia = '${fparse mesh_density*nelem_base_adia}'
nelem_cond = '${fparse mesh_density*nelem_base_cond}'

# Envelope thickness
t_env = 0.08e-2
# Liquid annulus thickness
t_ann = 0.07e-2
# Wick thickness
t_wick = 0.1e-2

# Radial geometry
# Envelope outer
R_hp_o = 1.05e-2
D_hp_o = '${fparse 2.0 * R_hp_o}'
# Inner Envelope/outer annulus
R_hp_i = '${fparse R_hp_o - t_env}'
D_hp_i = '${fparse 2.0 * R_hp_i}'
# Inner annulus/wick outer
R_wick_o = '${fparse R_hp_i - t_ann}'
D_wick_o = '${fparse 2.0 * R_wick_o}'
# Inner wick/vapor core outer
R_wick_i = '${fparse R_wick_o - t_wick}'
D_wick_i = '${fparse 2.0 * R_wick_i}'

# BCs for condenser
T_ext_cond = 800.
htc_ext_cond = 1.0e6

# Evaporator parameters
S_evap = '${fparse pi * D_hp_o * length_evap}'
q_evap = '${fparse Q_hp / S_evap}'

[GlobalParams]
  scaling_factor_temperature = 1e-2
  fp_2phase = fp_2phase
[]

[FluidProperties]
  [fp_2phase]
    type = PotassiumTwoPhaseFluidProperties
    emit_on_nan = none
  []
[]

[Components]
  [hp]
    type = HeatPipeConduction

    # Common to both HeatPipe2Phase and HeatPipeBlackbox
    position = '0 0 0'
    orientation = '0 0 1'
    length = '${length_evap} ${length_adia} ${length_cond}'
    n_elems = '${nelem_evap} ${nelem_adia} ${nelem_cond}'
    gravity_vector = '0 0 -9.8'
    D_wick_i = ${D_wick_i}
    D_wick_o = ${D_wick_o}
    R_pore = ${R_pore}
    porosity = ${porosity}
    permeability = ${permeability}

    # HeatPipeConduction only
    # Axial dimensions (for heat transfer & analytic limits)
    axial_region_names = 'evap adia cond'
    L_evap = ${length_evap}
    L_adia = ${length_adia}
    L_cond = ${length_cond}
    # Radial dimensions, mesh, and materials for heat transfer problem
    D_clad_o = ${D_hp_o}
    D_clad_i = ${D_hp_i}
    D_h_pore = ${D_h_pore}
    # Mesh
    n_elems_clad = 4
    n_elems_wick = 8
    n_elems_core = 10

    # Thermal properties
    sp_vapor = sp_vapor
    sp_liquid = sp_wick_ann
    sp_wick = sp_wick_ann
    sp_clad = sp_clad
    k_core = ${k_vapor}
    k_eff = ${k_wick}

    fp_2phase = fp_2phase
    evaporator_at_start_end = true
    # Initial temperature of block
    initial_T = ${T_ext_cond}
    T_ref = T_inner_avg

    # To evaluate the constant properties
    T_ref_density = 1000

    make_pressure_corrections = true
  []

  [condenser_boundary]
    type = HSBoundaryAmbientConvection
    boundary = 'hp:cond:outer'
    hs = hp
    T_ambient = ${T_ext_cond}
    htc_ambient = ${htc_ext_cond} #large value to approach an effective DirichletBC
    scale = 1
  []
  [evaporator_boundary]
    type = HSBoundaryExternalAppConvection
    boundary = 'hp:evap:outer'
    hs = hp
    T_ext = virtual_Text
    htc_ext = virtual_htc
  []
[]

[SolidProperties]
  [sp_vapor]
    type = ThermalFunctionSolidProperties
    rho = ${rho_vapor}
    cp = ${cp_vapor}
    k = ${k_vapor}
  []
  [sp_wick_ann]
    type = ThermalFunctionSolidProperties
    rho = ${rho_wick}
    cp = ${cp_wick}
    k = ${k_wick}
  []
  [sp_clad]
    type = ThermalFunctionSolidProperties
    rho = ${rho_env}
    cp = ${cp_env}
    k = ${k_env}
  []
[]

[UserObjects]
  [surf_T]
    type = LayeredSideAverage
    direction = z
    num_layers = 100
    variable = T_solid
    boundary = 'hp:evap:outer'
  []
[]

[AuxKernels]
  [hp_var]
    type = SpatialUserObjectAux
    variable = hp_temp_aux
    user_object = surf_T
  []
  [virtual_Text]
    type = ParsedAux
    variable = virtual_Text
    coupled_variables = 'T_solid master_flux virtual_htc'
    expression = 'master_flux/virtual_htc + T_solid'
  []
[]

[Functions]
  [scale_fcn]
    type = ParsedFunction
    symbol_names = 'catastrophic_pp recoverable_pp operational_pp'
    symbol_values = 'catastrophic_pp recoverable_pp operational_pp'
    expression = 'catastrophic_pp*recoverable_pp*operational_pp'
  []
[]

[AuxVariables]
  [T_wall_var]
    initial_condition = ${T_ext_cond}
  []
  [operational_aux]
    initial_condition = 1
  []
  [master_flux]
    initial_condition = ${q_evap}
  []
  [hp_temp_aux]
    initial_condition = ${T_ext_cond}
  []
  [virtual_Text]
    initial_condition = ${T_ext_cond}
  []
  [virtual_htc]
    initial_condition = 1.0
  []
[]

[Postprocessors]
  [Integral_BC_Total]
    type = SumPostprocessor
    values = 'condenser_boundary_integral evaporator_boundary_integral'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [ZeroPP]
    type = EmptyPostprocessor
  []
  [Integral_BC_Cond]
    type = DifferencePostprocessor
    value1 = ZeroPP
    value2 = condenser_boundary_integral
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Integral_BC_RelErr]
    type = RelativeDifferencePostprocessor
    value1 = evaporator_boundary_integral
    value2 = Integral_BC_Cond
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [operational_pp]
    type = ElementAverageValue
    variable = operational_aux
    execute_on = 'initial timestep_begin TIMESTEP_END'
  []
  # set `catastrophic_pp` as if it is recoverable for a solving to steady-state simulation
  # this MUST be changed back to catestrophic in transient simulations
  [catastrophic_pp]
    type = HeatRemovalRateLimitScale
    heat_addition_pps = 'evaporator_boundary_integral'
    limit_condenser_side = false
    recoverable_heat_removal_limit_pps = 'hp_boiling_limit hp_capillary_limit hp_entrainment_limit'
    catastrophic_heat_removal_limit_pps = ''
    T = T_inner_avg
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [recoverable_pp]
    type = HeatRemovalRateLimitScale
    heat_addition_pps = 'evaporator_boundary_integral'
    limit_condenser_side = false
    catastrophic_heat_removal_limit_pps = ''
    recoverable_heat_removal_limit_pps = 'hp_sonic_limit hp_viscous_limit'
    T = T_inner_avg
    execute_on = 'INITIAL linear nonlinear TIMESTEP_END'
  []
  [T_evap_inner]
    type = SideAverageValue
    boundary = hp:evap:inner
    variable = T_solid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_cond_inner]
    type = SideAverageValue
    boundary = hp:cond:inner
    variable = T_solid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_evap_outer]
    type = SideAverageValue
    boundary = hp:evap:outer
    variable = T_solid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_cond_outer]
    type = SideAverageValue
    boundary = hp:cond:outer
    variable = T_solid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_wall_var_avg]
    type = ElementAverageValue
    variable = T_wall_var
    execute_on = 'Initial timestep_end'
  []
  [T_inner_avg]
    type = SideAverageValue
    variable = T_solid
    boundary = hp:inner
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_inner_max]
    type = NodalExtremeValue
    variable = T_solid
    boundary = hp:inner
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_inner_min]
    type = NodalExtremeValue
    variable = T_solid
    boundary = hp:inner
    execute_on = 'INITIAL TIMESTEP_END'
    value_type = min
  []
  [DT_outer]
    type = DifferencePostprocessor
    value1 = T_evap_outer
    value2 = T_cond_outer
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [DT_inner]
    type = DifferencePostprocessor
    value1 = T_evap_inner
    value2 = T_cond_inner
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [A_int_master_flux]
    type = SideIntegralVariablePostprocessor
    variable = master_flux
    boundary = 'hp:evap:inner'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [A_int_T_ext]
    type = SideIntegralVariablePostprocessor
    variable = virtual_Text
    boundary = 'hp:evap:inner'
    execute_on = 'INITIAL LINEAR'
  []
  [A_avg_T_aux]
    type = AverageNodalVariableValue
    variable = hp_temp_aux
    boundary = 'hp:evap:inner'
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[VectorPostprocessors]
  [env_vpp]
    type = NodalValueSampler
    variable = T_solid
    block = 'hp:clad'
    sort_by = x
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

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

[Executioner]
  type = Transient

  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  scheme = bdf2
  line_search = none

  # ensure nl_abs_tol >> nl_rel_tol
  nl_abs_tol = 1e-6
  nl_rel_tol = 1e-8
  nl_max_its = 100

  l_tol = 1e-3
  l_max_its = 100

  start_time = -5e4 # negative start time so we can start running from t = 0
  end_time = 0
  dtmin = 1
  dt = 5000
[]

[Outputs]
  [console]
    type = Console
    max_rows = 5
    execute_postprocessors_on = 'INITIAL TIMESTEP_END FINAL FAILED'
  []
  [csv]
    type = CSV
    execute_on = 'INITIAL TIMESTEP_END FINAL FAILED'
    execute_vector_postprocessors_on = 'FINAL FAILED'
  []
[]

(microreactors/mrad/mesh/HPMR_OneSixth_Core_meshgenerator_tri.i)

################################################################
# This is the MOOSE input file to generate the mesh that can be
# used for the 1/6 core Heat-Pipe Microreactor Multiphysics
# simulations (Griffin neutronics simulation).
# Running this input requires MOOSE Reactor Module Objects
# Users should use
# --mesh-only HPMR_OneSixth_Core_meshgenerator_tri_rotate_bdry.e
# command line argument to generate
# exodus file for further Multiphysics simulations.
################################################################
[GlobalParams]
  create_outward_interface_boundaries = false
[]

[Mesh]
  [BatchMeshGeneratorAction]
    # We generate each control drum in 2 steps (concentric circles generation and azimuthal absorber/reflector splitting),
    # this syntax allows us to input the control drum (and reflector) generation in a very compact manner.
    [control_drums_ccg]
      mesh_generator_type = 'HexagonConcentricCircleAdaptiveBoundaryMeshGenerator'
      mesh_name_prefix = 'cd0'
      multi_batch_params_method = 'corresponding'
      # The main different parameters between the control drums are how the meshes are adapted to the reactor core
      # The batch parameters here are used to account for these differences
      batch_vector_input_param_names = 'meshes_to_adapt_to sides_to_adapt'
      batch_vector_input_param_types = 'MGNAME UINT'
      batch_vector_input_param_values = 'Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned|
                                         2 3 4;2 3;1 2 3;1 2;0 1 2;0 1;
                                         0 1 5;0 5;0 4 5;4 5;3 4 5;3 4'
      # The rest of the control drum parameters are the same for all control drums
      fixed_vector_input_param_names = 'num_sectors_per_side background_block_ids ring_radii ring_intervals ring_block_ids'
      fixed_vector_input_param_types = 'UINT USHORT REAL UINT USHORT'
      fixed_vector_input_param_values = '4 4 4 4 4 4;
                                         504;
                                         12.25 13.25;
                                         2 1;
                                         500 501 502'
      fixed_scalar_input_param_names = 'hexagon_size background_intervals preserve_volumes is_control_drum create_outward_interface_boundaries'
      fixed_scalar_input_param_types = 'REAL UINT BOOL BOOL BOOL'
      fixed_scalar_input_param_values = '13.376 2 true true false'
    []
    [control_drums_split]
      mesh_generator_type = 'AzimuthalBlockSplitGenerator'
      mesh_name_prefix = 'cd'
      multi_batch_params_method = 'corresponding'
      # Based on the locations of the control drums, we need to split the blocks in different ways
      batch_scalar_input_param_names = 'input start_angle'
      batch_scalar_input_param_types = 'MGNAME REAL'
      batch_scalar_input_param_values = 'cd0_0 cd0_1 cd0_2 cd0_3 cd0_4 cd0_5 cd0_6 cd0_7 cd0_8 cd0_9 cd0_10 cd0_11;
                                         15 345 315 285 255 225 195 165 135 105 75 45'
      fixed_vector_input_param_names = 'old_blocks new_block_ids'
      fixed_vector_input_param_types = 'SDNAME USHORT'
      fixed_vector_input_param_values = '502;503'
      fixed_scalar_input_param_names = 'angle_range'
      fixed_scalar_input_param_types = 'REAL'
      fixed_scalar_input_param_values = '90'
    []
    [reflector_blocks]
      mesh_generator_type = 'HexagonConcentricCircleAdaptiveBoundaryMeshGenerator'
      mesh_name_prefix = 'ref'
      multi_batch_params_method = 'corresponding'
      # Similar to the control drums, we need to adapt the reflector meshes to the reactor core based on the locations
      batch_vector_input_param_names = 'sides_to_adapt'
      batch_vector_input_param_types = 'UINT'
      batch_vector_input_param_values = '0;1;2;3;4;5'
      fixed_vector_input_param_names = 'meshes_to_adapt_to num_sectors_per_side background_block_ids'
      fixed_vector_input_param_types = 'MGNAME UINT USHORT'
      fixed_vector_input_param_values = 'Patterned;
                                         4 4 4 4 4 4;
                                         400 401'
      fixed_scalar_input_param_names = 'hexagon_size background_intervals create_outward_interface_boundaries'
      fixed_scalar_input_param_types = 'REAL UINT BOOL'
      fixed_scalar_input_param_values = '13.376 2 false'
    []
  []
  # Moderator pins
  [Mod_hex]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2 '
    background_intervals = 1
    background_block_ids = '10'
    polygon_size = 1.15
    polygon_size_style ='apothem'
    ring_radii = '0.825 0.92'
    ring_intervals = '2 1'
    ring_block_ids = '103 100 101' # 103 is tri mesh
    preserve_volumes = on
    quad_center_elements = false
  []
  # Heat Pipe pins
  [HP_hex]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2 '
    background_intervals = 1
    background_block_ids = '10'
    polygon_size = 1.15
    polygon_size_style ='apothem'
    ring_radii = '0.97 1.07'
    ring_intervals = '2 1'
    ring_block_ids = '203 200 201' # 203 is tri mesh
    preserve_volumes = on
    quad_center_elements = false
  []
  # Fuel pins
  [Fuel_hex]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2 '
    background_intervals = 1
    background_block_ids = '10'
    polygon_size = 1.15
    polygon_size_style ='apothem'
    ring_radii = '1'
    ring_intervals = '2'
    ring_block_ids = '303 301'  # 303 is tri mesh
    preserve_volumes = on
    quad_center_elements = false
  []
  # Assembly mesh with all the three types of pins
  [Patterned]
    type = PatternedHexMeshGenerator
    inputs = 'Fuel_hex HP_hex Mod_hex'
    hexagon_size = 13.376
    background_block_id = 10
    background_intervals = 1
    pattern = '1 0 1 0 1 0 1;
              0 2 0 2 0 2 0 0;
             1 0 1 0 1 0 1 2 1;
            0 2 0 2 0 2 0 0 0 0;
           1 0 1 0 1 0 1 2 1 2 1;
          0 2 0 2 0 2 0 0 0 0 0 0;
         1 0 1 0 1 0 1 2 1 2 1 2 1;
          0 2 0 2 0 2 0 0 0 0 0 0;
           1 0 1 0 1 0 1 2 1 2 1;
            0 2 0 2 0 2 0 0 0 0;
             1 0 1 0 1 0 1 2 1;
              0 2 0 2 0 2 0 0;
               1 0 1 0 1 0 1'
  []
  # Central void region
  [air_center]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    num_sectors_per_side= '4 4 4 4 4 4'
    meshes_to_adapt_to = 'Patterned Patterned Patterned Patterned Patterned Patterned'
    sides_to_adapt = '0 1 2 3 4 5'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = '600 601'
  []
  # Dummy assembly to help form a regular hexagonal core pattern
  # This will be deleted later
  [dummy]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    num_sectors_per_side= '4 4 4 4 4 4'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = '700 701'
       # external_boundary_id = 9998
  []
  # Stitching assemblies together to form the core mesh
  [core]
    type = PatternedHexMeshGenerator
    inputs = 'Patterned cd_0 cd_1 cd_2 cd_3 cd_4 cd_5 cd_6 cd_7 cd_8 cd_9 cd_10 cd_11 ref_0 ref_1 ref_2 ref_3 ref_4 ref_5 dummy air_center'
    # Pattern ID  0       1    2    3    4    5    6    7    8    9    10    11    12    13    14    15    16    17    18   19       20
    pattern_boundary = none
    generate_core_metadata = true
    pattern = '19 17 12 16 19;
             17 11  0  0  1 16;
           10  0  0  0  0  0  2;
          18 0  0   0  0  0  0 15;
        19  9  0  0  20  0  0  3 19;
          18 0  0   0  0  0  0 15;
            8  0  0  0  0  0  4;
             13  7  0  0  5 14;
               19 13  6 14 19'
    rotate_angle = 60
    assign_control_drum_id = true
  []
  # Delete the dummy assemblies
  [del_dummy]
    type = BlockDeletionGenerator
    block = '700 701'
    input = core
    new_boundary = 10000
  []
  # Add peripheral circular region
  [add_outer_shield]
    type = PeripheralRingMeshGenerator
    input = del_dummy
    peripheral_layer_num = 1
    peripheral_ring_radius = 115.0
    input_mesh_external_boundary = 10000
    peripheral_ring_block_id = 250
    peripheral_ring_block_name = outer_shield
 []
 # trim the full core to 1/6 core
 [del_1]
    type = PlaneDeletionGenerator
    point = '0 0 0'
    normal = '10 17.32 0'
    input = add_outer_shield
    new_boundary = 147
  []
  [del_2]
    type = PlaneDeletionGenerator
    point = '0 0 0'
    normal = '10 -17.32 0'
    input = del_1
    new_boundary = 147
  []
  # extrusion
  [extrude]
     type = AdvancedExtruderGenerator
     input = del_2
     heights = '20 160 20'
     # Use `num_layers = '6 16 6'` for BISON mesh
     num_layers = '1 8 1'
     # biased upper and lower reflector mesh, uncomment for BISON mesh
     # biases = '1.6 1.0 0.625'
     direction = '0 0 1'
     top_boundary = 2000
     bottom_boundary = 3000
  []
  # Define some special reflector blocks
  [reflector_bottom_quad]
    type = ParsedSubdomainMeshGenerator
    input = extrude
    combinatorial_geometry = 'z<=20'
    block_id = 1000
    excluded_subdomains = '103 203 303 250 400 401 500 501 502 503 504 600 601'
  []
  [reflector_bottom_tri]
    type = ParsedSubdomainMeshGenerator
    input = reflector_bottom_quad
    combinatorial_geometry = 'z<=20'
    block_id = 1003
    excluded_subdomains = '1000 250 400 401 500 501 502 503 504 600 601'
  []
  [reflector_top_quad]
    type = ParsedSubdomainMeshGenerator
    input = reflector_bottom_tri
    combinatorial_geometry = 'z>=180'
    block_id = 1000
    excluded_subdomains = '103 203 303 200 201 250 400 401 500 501 502 503 504 600 601'
  []
  [reflector_top_tri]
    type = ParsedSubdomainMeshGenerator
    input = reflector_top_quad
    combinatorial_geometry = 'z>=180'
    block_id = 1003
    excluded_subdomains = '1000 200 201 203 250 400 401 500 501 502 503 504 600 601'
  []
  # Assgin block names
  [rename_blocks]
    type = RenameBlockGenerator
    old_block = '     101       100         103          201        200           203            301       303        400            401         10      503      600         601        504          500                501          502         1000           1003'
    new_block = ' mod_ss moderator_quad moderator_tri hp_ss  heat_pipes_quad heat_pipes_tri fuel_quad fuel_tri reflector_tri reflector_quad monolith B4C air_gap_tri air_gap_quad air_gap_quad reflector_tri  reflector_quad reflector_quad reflector_quad reflector_tri'
    input = reflector_top_tri
  []
  # Assign boundary names
  [rename_boundaries]
    type = RenameBoundaryGenerator
    input = rename_blocks
    old_boundary = ' 10000 2000 3000 147'
    new_boundary = ' side top bottom side_mirror'
  []
  # Scale to m unit for general MOOSE application
  # Note that is_meter needs to be true in Griffin
  [scale]
    type = TransformGenerator
    input = rename_boundaries
    transform = SCALE
    vector_value = '1e-2 1e-2 1e-2'
  []
  # Rotate to the desired orientation for simulation
  [rotate_end]
    type = TransformGenerator
    input = scale
    transform = ROTATE
    vector_value = '0 0 -120'
  []
[]

(microreactors/mrad/mesh/HPMR_OneSixth_Core_meshgenerator_tri.i)

################################################################
# This is the MOOSE input file to generate the mesh that can be
# used for the 1/6 core Heat-Pipe Microreactor Multiphysics
# simulations (Griffin neutronics simulation).
# Running this input requires MOOSE Reactor Module Objects
# Users should use
# --mesh-only HPMR_OneSixth_Core_meshgenerator_tri_rotate_bdry.e
# command line argument to generate
# exodus file for further Multiphysics simulations.
################################################################
[GlobalParams]
  create_outward_interface_boundaries = false
[]

[Mesh]
  [BatchMeshGeneratorAction]
    # We generate each control drum in 2 steps (concentric circles generation and azimuthal absorber/reflector splitting),
    # this syntax allows us to input the control drum (and reflector) generation in a very compact manner.
    [control_drums_ccg]
      mesh_generator_type = 'HexagonConcentricCircleAdaptiveBoundaryMeshGenerator'
      mesh_name_prefix = 'cd0'
      multi_batch_params_method = 'corresponding'
      # The main different parameters between the control drums are how the meshes are adapted to the reactor core
      # The batch parameters here are used to account for these differences
      batch_vector_input_param_names = 'meshes_to_adapt_to sides_to_adapt'
      batch_vector_input_param_types = 'MGNAME UINT'
      batch_vector_input_param_values = 'Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned;
                                         Patterned Patterned Patterned;
                                         Patterned Patterned|
                                         2 3 4;2 3;1 2 3;1 2;0 1 2;0 1;
                                         0 1 5;0 5;0 4 5;4 5;3 4 5;3 4'
      # The rest of the control drum parameters are the same for all control drums
      fixed_vector_input_param_names = 'num_sectors_per_side background_block_ids ring_radii ring_intervals ring_block_ids'
      fixed_vector_input_param_types = 'UINT USHORT REAL UINT USHORT'
      fixed_vector_input_param_values = '4 4 4 4 4 4;
                                         504;
                                         12.25 13.25;
                                         2 1;
                                         500 501 502'
      fixed_scalar_input_param_names = 'hexagon_size background_intervals preserve_volumes is_control_drum create_outward_interface_boundaries'
      fixed_scalar_input_param_types = 'REAL UINT BOOL BOOL BOOL'
      fixed_scalar_input_param_values = '13.376 2 true true false'
    []
    [control_drums_split]
      mesh_generator_type = 'AzimuthalBlockSplitGenerator'
      mesh_name_prefix = 'cd'
      multi_batch_params_method = 'corresponding'
      # Based on the locations of the control drums, we need to split the blocks in different ways
      batch_scalar_input_param_names = 'input start_angle'
      batch_scalar_input_param_types = 'MGNAME REAL'
      batch_scalar_input_param_values = 'cd0_0 cd0_1 cd0_2 cd0_3 cd0_4 cd0_5 cd0_6 cd0_7 cd0_8 cd0_9 cd0_10 cd0_11;
                                         15 345 315 285 255 225 195 165 135 105 75 45'
      fixed_vector_input_param_names = 'old_blocks new_block_ids'
      fixed_vector_input_param_types = 'SDNAME USHORT'
      fixed_vector_input_param_values = '502;503'
      fixed_scalar_input_param_names = 'angle_range'
      fixed_scalar_input_param_types = 'REAL'
      fixed_scalar_input_param_values = '90'
    []
    [reflector_blocks]
      mesh_generator_type = 'HexagonConcentricCircleAdaptiveBoundaryMeshGenerator'
      mesh_name_prefix = 'ref'
      multi_batch_params_method = 'corresponding'
      # Similar to the control drums, we need to adapt the reflector meshes to the reactor core based on the locations
      batch_vector_input_param_names = 'sides_to_adapt'
      batch_vector_input_param_types = 'UINT'
      batch_vector_input_param_values = '0;1;2;3;4;5'
      fixed_vector_input_param_names = 'meshes_to_adapt_to num_sectors_per_side background_block_ids'
      fixed_vector_input_param_types = 'MGNAME UINT USHORT'
      fixed_vector_input_param_values = 'Patterned;
                                         4 4 4 4 4 4;
                                         400 401'
      fixed_scalar_input_param_names = 'hexagon_size background_intervals create_outward_interface_boundaries'
      fixed_scalar_input_param_types = 'REAL UINT BOOL'
      fixed_scalar_input_param_values = '13.376 2 false'
    []
  []
  # Moderator pins
  [Mod_hex]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2 '
    background_intervals = 1
    background_block_ids = '10'
    polygon_size = 1.15
    polygon_size_style ='apothem'
    ring_radii = '0.825 0.92'
    ring_intervals = '2 1'
    ring_block_ids = '103 100 101' # 103 is tri mesh
    preserve_volumes = on
    quad_center_elements = false
  []
  # Heat Pipe pins
  [HP_hex]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2 '
    background_intervals = 1
    background_block_ids = '10'
    polygon_size = 1.15
    polygon_size_style ='apothem'
    ring_radii = '0.97 1.07'
    ring_intervals = '2 1'
    ring_block_ids = '203 200 201' # 203 is tri mesh
    preserve_volumes = on
    quad_center_elements = false
  []
  # Fuel pins
  [Fuel_hex]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2 '
    background_intervals = 1
    background_block_ids = '10'
    polygon_size = 1.15
    polygon_size_style ='apothem'
    ring_radii = '1'
    ring_intervals = '2'
    ring_block_ids = '303 301'  # 303 is tri mesh
    preserve_volumes = on
    quad_center_elements = false
  []
  # Assembly mesh with all the three types of pins
  [Patterned]
    type = PatternedHexMeshGenerator
    inputs = 'Fuel_hex HP_hex Mod_hex'
    hexagon_size = 13.376
    background_block_id = 10
    background_intervals = 1
    pattern = '1 0 1 0 1 0 1;
              0 2 0 2 0 2 0 0;
             1 0 1 0 1 0 1 2 1;
            0 2 0 2 0 2 0 0 0 0;
           1 0 1 0 1 0 1 2 1 2 1;
          0 2 0 2 0 2 0 0 0 0 0 0;
         1 0 1 0 1 0 1 2 1 2 1 2 1;
          0 2 0 2 0 2 0 0 0 0 0 0;
           1 0 1 0 1 0 1 2 1 2 1;
            0 2 0 2 0 2 0 0 0 0;
             1 0 1 0 1 0 1 2 1;
              0 2 0 2 0 2 0 0;
               1 0 1 0 1 0 1'
  []
  # Central void region
  [air_center]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    num_sectors_per_side= '4 4 4 4 4 4'
    meshes_to_adapt_to = 'Patterned Patterned Patterned Patterned Patterned Patterned'
    sides_to_adapt = '0 1 2 3 4 5'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = '600 601'
  []
  # Dummy assembly to help form a regular hexagonal core pattern
  # This will be deleted later
  [dummy]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    num_sectors_per_side= '4 4 4 4 4 4'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = '700 701'
       # external_boundary_id = 9998
  []
  # Stitching assemblies together to form the core mesh
  [core]
    type = PatternedHexMeshGenerator
    inputs = 'Patterned cd_0 cd_1 cd_2 cd_3 cd_4 cd_5 cd_6 cd_7 cd_8 cd_9 cd_10 cd_11 ref_0 ref_1 ref_2 ref_3 ref_4 ref_5 dummy air_center'
    # Pattern ID  0       1    2    3    4    5    6    7    8    9    10    11    12    13    14    15    16    17    18   19       20
    pattern_boundary = none
    generate_core_metadata = true
    pattern = '19 17 12 16 19;
             17 11  0  0  1 16;
           10  0  0  0  0  0  2;
          18 0  0   0  0  0  0 15;
        19  9  0  0  20  0  0  3 19;
          18 0  0   0  0  0  0 15;
            8  0  0  0  0  0  4;
             13  7  0  0  5 14;
               19 13  6 14 19'
    rotate_angle = 60
    assign_control_drum_id = true
  []
  # Delete the dummy assemblies
  [del_dummy]
    type = BlockDeletionGenerator
    block = '700 701'
    input = core
    new_boundary = 10000
  []
  # Add peripheral circular region
  [add_outer_shield]
    type = PeripheralRingMeshGenerator
    input = del_dummy
    peripheral_layer_num = 1
    peripheral_ring_radius = 115.0
    input_mesh_external_boundary = 10000
    peripheral_ring_block_id = 250
    peripheral_ring_block_name = outer_shield
 []
 # trim the full core to 1/6 core
 [del_1]
    type = PlaneDeletionGenerator
    point = '0 0 0'
    normal = '10 17.32 0'
    input = add_outer_shield
    new_boundary = 147
  []
  [del_2]
    type = PlaneDeletionGenerator
    point = '0 0 0'
    normal = '10 -17.32 0'
    input = del_1
    new_boundary = 147
  []
  # extrusion
  [extrude]
     type = AdvancedExtruderGenerator
     input = del_2
     heights = '20 160 20'
     # Use `num_layers = '6 16 6'` for BISON mesh
     num_layers = '1 8 1'
     # biased upper and lower reflector mesh, uncomment for BISON mesh
     # biases = '1.6 1.0 0.625'
     direction = '0 0 1'
     top_boundary = 2000
     bottom_boundary = 3000
  []
  # Define some special reflector blocks
  [reflector_bottom_quad]
    type = ParsedSubdomainMeshGenerator
    input = extrude
    combinatorial_geometry = 'z<=20'
    block_id = 1000
    excluded_subdomains = '103 203 303 250 400 401 500 501 502 503 504 600 601'
  []
  [reflector_bottom_tri]
    type = ParsedSubdomainMeshGenerator
    input = reflector_bottom_quad
    combinatorial_geometry = 'z<=20'
    block_id = 1003
    excluded_subdomains = '1000 250 400 401 500 501 502 503 504 600 601'
  []
  [reflector_top_quad]
    type = ParsedSubdomainMeshGenerator
    input = reflector_bottom_tri
    combinatorial_geometry = 'z>=180'
    block_id = 1000
    excluded_subdomains = '103 203 303 200 201 250 400 401 500 501 502 503 504 600 601'
  []
  [reflector_top_tri]
    type = ParsedSubdomainMeshGenerator
    input = reflector_top_quad
    combinatorial_geometry = 'z>=180'
    block_id = 1003
    excluded_subdomains = '1000 200 201 203 250 400 401 500 501 502 503 504 600 601'
  []
  # Assgin block names
  [rename_blocks]
    type = RenameBlockGenerator
    old_block = '     101       100         103          201        200           203            301       303        400            401         10      503      600         601        504          500                501          502         1000           1003'
    new_block = ' mod_ss moderator_quad moderator_tri hp_ss  heat_pipes_quad heat_pipes_tri fuel_quad fuel_tri reflector_tri reflector_quad monolith B4C air_gap_tri air_gap_quad air_gap_quad reflector_tri  reflector_quad reflector_quad reflector_quad reflector_tri'
    input = reflector_top_tri
  []
  # Assign boundary names
  [rename_boundaries]
    type = RenameBoundaryGenerator
    input = rename_blocks
    old_boundary = ' 10000 2000 3000 147'
    new_boundary = ' side top bottom side_mirror'
  []
  # Scale to m unit for general MOOSE application
  # Note that is_meter needs to be true in Griffin
  [scale]
    type = TransformGenerator
    input = rename_boundaries
    transform = SCALE
    vector_value = '1e-2 1e-2 1e-2'
  []
  # Rotate to the desired orientation for simulation
  [rotate_end]
    type = TransformGenerator
    input = scale
    transform = ROTATE
    vector_value = '0 0 -120'
  []
[]

(microreactors/mrad/mesh/HPMR_OneSixth_Core_meshgenerator_tri_fine.i)

#####################################################################
# This is the MOOSE input file to generate the mesh that can be
# used for the 1/6 core Heat-Pipe Microreactor Multiphysics
# simulations (BISON thermal simulation).
# Running this input requires MOOSE Reactor Module Objects
# Users should use
# --mesh-only HPMR_OneSixth_Core_meshgenerator_tri_rotate_bdry_fine.e
# command line argument to generate
# exodus file for further Multiphysics simulations.
#####################################################################
[GlobalParams]
  # we use ':=' to override the values already defined in the !include file
  create_outward_interface_boundaries := false
[]

[Mesh]
  [HP_hex]
    # inner background boundary layer, only for BISON mesh
    background_inner_boundary_layer_bias = 1.5
    background_inner_boundary_layer_intervals = 3
    background_inner_boundary_layer_width = 0.03
  []
  [add_outer_shield]
    peripheral_layer_num := 2
    # Peripheral boundary layer, only for BISON mesh
    peripheral_outer_boundary_layer_bias = 0.625
    peripheral_outer_boundary_layer_intervals = 3
    peripheral_outer_boundary_layer_width = 2
  []
  [extrude]
    num_layers := '6 16 6'
    # biased upper and lower reflector mesh, only for BISON mesh
    biases = '1.6 1.0 0.625'
  []
  [split_hp_ss]
    type = SubdomainBoundingBoxGenerator
    input = rotate_end
    block_id = 1201
    block_name = 'hp_ss_up'
    restricted_subdomains = 'hp_ss'
    bottom_left = '-100 -100 1.8'
    top_right = '100 100 3.0'
  []
  [add_exterior_ht_low]
    type = SideSetsBetweenSubdomainsGenerator
    input = split_hp_ss
    paired_block = 'hp_ss'
    primary_block = 'monolith'
    new_boundary = 'heat_pipe_ht_surf_low'
  []
  [add_exterior_ht_up]
    type = SideSetsBetweenSubdomainsGenerator
    input = add_exterior_ht_low
    paired_block = 'hp_ss_up'
    primary_block = 'reflector_quad'
    new_boundary = 'heat_pipe_ht_surf_up'
  []
  [add_exterior_bot]
    type = SideSetsBetweenSubdomainsGenerator
    input = add_exterior_ht_up
    paired_block = 'heat_pipes_quad heat_pipes_tri hp_ss'
    primary_block = 'reflector_quad reflector_tri'
    new_boundary = 'heat_pipe_ht_surf_bot'
  []
  [merge_hp_surf]
    type = RenameBoundaryGenerator
    input = add_exterior_bot
    old_boundary = 'heat_pipe_ht_surf_low heat_pipe_ht_surf_up'
    new_boundary = 'heat_pipe_ht_surf heat_pipe_ht_surf'
  []
  [remove_hp]
    type = BlockDeletionGenerator
    input = merge_hp_surf
    block = 'hp_ss hp_ss_up heat_pipes_quad heat_pipes_tri'
  []

  # remove extra nodesets to limit the size of the mesh
  [clean_up]
    type = BoundaryDeletionGenerator
    input = remove_hp
    boundary_names = '1 3'
  []
[]

# This input defines the mesh generation for the Griffin mesh
!include HPMR_OneSixth_Core_meshgenerator_tri.i

(microreactors/mrad/steady/HPMR_dfem_griffin_ss.i)

################################################################################
## NEAMS Micro-Reactor Application Driver                                     ##
## Heat Pipe Microreactor Steady State                                        ##
## Griffin Main Application input file                                        ##
## DFEM-SN (1, 3) with CMFD acceleration                                      ##
################################################################################

[Mesh]
  [fmg]
    type = FileMeshGenerator
    file = '../mesh/gold/HPMR_OneSixth_Core_meshgenerator_tri_rotate_bdry.e'
  []
  [fmg_id]
    type = SubdomainExtraElementIDGenerator
    input = fmg
    subdomains = '200 203 100 103 301 303 10 503 600 601 201 101 400 401 250'
    extra_element_id_names = 'material_id equivalence_id'
    extra_element_ids = '815 815 802 802 801 801 803 811 820 820 817 816 805 805 820;
                         815 815 802 802 801 801 803 811 820 820 817 816 805 805 820'
  []
  [coarse_mesh]
    type = GeneratedMeshGenerator
    dim = 3
    nx = 10
    ny = 10
    nz = 10
    xmin = -0.1
    xmax = 1.1
    ymin = -0.1
    ymax = 1.2
    zmin = -0.0
    zmax = 2.1
  []
  [assign_coarse_id]
    type = CoarseMeshExtraElementIDGenerator
    input = fmg_id
    coarse_mesh = coarse_mesh
    extra_element_id_name = coarse_element_id
  []
  uniform_refine = 0
[]

[Executioner]
  type = SweepUpdate

  richardson_abs_tol = 1e-6 #-8
  richardson_rel_tol = 1e-8
  richardson_value = eigenvalue

  richardson_max_its = 1000
  inner_solve_type = GMRes
  max_inner_its = 20

  fixed_point_max_its = 1
  custom_pp = eigenvalue
  custom_rel_tol = 1e-6
  force_fixed_point_solve = true

  cmfd_acceleration = true
  coarse_element_id = coarse_element_id
  prolongation_type = multiplicative
  max_diffusion_coefficient = 1
[]

[TransportSystems]
  particle = neutron
  equation_type = eigenvalue

  G = 11
  VacuumBoundary = '10000 2000 3000'
  ReflectingBoundary = '147'

  [sn]
    scheme = DFEM-SN
    family = MONOMIAL
    order = first
    AQtype = Gauss-Chebyshev
    NPolar = 1
    NAzmthl = 3
    NA = 2
    sweep_type = asynchronous_parallel_sweeper
    using_array_variable = true
    collapse_scattering = true
    n_delay_groups = 6
  []
[]

[AuxVariables]
  [Tf]
    initial_condition = 873.15
    order = CONSTANT
    family = MONOMIAL
  []
  [Tm]
    initial_condition = 873.15
    order = CONSTANT
    family = MONOMIAL
  []
[]

[GlobalParams]
  library_file = '../isoxml/fullcore_xml_G11_endfb8_ss_tr.xml'
  library_name = fullcore_xml_G11_endfb8_ss_tr
  isotopes = 'pseudo'
  densities = 1.0
  is_meter = true
  plus = true
  dbgmat = false
  grid_names = 'Tfuel Tmod'
  grid_variables = 'Tf Tm'
[]

[Materials]
  [mod]
    type = CoupledFeedbackMatIDNeutronicsMaterial
    block = '200 203 100 103 301 303 10 503 600 601 201 101 400 401 250'
  []
[]

[PowerDensity]
  power = 345.6e3
  power_density_variable = power_density
  integrated_power_postprocessor = integrated_power
[]

[MultiApps]
  [bison]
    type = FullSolveMultiApp
    app_type = BisonApp
    positions = '0 0 0'
    input_files = HPMR_thermo_ss.i
    execute_on = 'timestep_end'
    keep_solution_during_restore = true
  []
[]

[Transfers]
  [to_sub_power_density]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    to_multi_app = bison
    source_variable = power_density
    variable = power_density
    from_postprocessors_to_be_preserved = integrated_power
    to_postprocessors_to_be_preserved = power

    search_value_conflicts = false
  []
  [from_sub_temp_fuel]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = bison
    variable = Tf
    source_variable = Tfuel
    execute_on = 'initial timestep_end'
    displaced_source_mesh = false
    displaced_target_mesh = false
    use_displaced_mesh = false

    search_value_conflicts = false
  []
  [from_sub_temp_mod]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = bison
    variable = Tm
    source_variable = Tmod
    execute_on = 'initial timestep_end'
    displaced_source_mesh = false
    displaced_target_mesh = false
    use_displaced_mesh = false

    search_value_conflicts = false
  []
[]

[UserObjects]
  [ss]
    type = TransportSolutionVectorFile
    transport_system = sn
    writing = true
    execute_on = final
  []
[]

[Postprocessors]
  [scaled_power_avg]
    type = ElementAverageValue
    block = 'fuel_quad fuel_tri'
    variable = power_density
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_avg]
    type = ElementAverageValue
    variable = Tf
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tf
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_min]
    type = ElementExtremeValue
    value_type = min
    variable = Tf
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [mod_temp_avg]
    type = ElementAverageValue
    variable = Tm
    block = 'moderator_quad moderator_tri'
    execute_on = 'initial timestep_end'
  []
  [mod_temp_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tm
    block = 'moderator_quad moderator_tri'
    execute_on = 'initial timestep_end'
  []
  [mod_temp_min]
    type = ElementExtremeValue
    value_type = min
    variable = Tm
    block = 'moderator_quad moderator_tri'
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  csv = true
  exodus = true
  perf_graph = true
[]

(microreactors/mrad/transient_null/HPMR_dfem_griffin_trN.i)

################################################################################
## NEAMS Micro-Reactor Application Driver                                     ##
## Heat Pipe Microreactor Null Transient                                      ##
## Griffin Main Application input file                                        ##
## DFEM-SN (1, 3) with CMFD acceleration                                      ##
################################################################################

[Mesh]
  [fmg]
    type = FileMeshGenerator
    file = '../mesh/gold/HPMR_OneSixth_Core_meshgenerator_tri_rotate_bdry.e'
  []
  [fmg_id]
    type = SubdomainExtraElementIDGenerator
    input = fmg
    subdomains = '200 203 100 103 301 303 10 503 600 601 201 101 400 401 250'
    extra_element_id_names = 'material_id equivalence_id'
    extra_element_ids = '815 815 802 802 801 801 803 811 820 820 817 816 805 805 820;
                         815 815 802 802 801 801 803 811 820 820 817 816 805 805 820'
  []
  [coarse_mesh]
    type = GeneratedMeshGenerator
    dim= 3
    nx = 10
    ny = 10
    nz = 10
    xmin = -0.1
    xmax =  1.1
    ymin = -0.1
    ymax =  1.2
    zmin = -0.0
    zmax =  2.1
  []
  [assign_coarse_id]
    type = CoarseMeshExtraElementIDGenerator
    input = fmg_id
    coarse_mesh = coarse_mesh
    extra_element_id_name = coarse_element_id
  []
  uniform_refine = 0
[]

[Executioner]
  type = TransientSweepUpdate

  richardson_abs_tol = 1e-6
  richardson_rel_tol = 1e-8
  richardson_max_its = 1000
  richardson_value = eigenvalue

  inner_solve_type = GMRes
  max_inner_its = 20

  fixed_point_max_its = 1
  # force_fixed_point_solve = true

  cmfd_acceleration = true
  coarse_element_id = coarse_element_id
  prolongation_type = multiplicative
  max_diffusion_coefficient = 1

  end_time = 200
  dt = 1
[]

[TransportSystems]
  particle = neutron
  equation_type = transient

  G = 11
  VacuumBoundary = '10000 2000 3000'
  ReflectingBoundary = '147'

  [sn]
    scheme = DFEM-SN
    family = MONOMIAL
    order = first
    AQtype = Gauss-Chebyshev
    NPolar = 1
    NAzmthl = 3
    NA = 2
    sweep_type = asynchronous_parallel_sweeper
    using_array_variable = true
    collapse_scattering  = true
    n_delay_groups = 6
  []
[]

[AuxVariables]
  [Tf]
    order = CONSTANT
    family = MONOMIAL
  []
  [Tm]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[GlobalParams]
  library_file = '../isoxml/fullcore_xml_G11_endfb8_ss_tr.xml'
  library_name = fullcore_xml_G11_endfb8_ss_tr
  isotopes = 'pseudo'
  densities = 1.0
  is_meter = true
  plus = true
  dbgmat = false
  grid_names = 'Tfuel Tmod'
  grid_variables = 'Tf Tm'
[]

[Materials]
  [mod]
    type = CoupledFeedbackMatIDNeutronicsMaterial
    block = '200 203 100 103 301 303 10 503 600 601 201 101 400 401 250'
  []
[]

[PowerDensity]
  power = 345.6e3
  power_density_variable = power_density
  integrated_power_postprocessor = integrated_power
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    positions = '0 0 0'
    input_files = HPMR_thermo_trN.i
    execute_on = 'timestep_end'
  []
[]


[Transfers]
  [to_sub_power_density]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    to_multi_app = bison
    source_variable = power_density
    variable = power_density
    from_postprocessors_to_be_preserved = integrated_power
    to_postprocessors_to_be_preserved = power

    search_value_conflicts = false
  []
  [from_sub_temp_fuel]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = bison
    variable = Tf
    source_variable = Tfuel
    execute_on = 'initial timestep_end'
    displaced_source_mesh = false
    displaced_target_mesh = false
    use_displaced_mesh = false

    search_value_conflicts = false
  []
  [from_sub_temp_mod]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = bison
    variable = Tm
    source_variable = Tmod
    execute_on = 'initial timestep_end'
    displaced_source_mesh = false
    displaced_target_mesh = false
    use_displaced_mesh = false

    search_value_conflicts = false
  []
[]

[Postprocessors]
  [scaled_power_avg]
    type = ElementAverageValue
    block = 'fuel_quad fuel_tri'
    variable = power_density
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_avg]
    type = ElementAverageValue
    variable = Tf
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tf
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_min]
    type = ElementExtremeValue
    value_type = min
    variable = Tf
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [mod_temp_avg]
    type = ElementAverageValue
    variable = Tm
    block = 'moderator_quad moderator_tri'
    execute_on = 'initial timestep_end'
  []
  [mod_temp_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tm
    block = 'moderator_quad moderator_tri'
    execute_on = 'initial timestep_end'
  []
  [mod_temp_min]
    type = ElementExtremeValue
    value_type = min
    variable = Tm
    block = 'moderator_quad moderator_tri'
    execute_on = 'initial timestep_end'
  []
[]

[UserObjects]
  [ss]
    type = TransportSolutionVectorFile
    transport_system = sn
    writing = false
    execute_on = initial
  []
[]

[Outputs]
  csv = true
  exodus = false
  perf_graph=true
  checkpoint = true
[]

(microreactors/mrad/heat_pipe_failure/HPMR_dfem_griffin_tr.i)

################################################################################
## NEAMS Micro-Reactor Application Driver                                     ##
## Heat Pipe Microreactor Heat Pipe Failure Transient                         ##
## Griffin Main Application input file                                        ##
## DFEM-SN (1, 3) with CMFD acceleration                                      ##
################################################################################

[Mesh]
  [fmg]
    type = FileMeshGenerator
    file = '../mesh/gold/HPMR_OneSixth_Core_meshgenerator_tri_rotate_bdry.e'
  []
  [fmg_id]
    type = SubdomainExtraElementIDGenerator
    input = fmg
    subdomains = '200 203 100 103 301 303 10 503 600 601 201 101 400 401 250'
    extra_element_id_names = 'material_id equivalence_id'
    extra_element_ids = '815 815 802 802 801 801 803 811 820 820 817 816 805 805 820;
                         815 815 802 802 801 801 803 811 820 820 817 816 805 805 820'
  []
  [coarse_mesh]
    type = GeneratedMeshGenerator
    dim= 3
    nx = 10
    ny = 10
    nz = 10
    xmin = -0.1
    xmax =  1.1
    ymin = -0.1
    ymax =  1.2
    zmin = -0.0
    zmax =  2.1
  []
  [assign_coarse_id]
    type = CoarseMeshExtraElementIDGenerator
    input = fmg_id
    coarse_mesh = coarse_mesh
    extra_element_id_name = coarse_element_id
  []
  uniform_refine = 0
[]

[Executioner]
  type = TransientSweepUpdate

  richardson_abs_tol = 1e-6
  richardson_rel_tol = 1e-8
  richardson_max_its = 1000
  richardson_value = eigenvalue

  inner_solve_type = GMRes
  max_inner_its = 20

  fixed_point_max_its = 1
  # force_fixed_point_solve = true

  cmfd_acceleration = true
  coarse_element_id = coarse_element_id
  prolongation_type = multiplicative
  max_diffusion_coefficient = 1

  end_time = 2000
  dt = 1
[]

[TransportSystems]
  particle = neutron
  equation_type = transient

  G = 11
  VacuumBoundary = '10000 2000 3000'
  ReflectingBoundary = '147'

  [sn]
    scheme = DFEM-SN
    family = MONOMIAL
    order = first
    AQtype = Gauss-Chebyshev
    NPolar = 1
    NAzmthl = 3
    NA = 2
    sweep_type = asynchronous_parallel_sweeper
    using_array_variable = true
    collapse_scattering  = true
    n_delay_groups = 6
  []
[]

[AuxVariables]
  [Tf]
    order = CONSTANT
    family = MONOMIAL
  []
  [Tm]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[GlobalParams]
  library_file = '../isoxml/fullcore_xml_G11_endfb8_ss_tr.xml'
  library_name = fullcore_xml_G11_endfb8_ss_tr
  isotopes = 'pseudo'
  densities = 1.0
  is_meter = true
  plus = true
  dbgmat = false
  grid_names = 'Tfuel Tmod'
  grid_variables = 'Tf Tm'
[]

[Materials]
  [mod]
    type = CoupledFeedbackMatIDNeutronicsMaterial
    block = '200 203 100 103 301 303 10 503 600 601 201 101 400 401 250'
  []
[]

[PowerDensity]
  power = 345.6e3
  power_density_variable = power_density
  integrated_power_postprocessor = integrated_power
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    positions = '0 0 0'
    input_files = HPMR_thermo_tr.i
    execute_on = 'timestep_end'
  []
[]


[Transfers]
  [to_sub_power_density]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    to_multi_app = bison
    source_variable = power_density
    variable = power_density
    from_postprocessors_to_be_preserved = integrated_power
    to_postprocessors_to_be_preserved = power

    search_value_conflicts = false
  []
  [from_sub_temp_fuel]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = bison
    variable = Tf
    source_variable = Tfuel
    execute_on = 'initial timestep_end'
    displaced_source_mesh = false
    displaced_target_mesh = false
    use_displaced_mesh = false

    search_value_conflicts = false
  []
  [from_sub_temp_mod]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = bison
    variable = Tm
    source_variable = Tmod
    execute_on = 'initial timestep_end'
    displaced_source_mesh = false
    displaced_target_mesh = false
    use_displaced_mesh = false

    search_value_conflicts = false
  []
[]

[Postprocessors]
  [scaled_power_avg]
    type = ElementAverageValue
    block = 'fuel_quad fuel_tri'
    variable = power_density
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_avg]
    type = ElementAverageValue
    variable = Tf
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tf
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_min]
    type = ElementExtremeValue
    value_type = min
    variable = Tf
    block = 'fuel_quad fuel_tri'
    execute_on = 'initial timestep_end'
  []
  [mod_temp_avg]
    type = ElementAverageValue
    variable = Tm
    block = 'moderator_quad moderator_tri'
    execute_on = 'initial timestep_end'
  []
  [mod_temp_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tm
    block = 'moderator_quad moderator_tri'
    execute_on = 'initial timestep_end'
  []
  [mod_temp_min]
    type = ElementExtremeValue
    value_type = min
    variable = Tm
    block = 'moderator_quad moderator_tri'
    execute_on = 'initial timestep_end'
  []
[]

[UserObjects]
  [ss]
    type = TransportSolutionVectorFile
    transport_system = sn
    writing = false
    execute_on = initial
  []
[]

[Outputs]
  csv = true
  exodus = false
  perf_graph=true
  checkpoint = true
[]

Mesh File
Multiphysics Model Setup
Simulation Cases
Guidance of Running Different Simulation Cases
References