Multiphysics for a TRISO Gas-Cooled Compact

In this tutorial, you will learn how to:

Couple OpenMC, NekRS/THM, and MOOSE together for multiphysics modeling of a TRISO compact
Use two different MultiApp hierarchies to achieve different data transfers
Use triggers to automatically terminate the OpenMC active batches once reaching the desired statistical uncertainty
Automatically detect steady state

To access this tutorial,

cd cardinal/tutorials/gas_compact_multiphysics

This tutorial also requires you to download mesh files and a NekRS restart file from Box. Please download the files from the gas_compact_multiphysics folder here and place these within the same directory structure in tutorials/gas_compact_multiphysics.

note:Computing Needs

This tutorial requires HPC resources to run the NekRS cases. You will be able to run the OpenMC-THM-MOOSE files without any special resources.

In this tutorial, we couple OpenMC to the MOOSE heat transfer module, with fluid feedback provided by either NekRS or THM,

NekRS: wall-resolved $k$ - $\tau$ Reynolds Averaged Navier-Stokes (RANS) equations
Thermal Hydraulics Module (THM): 1-D area-averaged Navier-Stokes equations

Two different multiapp hierarchies will be used in order to demonstrate the flexibility of the MultiApp system. The same OpenMC model can be used to provide feedback to different combinations of MOOSE applications.

In this tutorial, OpenMC receives temperature feedback from the MOOSE heat conduction module (for the solid regions) and NekRS/THM (for the fluid regions). Density feedback is provided by NekRS/THM for the fluid regions. This tutorial models a partial-height TRISO-fueled unit cell of a prismatic gas reactor assembly, and is a continuation of the conjugate heat transfer tutorial (where we coupled NekRS and MOOSE heat conduction) and the OpenMC-heat conduction tutorial (where we coupled OpenMC and MOOSE heat conduction) for this geometry.

This tutorial was developed with support from the NEAMS Thermal Fluids Center of Excellence and is described in more detail in our journal article Novak et al. (2022).

Geometry and Computational Model

The geometry consists of a TRISO-fueled gas reactor compact unit cell INL (2016). A top-down view of the geometry is shown in Figure 1. The fuel is cooled by helium flowing in a cylindrical channel of diameter $d_c$ . Cylindrical fuel compacts containing randomly-dispersed TRISO particles at 15% packing fraction are arranged around the coolant channel in a triangular lattice. The TRISO particles use a conventional design that consists of a central fissile uranium oxycarbide kernel enclosed in a carbon buffer, an inner PyC layer, a silicon carbide layer, and finally an outer PyC layer. The geometric specifications are summarized in Table 1. Heat is produced in the TRISO particles to yield a total power of 38 kW.

Figure 1: TRISO-fueled gas reactor compact unit cell

Table 1: Geometric specifications for a TRISO-fueled gas reactor compact

Parameter	Value (cm)
Coolant channel diameter, $d_c$	1.6
Fuel compact diameter, $d_f$	1.27
Fuel-to-coolant center distance, $p_{cf}$	1.88
Height	160
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

Two different MultiApp hierarchies are used:

A "single stack" design where each application has either a single "parent" application or a single "child" application
A "tree" design where each application has a single "parent" application, but multiple "child" applications

Figure 2 shows a conceptual depiction of these hierarchies. We will describe these in greater detail later, but introduce them here to assist with describing a few aspects of the single-physics models. The circled numbers indicate the order in which the applications run.

Solid lines depict transfers that occur directly from application $A$ to application $B$ , or between the source and receiver of that field. Dashed lines, on the other hand, depict transfers that do not occur directly between the source and receiver of the field - for instance, in the "single stack" hierarchy, the NekRS application can only communicate data with it's immediate parent application. Therefore, to send fluid density and temperature from NekRS to OpenMC, there are actually two transfers - 1) sending fluid density and temperature from NekRS to the MOOSE heat transfer module, and 2) sending fluid density and temperature from the MOOSE heat transfer module to OpenMC.

Figure 2: MultiApp hierarchies used in this tutorial; data transfers are shown with solid and dashed lines. Solid lines indicate transfers that occur directly from application A to application B, while dashed lines show transfers that have to first pass through an intermediate application to get to the eventual target application.

Conversely, with the "tree" hierarhcy, MOOSE MultiApps communicate with parent/child applications. Therefore, all data communicated between the MOOSE heat transfer module and THM actually has to first pass through their common parent application before reaching the desired target application.

note

In the time since we originally developed this tutorial, MOOSE has been extended to support sibling transfers, which would allow MOOSE and THM to communicate data directly to one another (in the "tree" hierarchy shown in Figure 2).

OpenMC Model

The OpenMC model is built using CSG. The TRISO positions are sampled using the RSA algorithm in OpenMC. OpenMC's Python API is used to create the model with the script shown below. First, we define materials. 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 pack uniform-radius spheres into a cylindrical region representing a fuel compact, setting each sphere to be filled with the TRISO universe.

#!/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):
    """
    THIS IS ONLY USED FOR SETTING AN INITIAL CONDITION IN OPENMC's XML FILES -
    the coolant temperature will be applied from MOOSE, we just set an initial
    value here in case you want to run these files in standalone mode (i.e. with
    the "openmc" executable).

    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
    distance_from_inlet = specs.unit_cell_height - z
    Qi = (l - l * np.cos(math.pi * distance_from_inlet / l)) / math.pi

    t = t_in + Qi / Q * dT

    return t

def coolant_density(t):
  """
    THIS IS ONLY USED FOR SETTING AN INITIAL CONDITION IN OPENMC's XML FILES -
    the coolant density will be applied from MOOSE, we just set an initial
    value here in case you want to run these files in standalone mode (i.e. with
    the "openmc" executable).

  Computes the helium density (kg/m3) 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
# -------------------------------------------

### RADIANT UNIT CELL SPECS (INFERRED FROM REPORTS) ###

unit_cell_mdot = specs.mdot / (specs.n_bundles * specs.n_coolant_channels_per_block)
unit_cell_power = specs.power / (specs.n_bundles * specs.n_coolant_channels_per_block) * (specs.unit_cell_height / specs.height)

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

# geometry
coolant_channel_diam = specs.channel_diameter * m
reactor_bottom = 0.0
reactor_height = specs.unit_cell_height * m
reactor_top = reactor_bottom + reactor_height

### ARBITRARILY DETERMINED PARAMETERS ###

cell_pitch = specs.fuel_to_coolant_distance * m
fuel_channel_diam = specs.compact_diameter * m

hex_orientation = 'x'

def unit_cell(n_ax_zones, n_inactive, n_active):
    axial_section_height = reactor_height / n_ax_zones

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

    lattice_orientation = 'x'
    cell_edge_length = cell_pitch

    model = openmc.model.Model()

    ### Geometry ###

    # 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)

    # 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, seed=1.0)
    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)

    # create a hexagonal lattice for the coolant and fuel channels
    fuel_univ = openmc.Universe(cells=[openmc.Cell(region=-fuel_cyl, fill=triso_lattice),
                                    openmc.Cell(region=+fuel_cyl, fill=mats.m_graphite_matrix)])

    # extract the coolant cell and set temperatures based on the axial profile
    coolant_cell = openmc.Cell(region=-coolant_cyl, fill=mats.m_coolant)
    # set the coolant temperature on the cell to approximately match the expected
    # temperature profile
    axial_coords = np.linspace(reactor_bottom, reactor_top, n_ax_zones + 1)
    lattice_univs = []

    fuel_ch_cells = []

    i = 0
    for z_min, z_max in zip(axial_coords[0:-1], axial_coords[1:]):
        # create a new coolant universe for each axial zone in the coolant channel;
        # this generates a new material as well (we only need to do this for all
        # cells except the first cell)
        if (i == 0):
          c_cell = coolant_cell
        else:
          c_cell = coolant_cell.clone()

        i += 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)

        c_cell.temperature = t
        c_cell.fill.set_density('kg/m3', coolant_density(t))

        # set the solid cells and their temperatures
        graphite_cell = openmc.Cell(region=+coolant_cyl, fill=mats.m_graphite_matrix)
        fuel_ch_cell = openmc.Cell(region=-fuel_cyl, fill=triso_lattice)
        fuel_ch_matrix_cell = openmc.Cell(region=+fuel_cyl, fill=mats.m_graphite_matrix)

        graphite_cell.temperature = t
        fuel_ch_cell.temperature = t
        fuel_ch_matrix_cell.temperature = t

        fuel_ch_cells.append(fuel_ch_cell)
        fuel_u = openmc.Universe(cells=[fuel_ch_cell, fuel_ch_matrix_cell])
        coolant_u = openmc.Universe(cells=[c_cell, graphite_cell])
        lattice_univs.append([[fuel_u] * 6, [coolant_u]])

    # create a hexagonal lattice used in each axial zone to represent the cell
    hex_lattice = openmc.HexLattice(name="Unit cell lattice")
    hex_lattice.orientation = lattice_orientation
    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

    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

    # hexagonal bounding cell
    hex = openmc.model.HexagonalPrism(cell_edge_length, hex_orientation, boundary_type='periodic')

    hex_cell_vol = 6.0 * (math.sqrt(3) / 4.0) * cell_edge_length**2 * reactor_height

    # create additional axial regions
    axial_planes = [openmc.ZPlane(z0=coord) for coord in axial_coords]
    # axial planes
    min_z = axial_planes[0]
    min_z.boundary_type = 'vacuum'
    max_z = axial_planes[-1]
    max_z.boundary_type = 'vacuum'

    # fill the unit cell with the hex lattice
    hex_cell = openmc.Cell(region=-hex & +min_z & -max_z, fill=hex_lattice)

    model.geometry = openmc.Geometry([hex_cell])

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

    settings.particles = 150000
    settings.inactive = n_inactive
    settings.batches = settings.inactive + n_active

    # the only reason we use 'nearest' here is to be sure we have a robust test for CI;
    # otherwise, 1e-16 differences in temperature (due to numerical roundoff when using
    # different MPI ranks) do change the tracking do to the stochastic interpolation
    settings.temperature['method'] = 'nearest'
    settings.temperature['range'] = (294.0, 1500.0)
    settings.temperature['tolerance'] = 200.0

    hexagon_half_flat = math.sqrt(3.0) / 2.0 * cell_edge_length
    lower_left = (-cell_edge_length, -hexagon_half_flat, reactor_bottom)
    upper_right = (cell_edge_length, hexagon_half_flat, reactor_top)
    source_dist = openmc.stats.Box(lower_left, upper_right)
    source = openmc.IndependentSource(space=source_dist)
    settings.source = source

    model.settings = settings

    m_colors = {}
    m_colors[mats.m_coolant] = 'royalblue'
    m_colors[mats.m_fuel] = 'red'
    m_colors[mats.m_graphite_c_buffer] = 'black'
    m_colors[mats.m_graphite_pyc] = 'orange'
    m_colors[mats.m_sic] = 'yellow'
    m_colors[mats.m_graphite_matrix] = 'silver'

    plot1          = openmc.Plot()
    plot1.filename = 'plot1'
    plot1.width    = (2 * cell_pitch, 4 * axial_section_height)
    plot1.basis    = 'xz'
    plot1.origin   = (0.0, 0.0, reactor_height/2.0)
    plot1.pixels   = (int(800 * 2 * cell_pitch), int(800 * 4 * axial_section_height))
    plot1.color_by = 'cell'

    plot2          = openmc.Plot()
    plot2.filename = 'plot2'
    plot2.width    = (3 * cell_pitch, 3 * cell_pitch)
    plot2.basis    = 'xy'
    plot2.origin   = (0.0, 0.0, axial_section_height / 2.0)
    plot2.pixels   = (int(800 * cell_pitch), int(800 * cell_pitch))
    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('-i', dest='n_inactive', type=int, default=20,
                    help='Number of inactive cycles')
    ap.add_argument('-a', dest='n_active', type=int, default=45,
                    help='Number of active cycles')

    args = ap.parse_args()

    model = unit_cell(args.n_axial, args.n_inactive, args.n_active)
    model.export_to_xml()

if __name__ == "__main__":
    main()

Finally, we loop over $n_l$ axial layers and create unique cells for each of the six compacts, the graphite block, and the coolant. Recall that we need unique cells in order for each region to obtain a a unique temperature from MOOSE. The level on which we will apply feedback from MOOSE is set to 1 because each layer is a component in a lattice nested once with respect to the highest level. 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 cell ID, is shown in Figure 3. The lateral faces of the unit cell are periodic, while the top and bottom boundaries are vacuum. The Cartesian search lattice in the fuel compact regions is also visible.

Figure 3: OpenMC model, colored by cell ID

For the "single-stack" MultiApp hierarchy, OpenMC runs first, so the initial temperature is set to uniform in the radial direction and given by a linear variation between the inlet and outlet fluid temperatures. The fluid density is then set using the ideal gas Equation of State (EOS) with pressure taken as the fixed outlet of 7.1 MPa given the temperature, i.e. $\rho_f(P, T)$ . For the "tree" MultiApp hierarchy, OpenMC instead runs after the MOOSE heat transfer module, but before THM. For this structure, initial conditions are only required for fluid temperature and density, which are taken as the same initial conditions as for the "single-stack" case.

To create the XML files required to run OpenMC, run the script:

python unit_cell.py

You can also use the XML files checked in to the tutorials/gas_compact_multiphysics directory.

Heat Conduction Model

The MOOSE heat transfer module is used to solve for energy conservation in the solid. The solid mesh is shown in Figure 4; the only sideset defined in the domain is the coolant channel surface. The solid geometry uses a length unit of meters.

Figure 4: Mesh for the solid heat conduction model

This mesh is generated using MOOSE mesh generators in the solid_mesh.i file.

[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 = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    ring_block_ids = '2 2'
    ring_block_names = 'compacts compacts'
    background_block_names = 'graphite'
    background_intervals = 4
  []
  [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 = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    ring_block_ids = '101 101'
    ring_block_names = 'coolant coolant'
    background_block_names = 'graphite'
    interface_boundary_id_shift = 100
    background_intervals = 3
  []
  [bundle]
    type = PatternedHexMeshGenerator
    inputs = 'fuel_pin coolant_pin'
    hexagon_size = '${fparse 2.0 * fuel_to_coolant_distance}'
    pattern = '0 0;
              0 1 0;
               0 0'
  []
  [trim]
    type = HexagonMeshTrimmer
    input = bundle
    trim_peripheral_region = '1 1 1 1 1 1'
    peripheral_trimming_section_boundary = peripheral_section
  []
  [rotate]
    type = TransformGenerator
    input = trim
    transform = rotate
    vector_value = '30.0 0.0 0.0'
  []
  [extrude]
    type = AdvancedExtruderGenerator
    input = rotate
    heights = ${height}
    num_layers = ${n_layers}
    direction = '0 0 1'
  []
  [fluid_solid_interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = extrude
    primary_block = 'graphite'
    paired_block = 'coolant'
    new_boundary = 'fluid_solid_interface'
  []
[]

We first create a full 7-pin bundle, and then apply a trimming operation to split the compacts. Because MOOSE does not support multiple element types (e.g. tets, hexes) on the same block ID, the trimmer automatically creates an additional block (compacts_trimmer_tri) to represent the triangular prism elements formed in the compacts. Note that within this mesh, we include the fluid region - for the "single stack" MultiApp hierarchy, we will need somewhere for NekRS to write the fluid temperature solution. So, while this block does not participate in the solid solve, we include it in the mesh just for data transfers. You can generate this mesh by running

cardinal-opt -i solid_mesh.i --mesh-only

which will create the mesh, named solid_mesh_in.e.

On the coolant channel surface, a Dirichlet temperature is provided by NekRS/THM. All other boundaries are insulated. The volumetric power density is provided by OpenMC, with normalization to ensure the total specified power. When using the "single stack" hierarchy, MOOSE runs after OpenMC but before NekRS, and an initial condition is only required for the wall temperature, which is set to a linear variation from inlet to outlet fluid temperature. When using the "tree" hierarchy, MOOSE runs first, in which case the initial wall temperature is taken as the same linear variation, while the power is taken as uniform.

NekRS Model

NekRS is used to solve the incompressible k-tau RANS model. The inlet mass flowrate is 0.0905 kg/s; with the channel diameter of 1.6 cm and material properties of helium, this results in a Reynolds number of 223214 and a Prandtl number of 0.655. This highly-turbulent flow results in extremely thin momentum and thermal boundary layers on the no-slip surfaces forming the periphery of the coolant channel. In order to resolve the near-wall behavior with a wall-resolved model, an extremely fine mesh is required in the NekRS simulation. To accelerate the overall coupled solve that is of interest in this tutorial, the NekRS model is split into a series of calculations:

We first run a partial-height, periodic flow-only case to obtain converged $P$ , $\vec{u}$ , and $\mu_T$ distributions.
Then, we extrapolate the $\vec{u}$ and $\mu_T$ to the full-height case.
We use the converged, full-height $\vec{u}$ and $\mu_T$ distributions to transport a temperature passive scalar in a Conjugate Heat Transfer (CHT) calculation with MOOSE.
Finally, we use the converged CHT case as an initial condition for the multiphysics simulation with OpenMC and MOOSE feedback.

Steps 1-3 were performed in an earlier tutorial - for brevity, we skip repeating the discussion of steps 1-3.

For the multiphysics case, we will load the restart file produced from step 3, compute $k_T$ from the loaded solutions for $k$ and $\tau$ , and then transport temperature with coupling to MOOSE heat conduction and OpenMC particle transport. Let's now describe the NekRS input files needed for the passive scalar solve:

ranstube.re2: NekRS mesh
ranstube.par: High-level settings for the solver, boundary condition mappings to sidesets, and the equations to solve
ranstube.udf: User-defined C++ functions for on-line postprocessing and model setup
ranstube.oudf: User-defined Open Concurrent Compute Abstraction (OCCA) kernels for boundary conditions and source terms

A detailed description of all of the available parameters, settings, and use cases for these input files is available on the NekRS documentation website. Because the purpose of this analysis is to demonstrate Cardinal's capabilities, only the aspects of NekRS required to understand the present case will be covered. First, the NekRS mesh is shown in Figure 5. Boundary 1 is the inlet, boundary 2 is the outlet, and boundary 3 is the wall. The same mesh was used for the periodic flow solve, except with a shorter height.

Figure 5: Mesh for the NekRS RANS model

Next, the .par file contains problem setup information. This input sets up a nondimensional passive scalar solution, loading $P$ , $\vec{u}$ , $k$ , and $\tau$ from a restart file. We "freeze" the flow by setting solver = none in the [VELOCITY], [SCALAR01] ( $k$ passive scalar), and [SCALAR02] ( $\tau$ passive scalar) blocks. In the nondimensional formulation, the "viscosity" becomes $1/Re$ , where $Re$ is the Reynolds number, while the "thermal conductivity" becomes $1/Pe$ , where $Pe$ is the Peclet number. These nondimensional numbers are used to set various diffusion coefficients in the governing equations with syntax like -223214, which is equivalent in NekRS syntax to $\frac{1}{223214}$ . The only equation that NekRS will solve is for temperature.

# Model of a turbulent channel flow.
#
# L_ref = 0.016 m          coolant channel diameter
# T_ref = 598 K            coolant inlet temperature
# dT_ref = 82.89 K         coolant nominal temperature rise in the unit cell
# U_ref = 81.089           coolant inlet velocity
# length = 1.6 m           coolant channel flow length
#
# rho_0 = 5.5508 kg/m3     coolant density
# mu_0 = 3.22639e-5 Pa-s   coolant dynamic viscosity
# Cp_0 = 5189 J/kg/K       coolant isobaric heat capacity
# k_0 = 0.2556 W/m/K       coolant thermal conductivity
#
# Re = 223214
# Pr = 0.655
# Pe = 146205

[GENERAL]
  startFrom = converged_cht.fld

  polynomialOrder = 7
  dt = 6e-3
  timeStepper = BDF2

  # end the simulation at specified number of time steps
  stopAt = numSteps
  numSteps = 1000

  # write an output file every specified time steps
  writeControl = steps
  writeInterval = 1000

[PROBLEMTYPE]
  stressFormulation = true

[PRESSURE]
  residualTol = 1e-04

[VELOCITY]
  solver = none
  viscosity = -223214
  density = 1.0
  boundaryTypeMap = v, o, W
  residualTol = 1e-06

[TEMPERATURE]
  rhoCp = 1.0
  conductivity = -146205
  boundaryTypeMap = t, o, f
  residualTol = 1e-06

[SCALAR01] # k
  solver = none
  boundaryTypeMap = t, I, o
  residualTol = 1e-06
  rho = 1.0
  diffusivity = -223214

[SCALAR02] # tau
  solver = none
  boundaryTypeMap = t, I, o
  residualTol = 1e-06
  rho = 1.0
  diffusivity = -223214

Next, the .udf file is used to setup initial conditions and define how $k_T$ should be computed based on $Pr_T$ and the restart values of $k$ and $\tau$ . In turbulent_props, a user-defined function, we use $k_f$ from the input file in combination with the $Pr_T$ and $\mu_T$ (read from the restart file later in the .udf file) to adjust the total diffusion coefficient on temperature to $k_f+k_T$ according to the turbulent Prandtl number definition. This adjustment must happen on device, in a new GPU kernel we name scalarScaledAddKernel. This kernel will be defined in the .oudf file; we instruct the JIT compilation to compile this new kernel by calling udfBuildKernel.

Then, in UDF_Setup we store the value of $\mu_T$ computed in the restart file.

#include <math.h>
#include "udf.hpp"

occa::memory o_mu_t_from_restart;

#define length 100.0    // non-dimensional channel length
#define Pr_T 0.91       // turbulent Prandtl number

static occa::kernel scalarScaledAddKernel;

void turbulent_props(nrs_t *nrs, double time, occa::memory o_U, occa::memory o_S,
            occa::memory o_UProp, occa::memory o_SProp)
{
  // fetch the laminar thermal conductivity and specify the desired turbulent Pr
  dfloat k_laminar;
  platform->options.getArgs("SCALAR00 DIFFUSIVITY", k_laminar);

  // o_diff is the turbulent conductivity for all passive scalars, so we grab the first
  // slice, since temperature is the first passive scalar
  occa::memory o_mu = nrs->cds->o_diff + 0 * nrs->cds->fieldOffset[0] * sizeof(dfloat);

  scalarScaledAddKernel(0, nrs->cds->mesh[0]->Nlocal, k_laminar, 1.0 / Pr_T,
    o_mu_t_from_restart /* turbulent viscosity */, o_mu /* laminar viscosity */);
}

void UDF_LoadKernels(occa::properties & kernelInfo)
{
  scalarScaledAddKernel = oudfBuildKernel(kernelInfo, "scalarScaledAdd");
}

void UDF_Setup(nrs_t *nrs)
{
  mesh_t *mesh = nrs->cds->mesh[0];
  int n_gll_points = mesh->Np * mesh->Nelements;

  // allocate space for the k and tau that we will read from the file
  dfloat * mu_t = (dfloat *) calloc(n_gll_points, sizeof(dfloat));

  for (int n = 0; n < n_gll_points; ++n)
  {
    // fetch the restart value of mu_T, which is equal to k * tau, the second and third
    // passive scalars
    mu_t[n] = nrs->cds->S[n + 1 * nrs->cds->fieldOffset[0]] *
              nrs->cds->S[n + 2 * nrs->cds->fieldOffset[0]];
  }

  o_mu_t_from_restart = platform->device.malloc(n_gll_points * sizeof(dfloat), mu_t);

  udf.properties = &turbulent_props;
}

void UDF_ExecuteStep(nrs_t * nrs, double time, int tstep)
{
}

In the .oudf file, we define boundary conditions for temperature and also the form of the scalarScaledAdd kernel that we use to compute $k_T$ . The inlet boundary is set to a temperature of 0 (a dimensional temperature of $T_{ref}$ ), while the fluid-solid interface will receive a heat flux from MOOSE.

@kernel void scalarScaledAdd(const dlong offset, const dlong N,
                             const dfloat a,
                             const dfloat b,
                             @restrict const dfloat* X,
                             @restrict dfloat* Y)
{
  for (dlong n = 0; n < N; ++n; @tile(256,@outer,@inner))
    if (n < N)
      Y[n] = a + b * X[n + offset];
}

// Boundary conditions are only needed for temperature, since the solves of pressure,
// velocity, k, and tau are all frozen

void scalarDirichletConditions(bcData *bc)
{
  // inlet temperature
  bc->s = 0.0;
}

void scalarNeumannConditions(bcData *bc)
{
  // wall temperature
  // note: when running with Cardinal, Cardinal will allocate the usrwrk
  // array. If running with NekRS standalone (e.g. nrsmpi), you need to
  // replace the usrwrk with some other value or allocate it youself from
  // the udf and populate it with values.
  bc->flux = bc->usrwrk[bc->idM];
}

For this tutorial, NekRS runs last in the "single-stack" MultiApp hierarchy, so no initial conditions are required aside from the $T$ , $\vec{u}$ , and $\mu_T$ taken from the converged_cht.fld restart file on Box.

THM Model

THM is used to solve the 1-D area-averaged Navier-Stokes equations. The THM mesh contains 150 elements; the mesh is constucted automatically within THM. 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 150 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.

For this tutorial, THM runs last in the "tree" MultiApp hierarchy; because THM solves time-dependent equations, initial conditions are only required for the solution variables for which THM solves - pressure, fluid temperature, and velocity, all of which are set to uniform conditions.

Multiphysics Coupling

In this section, OpenMC, NekRS/THM, and MOOSE heat conduction are coupled for multiphysics modeling of the TRISO gas compact. Two separate simulations are performed here:

Coupling of OpenMC, NekRS, and MOOSE heat conduction in a "single-stack" MultiApp hierarchy
Coupling of OpenMC, THM, and MOOSE heat conduction in a "tree" MultiApp hierarchy

By individually describing the two setups, you will understand the customizability of the MultiApp system and the flexibility shared by all MOOSE applications for seamlessly exchanging tools of varying resolution for one another.

OpenMC-NekRS-MOOSE

In this section, we describe the coupling of OpenMC, NekRS, and MOOSE in the "single-stack" MultiApp hierarchy shown in Figure 2.

OpenMC Input Files

The neutronics physics is solved over the entire domain with OpenMC. The OpenMC wrapping used for the OpenMC-NekRS-MOOSE coupling is described in the openmc_nek.i input 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 THM-MOOSE, and on which OpenMC will write the fission heat source. Because the coupled applications use length units of meters, the mesh mirror must also be in units of meters. For simplicity, we just use the same mesh as generated with solid_mesh.i earlier, though this it not necessary.

!include common_input.i

# This input file runs coupled OpenMC Monte Carlo transport, MOOSE heat
# conduction, and NekRS fluid flow and heat transfer.
# This input should be run with:
#
# cardinal-opt -i openmc_nek.i

num_layers_for_THM = 150
density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3
N = 1000

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

Next, we define a number of auxiliary variables to query the OpenMC solution and set up multiphysics coupling. We add cell_temperature and cell_density in order to read the cell temperatures and densities directly from OpenMC in order to visualize the temperatures and densities ultimately applied to OpenMC's cells. In order to compute density using the ideal gas EOS given a temperature and a fixed pressure, we use a FluidDensityAux to set density.

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

[AuxKernels]
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
  [density]
    type = FluidDensityAux
    variable = density
    p = ${outlet_P}
    T = nek_temp
    fp = helium
    execute_on = 'timestep_begin'
  []
[]

[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 necessary initial conditions using functions; all temperatures (nek_temp and solid_temp are to be discussed shortly) are set to a linear variation from the inlet to outlet fluid temperatures.

[ICs]
  [temp]
    type = FunctionIC
    variable = nek_temp
    function = temp_ic
  []
  [solid_temp]
    type = FunctionIC
    variable = solid_temp
    function = temp_ic
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / (${mdot} / ${n_bundles} / ${n_coolant_channels_per_block}) / ${fluid_Cp}'
  []
[]

The wrapping of OpenMC is specified in the [Problem] block and the addition of tallies is done in the [Tallies] block. Here, we indicate that we will provide both temperature and density feedback to OpenMC. In order to visualize the tally standard deviation, we output the fission tally standard deviation using the output parameter. The heat source from OpenMC will be relaxed using Robbins-Monro relaxation.

By default, OpenMC will try to read temperature from the temp variable. However, in this case we have multiple applications (NekRS/THM for the fluids, MOOSE for the solids) which want to send temperatures into OpenMC. To be sure one data transfer does not overwrite the second, we need to tell OpenMC the names of the variables to read temperature from. Cardinal contains convenient syntax to automatically set up the necessary receiver variables for temperature, by using the temperature_variables and temperature_blocks parameters.

Finally, a number of "triggers" are used to automatically terminate OpenMC's active batches once reaching desired uncertainties in $k$ and the fission tally. The number of batches here is terminated once both of the following are satisfied:

Standard deviation in $k$ is less than 75 pcm
Maximum fission tally relative error is less than 1%

These criteria are checked every batch_interval, up to a maximum number of batches.

[Problem]
  type = OpenMCCellAverageProblem

  power = ${unit_cell_power}
  scaling = 100.0
  density_blocks = ${density_blocks}
  cell_level = 1

  relaxation = robbins_monro

  temperature_variables = 'solid_temp;   nek_temp'
  temperature_blocks = '${temperature_blocks}; ${density_blocks}'

  k_trigger = std_dev
  k_trigger_threshold = 7.5e-4
  batches = 40
  max_batches = 100
  batch_interval = 5

  [Tallies]
    [heat_source]
      type = CellTally
      blocks = ${fuel_blocks}
      name = heat_source

      check_equal_mapped_tally_volumes = true

      trigger = rel_err
      trigger_threshold = 1e-2

      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

Next, we define a transient executioner - while OpenMC is technically solving a steady $k$ -eigenvalue problem, using a time-dependent executioner with the notion of a "time step" will allow us to control the frequency with which OpenMC sends data to/from it's sub-app (MOOSE heat conduction). Here, we set the time step to be 1000 times the NekRS fluid time step. We will terminate the solution once reaching the specified steady_state_tolerance, which by setting check_aux = true will terminate the entire solution once reaching less than 1% change in all auxiliary variables in the OpenMC input file (which includes temperatures and power).

[Executioner]
  type = Transient
  dt = '${fparse N * nek_dt * t0}'

  # This is somewhat loose, and you should use a tighter tolerance for production runs;
  # we use 1e-2 to obtain a faster tutorial
  steady_state_detection = true
  check_aux = true
  steady_state_tolerance = 1e-2
[]

Next, we define the MOOSE heat conduction sub-application and data transfers to/from that application. Most important to note is that while the MOOSE heat transfer module does not itself compute fluid temperature, we see a transfer getting the fluid temperature that has been transferred to the MOOSE heat transfer module by the doubly-nested NekRS sub-application.

[MultiApps]
  [bison]
    type = TransientMultiApp
    input_files = 'solid_nek.i'
    execute_on = timestep_end
    sub_cycling = true
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    to_multi_app = bison
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = power
  []
  [temp_from_nek]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = nek_bulk_temp
    from_multi_app = bison
    variable = nek_temp
  []
[]

Next, we define several postprocessors for querying the solution. The heat_source postprocessor will be used to ensure conservation of power when sent to the MOOSE heat conduction application. All other postprocessors are used for general solution monitoring purposes.

[Postprocessors]
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_err]
    type = TallyRelativeError
    value_type = max
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KStandardDeviation
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = ${fuel_blocks}
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = ${fuel_blocks}
  []
[]

For postprocessing, we also compute the average power distribution in a number of layers using a LayeredAverage and output to CSV using a SpatialUserObjectVectorPostprocessor, in combination with a CSV output.

[UserObjects]
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers_for_plots}
    block = ${fuel_blocks}
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
[]

[Outputs]
  [out]
    type = Exodus
    hide = 'solid_temp nek_temp'
  []

  [csv]
    type = CSV
    file_base = 'csv_nek/openmc_nek'
  []
[]

Solid Input Files

The solid heat conduction physics is solved over the solid regions of the unit cell using the MOOSE heat transfer module. The input file for this portion of the physics is the solid_nek.i input. We begin by defining a number of constants and by setting up the mesh for solving heat conduction.

channel_diameter = 0.016                 # diameter of the coolant channels (m)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected
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)
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)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
height = 6.343                           # height of the full core (m)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
triso_pf = 0.15                          # TRISO packing fraction (%)

density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}
unit_cell_mdot = ${fparse mdot / (n_bundles * n_coolant_channels_per_block)}

# 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}

N = 50

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

Because we have a block in the problem that we don't need to define any material properties on, we technically need to turn off a material coverage check, or else we're going to get an error from MOOSE. FEProblem is just the default problem, which we need to list in order to turn off the material coverage check.

[Problem]
  type = FEProblem
  material_coverage_check = false
[]

Next, we define the nonlinear variable that this application will solve for (T), the solid temperature. On the solid blocks, we solve the heat equation, but on the fluid blocks that exist exclusively for transferring data, we add a NullKernel to skip the solve in those regions. On the channel wall, the temperature from NekRS is applied as a Dirichlet condition.

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

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
    block = ${temperature_blocks}
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = ${fuel_blocks}
  []
  [null] # The value of T on the fluid blocks is meaningless
    type = NullKernel
    variable = T
    block = ${density_blocks}
  []
[]

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

[ICs]
  [nek_temp]
    type = FunctionIC
    variable = nek_temp
    function = axial_fluid_temp
  []
[]

Next, we define a number of functions to set the solid material properties and define an initial condition for the wall temperature. The solid material properties are then applied with a HeatConductionMaterial.

[ICs]
  [nek_temp]
    type = FunctionIC
    variable = nek_temp
    function = axial_fluid_temp
  []
[]

[Functions]
  [k_graphite]
    type = ParsedFunction
    expression = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    expression = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    expression = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    symbol_names = 'k_TRISO k_graphite'
    symbol_values = 'k_TRISO k_graphite'
  []
  [axial_fluid_temp]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / ${unit_cell_mdot} / ${fluid_Cp}'
  []
[]

[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 = ${fuel_blocks}
  []
[]

[AuxVariables]
  [nek_temp]
    block = ${density_blocks}
  []
  [nek_bulk_temp]
    block = ${density_blocks}
  []
  [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'
  []
[]

Next, we define a number of auxiliary variables - power will receive the heat source from OpenMC, nek_temp will receive the wall temperature from NekRS, nek_bulk_temp will receive the volumetric fluid temperature from NekRS, and finally flux will be used to compute the wall heat flux to send to NekRS.

[AuxVariables]
  [nek_temp]
    block = ${density_blocks}
  []
  [nek_bulk_temp]
    block = ${density_blocks}
  []
  [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 then add a NekRS sub-application and define the transfers to/from NekRS, including the transfer of fluid temperature from the NekRS volume to the dummy coolant blocks in the MOOSE solid model. Because the NekRS mesh mirror will be a volume mirror (in order to extract volumetric temperatures for the neutronics feedback), a significant cost savings can be obtained by using the "minimal transfer" feature of NekRSProblem (which requires sending a dummy Receiver postprocessor, here named synchronization_to_nek, to indicate when data is to be exchanged). For more information, please consult the documentation for NekRSProblem.

[MultiApps]
  [nek]
    type = TransientMultiApp
    input_files = 'nek.i'
    sub_cycling = true
    execute_on = timestep_end
  []
[]

[Transfers]
  [heat_flux_to_nek]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = flux
    variable = avg_flux
    from_boundaries = 'fluid_solid_interface'
    to_boundaries = '3'
    to_multi_app = nek
  []
  [flux_integral_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = flux_integral
    from_postprocessor = flux_integral
    to_multi_app = nek
  []
  [temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_temp
    from_boundaries = '3'
    to_boundaries = 'fluid_solid_interface'
    from_multi_app = nek
  []
  [bulk_temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_bulk_temp
    from_multi_app = nek
  []
  [synchronization_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = transfer_in
    from_postprocessor = synchronization_to_nek
    to_multi_app = nek
  []
[]

We add several postprocessors to facilitate the data transfers as well as to query the solution.

[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 = ${fuel_blocks}
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [synchronization_to_nek]
    type = Receiver
    default = 1.0
  []
  [power]
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = ${fuel_blocks}
    execute_on = transfer
  []
[]

We add a transient executioner - again, even though the MOOSE heat conduction module in this tutorial solves steady equations, a transient executioner allows us to control the frequency with which MOOSE and NekRS iterate the CHT physics. Here, we use a time step that is 50 times bigger than the NekRS time step.

[Executioner]
  type = Transient
  dt = '${fparse N * nek_dt * t0}'
  nl_abs_tol = 1e-5
  nl_rel_tol = 1e-16
  petsc_options_value = 'hypre boomeramg'
  petsc_options_iname = '-pc_type -pc_hypre_type'
[]

note

The ability to use different data transfer "frequencies" is a key advantage of the "single-stack" MultiApp hierarchy. So far, we have shown that OpenMC will exchange data every 1000 time steps with MOOSE, but MOOSE exchanges data every 50 times steps with NekRS. In other words, for every update of the OpenMC fission distribution, we perform 20 sub-iterations of the CHT physics. Depending on the problem, begin able to iterate the thermal-fluid physics on a finer granularity than the neutronics feedback can be essential to obtaining a stable solution without an inordinately high number of Monte Carlo solves.

Finally, we add a number of LayeredAverage user objects to compute averages of the fuel and graphite temperatures in the axial direction, which are output to CSV using a SpatialUserObjectVectorPostprocessor, in combination with a CSV output.

[UserObjects]
  [average_fuel_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = ${fuel_blocks}
  []
  [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
  hide = 'synchronization_to_nek'

  [csv]
    file_base = 'csv/solid_nek'
    type = CSV
  []
[]

Fluid Input Files

The fluid mass, momentum, and energy transport physics are solved using NekRS. The input file for this portion of the physics is the nek.i input. We begin by defining a number of constants and by setting up the NekRSMesh mesh mirror. Because we are coupling via boundary CHT to MOOSE, we set boundary = '3' so that we will be able to extract the boundary temperature from boundary 3 (the wall). We are also coupling via volumes to OpenMC higher in the MultiApp hierarchy. In order to extract volume representations of the fluid temperature, we also set volume = true.

# copy-pasta from common_input.i

channel_diameter = 0.016                 # diameter of the coolant channels (m)
height = 6.343                           # height of the full core (m)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected

[Mesh]
  type = NekRSMesh
  boundary = '3'
  volume = true
  scaling = ${channel_diameter}
[]

The bulk of the NekRS wrapping occurs in the [Problem] block with NekRSProblem. The NekRS input files are in non-dimensional form, whereas all other coupled applications use dimensional units. The various *_ref and *_0 parameters define the characteristic scales that were used to non-dimensionalize the NekRS input. In order to simplify the input file, we know a priori that OpenMC will not be sending a heat source to NekRS, so we set has_heat_source = false so that we don't need to add a dummy heat source kernel to the ranstube.oudf file. Finally, we indicate that we will be minimizing the data transfers in/out of NekRS unless new data is actually available from the MOOSE heat transfer module with the synchronization_interval = parent_app parameter.

[Problem]
  type = NekRSProblem
  casename = 'ranstube'
  has_heat_source = false
  n_usrwrk_slots = 2

  [Dimensionalize]
    U = '${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}'
    T = ${inlet_T}
    dT = '${fparse power / mdot / fluid_Cp * unit_cell_height / height}'
    L = ${channel_diameter}
    rho = ${fluid_density}
    Cp = ${fluid_Cp}
  []

  synchronization_interval = parent_app
[]

Next, we will allow NekRS to select its own time step using the NekTimeStepper, combined with a transient executioner.

[Executioner]
  type = Transient

  [TimeStepper]
    type = NekTimeStepper
  []
[]

We also add a number of postprocessors to query the Nek solution.

[Postprocessors]
  [inlet_T]
    type = NekSideAverage
    field = temperature
    boundary = '1'
  []
  [outlet_T]
    type = NekSideAverage
    field = temperature
    boundary = '2'
  []
  [max_T]
    type = NekVolumeExtremeValue
    field = temperature
    value_type = max
  []
[]

Finally, we define the output formats and hide the automatically-created flux_integral and transfer_in postprocessors from the screen (console) to have neater output.

[Outputs]
  exodus = true

  [screen]
    type = Console
    hide = 'flux_integral transfer_in'
  []

  [csv]
    file_base = 'csv/nek'
    type = CSV
  []
[]

OpenMC-THM-MOOSE

In this section, we describe the coupling of OpenMC, THM, and MOOSE in the "tree" MultiApp hierarchy shown in Figure 2. For the most part, the OpenMC and MOOSE heat conduction input files only have small differences from those presented earlier for the OpenMC-NekRS-MOOSE coupling. Therefore, for the OpenMC and MOOSE heat conduction cases, we only point out the aspects that differ from the earlier presentation. Differences will exist in several areas:

Due to the use of the "tree" MultiApp hierarchy, different initial conditions are required because the applications execute in a different order
Due to the use of the "tree" MultiApp hierarhcy, the wall heat flux and wall temperature exchanged between MOOSE and THM will need to go "up a level" to their common parent application, while the thermal-fluid application can now directly send fluid temperature to OpenMC
THM requires a slightly different data transfer than NekRS due to its 1-D representation of the flow channels, which requires us to compute an average of the wall heat flux along the wall channel

OpenMC Input Files

The neutronics physics is solved over the entire domain with OpenMC. The OpenMC wrapping used for the OpenMC-THM-MOOSE coupling is described in the openm_thm.i input file. Relative to the OpenMC input files used for the OpenMC-NekRS-MOOSE coupling shown previously, we now need to add variables to hold the wall temperature, thm_temp_wall, that THM passes to the MOOSE heat transfer module, and the wall heat flux, flux, that the MOOSE heat transfer module passes to THM. Both of these are simply receiver variables that are written into by the two sub-applications (MOOSE heat conduction and THM).

[AuxVariables]
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_density]
    family = MONOMIAL
    order = CONSTANT
  []
  [thm_temp_wall]
    family = MONOMIAL
    order = CONSTANT
    block = ${density_blocks}
  []
  [flux]
  []

  [q_wall]
  []
[]

Next, when using the "tree" MultiApp hierarchy, OpenMC runs after the MOOSE heat transfer module, but before THM. Therefore, initial conditions are required for the fluid temperature (for OpenMC) as well as for the fluid wall temperature (which will be sent to the MOOSE heat transfer module on the first time step). Because OpenMC sends its heat source to the MOOSE heat transfer module at the beginning of a time step, we also need to set an initial condition on the heat source.

[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 = ${fuel_blocks}
    value = ${fparse unit_cell_power / (2.0 * pi * compact_diameter * compact_diameter / 4.0 * unit_cell_height)}
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / (${mdot} / ${n_bundles} / ${n_coolant_channels_per_block}) / ${fluid_Cp}'
  []
[]

The [Problem] block is identical to that shown earlier for the OpenMC-NekRS-MOOSE coupling. Next, we define the two sub-applications and the transfers. The transfers of wall heat flux and wall temperature between MOOSE and THM first pass up through OpenMC to one of the sub-applications, giving four additional transfers in the main application's input file than we saw for the "single-stack" hierarchy.

[MultiApps]
  [bison]
    type = TransientMultiApp
    input_files = 'solid_thm.i'
    execute_on = timestep_begin
  []
  [thm]
    type = FullSolveMultiApp
    input_files = 'thm.i'
    execute_on = timestep_end
    max_procs_per_app = 1
    bounding_box_padding = '0.1 0.1 0'
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [heat_flux_to_openmc]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = flux
    variable = flux
    from_multi_app = bison
    from_boundaries = 'fluid_solid_interface'
    to_boundaries = '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
    to_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 = fluid_temp
    to_multi_app = bison
  []

  [q_wall_to_thm]
    type = MultiAppGeneralFieldUserObjectTransfer
    variable = q_wall
    to_multi_app = thm
    source_user_object = q_wall_avg
  []
  [T_wall_from_thm]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = T_wall
    from_multi_app = thm
    variable = thm_temp_wall
  []
  [T_bulk_from_thm]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = T
    from_multi_app = thm
    variable = thm_temp
  []
[]

note

The advantage of the "tree" MultiApp hierarhcy is in the simpler solid input file that we will see in Solid Input Files , which will not need to have a dummy part of the mesh for the fluid region or any NullKernels.

Finally, because we are coupling to THM, we need to average the wall heat flux around the coolant channel into a number of layers to transfer from the 3-D MOOSE heat conduction model to the 1-D fluid flow model. Heat flux is averaged using a LayeredSideAverage.

[UserObjects]
  [q_wall_avg]
    type = LayeredSideAverage
    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}

    direction_min = 0.0
    direction_max = ${unit_cell_height}
  []
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers}
    block = ${fuel_blocks}
  []
[]

All other aspects of the input file are the same as for the OpenMC-NekRS-MOOSE case.

Solid Input Files

The solid heat conduction input file is the solid_thm.i input. This input file is identical to the solid_nek.i input file for the OpenMC-NekRS-MOOSE case except for:

There are no MultiApps, and therefore no transfers, to a thermal-fluid code
There are no initial conditions in the input file, since OpenMC sends the MOOSE heat conduction file all necessary initial conditions for wall temperature and heat source
There is no need to explicitly turn off the material coverage check or add a NullKernel because there is not a dummy fluid block to receive fluid temperature from a sub-application

Because all of these differences are simply omissions, this concludes the discussion of the solid input file. For reference, the full file is below.

# 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}

fuel_blocks = 'compacts compacts_trimmer_tri'

[Mesh]
  [mesh]
   type = FileMeshGenerator
   file = solid_mesh_in.e
  []
  [delete_coolant]
    type = BlockDeletionGenerator
    input = mesh
    block = 'coolant'
  []
[]

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

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = ${fuel_blocks}
  []
[]

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

[Functions]
  [k_graphite]
    type = ParsedFunction
    expression = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    expression = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    expression = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    symbol_names = 'k_TRISO k_graphite'
    symbol_values = 'k_TRISO k_graphite'
  []
[]

[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 = ${fuel_blocks}
  []
[]

[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 = ${fuel_blocks}
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [power]
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = ${fuel_blocks}
    execute_on = 'transfer'
  []
[]

[AuxVariables]
  [fluid_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 = ${fuel_blocks}
  []
  [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
  print_linear_residuals = false

  [csv]
    type = CSV
    file_base = 'csv_thm/solid_out'
  []
[]

Fluid Input Files

The fluid mass, momentum, and energy transport physics are solved using THM. The input file for this portion of the physics is the thm.i input. The THM input file is built using syntax specific to THM - we will only briefly cover the syntax, and instead refer users to the THM documentation 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.

# 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)
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected

num_layers_for_THM = 150

[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 condition, and heat transfer to the pipe wall.

[Components]
  [channel]
    type = FlowChannel1Phase
    position = '0 0 0'
    orientation = '0 0 1'

    A = '${fparse pi * channel_diameter * channel_diameter / 4}'
    D_h = ${channel_diameter}
    length = ${unit_cell_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. For the Dittus-Boelter correlation, we use a corrected version of the closure (with the leading coefficient changed from 0.023 to 0.021) based on the NekRS simulations described in an earlier tutorial. 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 traditioanlly 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
    expression = '0.021 * pow(Re, 0.8) * pow(Pr, 0.4)'
    property_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 transfer 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. We will run THM to convergence based on a tight steady state relative tolerance of $10^{-8}$ .

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

[Executioner]
  type = Transient
  dt = 0.1

  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
  []

  [csv]
    type = CSV
    file_base = 'csv_thm/thm'
  []
[]

Execution and Postprocessing

To run the coupled OpenMC-NekRS-MOOSE calculation, run the following:

mpiexec -np 500 cardinal-opt -i openmc_nek.i

This will run with 500 MPI processes (you may run with other parallel configurations as needed, but you will find that the NekRS simulation requires HPC resources due to the large mesh). To run the coupled OpenMC-THM-MOOSE calculation, run the following:

mpiexec -np 2 cardinal-opt -i openmc_thm.i --n-threads=36

which will run with 2 MPI ranks with 36 threads each (again, these parallel resource choices are only examples). When the two simulations have completed, you will have created a number of different output files. For the OpenMC-NekRS-MOOSE calculation:

openmc_nek_out.e, an Exodus file with the OpenMC solution
openmc_nek_out_bison0.e, an Exodus file with the MOOSE heat conduction solution
openmc_nek_out_bison0_nek0.e, an Exodus file with the NekRS solution
csv_nek/*, CSV files with various postprocessors and userobjects

And for the OpenMC-THM-MOOSE calculation:

openmc_thm_out.e, an Exodus file with the OpenMC solution
openmc_thm_out_bison0.e, an Exodus file with the MOOSE heat conduction solution
openmc_thm_out_thm0.e, an Exodus file with the THM solution
csv_thm/*, CSV files with various postprocessors and userobjects

We now briefly present the coupled physics predictions. Figure 6 shows the fission power predicted by OpenMC with thermal-fluid feedback from NekRS-MOOSE and THM-MOOSE. Also shown is the difference between the two predictions (NekRS case minus the THM case). For the six compacts with 50 axial cell layers, the OpenMC model contains a total of 300 tallies; the maximum tally relative error of 1% does introduce some slight asymmetries among the six compacts, which are easiest to discern in the "difference" image in Figure 6.

Figure 6: Heat source predicted by OpenMC with thermal feedback from either NekRS-MOOSE or THM-MOOSE

The lack of reflectors results in a relatively high leakage percentage of about 18.8%. For the unit cell, recall that the fluid temperature rise is only 82 K. These two effects combine to give a power distribution that is very nearly symmetric in the axial direction - while the lower solid temperatures near the inlet do cause power to shift slightly downwards, the magnitude of the shift is moderated by the effect of pushing the fission source closer to external boundaries, where those neutrons would be more likely to exit the domain.

In the earlier conjugate heat transfer tutorial, we showed that CHT calculations with NekRS-MOOSE and THM-MOOSE agree very well with one another in terms of fluid temperature, solid temperature, and fluid density. Especially when considering that the neutronics feedback due to fluid temperature and density are very small for gas-cooled systems, the excellent agreement in OpenMC's fission distribution shown in Figure 6 is as expected. When cast in terms of a percent difference (as opposed to the absolute difference shown in Figure 6), the NekRS-MOOSE-OpenMC and THM-MOOSE-OpenMC cases agree to within 1%, which is within the range of the uncertainty in the fission tally itself.

Figure 7 shows the solid temperature predicted by the MOOSE heat transfer module, with physics feedback provided by either NekRS-OpenMC or THM-OpenMC. The solid temperature peaks slightly downstream of the maximum OpenMC power due to the combined effects of the power distribution and convective heat transfer at the wall. Figure 7 helps explain the trend in the power difference shown in Figure 6. Near the inlet, NekRS predicts a lower temperature than THM, while near the outlet NekRS predicts a higher temperature than THM. Due to the negative temperature reactivity coefficient of the unit cell, this causes the NekRS-based coupled model to predict a higher power than THM near the inlet, but a lower power than THM near the outlet.

Figure 7: Solid temperature predicted by the MOOSE heat transfer module with physics feedback from either OpenMC-NekRS or OpenMC-THM.

Figure 8 shows the solid temperature predictions on the axial midplane with physics feedback provided by either NekRS-OpenMC or THM-OpenMC. Also shown is the difference between the two (NekRS-based case minus the THM-based case). The THM simulations are unable to capture the radial variation in heat flux along the channel wall, causing THM to underpredict temperatures near the channel wall close to the fuel compacts. However, an additional source of difference is now also present - small radial asymmetries in OpenMC's fission distribution contribute a small asymmetry in solid temperature on the order of 2 K. The tolerance on the fission power uncertainty can simply be increased in order to push down this contribution.

Figure 8: Solid temperature predicted by the MOOSE heat transfer module with physics feedback from either OpenMC-NekRS or OpenMC-THM on the axial mid-plane.

On each axial layer, the OpenMC model receives a total of 14 temperatures (six compacts plus six graphite regions around them, one graphite cell surrounding the coolant channel, and the coolant channel). Figure 9 shows the solid temperature predicted by the NekRS-MOOSE-OpenMC high-resolution simulations along with the actual cell temperature imposed in OpenMC for a half-height slice of the unit cell. OpenMC sets a volume-average temperature for each cell according to the mesh mirror element centroid mappings to OpenMC's cells.

Figure 9: Solid temperature predicted by the MOOSE heat transfer module with physics feedback from OpenMC-NekRS and the solid temperature actually imposed in OpenMC for a half-height portion of the unit cell

Figure 10 shows the fluid temperature predicted by NekRS and THM, with physics feedback provided by MOOSE-OpenMC. NekRS resolves the thermal boundary layer, whereas the THM model uses the Dittus-Boelter correlation to represent the temperature drop from the heated wall to the bulk. Therefore, the fluid temperature shown in Figure 10 for THM is the area-averaged fluid temperature. The wall temperature predicted by NekRS follows a very similar distribution as the heat flux, with a profile that peaks slightly downstream of the maximum power due to the combined effects of convective heat transfer and the fission power distribution.

Figure 10: Fluid temperature predicted for the multiphysics simulations for NekRS-MOOSE-OpenMC and THM-MOOSE-OpenMC.

Figure 11 shows the radially-averaged temperatures for the two thermal-fluid feedback options. The fluid bulk temperature increases along the flow direction, with a faster rate of increase where the power density is highest.

Figure 11: Radially-averaged temperatures for the NekRS-MOOSE-OpenMC and THM-MOOSE-OpenMC simulations.

References

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, J. Fang, H. Yuan, D. Shaver, E. Merzari, P.K. Romano, and R.C. Martineau. Coupled Monte Carlo and Thermal-Fluid Modeling of High Temperature Gas Reactors Using Cardinal. Annals of Nuclear Energy, 177:109310, 2022. doi:10.1016/j.anucene.2022.109310.

@Article{novak2022_cardinal,
    author = "Novak, A.J. and Andrs, D. and Shriwise, P. and Fang, J. and Yuan, H. and Shaver, D. and Merzari, E. and Romano, P.K. and Martineau, R.C.",
    title = "{Coupled {Monte} {Carlo} and Thermal-Fluid Modeling of High Temperature Gas Reactors Using {Cardinal}}",
    journal = "Annals of Nuclear Energy",
    year = "2022",
    volume = "177",
    pages = "109310",
    DOI = "10.1016/j.anucene.2022.109310"
}

(tutorials/gas_compact_multiphysics/unit_cell.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):
    """
    THIS IS ONLY USED FOR SETTING AN INITIAL CONDITION IN OPENMC's XML FILES -
    the coolant temperature will be applied from MOOSE, we just set an initial
    value here in case you want to run these files in standalone mode (i.e. with
    the "openmc" executable).

    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
    distance_from_inlet = specs.unit_cell_height - z
    Qi = (l - l * np.cos(math.pi * distance_from_inlet / l)) / math.pi

    t = t_in + Qi / Q * dT

    return t

def coolant_density(t):
  """
    THIS IS ONLY USED FOR SETTING AN INITIAL CONDITION IN OPENMC's XML FILES -
    the coolant density will be applied from MOOSE, we just set an initial
    value here in case you want to run these files in standalone mode (i.e. with
    the "openmc" executable).

  Computes the helium density (kg/m3) 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
# -------------------------------------------

### RADIANT UNIT CELL SPECS (INFERRED FROM REPORTS) ###

unit_cell_mdot = specs.mdot / (specs.n_bundles * specs.n_coolant_channels_per_block)
unit_cell_power = specs.power / (specs.n_bundles * specs.n_coolant_channels_per_block) * (specs.unit_cell_height / specs.height)

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

# geometry
coolant_channel_diam = specs.channel_diameter * m
reactor_bottom = 0.0
reactor_height = specs.unit_cell_height * m
reactor_top = reactor_bottom + reactor_height

### ARBITRARILY DETERMINED PARAMETERS ###

cell_pitch = specs.fuel_to_coolant_distance * m
fuel_channel_diam = specs.compact_diameter * m

hex_orientation = 'x'

def unit_cell(n_ax_zones, n_inactive, n_active):
    axial_section_height = reactor_height / n_ax_zones

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

    lattice_orientation = 'x'
    cell_edge_length = cell_pitch

    model = openmc.model.Model()

    ### Geometry ###

    # 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)

    # 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, seed=1.0)
    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)

    # create a hexagonal lattice for the coolant and fuel channels
    fuel_univ = openmc.Universe(cells=[openmc.Cell(region=-fuel_cyl, fill=triso_lattice),
                                    openmc.Cell(region=+fuel_cyl, fill=mats.m_graphite_matrix)])

    # extract the coolant cell and set temperatures based on the axial profile
    coolant_cell = openmc.Cell(region=-coolant_cyl, fill=mats.m_coolant)
    # set the coolant temperature on the cell to approximately match the expected
    # temperature profile
    axial_coords = np.linspace(reactor_bottom, reactor_top, n_ax_zones + 1)
    lattice_univs = []

    fuel_ch_cells = []

    i = 0
    for z_min, z_max in zip(axial_coords[0:-1], axial_coords[1:]):
        # create a new coolant universe for each axial zone in the coolant channel;
        # this generates a new material as well (we only need to do this for all
        # cells except the first cell)
        if (i == 0):
          c_cell = coolant_cell
        else:
          c_cell = coolant_cell.clone()

        i += 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)

        c_cell.temperature = t
        c_cell.fill.set_density('kg/m3', coolant_density(t))

        # set the solid cells and their temperatures
        graphite_cell = openmc.Cell(region=+coolant_cyl, fill=mats.m_graphite_matrix)
        fuel_ch_cell = openmc.Cell(region=-fuel_cyl, fill=triso_lattice)
        fuel_ch_matrix_cell = openmc.Cell(region=+fuel_cyl, fill=mats.m_graphite_matrix)

        graphite_cell.temperature = t
        fuel_ch_cell.temperature = t
        fuel_ch_matrix_cell.temperature = t

        fuel_ch_cells.append(fuel_ch_cell)
        fuel_u = openmc.Universe(cells=[fuel_ch_cell, fuel_ch_matrix_cell])
        coolant_u = openmc.Universe(cells=[c_cell, graphite_cell])
        lattice_univs.append([[fuel_u] * 6, [coolant_u]])

    # create a hexagonal lattice used in each axial zone to represent the cell
    hex_lattice = openmc.HexLattice(name="Unit cell lattice")
    hex_lattice.orientation = lattice_orientation
    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

    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

    # hexagonal bounding cell
    hex = openmc.model.HexagonalPrism(cell_edge_length, hex_orientation, boundary_type='periodic')

    hex_cell_vol = 6.0 * (math.sqrt(3) / 4.0) * cell_edge_length**2 * reactor_height

    # create additional axial regions
    axial_planes = [openmc.ZPlane(z0=coord) for coord in axial_coords]
    # axial planes
    min_z = axial_planes[0]
    min_z.boundary_type = 'vacuum'
    max_z = axial_planes[-1]
    max_z.boundary_type = 'vacuum'

    # fill the unit cell with the hex lattice
    hex_cell = openmc.Cell(region=-hex & +min_z & -max_z, fill=hex_lattice)

    model.geometry = openmc.Geometry([hex_cell])

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

    settings.particles = 150000
    settings.inactive = n_inactive
    settings.batches = settings.inactive + n_active

    # the only reason we use 'nearest' here is to be sure we have a robust test for CI;
    # otherwise, 1e-16 differences in temperature (due to numerical roundoff when using
    # different MPI ranks) do change the tracking do to the stochastic interpolation
    settings.temperature['method'] = 'nearest'
    settings.temperature['range'] = (294.0, 1500.0)
    settings.temperature['tolerance'] = 200.0

    hexagon_half_flat = math.sqrt(3.0) / 2.0 * cell_edge_length
    lower_left = (-cell_edge_length, -hexagon_half_flat, reactor_bottom)
    upper_right = (cell_edge_length, hexagon_half_flat, reactor_top)
    source_dist = openmc.stats.Box(lower_left, upper_right)
    source = openmc.IndependentSource(space=source_dist)
    settings.source = source

    model.settings = settings

    m_colors = {}
    m_colors[mats.m_coolant] = 'royalblue'
    m_colors[mats.m_fuel] = 'red'
    m_colors[mats.m_graphite_c_buffer] = 'black'
    m_colors[mats.m_graphite_pyc] = 'orange'
    m_colors[mats.m_sic] = 'yellow'
    m_colors[mats.m_graphite_matrix] = 'silver'

    plot1          = openmc.Plot()
    plot1.filename = 'plot1'
    plot1.width    = (2 * cell_pitch, 4 * axial_section_height)
    plot1.basis    = 'xz'
    plot1.origin   = (0.0, 0.0, reactor_height/2.0)
    plot1.pixels   = (int(800 * 2 * cell_pitch), int(800 * 4 * axial_section_height))
    plot1.color_by = 'cell'

    plot2          = openmc.Plot()
    plot2.filename = 'plot2'
    plot2.width    = (3 * cell_pitch, 3 * cell_pitch)
    plot2.basis    = 'xy'
    plot2.origin   = (0.0, 0.0, axial_section_height / 2.0)
    plot2.pixels   = (int(800 * cell_pitch), int(800 * cell_pitch))
    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('-i', dest='n_inactive', type=int, default=20,
                    help='Number of inactive cycles')
    ap.add_argument('-a', dest='n_active', type=int, default=45,
                    help='Number of active cycles')

    args = ap.parse_args()

    model = unit_cell(args.n_axial, args.n_inactive, args.n_active)
    model.export_to_xml()

if __name__ == "__main__":
    main()

(tutorials/gas_compact_multiphysics/solid_mesh.i)

fuel_to_coolant_distance = 1.88e-2
compact_diameter = 1.27e-2
channel_diameter = 1.6e-2
height = 160e-2
n_layers = 50
ns = 8

[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 = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    ring_block_ids = '2 2'
    ring_block_names = 'compacts compacts'
    background_block_names = 'graphite'
    background_intervals = 4
  []
  [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 = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    ring_block_ids = '101 101'
    ring_block_names = 'coolant coolant'
    background_block_names = 'graphite'
    interface_boundary_id_shift = 100
    background_intervals = 3
  []
  [bundle]
    type = PatternedHexMeshGenerator
    inputs = 'fuel_pin coolant_pin'
    hexagon_size = ${fparse 2.0 * fuel_to_coolant_distance}
    pattern = '0 0;
              0 1 0;
               0 0'
  []
  [trim]
    type = HexagonMeshTrimmer
    input = bundle
    trim_peripheral_region = '1 1 1 1 1 1'
    peripheral_trimming_section_boundary = peripheral_section
  []
  [rotate]
    type = TransformGenerator
    input = trim
    transform = rotate
    vector_value = '30.0 0.0 0.0'
  []
  [extrude]
    type = AdvancedExtruderGenerator
    input = rotate
    heights = ${height}
    num_layers = ${n_layers}
    direction = '0 0 1'
  []
  [fluid_solid_interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = extrude
    primary_block = 'graphite'
    paired_block = 'coolant'
    new_boundary = 'fluid_solid_interface'
  []
[]

# The following content is adding postprocessor(s) to check sideset areas.
# The reactor module is unfortunately quite brittle in its assignment of sideset
# IDs, so we want to be extra sure that any changes to sideset numbering are detected
# in our test suite.
[Problem]
  type = FEProblem
  solve = false
[]

[Postprocessors]
  [area_walls] # should approximate 0.08042477193189872
    type = AreaPostprocessor
    boundary = 'fluid_solid_interface'
  []
  [volume_fuel] # should approximate 0.0004053659832779982
    type = VolumePostprocessor
    block = 'compacts compacts_trimmer_tri'
  []
  [volume_coolant] # should approximate 0.0003216990877275948
    type = VolumePostprocessor
    block = 'coolant'
  []
[]

[Executioner]
  type = Steady
[]

[Outputs]
  csv = true
[]

(tutorials/gas_compact_multiphysics/ranstube.par)

# Model of a turbulent channel flow.
#
# L_ref = 0.016 m          coolant channel diameter
# T_ref = 598 K            coolant inlet temperature
# dT_ref = 82.89 K         coolant nominal temperature rise in the unit cell
# U_ref = 81.089           coolant inlet velocity
# length = 1.6 m           coolant channel flow length
#
# rho_0 = 5.5508 kg/m3     coolant density
# mu_0 = 3.22639e-5 Pa-s   coolant dynamic viscosity
# Cp_0 = 5189 J/kg/K       coolant isobaric heat capacity
# k_0 = 0.2556 W/m/K       coolant thermal conductivity
#
# Re = 223214
# Pr = 0.655
# Pe = 146205

[GENERAL]
  startFrom = converged_cht.fld

  polynomialOrder = 7
  dt = 6e-3
  timeStepper = BDF2

  # end the simulation at specified number of time steps
  stopAt = numSteps
  numSteps = 1000

  # write an output file every specified time steps
  writeControl = steps
  writeInterval = 1000

[PROBLEMTYPE]
  stressFormulation = true

[PRESSURE]
  residualTol = 1e-04

[VELOCITY]
  solver = none
  viscosity = -223214
  density = 1.0
  boundaryTypeMap = v, o, W
  residualTol = 1e-06

[TEMPERATURE]
  rhoCp = 1.0
  conductivity = -146205
  boundaryTypeMap = t, o, f
  residualTol = 1e-06

[SCALAR01] # k
  solver = none
  boundaryTypeMap = t, I, o
  residualTol = 1e-06
  rho = 1.0
  diffusivity = -223214

[SCALAR02] # tau
  solver = none
  boundaryTypeMap = t, I, o
  residualTol = 1e-06
  rho = 1.0
  diffusivity = -223214

(tutorials/gas_compact_multiphysics/ranstube.udf)

#include <math.h>
#include "udf.hpp"

occa::memory o_mu_t_from_restart;

#define length 100.0    // non-dimensional channel length
#define Pr_T 0.91       // turbulent Prandtl number

static occa::kernel scalarScaledAddKernel;

void turbulent_props(nrs_t *nrs, double time, occa::memory o_U, occa::memory o_S,
            occa::memory o_UProp, occa::memory o_SProp)
{
  // fetch the laminar thermal conductivity and specify the desired turbulent Pr
  dfloat k_laminar;
  platform->options.getArgs("SCALAR00 DIFFUSIVITY", k_laminar);

  // o_diff is the turbulent conductivity for all passive scalars, so we grab the first
  // slice, since temperature is the first passive scalar
  occa::memory o_mu = nrs->cds->o_diff + 0 * nrs->cds->fieldOffset[0] * sizeof(dfloat);

  scalarScaledAddKernel(0, nrs->cds->mesh[0]->Nlocal, k_laminar, 1.0 / Pr_T,
    o_mu_t_from_restart /* turbulent viscosity */, o_mu /* laminar viscosity */);
}

void UDF_LoadKernels(occa::properties & kernelInfo)
{
  scalarScaledAddKernel = oudfBuildKernel(kernelInfo, "scalarScaledAdd");
}

void UDF_Setup(nrs_t *nrs)
{
  mesh_t *mesh = nrs->cds->mesh[0];
  int n_gll_points = mesh->Np * mesh->Nelements;

  // allocate space for the k and tau that we will read from the file
  dfloat * mu_t = (dfloat *) calloc(n_gll_points, sizeof(dfloat));

  for (int n = 0; n < n_gll_points; ++n)
  {
    // fetch the restart value of mu_T, which is equal to k * tau, the second and third
    // passive scalars
    mu_t[n] = nrs->cds->S[n + 1 * nrs->cds->fieldOffset[0]] *
              nrs->cds->S[n + 2 * nrs->cds->fieldOffset[0]];
  }

  o_mu_t_from_restart = platform->device.malloc(n_gll_points * sizeof(dfloat), mu_t);

  udf.properties = &turbulent_props;
}

void UDF_ExecuteStep(nrs_t * nrs, double time, int tstep)
{
}

(tutorials/gas_compact_multiphysics/ranstube.oudf)

@kernel void scalarScaledAdd(const dlong offset, const dlong N,
                             const dfloat a,
                             const dfloat b,
                             @restrict const dfloat* X,
                             @restrict dfloat* Y)
{
  for (dlong n = 0; n < N; ++n; @tile(256,@outer,@inner))
    if (n < N)
      Y[n] = a + b * X[n + offset];
}

// Boundary conditions are only needed for temperature, since the solves of pressure,
// velocity, k, and tau are all frozen

void scalarDirichletConditions(bcData *bc)
{
  // inlet temperature
  bc->s = 0.0;
}

void scalarNeumannConditions(bcData *bc)
{
  // wall temperature
  // note: when running with Cardinal, Cardinal will allocate the usrwrk
  // array. If running with NekRS standalone (e.g. nrsmpi), you need to
  // replace the usrwrk with some other value or allocate it youself from
  // the udf and populate it with values.
  bc->flux = bc->usrwrk[bc->idM];
}

(tutorials/gas_compact_multiphysics/openmc_nek.i)

!include common_input.i

# This input file runs coupled OpenMC Monte Carlo transport, MOOSE heat
# conduction, and NekRS fluid flow and heat transfer.
# This input should be run with:
#
# cardinal-opt -i openmc_nek.i

num_layers_for_THM = 150
density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3
N = 1000

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

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

[AuxKernels]
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
  [density]
    type = FluidDensityAux
    variable = density
    p = ${outlet_P}
    T = nek_temp
    fp = helium
    execute_on = 'timestep_begin'
  []
[]

[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
  []
[]

[ICs]
  [temp]
    type = FunctionIC
    variable = nek_temp
    function = temp_ic
  []
  [solid_temp]
    type = FunctionIC
    variable = solid_temp
    function = temp_ic
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / (${mdot} / ${n_bundles} / ${n_coolant_channels_per_block}) / ${fluid_Cp}'
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  power = ${unit_cell_power}
  scaling = 100.0
  density_blocks = ${density_blocks}
  cell_level = 1

  relaxation = robbins_monro

  temperature_variables = 'solid_temp;   nek_temp'
  temperature_blocks = '${temperature_blocks}; ${density_blocks}'

  k_trigger = std_dev
  k_trigger_threshold = 7.5e-4
  batches = 40
  max_batches = 100
  batch_interval = 5

  [Tallies]
    [heat_source]
      type = CellTally
      blocks = ${fuel_blocks}
      name = heat_source

      check_equal_mapped_tally_volumes = true

      trigger = rel_err
      trigger_threshold = 1e-2

      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Executioner]
  type = Transient
  dt = ${fparse N * nek_dt * t0}

  # This is somewhat loose, and you should use a tighter tolerance for production runs;
  # we use 1e-2 to obtain a faster tutorial
  steady_state_detection = true
  check_aux = true
  steady_state_tolerance = 1e-2
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    input_files = 'solid_nek.i'
    execute_on = timestep_end
    sub_cycling = true
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    to_multi_app = bison
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = power
  []
  [temp_from_nek]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = nek_bulk_temp
    from_multi_app = bison
    variable = nek_temp
  []
[]

[Postprocessors]
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_err]
    type = TallyRelativeError
    value_type = max
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KStandardDeviation
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = ${fuel_blocks}
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = ${fuel_blocks}
  []
[]

[UserObjects]
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers_for_plots}
    block = ${fuel_blocks}
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
[]

[Outputs]
  [out]
    type = Exodus
    hide = 'solid_temp nek_temp'
  []

  [csv]
    type = CSV
    file_base = 'csv_nek/openmc_nek'
  []
[]

(tutorials/gas_compact_multiphysics/openmc_nek.i)

!include common_input.i

# This input file runs coupled OpenMC Monte Carlo transport, MOOSE heat
# conduction, and NekRS fluid flow and heat transfer.
# This input should be run with:
#
# cardinal-opt -i openmc_nek.i

num_layers_for_THM = 150
density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3
N = 1000

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

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

[AuxKernels]
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
  [density]
    type = FluidDensityAux
    variable = density
    p = ${outlet_P}
    T = nek_temp
    fp = helium
    execute_on = 'timestep_begin'
  []
[]

[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
  []
[]

[ICs]
  [temp]
    type = FunctionIC
    variable = nek_temp
    function = temp_ic
  []
  [solid_temp]
    type = FunctionIC
    variable = solid_temp
    function = temp_ic
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / (${mdot} / ${n_bundles} / ${n_coolant_channels_per_block}) / ${fluid_Cp}'
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  power = ${unit_cell_power}
  scaling = 100.0
  density_blocks = ${density_blocks}
  cell_level = 1

  relaxation = robbins_monro

  temperature_variables = 'solid_temp;   nek_temp'
  temperature_blocks = '${temperature_blocks}; ${density_blocks}'

  k_trigger = std_dev
  k_trigger_threshold = 7.5e-4
  batches = 40
  max_batches = 100
  batch_interval = 5

  [Tallies]
    [heat_source]
      type = CellTally
      blocks = ${fuel_blocks}
      name = heat_source

      check_equal_mapped_tally_volumes = true

      trigger = rel_err
      trigger_threshold = 1e-2

      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Executioner]
  type = Transient
  dt = ${fparse N * nek_dt * t0}

  # This is somewhat loose, and you should use a tighter tolerance for production runs;
  # we use 1e-2 to obtain a faster tutorial
  steady_state_detection = true
  check_aux = true
  steady_state_tolerance = 1e-2
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    input_files = 'solid_nek.i'
    execute_on = timestep_end
    sub_cycling = true
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    to_multi_app = bison
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = power
  []
  [temp_from_nek]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = nek_bulk_temp
    from_multi_app = bison
    variable = nek_temp
  []
[]

[Postprocessors]
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_err]
    type = TallyRelativeError
    value_type = max
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KStandardDeviation
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = ${fuel_blocks}
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = ${fuel_blocks}
  []
[]

[UserObjects]
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers_for_plots}
    block = ${fuel_blocks}
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
[]

[Outputs]
  [out]
    type = Exodus
    hide = 'solid_temp nek_temp'
  []

  [csv]
    type = CSV
    file_base = 'csv_nek/openmc_nek'
  []
[]

(tutorials/gas_compact_multiphysics/openmc_nek.i)

!include common_input.i

# This input file runs coupled OpenMC Monte Carlo transport, MOOSE heat
# conduction, and NekRS fluid flow and heat transfer.
# This input should be run with:
#
# cardinal-opt -i openmc_nek.i

num_layers_for_THM = 150
density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3
N = 1000

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

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

[AuxKernels]
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
  [density]
    type = FluidDensityAux
    variable = density
    p = ${outlet_P}
    T = nek_temp
    fp = helium
    execute_on = 'timestep_begin'
  []
[]

[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
  []
[]

[ICs]
  [temp]
    type = FunctionIC
    variable = nek_temp
    function = temp_ic
  []
  [solid_temp]
    type = FunctionIC
    variable = solid_temp
    function = temp_ic
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / (${mdot} / ${n_bundles} / ${n_coolant_channels_per_block}) / ${fluid_Cp}'
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  power = ${unit_cell_power}
  scaling = 100.0
  density_blocks = ${density_blocks}
  cell_level = 1

  relaxation = robbins_monro

  temperature_variables = 'solid_temp;   nek_temp'
  temperature_blocks = '${temperature_blocks}; ${density_blocks}'

  k_trigger = std_dev
  k_trigger_threshold = 7.5e-4
  batches = 40
  max_batches = 100
  batch_interval = 5

  [Tallies]
    [heat_source]
      type = CellTally
      blocks = ${fuel_blocks}
      name = heat_source

      check_equal_mapped_tally_volumes = true

      trigger = rel_err
      trigger_threshold = 1e-2

      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Executioner]
  type = Transient
  dt = ${fparse N * nek_dt * t0}

  # This is somewhat loose, and you should use a tighter tolerance for production runs;
  # we use 1e-2 to obtain a faster tutorial
  steady_state_detection = true
  check_aux = true
  steady_state_tolerance = 1e-2
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    input_files = 'solid_nek.i'
    execute_on = timestep_end
    sub_cycling = true
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    to_multi_app = bison
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = power
  []
  [temp_from_nek]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = nek_bulk_temp
    from_multi_app = bison
    variable = nek_temp
  []
[]

[Postprocessors]
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_err]
    type = TallyRelativeError
    value_type = max
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KStandardDeviation
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = ${fuel_blocks}
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = ${fuel_blocks}
  []
[]

[UserObjects]
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers_for_plots}
    block = ${fuel_blocks}
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
[]

[Outputs]
  [out]
    type = Exodus
    hide = 'solid_temp nek_temp'
  []

  [csv]
    type = CSV
    file_base = 'csv_nek/openmc_nek'
  []
[]

(tutorials/gas_compact_multiphysics/openmc_nek.i)

!include common_input.i

# This input file runs coupled OpenMC Monte Carlo transport, MOOSE heat
# conduction, and NekRS fluid flow and heat transfer.
# This input should be run with:
#
# cardinal-opt -i openmc_nek.i

num_layers_for_THM = 150
density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3
N = 1000

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

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

[AuxKernels]
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
  [density]
    type = FluidDensityAux
    variable = density
    p = ${outlet_P}
    T = nek_temp
    fp = helium
    execute_on = 'timestep_begin'
  []
[]

[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
  []
[]

[ICs]
  [temp]
    type = FunctionIC
    variable = nek_temp
    function = temp_ic
  []
  [solid_temp]
    type = FunctionIC
    variable = solid_temp
    function = temp_ic
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / (${mdot} / ${n_bundles} / ${n_coolant_channels_per_block}) / ${fluid_Cp}'
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  power = ${unit_cell_power}
  scaling = 100.0
  density_blocks = ${density_blocks}
  cell_level = 1

  relaxation = robbins_monro

  temperature_variables = 'solid_temp;   nek_temp'
  temperature_blocks = '${temperature_blocks}; ${density_blocks}'

  k_trigger = std_dev
  k_trigger_threshold = 7.5e-4
  batches = 40
  max_batches = 100
  batch_interval = 5

  [Tallies]
    [heat_source]
      type = CellTally
      blocks = ${fuel_blocks}
      name = heat_source

      check_equal_mapped_tally_volumes = true

      trigger = rel_err
      trigger_threshold = 1e-2

      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Executioner]
  type = Transient
  dt = ${fparse N * nek_dt * t0}

  # This is somewhat loose, and you should use a tighter tolerance for production runs;
  # we use 1e-2 to obtain a faster tutorial
  steady_state_detection = true
  check_aux = true
  steady_state_tolerance = 1e-2
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    input_files = 'solid_nek.i'
    execute_on = timestep_end
    sub_cycling = true
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    to_multi_app = bison
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = power
  []
  [temp_from_nek]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = nek_bulk_temp
    from_multi_app = bison
    variable = nek_temp
  []
[]

[Postprocessors]
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_err]
    type = TallyRelativeError
    value_type = max
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KStandardDeviation
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = ${fuel_blocks}
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = ${fuel_blocks}
  []
[]

[UserObjects]
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers_for_plots}
    block = ${fuel_blocks}
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
[]

[Outputs]
  [out]
    type = Exodus
    hide = 'solid_temp nek_temp'
  []

  [csv]
    type = CSV
    file_base = 'csv_nek/openmc_nek'
  []
[]

(tutorials/gas_compact_multiphysics/openmc_nek.i)

!include common_input.i

# This input file runs coupled OpenMC Monte Carlo transport, MOOSE heat
# conduction, and NekRS fluid flow and heat transfer.
# This input should be run with:
#
# cardinal-opt -i openmc_nek.i

num_layers_for_THM = 150
density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3
N = 1000

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

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

[AuxKernels]
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
  [density]
    type = FluidDensityAux
    variable = density
    p = ${outlet_P}
    T = nek_temp
    fp = helium
    execute_on = 'timestep_begin'
  []
[]

[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
  []
[]

[ICs]
  [temp]
    type = FunctionIC
    variable = nek_temp
    function = temp_ic
  []
  [solid_temp]
    type = FunctionIC
    variable = solid_temp
    function = temp_ic
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / (${mdot} / ${n_bundles} / ${n_coolant_channels_per_block}) / ${fluid_Cp}'
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  power = ${unit_cell_power}
  scaling = 100.0
  density_blocks = ${density_blocks}
  cell_level = 1

  relaxation = robbins_monro

  temperature_variables = 'solid_temp;   nek_temp'
  temperature_blocks = '${temperature_blocks}; ${density_blocks}'

  k_trigger = std_dev
  k_trigger_threshold = 7.5e-4
  batches = 40
  max_batches = 100
  batch_interval = 5

  [Tallies]
    [heat_source]
      type = CellTally
      blocks = ${fuel_blocks}
      name = heat_source

      check_equal_mapped_tally_volumes = true

      trigger = rel_err
      trigger_threshold = 1e-2

      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Executioner]
  type = Transient
  dt = ${fparse N * nek_dt * t0}

  # This is somewhat loose, and you should use a tighter tolerance for production runs;
  # we use 1e-2 to obtain a faster tutorial
  steady_state_detection = true
  check_aux = true
  steady_state_tolerance = 1e-2
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    input_files = 'solid_nek.i'
    execute_on = timestep_end
    sub_cycling = true
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    to_multi_app = bison
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = power
  []
  [temp_from_nek]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = nek_bulk_temp
    from_multi_app = bison
    variable = nek_temp
  []
[]

[Postprocessors]
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_err]
    type = TallyRelativeError
    value_type = max
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KStandardDeviation
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = ${fuel_blocks}
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = ${fuel_blocks}
  []
[]

[UserObjects]
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers_for_plots}
    block = ${fuel_blocks}
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
[]

[Outputs]
  [out]
    type = Exodus
    hide = 'solid_temp nek_temp'
  []

  [csv]
    type = CSV
    file_base = 'csv_nek/openmc_nek'
  []
[]

(tutorials/gas_compact_multiphysics/openmc_nek.i)

!include common_input.i

# This input file runs coupled OpenMC Monte Carlo transport, MOOSE heat
# conduction, and NekRS fluid flow and heat transfer.
# This input should be run with:
#
# cardinal-opt -i openmc_nek.i

num_layers_for_THM = 150
density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3
N = 1000

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

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

[AuxKernels]
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
  [density]
    type = FluidDensityAux
    variable = density
    p = ${outlet_P}
    T = nek_temp
    fp = helium
    execute_on = 'timestep_begin'
  []
[]

[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
  []
[]

[ICs]
  [temp]
    type = FunctionIC
    variable = nek_temp
    function = temp_ic
  []
  [solid_temp]
    type = FunctionIC
    variable = solid_temp
    function = temp_ic
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / (${mdot} / ${n_bundles} / ${n_coolant_channels_per_block}) / ${fluid_Cp}'
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  power = ${unit_cell_power}
  scaling = 100.0
  density_blocks = ${density_blocks}
  cell_level = 1

  relaxation = robbins_monro

  temperature_variables = 'solid_temp;   nek_temp'
  temperature_blocks = '${temperature_blocks}; ${density_blocks}'

  k_trigger = std_dev
  k_trigger_threshold = 7.5e-4
  batches = 40
  max_batches = 100
  batch_interval = 5

  [Tallies]
    [heat_source]
      type = CellTally
      blocks = ${fuel_blocks}
      name = heat_source

      check_equal_mapped_tally_volumes = true

      trigger = rel_err
      trigger_threshold = 1e-2

      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Executioner]
  type = Transient
  dt = ${fparse N * nek_dt * t0}

  # This is somewhat loose, and you should use a tighter tolerance for production runs;
  # we use 1e-2 to obtain a faster tutorial
  steady_state_detection = true
  check_aux = true
  steady_state_tolerance = 1e-2
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    input_files = 'solid_nek.i'
    execute_on = timestep_end
    sub_cycling = true
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    to_multi_app = bison
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = power
  []
  [temp_from_nek]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = nek_bulk_temp
    from_multi_app = bison
    variable = nek_temp
  []
[]

[Postprocessors]
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_err]
    type = TallyRelativeError
    value_type = max
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KStandardDeviation
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = ${fuel_blocks}
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = ${fuel_blocks}
  []
[]

[UserObjects]
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers_for_plots}
    block = ${fuel_blocks}
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
[]

[Outputs]
  [out]
    type = Exodus
    hide = 'solid_temp nek_temp'
  []

  [csv]
    type = CSV
    file_base = 'csv_nek/openmc_nek'
  []
[]

(tutorials/gas_compact_multiphysics/openmc_nek.i)

!include common_input.i

# This input file runs coupled OpenMC Monte Carlo transport, MOOSE heat
# conduction, and NekRS fluid flow and heat transfer.
# This input should be run with:
#
# cardinal-opt -i openmc_nek.i

num_layers_for_THM = 150
density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3
N = 1000

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

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

[AuxKernels]
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
  [density]
    type = FluidDensityAux
    variable = density
    p = ${outlet_P}
    T = nek_temp
    fp = helium
    execute_on = 'timestep_begin'
  []
[]

[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
  []
[]

[ICs]
  [temp]
    type = FunctionIC
    variable = nek_temp
    function = temp_ic
  []
  [solid_temp]
    type = FunctionIC
    variable = solid_temp
    function = temp_ic
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / (${mdot} / ${n_bundles} / ${n_coolant_channels_per_block}) / ${fluid_Cp}'
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  power = ${unit_cell_power}
  scaling = 100.0
  density_blocks = ${density_blocks}
  cell_level = 1

  relaxation = robbins_monro

  temperature_variables = 'solid_temp;   nek_temp'
  temperature_blocks = '${temperature_blocks}; ${density_blocks}'

  k_trigger = std_dev
  k_trigger_threshold = 7.5e-4
  batches = 40
  max_batches = 100
  batch_interval = 5

  [Tallies]
    [heat_source]
      type = CellTally
      blocks = ${fuel_blocks}
      name = heat_source

      check_equal_mapped_tally_volumes = true

      trigger = rel_err
      trigger_threshold = 1e-2

      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Executioner]
  type = Transient
  dt = ${fparse N * nek_dt * t0}

  # This is somewhat loose, and you should use a tighter tolerance for production runs;
  # we use 1e-2 to obtain a faster tutorial
  steady_state_detection = true
  check_aux = true
  steady_state_tolerance = 1e-2
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    input_files = 'solid_nek.i'
    execute_on = timestep_end
    sub_cycling = true
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    to_multi_app = bison
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = power
  []
  [temp_from_nek]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = nek_bulk_temp
    from_multi_app = bison
    variable = nek_temp
  []
[]

[Postprocessors]
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_err]
    type = TallyRelativeError
    value_type = max
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KStandardDeviation
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = ${fuel_blocks}
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = ${fuel_blocks}
  []
[]

[UserObjects]
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers_for_plots}
    block = ${fuel_blocks}
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
[]

[Outputs]
  [out]
    type = Exodus
    hide = 'solid_temp nek_temp'
  []

  [csv]
    type = CSV
    file_base = 'csv_nek/openmc_nek'
  []
[]

(tutorials/gas_compact_multiphysics/openmc_nek.i)

!include common_input.i

# This input file runs coupled OpenMC Monte Carlo transport, MOOSE heat
# conduction, and NekRS fluid flow and heat transfer.
# This input should be run with:
#
# cardinal-opt -i openmc_nek.i

num_layers_for_THM = 150
density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3
N = 1000

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

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

[AuxKernels]
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
  [density]
    type = FluidDensityAux
    variable = density
    p = ${outlet_P}
    T = nek_temp
    fp = helium
    execute_on = 'timestep_begin'
  []
[]

[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
  []
[]

[ICs]
  [temp]
    type = FunctionIC
    variable = nek_temp
    function = temp_ic
  []
  [solid_temp]
    type = FunctionIC
    variable = solid_temp
    function = temp_ic
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / (${mdot} / ${n_bundles} / ${n_coolant_channels_per_block}) / ${fluid_Cp}'
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  power = ${unit_cell_power}
  scaling = 100.0
  density_blocks = ${density_blocks}
  cell_level = 1

  relaxation = robbins_monro

  temperature_variables = 'solid_temp;   nek_temp'
  temperature_blocks = '${temperature_blocks}; ${density_blocks}'

  k_trigger = std_dev
  k_trigger_threshold = 7.5e-4
  batches = 40
  max_batches = 100
  batch_interval = 5

  [Tallies]
    [heat_source]
      type = CellTally
      blocks = ${fuel_blocks}
      name = heat_source

      check_equal_mapped_tally_volumes = true

      trigger = rel_err
      trigger_threshold = 1e-2

      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Executioner]
  type = Transient
  dt = ${fparse N * nek_dt * t0}

  # This is somewhat loose, and you should use a tighter tolerance for production runs;
  # we use 1e-2 to obtain a faster tutorial
  steady_state_detection = true
  check_aux = true
  steady_state_tolerance = 1e-2
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    input_files = 'solid_nek.i'
    execute_on = timestep_end
    sub_cycling = true
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    to_multi_app = bison
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = power
  []
  [temp_from_nek]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = nek_bulk_temp
    from_multi_app = bison
    variable = nek_temp
  []
[]

[Postprocessors]
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_err]
    type = TallyRelativeError
    value_type = max
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KStandardDeviation
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = ${fuel_blocks}
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = ${fuel_blocks}
  []
[]

[UserObjects]
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers_for_plots}
    block = ${fuel_blocks}
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
[]

[Outputs]
  [out]
    type = Exodus
    hide = 'solid_temp nek_temp'
  []

  [csv]
    type = CSV
    file_base = 'csv_nek/openmc_nek'
  []
[]

(tutorials/gas_compact_multiphysics/solid_nek.i)

channel_diameter = 0.016                 # diameter of the coolant channels (m)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected
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)
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)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
height = 6.343                           # height of the full core (m)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
triso_pf = 0.15                          # TRISO packing fraction (%)

density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}
unit_cell_mdot = ${fparse mdot / (n_bundles * n_coolant_channels_per_block)}

# 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}

N = 50

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

[Problem]
  type = FEProblem
  material_coverage_check = false
[]

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

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
    block = ${temperature_blocks}
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = ${fuel_blocks}
  []
  [null] # The value of T on the fluid blocks is meaningless
    type = NullKernel
    variable = T
    block = ${density_blocks}
  []
[]

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

[ICs]
  [nek_temp]
    type = FunctionIC
    variable = nek_temp
    function = axial_fluid_temp
  []
[]

[Functions]
  [k_graphite]
    type = ParsedFunction
    expression = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    expression = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    expression = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    symbol_names = 'k_TRISO k_graphite'
    symbol_values = 'k_TRISO k_graphite'
  []
  [axial_fluid_temp]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / ${unit_cell_mdot} / ${fluid_Cp}'
  []
[]

[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 = ${fuel_blocks}
  []
[]

[AuxVariables]
  [nek_temp]
    block = ${density_blocks}
  []
  [nek_bulk_temp]
    block = ${density_blocks}
  []
  [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'
  []
[]

[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 = ${fuel_blocks}
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [synchronization_to_nek]
    type = Receiver
    default = 1.0
  []
  [power]
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = ${fuel_blocks}
    execute_on = transfer
  []
[]

[MultiApps]
  [nek]
    type = TransientMultiApp
    input_files = 'nek.i'
    sub_cycling = true
    execute_on = timestep_end
  []
[]

[Transfers]
  [heat_flux_to_nek]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = flux
    variable = avg_flux
    from_boundaries = 'fluid_solid_interface'
    to_boundaries = '3'
    to_multi_app = nek
  []
  [flux_integral_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = flux_integral
    from_postprocessor = flux_integral
    to_multi_app = nek
  []
  [temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_temp
    from_boundaries = '3'
    to_boundaries = 'fluid_solid_interface'
    from_multi_app = nek
  []
  [bulk_temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_bulk_temp
    from_multi_app = nek
  []
  [synchronization_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = transfer_in
    from_postprocessor = synchronization_to_nek
    to_multi_app = nek
  []
[]

[Executioner]
  type = Transient
  dt = ${fparse N * nek_dt * t0}
  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 = ${fuel_blocks}
  []
  [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
  hide = 'synchronization_to_nek'

  [csv]
    file_base = 'csv/solid_nek'
    type = CSV
  []
[]

(tutorials/gas_compact_multiphysics/solid_nek.i)

channel_diameter = 0.016                 # diameter of the coolant channels (m)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected
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)
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)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
height = 6.343                           # height of the full core (m)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
triso_pf = 0.15                          # TRISO packing fraction (%)

density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}
unit_cell_mdot = ${fparse mdot / (n_bundles * n_coolant_channels_per_block)}

# 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}

N = 50

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

[Problem]
  type = FEProblem
  material_coverage_check = false
[]

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

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
    block = ${temperature_blocks}
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = ${fuel_blocks}
  []
  [null] # The value of T on the fluid blocks is meaningless
    type = NullKernel
    variable = T
    block = ${density_blocks}
  []
[]

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

[ICs]
  [nek_temp]
    type = FunctionIC
    variable = nek_temp
    function = axial_fluid_temp
  []
[]

[Functions]
  [k_graphite]
    type = ParsedFunction
    expression = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    expression = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    expression = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    symbol_names = 'k_TRISO k_graphite'
    symbol_values = 'k_TRISO k_graphite'
  []
  [axial_fluid_temp]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / ${unit_cell_mdot} / ${fluid_Cp}'
  []
[]

[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 = ${fuel_blocks}
  []
[]

[AuxVariables]
  [nek_temp]
    block = ${density_blocks}
  []
  [nek_bulk_temp]
    block = ${density_blocks}
  []
  [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'
  []
[]

[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 = ${fuel_blocks}
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [synchronization_to_nek]
    type = Receiver
    default = 1.0
  []
  [power]
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = ${fuel_blocks}
    execute_on = transfer
  []
[]

[MultiApps]
  [nek]
    type = TransientMultiApp
    input_files = 'nek.i'
    sub_cycling = true
    execute_on = timestep_end
  []
[]

[Transfers]
  [heat_flux_to_nek]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = flux
    variable = avg_flux
    from_boundaries = 'fluid_solid_interface'
    to_boundaries = '3'
    to_multi_app = nek
  []
  [flux_integral_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = flux_integral
    from_postprocessor = flux_integral
    to_multi_app = nek
  []
  [temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_temp
    from_boundaries = '3'
    to_boundaries = 'fluid_solid_interface'
    from_multi_app = nek
  []
  [bulk_temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_bulk_temp
    from_multi_app = nek
  []
  [synchronization_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = transfer_in
    from_postprocessor = synchronization_to_nek
    to_multi_app = nek
  []
[]

[Executioner]
  type = Transient
  dt = ${fparse N * nek_dt * t0}
  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 = ${fuel_blocks}
  []
  [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
  hide = 'synchronization_to_nek'

  [csv]
    file_base = 'csv/solid_nek'
    type = CSV
  []
[]

(tutorials/gas_compact_multiphysics/solid_nek.i)

channel_diameter = 0.016                 # diameter of the coolant channels (m)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected
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)
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)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
height = 6.343                           # height of the full core (m)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
triso_pf = 0.15                          # TRISO packing fraction (%)

density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}
unit_cell_mdot = ${fparse mdot / (n_bundles * n_coolant_channels_per_block)}

# 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}

N = 50

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

[Problem]
  type = FEProblem
  material_coverage_check = false
[]

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

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
    block = ${temperature_blocks}
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = ${fuel_blocks}
  []
  [null] # The value of T on the fluid blocks is meaningless
    type = NullKernel
    variable = T
    block = ${density_blocks}
  []
[]

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

[ICs]
  [nek_temp]
    type = FunctionIC
    variable = nek_temp
    function = axial_fluid_temp
  []
[]

[Functions]
  [k_graphite]
    type = ParsedFunction
    expression = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    expression = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    expression = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    symbol_names = 'k_TRISO k_graphite'
    symbol_values = 'k_TRISO k_graphite'
  []
  [axial_fluid_temp]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / ${unit_cell_mdot} / ${fluid_Cp}'
  []
[]

[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 = ${fuel_blocks}
  []
[]

[AuxVariables]
  [nek_temp]
    block = ${density_blocks}
  []
  [nek_bulk_temp]
    block = ${density_blocks}
  []
  [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'
  []
[]

[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 = ${fuel_blocks}
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [synchronization_to_nek]
    type = Receiver
    default = 1.0
  []
  [power]
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = ${fuel_blocks}
    execute_on = transfer
  []
[]

[MultiApps]
  [nek]
    type = TransientMultiApp
    input_files = 'nek.i'
    sub_cycling = true
    execute_on = timestep_end
  []
[]

[Transfers]
  [heat_flux_to_nek]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = flux
    variable = avg_flux
    from_boundaries = 'fluid_solid_interface'
    to_boundaries = '3'
    to_multi_app = nek
  []
  [flux_integral_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = flux_integral
    from_postprocessor = flux_integral
    to_multi_app = nek
  []
  [temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_temp
    from_boundaries = '3'
    to_boundaries = 'fluid_solid_interface'
    from_multi_app = nek
  []
  [bulk_temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_bulk_temp
    from_multi_app = nek
  []
  [synchronization_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = transfer_in
    from_postprocessor = synchronization_to_nek
    to_multi_app = nek
  []
[]

[Executioner]
  type = Transient
  dt = ${fparse N * nek_dt * t0}
  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 = ${fuel_blocks}
  []
  [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
  hide = 'synchronization_to_nek'

  [csv]
    file_base = 'csv/solid_nek'
    type = CSV
  []
[]

(tutorials/gas_compact_multiphysics/solid_nek.i)

channel_diameter = 0.016                 # diameter of the coolant channels (m)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected
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)
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)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
height = 6.343                           # height of the full core (m)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
triso_pf = 0.15                          # TRISO packing fraction (%)

density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}
unit_cell_mdot = ${fparse mdot / (n_bundles * n_coolant_channels_per_block)}

# 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}

N = 50

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

[Problem]
  type = FEProblem
  material_coverage_check = false
[]

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

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
    block = ${temperature_blocks}
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = ${fuel_blocks}
  []
  [null] # The value of T on the fluid blocks is meaningless
    type = NullKernel
    variable = T
    block = ${density_blocks}
  []
[]

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

[ICs]
  [nek_temp]
    type = FunctionIC
    variable = nek_temp
    function = axial_fluid_temp
  []
[]

[Functions]
  [k_graphite]
    type = ParsedFunction
    expression = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    expression = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    expression = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    symbol_names = 'k_TRISO k_graphite'
    symbol_values = 'k_TRISO k_graphite'
  []
  [axial_fluid_temp]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / ${unit_cell_mdot} / ${fluid_Cp}'
  []
[]

[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 = ${fuel_blocks}
  []
[]

[AuxVariables]
  [nek_temp]
    block = ${density_blocks}
  []
  [nek_bulk_temp]
    block = ${density_blocks}
  []
  [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'
  []
[]

[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 = ${fuel_blocks}
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [synchronization_to_nek]
    type = Receiver
    default = 1.0
  []
  [power]
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = ${fuel_blocks}
    execute_on = transfer
  []
[]

[MultiApps]
  [nek]
    type = TransientMultiApp
    input_files = 'nek.i'
    sub_cycling = true
    execute_on = timestep_end
  []
[]

[Transfers]
  [heat_flux_to_nek]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = flux
    variable = avg_flux
    from_boundaries = 'fluid_solid_interface'
    to_boundaries = '3'
    to_multi_app = nek
  []
  [flux_integral_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = flux_integral
    from_postprocessor = flux_integral
    to_multi_app = nek
  []
  [temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_temp
    from_boundaries = '3'
    to_boundaries = 'fluid_solid_interface'
    from_multi_app = nek
  []
  [bulk_temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_bulk_temp
    from_multi_app = nek
  []
  [synchronization_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = transfer_in
    from_postprocessor = synchronization_to_nek
    to_multi_app = nek
  []
[]

[Executioner]
  type = Transient
  dt = ${fparse N * nek_dt * t0}
  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 = ${fuel_blocks}
  []
  [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
  hide = 'synchronization_to_nek'

  [csv]
    file_base = 'csv/solid_nek'
    type = CSV
  []
[]

(tutorials/gas_compact_multiphysics/solid_nek.i)

channel_diameter = 0.016                 # diameter of the coolant channels (m)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected
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)
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)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
height = 6.343                           # height of the full core (m)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
triso_pf = 0.15                          # TRISO packing fraction (%)

density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}
unit_cell_mdot = ${fparse mdot / (n_bundles * n_coolant_channels_per_block)}

# 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}

N = 50

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

[Problem]
  type = FEProblem
  material_coverage_check = false
[]

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

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
    block = ${temperature_blocks}
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = ${fuel_blocks}
  []
  [null] # The value of T on the fluid blocks is meaningless
    type = NullKernel
    variable = T
    block = ${density_blocks}
  []
[]

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

[ICs]
  [nek_temp]
    type = FunctionIC
    variable = nek_temp
    function = axial_fluid_temp
  []
[]

[Functions]
  [k_graphite]
    type = ParsedFunction
    expression = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    expression = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    expression = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    symbol_names = 'k_TRISO k_graphite'
    symbol_values = 'k_TRISO k_graphite'
  []
  [axial_fluid_temp]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / ${unit_cell_mdot} / ${fluid_Cp}'
  []
[]

[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 = ${fuel_blocks}
  []
[]

[AuxVariables]
  [nek_temp]
    block = ${density_blocks}
  []
  [nek_bulk_temp]
    block = ${density_blocks}
  []
  [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'
  []
[]

[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 = ${fuel_blocks}
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [synchronization_to_nek]
    type = Receiver
    default = 1.0
  []
  [power]
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = ${fuel_blocks}
    execute_on = transfer
  []
[]

[MultiApps]
  [nek]
    type = TransientMultiApp
    input_files = 'nek.i'
    sub_cycling = true
    execute_on = timestep_end
  []
[]

[Transfers]
  [heat_flux_to_nek]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = flux
    variable = avg_flux
    from_boundaries = 'fluid_solid_interface'
    to_boundaries = '3'
    to_multi_app = nek
  []
  [flux_integral_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = flux_integral
    from_postprocessor = flux_integral
    to_multi_app = nek
  []
  [temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_temp
    from_boundaries = '3'
    to_boundaries = 'fluid_solid_interface'
    from_multi_app = nek
  []
  [bulk_temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_bulk_temp
    from_multi_app = nek
  []
  [synchronization_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = transfer_in
    from_postprocessor = synchronization_to_nek
    to_multi_app = nek
  []
[]

[Executioner]
  type = Transient
  dt = ${fparse N * nek_dt * t0}
  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 = ${fuel_blocks}
  []
  [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
  hide = 'synchronization_to_nek'

  [csv]
    file_base = 'csv/solid_nek'
    type = CSV
  []
[]

(tutorials/gas_compact_multiphysics/solid_nek.i)

channel_diameter = 0.016                 # diameter of the coolant channels (m)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected
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)
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)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
height = 6.343                           # height of the full core (m)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
triso_pf = 0.15                          # TRISO packing fraction (%)

density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}
unit_cell_mdot = ${fparse mdot / (n_bundles * n_coolant_channels_per_block)}

# 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}

N = 50

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

[Problem]
  type = FEProblem
  material_coverage_check = false
[]

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

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
    block = ${temperature_blocks}
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = ${fuel_blocks}
  []
  [null] # The value of T on the fluid blocks is meaningless
    type = NullKernel
    variable = T
    block = ${density_blocks}
  []
[]

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

[ICs]
  [nek_temp]
    type = FunctionIC
    variable = nek_temp
    function = axial_fluid_temp
  []
[]

[Functions]
  [k_graphite]
    type = ParsedFunction
    expression = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    expression = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    expression = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    symbol_names = 'k_TRISO k_graphite'
    symbol_values = 'k_TRISO k_graphite'
  []
  [axial_fluid_temp]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / ${unit_cell_mdot} / ${fluid_Cp}'
  []
[]

[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 = ${fuel_blocks}
  []
[]

[AuxVariables]
  [nek_temp]
    block = ${density_blocks}
  []
  [nek_bulk_temp]
    block = ${density_blocks}
  []
  [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'
  []
[]

[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 = ${fuel_blocks}
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [synchronization_to_nek]
    type = Receiver
    default = 1.0
  []
  [power]
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = ${fuel_blocks}
    execute_on = transfer
  []
[]

[MultiApps]
  [nek]
    type = TransientMultiApp
    input_files = 'nek.i'
    sub_cycling = true
    execute_on = timestep_end
  []
[]

[Transfers]
  [heat_flux_to_nek]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = flux
    variable = avg_flux
    from_boundaries = 'fluid_solid_interface'
    to_boundaries = '3'
    to_multi_app = nek
  []
  [flux_integral_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = flux_integral
    from_postprocessor = flux_integral
    to_multi_app = nek
  []
  [temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_temp
    from_boundaries = '3'
    to_boundaries = 'fluid_solid_interface'
    from_multi_app = nek
  []
  [bulk_temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_bulk_temp
    from_multi_app = nek
  []
  [synchronization_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = transfer_in
    from_postprocessor = synchronization_to_nek
    to_multi_app = nek
  []
[]

[Executioner]
  type = Transient
  dt = ${fparse N * nek_dt * t0}
  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 = ${fuel_blocks}
  []
  [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
  hide = 'synchronization_to_nek'

  [csv]
    file_base = 'csv/solid_nek'
    type = CSV
  []
[]

(tutorials/gas_compact_multiphysics/solid_nek.i)

channel_diameter = 0.016                 # diameter of the coolant channels (m)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected
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)
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)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
height = 6.343                           # height of the full core (m)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
triso_pf = 0.15                          # TRISO packing fraction (%)

density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}
unit_cell_mdot = ${fparse mdot / (n_bundles * n_coolant_channels_per_block)}

# 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}

N = 50

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

[Problem]
  type = FEProblem
  material_coverage_check = false
[]

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

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
    block = ${temperature_blocks}
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = ${fuel_blocks}
  []
  [null] # The value of T on the fluid blocks is meaningless
    type = NullKernel
    variable = T
    block = ${density_blocks}
  []
[]

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

[ICs]
  [nek_temp]
    type = FunctionIC
    variable = nek_temp
    function = axial_fluid_temp
  []
[]

[Functions]
  [k_graphite]
    type = ParsedFunction
    expression = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    expression = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    expression = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    symbol_names = 'k_TRISO k_graphite'
    symbol_values = 'k_TRISO k_graphite'
  []
  [axial_fluid_temp]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / ${unit_cell_mdot} / ${fluid_Cp}'
  []
[]

[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 = ${fuel_blocks}
  []
[]

[AuxVariables]
  [nek_temp]
    block = ${density_blocks}
  []
  [nek_bulk_temp]
    block = ${density_blocks}
  []
  [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'
  []
[]

[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 = ${fuel_blocks}
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [synchronization_to_nek]
    type = Receiver
    default = 1.0
  []
  [power]
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = ${fuel_blocks}
    execute_on = transfer
  []
[]

[MultiApps]
  [nek]
    type = TransientMultiApp
    input_files = 'nek.i'
    sub_cycling = true
    execute_on = timestep_end
  []
[]

[Transfers]
  [heat_flux_to_nek]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = flux
    variable = avg_flux
    from_boundaries = 'fluid_solid_interface'
    to_boundaries = '3'
    to_multi_app = nek
  []
  [flux_integral_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = flux_integral
    from_postprocessor = flux_integral
    to_multi_app = nek
  []
  [temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_temp
    from_boundaries = '3'
    to_boundaries = 'fluid_solid_interface'
    from_multi_app = nek
  []
  [bulk_temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_bulk_temp
    from_multi_app = nek
  []
  [synchronization_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = transfer_in
    from_postprocessor = synchronization_to_nek
    to_multi_app = nek
  []
[]

[Executioner]
  type = Transient
  dt = ${fparse N * nek_dt * t0}
  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 = ${fuel_blocks}
  []
  [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
  hide = 'synchronization_to_nek'

  [csv]
    file_base = 'csv/solid_nek'
    type = CSV
  []
[]

(tutorials/gas_compact_multiphysics/solid_nek.i)

channel_diameter = 0.016                 # diameter of the coolant channels (m)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected
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)
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)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
height = 6.343                           # height of the full core (m)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
triso_pf = 0.15                          # TRISO packing fraction (%)

density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}
unit_cell_mdot = ${fparse mdot / (n_bundles * n_coolant_channels_per_block)}

# 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}

N = 50

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

[Problem]
  type = FEProblem
  material_coverage_check = false
[]

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

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
    block = ${temperature_blocks}
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = ${fuel_blocks}
  []
  [null] # The value of T on the fluid blocks is meaningless
    type = NullKernel
    variable = T
    block = ${density_blocks}
  []
[]

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

[ICs]
  [nek_temp]
    type = FunctionIC
    variable = nek_temp
    function = axial_fluid_temp
  []
[]

[Functions]
  [k_graphite]
    type = ParsedFunction
    expression = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    expression = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    expression = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    symbol_names = 'k_TRISO k_graphite'
    symbol_values = 'k_TRISO k_graphite'
  []
  [axial_fluid_temp]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / ${unit_cell_mdot} / ${fluid_Cp}'
  []
[]

[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 = ${fuel_blocks}
  []
[]

[AuxVariables]
  [nek_temp]
    block = ${density_blocks}
  []
  [nek_bulk_temp]
    block = ${density_blocks}
  []
  [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'
  []
[]

[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 = ${fuel_blocks}
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [synchronization_to_nek]
    type = Receiver
    default = 1.0
  []
  [power]
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = ${fuel_blocks}
    execute_on = transfer
  []
[]

[MultiApps]
  [nek]
    type = TransientMultiApp
    input_files = 'nek.i'
    sub_cycling = true
    execute_on = timestep_end
  []
[]

[Transfers]
  [heat_flux_to_nek]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = flux
    variable = avg_flux
    from_boundaries = 'fluid_solid_interface'
    to_boundaries = '3'
    to_multi_app = nek
  []
  [flux_integral_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = flux_integral
    from_postprocessor = flux_integral
    to_multi_app = nek
  []
  [temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_temp
    from_boundaries = '3'
    to_boundaries = 'fluid_solid_interface'
    from_multi_app = nek
  []
  [bulk_temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_bulk_temp
    from_multi_app = nek
  []
  [synchronization_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = transfer_in
    from_postprocessor = synchronization_to_nek
    to_multi_app = nek
  []
[]

[Executioner]
  type = Transient
  dt = ${fparse N * nek_dt * t0}
  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 = ${fuel_blocks}
  []
  [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
  hide = 'synchronization_to_nek'

  [csv]
    file_base = 'csv/solid_nek'
    type = CSV
  []
[]

(tutorials/gas_compact_multiphysics/solid_nek.i)

channel_diameter = 0.016                 # diameter of the coolant channels (m)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected
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)
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)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
height = 6.343                           # height of the full core (m)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
triso_pf = 0.15                          # TRISO packing fraction (%)

density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}
unit_cell_mdot = ${fparse mdot / (n_bundles * n_coolant_channels_per_block)}

# 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}

N = 50

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

[Problem]
  type = FEProblem
  material_coverage_check = false
[]

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

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
    block = ${temperature_blocks}
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = ${fuel_blocks}
  []
  [null] # The value of T on the fluid blocks is meaningless
    type = NullKernel
    variable = T
    block = ${density_blocks}
  []
[]

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

[ICs]
  [nek_temp]
    type = FunctionIC
    variable = nek_temp
    function = axial_fluid_temp
  []
[]

[Functions]
  [k_graphite]
    type = ParsedFunction
    expression = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    expression = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    expression = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    symbol_names = 'k_TRISO k_graphite'
    symbol_values = 'k_TRISO k_graphite'
  []
  [axial_fluid_temp]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / ${unit_cell_mdot} / ${fluid_Cp}'
  []
[]

[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 = ${fuel_blocks}
  []
[]

[AuxVariables]
  [nek_temp]
    block = ${density_blocks}
  []
  [nek_bulk_temp]
    block = ${density_blocks}
  []
  [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'
  []
[]

[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 = ${fuel_blocks}
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [synchronization_to_nek]
    type = Receiver
    default = 1.0
  []
  [power]
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = ${fuel_blocks}
    execute_on = transfer
  []
[]

[MultiApps]
  [nek]
    type = TransientMultiApp
    input_files = 'nek.i'
    sub_cycling = true
    execute_on = timestep_end
  []
[]

[Transfers]
  [heat_flux_to_nek]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = flux
    variable = avg_flux
    from_boundaries = 'fluid_solid_interface'
    to_boundaries = '3'
    to_multi_app = nek
  []
  [flux_integral_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = flux_integral
    from_postprocessor = flux_integral
    to_multi_app = nek
  []
  [temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_temp
    from_boundaries = '3'
    to_boundaries = 'fluid_solid_interface'
    from_multi_app = nek
  []
  [bulk_temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = nek_bulk_temp
    from_multi_app = nek
  []
  [synchronization_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = transfer_in
    from_postprocessor = synchronization_to_nek
    to_multi_app = nek
  []
[]

[Executioner]
  type = Transient
  dt = ${fparse N * nek_dt * t0}
  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 = ${fuel_blocks}
  []
  [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
  hide = 'synchronization_to_nek'

  [csv]
    file_base = 'csv/solid_nek'
    type = CSV
  []
[]

(tutorials/gas_compact_multiphysics/nek.i)

# copy-pasta from common_input.i

channel_diameter = 0.016                 # diameter of the coolant channels (m)
height = 6.343                           # height of the full core (m)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected

[Mesh]
  type = NekRSMesh
  boundary = '3'
  volume = true
  scaling = ${channel_diameter}
[]

[Problem]
  type = NekRSProblem
  casename = 'ranstube'
  has_heat_source = false
  n_usrwrk_slots = 2

  [Dimensionalize]
    U = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
    T = ${inlet_T}
    dT = ${fparse power / mdot / fluid_Cp * unit_cell_height / height}
    L = ${channel_diameter}
    rho = ${fluid_density}
    Cp = ${fluid_Cp}
  []

  synchronization_interval = parent_app
[]

[Executioner]
  type = Transient

  [TimeStepper]
    type = NekTimeStepper
  []
[]

[Outputs]
  exodus = true

  [screen]
    type = Console
    hide = 'flux_integral transfer_in'
  []

  [csv]
    file_base = 'csv/nek'
    type = CSV
  []
[]

[Postprocessors]
  [inlet_T]
    type = NekSideAverage
    field = temperature
    boundary = '1'
  []
  [outlet_T]
    type = NekSideAverage
    field = temperature
    boundary = '2'
  []
  [max_T]
    type = NekVolumeExtremeValue
    field = temperature
    value_type = max
  []
[]

(tutorials/gas_compact_multiphysics/nek.i)

# copy-pasta from common_input.i

channel_diameter = 0.016                 # diameter of the coolant channels (m)
height = 6.343                           # height of the full core (m)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected

[Mesh]
  type = NekRSMesh
  boundary = '3'
  volume = true
  scaling = ${channel_diameter}
[]

[Problem]
  type = NekRSProblem
  casename = 'ranstube'
  has_heat_source = false
  n_usrwrk_slots = 2

  [Dimensionalize]
    U = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
    T = ${inlet_T}
    dT = ${fparse power / mdot / fluid_Cp * unit_cell_height / height}
    L = ${channel_diameter}
    rho = ${fluid_density}
    Cp = ${fluid_Cp}
  []

  synchronization_interval = parent_app
[]

[Executioner]
  type = Transient

  [TimeStepper]
    type = NekTimeStepper
  []
[]

[Outputs]
  exodus = true

  [screen]
    type = Console
    hide = 'flux_integral transfer_in'
  []

  [csv]
    file_base = 'csv/nek'
    type = CSV
  []
[]

[Postprocessors]
  [inlet_T]
    type = NekSideAverage
    field = temperature
    boundary = '1'
  []
  [outlet_T]
    type = NekSideAverage
    field = temperature
    boundary = '2'
  []
  [max_T]
    type = NekVolumeExtremeValue
    field = temperature
    value_type = max
  []
[]

(tutorials/gas_compact_multiphysics/nek.i)

# copy-pasta from common_input.i

channel_diameter = 0.016                 # diameter of the coolant channels (m)
height = 6.343                           # height of the full core (m)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected

[Mesh]
  type = NekRSMesh
  boundary = '3'
  volume = true
  scaling = ${channel_diameter}
[]

[Problem]
  type = NekRSProblem
  casename = 'ranstube'
  has_heat_source = false
  n_usrwrk_slots = 2

  [Dimensionalize]
    U = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
    T = ${inlet_T}
    dT = ${fparse power / mdot / fluid_Cp * unit_cell_height / height}
    L = ${channel_diameter}
    rho = ${fluid_density}
    Cp = ${fluid_Cp}
  []

  synchronization_interval = parent_app
[]

[Executioner]
  type = Transient

  [TimeStepper]
    type = NekTimeStepper
  []
[]

[Outputs]
  exodus = true

  [screen]
    type = Console
    hide = 'flux_integral transfer_in'
  []

  [csv]
    file_base = 'csv/nek'
    type = CSV
  []
[]

[Postprocessors]
  [inlet_T]
    type = NekSideAverage
    field = temperature
    boundary = '1'
  []
  [outlet_T]
    type = NekSideAverage
    field = temperature
    boundary = '2'
  []
  [max_T]
    type = NekVolumeExtremeValue
    field = temperature
    value_type = max
  []
[]

(tutorials/gas_compact_multiphysics/nek.i)

# copy-pasta from common_input.i

channel_diameter = 0.016                 # diameter of the coolant channels (m)
height = 6.343                           # height of the full core (m)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected

[Mesh]
  type = NekRSMesh
  boundary = '3'
  volume = true
  scaling = ${channel_diameter}
[]

[Problem]
  type = NekRSProblem
  casename = 'ranstube'
  has_heat_source = false
  n_usrwrk_slots = 2

  [Dimensionalize]
    U = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
    T = ${inlet_T}
    dT = ${fparse power / mdot / fluid_Cp * unit_cell_height / height}
    L = ${channel_diameter}
    rho = ${fluid_density}
    Cp = ${fluid_Cp}
  []

  synchronization_interval = parent_app
[]

[Executioner]
  type = Transient

  [TimeStepper]
    type = NekTimeStepper
  []
[]

[Outputs]
  exodus = true

  [screen]
    type = Console
    hide = 'flux_integral transfer_in'
  []

  [csv]
    file_base = 'csv/nek'
    type = CSV
  []
[]

[Postprocessors]
  [inlet_T]
    type = NekSideAverage
    field = temperature
    boundary = '1'
  []
  [outlet_T]
    type = NekSideAverage
    field = temperature
    boundary = '2'
  []
  [max_T]
    type = NekVolumeExtremeValue
    field = temperature
    value_type = max
  []
[]

(tutorials/gas_compact_multiphysics/nek.i)

# copy-pasta from common_input.i

channel_diameter = 0.016                 # diameter of the coolant channels (m)
height = 6.343                           # height of the full core (m)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected

[Mesh]
  type = NekRSMesh
  boundary = '3'
  volume = true
  scaling = ${channel_diameter}
[]

[Problem]
  type = NekRSProblem
  casename = 'ranstube'
  has_heat_source = false
  n_usrwrk_slots = 2

  [Dimensionalize]
    U = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
    T = ${inlet_T}
    dT = ${fparse power / mdot / fluid_Cp * unit_cell_height / height}
    L = ${channel_diameter}
    rho = ${fluid_density}
    Cp = ${fluid_Cp}
  []

  synchronization_interval = parent_app
[]

[Executioner]
  type = Transient

  [TimeStepper]
    type = NekTimeStepper
  []
[]

[Outputs]
  exodus = true

  [screen]
    type = Console
    hide = 'flux_integral transfer_in'
  []

  [csv]
    file_base = 'csv/nek'
    type = CSV
  []
[]

[Postprocessors]
  [inlet_T]
    type = NekSideAverage
    field = temperature
    boundary = '1'
  []
  [outlet_T]
    type = NekSideAverage
    field = temperature
    boundary = '2'
  []
  [max_T]
    type = NekVolumeExtremeValue
    field = temperature
    value_type = max
  []
[]

(tutorials/gas_compact_multiphysics/openmc_thm.i)

!include common_input.i

# This input file runs coupled OpenMC Monte Carlo transport, MOOSE heat
# conduction, and THM fluid flow and heat transfer.
# This input should be run with:
#
# cardinal-opt -i openmc_thm.i

num_layers_for_THM = 150
num_layers = 50
density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

[AuxVariables]
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_density]
    family = MONOMIAL
    order = CONSTANT
  []
  [thm_temp_wall]
    family = MONOMIAL
    order = CONSTANT
    block = ${density_blocks}
  []
  [flux]
  []
[]

[AuxKernels]
  [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'
  []
[]

[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
  []
[]

[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 = ${fuel_blocks}
    value = ${fparse unit_cell_power / (2.0 * pi * compact_diameter * compact_diameter / 4.0 * unit_cell_height)}
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / (${mdot} / ${n_bundles} / ${n_coolant_channels_per_block}) / ${fluid_Cp}'
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  power = ${unit_cell_power}
  scaling = 100.0
  density_blocks = ${density_blocks}
  cell_level = 1

  relaxation = robbins_monro

  temperature_variables = 'solid_temp; thm_temp'
  temperature_blocks = '${temperature_blocks}; ${density_blocks}'

  k_trigger = std_dev
  k_trigger_threshold = 7.5e-4
  batches = 40
  max_batches = 100
  batch_interval = 5

  [Tallies]
    [heat_source]
      type = CellTally
      blocks = ${fuel_blocks}
      name = heat_source

      check_equal_mapped_tally_volumes = true

      trigger = rel_err
      trigger_threshold = 1e-2

      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Executioner]
  type = Transient

  # We use a fairly loose tolerance here for a tutorial; you may consider increasing this
  # for production runs
  steady_state_detection = true
  check_aux = true
  steady_state_tolerance = 5e-3
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    input_files = 'solid_thm.i'
    execute_on = timestep_begin
  []
  [thm]
    type = FullSolveMultiApp
    input_files = 'thm.i'
    execute_on = timestep_end
    max_procs_per_app = 1
    bounding_box_padding = '0.1 0.1 0'
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [heat_flux_to_openmc]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = flux
    variable = flux
    from_multi_app = bison
    from_boundaries = 'fluid_solid_interface'
    to_boundaries = '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
    to_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 = fluid_temp
    to_multi_app = bison
  []

  [q_wall_to_thm]
    type = MultiAppGeneralFieldUserObjectTransfer
    variable = q_wall
    to_multi_app = thm
    source_user_object = q_wall_avg
  []
  [T_wall_from_thm]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = T_wall
    from_multi_app = thm
    variable = thm_temp_wall
  []
  [T_bulk_from_thm]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = T
    from_multi_app = thm
    variable = thm_temp
  []
[]

[Postprocessors]
  [flux_integral]
    type = SideIntegralVariablePostprocessor
    variable = flux
    boundary = 'fluid_solid_interface'
    execute_on = 'transfer'
  []
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_err]
    type = TallyRelativeError
    value_type = max
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KStandardDeviation
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = ${fuel_blocks}
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = ${fuel_blocks}
  []
[]


[AuxVariables]
  [q_wall]
  []
[]

[AuxKernels]
  [q_wall]
    type = SpatialUserObjectAux
    variable = q_wall
    user_object = q_wall_avg
  []
[]

[UserObjects]
  [q_wall_avg]
    type = LayeredSideAverage
    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}

    direction_min = 0.0
    direction_max = ${unit_cell_height}
  []
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers}
    block = ${fuel_blocks}
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
[]

[Outputs]
  exodus = true

  [csv]
    type = CSV
    file_base = 'csv_thm/openmc_thm'
  []
[]

(tutorials/gas_compact_multiphysics/openmc_thm.i)

!include common_input.i

# This input file runs coupled OpenMC Monte Carlo transport, MOOSE heat
# conduction, and THM fluid flow and heat transfer.
# This input should be run with:
#
# cardinal-opt -i openmc_thm.i

num_layers_for_THM = 150
num_layers = 50
density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

[AuxVariables]
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_density]
    family = MONOMIAL
    order = CONSTANT
  []
  [thm_temp_wall]
    family = MONOMIAL
    order = CONSTANT
    block = ${density_blocks}
  []
  [flux]
  []
[]

[AuxKernels]
  [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'
  []
[]

[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
  []
[]

[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 = ${fuel_blocks}
    value = ${fparse unit_cell_power / (2.0 * pi * compact_diameter * compact_diameter / 4.0 * unit_cell_height)}
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / (${mdot} / ${n_bundles} / ${n_coolant_channels_per_block}) / ${fluid_Cp}'
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  power = ${unit_cell_power}
  scaling = 100.0
  density_blocks = ${density_blocks}
  cell_level = 1

  relaxation = robbins_monro

  temperature_variables = 'solid_temp; thm_temp'
  temperature_blocks = '${temperature_blocks}; ${density_blocks}'

  k_trigger = std_dev
  k_trigger_threshold = 7.5e-4
  batches = 40
  max_batches = 100
  batch_interval = 5

  [Tallies]
    [heat_source]
      type = CellTally
      blocks = ${fuel_blocks}
      name = heat_source

      check_equal_mapped_tally_volumes = true

      trigger = rel_err
      trigger_threshold = 1e-2

      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Executioner]
  type = Transient

  # We use a fairly loose tolerance here for a tutorial; you may consider increasing this
  # for production runs
  steady_state_detection = true
  check_aux = true
  steady_state_tolerance = 5e-3
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    input_files = 'solid_thm.i'
    execute_on = timestep_begin
  []
  [thm]
    type = FullSolveMultiApp
    input_files = 'thm.i'
    execute_on = timestep_end
    max_procs_per_app = 1
    bounding_box_padding = '0.1 0.1 0'
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [heat_flux_to_openmc]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = flux
    variable = flux
    from_multi_app = bison
    from_boundaries = 'fluid_solid_interface'
    to_boundaries = '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
    to_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 = fluid_temp
    to_multi_app = bison
  []

  [q_wall_to_thm]
    type = MultiAppGeneralFieldUserObjectTransfer
    variable = q_wall
    to_multi_app = thm
    source_user_object = q_wall_avg
  []
  [T_wall_from_thm]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = T_wall
    from_multi_app = thm
    variable = thm_temp_wall
  []
  [T_bulk_from_thm]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = T
    from_multi_app = thm
    variable = thm_temp
  []
[]

[Postprocessors]
  [flux_integral]
    type = SideIntegralVariablePostprocessor
    variable = flux
    boundary = 'fluid_solid_interface'
    execute_on = 'transfer'
  []
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_err]
    type = TallyRelativeError
    value_type = max
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KStandardDeviation
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = ${fuel_blocks}
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = ${fuel_blocks}
  []
[]


[AuxVariables]
  [q_wall]
  []
[]

[AuxKernels]
  [q_wall]
    type = SpatialUserObjectAux
    variable = q_wall
    user_object = q_wall_avg
  []
[]

[UserObjects]
  [q_wall_avg]
    type = LayeredSideAverage
    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}

    direction_min = 0.0
    direction_max = ${unit_cell_height}
  []
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers}
    block = ${fuel_blocks}
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
[]

[Outputs]
  exodus = true

  [csv]
    type = CSV
    file_base = 'csv_thm/openmc_thm'
  []
[]

(tutorials/gas_compact_multiphysics/openmc_thm.i)

!include common_input.i

# This input file runs coupled OpenMC Monte Carlo transport, MOOSE heat
# conduction, and THM fluid flow and heat transfer.
# This input should be run with:
#
# cardinal-opt -i openmc_thm.i

num_layers_for_THM = 150
num_layers = 50
density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

[AuxVariables]
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_density]
    family = MONOMIAL
    order = CONSTANT
  []
  [thm_temp_wall]
    family = MONOMIAL
    order = CONSTANT
    block = ${density_blocks}
  []
  [flux]
  []
[]

[AuxKernels]
  [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'
  []
[]

[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
  []
[]

[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 = ${fuel_blocks}
    value = ${fparse unit_cell_power / (2.0 * pi * compact_diameter * compact_diameter / 4.0 * unit_cell_height)}
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / (${mdot} / ${n_bundles} / ${n_coolant_channels_per_block}) / ${fluid_Cp}'
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  power = ${unit_cell_power}
  scaling = 100.0
  density_blocks = ${density_blocks}
  cell_level = 1

  relaxation = robbins_monro

  temperature_variables = 'solid_temp; thm_temp'
  temperature_blocks = '${temperature_blocks}; ${density_blocks}'

  k_trigger = std_dev
  k_trigger_threshold = 7.5e-4
  batches = 40
  max_batches = 100
  batch_interval = 5

  [Tallies]
    [heat_source]
      type = CellTally
      blocks = ${fuel_blocks}
      name = heat_source

      check_equal_mapped_tally_volumes = true

      trigger = rel_err
      trigger_threshold = 1e-2

      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Executioner]
  type = Transient

  # We use a fairly loose tolerance here for a tutorial; you may consider increasing this
  # for production runs
  steady_state_detection = true
  check_aux = true
  steady_state_tolerance = 5e-3
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    input_files = 'solid_thm.i'
    execute_on = timestep_begin
  []
  [thm]
    type = FullSolveMultiApp
    input_files = 'thm.i'
    execute_on = timestep_end
    max_procs_per_app = 1
    bounding_box_padding = '0.1 0.1 0'
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [heat_flux_to_openmc]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = flux
    variable = flux
    from_multi_app = bison
    from_boundaries = 'fluid_solid_interface'
    to_boundaries = '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
    to_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 = fluid_temp
    to_multi_app = bison
  []

  [q_wall_to_thm]
    type = MultiAppGeneralFieldUserObjectTransfer
    variable = q_wall
    to_multi_app = thm
    source_user_object = q_wall_avg
  []
  [T_wall_from_thm]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = T_wall
    from_multi_app = thm
    variable = thm_temp_wall
  []
  [T_bulk_from_thm]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = T
    from_multi_app = thm
    variable = thm_temp
  []
[]

[Postprocessors]
  [flux_integral]
    type = SideIntegralVariablePostprocessor
    variable = flux
    boundary = 'fluid_solid_interface'
    execute_on = 'transfer'
  []
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_err]
    type = TallyRelativeError
    value_type = max
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KStandardDeviation
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = ${fuel_blocks}
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = ${fuel_blocks}
  []
[]


[AuxVariables]
  [q_wall]
  []
[]

[AuxKernels]
  [q_wall]
    type = SpatialUserObjectAux
    variable = q_wall
    user_object = q_wall_avg
  []
[]

[UserObjects]
  [q_wall_avg]
    type = LayeredSideAverage
    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}

    direction_min = 0.0
    direction_max = ${unit_cell_height}
  []
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers}
    block = ${fuel_blocks}
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
[]

[Outputs]
  exodus = true

  [csv]
    type = CSV
    file_base = 'csv_thm/openmc_thm'
  []
[]

(tutorials/gas_compact_multiphysics/openmc_thm.i)

!include common_input.i

# This input file runs coupled OpenMC Monte Carlo transport, MOOSE heat
# conduction, and THM fluid flow and heat transfer.
# This input should be run with:
#
# cardinal-opt -i openmc_thm.i

num_layers_for_THM = 150
num_layers = 50
density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

[AuxVariables]
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_density]
    family = MONOMIAL
    order = CONSTANT
  []
  [thm_temp_wall]
    family = MONOMIAL
    order = CONSTANT
    block = ${density_blocks}
  []
  [flux]
  []
[]

[AuxKernels]
  [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'
  []
[]

[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
  []
[]

[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 = ${fuel_blocks}
    value = ${fparse unit_cell_power / (2.0 * pi * compact_diameter * compact_diameter / 4.0 * unit_cell_height)}
  []
[]

[Functions]
  [temp_ic]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / (${mdot} / ${n_bundles} / ${n_coolant_channels_per_block}) / ${fluid_Cp}'
  []
[]

[Problem]
  type = OpenMCCellAverageProblem

  power = ${unit_cell_power}
  scaling = 100.0
  density_blocks = ${density_blocks}
  cell_level = 1

  relaxation = robbins_monro

  temperature_variables = 'solid_temp; thm_temp'
  temperature_blocks = '${temperature_blocks}; ${density_blocks}'

  k_trigger = std_dev
  k_trigger_threshold = 7.5e-4
  batches = 40
  max_batches = 100
  batch_interval = 5

  [Tallies]
    [heat_source]
      type = CellTally
      blocks = ${fuel_blocks}
      name = heat_source

      check_equal_mapped_tally_volumes = true

      trigger = rel_err
      trigger_threshold = 1e-2

      output = 'unrelaxed_tally_std_dev'
    []
  []
[]

[Executioner]
  type = Transient

  # We use a fairly loose tolerance here for a tutorial; you may consider increasing this
  # for production runs
  steady_state_detection = true
  check_aux = true
  steady_state_tolerance = 5e-3
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    input_files = 'solid_thm.i'
    execute_on = timestep_begin
  []
  [thm]
    type = FullSolveMultiApp
    input_files = 'thm.i'
    execute_on = timestep_end
    max_procs_per_app = 1
    bounding_box_padding = '0.1 0.1 0'
  []
[]

[Transfers]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer
    source_variable = T
    variable = solid_temp
    from_multi_app = bison
  []
  [heat_flux_to_openmc]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = flux
    variable = flux
    from_multi_app = bison
    from_boundaries = 'fluid_solid_interface'
    to_boundaries = '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
    to_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 = fluid_temp
    to_multi_app = bison
  []

  [q_wall_to_thm]
    type = MultiAppGeneralFieldUserObjectTransfer
    variable = q_wall
    to_multi_app = thm
    source_user_object = q_wall_avg
  []
  [T_wall_from_thm]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = T_wall
    from_multi_app = thm
    variable = thm_temp_wall
  []
  [T_bulk_from_thm]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = T
    from_multi_app = thm
    variable = thm_temp
  []
[]

[Postprocessors]
  [flux_integral]
    type = SideIntegralVariablePostprocessor
    variable = flux
    boundary = 'fluid_solid_interface'
    execute_on = 'transfer'
  []
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = 'transfer initial timestep_end'
  []
  [max_tally_err]
    type = TallyRelativeError
    value_type = max
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KStandardDeviation
  []
  [min_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = min
    block = ${fuel_blocks}
  []
  [max_power]
    type = ElementExtremeValue
    variable = heat_source
    value_type = max
    block = ${fuel_blocks}
  []
[]


[AuxVariables]
  [q_wall]
  []
[]

[AuxKernels]
  [q_wall]
    type = SpatialUserObjectAux
    variable = q_wall
    user_object = q_wall_avg
  []
[]

[UserObjects]
  [q_wall_avg]
    type = LayeredSideAverage
    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}

    direction_min = 0.0
    direction_max = ${unit_cell_height}
  []
  [average_power_axial]
    type = LayeredAverage
    variable = heat_source
    direction = z
    num_layers = ${num_layers}
    block = ${fuel_blocks}
  []
[]

[VectorPostprocessors]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_power_axial
  []
[]

[Outputs]
  exodus = true

  [csv]
    type = CSV
    file_base = 'csv_thm/openmc_thm'
  []
[]

(tutorials/gas_compact_multiphysics/solid_thm.i)

# 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}

fuel_blocks = 'compacts compacts_trimmer_tri'

[Mesh]
  [mesh]
   type = FileMeshGenerator
   file = solid_mesh_in.e
  []
  [delete_coolant]
    type = BlockDeletionGenerator
    input = mesh
    block = 'coolant'
  []
[]

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

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = ${fuel_blocks}
  []
[]

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

[Functions]
  [k_graphite]
    type = ParsedFunction
    expression = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    expression = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    expression = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    symbol_names = 'k_TRISO k_graphite'
    symbol_values = 'k_TRISO k_graphite'
  []
[]

[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 = ${fuel_blocks}
  []
[]

[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 = ${fuel_blocks}
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [power]
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = ${fuel_blocks}
    execute_on = 'transfer'
  []
[]

[AuxVariables]
  [fluid_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 = ${fuel_blocks}
  []
  [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
  print_linear_residuals = false

  [csv]
    type = CSV
    file_base = 'csv_thm/solid_out'
  []
[]

(tutorials/gas_compact_multiphysics/thm.i)

# 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)
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected

num_layers_for_THM = 150

[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
    expression = '0.021 * pow(Re, 0.8) * pow(Pr, 0.4)'
    property_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 0'
    orientation = '0 0 1'

    A = ${fparse pi * channel_diameter * channel_diameter / 4}
    D_h = ${channel_diameter}
    length = ${unit_cell_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

  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
  []

  [csv]
    type = CSV
    file_base = 'csv_thm/thm'
  []
[]

(tutorials/gas_compact_multiphysics/thm.i)

# 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)
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected

num_layers_for_THM = 150

[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
    expression = '0.021 * pow(Re, 0.8) * pow(Pr, 0.4)'
    property_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 0'
    orientation = '0 0 1'

    A = ${fparse pi * channel_diameter * channel_diameter / 4}
    D_h = ${channel_diameter}
    length = ${unit_cell_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

  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
  []

  [csv]
    type = CSV
    file_base = 'csv_thm/thm'
  []
[]

(tutorials/gas_compact_multiphysics/thm.i)

# 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)
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected

num_layers_for_THM = 150

[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
    expression = '0.021 * pow(Re, 0.8) * pow(Pr, 0.4)'
    property_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 0'
    orientation = '0 0 1'

    A = ${fparse pi * channel_diameter * channel_diameter / 4}
    D_h = ${channel_diameter}
    length = ${unit_cell_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

  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
  []

  [csv]
    type = CSV
    file_base = 'csv_thm/thm'
  []
[]

(tutorials/gas_compact_multiphysics/thm.i)

# 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)
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected

num_layers_for_THM = 150

[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
    expression = '0.021 * pow(Re, 0.8) * pow(Pr, 0.4)'
    property_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 0'
    orientation = '0 0 1'

    A = ${fparse pi * channel_diameter * channel_diameter / 4}
    D_h = ${channel_diameter}
    length = ${unit_cell_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

  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
  []

  [csv]
    type = CSV
    file_base = 'csv_thm/thm'
  []
[]

(tutorials/gas_compact_multiphysics/thm.i)

# 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)
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected

num_layers_for_THM = 150

[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
    expression = '0.021 * pow(Re, 0.8) * pow(Pr, 0.4)'
    property_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 0'
    orientation = '0 0 1'

    A = ${fparse pi * channel_diameter * channel_diameter / 4}
    D_h = ${channel_diameter}
    length = ${unit_cell_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

  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
  []

  [csv]
    type = CSV
    file_base = 'csv_thm/thm'
  []
[]

Geometry and Computational Model
Multiphysics Coupling
Execution and Postprocessing
References