Multiphysics Coupling of OpenMC, MOOSE, and THM for a prismatic HTGR

Contact: April Novak, anovak.at.anl.gov

In this tutorial, we couple OpenMC Monte Carlo transport to the MOOSE heat transfer module and the MOOSE Thermal-Hydraulics Module (THM), a set of 1-D systems-level thermal-hydraulics kernels in MOOSE Berry et al. (2014), for application to a generic Tri-Structural Isotropic (TRISO)-fueled High Temperature Gas Reactor (HTGR) assembly. This coupling will be performed with Cardinal, an open-source application available on GitHub. OpenMC will receive temperature feedback from both the MOOSE heat conduction module (for the solid regions) and from THM (for the fluid regions). Density feedback will be provided by THM for the fluid regions.

This tutorial was developed with support from the NEAMS Thermal Fluids Center of Excellence. A paper Novak et al. (2021) describing the physics models and mesh refinement studies provides additional context beyond the scope of this tutorial.

Geometry and Computational Model

The geometry consists of a TRISO-fueled gas reactor assembly, loosely based on a point design available in the literature INL (2016). A top-down view of the geometry is shown in Figure 1. The assembly is a graphite prismatic hexagonal block with 108 helium coolant channels, 210 fuel compacts, and 6 poison compacts. Each fuel compact contains TRISO particles dispersed in a graphite matrix with a packing fraction of 15%. All channels and compacts are ordered in a triangular lattice with pitch $var element = document.getElementById("moose-equation-8d69a465-cf2c-46fb-918d-66a919ed9dcc");katex.render("p_{cf}", element, {displayMode:false,throwOnError:false});$ . Due to irradiation- and temperature-induced swelling of graphite, small helium gaps exist between assemblies. In this tutorial, rather than explicitly model the inter-assembly flow, we treat the gap regions as solid graphite. There are also graphite reflectors above and below the assembly.

Figure 1: TRISO-fueled gas reactor fuel assembly

The TRISO particles use a conventional design that consists of a central fissile uranium oxycarbide kernel enclosed in a carbon buffer, an inner Pyrolitic Carbon (PyC) layer, a silicon carbide layer, and finally an outer PyC layer. The geometric specifications for the assembly dimensions are shown in Figure 1, while dimensions for the TRISO particles are summarized in Table 1.

Table 1: Geometric specifications for the TRISO particles

Parameter	Value (cm)
TRISO kernel radius	214.85e-4
Buffer layer radius	314.85e-4
Inner PyC layer radius	354.85e-4
Silicon carbide layer radius	389.85e-4
Outer PyC layer radius	429.85e-4

Heat is produced in the TRISO particles to yield a total power of 16.7 MWth. This heat is removed by helium flowing downwards through the coolant channels with a total mass flowrate of 9.775 kg/s, which is assumed to be uniformly distributed among the coolant channels. The outlet pressure is 7.1 MPa.

Heat Conduction Model

The MOOSE heat conduction module is used to solve for steady-state heat conduction,

$(1)var element = document.getElementById("moose-equation-1a69e08f-8ae5-4c9f-9e52-5ad71110ea9b");katex.render(" -\\nabla\\cdot\\left(k_s\\nabla T_s\\right)=\\dot{q}_s", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-4ccfdda6-9261-4849-86fa-522e4e85a05a");katex.render("k_s", element, {displayMode:false,throwOnError:false});$ is the solid thermal conductivity, $var element = document.getElementById("moose-equation-9a184fd5-0628-46b8-b1f9-709d9cc56951");katex.render("T_s", element, {displayMode:false,throwOnError:false});$ is the solid temperature, and $var element = document.getElementById("moose-equation-59f4df74-62d4-446a-a1f4-9c548a74cda8");katex.render("\\dot{q}_s", element, {displayMode:false,throwOnError:false});$ is a volumetric heat source in the solid.

To greatly reduce meshing requirements, the TRISO particles are homogenized into the compact regions by volume-averaging material properties. The solid mesh is shown in Figure 2. The only sideset in the domain is the coolant channel surface, which is named fluid_solid_interface. To simplify the specification of material properties, the solid geometry uses a length unit of meters. The solid mesh is created using the MOOSE reactor module Shemon et al. (2022), which provides easy-to-use mesh generators to programmatically construct reactor core meshes as building blocks of bundle and pincell meshes.

Figure 2: Mesh for the solid heat conduction model

The file used to generate the solid mesh is shown below. The mesh is created by first building pincell meshes for a fuel pin, a coolant pin, a poison pin, and a graphite "pin" (to represent the central graphite region). The pin meshes are then combined together into a bundle pattern and extruded into the $var element = document.getElementById("moose-equation-403be384-942d-4c15-b090-182f6354c9d0");katex.render("z", element, {displayMode:false,throwOnError:false});$ direction.

# HTGR assembly meshing file
# Application: MOOSE reactor module, Griffin, BlueCRAB
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html

# This input file is used to create the mesh for the solid phase.
# This should be run with:
#
# cardinal-opt -i common_input.i solid.i --mesh-only

n_layers = 100

[GlobalParams]
  quad_center_elements = true
[]

[Mesh]
  [fuel_pin]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    polygon_size = '${fparse fuel_to_coolant_distance / 2.0}'
    ring_radii = '${fparse 0.8 * compact_diameter / 2.0} ${fparse compact_diameter / 2.0}'
    ring_intervals = '1 1'
    num_sectors_per_side = '4 4 4 4 4 4'
    ring_block_ids = '2 2'
    ring_block_names = 'compacts compacts'
    background_block_ids = '1'
    background_block_names = 'graphite'
    background_intervals = 2
  []
  [coolant_pin]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    polygon_size = '${fparse fuel_to_coolant_distance / 2.0}'
    ring_radii = '${fparse channel_diameter / 2.0}'
    ring_intervals = '2'
    num_sectors_per_side = '4 4 4 4 4 4'
    ring_block_ids = '101 101'
    ring_block_names = 'coolant coolant'
    background_block_ids = '1'
    background_block_names = 'graphite'
    interface_boundary_id_shift = 100
    background_intervals = 2
    create_inward_interface_boundaries = true
  []
  [poison_pin]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    polygon_size = '${fparse fuel_to_coolant_distance / 2.0}'
    ring_radii = '${fparse compact_diameter / 2.0}'
    ring_intervals = '2'
    num_sectors_per_side = '4 4 4 4 4 4'
    ring_block_ids = '4 4'
    ring_block_names = 'poison poison'
    background_block_ids = '1'
    background_block_names = 'graphite'
    background_intervals = 2
  []
  [graphite_pin]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    polygon_size = '${fparse fuel_to_coolant_distance / 2.0}'
    ring_radii = '${fparse compact_diameter / 2.0}'
    ring_intervals = '2'
    num_sectors_per_side = '4 4 4 4 4 4'
    ring_block_ids = '1 1'
    ring_block_names = 'graphite graphite'
    background_block_ids = '1'
    background_block_names = 'graphite'
  []
  [bundle]
    type = PatternedHexMeshGenerator
    inputs = 'fuel_pin coolant_pin poison_pin graphite_pin'
    hexagon_size = '${fparse bundle_flat_to_flat / 2.0 + bundle_gap_width / 2.0}'
    pattern = '2 0 1 0 0 1 0 0 1 0 2;
              0 1 0 0 1 0 0 1 0 0 1 0;
             1 0 0 1 0 0 1 0 0 1 0 0 1;
            0 0 1 0 0 1 0 0 1 0 0 1 0 0;
           0 1 0 0 1 0 0 1 0 0 1 0 0 1 0;
          1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1;
         0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0;
        0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0;
       1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1;
      0 0 1 0 0 1 0 0 1 3 3 1 0 0 1 0 0 1 0 0;
     2 1 0 0 1 0 0 1 0 3 3 3 0 1 0 0 1 0 0 1 2;
      0 0 1 0 0 1 0 0 1 3 3 1 0 0 1 0 0 1 0 0;
       1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1;
        0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0;
         0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0;
          1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1;
           0 1 0 0 1 0 0 1 0 0 1 0 0 1 0;
            0 0 1 0 0 1 0 0 1 0 0 1 0 0;
             1 0 0 1 0 0 1 0 0 1 0 0 1;
              0 1 0 0 1 0 0 1 0 0 1 0;
               2 0 1 0 0 1 0 0 1 0 2'
    rotate_angle = 0

    background_block_id = '1'
    background_block_names = 'graphite'
  []
  [extrude]
    type = AdvancedExtruderGenerator
    input = bundle
    heights = ${height}
    num_layers = ${n_layers}
    direction = '0 0 1'
  []
  [delete_coolant]
    type = BlockDeletionGenerator
    input = extrude
    block = '101'
  []
  [rename_coolant_sideset]
    type = RenameBoundaryGenerator
    input = delete_coolant
    old_boundary = 102
    new_boundary = 'fluid_solid_interface'
  []

  construct_side_list_from_node_list = true
[]

You can create this mesh by running:


moose-opt -i common_input.i solid_mesh.i --mesh-only

which will create a mesh named solid_mesh_in.e. Note that the above command takes advantage of a MOOSE feature for combining input files together by placing some common parameters used by the other applications into a file named common_input.i.

The temperature on the fluid-solid interface is provided by THM, while the heat source is provided by OpenMC. Because MOOSE heat conduction will run first in the coupled case, the initial fluid temperature is set to an axial distribution given by bulk energy conservation ( $var element = document.getElementById("moose-equation-edfc0bf8-b67b-4147-87cf-d170805bd60f");katex.render("q=\\dot{m}C_{p,f}\\left(T_f-T_{inlet}\\right)", element, {displayMode:false,throwOnError:false});$ ) given the inlet temperature $var element = document.getElementById("moose-equation-91ef32df-cfd5-49a5-b3e4-a6439185f5f4");katex.render("T_{inlet}", element, {displayMode:false,throwOnError:false});$ , mass flowrate $var element = document.getElementById("moose-equation-22e00e1b-66db-473a-8b0a-55ff02f8da6e");katex.render("\\dot{m}", element, {displayMode:false,throwOnError:false});$ , and fluid isobaric specific heat $var element = document.getElementById("moose-equation-2b1a7cfb-a7a6-41c6-86a5-ce22d72b1b37");katex.render("C_{p,f}", element, {displayMode:false,throwOnError:false});$ . The initial heat source distribution is assumed uniform in the radial direction with a sinusoidal dependence in the axial direction.

OpenMC Model

The OpenMC model is built using Constructive Solid Geometry (CSG). The TRISO positions are sampled using the Random Sequential Addition (RSA) algorithm in OpenMC. OpenMC's Python Application Programming Interface (API) is used to create the model with the script shown below. First, we define materials for the various regions. Next, we create a single TRISO particle universe consisting of the five layers of the particle and an infinite extent of graphite filling all other space. We then pack uniform-radius spheres into a cylindrical region representing a fuel compact, setting each sphere to be filled with the TRISO universe.

Next, we loop over 50 axial layers and create a unique hexagonal lattice for each layer. This hexagonal lattice defines the fuel assembly structure, and consists of four different universes:

A fuel pin plus surrounding matrix (f),
A coolant channel plus surrounding matrix (c),
A boron carbide poision pin plus surrounding matrix (p), and
A homogeneous graphite hexagonal pincell to fill the "boundaries" and centermost region (g).

In each layer we set up the lattice structure by listing the universes in each "ring" of the lattice, with ring0 being the centermost ring and ring11 being the outermost ring.

Temperatures in OpenMC can be set directly on the cell, but fluid densities can only be set on materials. For this reason, we need to create 108 unique coolant materials for each axial plane if we want to be able to set unique densities in each coolant channel region. Rather than creating 108 materials in a loop, we use the clone() feature in OpenMC to clone an existing coolant material 108 times per layer. This duplicates the material properties (densities and isotopic composition), but assigns a new ID that allows individual tracking of density. The Python script used to create the OpenMC model is shown below.

#!/bin/env python

from argparse import ArgumentParser
import math
import numpy as np
import matplotlib.pyplot as plt

import openmc
import sys
import os

# Get common input parameters shared by other physics
script_dir = os.path.dirname(__file__)
sys.path.append(script_dir)
import common_input as specs
import materials as mats

def coolant_temp(t_in, t_out, l, z):
    """
    Computes the coolant temperature based on an expected cosine power distribution
    for a specified temperature rise. The total core temperature rise is governed
    by energy conservation as dT = Q / m / Cp, where dT is the total core temperature
    rise, Q is the total core power, m is the mass flowrate, and Cp is the fluid
    isobaric specific heat. If you neglect axial heat conduction and assume steady
    state, then the temperature rise in a layer of fluid i can be related to the
    ratio of the power in that layer to the total power,
    dT_i / dT = Q_i / Q. We assume here a sinusoidal power distribution to get
    a reasonable estimate of an initial coolant temperature distribution.

    Parameters
    ----------

    t_in : float
        Inlet temperature of the channel
    t_out : float
        Outlet temperature of the channel
    l : float
        Length of the channel
    z : float or 1-D numpy.array
        Axial position where the temperature will be computed

    Returns
    -------
        float or 1-D numpy array of float depending on z
    """
    dT = t_out - t_in
    Q = 2 * l / math.pi
    Qi = (l - l * np.cos(math.pi * z / l)) / math.pi

    t = t_in + Qi / Q * dT

    return t

def coolant_density(t):
  """
  Computes the helium density from temperature assuming a fixed operating pressure.

  Parameters
  ----------

  t : float
    Fluid temperature

  Returns
  _______
    float or 1-D numpy array of float depending on t
  """

  p_in_bar = specs.outlet_P * 1.0e-5;
  return 48.14 * p_in_bar / (t + 0.4446 * p_in_bar / math.pow(t, 0.2));

# -------------- Unit Conversions: OpenMC requires cm -----------
m = 100.0

# estimate the outlet temperature using bulk energy conservation for steady state
coolant_outlet_temp = specs.power / specs.mdot / specs.fluid_Cp + specs.inlet_T

# geometry parameters
coolant_channel_diam = specs.channel_diameter * m
reactor_bottom = 0.0
reactor_height = specs.height * m
reactor_top = reactor_bottom + reactor_height
bundle_pitch = specs.bundle_flat_to_flat * m + specs.bundle_gap_width * m
cell_pitch = specs.fuel_to_coolant_distance * m
fuel_channel_diam = specs.compact_diameter * m
top_reflector_height = specs.top_reflector_thickness * m
bottom_reflector_height = specs.bottom_reflector_thickness * m

def assembly(n_ax_zones, n_inactive, n_active, add_entropy_mesh=False):
    axial_section_height = reactor_height / n_ax_zones

    # superimposed search lattice
    triso_lattice_shape = (4, 4, int(axial_section_height / 0.5))

    model = openmc.model.Model()

    m_b4c = openmc.Material(name='B4C')
    enrichment_10 = specs.B10_enrichment
    mass_10 = openmc.data.atomic_mass('B10')
    mass_11 = openmc.data.atomic_mass('B11')

    # number of atoms in one gram of boron mixture
    n_10 = enrichment_10 / mass_10
    n_11 = (1.0 - enrichment_10) / mass_11
    total_n = n_10 + n_11
    grams_10 = n_10 / total_n
    grams_11 = n_11 / total_n

    # now, figure out how much carbon needs to be in the poison to get
    # an overall specified B10 weight percent
    total_b10_weight_percent = specs.total_B10_wt_percent
    total_mass = grams_10 / total_b10_weight_percent
    carbon_mass = total_mass - grams_10 - grams_11

    m_b4c.add_nuclide('B10', grams_10 / total_mass, 'wo')
    m_b4c.add_nuclide('B11', grams_11 / total_mass, 'wo')
    m_b4c.add_element('C', carbon_mass / total_mass, 'wo')
    m_b4c.set_density('kg/m3', specs.B4C_density)

    # reflector is 40 percent helium by volume (arbitrary assumption), with helium
    # at the inlet conditions
    rho = coolant_density(specs.inlet_T)
    matrix_density = 1700
    reflector_porosity = 0.40
    n_helium = reflector_porosity * rho / 4.002602
    n_carbon = (1.0 - reflector_porosity) * matrix_density / 12.0107
    combined_density = rho * reflector_porosity + matrix_density * (1.0 - reflector_porosity)
    m_reflector = openmc.Material(name='reflector')
    m_reflector.add_element('He', n_helium / (n_helium + n_carbon))
    m_reflector.add_element('C', n_carbon / (n_helium + n_carbon))
    m_reflector.set_density('kg/m3', combined_density)

    # TRISO particle
    radius_pyc_outer   = specs.oPyC_radius * m

    s_fuel             = openmc.Sphere(r=specs.kernel_radius*m)
    s_c_buffer         = openmc.Sphere(r=specs.buffer_radius*m)
    s_pyc_inner        = openmc.Sphere(r=specs.iPyC_radius*m)
    s_sic              = openmc.Sphere(r=specs.SiC_radius*m)
    s_pyc_outer        = openmc.Sphere(r=radius_pyc_outer)
    c_triso_fuel       = openmc.Cell(name='c_triso_fuel'     , fill=mats.m_fuel,              region=-s_fuel)
    c_triso_c_buffer   = openmc.Cell(name='c_triso_c_buffer' , fill=mats.m_graphite_c_buffer, region=+s_fuel      & -s_c_buffer)
    c_triso_pyc_inner  = openmc.Cell(name='c_triso_pyc_inner', fill=mats.m_graphite_pyc,      region=+s_c_buffer  & -s_pyc_inner)
    c_triso_sic        = openmc.Cell(name='c_triso_sic'      , fill=mats.m_sic,               region=+s_pyc_inner & -s_sic)
    c_triso_pyc_outer  = openmc.Cell(name='c_triso_pyc_outer', fill=mats.m_graphite_pyc,      region=+s_sic       & -s_pyc_outer)
    c_triso_matrix     = openmc.Cell(name='c_triso_matrix'   , fill=mats.m_graphite_matrix,   region=+s_pyc_outer)
    u_triso            = openmc.Universe(cells=[c_triso_fuel, c_triso_c_buffer, c_triso_pyc_inner, c_triso_sic, c_triso_pyc_outer, c_triso_matrix])

    # Channel surfaces
    fuel_cyl = openmc.ZCylinder(r=0.5 * fuel_channel_diam)
    coolant_cyl = openmc.ZCylinder(r=0.5 * coolant_channel_diam)
    poison_cyl = openmc.ZCylinder(r=0.5 * fuel_channel_diam)
    graphite_cyl = openmc.ZCylinder(r=0.5 * fuel_channel_diam)

    # create a TRISO lattice for one axial section (to be used in the rest of the axial zones)
    # center the TRISO region on the origin so it fills lattice cells appropriately
    min_z = openmc.ZPlane(z0=-0.5 * axial_section_height)
    max_z = openmc.ZPlane(z0=0.5 * axial_section_height)

    # region in which TRISOs are generated
    r_triso = -fuel_cyl & +min_z & -max_z

    rand_spheres = openmc.model.pack_spheres(radius=radius_pyc_outer, region=r_triso, pf=specs.triso_pf)
    random_trisos = [openmc.model.TRISO(radius_pyc_outer, u_triso, i) for i in rand_spheres]

    llc, urc = r_triso.bounding_box
    pitch = (urc - llc) / triso_lattice_shape

    # insert TRISOs into a lattice to accelerate point location queries
    triso_lattice = openmc.model.create_triso_lattice(random_trisos, llc, pitch, triso_lattice_shape, mats.m_graphite_matrix)

    axial_coords = np.linspace(reactor_bottom, reactor_top, n_ax_zones + 1)
    lattice_univs = []
    bundle_univs = []

    m_colors = {}

    for z_min, z_max in zip(axial_coords[0:-1], axial_coords[1:]):
        # use the middle of the axial section to compute the temperature and density
        ax_pos = 0.5 * (z_min + z_max)
        T = coolant_temp(specs.inlet_T, coolant_outlet_temp, reactor_height, ax_pos)

        # create solid cells, which don't require us to clone materials in order to set temperatures
        fuel_ch_cell = openmc.Cell(region=-fuel_cyl, fill=triso_lattice)
        fuel_ch_cell.temperature = T

        fuel_ch_matrix_cell = openmc.Cell(region=+fuel_cyl, fill=mats.m_graphite_matrix)
        fuel_ch_matrix_cell.temperature = T

        poison_cell = openmc.Cell(region=-poison_cyl, fill=m_b4c)
        poison_cell.temperature = T

        poison_matrix_cell = openmc.Cell(region=+poison_cyl, fill=mats.m_graphite_matrix)
        poison_matrix_cell.temperature = T

        graphite_cell = openmc.Cell(fill=mats.m_graphite_matrix)
        graphite_cell.temperature = T

        coolant_matrix_cell = openmc.Cell(region=+coolant_cyl, fill=mats.m_graphite_matrix)
        coolant_matrix_cell.temperature = T

        # create fluid cells, which require us to clone the material in order to be able to
        # set unique densities
        coolant_cell = openmc.Cell(region=-coolant_cyl, fill=mats.m_coolant)
        coolant_cell.fill = [mats.m_coolant.clone() for i in range(specs.n_coolant_channels_per_block)]

        for mat in range(len(coolant_cell.fill)):
          m_colors[coolant_cell.fill[mat]] = 'red'

        # Define a universe for each type of solid-only pin (fuel, poison, and graphite)
        f = openmc.Universe(cells=[fuel_ch_cell, fuel_ch_matrix_cell])
        p = openmc.Universe(cells=[poison_cell, poison_matrix_cell])
        c = openmc.Universe(cells=[coolant_cell, coolant_matrix_cell])
        g = openmc.Universe(cells=[graphite_cell])

        d = [f] * 2

        ring2 = ([f] + [c]) * 6
        ring3 = ([c] + d) * 6
        ring4 = (d + [c] + [f]) * 6
        ring5 = ([f] + [c] + d + [c]) * 6
        ring6 = ([c] + d + [c] + d) * 6
        ring7 = (d + [c] + d + [c] + [f]) * 6
        ring8 = ([f] + [c] + d + [c] + d + [c]) * 6
        ring9 = ([c] + d + [c] + d + [c] + d) * 6
        ring10 = ([p] + [f] + [c] + d + [c] + d + [c] + [f]) * 6
        ring11 = [g] * 66

        # inner two rings where there aren't any fuel/compact/poison pins
        ring1 = [g] * 6
        ring0 = [g]

        lattice_univs.append([ring11, ring10, ring9, ring8, ring7, ring6, ring5, ring4, ring3, ring2, ring1, ring0])

    # create a hexagonal lattice used in each axial zone to represent the cell
    hex_lattice = openmc.HexLattice(name="Bundle cell lattice")
    hex_lattice.orientation = 'x'
    hex_lattice.center = (0.0, 0.0, 0.5 * (reactor_bottom + reactor_top))
    hex_lattice.pitch = (cell_pitch, axial_section_height)
    hex_lattice.universes = lattice_univs

    hexagon_volume = reactor_height * math.sqrt(3) / 2.0 * bundle_pitch**2
    coolant_channel_volume = math.pi * coolant_channel_diam**2 / 4.0 * reactor_height

    graphite_outer_cell = openmc.Cell(fill=mats.m_graphite_matrix)
    graphite_outer_cell.temperature = T
    inf_graphite_univ = openmc.Universe(cells=[graphite_outer_cell])
    hex_lattice.outer = inf_graphite_univ

    # create additional axial regions
    axial_planes = [openmc.ZPlane(z0=coord) for coord in axial_coords]

    # axial planes
    min_z = axial_planes[0]
    max_z = axial_planes[-1]

    # fill the unit cell with the hex lattice
    hex_prism = openmc.hexagonal_prism(bundle_pitch / math.sqrt(3.0), 'x', boundary_type='periodic')
    outer_cell = openmc.Cell(region=hex_prism & +min_z & -max_z, fill=hex_lattice)

    # add the top and bottom reflector
    top_refl_z = reactor_height + top_reflector_height
    top_refl = openmc.ZPlane(z0=top_refl_z, boundary_type='vacuum')
    bottom_refl_z = -bottom_reflector_height
    bottom_refl = openmc.ZPlane(z0=bottom_refl_z, boundary_type='vacuum')

    top_refl_cell = openmc.Cell(region=hex_prism & +max_z & -top_refl, fill=m_reflector)
    bottom_refl_cell = openmc.Cell(region=hex_prism & -min_z & +bottom_refl, fill=m_reflector)
    top_refl_cell.temperature = specs.inlet_T
    bottom_refl_cell.temperature = coolant_outlet_temp

    model.geometry = openmc.Geometry([outer_cell, top_refl_cell, bottom_refl_cell])

    ### Settings ###
    settings = openmc.Settings()

    settings.particles = 10000
    settings.inactive = n_inactive
    settings.batches = settings.inactive + n_active
    settings.temperature['method'] = 'interpolation'
    settings.temperature['range'] = (294.0, 1500.0)

    l = bundle_pitch / math.sqrt(3.0)
    lower_left = (-l, -l, reactor_bottom)
    upper_right = (l, l, reactor_top)
    source_dist = openmc.stats.Box(lower_left, upper_right, only_fissionable=True)
    source = openmc.Source(space=source_dist)
    settings.source = source

    if (add_entropy_mesh):
        entropy_mesh = openmc.RegularMesh()
        entropy_mesh.lower_left = lower_left
        entropy_mesh.upper_right = upper_right
        entropy_mesh.dimension = (6, 6, 20)
        settings.entropy_mesh = entropy_mesh

    vol_calc = openmc.VolumeCalculation([outer_cell],
                                        100_000_000, lower_left, upper_right)
    settings.volume_calculations = [vol_calc]

    model.settings = settings

    m_colors[mats.m_fuel] = 'palegreen'
    m_colors[mats.m_graphite_c_buffer] = 'sandybrown'
    m_colors[mats.m_graphite_pyc] = 'orange'
    m_colors[mats.m_sic] = 'yellow'
    m_colors[mats.m_graphite_matrix] = 'darkblue'
    m_colors[m_b4c] = 'lightskyblue'

    bundle_p_rounded = int(bundle_pitch)

    plot1          = openmc.Plot()
    plot1.filename = 'plot1'
    plot1.width    = (2 * bundle_pitch, 2 * axial_section_height)
    plot1.basis    = 'xz'
    plot1.origin   = (0.0, 0.0, reactor_height / 2.0)
    plot1.pixels   = (100 * 2 * bundle_p_rounded, int(100 * 3 * axial_section_height))
    plot1.color_by = 'cell'

    plot2          = openmc.Plot()
    plot2.filename = 'plot2'
    plot2.width    = (1.5 * bundle_pitch, 1.5 * bundle_pitch)
    plot2.basis    = 'xy'
    plot2.origin   = (0.0, 0.0, axial_section_height / 4.0)
    plot2.pixels   = (500 * bundle_p_rounded, 500 * bundle_p_rounded)
    plot2.color_by = 'material'
    plot2.colors   = m_colors

    plot3          = openmc.Plot()
    plot3.filename = 'plot3'
    plot3.width    = plot2.width
    plot3.basis    = plot2.basis
    plot3.origin   = plot2.origin
    plot3.pixels   = plot2.pixels
    plot3.color_by = 'cell'

    model.plots = openmc.Plots([plot1, plot2, plot3])

    return model

def main():

    ap = ArgumentParser()
    ap.add_argument('-n', dest='n_axial', type=int, default=50,
                    help='Number of axial cell divisions')
    ap.add_argument('-s', '--entropy', action='store_true',
                    help='Whether to add a Shannon entropy mesh')
    ap.add_argument('-i', dest='n_inactive', type=int, default=500,
                    help='Number of inactive cycles')
    ap.add_argument('-a', dest='n_active', type=int, default=2000,
                    help='Number of active cycles')

    args = ap.parse_args()

    model = assembly(args.n_axial, args.n_inactive, args.n_active, args.entropy)
    model.export_to_xml()

if __name__ == "__main__":
    main()

Cardinal contains several options for the "level" in the geometry at which to apply temperature feedback. Because we are not resolving the TRISO particles in the heat conduction model (and we also are not using a multiscale solid model that might otherwise predict individual temperatures for the different TRISO layers), we will set each TRISO compact to one temperature. That is, all the TRISO particles in cell $var element = document.getElementById("moose-equation-f93fb719-44e8-42b6-b4f0-9feb381be054");katex.render("i", element, {displayMode:false,throwOnError:false});$ , plus the surrounding matrix, will be set to $var element = document.getElementById("moose-equation-51d90257-4677-472d-8d9b-27adb706a5ff");katex.render("T_i", element, {displayMode:false,throwOnError:false});$ . From this selection, the level on which we will apply feedback from MOOSE is 1, because the geometry consists of a hexagonal lattice (level 0), and we want to apply individual cell feedback within that lattice (level 1).

Finally, to accelerate the particle tracking, we:

Repeat the same TRISO universe in each axial layer and within each compact
Superimpose a Cartesian search lattice in the fuel channel regions.

The OpenMC geometry, colored by either cell ID or instance, is shown in Figure 3. Not shown are the axial reflectors on the top and bottom of the assembly. The lateral faces are periodic, while the top and bottom boundaries are vacuum. The Cartesian search lattice in the fuel compact regions is also visible in Figure 3.

Figure 3: OpenMC model, colored by cell ID or instance

Cardinal applies uniform temperature and density feedback to OpenMC for each unique cell ID $var element = document.getElementById("moose-equation-4665b492-685f-451b-8593-8c36f2d34b23");katex.render("+", element, {displayMode:false,throwOnError:false});$ instance combination. For this setup, OpenMC receives on each axial plane a total of 721 temperatures and 108 densities (one density per coolant channel). With references to the colors shown in Figure 3, the 721 cell temperatures correspond to:

$var element = document.getElementById("moose-equation-9844866d-0e8f-4d87-8b1b-5cee025e0e65");katex.render("210\\underbrace{\\text{ fuel compacts}}_{\\substack{\\textup{1 TRISO compact (\\textit{rainbow})}\\\\\\textup{1 matrix region (\\textit{purple})}}}\\ \\ +\\ \\ \\ 108\\underbrace{\\text{ coolant channels}}_{\\substack{\\textup{1 coolant region (\\textit{various})}\\\\\\textup{1 matrix region (\\textit{various})}}}\\ \\ +\\ \\ \\ 6\\underbrace{\\text{ poison compacts}}_{\\substack{\\textup{1 poison region (\\textit{brown})}\\\\\\textup{1 matrix region (\\textit{blue})}}}\\ \\ +\\ \\ \\ 73\\underbrace{\\text{ graphite fillers}}_\\text{1 matrix region (\\textit{mustard})}", element, {displayMode:true,throwOnError:false});$

The solid temperature is provided by the MOOSE heat conduction module, while the fluid temperature and density are provided by THM. Because we will run OpenMC second, the initial fluid temperature is set to the same initial condition that is imposed in the MOOSE solid model. The fluid density is then set using the ideal gas Equation of State (EOS) at a fixed pressure of 7.1 MPa given the imposed temperature, i.e. $var element = document.getElementById("moose-equation-7cc763a9-792f-46d9-bd95-864ad3e5c037");katex.render("\\rho_f(P,T_f)", element, {displayMode:false,throwOnError:false});$ .

To create the XML files required to run OpenMC, if you have the Python API installed, you can run the script:


$ python assembly.py

You can also use the XML files checked in to the assembly directory; if you use these already-existing files.

THM Model

THM solves for conservation of mass, momentum, and energy with 1-D area averages of the Navier-Stokes equations,

$(2)var element = document.getElementById("moose-equation-f0db4487-7cfe-47f4-878f-af8e6e1fe339");katex.render(" \\frac{\\partial}{\\partial t}\\left(A\\rho_f\\right)+\\frac{\\partial}{\\partial x}\\left(A\\rho_fu\\right)=0\\ ,", element, {displayMode:true,throwOnError:false});$

$(3)var element = document.getElementById("moose-equation-8c796cc0-9f8b-4e6a-8914-ca911e0f6e34");katex.render(" \\frac{\\partial}{\\partial t}\\left(A\\rho_fu\\right)+\\frac{\\partial}{\\partial x}\\left(A\\rho_fu^2+AP\\right)=\\tilde{P}\\frac{\\partial A}{\\partial x}-\\frac{f}{2D_h}\\rho_fu|u|A", element, {displayMode:true,throwOnError:false});$

$(4)var element = document.getElementById("moose-equation-2be17bc7-ff70-4683-b174-39fbb260dcdf");katex.render(" \\frac{\\partial}{\\partial t}\\left(A\\rho_f E_f\\right)+\\frac{\\partial}{\\partial x}\\left\\lbrack Au\\left(\\rho_fE_f+P\\right)\\right\\rbrack=H_wa_w\\left(T_\\text{wall}-T_\\text{bulk}\\right)A", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-c3097f74-a1fe-485c-a75d-5d0c6a513f70");katex.render("x", element, {displayMode:false,throwOnError:false});$ is the coordinate along the flow length, $var element = document.getElementById("moose-equation-2f87c3a8-5d06-451f-b35d-8255ff954d2b");katex.render("A", element, {displayMode:false,throwOnError:false});$ is the channel cross-sectional area, $var element = document.getElementById("moose-equation-4b1b39aa-0782-4f81-ac77-9a204b661369");katex.render("\\rho_f", element, {displayMode:false,throwOnError:false});$ is the fluid density, $var element = document.getElementById("moose-equation-a6fdd6b9-2cdb-4289-bb98-f447d7473694");katex.render("u", element, {displayMode:false,throwOnError:false});$ is the $var element = document.getElementById("moose-equation-41ba5e7a-33bb-426e-ab85-4a57fe26d996");katex.render("x", element, {displayMode:false,throwOnError:false});$ -component of velocity, $var element = document.getElementById("moose-equation-c793ffca-536a-4157-8ee7-84494b454e30");katex.render("\\tilde{P}", element, {displayMode:false,throwOnError:false});$ is the average pressure on the curve boundary, $var element = document.getElementById("moose-equation-5ed632e1-50d3-488a-87e9-ee35215cedd3");katex.render("E_f", element, {displayMode:false,throwOnError:false});$ is the fluid total energy, $var element = document.getElementById("moose-equation-0ed30771-cc1d-49c5-9349-71698cc4c5e1");katex.render("f", element, {displayMode:false,throwOnError:false});$ is the friction factor, $var element = document.getElementById("moose-equation-15a9ce7e-ae4e-4ac1-b52c-40c4f6d111e7");katex.render("H_w", element, {displayMode:false,throwOnError:false});$ is the wall heat transfer coefficient, $var element = document.getElementById("moose-equation-d92b6506-7cc4-4014-b265-8a7f51152b30");katex.render("a_w", element, {displayMode:false,throwOnError:false});$ is the heat transfer area density, $var element = document.getElementById("moose-equation-8a5f06bf-e4e8-4167-a62e-ef93052613e5");katex.render("T_\\text{wall}", element, {displayMode:false,throwOnError:false});$ is the wall temperature, and $var element = document.getElementById("moose-equation-16340234-08bb-4431-814e-8e5411569110");katex.render("T_\\text{bulk}", element, {displayMode:false,throwOnError:false});$ is the area average bulk fluid temperature. The Churchill correlation is used for $var element = document.getElementById("moose-equation-196238a3-93c7-403e-9d2c-3ba1d38a52b1");katex.render("f", element, {displayMode:false,throwOnError:false});$ and the Dittus-Boelter correlation is used for $var element = document.getElementById("moose-equation-7a7ba75a-e2d3-45ec-84a7-cea4ff80b547");katex.render("H_w", element, {displayMode:false,throwOnError:false});$ Berry et al. (2014).

The THM mesh for each flow channel contains 50 elements. The mesh is constructed automatically within THM. To simplify the specification of material properties, the fluid geometry uses a length unit of meters. The heat flux imposed in the THM elements is obtained by area averaging the heat flux from the heat conduction model in 50 layers along the fluid-solid interface. For the reverse transfer, the wall temperature sent to MOOSE heat conduction is set to a uniform value along the fluid-solid interface according to a nearest-node mapping to the THM elements.

Because THM will run last in the coupled case, initial conditions are only required for pressure, fluid temperature, and velocity, which are set to uniform distributions.

Multiphysics Coupling

In this section, OpenMC, MOOSE, and THM are coupled for heat source, temperature, and density feedback for the fluid and solid regions of a TRISO-fueled gas reactor assembly. For this system, fluid density and temperature reactivity feedback is quite weak, but is included anyways to be comprehensive. A summary of the data transfers between the three applications is shown in Figure 4. The inset describes the 1-D/3-D data transfers with THM.

Figure 4: Summary of data transfers between OpenMC, MOOSE, and THM

The following sub-sections describe the input files.

Solid Input Files

The solid phase is solved with the MOOSE heat conduction module, and is described in the solid.i input. We define a number of constants at the beginning of the file and set up the mesh from a file.

## HTGR assembly heat conduction input
## Application: MOOSE heat conduction module
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

# copy-pasta from common_input.i
inlet_T = 598.0                          # inlet fluid temperature (K)
buffer_k = 0.5                           # buffer thermal conductivity (W/m/K)
PyC_k = 4.0                              # PyC thermal conductivity (W/m/K)
SiC_k = 13.9                             # SiC thermal conductivity (W/m/K)
kernel_k = 3.5                           # fissil kernel thermal conductivity (W/m/K)
matrix_k = 15.0                          # graphite matrix thermal conductivity (W/m/K)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
triso_pf = 0.15                          # TRISO packing fraction (%)
kernel_radius = 214.85e-6                # fissile kernel outer radius (m)
buffer_radius = 314.85e-6                # buffer outer radius (m)
iPyC_radius = 354.85e-6                  # inner PyC outer radius (m)
SiC_radius = 389.85e-6                   # SiC outer radius (m)
oPyC_radius = 429.85e-6                  # outer PyC outer radius (m)

# compute the volume fraction of each TRISO layer in a TRISO particle
# for use in computing average thermophysical properties
kernel_fraction = ${fparse kernel_radius^3 / oPyC_radius^3}
buffer_fraction = ${fparse (buffer_radius^3 - kernel_radius^3) / oPyC_radius^3}
ipyc_fraction = ${fparse (iPyC_radius^3 - buffer_radius^3) / oPyC_radius^3}
sic_fraction = ${fparse (SiC_radius^3 - iPyC_radius^3) / oPyC_radius^3}
opyc_fraction = ${fparse (oPyC_radius^3 - SiC_radius^3) / oPyC_radius^3}

[Mesh]
  type = FileMesh
  file = meshes/solid_in.e
[]

Next, we define the temperature variable, T, and specify the governing equations and boundary conditions we will apply.

[Variables]
  [T]
    initial_condition = ${inlet_T}
  []
[]

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

[BCs]
  [pin_outer]
    type = MatchedValueBC
    variable = T
    v = thm_temp
    boundary = 'fluid_solid_interface'
  []
[]

The MOOSE heat conduction module will receive power from OpenMC in the form of an AuxVariable, so we define a receiver variable for the fission power, as power. The MOOSE heat conduction module will also receive a fluid wall temperature from THM as another AuxVariable which we name thm_temp. Finally, the MOOSE heat conduction module will send the heat flux to THM, so we add a variable named flux that we will use to compute the heat flux using the DiffusionFluxAux auxiliary kernel.

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

[AuxKernels]
  [flux]
    type = DiffusionFluxAux
    diffusion_variable = T
    component = normal
    diffusivity = thermal_conductivity
    variable = flux
    boundary = 'fluid_solid_interface'
  []
[]

We use functions to define the thermal conductivities. The material properties for the TRISO compacts are taken as volume averages of the various constituent materials. We will evaluate the thermal conductivity for the boron carbide as a function of temperature by using t (which usually is interpeted as time) as a variable to represent temperature. This is syntax supported by the HeatConductionMaterials used to apply these functions to the thermal conductivity.

[Functions]
  [k_graphite]
    type = ParsedFunction
    value = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    value = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    value = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    vars = 'k_TRISO k_graphite'
    vals = 'k_TRISO k_graphite'
  []
  [k_b4c]
    type = ParsedFunction
    value = '5.096154e-6 * t - 1.952360e-2 * t + 2.558435e1'
  []
[]

[Materials]
  [graphite]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_graphite
    temp = T
    block = 'graphite'
  []
  [compacts]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_compacts
    temp = T
    block = 'compacts'
  []
  [poison]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_b4c
    temp = T
    block = 'poison'
  []
[]

We define a number of postprocessors for querying the solution as well as for normalizing the fission power and heat flux, to be described at greater length in Neutronics Input Files .

[Postprocessors]
  [flux_integral]
    # evaluate the total heat flux for normalization
    type = SideDiffusiveFluxIntegral
    diffusivity = thermal_conductivity
    variable = T
    boundary = 'fluid_solid_interface'
  []
  [max_fuel_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'compacts'
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [power]
    # evaluate the total power for normalization
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = 'compacts'
    execute_on = 'transfer'
  []
[]

For visualization purposes only, we add LayeredAverages for the fuel and block temperatures. These will average the temperature in layers oriented in the $var element = document.getElementById("moose-equation-e737b3a1-a0e2-4134-b583-8bb0b7ae3f78");katex.render("z", element, {displayMode:false,throwOnError:false});$ direction, which we will use for plotting axial temperature distributions. We output the results of these userobjects to CSV using SpatialUserObjectVectorPostprocessors and by setting csv = true in the output.

[UserObjects]
  [average_fuel_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_block_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'graphite'
  []
[]

[VectorPostprocessors]
  [fuel_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fuel_axial
  []
  [block_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_block_axial
  []
[]

[Outputs]
  exodus = true
  csv = true
  print_linear_residuals = false
[]

Finally, we specify a Transient executioner. Because there are no time-dependent kernels in this input file, this is equivalent in practice to using a Steady executioner, but allows you to potentially sub-cycle the MOOSE heat conduction solve relative to the OpenMC solve (such as if you wanted to converge the Conjugate Heat Transfer (CHT) fully inbetween data exchanges with OpenMC).

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

Fluid Input Files

The fluid phase is solved with THM, and is described in the thm.i input. This input file is built using syntax specific to THM - we will only briefly cover this syntax, and instead refer users to the THM manuals for more information. First we define a number of constants at the beginning of the file and apply some global settings. We set the initial conditions for pressure, velocity, and temperature and indicate the fluid EOS object using IdealGasFluidProperties.

## HTGR assembly thermal hydraulics simulation
## Application: MOOSE thermal hydraulics module
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

# copy-pasta from common.i
inlet_T = 598.0 # inlet fluid temperature (K)
mdot = '${fparse 117.3 / 12 / 108}' # fluid mass flowrate (kg/s)
outlet_P = 7.1e6 # fluid outlet pressure (Pa)
channel_diameter = 0.016 # diameter of the coolant channels (m)
height = 6.343 # height of the assembly (m)

num_layers_for_THM = 50 # number of elements in the THM model

[GlobalParams]
  initial_p = ${outlet_P}
  initial_T = ${inlet_T}
  initial_vel = '${fparse mdot / outlet_P / 8.3144598 * 4.0e-3 / inlet_T / (pi * channel_diameter * channel_diameter / 4.0)}'

  rdg_slope_reconstruction = full
  closures = none
  fp = helium
[]

[Closures]
  [none]
    type = Closures1PhaseNone
  []
[]

[FluidProperties]
  [helium]
    type = IdealGasFluidProperties
    molar_mass = 4e-3
    gamma = 1.668282 # should correspond to  Cp = 5189 J/kg/K
    k = 0.2556
    mu = 3.22639e-5
  []
[]

Next, we define the "components" in the domain. These components essentially consist of the physics equations and boundary conditions solved by THM, but expressed in THM-specific syntax. These components define single-phase flow in a pipe, an inlet mass flowrate boundary condition, an outlet pressure boundary condition, and heat transfer to the pipe wall.

[Components]
  [channel]
    type = FlowChannel1Phase
    position = '0 0 ${height}'
    orientation = '0 0 -1'

    A = '${fparse pi * channel_diameter * channel_diameter / 4}'
    D_h = ${channel_diameter}
    length = ${height}
    n_elems = ${num_layers_for_THM}
  []

  [inlet]
    type = InletMassFlowRateTemperature1Phase
    input = 'channel:in'
    m_dot = ${mdot}
    T = ${inlet_T}
  []

  [outlet]
    type = Outlet1Phase
    input = 'channel:out'
    p = ${outlet_P}
  []

  [ht_ext]
    type = HeatTransferFromExternalAppHeatFlux1Phase
    flow_channel = channel
    P_hf = '${fparse channel_diameter * pi}'
  []
[]

Associated with these components are a number of closures, defined as materials. We set up the Churchill correlation for the friction factor and the Dittus-Boelter correlation for the convective heat transfer coefficient. Additional materials are created to represent dimensionless numbers and other auxiliary terms, such as the wall temperature. As can be seen here, the Material system is not always used to represent quantities traditionally thought of as "material properties."

[Materials]
  # wall friction closure
  [f_mat]
    type = ADWallFrictionChurchillMaterial
    block = channel
    D_h = D_h
    f_D = f_D
    vel = vel
    rho = rho
    mu = mu
  []

  # Wall heat transfer closure (all important is in Nu_mat)
  [Re_mat]
    type = ADReynoldsNumberMaterial
    block = channel
    Re = Re
    D_h = D_h
    mu = mu
    vel = vel
    rho = rho
  []
  [Pr_mat]
    type = ADPrandtlNumberMaterial
    block = channel
    cp = cp
    mu = mu
    k = k
  []
  [Nu_mat]
    type = ADParsedMaterial
    block = channel
    # Dittus-Boelter
    function = '0.022 * pow(Re, 0.8) * pow(Pr, 0.4)'
    f_name = 'Nu'
    material_property_names = 'Re Pr'
  []
  [Hw_mat]
    type = ADConvectiveHeatTransferCoefficientMaterial
    block = channel
    D_h = D_h
    Nu = Nu
    k = k
  []
  [T_wall]
    type = ADTemperatureWall3EqnMaterial
    Hw = Hw
    T = T
    q_wall = q_wall
  []
[]

THM computes the wall temperature to apply a boundary condition in the MOOSE heat conduction module. To convert the T_wall material into an auxiliary variable, we use the ADMaterialRealAux.

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

[AuxKernels]
  [Tw_aux]
    type = ADMaterialRealAux
    block = channel
    variable = T_wall
    property = T_wall
  []
[]

Finally, we set the preconditioner, a Transient executioner, and set an Exodus output. The steady_state_detection and steady_state_tolerance parameters will automatically terminate the THM solution once the relative change in the solution is less than $var element = document.getElementById("moose-equation-c225d78b-3217-40b3-97f4-5dbdc290bb4b");katex.render("10^{-8}", element, {displayMode:false,throwOnError:false});$ .

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

[Executioner]
  type = Transient
  dt = 0.1
  start_time = 0

  steady_state_detection = true
  steady_state_tolerance = 1e-08

  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-6

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  solve_type = NEWTON
  line_search = basic
[]

[Outputs]
  exodus = true
  print_linear_residuals = false

  [console]
    type = Console
    outlier_variable_norms = false
  []
[]

As you may notice, this THM input file only models a single coolant flow channel. We will leverage a feature in MOOSE that allows a single application to be repeated multiple times throughout a master application without having to merge input files or perform other transformations. We will run OpenMC as the main application; the syntax needed to achieve this setup is covered next.

Neutronics Input Files

The neutronics physics is solved with OpenMC over the entire domain. The OpenMC wrapping is described in the openmc.i file. We begin by defining a number of constants and by setting up the mesh mirror on which OpenMC will receive temperature and density from the coupled applications, and on which OpenMC will write the fission heat source. Because the coupled MOOSE applications use length units of meters, the mesh mirror must also be in units of meters in order to obtain correct data transfers. For simplicity, we use the same solid mesh as used by the solid heat conduction solution, but this is not required. For the fluid region, we use MOOSE mesh generators to construct a mesh for a single coolant channel, and then repeat it for the 108 coolant channels.

## HTGR assembly neutronics simulation
## Application: Cardinal
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

num_layers_for_THM = 50 # number of elements in the THM model; for the converged
# case, we set this to 150

[GlobalParams]
  # Reduces transfers efficiency for now, can be removed once transferred fields are checked
  bbox_factor = 10
[]

[Mesh]
  # mesh mirror for the solid regions
  [solid]
    type = FileMeshGenerator
    file = meshes/solid_in.e
  []

  # create a mesh for a single coolant channel; because we will receive uniform
  # temperatures and densities from THM on each x-y plane, we can use a very coarse
  # mesh in the radial direction
  [coolant_face]
    type = AnnularMeshGenerator
    nr = 1
    nt = 8
    rmin = 0.0
    rmax = '${fparse channel_diameter / 2.0}'
  []
  [extrude]
    type = AdvancedExtruderGenerator
    input = coolant_face
    num_layers = ${num_layers_for_THM}
    direction = '0 0 1'
    heights = '${height}'
    top_boundary = '300' # inlet
    bottom_boundary = '400' # outlet
  []
  [rename]
    type = RenameBlockGenerator
    input = extrude
    old_block = '1'
    new_block = '101'
  []

  # repeat the coolant channels and then combine together to get a combined mesh mirror
  [repeat]
    type = CombinerGenerator
    inputs = rename
    positions_file = coolant_channel_positions.txt
  []
  [add]
    type = CombinerGenerator
    inputs = 'solid repeat'
  []
[]

Before progressing further, we first need to describe the multiapp structure connecting OpenMC, MOOSE heat conduction, and THM. We set the main application to OpenMC, and will have both MOOSE heat conduction and THM as "sibling" sub-applications. Therefore, we will see in the OpenMC input file information related to data transfers between MOOSE heat conduction and THM. The auxiliary variables defined for the OpenMC model are shown below.

[AuxVariables]
  [material_id]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_density]
    family = MONOMIAL
    order = CONSTANT
  []
  [thm_temp_wall]
    block = '101'
  []
  [flux]
  []

  # just for postprocessing purposes
  [thm_pressure]
    block = '101'
  []
  [thm_velocity]
    block = '101'
  []
  [z]
    family = MONOMIAL
    order = CONSTANT
    block = 'compacts'
  []
[]

[AuxKernels]
  [material_id]
    type = CellMaterialIDAux
    variable = material_id
  []
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
  [density]
    type = FluidDensityAux
    variable = density
    p = ${outlet_P}
    T = thm_temp
    fp = helium
    execute_on = 'timestep_begin linear'
  []
  [z]
    type = ParsedAux
    variable = z
    use_xyzt = true
    expression = 'z'
  []
[]

For visualization purposes, we use a CellTemperatureAux to view the temperature set in each OpenMC cell and a CellDensityAux to view the density set in each fluid OpenMC cell. To understand how the OpenMC model maps to the [Mesh], we also include CellMaterialIDAux, CellIDAux, and CellInstanceAux auxiliary kernels. Next, we add a receiver flux variable that will hold the heat flux received from MOOSE (and sent to THM) and another receiver variable thm_temp_wall that will hold the wall temperature received from THM (and sent to MOOSE). The OpenMC wrapping in Cardinal automatically adds auxiliary variables named temp and density when receiving feedback from coupled applications.

Finally, to reduce the number of transfers from THM, we will receive fluid temperature from THM, but re-compute the density locally in the OpenMC wrapping using a FluidDensityAux with the same EOS as used in the THM input files.

[FluidProperties]
  [helium]
    type = IdealGasFluidProperties
    molar_mass = 4e-3
    gamma = 1.668282 # should correspond to  Cp = 5189 J/kg/K
    k = 0.2556
    mu = 3.22639e-5
  []
[]

Next, we set initial conditions for the fluid wall temperature, the fluid bulk temperature, and the heat source. We set these initial conditions in the OpenMC wrapper because the very first time that the transfers to the MOOSE heat conduction module and to THM occur, these initial conditions will be passed.

[ICs]
  [fluid_temp_wall]
    type = FunctionIC
    variable = thm_temp_wall
    function = temp_ic
  []
  [fluid_temp]
    type = FunctionIC
    variable = thm_temp
    function = temp_ic
  []
  [heat_source]
    type = ConstantIC
    variable = heat_source
    block = 'compacts'
    value = '${fparse power / (n_bundles * n_fuel_compacts_per_block) / (pi * compact_diameter * compact_diameter / 4.0 * height)}'
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + (${height} - z) / ${height} * ${power} / ${mdot} / ${fluid_Cp}'
  []
[]

The [Problem] block is then used to specify settings for the OpenMC wrapping. We define the total power for normalization, indicate that blocks 1, 2, and 4 are solid (graphite, compacts, and poison) while block 101 is fluid. We automatically add tallies to block 2, the fuel compacts. Because OpenMC solves in units of centimeters, we specify a scaling of 100, i.e. a multiplicative factor to apply to the [Mesh] to get into OpenMC's centimeter units.

[Problem]
  type = OpenMCCellAverageProblem

  # optimizations to increase tracking rate by assuring that the tallies
  # are spatially independent
  check_tally_sum = false
  normalize_by_global_tally = false
  assume_separate_tallies = true

  power = '${fparse power / n_bundles}'
  scaling = 100.0

  relaxation = constant
  relaxation_factor = 0.5

  inactive_batches = 200
  batches = 1000

  # we will collate temperature from THM (for the fluid) and MOOSE (for the solid)
  # into variables we name as 'solid_temp' and 'thm_temp'. This syntax will automatically
  # create those variabes for us
  temperature_variables = 'solid_temp; solid_temp; solid_temp; thm_temp'
  temperature_blocks = 'graphite; compacts; poison; 101'
  cell_level = 0

  density_blocks = '101'

  [Tallies]
    [heat_source]
      type = CellTally
      block = '2'
      name = heat_source
      score = heating

      trigger = rel_err
      trigger_threshold = 2e-2

      check_equal_mapped_tally_volumes = true
      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

Other features we use include an output of the fission tally standard deviation in units of W/m $var element = document.getElementById("moose-equation-d068b1f3-136b-4cd7-8e9e-6d51698295ce");katex.render("^3", element, {displayMode:false,throwOnError:false});$ to the [Mesh] by setting output = 'fission_tally_std_dev'. This is used to obtain uncertainty estimates of the heat source distribution from OpenMC in the same units as the heat source. We also leverage a helper utility in Cardinal by setting check_equal_mapped_tally_volumes = true. This parameter will throw an error if the tallied OpenMC cells map to different volumes in the MOOSE domain. Because we know a priori that the equal-volume OpenMC tally cells should all map to equal volumes, this will help ensure that the volumes used for heat source normalization are also all equal. For further discussion of this setting and a pictorial description of the possible effect of non-equal mapped volumes, please see the OpenMCCellAverageProblem documentation.

Because the blocks in the OpenMC mesh mirror receive temperatures from different applications, we use the temperature_variables and temperature_blocks parameters of OpenMCCellAverageProblem to automatically create separate variables to hold the temperatures from THM and MOOSE (thm_temp, solid_temp) and create several SelfAux auxiliary kernels to write into temp. The temperature_blocks and temperature_variables parameters simply allow shorter input syntax by creating the variables and their auxiliary kernels for you.

Finally, we apply a constant relaxation model to the heat source. A constant relaxation will compute the heat source in iteration $var element = document.getElementById("moose-equation-7a000d8d-7a4d-466c-8df8-b1902b1b4e80");katex.render("i+1", element, {displayMode:false,throwOnError:false});$ as an average of the heat source from iteration $var element = document.getElementById("moose-equation-f83ec019-b881-4e06-a904-4ff3ce9c769f");katex.render("i", element, {displayMode:false,throwOnError:false});$ and the most-recently-computed heat source, indicated here as a generic operator $var element = document.getElementById("moose-equation-b4e14b22-719c-4a50-b325-ab12c83881fb");katex.render("\\Phi", element, {displayMode:false,throwOnError:false});$ which represents the Monte Carlo solve:

$(5)var element = document.getElementById("moose-equation-09081ddb-d9ff-48f6-aa5b-4729abf03ade");katex.render(" q^{i+1}=\\left(1-\\alpha\\right)q^i+\\alpha\\Phi", element, {displayMode:true,throwOnError:false});$

In Eq. (5), $var element = document.getElementById("moose-equation-798494ee-b880-4fcc-8361-578f39611038");katex.render("\\alpha", element, {displayMode:false,throwOnError:false});$ is the relaxation factor, which we set here to 0.5. In other words, the heat source in iteration $var element = document.getElementById("moose-equation-af46d147-0950-4bf1-9a11-8c834fa2b696");katex.render("i+1", element, {displayMode:false,throwOnError:false});$ is an average of the most recent Monte Carlo solution and the previous iterate. This is necessary to accelerate the fixed point iterations and achieve convergence in a reasonable time - otherwise oscillations can occur in the coupled physics.

We run OpenMC as the main application, so we next need to define MultiApps to run the solid heat conduction model and the THM fluid model as the sub-applications. We also require a number of transfers both for 1) sending necessary coupling data between the three applications and 2) visualizing the combined THM output. To couple OpenMC to MOOSE heat conduction, we use four transfers:

MultiAppShapeEvaluationTransfer to transfer power from OpenMC to MOOSE (with conservation of total power)
MultiAppInterpolationTransfer to transfer:
- solid temperature from MOOSE to OpenMC
- wall temperature from OpenMC (which doesn't directly compute the wall temperature, but instead receives it from THM through a separate transfer) to MOOSE
MultiAppNearestNodeTransfer to transfer and conserve heat flux from MOOSE to OpenMC (which isn't used directly in OpenMC, but instead gets sent later to THM through a separate transfer)

To couple OpenMC to THM, we require three transfers:

MultiAppUserObjectTransfer to send the layer-averaged wall heat flux from OpenMC (which computes the layered-average heat flux from the heat flux received from MOOSE heat conduction) to THM
MultiAppNearestNodeTransfer to transfer:
- fluid wall temperature from THM to OpenMC (which isn't used directly in OpenMC, but instead gets sent to MOOSE heat conduction in a separate transfer)
- fluid bulk temperature from THM to OpenMC

For visualization purposes, we also send the pressure and velocity computed by THM to the OpenMC mesh mirror.

[MultiApps]
  [bison]
    type = TransientMultiApp
    app_type = CardinalApp
    input_files = 'solid.i'
    execute_on = timestep_begin
  []
  [thm]
    type = FullSolveMultiApp
    app_type = CardinalApp
    input_files = 'thm.i'
    execute_on = timestep_end
    max_procs_per_app = 1
    bounding_box_padding = '0.1 0.1 0'
    positions_file = coolant_channel_positions.txt
    output_in_position = true
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [heat_flux_to_openmc]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = flux
    variable = flux
    from_multi_app = bison
    source_boundary = 'fluid_solid_interface'
    target_boundary = 'fluid_solid_interface'
    from_postprocessors_to_be_preserved = flux_integral
    to_postprocessors_to_be_preserved = flux_integral
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    direction = to_multiapp
    multi_app = bison
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = power
  []
  [thm_temp_to_bison]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = thm_temp_wall
    variable = thm_temp
    direction = to_multiapp
    multi_app = bison
  []

  [q_wall_to_thm]
    type = MultiAppGeneralFieldUserObjectTransfer
    variable = q_wall
    direction = to_multiapp
    multi_app = thm
    source_user_object = q_wall_avg
  []
  [T_wall_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = T_wall
    direction = from_multiapp
    multi_app = thm
    variable = thm_temp_wall
    target_boundary = 'fluid_solid_interface'
  []
  [T_bulk_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = T
    direction = from_multiapp
    multi_app = thm
    variable = thm_temp
  []

  # just for postprocessing purposes
  [pressure_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = p
    direction = from_multiapp
    multi_app = thm
    variable = thm_pressure
  []
  [velocity_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = vel_z
    direction = from_multiapp
    multi_app = thm
    variable = thm_velocity
  []
[]

To compute the layer-averaged heat flux on the surface of each coolant channel (which is used as a boundary condition in THM), we use a NearestPointLayeredSideAverage user object, where by providing the center points of each of the coolant channels, we can get a unique heat flux along each channel wall. We also add several LayeredAverage user objects in order to compute radially-averaged power, temperatures, pressures, and velocities that we will use later in making axial plots of the solution. We can automatically output these user objects into CSV format by translating the user objects into SpatialUserObjectVectorPostprocessors.

[UserObjects]
  [q_wall_avg]
    type = NearestPointLayeredSideAverage
    boundary = fluid_solid_interface
    variable = flux

    # Note: make this to match the num_elems in the channel
    direction = z
    num_layers = ${num_layers_for_THM}
    points_file = coolant_channel_positions.txt

    direction_min = 0.0
    direction_max = ${height}
  []
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_fluid_axial]
    type = LayeredAverage
    variable = thm_temp
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
  [average_pressure]
    type = LayeredAverage
    variable = thm_pressure
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
  [average_axial_velocity]
    type = LayeredAverage
    variable = thm_velocity
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
  [fluid_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fluid_axial
  []
  [pressure_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_pressure
  []
  [velocity_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_axial_velocity
  []
[]

Finally, we use a Transient executioner and specify Exodus and CSV output formats. Note that the time step size is inconsequential in this case, but instead represents the Picard iteration. We will run for 10 "time steps," which represent fixed point iterations.

[Executioner]
  type = Transient
  num_steps = 10
[]

[Outputs]
  exodus = true
  csv = true
  hide = 'P_in flux_integral z'
[]

Execution and Postprocessing

To run the coupled calculation, run the following:


$ mpiexec -np 6 cardinal-opt -i common_input.i openmc.i --n-threads=12

This will run with 6 Message Passing Interface (MPI) processes and 12 OpenMP threads (you may use other parallel configurations as needed). This tutorial uses quite large meshes due to the 6 meter height of the domain - if you wish to run this tutorial with fewer computational resources, you can reduce the height and the various mesh parameters (number of extruded layers and elements in the THM domain) and then recreate the OpenMC model and meshes.

When the simulation has completed, you will have created a number of different output files:

openmc_out.e, an Exodus file with the OpenMC solution and the data that was ultimately transferred in/out of OpenMC
openmc_out_bison0.e, an Exodus file with the solid solution
openmc_out_thm<n>.e, Exodus files with each of the <n> THM solutions
openmc_out.csv, a CSV file with the results of the postprocessors in the OpenMC wrapping input file for each time step
openmc_out_bison0.csv, a CSV file with the results of the postprocessors in the solid heat conduction input file for each time step
openmc_out_bison0_block_axial_avg_<n>.csv, a CSV file with the layer-averaged block temperature at time step <n>
openmc_out_bison0_flux_axial_avg_<n>.csv, a CSV file with the layer-averaged fluid-solid interface heat flux at time step <n>
openmc_out_bison0_fuel_axial_avg_<n>.csv, a CSV file with the layer-averaged compact temperature at time step <n>
openmc_out_fluid_avg_<n>.csv, a CSV file with the layer-averaged fluid bulk temperature at time step <n>
openmc_out_power_avg_<n>.csv, a CSV file with the layer-averaged heat source at time step <n>
openmc_out_pressure_avg_<n>.csv, a CSV file with the layer-averaged pressure at time step <n>
openmc_out_velocity_avg_<n>.csv, a CSV file with the layer-average axial fluid velocity at time step <n>

Coupled convergence is defined at the point when 1) the $var element = document.getElementById("moose-equation-ba3660a1-97d1-4a3e-b6c5-985b62994d22");katex.render("k", element, {displayMode:false,throwOnError:false});$ eigenvalue is within the uncertainty band of the previous iteration and 2) there is less than 2 K absolute change in maximum fuel, block, and fluid temperatures. Convergence is obtained in 6 fixed point iterations. Figure 5 shows the maximum fuel, compact, and fluid temperatures, along with the $var element = document.getElementById("moose-equation-d0137971-af72-4b14-b1e9-b45e375cb2ee");katex.render("k", element, {displayMode:false,throwOnError:false});$ eigenvalue, as a function of iteration number for all 10 iterations run. For each iteration, we first run MOOSE heat conduction (which we indicate as iterations A1, A2, and so on), followed by OpenMC (which we indicate as iterations B1, B2, and so on), and finally THM (which we indicate as iterations C1, C2, and so on). In other words, in iteration $var element = document.getElementById("moose-equation-0157a84b-86b9-49aa-b40c-db02fefc73c7");katex.render("i", element, {displayMode:false,throwOnError:false});$ :

Iteration A $var element = document.getElementById("moose-equation-826e4971-1ae1-43de-8ea6-d69f00e14053");katex.render("i", element, {displayMode:false,throwOnError:false});$ represents a MOOSE heat conduction solve using the power and fluid-solid wall temperature from iteration $var element = document.getElementById("moose-equation-31369fb7-8858-4631-af26-ce9fc9b60c71");katex.render("i-1", element, {displayMode:false,throwOnError:false});$
Iteration B $var element = document.getElementById("moose-equation-fbffa8b5-b53a-43ed-83bc-1f3c5d18187f");katex.render("i", element, {displayMode:false,throwOnError:false});$ represents an OpenMC solve using the solid temperature from iteration $var element = document.getElementById("moose-equation-573d719e-211a-44e8-bdb9-0e4bdb379cad");katex.render("i", element, {displayMode:false,throwOnError:false});$ and the fluid temperature and density from iteration $var element = document.getElementById("moose-equation-7a830f15-0cc2-411e-ab72-a59075cca8e5");katex.render("i-1", element, {displayMode:false,throwOnError:false});$
Iteration C $var element = document.getElementById("moose-equation-75f9903a-1785-4625-9f57-da82b83067df");katex.render("i", element, {displayMode:false,throwOnError:false});$ represents a THM solve using the fluid-solid wall heat flux from iteration $var element = document.getElementById("moose-equation-e35f7b3f-e42c-4fcd-830c-c28756995ce4");katex.render("i", element, {displayMode:false,throwOnError:false});$

Figure 5: Convergence in various metrics as a function of iteration number

Figure 6 shows the radially-averaged power from OpenMC as a function of iteration number. There is essentially no change in the axial distribution beyond 6 fixed point iterations, which further confirms that we have obtained a converged solution. The remainder of the depicted results correpond to iteration 6.

Figure 6: Radially-average fission power distribution as a function of iteration number $var element = document.getElementById("moose-equation-7d7086cd-4f3f-4dd5-8084-805bf2654a39");katex.render("i", element, {displayMode:false,throwOnError:false});$

Figure 7 shows the fission power distribution computed by OpenMC. The inset shows the power distribution on an $var element = document.getElementById("moose-equation-f5c27d86-1b97-4885-8349-a5211972362b");katex.render("x-y", element, {displayMode:false,throwOnError:false});$ plane 2.15 meters from the inlet, where the maximum power occurs. Slight azimuthal asymmetries exist due to the finite tally uncertainty. Neutrons reflecting back into the fuel region from the axial reflectors cause local power peaking at the ends of the assembly, while the negative fuel temperature coefficient causes the power distribution to shift towards the reactor inlet where temperatures are lower. The six poison compacts in the corners result in local power depression such that the highest compact powers occur near the center of the assembly.

Figure 7: Fission power predicted by OpenMC; note the use of a separate color scale for the inset.

Figure 8 shows the solid temperature (left) computed by MOOSE and the cell temperature imposed in OpenMC (right). The bottom row shows the temperature on an $var element = document.getElementById("moose-equation-93e57011-a00d-4c25-95d4-261efbf85010");katex.render("x-y", element, {displayMode:false,throwOnError:false});$ plane 2.15 meters from the inlet. The insulated boundary conditions, combined with a lower "density" of coolant channels near the lateral faces, result in higher compact temperatures near the assembly peripheries, except in the vicinity of the poison pins where the increased absorption reduces the local power densities. Each OpenMC cell is set to the volume-average temperature from the mesh mirror elements whose centroids mapped to each cell; a similar procedure is used to set cell temperatures and densities for the fluid cells.

Figure 8: Solid temperature predicted by MOOSE (left) and volume-average temperature imposed in OpenMC (right). Note the use of a separate color scale for the insets.

Figure 9 shows the solid temperature computed by MOOSE on several $var element = document.getElementById("moose-equation-c4f091c7-1cff-426a-99f6-1d617152ad3b");katex.render("x-y", element, {displayMode:false,throwOnError:false});$ planes with the fluid temperature computed by THM shown as tubes. An inset shows the fluid temperature on the outlet plane. The absence of compacts in the center region results in the lowest fluid temperatures in this region, while the highest fluid temperatures are observed for channels surrounded by 6 compacts that are sufficiently close to the periphery to be affected by the lateral insulated boundary conditions.

Figure 9: Fluid temperature predicted by THM (tubes and inset) and solid temperature predicted by MOOSE (five slices). Note the use of three separate color scales.

Finally, Fig. Figure 10 shows the radially-averaged fission distribution and fluid, compact, and graphite temperatures (left); and velocity and pressure (right) as a function of axial position. The negative temperature feedback results in a top-peaked power distribution. The fuel temperature peaks near the mid-plane due to the combined effects of the relatively high power density and the continually-increasing fluid temperature with distance from the inlet. The pressure gradient is nearly constant with axial position. Due to mass conservation, the heating of the fluid results in the velocity increasing with distance from the inlet.

Figure 10: Radially-averaged temperatures, power, pressure, and velocity as a function of position.

References

R. A. Berry, J. W. Peterson, H. Zhang, R. C. Martineau, H. Zhao, L. Zou, and D. Andrš. RELAP-7 Theory Manual. Technical Report INL/EXT-14-31366, Idaho National Laboratory, 2014.

@TechReport{relap7,
    author = "Berry, R. A. and Peterson, J. W. and Zhang, H. and Martineau, R. C. and Zhao, H. and Zou, L. and Andr\v{s}, D.",
    title = "{{RELAP-7} Theory Manual}",
    institution = "Idaho National Laboratory",
    year = "2014",
    number = "INL/EXT-14-31366"
}

INL. High-Temperature Gas-Cooled Test Reactor Point Design. Technical Report INL/EXT-16-38296, Idaho National Laboratory, 2016. URL: https://tinyurl.com/v9a4hym5.

@TechReport{sterbentz,
    author = "INL",
    title = "{High-Temperature Gas-Cooled Test Reactor Point Design}",
    institution = "Idaho National Laboratory",
    year = "2016",
    number = "INL/EXT-16-38296",
    url = "https://tinyurl.com/v9a4hym5"
}

A.J. Novak, D. Andrs, P. Shriwise, D. Shaver, P.K. Romano, E. Merzari, and P. Keutelian. Coupled Monte Carlo and Thermal-Hydraulics Modeling of a Prismatic Gas Reactor Fuel Assembly Using Cardinal. In Proceedings of Physor. 2021.

@InProceedings{novak_2021c,
    author = "Novak, A.J. and Andrs, D. and Shriwise, P. and Shaver, D. and Romano, P.K. and Merzari, E. and Keutelian, P.",
    title = "{Coupled {Monte} {Carlo} and Thermal-Hydraulics Modeling of a Prismatic Gas Reactor Fuel Assembly Using {Cardinal}}",
    booktitle = "{Proceedings of Physor}",
    year = "2021"
}

E. Shemon, K. Mo, Y. Miao, Y. Jung, S. Richards, A. Oaks, and S. Kumar. MOOSE Framework Enhancements for Meshing Reactor Geometries. In Proceedings of Physor. 2022.

@InProceedings{shemon_2021,
    author = "Shemon, E. and Mo, K. and Miao, Y. and Jung, Y. and Richards, S. and Oaks, A. and Kumar, S.",
    title = "{{MOOSE} Framework Enhancements for Meshing Reactor Geometries}",
    booktitle = "{Proceedings of Physor}",
    year = "2022"
}

(htgr/assembly/meshes/solid.i)

# HTGR assembly meshing file
# Application: MOOSE reactor module, Griffin, BlueCRAB
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html

# This input file is used to create the mesh for the solid phase.
# This should be run with:
#
# cardinal-opt -i common_input.i solid.i --mesh-only

n_layers = 100

[GlobalParams]
  quad_center_elements = true
[]

[Mesh]
  [fuel_pin]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    polygon_size = '${fparse fuel_to_coolant_distance / 2.0}'
    ring_radii = '${fparse 0.8 * compact_diameter / 2.0} ${fparse compact_diameter / 2.0}'
    ring_intervals = '1 1'
    num_sectors_per_side = '4 4 4 4 4 4'
    ring_block_ids = '2 2'
    ring_block_names = 'compacts compacts'
    background_block_ids = '1'
    background_block_names = 'graphite'
    background_intervals = 2
  []
  [coolant_pin]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    polygon_size = '${fparse fuel_to_coolant_distance / 2.0}'
    ring_radii = '${fparse channel_diameter / 2.0}'
    ring_intervals = '2'
    num_sectors_per_side = '4 4 4 4 4 4'
    ring_block_ids = '101 101'
    ring_block_names = 'coolant coolant'
    background_block_ids = '1'
    background_block_names = 'graphite'
    interface_boundary_id_shift = 100
    background_intervals = 2
    create_inward_interface_boundaries = true
  []
  [poison_pin]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    polygon_size = '${fparse fuel_to_coolant_distance / 2.0}'
    ring_radii = '${fparse compact_diameter / 2.0}'
    ring_intervals = '2'
    num_sectors_per_side = '4 4 4 4 4 4'
    ring_block_ids = '4 4'
    ring_block_names = 'poison poison'
    background_block_ids = '1'
    background_block_names = 'graphite'
    background_intervals = 2
  []
  [graphite_pin]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    polygon_size = '${fparse fuel_to_coolant_distance / 2.0}'
    ring_radii = '${fparse compact_diameter / 2.0}'
    ring_intervals = '2'
    num_sectors_per_side = '4 4 4 4 4 4'
    ring_block_ids = '1 1'
    ring_block_names = 'graphite graphite'
    background_block_ids = '1'
    background_block_names = 'graphite'
  []
  [bundle]
    type = PatternedHexMeshGenerator
    inputs = 'fuel_pin coolant_pin poison_pin graphite_pin'
    hexagon_size = '${fparse bundle_flat_to_flat / 2.0 + bundle_gap_width / 2.0}'
    pattern = '2 0 1 0 0 1 0 0 1 0 2;
              0 1 0 0 1 0 0 1 0 0 1 0;
             1 0 0 1 0 0 1 0 0 1 0 0 1;
            0 0 1 0 0 1 0 0 1 0 0 1 0 0;
           0 1 0 0 1 0 0 1 0 0 1 0 0 1 0;
          1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1;
         0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0;
        0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0;
       1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1;
      0 0 1 0 0 1 0 0 1 3 3 1 0 0 1 0 0 1 0 0;
     2 1 0 0 1 0 0 1 0 3 3 3 0 1 0 0 1 0 0 1 2;
      0 0 1 0 0 1 0 0 1 3 3 1 0 0 1 0 0 1 0 0;
       1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1;
        0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0;
         0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0;
          1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1;
           0 1 0 0 1 0 0 1 0 0 1 0 0 1 0;
            0 0 1 0 0 1 0 0 1 0 0 1 0 0;
             1 0 0 1 0 0 1 0 0 1 0 0 1;
              0 1 0 0 1 0 0 1 0 0 1 0;
               2 0 1 0 0 1 0 0 1 0 2'
    rotate_angle = 0

    background_block_id = '1'
    background_block_names = 'graphite'
  []
  [extrude]
    type = AdvancedExtruderGenerator
    input = bundle
    heights = ${height}
    num_layers = ${n_layers}
    direction = '0 0 1'
  []
  [delete_coolant]
    type = BlockDeletionGenerator
    input = extrude
    block = '101'
  []
  [rename_coolant_sideset]
    type = RenameBoundaryGenerator
    input = delete_coolant
    old_boundary = 102
    new_boundary = 'fluid_solid_interface'
  []

  construct_side_list_from_node_list = true
[]

(htgr/assembly/assembly.py)

#!/bin/env python

from argparse import ArgumentParser
import math
import numpy as np
import matplotlib.pyplot as plt

import openmc
import sys
import os

# Get common input parameters shared by other physics
script_dir = os.path.dirname(__file__)
sys.path.append(script_dir)
import common_input as specs
import materials as mats

def coolant_temp(t_in, t_out, l, z):
    """
    Computes the coolant temperature based on an expected cosine power distribution
    for a specified temperature rise. The total core temperature rise is governed
    by energy conservation as dT = Q / m / Cp, where dT is the total core temperature
    rise, Q is the total core power, m is the mass flowrate, and Cp is the fluid
    isobaric specific heat. If you neglect axial heat conduction and assume steady
    state, then the temperature rise in a layer of fluid i can be related to the
    ratio of the power in that layer to the total power,
    dT_i / dT = Q_i / Q. We assume here a sinusoidal power distribution to get
    a reasonable estimate of an initial coolant temperature distribution.

    Parameters
    ----------

    t_in : float
        Inlet temperature of the channel
    t_out : float
        Outlet temperature of the channel
    l : float
        Length of the channel
    z : float or 1-D numpy.array
        Axial position where the temperature will be computed

    Returns
    -------
        float or 1-D numpy array of float depending on z
    """
    dT = t_out - t_in
    Q = 2 * l / math.pi
    Qi = (l - l * np.cos(math.pi * z / l)) / math.pi

    t = t_in + Qi / Q * dT

    return t

def coolant_density(t):
  """
  Computes the helium density from temperature assuming a fixed operating pressure.

  Parameters
  ----------

  t : float
    Fluid temperature

  Returns
  _______
    float or 1-D numpy array of float depending on t
  """

  p_in_bar = specs.outlet_P * 1.0e-5;
  return 48.14 * p_in_bar / (t + 0.4446 * p_in_bar / math.pow(t, 0.2));

# -------------- Unit Conversions: OpenMC requires cm -----------
m = 100.0

# estimate the outlet temperature using bulk energy conservation for steady state
coolant_outlet_temp = specs.power / specs.mdot / specs.fluid_Cp + specs.inlet_T

# geometry parameters
coolant_channel_diam = specs.channel_diameter * m
reactor_bottom = 0.0
reactor_height = specs.height * m
reactor_top = reactor_bottom + reactor_height
bundle_pitch = specs.bundle_flat_to_flat * m + specs.bundle_gap_width * m
cell_pitch = specs.fuel_to_coolant_distance * m
fuel_channel_diam = specs.compact_diameter * m
top_reflector_height = specs.top_reflector_thickness * m
bottom_reflector_height = specs.bottom_reflector_thickness * m

def assembly(n_ax_zones, n_inactive, n_active, add_entropy_mesh=False):
    axial_section_height = reactor_height / n_ax_zones

    # superimposed search lattice
    triso_lattice_shape = (4, 4, int(axial_section_height / 0.5))

    model = openmc.model.Model()

    m_b4c = openmc.Material(name='B4C')
    enrichment_10 = specs.B10_enrichment
    mass_10 = openmc.data.atomic_mass('B10')
    mass_11 = openmc.data.atomic_mass('B11')

    # number of atoms in one gram of boron mixture
    n_10 = enrichment_10 / mass_10
    n_11 = (1.0 - enrichment_10) / mass_11
    total_n = n_10 + n_11
    grams_10 = n_10 / total_n
    grams_11 = n_11 / total_n

    # now, figure out how much carbon needs to be in the poison to get
    # an overall specified B10 weight percent
    total_b10_weight_percent = specs.total_B10_wt_percent
    total_mass = grams_10 / total_b10_weight_percent
    carbon_mass = total_mass - grams_10 - grams_11

    m_b4c.add_nuclide('B10', grams_10 / total_mass, 'wo')
    m_b4c.add_nuclide('B11', grams_11 / total_mass, 'wo')
    m_b4c.add_element('C', carbon_mass / total_mass, 'wo')
    m_b4c.set_density('kg/m3', specs.B4C_density)

    # reflector is 40 percent helium by volume (arbitrary assumption), with helium
    # at the inlet conditions
    rho = coolant_density(specs.inlet_T)
    matrix_density = 1700
    reflector_porosity = 0.40
    n_helium = reflector_porosity * rho / 4.002602
    n_carbon = (1.0 - reflector_porosity) * matrix_density / 12.0107
    combined_density = rho * reflector_porosity + matrix_density * (1.0 - reflector_porosity)
    m_reflector = openmc.Material(name='reflector')
    m_reflector.add_element('He', n_helium / (n_helium + n_carbon))
    m_reflector.add_element('C', n_carbon / (n_helium + n_carbon))
    m_reflector.set_density('kg/m3', combined_density)

    # TRISO particle
    radius_pyc_outer   = specs.oPyC_radius * m

    s_fuel             = openmc.Sphere(r=specs.kernel_radius*m)
    s_c_buffer         = openmc.Sphere(r=specs.buffer_radius*m)
    s_pyc_inner        = openmc.Sphere(r=specs.iPyC_radius*m)
    s_sic              = openmc.Sphere(r=specs.SiC_radius*m)
    s_pyc_outer        = openmc.Sphere(r=radius_pyc_outer)
    c_triso_fuel       = openmc.Cell(name='c_triso_fuel'     , fill=mats.m_fuel,              region=-s_fuel)
    c_triso_c_buffer   = openmc.Cell(name='c_triso_c_buffer' , fill=mats.m_graphite_c_buffer, region=+s_fuel      & -s_c_buffer)
    c_triso_pyc_inner  = openmc.Cell(name='c_triso_pyc_inner', fill=mats.m_graphite_pyc,      region=+s_c_buffer  & -s_pyc_inner)
    c_triso_sic        = openmc.Cell(name='c_triso_sic'      , fill=mats.m_sic,               region=+s_pyc_inner & -s_sic)
    c_triso_pyc_outer  = openmc.Cell(name='c_triso_pyc_outer', fill=mats.m_graphite_pyc,      region=+s_sic       & -s_pyc_outer)
    c_triso_matrix     = openmc.Cell(name='c_triso_matrix'   , fill=mats.m_graphite_matrix,   region=+s_pyc_outer)
    u_triso            = openmc.Universe(cells=[c_triso_fuel, c_triso_c_buffer, c_triso_pyc_inner, c_triso_sic, c_triso_pyc_outer, c_triso_matrix])

    # Channel surfaces
    fuel_cyl = openmc.ZCylinder(r=0.5 * fuel_channel_diam)
    coolant_cyl = openmc.ZCylinder(r=0.5 * coolant_channel_diam)
    poison_cyl = openmc.ZCylinder(r=0.5 * fuel_channel_diam)
    graphite_cyl = openmc.ZCylinder(r=0.5 * fuel_channel_diam)

    # create a TRISO lattice for one axial section (to be used in the rest of the axial zones)
    # center the TRISO region on the origin so it fills lattice cells appropriately
    min_z = openmc.ZPlane(z0=-0.5 * axial_section_height)
    max_z = openmc.ZPlane(z0=0.5 * axial_section_height)

    # region in which TRISOs are generated
    r_triso = -fuel_cyl & +min_z & -max_z

    rand_spheres = openmc.model.pack_spheres(radius=radius_pyc_outer, region=r_triso, pf=specs.triso_pf)
    random_trisos = [openmc.model.TRISO(radius_pyc_outer, u_triso, i) for i in rand_spheres]

    llc, urc = r_triso.bounding_box
    pitch = (urc - llc) / triso_lattice_shape

    # insert TRISOs into a lattice to accelerate point location queries
    triso_lattice = openmc.model.create_triso_lattice(random_trisos, llc, pitch, triso_lattice_shape, mats.m_graphite_matrix)

    axial_coords = np.linspace(reactor_bottom, reactor_top, n_ax_zones + 1)
    lattice_univs = []
    bundle_univs = []

    m_colors = {}

    for z_min, z_max in zip(axial_coords[0:-1], axial_coords[1:]):
        # use the middle of the axial section to compute the temperature and density
        ax_pos = 0.5 * (z_min + z_max)
        T = coolant_temp(specs.inlet_T, coolant_outlet_temp, reactor_height, ax_pos)

        # create solid cells, which don't require us to clone materials in order to set temperatures
        fuel_ch_cell = openmc.Cell(region=-fuel_cyl, fill=triso_lattice)
        fuel_ch_cell.temperature = T

        fuel_ch_matrix_cell = openmc.Cell(region=+fuel_cyl, fill=mats.m_graphite_matrix)
        fuel_ch_matrix_cell.temperature = T

        poison_cell = openmc.Cell(region=-poison_cyl, fill=m_b4c)
        poison_cell.temperature = T

        poison_matrix_cell = openmc.Cell(region=+poison_cyl, fill=mats.m_graphite_matrix)
        poison_matrix_cell.temperature = T

        graphite_cell = openmc.Cell(fill=mats.m_graphite_matrix)
        graphite_cell.temperature = T

        coolant_matrix_cell = openmc.Cell(region=+coolant_cyl, fill=mats.m_graphite_matrix)
        coolant_matrix_cell.temperature = T

        # create fluid cells, which require us to clone the material in order to be able to
        # set unique densities
        coolant_cell = openmc.Cell(region=-coolant_cyl, fill=mats.m_coolant)
        coolant_cell.fill = [mats.m_coolant.clone() for i in range(specs.n_coolant_channels_per_block)]

        for mat in range(len(coolant_cell.fill)):
          m_colors[coolant_cell.fill[mat]] = 'red'

        # Define a universe for each type of solid-only pin (fuel, poison, and graphite)
        f = openmc.Universe(cells=[fuel_ch_cell, fuel_ch_matrix_cell])
        p = openmc.Universe(cells=[poison_cell, poison_matrix_cell])
        c = openmc.Universe(cells=[coolant_cell, coolant_matrix_cell])
        g = openmc.Universe(cells=[graphite_cell])

        d = [f] * 2

        ring2 = ([f] + [c]) * 6
        ring3 = ([c] + d) * 6
        ring4 = (d + [c] + [f]) * 6
        ring5 = ([f] + [c] + d + [c]) * 6
        ring6 = ([c] + d + [c] + d) * 6
        ring7 = (d + [c] + d + [c] + [f]) * 6
        ring8 = ([f] + [c] + d + [c] + d + [c]) * 6
        ring9 = ([c] + d + [c] + d + [c] + d) * 6
        ring10 = ([p] + [f] + [c] + d + [c] + d + [c] + [f]) * 6
        ring11 = [g] * 66

        # inner two rings where there aren't any fuel/compact/poison pins
        ring1 = [g] * 6
        ring0 = [g]

        lattice_univs.append([ring11, ring10, ring9, ring8, ring7, ring6, ring5, ring4, ring3, ring2, ring1, ring0])

    # create a hexagonal lattice used in each axial zone to represent the cell
    hex_lattice = openmc.HexLattice(name="Bundle cell lattice")
    hex_lattice.orientation = 'x'
    hex_lattice.center = (0.0, 0.0, 0.5 * (reactor_bottom + reactor_top))
    hex_lattice.pitch = (cell_pitch, axial_section_height)
    hex_lattice.universes = lattice_univs

    hexagon_volume = reactor_height * math.sqrt(3) / 2.0 * bundle_pitch**2
    coolant_channel_volume = math.pi * coolant_channel_diam**2 / 4.0 * reactor_height

    graphite_outer_cell = openmc.Cell(fill=mats.m_graphite_matrix)
    graphite_outer_cell.temperature = T
    inf_graphite_univ = openmc.Universe(cells=[graphite_outer_cell])
    hex_lattice.outer = inf_graphite_univ

    # create additional axial regions
    axial_planes = [openmc.ZPlane(z0=coord) for coord in axial_coords]

    # axial planes
    min_z = axial_planes[0]
    max_z = axial_planes[-1]

    # fill the unit cell with the hex lattice
    hex_prism = openmc.hexagonal_prism(bundle_pitch / math.sqrt(3.0), 'x', boundary_type='periodic')
    outer_cell = openmc.Cell(region=hex_prism & +min_z & -max_z, fill=hex_lattice)

    # add the top and bottom reflector
    top_refl_z = reactor_height + top_reflector_height
    top_refl = openmc.ZPlane(z0=top_refl_z, boundary_type='vacuum')
    bottom_refl_z = -bottom_reflector_height
    bottom_refl = openmc.ZPlane(z0=bottom_refl_z, boundary_type='vacuum')

    top_refl_cell = openmc.Cell(region=hex_prism & +max_z & -top_refl, fill=m_reflector)
    bottom_refl_cell = openmc.Cell(region=hex_prism & -min_z & +bottom_refl, fill=m_reflector)
    top_refl_cell.temperature = specs.inlet_T
    bottom_refl_cell.temperature = coolant_outlet_temp

    model.geometry = openmc.Geometry([outer_cell, top_refl_cell, bottom_refl_cell])

    ### Settings ###
    settings = openmc.Settings()

    settings.particles = 10000
    settings.inactive = n_inactive
    settings.batches = settings.inactive + n_active
    settings.temperature['method'] = 'interpolation'
    settings.temperature['range'] = (294.0, 1500.0)

    l = bundle_pitch / math.sqrt(3.0)
    lower_left = (-l, -l, reactor_bottom)
    upper_right = (l, l, reactor_top)
    source_dist = openmc.stats.Box(lower_left, upper_right, only_fissionable=True)
    source = openmc.Source(space=source_dist)
    settings.source = source

    if (add_entropy_mesh):
        entropy_mesh = openmc.RegularMesh()
        entropy_mesh.lower_left = lower_left
        entropy_mesh.upper_right = upper_right
        entropy_mesh.dimension = (6, 6, 20)
        settings.entropy_mesh = entropy_mesh

    vol_calc = openmc.VolumeCalculation([outer_cell],
                                        100_000_000, lower_left, upper_right)
    settings.volume_calculations = [vol_calc]

    model.settings = settings

    m_colors[mats.m_fuel] = 'palegreen'
    m_colors[mats.m_graphite_c_buffer] = 'sandybrown'
    m_colors[mats.m_graphite_pyc] = 'orange'
    m_colors[mats.m_sic] = 'yellow'
    m_colors[mats.m_graphite_matrix] = 'darkblue'
    m_colors[m_b4c] = 'lightskyblue'

    bundle_p_rounded = int(bundle_pitch)

    plot1          = openmc.Plot()
    plot1.filename = 'plot1'
    plot1.width    = (2 * bundle_pitch, 2 * axial_section_height)
    plot1.basis    = 'xz'
    plot1.origin   = (0.0, 0.0, reactor_height / 2.0)
    plot1.pixels   = (100 * 2 * bundle_p_rounded, int(100 * 3 * axial_section_height))
    plot1.color_by = 'cell'

    plot2          = openmc.Plot()
    plot2.filename = 'plot2'
    plot2.width    = (1.5 * bundle_pitch, 1.5 * bundle_pitch)
    plot2.basis    = 'xy'
    plot2.origin   = (0.0, 0.0, axial_section_height / 4.0)
    plot2.pixels   = (500 * bundle_p_rounded, 500 * bundle_p_rounded)
    plot2.color_by = 'material'
    plot2.colors   = m_colors

    plot3          = openmc.Plot()
    plot3.filename = 'plot3'
    plot3.width    = plot2.width
    plot3.basis    = plot2.basis
    plot3.origin   = plot2.origin
    plot3.pixels   = plot2.pixels
    plot3.color_by = 'cell'

    model.plots = openmc.Plots([plot1, plot2, plot3])

    return model


def main():

    ap = ArgumentParser()
    ap.add_argument('-n', dest='n_axial', type=int, default=50,
                    help='Number of axial cell divisions')
    ap.add_argument('-s', '--entropy', action='store_true',
                    help='Whether to add a Shannon entropy mesh')
    ap.add_argument('-i', dest='n_inactive', type=int, default=500,
                    help='Number of inactive cycles')
    ap.add_argument('-a', dest='n_active', type=int, default=2000,
                    help='Number of active cycles')

    args = ap.parse_args()

    model = assembly(args.n_axial, args.n_inactive, args.n_active, args.entropy)
    model.export_to_xml()

if __name__ == "__main__":
    main()

(htgr/assembly/solid.i)

## HTGR assembly heat conduction input
## Application: MOOSE heat conduction module
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

# copy-pasta from common_input.i
inlet_T = 598.0                          # inlet fluid temperature (K)
buffer_k = 0.5                           # buffer thermal conductivity (W/m/K)
PyC_k = 4.0                              # PyC thermal conductivity (W/m/K)
SiC_k = 13.9                             # SiC thermal conductivity (W/m/K)
kernel_k = 3.5                           # fissil kernel thermal conductivity (W/m/K)
matrix_k = 15.0                          # graphite matrix thermal conductivity (W/m/K)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
triso_pf = 0.15                          # TRISO packing fraction (%)
kernel_radius = 214.85e-6                # fissile kernel outer radius (m)
buffer_radius = 314.85e-6                # buffer outer radius (m)
iPyC_radius = 354.85e-6                  # inner PyC outer radius (m)
SiC_radius = 389.85e-6                   # SiC outer radius (m)
oPyC_radius = 429.85e-6                  # outer PyC outer radius (m)

# compute the volume fraction of each TRISO layer in a TRISO particle
# for use in computing average thermophysical properties
kernel_fraction = ${fparse kernel_radius^3 / oPyC_radius^3}
buffer_fraction = ${fparse (buffer_radius^3 - kernel_radius^3) / oPyC_radius^3}
ipyc_fraction = ${fparse (iPyC_radius^3 - buffer_radius^3) / oPyC_radius^3}
sic_fraction = ${fparse (SiC_radius^3 - iPyC_radius^3) / oPyC_radius^3}
opyc_fraction = ${fparse (oPyC_radius^3 - SiC_radius^3) / oPyC_radius^3}

[Mesh]
  type = FileMesh
  file = meshes/solid_in.e
[]

[Variables]
  [T]
    initial_condition = ${inlet_T}
  []
[]

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

[BCs]
  [pin_outer]
    type = MatchedValueBC
    variable = T
    v = thm_temp
    boundary = 'fluid_solid_interface'
  []
[]

[Functions]
  [k_graphite]
    type = ParsedFunction
    value = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    value = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    value = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    vars = 'k_TRISO k_graphite'
    vals = 'k_TRISO k_graphite'
  []
  [k_b4c]
    type = ParsedFunction
    value = '5.096154e-6 * t - 1.952360e-2 * t + 2.558435e1'
  []
[]

[Materials]
  [graphite]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_graphite
    temp = T
    block = 'graphite'
  []
  [compacts]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_compacts
    temp = T
    block = 'compacts'
  []
  [poison]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_b4c
    temp = T
    block = 'poison'
  []
[]

[Postprocessors]
  [flux_integral] # evaluate the total heat flux for normalization
    type = SideDiffusiveFluxIntegral
    diffusivity = thermal_conductivity
    variable = T
    boundary = 'fluid_solid_interface'
  []
  [max_fuel_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'compacts'
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [power] # evaluate the total power for normalization
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = 'compacts'
    execute_on = 'transfer'
  []
[]

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

[AuxKernels]
  [flux]
    type = DiffusionFluxAux
    diffusion_variable = T
    component = normal
    diffusivity = thermal_conductivity
    variable = flux
    boundary = 'fluid_solid_interface'
  []
[]

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

[UserObjects]
  [average_fuel_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_block_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'graphite'
  []
[]

[VectorPostprocessors]
  [fuel_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fuel_axial
  []
  [block_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_block_axial
  []
[]

[Outputs]
  exodus = true
  csv = true
  print_linear_residuals = false
[]

(htgr/assembly/solid.i)

## HTGR assembly heat conduction input
## Application: MOOSE heat conduction module
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

# copy-pasta from common_input.i
inlet_T = 598.0                          # inlet fluid temperature (K)
buffer_k = 0.5                           # buffer thermal conductivity (W/m/K)
PyC_k = 4.0                              # PyC thermal conductivity (W/m/K)
SiC_k = 13.9                             # SiC thermal conductivity (W/m/K)
kernel_k = 3.5                           # fissil kernel thermal conductivity (W/m/K)
matrix_k = 15.0                          # graphite matrix thermal conductivity (W/m/K)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
triso_pf = 0.15                          # TRISO packing fraction (%)
kernel_radius = 214.85e-6                # fissile kernel outer radius (m)
buffer_radius = 314.85e-6                # buffer outer radius (m)
iPyC_radius = 354.85e-6                  # inner PyC outer radius (m)
SiC_radius = 389.85e-6                   # SiC outer radius (m)
oPyC_radius = 429.85e-6                  # outer PyC outer radius (m)

# compute the volume fraction of each TRISO layer in a TRISO particle
# for use in computing average thermophysical properties
kernel_fraction = ${fparse kernel_radius^3 / oPyC_radius^3}
buffer_fraction = ${fparse (buffer_radius^3 - kernel_radius^3) / oPyC_radius^3}
ipyc_fraction = ${fparse (iPyC_radius^3 - buffer_radius^3) / oPyC_radius^3}
sic_fraction = ${fparse (SiC_radius^3 - iPyC_radius^3) / oPyC_radius^3}
opyc_fraction = ${fparse (oPyC_radius^3 - SiC_radius^3) / oPyC_radius^3}

[Mesh]
  type = FileMesh
  file = meshes/solid_in.e
[]

[Variables]
  [T]
    initial_condition = ${inlet_T}
  []
[]

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

[BCs]
  [pin_outer]
    type = MatchedValueBC
    variable = T
    v = thm_temp
    boundary = 'fluid_solid_interface'
  []
[]

[Functions]
  [k_graphite]
    type = ParsedFunction
    value = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    value = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    value = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    vars = 'k_TRISO k_graphite'
    vals = 'k_TRISO k_graphite'
  []
  [k_b4c]
    type = ParsedFunction
    value = '5.096154e-6 * t - 1.952360e-2 * t + 2.558435e1'
  []
[]

[Materials]
  [graphite]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_graphite
    temp = T
    block = 'graphite'
  []
  [compacts]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_compacts
    temp = T
    block = 'compacts'
  []
  [poison]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_b4c
    temp = T
    block = 'poison'
  []
[]

[Postprocessors]
  [flux_integral] # evaluate the total heat flux for normalization
    type = SideDiffusiveFluxIntegral
    diffusivity = thermal_conductivity
    variable = T
    boundary = 'fluid_solid_interface'
  []
  [max_fuel_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'compacts'
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [power] # evaluate the total power for normalization
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = 'compacts'
    execute_on = 'transfer'
  []
[]

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

[AuxKernels]
  [flux]
    type = DiffusionFluxAux
    diffusion_variable = T
    component = normal
    diffusivity = thermal_conductivity
    variable = flux
    boundary = 'fluid_solid_interface'
  []
[]

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

[UserObjects]
  [average_fuel_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_block_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'graphite'
  []
[]

[VectorPostprocessors]
  [fuel_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fuel_axial
  []
  [block_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_block_axial
  []
[]

[Outputs]
  exodus = true
  csv = true
  print_linear_residuals = false
[]

(htgr/assembly/solid.i)

## HTGR assembly heat conduction input
## Application: MOOSE heat conduction module
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

# copy-pasta from common_input.i
inlet_T = 598.0                          # inlet fluid temperature (K)
buffer_k = 0.5                           # buffer thermal conductivity (W/m/K)
PyC_k = 4.0                              # PyC thermal conductivity (W/m/K)
SiC_k = 13.9                             # SiC thermal conductivity (W/m/K)
kernel_k = 3.5                           # fissil kernel thermal conductivity (W/m/K)
matrix_k = 15.0                          # graphite matrix thermal conductivity (W/m/K)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
triso_pf = 0.15                          # TRISO packing fraction (%)
kernel_radius = 214.85e-6                # fissile kernel outer radius (m)
buffer_radius = 314.85e-6                # buffer outer radius (m)
iPyC_radius = 354.85e-6                  # inner PyC outer radius (m)
SiC_radius = 389.85e-6                   # SiC outer radius (m)
oPyC_radius = 429.85e-6                  # outer PyC outer radius (m)

# compute the volume fraction of each TRISO layer in a TRISO particle
# for use in computing average thermophysical properties
kernel_fraction = ${fparse kernel_radius^3 / oPyC_radius^3}
buffer_fraction = ${fparse (buffer_radius^3 - kernel_radius^3) / oPyC_radius^3}
ipyc_fraction = ${fparse (iPyC_radius^3 - buffer_radius^3) / oPyC_radius^3}
sic_fraction = ${fparse (SiC_radius^3 - iPyC_radius^3) / oPyC_radius^3}
opyc_fraction = ${fparse (oPyC_radius^3 - SiC_radius^3) / oPyC_radius^3}

[Mesh]
  type = FileMesh
  file = meshes/solid_in.e
[]

[Variables]
  [T]
    initial_condition = ${inlet_T}
  []
[]

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

[BCs]
  [pin_outer]
    type = MatchedValueBC
    variable = T
    v = thm_temp
    boundary = 'fluid_solid_interface'
  []
[]

[Functions]
  [k_graphite]
    type = ParsedFunction
    value = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    value = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    value = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    vars = 'k_TRISO k_graphite'
    vals = 'k_TRISO k_graphite'
  []
  [k_b4c]
    type = ParsedFunction
    value = '5.096154e-6 * t - 1.952360e-2 * t + 2.558435e1'
  []
[]

[Materials]
  [graphite]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_graphite
    temp = T
    block = 'graphite'
  []
  [compacts]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_compacts
    temp = T
    block = 'compacts'
  []
  [poison]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_b4c
    temp = T
    block = 'poison'
  []
[]

[Postprocessors]
  [flux_integral] # evaluate the total heat flux for normalization
    type = SideDiffusiveFluxIntegral
    diffusivity = thermal_conductivity
    variable = T
    boundary = 'fluid_solid_interface'
  []
  [max_fuel_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'compacts'
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [power] # evaluate the total power for normalization
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = 'compacts'
    execute_on = 'transfer'
  []
[]

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

[AuxKernels]
  [flux]
    type = DiffusionFluxAux
    diffusion_variable = T
    component = normal
    diffusivity = thermal_conductivity
    variable = flux
    boundary = 'fluid_solid_interface'
  []
[]

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

[UserObjects]
  [average_fuel_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_block_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'graphite'
  []
[]

[VectorPostprocessors]
  [fuel_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fuel_axial
  []
  [block_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_block_axial
  []
[]

[Outputs]
  exodus = true
  csv = true
  print_linear_residuals = false
[]

(htgr/assembly/solid.i)

## HTGR assembly heat conduction input
## Application: MOOSE heat conduction module
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

# copy-pasta from common_input.i
inlet_T = 598.0                          # inlet fluid temperature (K)
buffer_k = 0.5                           # buffer thermal conductivity (W/m/K)
PyC_k = 4.0                              # PyC thermal conductivity (W/m/K)
SiC_k = 13.9                             # SiC thermal conductivity (W/m/K)
kernel_k = 3.5                           # fissil kernel thermal conductivity (W/m/K)
matrix_k = 15.0                          # graphite matrix thermal conductivity (W/m/K)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
triso_pf = 0.15                          # TRISO packing fraction (%)
kernel_radius = 214.85e-6                # fissile kernel outer radius (m)
buffer_radius = 314.85e-6                # buffer outer radius (m)
iPyC_radius = 354.85e-6                  # inner PyC outer radius (m)
SiC_radius = 389.85e-6                   # SiC outer radius (m)
oPyC_radius = 429.85e-6                  # outer PyC outer radius (m)

# compute the volume fraction of each TRISO layer in a TRISO particle
# for use in computing average thermophysical properties
kernel_fraction = ${fparse kernel_radius^3 / oPyC_radius^3}
buffer_fraction = ${fparse (buffer_radius^3 - kernel_radius^3) / oPyC_radius^3}
ipyc_fraction = ${fparse (iPyC_radius^3 - buffer_radius^3) / oPyC_radius^3}
sic_fraction = ${fparse (SiC_radius^3 - iPyC_radius^3) / oPyC_radius^3}
opyc_fraction = ${fparse (oPyC_radius^3 - SiC_radius^3) / oPyC_radius^3}

[Mesh]
  type = FileMesh
  file = meshes/solid_in.e
[]

[Variables]
  [T]
    initial_condition = ${inlet_T}
  []
[]

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

[BCs]
  [pin_outer]
    type = MatchedValueBC
    variable = T
    v = thm_temp
    boundary = 'fluid_solid_interface'
  []
[]

[Functions]
  [k_graphite]
    type = ParsedFunction
    value = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    value = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    value = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    vars = 'k_TRISO k_graphite'
    vals = 'k_TRISO k_graphite'
  []
  [k_b4c]
    type = ParsedFunction
    value = '5.096154e-6 * t - 1.952360e-2 * t + 2.558435e1'
  []
[]

[Materials]
  [graphite]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_graphite
    temp = T
    block = 'graphite'
  []
  [compacts]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_compacts
    temp = T
    block = 'compacts'
  []
  [poison]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_b4c
    temp = T
    block = 'poison'
  []
[]

[Postprocessors]
  [flux_integral] # evaluate the total heat flux for normalization
    type = SideDiffusiveFluxIntegral
    diffusivity = thermal_conductivity
    variable = T
    boundary = 'fluid_solid_interface'
  []
  [max_fuel_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'compacts'
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [power] # evaluate the total power for normalization
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = 'compacts'
    execute_on = 'transfer'
  []
[]

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

[AuxKernels]
  [flux]
    type = DiffusionFluxAux
    diffusion_variable = T
    component = normal
    diffusivity = thermal_conductivity
    variable = flux
    boundary = 'fluid_solid_interface'
  []
[]

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

[UserObjects]
  [average_fuel_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_block_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'graphite'
  []
[]

[VectorPostprocessors]
  [fuel_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fuel_axial
  []
  [block_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_block_axial
  []
[]

[Outputs]
  exodus = true
  csv = true
  print_linear_residuals = false
[]

(htgr/assembly/solid.i)

## HTGR assembly heat conduction input
## Application: MOOSE heat conduction module
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

# copy-pasta from common_input.i
inlet_T = 598.0                          # inlet fluid temperature (K)
buffer_k = 0.5                           # buffer thermal conductivity (W/m/K)
PyC_k = 4.0                              # PyC thermal conductivity (W/m/K)
SiC_k = 13.9                             # SiC thermal conductivity (W/m/K)
kernel_k = 3.5                           # fissil kernel thermal conductivity (W/m/K)
matrix_k = 15.0                          # graphite matrix thermal conductivity (W/m/K)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
triso_pf = 0.15                          # TRISO packing fraction (%)
kernel_radius = 214.85e-6                # fissile kernel outer radius (m)
buffer_radius = 314.85e-6                # buffer outer radius (m)
iPyC_radius = 354.85e-6                  # inner PyC outer radius (m)
SiC_radius = 389.85e-6                   # SiC outer radius (m)
oPyC_radius = 429.85e-6                  # outer PyC outer radius (m)

# compute the volume fraction of each TRISO layer in a TRISO particle
# for use in computing average thermophysical properties
kernel_fraction = ${fparse kernel_radius^3 / oPyC_radius^3}
buffer_fraction = ${fparse (buffer_radius^3 - kernel_radius^3) / oPyC_radius^3}
ipyc_fraction = ${fparse (iPyC_radius^3 - buffer_radius^3) / oPyC_radius^3}
sic_fraction = ${fparse (SiC_radius^3 - iPyC_radius^3) / oPyC_radius^3}
opyc_fraction = ${fparse (oPyC_radius^3 - SiC_radius^3) / oPyC_radius^3}

[Mesh]
  type = FileMesh
  file = meshes/solid_in.e
[]

[Variables]
  [T]
    initial_condition = ${inlet_T}
  []
[]

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

[BCs]
  [pin_outer]
    type = MatchedValueBC
    variable = T
    v = thm_temp
    boundary = 'fluid_solid_interface'
  []
[]

[Functions]
  [k_graphite]
    type = ParsedFunction
    value = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    value = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    value = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    vars = 'k_TRISO k_graphite'
    vals = 'k_TRISO k_graphite'
  []
  [k_b4c]
    type = ParsedFunction
    value = '5.096154e-6 * t - 1.952360e-2 * t + 2.558435e1'
  []
[]

[Materials]
  [graphite]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_graphite
    temp = T
    block = 'graphite'
  []
  [compacts]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_compacts
    temp = T
    block = 'compacts'
  []
  [poison]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_b4c
    temp = T
    block = 'poison'
  []
[]

[Postprocessors]
  [flux_integral] # evaluate the total heat flux for normalization
    type = SideDiffusiveFluxIntegral
    diffusivity = thermal_conductivity
    variable = T
    boundary = 'fluid_solid_interface'
  []
  [max_fuel_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'compacts'
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [power] # evaluate the total power for normalization
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = 'compacts'
    execute_on = 'transfer'
  []
[]

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

[AuxKernels]
  [flux]
    type = DiffusionFluxAux
    diffusion_variable = T
    component = normal
    diffusivity = thermal_conductivity
    variable = flux
    boundary = 'fluid_solid_interface'
  []
[]

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

[UserObjects]
  [average_fuel_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_block_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'graphite'
  []
[]

[VectorPostprocessors]
  [fuel_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fuel_axial
  []
  [block_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_block_axial
  []
[]

[Outputs]
  exodus = true
  csv = true
  print_linear_residuals = false
[]

(htgr/assembly/solid.i)

## HTGR assembly heat conduction input
## Application: MOOSE heat conduction module
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

# copy-pasta from common_input.i
inlet_T = 598.0                          # inlet fluid temperature (K)
buffer_k = 0.5                           # buffer thermal conductivity (W/m/K)
PyC_k = 4.0                              # PyC thermal conductivity (W/m/K)
SiC_k = 13.9                             # SiC thermal conductivity (W/m/K)
kernel_k = 3.5                           # fissil kernel thermal conductivity (W/m/K)
matrix_k = 15.0                          # graphite matrix thermal conductivity (W/m/K)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
triso_pf = 0.15                          # TRISO packing fraction (%)
kernel_radius = 214.85e-6                # fissile kernel outer radius (m)
buffer_radius = 314.85e-6                # buffer outer radius (m)
iPyC_radius = 354.85e-6                  # inner PyC outer radius (m)
SiC_radius = 389.85e-6                   # SiC outer radius (m)
oPyC_radius = 429.85e-6                  # outer PyC outer radius (m)

# compute the volume fraction of each TRISO layer in a TRISO particle
# for use in computing average thermophysical properties
kernel_fraction = ${fparse kernel_radius^3 / oPyC_radius^3}
buffer_fraction = ${fparse (buffer_radius^3 - kernel_radius^3) / oPyC_radius^3}
ipyc_fraction = ${fparse (iPyC_radius^3 - buffer_radius^3) / oPyC_radius^3}
sic_fraction = ${fparse (SiC_radius^3 - iPyC_radius^3) / oPyC_radius^3}
opyc_fraction = ${fparse (oPyC_radius^3 - SiC_radius^3) / oPyC_radius^3}

[Mesh]
  type = FileMesh
  file = meshes/solid_in.e
[]

[Variables]
  [T]
    initial_condition = ${inlet_T}
  []
[]

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

[BCs]
  [pin_outer]
    type = MatchedValueBC
    variable = T
    v = thm_temp
    boundary = 'fluid_solid_interface'
  []
[]

[Functions]
  [k_graphite]
    type = ParsedFunction
    value = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    value = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    value = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    vars = 'k_TRISO k_graphite'
    vals = 'k_TRISO k_graphite'
  []
  [k_b4c]
    type = ParsedFunction
    value = '5.096154e-6 * t - 1.952360e-2 * t + 2.558435e1'
  []
[]

[Materials]
  [graphite]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_graphite
    temp = T
    block = 'graphite'
  []
  [compacts]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_compacts
    temp = T
    block = 'compacts'
  []
  [poison]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_b4c
    temp = T
    block = 'poison'
  []
[]

[Postprocessors]
  [flux_integral] # evaluate the total heat flux for normalization
    type = SideDiffusiveFluxIntegral
    diffusivity = thermal_conductivity
    variable = T
    boundary = 'fluid_solid_interface'
  []
  [max_fuel_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'compacts'
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [power] # evaluate the total power for normalization
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = 'compacts'
    execute_on = 'transfer'
  []
[]

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

[AuxKernels]
  [flux]
    type = DiffusionFluxAux
    diffusion_variable = T
    component = normal
    diffusivity = thermal_conductivity
    variable = flux
    boundary = 'fluid_solid_interface'
  []
[]

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

[UserObjects]
  [average_fuel_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_block_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'graphite'
  []
[]

[VectorPostprocessors]
  [fuel_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fuel_axial
  []
  [block_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_block_axial
  []
[]

[Outputs]
  exodus = true
  csv = true
  print_linear_residuals = false
[]

(htgr/assembly/solid.i)

## HTGR assembly heat conduction input
## Application: MOOSE heat conduction module
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

# copy-pasta from common_input.i
inlet_T = 598.0                          # inlet fluid temperature (K)
buffer_k = 0.5                           # buffer thermal conductivity (W/m/K)
PyC_k = 4.0                              # PyC thermal conductivity (W/m/K)
SiC_k = 13.9                             # SiC thermal conductivity (W/m/K)
kernel_k = 3.5                           # fissil kernel thermal conductivity (W/m/K)
matrix_k = 15.0                          # graphite matrix thermal conductivity (W/m/K)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
triso_pf = 0.15                          # TRISO packing fraction (%)
kernel_radius = 214.85e-6                # fissile kernel outer radius (m)
buffer_radius = 314.85e-6                # buffer outer radius (m)
iPyC_radius = 354.85e-6                  # inner PyC outer radius (m)
SiC_radius = 389.85e-6                   # SiC outer radius (m)
oPyC_radius = 429.85e-6                  # outer PyC outer radius (m)

# compute the volume fraction of each TRISO layer in a TRISO particle
# for use in computing average thermophysical properties
kernel_fraction = ${fparse kernel_radius^3 / oPyC_radius^3}
buffer_fraction = ${fparse (buffer_radius^3 - kernel_radius^3) / oPyC_radius^3}
ipyc_fraction = ${fparse (iPyC_radius^3 - buffer_radius^3) / oPyC_radius^3}
sic_fraction = ${fparse (SiC_radius^3 - iPyC_radius^3) / oPyC_radius^3}
opyc_fraction = ${fparse (oPyC_radius^3 - SiC_radius^3) / oPyC_radius^3}

[Mesh]
  type = FileMesh
  file = meshes/solid_in.e
[]

[Variables]
  [T]
    initial_condition = ${inlet_T}
  []
[]

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

[BCs]
  [pin_outer]
    type = MatchedValueBC
    variable = T
    v = thm_temp
    boundary = 'fluid_solid_interface'
  []
[]

[Functions]
  [k_graphite]
    type = ParsedFunction
    value = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    value = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    value = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    vars = 'k_TRISO k_graphite'
    vals = 'k_TRISO k_graphite'
  []
  [k_b4c]
    type = ParsedFunction
    value = '5.096154e-6 * t - 1.952360e-2 * t + 2.558435e1'
  []
[]

[Materials]
  [graphite]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_graphite
    temp = T
    block = 'graphite'
  []
  [compacts]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_compacts
    temp = T
    block = 'compacts'
  []
  [poison]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_b4c
    temp = T
    block = 'poison'
  []
[]

[Postprocessors]
  [flux_integral] # evaluate the total heat flux for normalization
    type = SideDiffusiveFluxIntegral
    diffusivity = thermal_conductivity
    variable = T
    boundary = 'fluid_solid_interface'
  []
  [max_fuel_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'compacts'
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [power] # evaluate the total power for normalization
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = 'compacts'
    execute_on = 'transfer'
  []
[]

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

[AuxKernels]
  [flux]
    type = DiffusionFluxAux
    diffusion_variable = T
    component = normal
    diffusivity = thermal_conductivity
    variable = flux
    boundary = 'fluid_solid_interface'
  []
[]

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

[UserObjects]
  [average_fuel_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_block_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'graphite'
  []
[]

[VectorPostprocessors]
  [fuel_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fuel_axial
  []
  [block_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_block_axial
  []
[]

[Outputs]
  exodus = true
  csv = true
  print_linear_residuals = false
[]

(htgr/assembly/thm.i)

## HTGR assembly thermal hydraulics simulation
## Application: MOOSE thermal hydraulics module
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

# copy-pasta from common.i
inlet_T = 598.0 # inlet fluid temperature (K)
mdot = '${fparse 117.3 / 12 / 108}' # fluid mass flowrate (kg/s)
outlet_P = 7.1e6 # fluid outlet pressure (Pa)
channel_diameter = 0.016 # diameter of the coolant channels (m)
height = 6.343 # height of the assembly (m)

num_layers_for_THM = 50 # number of elements in the THM model

[GlobalParams]
  initial_p = ${outlet_P}
  initial_T = ${inlet_T}
  initial_vel = '${fparse mdot / outlet_P / 8.3144598 * 4.0e-3 / inlet_T / (pi * channel_diameter * channel_diameter / 4.0)}'

  rdg_slope_reconstruction = full
  closures = none
  fp = helium
[]

[Closures]
  [none]
    type = Closures1PhaseNone
  []
[]

[FluidProperties]
  [helium]
    type = IdealGasFluidProperties
    molar_mass = 4e-3
    gamma = 1.668282 # should correspond to  Cp = 5189 J/kg/K
    k = 0.2556
    mu = 3.22639e-5
  []
[]

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

[AuxKernels]
  [Tw_aux]
    type = ADMaterialRealAux
    block = channel
    variable = T_wall
    property = T_wall
  []
[]

[Materials]
  # wall friction closure
  [f_mat]
    type = ADWallFrictionChurchillMaterial
    block = channel
    D_h = D_h
    f_D = f_D
    vel = vel
    rho = rho
    mu = mu
  []

  # Wall heat transfer closure (all important is in Nu_mat)
  [Re_mat]
    type = ADReynoldsNumberMaterial
    block = channel
    Re = Re
    D_h = D_h
    mu = mu
    vel = vel
    rho = rho
  []
  [Pr_mat]
    type = ADPrandtlNumberMaterial
    block = channel
    cp = cp
    mu = mu
    k = k
  []
  [Nu_mat]
    type = ADParsedMaterial
    block = channel
    # Dittus-Boelter
    function = '0.022 * pow(Re, 0.8) * pow(Pr, 0.4)'
    f_name = 'Nu'
    material_property_names = 'Re Pr'
  []
  [Hw_mat]
    type = ADConvectiveHeatTransferCoefficientMaterial
    block = channel
    D_h = D_h
    Nu = Nu
    k = k
  []
  [T_wall]
    type = ADTemperatureWall3EqnMaterial
    Hw = Hw
    T = T
    q_wall = q_wall
  []
[]

[Components]
  [channel]
    type = FlowChannel1Phase
    position = '0 0 ${height}'
    orientation = '0 0 -1'

    A = '${fparse pi * channel_diameter * channel_diameter / 4}'
    D_h = ${channel_diameter}
    length = ${height}
    n_elems = ${num_layers_for_THM}
  []

  [inlet]
    type = InletMassFlowRateTemperature1Phase
    input = 'channel:in'
    m_dot = ${mdot}
    T = ${inlet_T}
  []

  [outlet]
    type = Outlet1Phase
    input = 'channel:out'
    p = ${outlet_P}
  []

  [ht_ext]
    type = HeatTransferFromExternalAppHeatFlux1Phase
    flow_channel = channel
    P_hf = '${fparse channel_diameter * pi}'
  []
[]

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

[Executioner]
  type = Transient
  dt = 0.1
  start_time = 0

  steady_state_detection = true
  steady_state_tolerance = 1e-08

  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-6

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  solve_type = NEWTON
  line_search = basic
[]

[Outputs]
  exodus = true
  print_linear_residuals = false

  [console]
    type = Console
    outlier_variable_norms = false
  []
[]

(htgr/assembly/thm.i)

## HTGR assembly thermal hydraulics simulation
## Application: MOOSE thermal hydraulics module
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

# copy-pasta from common.i
inlet_T = 598.0 # inlet fluid temperature (K)
mdot = '${fparse 117.3 / 12 / 108}' # fluid mass flowrate (kg/s)
outlet_P = 7.1e6 # fluid outlet pressure (Pa)
channel_diameter = 0.016 # diameter of the coolant channels (m)
height = 6.343 # height of the assembly (m)

num_layers_for_THM = 50 # number of elements in the THM model

[GlobalParams]
  initial_p = ${outlet_P}
  initial_T = ${inlet_T}
  initial_vel = '${fparse mdot / outlet_P / 8.3144598 * 4.0e-3 / inlet_T / (pi * channel_diameter * channel_diameter / 4.0)}'

  rdg_slope_reconstruction = full
  closures = none
  fp = helium
[]

[Closures]
  [none]
    type = Closures1PhaseNone
  []
[]

[FluidProperties]
  [helium]
    type = IdealGasFluidProperties
    molar_mass = 4e-3
    gamma = 1.668282 # should correspond to  Cp = 5189 J/kg/K
    k = 0.2556
    mu = 3.22639e-5
  []
[]

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

[AuxKernels]
  [Tw_aux]
    type = ADMaterialRealAux
    block = channel
    variable = T_wall
    property = T_wall
  []
[]

[Materials]
  # wall friction closure
  [f_mat]
    type = ADWallFrictionChurchillMaterial
    block = channel
    D_h = D_h
    f_D = f_D
    vel = vel
    rho = rho
    mu = mu
  []

  # Wall heat transfer closure (all important is in Nu_mat)
  [Re_mat]
    type = ADReynoldsNumberMaterial
    block = channel
    Re = Re
    D_h = D_h
    mu = mu
    vel = vel
    rho = rho
  []
  [Pr_mat]
    type = ADPrandtlNumberMaterial
    block = channel
    cp = cp
    mu = mu
    k = k
  []
  [Nu_mat]
    type = ADParsedMaterial
    block = channel
    # Dittus-Boelter
    function = '0.022 * pow(Re, 0.8) * pow(Pr, 0.4)'
    f_name = 'Nu'
    material_property_names = 'Re Pr'
  []
  [Hw_mat]
    type = ADConvectiveHeatTransferCoefficientMaterial
    block = channel
    D_h = D_h
    Nu = Nu
    k = k
  []
  [T_wall]
    type = ADTemperatureWall3EqnMaterial
    Hw = Hw
    T = T
    q_wall = q_wall
  []
[]

[Components]
  [channel]
    type = FlowChannel1Phase
    position = '0 0 ${height}'
    orientation = '0 0 -1'

    A = '${fparse pi * channel_diameter * channel_diameter / 4}'
    D_h = ${channel_diameter}
    length = ${height}
    n_elems = ${num_layers_for_THM}
  []

  [inlet]
    type = InletMassFlowRateTemperature1Phase
    input = 'channel:in'
    m_dot = ${mdot}
    T = ${inlet_T}
  []

  [outlet]
    type = Outlet1Phase
    input = 'channel:out'
    p = ${outlet_P}
  []

  [ht_ext]
    type = HeatTransferFromExternalAppHeatFlux1Phase
    flow_channel = channel
    P_hf = '${fparse channel_diameter * pi}'
  []
[]

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

[Executioner]
  type = Transient
  dt = 0.1
  start_time = 0

  steady_state_detection = true
  steady_state_tolerance = 1e-08

  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-6

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  solve_type = NEWTON
  line_search = basic
[]

[Outputs]
  exodus = true
  print_linear_residuals = false

  [console]
    type = Console
    outlier_variable_norms = false
  []
[]

(htgr/assembly/thm.i)

## HTGR assembly thermal hydraulics simulation
## Application: MOOSE thermal hydraulics module
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

# copy-pasta from common.i
inlet_T = 598.0 # inlet fluid temperature (K)
mdot = '${fparse 117.3 / 12 / 108}' # fluid mass flowrate (kg/s)
outlet_P = 7.1e6 # fluid outlet pressure (Pa)
channel_diameter = 0.016 # diameter of the coolant channels (m)
height = 6.343 # height of the assembly (m)

num_layers_for_THM = 50 # number of elements in the THM model

[GlobalParams]
  initial_p = ${outlet_P}
  initial_T = ${inlet_T}
  initial_vel = '${fparse mdot / outlet_P / 8.3144598 * 4.0e-3 / inlet_T / (pi * channel_diameter * channel_diameter / 4.0)}'

  rdg_slope_reconstruction = full
  closures = none
  fp = helium
[]

[Closures]
  [none]
    type = Closures1PhaseNone
  []
[]

[FluidProperties]
  [helium]
    type = IdealGasFluidProperties
    molar_mass = 4e-3
    gamma = 1.668282 # should correspond to  Cp = 5189 J/kg/K
    k = 0.2556
    mu = 3.22639e-5
  []
[]

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

[AuxKernels]
  [Tw_aux]
    type = ADMaterialRealAux
    block = channel
    variable = T_wall
    property = T_wall
  []
[]

[Materials]
  # wall friction closure
  [f_mat]
    type = ADWallFrictionChurchillMaterial
    block = channel
    D_h = D_h
    f_D = f_D
    vel = vel
    rho = rho
    mu = mu
  []

  # Wall heat transfer closure (all important is in Nu_mat)
  [Re_mat]
    type = ADReynoldsNumberMaterial
    block = channel
    Re = Re
    D_h = D_h
    mu = mu
    vel = vel
    rho = rho
  []
  [Pr_mat]
    type = ADPrandtlNumberMaterial
    block = channel
    cp = cp
    mu = mu
    k = k
  []
  [Nu_mat]
    type = ADParsedMaterial
    block = channel
    # Dittus-Boelter
    function = '0.022 * pow(Re, 0.8) * pow(Pr, 0.4)'
    f_name = 'Nu'
    material_property_names = 'Re Pr'
  []
  [Hw_mat]
    type = ADConvectiveHeatTransferCoefficientMaterial
    block = channel
    D_h = D_h
    Nu = Nu
    k = k
  []
  [T_wall]
    type = ADTemperatureWall3EqnMaterial
    Hw = Hw
    T = T
    q_wall = q_wall
  []
[]

[Components]
  [channel]
    type = FlowChannel1Phase
    position = '0 0 ${height}'
    orientation = '0 0 -1'

    A = '${fparse pi * channel_diameter * channel_diameter / 4}'
    D_h = ${channel_diameter}
    length = ${height}
    n_elems = ${num_layers_for_THM}
  []

  [inlet]
    type = InletMassFlowRateTemperature1Phase
    input = 'channel:in'
    m_dot = ${mdot}
    T = ${inlet_T}
  []

  [outlet]
    type = Outlet1Phase
    input = 'channel:out'
    p = ${outlet_P}
  []

  [ht_ext]
    type = HeatTransferFromExternalAppHeatFlux1Phase
    flow_channel = channel
    P_hf = '${fparse channel_diameter * pi}'
  []
[]

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

[Executioner]
  type = Transient
  dt = 0.1
  start_time = 0

  steady_state_detection = true
  steady_state_tolerance = 1e-08

  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-6

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  solve_type = NEWTON
  line_search = basic
[]

[Outputs]
  exodus = true
  print_linear_residuals = false

  [console]
    type = Console
    outlier_variable_norms = false
  []
[]

(htgr/assembly/thm.i)

## HTGR assembly thermal hydraulics simulation
## Application: MOOSE thermal hydraulics module
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

# copy-pasta from common.i
inlet_T = 598.0 # inlet fluid temperature (K)
mdot = '${fparse 117.3 / 12 / 108}' # fluid mass flowrate (kg/s)
outlet_P = 7.1e6 # fluid outlet pressure (Pa)
channel_diameter = 0.016 # diameter of the coolant channels (m)
height = 6.343 # height of the assembly (m)

num_layers_for_THM = 50 # number of elements in the THM model

[GlobalParams]
  initial_p = ${outlet_P}
  initial_T = ${inlet_T}
  initial_vel = '${fparse mdot / outlet_P / 8.3144598 * 4.0e-3 / inlet_T / (pi * channel_diameter * channel_diameter / 4.0)}'

  rdg_slope_reconstruction = full
  closures = none
  fp = helium
[]

[Closures]
  [none]
    type = Closures1PhaseNone
  []
[]

[FluidProperties]
  [helium]
    type = IdealGasFluidProperties
    molar_mass = 4e-3
    gamma = 1.668282 # should correspond to  Cp = 5189 J/kg/K
    k = 0.2556
    mu = 3.22639e-5
  []
[]

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

[AuxKernels]
  [Tw_aux]
    type = ADMaterialRealAux
    block = channel
    variable = T_wall
    property = T_wall
  []
[]

[Materials]
  # wall friction closure
  [f_mat]
    type = ADWallFrictionChurchillMaterial
    block = channel
    D_h = D_h
    f_D = f_D
    vel = vel
    rho = rho
    mu = mu
  []

  # Wall heat transfer closure (all important is in Nu_mat)
  [Re_mat]
    type = ADReynoldsNumberMaterial
    block = channel
    Re = Re
    D_h = D_h
    mu = mu
    vel = vel
    rho = rho
  []
  [Pr_mat]
    type = ADPrandtlNumberMaterial
    block = channel
    cp = cp
    mu = mu
    k = k
  []
  [Nu_mat]
    type = ADParsedMaterial
    block = channel
    # Dittus-Boelter
    function = '0.022 * pow(Re, 0.8) * pow(Pr, 0.4)'
    f_name = 'Nu'
    material_property_names = 'Re Pr'
  []
  [Hw_mat]
    type = ADConvectiveHeatTransferCoefficientMaterial
    block = channel
    D_h = D_h
    Nu = Nu
    k = k
  []
  [T_wall]
    type = ADTemperatureWall3EqnMaterial
    Hw = Hw
    T = T
    q_wall = q_wall
  []
[]

[Components]
  [channel]
    type = FlowChannel1Phase
    position = '0 0 ${height}'
    orientation = '0 0 -1'

    A = '${fparse pi * channel_diameter * channel_diameter / 4}'
    D_h = ${channel_diameter}
    length = ${height}
    n_elems = ${num_layers_for_THM}
  []

  [inlet]
    type = InletMassFlowRateTemperature1Phase
    input = 'channel:in'
    m_dot = ${mdot}
    T = ${inlet_T}
  []

  [outlet]
    type = Outlet1Phase
    input = 'channel:out'
    p = ${outlet_P}
  []

  [ht_ext]
    type = HeatTransferFromExternalAppHeatFlux1Phase
    flow_channel = channel
    P_hf = '${fparse channel_diameter * pi}'
  []
[]

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

[Executioner]
  type = Transient
  dt = 0.1
  start_time = 0

  steady_state_detection = true
  steady_state_tolerance = 1e-08

  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-6

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  solve_type = NEWTON
  line_search = basic
[]

[Outputs]
  exodus = true
  print_linear_residuals = false

  [console]
    type = Console
    outlier_variable_norms = false
  []
[]

(htgr/assembly/thm.i)

## HTGR assembly thermal hydraulics simulation
## Application: MOOSE thermal hydraulics module
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

# copy-pasta from common.i
inlet_T = 598.0 # inlet fluid temperature (K)
mdot = '${fparse 117.3 / 12 / 108}' # fluid mass flowrate (kg/s)
outlet_P = 7.1e6 # fluid outlet pressure (Pa)
channel_diameter = 0.016 # diameter of the coolant channels (m)
height = 6.343 # height of the assembly (m)

num_layers_for_THM = 50 # number of elements in the THM model

[GlobalParams]
  initial_p = ${outlet_P}
  initial_T = ${inlet_T}
  initial_vel = '${fparse mdot / outlet_P / 8.3144598 * 4.0e-3 / inlet_T / (pi * channel_diameter * channel_diameter / 4.0)}'

  rdg_slope_reconstruction = full
  closures = none
  fp = helium
[]

[Closures]
  [none]
    type = Closures1PhaseNone
  []
[]

[FluidProperties]
  [helium]
    type = IdealGasFluidProperties
    molar_mass = 4e-3
    gamma = 1.668282 # should correspond to  Cp = 5189 J/kg/K
    k = 0.2556
    mu = 3.22639e-5
  []
[]

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

[AuxKernels]
  [Tw_aux]
    type = ADMaterialRealAux
    block = channel
    variable = T_wall
    property = T_wall
  []
[]

[Materials]
  # wall friction closure
  [f_mat]
    type = ADWallFrictionChurchillMaterial
    block = channel
    D_h = D_h
    f_D = f_D
    vel = vel
    rho = rho
    mu = mu
  []

  # Wall heat transfer closure (all important is in Nu_mat)
  [Re_mat]
    type = ADReynoldsNumberMaterial
    block = channel
    Re = Re
    D_h = D_h
    mu = mu
    vel = vel
    rho = rho
  []
  [Pr_mat]
    type = ADPrandtlNumberMaterial
    block = channel
    cp = cp
    mu = mu
    k = k
  []
  [Nu_mat]
    type = ADParsedMaterial
    block = channel
    # Dittus-Boelter
    function = '0.022 * pow(Re, 0.8) * pow(Pr, 0.4)'
    f_name = 'Nu'
    material_property_names = 'Re Pr'
  []
  [Hw_mat]
    type = ADConvectiveHeatTransferCoefficientMaterial
    block = channel
    D_h = D_h
    Nu = Nu
    k = k
  []
  [T_wall]
    type = ADTemperatureWall3EqnMaterial
    Hw = Hw
    T = T
    q_wall = q_wall
  []
[]

[Components]
  [channel]
    type = FlowChannel1Phase
    position = '0 0 ${height}'
    orientation = '0 0 -1'

    A = '${fparse pi * channel_diameter * channel_diameter / 4}'
    D_h = ${channel_diameter}
    length = ${height}
    n_elems = ${num_layers_for_THM}
  []

  [inlet]
    type = InletMassFlowRateTemperature1Phase
    input = 'channel:in'
    m_dot = ${mdot}
    T = ${inlet_T}
  []

  [outlet]
    type = Outlet1Phase
    input = 'channel:out'
    p = ${outlet_P}
  []

  [ht_ext]
    type = HeatTransferFromExternalAppHeatFlux1Phase
    flow_channel = channel
    P_hf = '${fparse channel_diameter * pi}'
  []
[]

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

[Executioner]
  type = Transient
  dt = 0.1
  start_time = 0

  steady_state_detection = true
  steady_state_tolerance = 1e-08

  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-6

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  solve_type = NEWTON
  line_search = basic
[]

[Outputs]
  exodus = true
  print_linear_residuals = false

  [console]
    type = Console
    outlier_variable_norms = false
  []
[]

(htgr/assembly/openmc.i)

## HTGR assembly neutronics simulation
## Application: Cardinal
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

num_layers_for_THM = 50 # number of elements in the THM model; for the converged
# case, we set this to 150

[GlobalParams]
  # Reduces transfers efficiency for now, can be removed once transferred fields are checked
  bbox_factor = 10
[]

[Mesh]
  # mesh mirror for the solid regions
  [solid]
    type = FileMeshGenerator
    file = meshes/solid_in.e
  []

  # create a mesh for a single coolant channel; because we will receive uniform
  # temperatures and densities from THM on each x-y plane, we can use a very coarse
  # mesh in the radial direction
  [coolant_face]
    type = AnnularMeshGenerator
    nr = 1
    nt = 8
    rmin = 0.0
    rmax = '${fparse channel_diameter / 2.0}'
  []
  [extrude]
    type = AdvancedExtruderGenerator
    input = coolant_face
    num_layers = ${num_layers_for_THM}
    direction = '0 0 1'
    heights = '${height}'
    top_boundary = '300' # inlet
    bottom_boundary = '400' # outlet
  []
  [rename]
    type = RenameBlockGenerator
    input = extrude
    old_block = '1'
    new_block = '101'
  []

  # repeat the coolant channels and then combine together to get a combined mesh mirror
  [repeat]
    type = CombinerGenerator
    inputs = rename
    positions_file = coolant_channel_positions.txt
  []
  [add]
    type = CombinerGenerator
    inputs = 'solid repeat'
  []
[]

[AuxVariables]
  [material_id]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_density]
    family = MONOMIAL
    order = CONSTANT
  []
  [thm_temp_wall]
    block = '101'
  []
  [flux]
  []

  # just for postprocessing purposes
  [thm_pressure]
    block = '101'
  []
  [thm_velocity]
    block = '101'
  []
  [z]
    family = MONOMIAL
    order = CONSTANT
    block = 'compacts'
  []
[]

[AuxKernels]
  [material_id]
    type = CellMaterialIDAux
    variable = material_id
  []
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
  [density]
    type = FluidDensityAux
    variable = density
    p = ${outlet_P}
    T = thm_temp
    fp = helium
    execute_on = 'timestep_begin linear'
  []
  [z]
    type = ParsedAux
    variable = z
    use_xyzt = true
    expression = 'z'
  []
[]

[ICs]
  [fluid_temp_wall]
    type = FunctionIC
    variable = thm_temp_wall
    function = temp_ic
  []
  [fluid_temp]
    type = FunctionIC
    variable = thm_temp
    function = temp_ic
  []
  [heat_source]
    type = ConstantIC
    variable = heat_source
    block = 'compacts'
    value = '${fparse power / (n_bundles * n_fuel_compacts_per_block) / (pi * compact_diameter * compact_diameter / 4.0 * height)}'
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + (${height} - z) / ${height} * ${power} / ${mdot} / ${fluid_Cp}'
  []
[]

[FluidProperties]
  [helium]
    type = IdealGasFluidProperties
    molar_mass = 4e-3
    gamma = 1.668282 # should correspond to  Cp = 5189 J/kg/K
    k = 0.2556
    mu = 3.22639e-5
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  # optimizations to increase tracking rate by assuring that the tallies
  # are spatially independent
  check_tally_sum = false
  normalize_by_global_tally = false
  assume_separate_tallies = true

  power = '${fparse power / n_bundles}'
  scaling = 100.0

  relaxation = constant
  relaxation_factor = 0.5

  inactive_batches = 200
  batches = 1000

  # we will collate temperature from THM (for the fluid) and MOOSE (for the solid)
  # into variables we name as 'solid_temp' and 'thm_temp'. This syntax will automatically
  # create those variabes for us
  temperature_variables = 'solid_temp; solid_temp; solid_temp; thm_temp'
  temperature_blocks = 'graphite; compacts; poison; 101'
  cell_level = 0

  density_blocks = '101'

  [Tallies]
    [heat_source]
      type = CellTally
      block = '2'
      name = heat_source
      score = heating

      trigger = rel_err
      trigger_threshold = 2e-2

      check_equal_mapped_tally_volumes = true
      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Postprocessors]
  [flux_integral]
    type = SideIntegralVariablePostprocessor
    variable = flux
    boundary = 'fluid_solid_interface'
    execute_on = 'transfer linear'
  []
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_rel_err]
    type = TallyRelativeError
    value_type = max
    tally_score = 'heating'
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KEigenvalue
    output = 'std_dev'
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = 'compacts'
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = 'compacts'
  []
  [z_max_power]
    type = ElementExtremeValue
    proxy_variable = heat_source
    variable = z
    block = 'compacts'
  []
  [max_Tf]
    type = ElementExtremeValue
    variable = thm_temp
    block = '101'
  []
  [P_in]
    type = SideAverageValue
    variable = thm_pressure
    boundary = '300'
  []
  [pressure_drop]
    type = LinearCombinationPostprocessor
    pp_names = 'P_in'
    pp_coefs = '1.0'
    b = '${fparse -outlet_P}'
  []
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    app_type = CardinalApp
    input_files = 'solid.i'
    execute_on = timestep_begin
  []
  [thm]
    type = FullSolveMultiApp
    app_type = CardinalApp
    input_files = 'thm.i'
    execute_on = timestep_end
    max_procs_per_app = 1
    bounding_box_padding = '0.1 0.1 0'
    positions_file = coolant_channel_positions.txt
    output_in_position = true
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [heat_flux_to_openmc]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = flux
    variable = flux
    from_multi_app = bison
    source_boundary = 'fluid_solid_interface'
    target_boundary = 'fluid_solid_interface'
    from_postprocessors_to_be_preserved = flux_integral
    to_postprocessors_to_be_preserved = flux_integral
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    direction = to_multiapp
    multi_app = bison
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = power
  []
  [thm_temp_to_bison]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = thm_temp_wall
    variable = thm_temp
    direction = to_multiapp
    multi_app = bison
  []

  [q_wall_to_thm]
    type = MultiAppGeneralFieldUserObjectTransfer
    variable = q_wall
    direction = to_multiapp
    multi_app = thm
    source_user_object = q_wall_avg
  []
  [T_wall_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = T_wall
    direction = from_multiapp
    multi_app = thm
    variable = thm_temp_wall
    target_boundary = 'fluid_solid_interface'
  []
  [T_bulk_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = T
    direction = from_multiapp
    multi_app = thm
    variable = thm_temp
  []

  # just for postprocessing purposes
  [pressure_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = p
    direction = from_multiapp
    multi_app = thm
    variable = thm_pressure
  []
  [velocity_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = vel_z
    direction = from_multiapp
    multi_app = thm
    variable = thm_velocity
  []
[]

[UserObjects]
  [q_wall_avg]
    type = NearestPointLayeredSideAverage
    boundary = fluid_solid_interface
    variable = flux

    # Note: make this to match the num_elems in the channel
    direction = z
    num_layers = ${num_layers_for_THM}
    points_file = coolant_channel_positions.txt

    direction_min = 0.0
    direction_max = ${height}
  []
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_fluid_axial]
    type = LayeredAverage
    variable = thm_temp
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
  [average_pressure]
    type = LayeredAverage
    variable = thm_pressure
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
  [average_axial_velocity]
    type = LayeredAverage
    variable = thm_velocity
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
  [fluid_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fluid_axial
  []
  [pressure_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_pressure
  []
  [velocity_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_axial_velocity
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Outputs]
  exodus = true
  csv = true
  hide = 'P_in flux_integral z'
[]

(htgr/assembly/openmc.i)

## HTGR assembly neutronics simulation
## Application: Cardinal
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

num_layers_for_THM = 50 # number of elements in the THM model; for the converged
# case, we set this to 150

[GlobalParams]
  # Reduces transfers efficiency for now, can be removed once transferred fields are checked
  bbox_factor = 10
[]

[Mesh]
  # mesh mirror for the solid regions
  [solid]
    type = FileMeshGenerator
    file = meshes/solid_in.e
  []

  # create a mesh for a single coolant channel; because we will receive uniform
  # temperatures and densities from THM on each x-y plane, we can use a very coarse
  # mesh in the radial direction
  [coolant_face]
    type = AnnularMeshGenerator
    nr = 1
    nt = 8
    rmin = 0.0
    rmax = '${fparse channel_diameter / 2.0}'
  []
  [extrude]
    type = AdvancedExtruderGenerator
    input = coolant_face
    num_layers = ${num_layers_for_THM}
    direction = '0 0 1'
    heights = '${height}'
    top_boundary = '300' # inlet
    bottom_boundary = '400' # outlet
  []
  [rename]
    type = RenameBlockGenerator
    input = extrude
    old_block = '1'
    new_block = '101'
  []

  # repeat the coolant channels and then combine together to get a combined mesh mirror
  [repeat]
    type = CombinerGenerator
    inputs = rename
    positions_file = coolant_channel_positions.txt
  []
  [add]
    type = CombinerGenerator
    inputs = 'solid repeat'
  []
[]

[AuxVariables]
  [material_id]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_density]
    family = MONOMIAL
    order = CONSTANT
  []
  [thm_temp_wall]
    block = '101'
  []
  [flux]
  []

  # just for postprocessing purposes
  [thm_pressure]
    block = '101'
  []
  [thm_velocity]
    block = '101'
  []
  [z]
    family = MONOMIAL
    order = CONSTANT
    block = 'compacts'
  []
[]

[AuxKernels]
  [material_id]
    type = CellMaterialIDAux
    variable = material_id
  []
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
  [density]
    type = FluidDensityAux
    variable = density
    p = ${outlet_P}
    T = thm_temp
    fp = helium
    execute_on = 'timestep_begin linear'
  []
  [z]
    type = ParsedAux
    variable = z
    use_xyzt = true
    expression = 'z'
  []
[]

[ICs]
  [fluid_temp_wall]
    type = FunctionIC
    variable = thm_temp_wall
    function = temp_ic
  []
  [fluid_temp]
    type = FunctionIC
    variable = thm_temp
    function = temp_ic
  []
  [heat_source]
    type = ConstantIC
    variable = heat_source
    block = 'compacts'
    value = '${fparse power / (n_bundles * n_fuel_compacts_per_block) / (pi * compact_diameter * compact_diameter / 4.0 * height)}'
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + (${height} - z) / ${height} * ${power} / ${mdot} / ${fluid_Cp}'
  []
[]

[FluidProperties]
  [helium]
    type = IdealGasFluidProperties
    molar_mass = 4e-3
    gamma = 1.668282 # should correspond to  Cp = 5189 J/kg/K
    k = 0.2556
    mu = 3.22639e-5
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  # optimizations to increase tracking rate by assuring that the tallies
  # are spatially independent
  check_tally_sum = false
  normalize_by_global_tally = false
  assume_separate_tallies = true

  power = '${fparse power / n_bundles}'
  scaling = 100.0

  relaxation = constant
  relaxation_factor = 0.5

  inactive_batches = 200
  batches = 1000

  # we will collate temperature from THM (for the fluid) and MOOSE (for the solid)
  # into variables we name as 'solid_temp' and 'thm_temp'. This syntax will automatically
  # create those variabes for us
  temperature_variables = 'solid_temp; solid_temp; solid_temp; thm_temp'
  temperature_blocks = 'graphite; compacts; poison; 101'
  cell_level = 0

  density_blocks = '101'

  [Tallies]
    [heat_source]
      type = CellTally
      block = '2'
      name = heat_source
      score = heating

      trigger = rel_err
      trigger_threshold = 2e-2

      check_equal_mapped_tally_volumes = true
      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Postprocessors]
  [flux_integral]
    type = SideIntegralVariablePostprocessor
    variable = flux
    boundary = 'fluid_solid_interface'
    execute_on = 'transfer linear'
  []
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_rel_err]
    type = TallyRelativeError
    value_type = max
    tally_score = 'heating'
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KEigenvalue
    output = 'std_dev'
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = 'compacts'
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = 'compacts'
  []
  [z_max_power]
    type = ElementExtremeValue
    proxy_variable = heat_source
    variable = z
    block = 'compacts'
  []
  [max_Tf]
    type = ElementExtremeValue
    variable = thm_temp
    block = '101'
  []
  [P_in]
    type = SideAverageValue
    variable = thm_pressure
    boundary = '300'
  []
  [pressure_drop]
    type = LinearCombinationPostprocessor
    pp_names = 'P_in'
    pp_coefs = '1.0'
    b = '${fparse -outlet_P}'
  []
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    app_type = CardinalApp
    input_files = 'solid.i'
    execute_on = timestep_begin
  []
  [thm]
    type = FullSolveMultiApp
    app_type = CardinalApp
    input_files = 'thm.i'
    execute_on = timestep_end
    max_procs_per_app = 1
    bounding_box_padding = '0.1 0.1 0'
    positions_file = coolant_channel_positions.txt
    output_in_position = true
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [heat_flux_to_openmc]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = flux
    variable = flux
    from_multi_app = bison
    source_boundary = 'fluid_solid_interface'
    target_boundary = 'fluid_solid_interface'
    from_postprocessors_to_be_preserved = flux_integral
    to_postprocessors_to_be_preserved = flux_integral
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    direction = to_multiapp
    multi_app = bison
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = power
  []
  [thm_temp_to_bison]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = thm_temp_wall
    variable = thm_temp
    direction = to_multiapp
    multi_app = bison
  []

  [q_wall_to_thm]
    type = MultiAppGeneralFieldUserObjectTransfer
    variable = q_wall
    direction = to_multiapp
    multi_app = thm
    source_user_object = q_wall_avg
  []
  [T_wall_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = T_wall
    direction = from_multiapp
    multi_app = thm
    variable = thm_temp_wall
    target_boundary = 'fluid_solid_interface'
  []
  [T_bulk_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = T
    direction = from_multiapp
    multi_app = thm
    variable = thm_temp
  []

  # just for postprocessing purposes
  [pressure_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = p
    direction = from_multiapp
    multi_app = thm
    variable = thm_pressure
  []
  [velocity_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = vel_z
    direction = from_multiapp
    multi_app = thm
    variable = thm_velocity
  []
[]

[UserObjects]
  [q_wall_avg]
    type = NearestPointLayeredSideAverage
    boundary = fluid_solid_interface
    variable = flux

    # Note: make this to match the num_elems in the channel
    direction = z
    num_layers = ${num_layers_for_THM}
    points_file = coolant_channel_positions.txt

    direction_min = 0.0
    direction_max = ${height}
  []
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_fluid_axial]
    type = LayeredAverage
    variable = thm_temp
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
  [average_pressure]
    type = LayeredAverage
    variable = thm_pressure
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
  [average_axial_velocity]
    type = LayeredAverage
    variable = thm_velocity
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
  [fluid_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fluid_axial
  []
  [pressure_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_pressure
  []
  [velocity_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_axial_velocity
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Outputs]
  exodus = true
  csv = true
  hide = 'P_in flux_integral z'
[]

(htgr/assembly/openmc.i)

## HTGR assembly neutronics simulation
## Application: Cardinal
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

num_layers_for_THM = 50 # number of elements in the THM model; for the converged
# case, we set this to 150

[GlobalParams]
  # Reduces transfers efficiency for now, can be removed once transferred fields are checked
  bbox_factor = 10
[]

[Mesh]
  # mesh mirror for the solid regions
  [solid]
    type = FileMeshGenerator
    file = meshes/solid_in.e
  []

  # create a mesh for a single coolant channel; because we will receive uniform
  # temperatures and densities from THM on each x-y plane, we can use a very coarse
  # mesh in the radial direction
  [coolant_face]
    type = AnnularMeshGenerator
    nr = 1
    nt = 8
    rmin = 0.0
    rmax = '${fparse channel_diameter / 2.0}'
  []
  [extrude]
    type = AdvancedExtruderGenerator
    input = coolant_face
    num_layers = ${num_layers_for_THM}
    direction = '0 0 1'
    heights = '${height}'
    top_boundary = '300' # inlet
    bottom_boundary = '400' # outlet
  []
  [rename]
    type = RenameBlockGenerator
    input = extrude
    old_block = '1'
    new_block = '101'
  []

  # repeat the coolant channels and then combine together to get a combined mesh mirror
  [repeat]
    type = CombinerGenerator
    inputs = rename
    positions_file = coolant_channel_positions.txt
  []
  [add]
    type = CombinerGenerator
    inputs = 'solid repeat'
  []
[]

[AuxVariables]
  [material_id]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_density]
    family = MONOMIAL
    order = CONSTANT
  []
  [thm_temp_wall]
    block = '101'
  []
  [flux]
  []

  # just for postprocessing purposes
  [thm_pressure]
    block = '101'
  []
  [thm_velocity]
    block = '101'
  []
  [z]
    family = MONOMIAL
    order = CONSTANT
    block = 'compacts'
  []
[]

[AuxKernels]
  [material_id]
    type = CellMaterialIDAux
    variable = material_id
  []
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
  [density]
    type = FluidDensityAux
    variable = density
    p = ${outlet_P}
    T = thm_temp
    fp = helium
    execute_on = 'timestep_begin linear'
  []
  [z]
    type = ParsedAux
    variable = z
    use_xyzt = true
    expression = 'z'
  []
[]

[ICs]
  [fluid_temp_wall]
    type = FunctionIC
    variable = thm_temp_wall
    function = temp_ic
  []
  [fluid_temp]
    type = FunctionIC
    variable = thm_temp
    function = temp_ic
  []
  [heat_source]
    type = ConstantIC
    variable = heat_source
    block = 'compacts'
    value = '${fparse power / (n_bundles * n_fuel_compacts_per_block) / (pi * compact_diameter * compact_diameter / 4.0 * height)}'
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + (${height} - z) / ${height} * ${power} / ${mdot} / ${fluid_Cp}'
  []
[]

[FluidProperties]
  [helium]
    type = IdealGasFluidProperties
    molar_mass = 4e-3
    gamma = 1.668282 # should correspond to  Cp = 5189 J/kg/K
    k = 0.2556
    mu = 3.22639e-5
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  # optimizations to increase tracking rate by assuring that the tallies
  # are spatially independent
  check_tally_sum = false
  normalize_by_global_tally = false
  assume_separate_tallies = true

  power = '${fparse power / n_bundles}'
  scaling = 100.0

  relaxation = constant
  relaxation_factor = 0.5

  inactive_batches = 200
  batches = 1000

  # we will collate temperature from THM (for the fluid) and MOOSE (for the solid)
  # into variables we name as 'solid_temp' and 'thm_temp'. This syntax will automatically
  # create those variabes for us
  temperature_variables = 'solid_temp; solid_temp; solid_temp; thm_temp'
  temperature_blocks = 'graphite; compacts; poison; 101'
  cell_level = 0

  density_blocks = '101'

  [Tallies]
    [heat_source]
      type = CellTally
      block = '2'
      name = heat_source
      score = heating

      trigger = rel_err
      trigger_threshold = 2e-2

      check_equal_mapped_tally_volumes = true
      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Postprocessors]
  [flux_integral]
    type = SideIntegralVariablePostprocessor
    variable = flux
    boundary = 'fluid_solid_interface'
    execute_on = 'transfer linear'
  []
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_rel_err]
    type = TallyRelativeError
    value_type = max
    tally_score = 'heating'
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KEigenvalue
    output = 'std_dev'
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = 'compacts'
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = 'compacts'
  []
  [z_max_power]
    type = ElementExtremeValue
    proxy_variable = heat_source
    variable = z
    block = 'compacts'
  []
  [max_Tf]
    type = ElementExtremeValue
    variable = thm_temp
    block = '101'
  []
  [P_in]
    type = SideAverageValue
    variable = thm_pressure
    boundary = '300'
  []
  [pressure_drop]
    type = LinearCombinationPostprocessor
    pp_names = 'P_in'
    pp_coefs = '1.0'
    b = '${fparse -outlet_P}'
  []
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    app_type = CardinalApp
    input_files = 'solid.i'
    execute_on = timestep_begin
  []
  [thm]
    type = FullSolveMultiApp
    app_type = CardinalApp
    input_files = 'thm.i'
    execute_on = timestep_end
    max_procs_per_app = 1
    bounding_box_padding = '0.1 0.1 0'
    positions_file = coolant_channel_positions.txt
    output_in_position = true
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [heat_flux_to_openmc]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = flux
    variable = flux
    from_multi_app = bison
    source_boundary = 'fluid_solid_interface'
    target_boundary = 'fluid_solid_interface'
    from_postprocessors_to_be_preserved = flux_integral
    to_postprocessors_to_be_preserved = flux_integral
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    direction = to_multiapp
    multi_app = bison
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = power
  []
  [thm_temp_to_bison]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = thm_temp_wall
    variable = thm_temp
    direction = to_multiapp
    multi_app = bison
  []

  [q_wall_to_thm]
    type = MultiAppGeneralFieldUserObjectTransfer
    variable = q_wall
    direction = to_multiapp
    multi_app = thm
    source_user_object = q_wall_avg
  []
  [T_wall_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = T_wall
    direction = from_multiapp
    multi_app = thm
    variable = thm_temp_wall
    target_boundary = 'fluid_solid_interface'
  []
  [T_bulk_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = T
    direction = from_multiapp
    multi_app = thm
    variable = thm_temp
  []

  # just for postprocessing purposes
  [pressure_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = p
    direction = from_multiapp
    multi_app = thm
    variable = thm_pressure
  []
  [velocity_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = vel_z
    direction = from_multiapp
    multi_app = thm
    variable = thm_velocity
  []
[]

[UserObjects]
  [q_wall_avg]
    type = NearestPointLayeredSideAverage
    boundary = fluid_solid_interface
    variable = flux

    # Note: make this to match the num_elems in the channel
    direction = z
    num_layers = ${num_layers_for_THM}
    points_file = coolant_channel_positions.txt

    direction_min = 0.0
    direction_max = ${height}
  []
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_fluid_axial]
    type = LayeredAverage
    variable = thm_temp
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
  [average_pressure]
    type = LayeredAverage
    variable = thm_pressure
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
  [average_axial_velocity]
    type = LayeredAverage
    variable = thm_velocity
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
  [fluid_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fluid_axial
  []
  [pressure_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_pressure
  []
  [velocity_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_axial_velocity
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Outputs]
  exodus = true
  csv = true
  hide = 'P_in flux_integral z'
[]

(htgr/assembly/openmc.i)

## HTGR assembly neutronics simulation
## Application: Cardinal
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

num_layers_for_THM = 50 # number of elements in the THM model; for the converged
# case, we set this to 150

[GlobalParams]
  # Reduces transfers efficiency for now, can be removed once transferred fields are checked
  bbox_factor = 10
[]

[Mesh]
  # mesh mirror for the solid regions
  [solid]
    type = FileMeshGenerator
    file = meshes/solid_in.e
  []

  # create a mesh for a single coolant channel; because we will receive uniform
  # temperatures and densities from THM on each x-y plane, we can use a very coarse
  # mesh in the radial direction
  [coolant_face]
    type = AnnularMeshGenerator
    nr = 1
    nt = 8
    rmin = 0.0
    rmax = '${fparse channel_diameter / 2.0}'
  []
  [extrude]
    type = AdvancedExtruderGenerator
    input = coolant_face
    num_layers = ${num_layers_for_THM}
    direction = '0 0 1'
    heights = '${height}'
    top_boundary = '300' # inlet
    bottom_boundary = '400' # outlet
  []
  [rename]
    type = RenameBlockGenerator
    input = extrude
    old_block = '1'
    new_block = '101'
  []

  # repeat the coolant channels and then combine together to get a combined mesh mirror
  [repeat]
    type = CombinerGenerator
    inputs = rename
    positions_file = coolant_channel_positions.txt
  []
  [add]
    type = CombinerGenerator
    inputs = 'solid repeat'
  []
[]

[AuxVariables]
  [material_id]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_density]
    family = MONOMIAL
    order = CONSTANT
  []
  [thm_temp_wall]
    block = '101'
  []
  [flux]
  []

  # just for postprocessing purposes
  [thm_pressure]
    block = '101'
  []
  [thm_velocity]
    block = '101'
  []
  [z]
    family = MONOMIAL
    order = CONSTANT
    block = 'compacts'
  []
[]

[AuxKernels]
  [material_id]
    type = CellMaterialIDAux
    variable = material_id
  []
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
  [density]
    type = FluidDensityAux
    variable = density
    p = ${outlet_P}
    T = thm_temp
    fp = helium
    execute_on = 'timestep_begin linear'
  []
  [z]
    type = ParsedAux
    variable = z
    use_xyzt = true
    expression = 'z'
  []
[]

[ICs]
  [fluid_temp_wall]
    type = FunctionIC
    variable = thm_temp_wall
    function = temp_ic
  []
  [fluid_temp]
    type = FunctionIC
    variable = thm_temp
    function = temp_ic
  []
  [heat_source]
    type = ConstantIC
    variable = heat_source
    block = 'compacts'
    value = '${fparse power / (n_bundles * n_fuel_compacts_per_block) / (pi * compact_diameter * compact_diameter / 4.0 * height)}'
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + (${height} - z) / ${height} * ${power} / ${mdot} / ${fluid_Cp}'
  []
[]

[FluidProperties]
  [helium]
    type = IdealGasFluidProperties
    molar_mass = 4e-3
    gamma = 1.668282 # should correspond to  Cp = 5189 J/kg/K
    k = 0.2556
    mu = 3.22639e-5
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  # optimizations to increase tracking rate by assuring that the tallies
  # are spatially independent
  check_tally_sum = false
  normalize_by_global_tally = false
  assume_separate_tallies = true

  power = '${fparse power / n_bundles}'
  scaling = 100.0

  relaxation = constant
  relaxation_factor = 0.5

  inactive_batches = 200
  batches = 1000

  # we will collate temperature from THM (for the fluid) and MOOSE (for the solid)
  # into variables we name as 'solid_temp' and 'thm_temp'. This syntax will automatically
  # create those variabes for us
  temperature_variables = 'solid_temp; solid_temp; solid_temp; thm_temp'
  temperature_blocks = 'graphite; compacts; poison; 101'
  cell_level = 0

  density_blocks = '101'

  [Tallies]
    [heat_source]
      type = CellTally
      block = '2'
      name = heat_source
      score = heating

      trigger = rel_err
      trigger_threshold = 2e-2

      check_equal_mapped_tally_volumes = true
      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Postprocessors]
  [flux_integral]
    type = SideIntegralVariablePostprocessor
    variable = flux
    boundary = 'fluid_solid_interface'
    execute_on = 'transfer linear'
  []
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_rel_err]
    type = TallyRelativeError
    value_type = max
    tally_score = 'heating'
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KEigenvalue
    output = 'std_dev'
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = 'compacts'
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = 'compacts'
  []
  [z_max_power]
    type = ElementExtremeValue
    proxy_variable = heat_source
    variable = z
    block = 'compacts'
  []
  [max_Tf]
    type = ElementExtremeValue
    variable = thm_temp
    block = '101'
  []
  [P_in]
    type = SideAverageValue
    variable = thm_pressure
    boundary = '300'
  []
  [pressure_drop]
    type = LinearCombinationPostprocessor
    pp_names = 'P_in'
    pp_coefs = '1.0'
    b = '${fparse -outlet_P}'
  []
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    app_type = CardinalApp
    input_files = 'solid.i'
    execute_on = timestep_begin
  []
  [thm]
    type = FullSolveMultiApp
    app_type = CardinalApp
    input_files = 'thm.i'
    execute_on = timestep_end
    max_procs_per_app = 1
    bounding_box_padding = '0.1 0.1 0'
    positions_file = coolant_channel_positions.txt
    output_in_position = true
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [heat_flux_to_openmc]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = flux
    variable = flux
    from_multi_app = bison
    source_boundary = 'fluid_solid_interface'
    target_boundary = 'fluid_solid_interface'
    from_postprocessors_to_be_preserved = flux_integral
    to_postprocessors_to_be_preserved = flux_integral
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    direction = to_multiapp
    multi_app = bison
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = power
  []
  [thm_temp_to_bison]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = thm_temp_wall
    variable = thm_temp
    direction = to_multiapp
    multi_app = bison
  []

  [q_wall_to_thm]
    type = MultiAppGeneralFieldUserObjectTransfer
    variable = q_wall
    direction = to_multiapp
    multi_app = thm
    source_user_object = q_wall_avg
  []
  [T_wall_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = T_wall
    direction = from_multiapp
    multi_app = thm
    variable = thm_temp_wall
    target_boundary = 'fluid_solid_interface'
  []
  [T_bulk_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = T
    direction = from_multiapp
    multi_app = thm
    variable = thm_temp
  []

  # just for postprocessing purposes
  [pressure_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = p
    direction = from_multiapp
    multi_app = thm
    variable = thm_pressure
  []
  [velocity_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = vel_z
    direction = from_multiapp
    multi_app = thm
    variable = thm_velocity
  []
[]

[UserObjects]
  [q_wall_avg]
    type = NearestPointLayeredSideAverage
    boundary = fluid_solid_interface
    variable = flux

    # Note: make this to match the num_elems in the channel
    direction = z
    num_layers = ${num_layers_for_THM}
    points_file = coolant_channel_positions.txt

    direction_min = 0.0
    direction_max = ${height}
  []
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_fluid_axial]
    type = LayeredAverage
    variable = thm_temp
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
  [average_pressure]
    type = LayeredAverage
    variable = thm_pressure
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
  [average_axial_velocity]
    type = LayeredAverage
    variable = thm_velocity
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
  [fluid_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fluid_axial
  []
  [pressure_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_pressure
  []
  [velocity_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_axial_velocity
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Outputs]
  exodus = true
  csv = true
  hide = 'P_in flux_integral z'
[]

(htgr/assembly/openmc.i)

## HTGR assembly neutronics simulation
## Application: Cardinal
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

num_layers_for_THM = 50 # number of elements in the THM model; for the converged
# case, we set this to 150

[GlobalParams]
  # Reduces transfers efficiency for now, can be removed once transferred fields are checked
  bbox_factor = 10
[]

[Mesh]
  # mesh mirror for the solid regions
  [solid]
    type = FileMeshGenerator
    file = meshes/solid_in.e
  []

  # create a mesh for a single coolant channel; because we will receive uniform
  # temperatures and densities from THM on each x-y plane, we can use a very coarse
  # mesh in the radial direction
  [coolant_face]
    type = AnnularMeshGenerator
    nr = 1
    nt = 8
    rmin = 0.0
    rmax = '${fparse channel_diameter / 2.0}'
  []
  [extrude]
    type = AdvancedExtruderGenerator
    input = coolant_face
    num_layers = ${num_layers_for_THM}
    direction = '0 0 1'
    heights = '${height}'
    top_boundary = '300' # inlet
    bottom_boundary = '400' # outlet
  []
  [rename]
    type = RenameBlockGenerator
    input = extrude
    old_block = '1'
    new_block = '101'
  []

  # repeat the coolant channels and then combine together to get a combined mesh mirror
  [repeat]
    type = CombinerGenerator
    inputs = rename
    positions_file = coolant_channel_positions.txt
  []
  [add]
    type = CombinerGenerator
    inputs = 'solid repeat'
  []
[]

[AuxVariables]
  [material_id]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_density]
    family = MONOMIAL
    order = CONSTANT
  []
  [thm_temp_wall]
    block = '101'
  []
  [flux]
  []

  # just for postprocessing purposes
  [thm_pressure]
    block = '101'
  []
  [thm_velocity]
    block = '101'
  []
  [z]
    family = MONOMIAL
    order = CONSTANT
    block = 'compacts'
  []
[]

[AuxKernels]
  [material_id]
    type = CellMaterialIDAux
    variable = material_id
  []
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
  [density]
    type = FluidDensityAux
    variable = density
    p = ${outlet_P}
    T = thm_temp
    fp = helium
    execute_on = 'timestep_begin linear'
  []
  [z]
    type = ParsedAux
    variable = z
    use_xyzt = true
    expression = 'z'
  []
[]

[ICs]
  [fluid_temp_wall]
    type = FunctionIC
    variable = thm_temp_wall
    function = temp_ic
  []
  [fluid_temp]
    type = FunctionIC
    variable = thm_temp
    function = temp_ic
  []
  [heat_source]
    type = ConstantIC
    variable = heat_source
    block = 'compacts'
    value = '${fparse power / (n_bundles * n_fuel_compacts_per_block) / (pi * compact_diameter * compact_diameter / 4.0 * height)}'
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + (${height} - z) / ${height} * ${power} / ${mdot} / ${fluid_Cp}'
  []
[]

[FluidProperties]
  [helium]
    type = IdealGasFluidProperties
    molar_mass = 4e-3
    gamma = 1.668282 # should correspond to  Cp = 5189 J/kg/K
    k = 0.2556
    mu = 3.22639e-5
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  # optimizations to increase tracking rate by assuring that the tallies
  # are spatially independent
  check_tally_sum = false
  normalize_by_global_tally = false
  assume_separate_tallies = true

  power = '${fparse power / n_bundles}'
  scaling = 100.0

  relaxation = constant
  relaxation_factor = 0.5

  inactive_batches = 200
  batches = 1000

  # we will collate temperature from THM (for the fluid) and MOOSE (for the solid)
  # into variables we name as 'solid_temp' and 'thm_temp'. This syntax will automatically
  # create those variabes for us
  temperature_variables = 'solid_temp; solid_temp; solid_temp; thm_temp'
  temperature_blocks = 'graphite; compacts; poison; 101'
  cell_level = 0

  density_blocks = '101'

  [Tallies]
    [heat_source]
      type = CellTally
      block = '2'
      name = heat_source
      score = heating

      trigger = rel_err
      trigger_threshold = 2e-2

      check_equal_mapped_tally_volumes = true
      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Postprocessors]
  [flux_integral]
    type = SideIntegralVariablePostprocessor
    variable = flux
    boundary = 'fluid_solid_interface'
    execute_on = 'transfer linear'
  []
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_rel_err]
    type = TallyRelativeError
    value_type = max
    tally_score = 'heating'
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KEigenvalue
    output = 'std_dev'
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = 'compacts'
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = 'compacts'
  []
  [z_max_power]
    type = ElementExtremeValue
    proxy_variable = heat_source
    variable = z
    block = 'compacts'
  []
  [max_Tf]
    type = ElementExtremeValue
    variable = thm_temp
    block = '101'
  []
  [P_in]
    type = SideAverageValue
    variable = thm_pressure
    boundary = '300'
  []
  [pressure_drop]
    type = LinearCombinationPostprocessor
    pp_names = 'P_in'
    pp_coefs = '1.0'
    b = '${fparse -outlet_P}'
  []
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    app_type = CardinalApp
    input_files = 'solid.i'
    execute_on = timestep_begin
  []
  [thm]
    type = FullSolveMultiApp
    app_type = CardinalApp
    input_files = 'thm.i'
    execute_on = timestep_end
    max_procs_per_app = 1
    bounding_box_padding = '0.1 0.1 0'
    positions_file = coolant_channel_positions.txt
    output_in_position = true
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [heat_flux_to_openmc]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = flux
    variable = flux
    from_multi_app = bison
    source_boundary = 'fluid_solid_interface'
    target_boundary = 'fluid_solid_interface'
    from_postprocessors_to_be_preserved = flux_integral
    to_postprocessors_to_be_preserved = flux_integral
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    direction = to_multiapp
    multi_app = bison
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = power
  []
  [thm_temp_to_bison]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = thm_temp_wall
    variable = thm_temp
    direction = to_multiapp
    multi_app = bison
  []

  [q_wall_to_thm]
    type = MultiAppGeneralFieldUserObjectTransfer
    variable = q_wall
    direction = to_multiapp
    multi_app = thm
    source_user_object = q_wall_avg
  []
  [T_wall_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = T_wall
    direction = from_multiapp
    multi_app = thm
    variable = thm_temp_wall
    target_boundary = 'fluid_solid_interface'
  []
  [T_bulk_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = T
    direction = from_multiapp
    multi_app = thm
    variable = thm_temp
  []

  # just for postprocessing purposes
  [pressure_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = p
    direction = from_multiapp
    multi_app = thm
    variable = thm_pressure
  []
  [velocity_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = vel_z
    direction = from_multiapp
    multi_app = thm
    variable = thm_velocity
  []
[]

[UserObjects]
  [q_wall_avg]
    type = NearestPointLayeredSideAverage
    boundary = fluid_solid_interface
    variable = flux

    # Note: make this to match the num_elems in the channel
    direction = z
    num_layers = ${num_layers_for_THM}
    points_file = coolant_channel_positions.txt

    direction_min = 0.0
    direction_max = ${height}
  []
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_fluid_axial]
    type = LayeredAverage
    variable = thm_temp
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
  [average_pressure]
    type = LayeredAverage
    variable = thm_pressure
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
  [average_axial_velocity]
    type = LayeredAverage
    variable = thm_velocity
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
  [fluid_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fluid_axial
  []
  [pressure_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_pressure
  []
  [velocity_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_axial_velocity
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Outputs]
  exodus = true
  csv = true
  hide = 'P_in flux_integral z'
[]

(htgr/assembly/openmc.i)

## HTGR assembly neutronics simulation
## Application: Cardinal
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

num_layers_for_THM = 50 # number of elements in the THM model; for the converged
# case, we set this to 150

[GlobalParams]
  # Reduces transfers efficiency for now, can be removed once transferred fields are checked
  bbox_factor = 10
[]

[Mesh]
  # mesh mirror for the solid regions
  [solid]
    type = FileMeshGenerator
    file = meshes/solid_in.e
  []

  # create a mesh for a single coolant channel; because we will receive uniform
  # temperatures and densities from THM on each x-y plane, we can use a very coarse
  # mesh in the radial direction
  [coolant_face]
    type = AnnularMeshGenerator
    nr = 1
    nt = 8
    rmin = 0.0
    rmax = '${fparse channel_diameter / 2.0}'
  []
  [extrude]
    type = AdvancedExtruderGenerator
    input = coolant_face
    num_layers = ${num_layers_for_THM}
    direction = '0 0 1'
    heights = '${height}'
    top_boundary = '300' # inlet
    bottom_boundary = '400' # outlet
  []
  [rename]
    type = RenameBlockGenerator
    input = extrude
    old_block = '1'
    new_block = '101'
  []

  # repeat the coolant channels and then combine together to get a combined mesh mirror
  [repeat]
    type = CombinerGenerator
    inputs = rename
    positions_file = coolant_channel_positions.txt
  []
  [add]
    type = CombinerGenerator
    inputs = 'solid repeat'
  []
[]

[AuxVariables]
  [material_id]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_density]
    family = MONOMIAL
    order = CONSTANT
  []
  [thm_temp_wall]
    block = '101'
  []
  [flux]
  []

  # just for postprocessing purposes
  [thm_pressure]
    block = '101'
  []
  [thm_velocity]
    block = '101'
  []
  [z]
    family = MONOMIAL
    order = CONSTANT
    block = 'compacts'
  []
[]

[AuxKernels]
  [material_id]
    type = CellMaterialIDAux
    variable = material_id
  []
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
  [density]
    type = FluidDensityAux
    variable = density
    p = ${outlet_P}
    T = thm_temp
    fp = helium
    execute_on = 'timestep_begin linear'
  []
  [z]
    type = ParsedAux
    variable = z
    use_xyzt = true
    expression = 'z'
  []
[]

[ICs]
  [fluid_temp_wall]
    type = FunctionIC
    variable = thm_temp_wall
    function = temp_ic
  []
  [fluid_temp]
    type = FunctionIC
    variable = thm_temp
    function = temp_ic
  []
  [heat_source]
    type = ConstantIC
    variable = heat_source
    block = 'compacts'
    value = '${fparse power / (n_bundles * n_fuel_compacts_per_block) / (pi * compact_diameter * compact_diameter / 4.0 * height)}'
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + (${height} - z) / ${height} * ${power} / ${mdot} / ${fluid_Cp}'
  []
[]

[FluidProperties]
  [helium]
    type = IdealGasFluidProperties
    molar_mass = 4e-3
    gamma = 1.668282 # should correspond to  Cp = 5189 J/kg/K
    k = 0.2556
    mu = 3.22639e-5
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  # optimizations to increase tracking rate by assuring that the tallies
  # are spatially independent
  check_tally_sum = false
  normalize_by_global_tally = false
  assume_separate_tallies = true

  power = '${fparse power / n_bundles}'
  scaling = 100.0

  relaxation = constant
  relaxation_factor = 0.5

  inactive_batches = 200
  batches = 1000

  # we will collate temperature from THM (for the fluid) and MOOSE (for the solid)
  # into variables we name as 'solid_temp' and 'thm_temp'. This syntax will automatically
  # create those variabes for us
  temperature_variables = 'solid_temp; solid_temp; solid_temp; thm_temp'
  temperature_blocks = 'graphite; compacts; poison; 101'
  cell_level = 0

  density_blocks = '101'

  [Tallies]
    [heat_source]
      type = CellTally
      block = '2'
      name = heat_source
      score = heating

      trigger = rel_err
      trigger_threshold = 2e-2

      check_equal_mapped_tally_volumes = true
      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Postprocessors]
  [flux_integral]
    type = SideIntegralVariablePostprocessor
    variable = flux
    boundary = 'fluid_solid_interface'
    execute_on = 'transfer linear'
  []
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_rel_err]
    type = TallyRelativeError
    value_type = max
    tally_score = 'heating'
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KEigenvalue
    output = 'std_dev'
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = 'compacts'
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = 'compacts'
  []
  [z_max_power]
    type = ElementExtremeValue
    proxy_variable = heat_source
    variable = z
    block = 'compacts'
  []
  [max_Tf]
    type = ElementExtremeValue
    variable = thm_temp
    block = '101'
  []
  [P_in]
    type = SideAverageValue
    variable = thm_pressure
    boundary = '300'
  []
  [pressure_drop]
    type = LinearCombinationPostprocessor
    pp_names = 'P_in'
    pp_coefs = '1.0'
    b = '${fparse -outlet_P}'
  []
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    app_type = CardinalApp
    input_files = 'solid.i'
    execute_on = timestep_begin
  []
  [thm]
    type = FullSolveMultiApp
    app_type = CardinalApp
    input_files = 'thm.i'
    execute_on = timestep_end
    max_procs_per_app = 1
    bounding_box_padding = '0.1 0.1 0'
    positions_file = coolant_channel_positions.txt
    output_in_position = true
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [heat_flux_to_openmc]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = flux
    variable = flux
    from_multi_app = bison
    source_boundary = 'fluid_solid_interface'
    target_boundary = 'fluid_solid_interface'
    from_postprocessors_to_be_preserved = flux_integral
    to_postprocessors_to_be_preserved = flux_integral
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    direction = to_multiapp
    multi_app = bison
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = power
  []
  [thm_temp_to_bison]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = thm_temp_wall
    variable = thm_temp
    direction = to_multiapp
    multi_app = bison
  []

  [q_wall_to_thm]
    type = MultiAppGeneralFieldUserObjectTransfer
    variable = q_wall
    direction = to_multiapp
    multi_app = thm
    source_user_object = q_wall_avg
  []
  [T_wall_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = T_wall
    direction = from_multiapp
    multi_app = thm
    variable = thm_temp_wall
    target_boundary = 'fluid_solid_interface'
  []
  [T_bulk_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = T
    direction = from_multiapp
    multi_app = thm
    variable = thm_temp
  []

  # just for postprocessing purposes
  [pressure_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = p
    direction = from_multiapp
    multi_app = thm
    variable = thm_pressure
  []
  [velocity_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = vel_z
    direction = from_multiapp
    multi_app = thm
    variable = thm_velocity
  []
[]

[UserObjects]
  [q_wall_avg]
    type = NearestPointLayeredSideAverage
    boundary = fluid_solid_interface
    variable = flux

    # Note: make this to match the num_elems in the channel
    direction = z
    num_layers = ${num_layers_for_THM}
    points_file = coolant_channel_positions.txt

    direction_min = 0.0
    direction_max = ${height}
  []
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_fluid_axial]
    type = LayeredAverage
    variable = thm_temp
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
  [average_pressure]
    type = LayeredAverage
    variable = thm_pressure
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
  [average_axial_velocity]
    type = LayeredAverage
    variable = thm_velocity
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
  [fluid_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fluid_axial
  []
  [pressure_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_pressure
  []
  [velocity_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_axial_velocity
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Outputs]
  exodus = true
  csv = true
  hide = 'P_in flux_integral z'
[]

(htgr/assembly/openmc.i)

## HTGR assembly neutronics simulation
## Application: Cardinal
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

num_layers_for_THM = 50 # number of elements in the THM model; for the converged
# case, we set this to 150

[GlobalParams]
  # Reduces transfers efficiency for now, can be removed once transferred fields are checked
  bbox_factor = 10
[]

[Mesh]
  # mesh mirror for the solid regions
  [solid]
    type = FileMeshGenerator
    file = meshes/solid_in.e
  []

  # create a mesh for a single coolant channel; because we will receive uniform
  # temperatures and densities from THM on each x-y plane, we can use a very coarse
  # mesh in the radial direction
  [coolant_face]
    type = AnnularMeshGenerator
    nr = 1
    nt = 8
    rmin = 0.0
    rmax = '${fparse channel_diameter / 2.0}'
  []
  [extrude]
    type = AdvancedExtruderGenerator
    input = coolant_face
    num_layers = ${num_layers_for_THM}
    direction = '0 0 1'
    heights = '${height}'
    top_boundary = '300' # inlet
    bottom_boundary = '400' # outlet
  []
  [rename]
    type = RenameBlockGenerator
    input = extrude
    old_block = '1'
    new_block = '101'
  []

  # repeat the coolant channels and then combine together to get a combined mesh mirror
  [repeat]
    type = CombinerGenerator
    inputs = rename
    positions_file = coolant_channel_positions.txt
  []
  [add]
    type = CombinerGenerator
    inputs = 'solid repeat'
  []
[]

[AuxVariables]
  [material_id]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_density]
    family = MONOMIAL
    order = CONSTANT
  []
  [thm_temp_wall]
    block = '101'
  []
  [flux]
  []

  # just for postprocessing purposes
  [thm_pressure]
    block = '101'
  []
  [thm_velocity]
    block = '101'
  []
  [z]
    family = MONOMIAL
    order = CONSTANT
    block = 'compacts'
  []
[]

[AuxKernels]
  [material_id]
    type = CellMaterialIDAux
    variable = material_id
  []
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
  [density]
    type = FluidDensityAux
    variable = density
    p = ${outlet_P}
    T = thm_temp
    fp = helium
    execute_on = 'timestep_begin linear'
  []
  [z]
    type = ParsedAux
    variable = z
    use_xyzt = true
    expression = 'z'
  []
[]

[ICs]
  [fluid_temp_wall]
    type = FunctionIC
    variable = thm_temp_wall
    function = temp_ic
  []
  [fluid_temp]
    type = FunctionIC
    variable = thm_temp
    function = temp_ic
  []
  [heat_source]
    type = ConstantIC
    variable = heat_source
    block = 'compacts'
    value = '${fparse power / (n_bundles * n_fuel_compacts_per_block) / (pi * compact_diameter * compact_diameter / 4.0 * height)}'
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + (${height} - z) / ${height} * ${power} / ${mdot} / ${fluid_Cp}'
  []
[]

[FluidProperties]
  [helium]
    type = IdealGasFluidProperties
    molar_mass = 4e-3
    gamma = 1.668282 # should correspond to  Cp = 5189 J/kg/K
    k = 0.2556
    mu = 3.22639e-5
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  # optimizations to increase tracking rate by assuring that the tallies
  # are spatially independent
  check_tally_sum = false
  normalize_by_global_tally = false
  assume_separate_tallies = true

  power = '${fparse power / n_bundles}'
  scaling = 100.0

  relaxation = constant
  relaxation_factor = 0.5

  inactive_batches = 200
  batches = 1000

  # we will collate temperature from THM (for the fluid) and MOOSE (for the solid)
  # into variables we name as 'solid_temp' and 'thm_temp'. This syntax will automatically
  # create those variabes for us
  temperature_variables = 'solid_temp; solid_temp; solid_temp; thm_temp'
  temperature_blocks = 'graphite; compacts; poison; 101'
  cell_level = 0

  density_blocks = '101'

  [Tallies]
    [heat_source]
      type = CellTally
      block = '2'
      name = heat_source
      score = heating

      trigger = rel_err
      trigger_threshold = 2e-2

      check_equal_mapped_tally_volumes = true
      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Postprocessors]
  [flux_integral]
    type = SideIntegralVariablePostprocessor
    variable = flux
    boundary = 'fluid_solid_interface'
    execute_on = 'transfer linear'
  []
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_rel_err]
    type = TallyRelativeError
    value_type = max
    tally_score = 'heating'
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KEigenvalue
    output = 'std_dev'
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = 'compacts'
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = 'compacts'
  []
  [z_max_power]
    type = ElementExtremeValue
    proxy_variable = heat_source
    variable = z
    block = 'compacts'
  []
  [max_Tf]
    type = ElementExtremeValue
    variable = thm_temp
    block = '101'
  []
  [P_in]
    type = SideAverageValue
    variable = thm_pressure
    boundary = '300'
  []
  [pressure_drop]
    type = LinearCombinationPostprocessor
    pp_names = 'P_in'
    pp_coefs = '1.0'
    b = '${fparse -outlet_P}'
  []
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    app_type = CardinalApp
    input_files = 'solid.i'
    execute_on = timestep_begin
  []
  [thm]
    type = FullSolveMultiApp
    app_type = CardinalApp
    input_files = 'thm.i'
    execute_on = timestep_end
    max_procs_per_app = 1
    bounding_box_padding = '0.1 0.1 0'
    positions_file = coolant_channel_positions.txt
    output_in_position = true
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [heat_flux_to_openmc]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = flux
    variable = flux
    from_multi_app = bison
    source_boundary = 'fluid_solid_interface'
    target_boundary = 'fluid_solid_interface'
    from_postprocessors_to_be_preserved = flux_integral
    to_postprocessors_to_be_preserved = flux_integral
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    direction = to_multiapp
    multi_app = bison
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = power
  []
  [thm_temp_to_bison]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = thm_temp_wall
    variable = thm_temp
    direction = to_multiapp
    multi_app = bison
  []

  [q_wall_to_thm]
    type = MultiAppGeneralFieldUserObjectTransfer
    variable = q_wall
    direction = to_multiapp
    multi_app = thm
    source_user_object = q_wall_avg
  []
  [T_wall_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = T_wall
    direction = from_multiapp
    multi_app = thm
    variable = thm_temp_wall
    target_boundary = 'fluid_solid_interface'
  []
  [T_bulk_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = T
    direction = from_multiapp
    multi_app = thm
    variable = thm_temp
  []

  # just for postprocessing purposes
  [pressure_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = p
    direction = from_multiapp
    multi_app = thm
    variable = thm_pressure
  []
  [velocity_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = vel_z
    direction = from_multiapp
    multi_app = thm
    variable = thm_velocity
  []
[]

[UserObjects]
  [q_wall_avg]
    type = NearestPointLayeredSideAverage
    boundary = fluid_solid_interface
    variable = flux

    # Note: make this to match the num_elems in the channel
    direction = z
    num_layers = ${num_layers_for_THM}
    points_file = coolant_channel_positions.txt

    direction_min = 0.0
    direction_max = ${height}
  []
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_fluid_axial]
    type = LayeredAverage
    variable = thm_temp
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
  [average_pressure]
    type = LayeredAverage
    variable = thm_pressure
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
  [average_axial_velocity]
    type = LayeredAverage
    variable = thm_velocity
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
  [fluid_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fluid_axial
  []
  [pressure_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_pressure
  []
  [velocity_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_axial_velocity
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Outputs]
  exodus = true
  csv = true
  hide = 'P_in flux_integral z'
[]

(htgr/assembly/openmc.i)

## HTGR assembly neutronics simulation
## Application: Cardinal
## POC: April Novak anovak at anl.gov
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

num_layers_for_THM = 50 # number of elements in the THM model; for the converged
# case, we set this to 150

[GlobalParams]
  # Reduces transfers efficiency for now, can be removed once transferred fields are checked
  bbox_factor = 10
[]

[Mesh]
  # mesh mirror for the solid regions
  [solid]
    type = FileMeshGenerator
    file = meshes/solid_in.e
  []

  # create a mesh for a single coolant channel; because we will receive uniform
  # temperatures and densities from THM on each x-y plane, we can use a very coarse
  # mesh in the radial direction
  [coolant_face]
    type = AnnularMeshGenerator
    nr = 1
    nt = 8
    rmin = 0.0
    rmax = '${fparse channel_diameter / 2.0}'
  []
  [extrude]
    type = AdvancedExtruderGenerator
    input = coolant_face
    num_layers = ${num_layers_for_THM}
    direction = '0 0 1'
    heights = '${height}'
    top_boundary = '300' # inlet
    bottom_boundary = '400' # outlet
  []
  [rename]
    type = RenameBlockGenerator
    input = extrude
    old_block = '1'
    new_block = '101'
  []

  # repeat the coolant channels and then combine together to get a combined mesh mirror
  [repeat]
    type = CombinerGenerator
    inputs = rename
    positions_file = coolant_channel_positions.txt
  []
  [add]
    type = CombinerGenerator
    inputs = 'solid repeat'
  []
[]

[AuxVariables]
  [material_id]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_density]
    family = MONOMIAL
    order = CONSTANT
  []
  [thm_temp_wall]
    block = '101'
  []
  [flux]
  []

  # just for postprocessing purposes
  [thm_pressure]
    block = '101'
  []
  [thm_velocity]
    block = '101'
  []
  [z]
    family = MONOMIAL
    order = CONSTANT
    block = 'compacts'
  []
[]

[AuxKernels]
  [material_id]
    type = CellMaterialIDAux
    variable = material_id
  []
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
  [density]
    type = FluidDensityAux
    variable = density
    p = ${outlet_P}
    T = thm_temp
    fp = helium
    execute_on = 'timestep_begin linear'
  []
  [z]
    type = ParsedAux
    variable = z
    use_xyzt = true
    expression = 'z'
  []
[]

[ICs]
  [fluid_temp_wall]
    type = FunctionIC
    variable = thm_temp_wall
    function = temp_ic
  []
  [fluid_temp]
    type = FunctionIC
    variable = thm_temp
    function = temp_ic
  []
  [heat_source]
    type = ConstantIC
    variable = heat_source
    block = 'compacts'
    value = '${fparse power / (n_bundles * n_fuel_compacts_per_block) / (pi * compact_diameter * compact_diameter / 4.0 * height)}'
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + (${height} - z) / ${height} * ${power} / ${mdot} / ${fluid_Cp}'
  []
[]

[FluidProperties]
  [helium]
    type = IdealGasFluidProperties
    molar_mass = 4e-3
    gamma = 1.668282 # should correspond to  Cp = 5189 J/kg/K
    k = 0.2556
    mu = 3.22639e-5
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  # optimizations to increase tracking rate by assuring that the tallies
  # are spatially independent
  check_tally_sum = false
  normalize_by_global_tally = false
  assume_separate_tallies = true

  power = '${fparse power / n_bundles}'
  scaling = 100.0

  relaxation = constant
  relaxation_factor = 0.5

  inactive_batches = 200
  batches = 1000

  # we will collate temperature from THM (for the fluid) and MOOSE (for the solid)
  # into variables we name as 'solid_temp' and 'thm_temp'. This syntax will automatically
  # create those variabes for us
  temperature_variables = 'solid_temp; solid_temp; solid_temp; thm_temp'
  temperature_blocks = 'graphite; compacts; poison; 101'
  cell_level = 0

  density_blocks = '101'

  [Tallies]
    [heat_source]
      type = CellTally
      block = '2'
      name = heat_source
      score = heating

      trigger = rel_err
      trigger_threshold = 2e-2

      check_equal_mapped_tally_volumes = true
      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Postprocessors]
  [flux_integral]
    type = SideIntegralVariablePostprocessor
    variable = flux
    boundary = 'fluid_solid_interface'
    execute_on = 'transfer linear'
  []
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_rel_err]
    type = TallyRelativeError
    value_type = max
    tally_score = 'heating'
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KEigenvalue
    output = 'std_dev'
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = 'compacts'
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = 'compacts'
  []
  [z_max_power]
    type = ElementExtremeValue
    proxy_variable = heat_source
    variable = z
    block = 'compacts'
  []
  [max_Tf]
    type = ElementExtremeValue
    variable = thm_temp
    block = '101'
  []
  [P_in]
    type = SideAverageValue
    variable = thm_pressure
    boundary = '300'
  []
  [pressure_drop]
    type = LinearCombinationPostprocessor
    pp_names = 'P_in'
    pp_coefs = '1.0'
    b = '${fparse -outlet_P}'
  []
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    app_type = CardinalApp
    input_files = 'solid.i'
    execute_on = timestep_begin
  []
  [thm]
    type = FullSolveMultiApp
    app_type = CardinalApp
    input_files = 'thm.i'
    execute_on = timestep_end
    max_procs_per_app = 1
    bounding_box_padding = '0.1 0.1 0'
    positions_file = coolant_channel_positions.txt
    output_in_position = true
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [heat_flux_to_openmc]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = flux
    variable = flux
    from_multi_app = bison
    source_boundary = 'fluid_solid_interface'
    target_boundary = 'fluid_solid_interface'
    from_postprocessors_to_be_preserved = flux_integral
    to_postprocessors_to_be_preserved = flux_integral
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    direction = to_multiapp
    multi_app = bison
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = power
  []
  [thm_temp_to_bison]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = thm_temp_wall
    variable = thm_temp
    direction = to_multiapp
    multi_app = bison
  []

  [q_wall_to_thm]
    type = MultiAppGeneralFieldUserObjectTransfer
    variable = q_wall
    direction = to_multiapp
    multi_app = thm
    source_user_object = q_wall_avg
  []
  [T_wall_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = T_wall
    direction = from_multiapp
    multi_app = thm
    variable = thm_temp_wall
    target_boundary = 'fluid_solid_interface'
  []
  [T_bulk_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = T
    direction = from_multiapp
    multi_app = thm
    variable = thm_temp
  []

  # just for postprocessing purposes
  [pressure_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = p
    direction = from_multiapp
    multi_app = thm
    variable = thm_pressure
  []
  [velocity_from_thm]
    type = MultiAppGeneralFieldNearestNodeTransfer
    source_variable = vel_z
    direction = from_multiapp
    multi_app = thm
    variable = thm_velocity
  []
[]

[UserObjects]
  [q_wall_avg]
    type = NearestPointLayeredSideAverage
    boundary = fluid_solid_interface
    variable = flux

    # Note: make this to match the num_elems in the channel
    direction = z
    num_layers = ${num_layers_for_THM}
    points_file = coolant_channel_positions.txt

    direction_min = 0.0
    direction_max = ${height}
  []
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_fluid_axial]
    type = LayeredAverage
    variable = thm_temp
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
  [average_pressure]
    type = LayeredAverage
    variable = thm_pressure
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
  [average_axial_velocity]
    type = LayeredAverage
    variable = thm_velocity
    direction = z
    num_layers = ${num_layers_for_plots}
    block = '101'
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
  [fluid_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fluid_axial
  []
  [pressure_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_pressure
  []
  [velocity_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_axial_velocity
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Outputs]
  exodus = true
  csv = true
  hide = 'P_in flux_integral z'
[]

Geometry and Computational Model
Multiphysics Coupling
Execution and Postprocessing
References