Thermal coupling of multiple SCM assemblies

Contact: Vasileios Kyriakopoulos, vasileios.kyriakopoulos.at.inl.gov

Model link: Multiple SCM assemblies in a wrapper/inter-wrapper configuration

We include two example cases of coupling multiple subchannel models (SCM) with a 3D heat conduction model (MOOSE) of the wrapper and inter-wrapper. The wrapper is the hexagonal duct that contains the subchannel assemblies and the inter-wrapper is the space in-between the ducts. A picture of the wrapper and inter-wrapper for the case of 19 assemblies is shown in Figure 1

Figure 1: Wrapper and inter-wrapper of 19 hexagonal assemblies

The index of the fuel-pins/subchannels and their relative location for the SCM model is shown in Figure 2

Figure 2: Index and location of fuel pin centers and subchannel centroids

Example Description

One example models 7 subchannel assemblies with a heat conduction model of the wrapper/inter-wrapper and the other one models 19 subchannel assemblies with a heat conduction model of the wrapper/inter-wrapper. In both models subchannel simulations are run as a MultiApp and heat conduction as the main app. The wrapper material is steel and the inter-wrapper is liquid sodium. Both materials are treated as solids and only heat conduction is considered. In both cases the central assembly has twice the power of its neighbors.

Coupling scheme

The coupling approach is that of domain decomposition. The main app calculates the temperature field in the wrapper and inter-wrapper regions. Based on the temperature field it calculates the heat-flux $var element = document.getElementById("moose-equation-25317d38-9696-4889-967b-42e880e39735");katex.render("(W/m^2)", element, {displayMode:false,throwOnError:false});$ on the inner surface of the wrapper (inner surface of the hexagonal ducts) and transfers that information to the equivalent duct model of SCM. The calculated heat-flux by the heat conduction model for the 7 assemblies case is shown in Figure 3

Figure 3: Heat-flux on the inner ducts for the 7 assembly case. Heat conduction model.

The heat-flux transfered to the duct model of SCM is shown in Figure 4.

Figure 4: Heat-flux on the inner ducts for the 7 assembly case. SCM model.

We show in Figure 4 two out of the seven SCM inner ducts. The center inner duct and one neighbor. In the 7 assembly case, the central assembly has twice the power of its neighbors. For this reason, heat flows from the high temperature, high power assembly through the wrapper/inter-wrapper domain into the neighboring colder, lower power assemblies.

SCM calculates the duct temperature at the inner duct surface and that information is transferred to the heat conduction model to be used as a boundary condition. All other inter-wrapper/wrapper surfaces assume a Neumann boundary condition. The SCM calculated duct temperature and the transferred BC are shown in Figure 5 and Figure 6.

Figure 5: Duct temperature calculated in two of the seven SCM assemblies.

Figure 6: Duct temperature transfered to the wrapper/inter-wrapper model.

Energy conservation

There is a postprocessor defined in the main app of type: "ADSideDiffusiveFluxIntegral" for each inner duct that calculates the total heat that flows through each duct.

[Postprocessors]
  # Total diffusive heat loss through central duct and the rest
  [center_duct_heat_loss_00]
    type = ADSideDiffusiveFluxIntegral
    variable = T_wrapper
    boundary = prsb_interface_00
    diffusivity = ${k_wrapper}
    execute_on = 'timestep_end'
  []
  [center_duct_heat_loss_01]
    type = ADSideDiffusiveFluxIntegral
    variable = T_wrapper
    boundary = prsb_interface_01
    diffusivity = ${k_wrapper}
    execute_on = 'timestep_end'
  []
  [center_duct_heat_loss_02]
    type = ADSideDiffusiveFluxIntegral
    variable = T_wrapper
    boundary = prsb_interface_02
    diffusivity = ${k_wrapper}
    execute_on = 'timestep_end'
  []
  [center_duct_heat_loss_03]
    type = ADSideDiffusiveFluxIntegral
    variable = T_wrapper
    boundary = prsb_interface_03
    diffusivity = ${k_wrapper}
    execute_on = 'timestep_end'
  []
  [center_duct_heat_loss_04]
    type = ADSideDiffusiveFluxIntegral
    variable = T_wrapper
    boundary = prsb_interface_04
    diffusivity = ${k_wrapper}
    execute_on = 'timestep_end'
  []
  [center_duct_heat_loss_05]
    type = ADSideDiffusiveFluxIntegral
    variable = T_wrapper
    boundary = prsb_interface_05
    diffusivity = ${k_wrapper}
    execute_on = 'timestep_end'
  []
  [center_duct_heat_loss_06]
    type = ADSideDiffusiveFluxIntegral
    variable = T_wrapper
    boundary = prsb_interface_06
    diffusivity = ${k_wrapper}
    execute_on = 'timestep_end'
  []
[]

There is also a postprocessor defined in the sub apps of type: "SCMDuctHeatRatePostprocessor" that calculates the total heat that flows through the SCM duct model.

[Postprocessors]
  [temp-1]
    type = SubChannelPointValue
    variable = T
    index = 527
    execute_on = 'TIMESTEP_END'
    height = 4.8
  []
  [temp-2]
    type = SubChannelPointValue
    variable = T
    index = 240
    execute_on = 'TIMESTEP_END'
    height = 4.8
  []
  [temp-3]
    type = SubChannelPointValue
    variable = T
    index = 112
    execute_on = 'TIMESTEP_END'
    height = 4.8
  []
  [temp-4]
    type = SubChannelPointValue
    variable = T
    index = 4
    execute_on = 'TIMESTEP_END'
    height = 4.8
  []
  [temp-5]
    type = SubChannelPointValue
    variable = T
    index = 1
    execute_on = 'TIMESTEP_END'
    height = 4.8
  []
  [temp-6]
    type = SubChannelPointValue
    variable = T
    index = 75
    execute_on = 'TIMESTEP_END'
    height = 4.8
  []
  [temp-7]
    type = SubChannelPointValue
    variable = T
    index = 258
    execute_on = 'TIMESTEP_END'
    height = 4.8
  []
  [temp-8]
    type = SubChannelPointValue
    variable = T
    index = 497
    execute_on = 'TIMESTEP_END'
    height = 4.8
  []
  ####### Assembly pressure drop
  [DP_SubchannelDelta]
    type = SubChannelDelta
    variable = P
    execute_on = 'TIMESTEP_END'
  []
  ##### Average temperature on plane
  [Mean_Temp]
    type = SCMPlanarMean
    variable = T
    height = 4.8
    execute_on = 'TIMESTEP_END'
  []
  #### Duct contribution to total power that goes into/out of the fluid
  [Total_Net_Power_Through_Duct]
    type = SCMDuctHeatRatePostprocessor
    execute_on ='transfer'
  []
  #### Total net power that goes into/out of the fluid
  [Total_coolant_power]
    type = SCMTHPowerPostprocessor
  []
[]

The "MultiAppGeneralFieldNearestLocationTransfer" uses these postprocessors as inputs to the parameters: "from_postprocessors_to_be_preserved, to_postprocessors_to_be_preserved" to ensure that energy is conserved during the transfer.

[Transfers]
  [q_duct]
    # send duct heat flux from heat conduction solve, to SCM
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = subchannel
    source_variable = duct_heat_flux
    variable = duct_heat_flux
    greedy_search = true
    from_boundaries = 'prsb_interface_00 prsb_interface_01  prsb_interface_02 prsb_interface_03
                       prsb_interface_04 prsb_interface_05 prsb_interface_06'
    execute_on = 'timestep_end'
    from_postprocessors_to_be_preserved = 'center_duct_heat_loss_00 center_duct_heat_loss_01 center_duct_heat_loss_02 center_duct_heat_loss_03 center_duct_heat_loss_04 center_duct_heat_loss_05 center_duct_heat_loss_06'
    to_postprocessors_to_be_preserved = 'Total_Net_Power_Through_Duct'
    allow_skipped_adjustment = true
  []

  [T_duct]
    # retrieve Tduct from SCM solve
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = subchannel
    source_variable = Tduct
    variable = duct_surface_temperature
    greedy_search = true
    to_boundaries = 'prsb_interface_00 prsb_interface_01  prsb_interface_02 prsb_interface_03
                     prsb_interface_04 prsb_interface_05 prsb_interface_06'
    execute_on = 'timestep_end'
  []
[]

Input files

The input files are the following:

7 assemblies

Input file to creates the wrapper/interwrapper model geometry. This only needs to be run once to generate the mesh.

# a MOOSE mesh for 7 ABR assemblies
# sqrt(3) / 2 is by how much flat to flat is smaller than corner to corner
f = '${fparse sqrt(3) / 2}'

# units are cm - do not forget to convert to meter
outer_duct_out = 15.8123 # outer size of the hexagonal duct (side to side)
outer_duct_in = 15.0191 # flat to flat
inter_wrapper_width = 0.4348
height = 480.2

# discretization
n_ax = 50
ns = 10
duct_intervals_perishperic = '4 4'

[Mesh]
  [XX00]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '12'
    polygon_size = '${fparse outer_duct_out / 2 + inter_wrapper_width / 2}'
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall inter_wrapper'
    outward_interface_boundary_names = 'wall_in_00 wall_out_00'
    interface_boundary_id_shift = 100
  []

  [XX01]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '13'
    polygon_size = '${fparse outer_duct_out / 2 + inter_wrapper_width / 2}'
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_01 wall_out_01'
    interface_boundary_id_shift = 200
  []

  [XX02]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '14'
    polygon_size = '${fparse outer_duct_out / 2 + inter_wrapper_width / 2}'
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_02 wall_out_02'
    interface_boundary_id_shift = 300
  []

  [XX03]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '15'
    polygon_size = '${fparse outer_duct_out / 2 + inter_wrapper_width / 2}'
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_03 wall_out_03'
    interface_boundary_id_shift = 400
  []

  [XX04]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '16'
    polygon_size = '${fparse outer_duct_out / 2 + inter_wrapper_width / 2}'
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_04 wall_out_04'
    interface_boundary_id_shift = 500
  []

  [XX05]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '17'
    polygon_size = '${fparse outer_duct_out / 2 + inter_wrapper_width / 2}'
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_05 wall_out_05'
    interface_boundary_id_shift = 600
  []

  [XX06]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '18'
    polygon_size = '${fparse outer_duct_out / 2 + inter_wrapper_width / 2}'
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_06 wall_out_06'
    interface_boundary_id_shift = 700
  []

  [pattern]
    type = PatternedHexMeshGenerator
    inputs = 'XX00 XX01 XX02 XX03 XX04 XX05 XX06'
    pattern = '5 4 ;
              6 0 3 ;
               1 2'
    pattern_boundary = none
  []

  [extrude]
    type = AdvancedExtruderGenerator
    direction = '0 0 1'
    input = pattern
    heights = '${height}'
    num_layers = '${n_ax}'
  []

  [rename]
    type = RenameBlockGenerator
    input = extrude
    old_block = '12             13             14             15             16             17             18'
    new_block = 'porous_flow_00 porous_flow_01 porous_flow_02 porous_flow_03 porous_flow_04 porous_flow_05 porous_flow_06'
  []

  [inlet_interwrapper]
    type = ParsedGenerateSideset
    input = rename
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'inter_wrapper duct_wall'
    normal = '0 0 -1'
    new_sideset_name = inlet_interwrapper
  []

  [inlet_porous_flow_hfd]
    type = ParsedGenerateSideset
    input = inlet_interwrapper
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_01
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_01
  []

  [inlet_porous_flow_p]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_hfd
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_02
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_02
  []

  [inlet_porous_flow_d1]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_p
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_03
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_03
  []

  [inlet_porous_flow_d2]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_d1
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_04
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_04
  []

  [inlet_porous_flow_k011]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_d2
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_05
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_05
  []

  [inlet_porous_flow_x402]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_k011
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_06
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_06
  []

  [inlet_central_assembly]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_x402
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_00'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_00
  []

  [outlet_interwrapper]
    type = ParsedGenerateSideset
    input = inlet_central_assembly
    included_subdomains = 'inter_wrapper duct_wall'
    combinatorial_geometry = 'abs(z - ${fparse height}) < 1e-6'
    normal = '0 0 1'
    new_sideset_name = outlet_interwrapper
  []

  [outlet_porous_flow]
    type = ParsedGenerateSideset
    input = outlet_interwrapper
    included_subdomains = 'porous_flow_01 porous_flow_02 porous_flow_03 porous_flow_04 porous_flow_05 porous_flow_06'
    combinatorial_geometry = 'abs(z - ${fparse height}) < 1e-6'
    normal = '0 0 1'
    new_sideset_name = outlet_porous_flow
  []

  [outlet_central_assembly]
    type = ParsedGenerateSideset
    input = outlet_porous_flow
    included_subdomains = 'porous_flow_00'
    combinatorial_geometry = 'abs(z - ${fparse height}) < 1e-6'
    normal = '0 0 1'
    new_sideset_name = outlet_porous_flow_00
  []

  [rotate]
    type = TransformGenerator
    input = outlet_central_assembly
    transform = ROTATE
    vector_value = '0 0 0'
  []

  # turn into meters
  [scale]
    type = TransformGenerator
    vector_value = '0.01 0.01 0.01'
    transform = SCALE
    input = rotate
  []

  [new_wall_boundary_00]
    type = SideSetsBetweenSubdomainsGenerator
    input = scale
    new_boundary = 'prsb_interface_00'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_00'
  []

  [new_wall_boundary_01]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_00
    new_boundary = 'prsb_interface_01'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_01'
  []

  [new_wall_boundary_02]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_01
    new_boundary = 'prsb_interface_02'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_02'
  []

  [new_wall_boundary_03]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_02
    new_boundary = 'prsb_interface_03'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_03'
  []

  [new_wall_boundary_04]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_03
    new_boundary = 'prsb_interface_04'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_04'
  []

  [new_wall_boundary_05]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_04
    new_boundary = 'prsb_interface_05'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_05'
  []

  [new_wall_boundary_06]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_05
    new_boundary = 'prsb_interface_06'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_06'
  []

  [delete_assemblies]
    type = BlockDeletionGenerator
    block = 'porous_flow_00 porous_flow_01 porous_flow_02 porous_flow_03 porous_flow_04 porous_flow_05 porous_flow_06'
    input = 'new_wall_boundary_06'
  []
[]

File that runs the main App, the heat conduction solve in the wrapper/interwrapper.

# This heat conduction model (main app) drives a multiapp setup where subchannel hexagonal ducts (solid wrappers)
# and an interwrapper (liquid sodium interwrapper between the ducts) are modeled and the temperature field is retrieved.
# This model treats the liquid sodium interwrapper as a solid and only heat conduction is considered.
# The calculated inner duct_wall axial heat flux [W/m] is transfered to 7 SCM models that run in parallel as multiapps.
# The method of coupling is domain decomposition:
# As mentioned above the heat conduction solution sends linear heat rate calculated on the inner duct surfaces to SCM
# and retrieves inner duct surface temperature from SCM. All other surfaces of the model have a Neuman BC.
##################### Instructions ###########################
# Create mesh file: abr_7assemblies_in.e
# run with : mpiexec -n N location-of-executable -i HC_master.i --keep-cout -w
# Units are SI

#### physics parameters
inlet_temperature = 653.15

####  Fluid properties of Sodium
####  Density #####
# A12 = 1.00423e3
# A13 = -0.21390
# A14 = -1.1046e-5
# rho = ${fparse A12 + A13 * inlet_temperature + A14 * inlet_temperature * inlet_temperature}
# #### Viscosity
# A52 = 3.6522e-5
# A53 = 0.16626
# A54 = -4.56877e1
# A55 = 2.8733e4
# mu = ${fparse A52 + A53 / inlet_temperature + A54 / inlet_temperature / inlet_temperature +
#         A55 / (inlet_temperature * inlet_temperature * inlet_temperature)}
#### Specific heat at constant pressure
# A28 = 7.3898e5
# A29 = 3.154e5
# A30 = 1.1340e3
# A31 = -2.2153e-1
# A32 = 1.1156e-4
# dt = ${fparse 2503.3 - inlet_temperature}
# cp = ${fparse A28 / dt / dt + A29 / dt + A30 + A31 * dt + A32 * dt * dt}
#### Heat conduction coefficient
A48 = 1.1045e2
A49 = -6.5112e-2
A50 = 1.5430e-5
A51 = -2.4617e-9
k_sodium = '${fparse A48 + A49 * inlet_temperature + A50 * inlet_temperature * inlet_temperature +
        A51 * inlet_temperature * inlet_temperature * inlet_temperature}'

# wrapper properties
k_wrapper = 15

wrapper_blocks = 'duct_wall'
inter_wrapper_blocks = 'inter_wrapper'

[Mesh]
  [fmesh]
    type = FileMeshGenerator
    file = 'abr_7assemblies_in.e'
  []
  coord_type = XYZ
[]

[Materials]
  [wall_heat_conductor]
    type = HeatConductionMaterial
    thermal_conductivity = ${k_wrapper}
    block = ${wrapper_blocks}
  []
  [sodium_heat_conductor]
    type = HeatConductionMaterial
    thermal_conductivity = ${k_sodium}
    block = ${inter_wrapper_blocks}
  []
[]

[Variables]
  [T_wrapper]
    order = FIRST
    family = LAGRANGE
  []
[]

[Kernels]
  [heat_conduction]
    type = HeatConduction
    variable = T_wrapper
    block = '${wrapper_blocks} ${inter_wrapper_blocks}'
  []
[]

[BCs]
  [T_wrapper_inside_wall]
    type = MatchedValueBC
    variable = T_wrapper
    boundary = 'prsb_interface_00 prsb_interface_01 prsb_interface_02 prsb_interface_03
                prsb_interface_04 prsb_interface_05 prsb_interface_06'
    v = duct_surface_temperature
  []
  [outside_bc]
   type = NeumannBC
   variable = T_wrapper
   boundary = '10000 inlet_interwrapper outlet_interwrapper'
  []
[]

[AuxVariables]
  [duct_surface_temperature]
    initial_condition = ${inlet_temperature}
  []

  [duct_heat_flux]
    order = CONSTANT
    family = MONOMIAL
    block = ${wrapper_blocks}
    initial_condition = 0
  []
[]

[AuxKernels]
  [HEAT_FLUX]
    type = DiffusionFluxAux
    diffusivity = ${k_wrapper}
    variable = duct_heat_flux
    diffusion_variable = T_wrapper
    component = normal
    boundary = 'prsb_interface_00 prsb_interface_01 prsb_interface_02 prsb_interface_03
                prsb_interface_04 prsb_interface_05 prsb_interface_06'
    execute_on = 'initial timestep_end'
  []
[]

[Executioner]
  type = Steady
  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  fixed_point_max_its = 5
  fixed_point_min_its = 4
  fixed_point_rel_tol = 1e-3
  fixed_point_abs_tol = 1e-3

  [Quadrature]
    order = THIRD
    side_order = FOURTH
  []
[]

[Postprocessors]
  # Total diffusive heat loss through central duct and the rest
  [center_duct_heat_loss_00]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_00
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_01]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_01
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_02]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_02
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_03]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_03
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_04]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_04
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_05]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_05
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_06]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_06
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
[]

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

################################################################################
######    A multiapp that couples heat conduction to subchannel via duct interface
################################################################################
[MultiApps]
  [subchannel]
    type = FullSolveMultiApp
    input_files = 'fuel_assembly_center.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i'
    execute_on = 'timestep_end'
    positions = '0         0           0
                 0.140704 -0.0812355   0
                 0.140704 0.0812355    0
                 0  0.162471           0
                 -0.140704 0.0812355   0
                 -0.140704 -0.0812355  0
                 0 -0.162471           0'
    max_procs_per_app = 1
    output_in_position = true
    bounding_box_padding = '0 0 0.1'
  []
[]

[Transfers]
  [q_duct] # send duct heat flux from heat conduction solve, to SCM
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = subchannel
    source_variable = duct_heat_flux
    variable = duct_heat_flux
    greedy_search = true
    from_boundaries = 'prsb_interface_00 prsb_interface_01  prsb_interface_02 prsb_interface_03
                       prsb_interface_04 prsb_interface_05 prsb_interface_06'
    execute_on = 'timestep_end'
    from_postprocessors_to_be_preserved = 'center_duct_heat_loss_00 center_duct_heat_loss_01 center_duct_heat_loss_02 center_duct_heat_loss_03 center_duct_heat_loss_04 center_duct_heat_loss_05 center_duct_heat_loss_06'
    to_postprocessors_to_be_preserved = 'Total_Net_Power_Through_Duct'
    allow_skipped_adjustment = true
  []

  [T_duct] # retrieve Tduct from SCM solve
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = subchannel
    source_variable = Tduct
    variable = duct_surface_temperature
    greedy_search = true
    to_boundaries = 'prsb_interface_00 prsb_interface_01  prsb_interface_02 prsb_interface_03
                     prsb_interface_04 prsb_interface_05 prsb_interface_06'
    execute_on = 'timestep_end'
  []
[]

Files that run the MultiApp.

###################################################
# Thermal-hydraulics parameters
###################################################
T_in = 653.15
P_out = 758423 # Pa
mdot_core = 5024
n_fuel_assemblies = 180
flow_area = 0.00650277
mass_flux_in = '${fparse mdot_core/n_fuel_assemblies/flow_area}' # kg/(m2.s)

###################################################
# Geometric parameters
###################################################
#f = ${fparse sqrt(3) / 2}

# units are cm - do not forget to convert to meter
scale_factor = 0.01
fuel_element_pitch = '${fparse 16.2471*scale_factor}'
inter_assembly_gap = '${fparse 0.4348*scale_factor}'
duct_thickness = '${fparse 0.3966*scale_factor}'
fuel_pin_pitch = '${fparse 0.8966*scale_factor}'
fuel_pin_diameter = '${fparse 0.7714*scale_factor}'
wire_z_spacing = '${fparse 20.43*scale_factor}' ###
wire_diameter = '${fparse 0.1307*scale_factor}' #check
n_rings = 10
length_entry_fuel = '${fparse 160.92*scale_factor}'
length_heated_fuel = '${fparse 85.82*scale_factor}'
length_outlet_fuel =  '${fparse 233.46*scale_factor}'
orifice_plate_height = '${fparse 5*scale_factor}'
duct_outside = '${fparse fuel_element_pitch - inter_assembly_gap}'
duct_inside = '${fparse duct_outside - 2 * duct_thickness}'
###################################################

[TriSubChannelMesh]
  [subchannel]
    type = SCMTriSubChannelMeshGenerator
    nrings = '${fparse n_rings}'
    n_cells = 50
    flat_to_flat = '${fparse duct_inside}'
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pin_diameter = '${fparse fuel_pin_diameter}'
    pitch = '${fparse fuel_pin_pitch}'
    dwire = '${fparse wire_diameter}'
    hwire = '${fparse wire_z_spacing}'
    spacer_z = '${fparse orifice_plate_height} ${fparse length_entry_fuel}'
    spacer_k = '0.5 0.5'
  []

  [fuel_pins]
    type = SCMTriPinMeshGenerator
    input = subchannel
    nrings = '${fparse n_rings}'
    n_cells = 50 #100
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pitch = '${fparse fuel_pin_pitch}'
  []

  [duct]
    type = SCMTriDuctMeshGenerator
    input = fuel_pins
    nrings = '${fparse n_rings}'
    n_cells = 50 #100
    flat_to_flat = '${fparse duct_inside}'
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pitch = '${fparse fuel_pin_pitch}'
  []
[]
# The function is evaluated with the origin z at unheated_length_entry
[Functions]
  [axial_heat_rate]
    type = ParsedFunction
    expression = '(pi/2)*sin(pi*z/L)'
    symbol_names = 'L'
    symbol_values = '${length_heated_fuel}'
  []
[]

[AuxVariables]
  [Tduct]
    block = duct
    initial_condition = 653.15
  []
[]

[FluidProperties]
  [sodium]
    type = PBSodiumFluidProperties
  []
[]

[SubChannel]
  type = TriSubChannel1PhaseProblem
  fp = sodium
  n_blocks = 1
  P_out = ${P_out}
  CT = 1.0
  compute_density = true
  compute_viscosity = true
  compute_power = true
  P_tol = 1.0e-4
  T_tol = 1.0e-3
  implicit = true
  segregated = false
  verbose_multiapps = true
  verbose_subchannel = true
[]

[ICs]
  [S_IC]
    type = SCMTriFlowAreaIC
    variable = S
  []

  [w_perim_IC]
    type = SCMTriWettedPerimIC
    variable = w_perim
  []

  [q_prime_IC]
    type = SCMTriPowerIC
    variable = q_prime
    power = 1000000.0 #W
    filename = "pin_power_profile_uni.txt"
    axial_heat_rate = axial_heat_rate
  []

  [T_ic]
    type = ConstantIC
    variable = T
    value = ${T_in}
  []

  [P_ic]
    type = ConstantIC
    variable = P
    value = 0.0 #always set to zero
  []

  [DP_ic]
    type = ConstantIC
    variable = DP
    value = 0.0
  []

  [Dpin_ic]
    type = ConstantIC
    variable = Dpin
    value = ${fuel_pin_diameter}
  []

  [Viscosity_ic]
    type = ViscosityIC
    variable = mu
    p = ${P_out}
    T = T
    fp = sodium
  []

  [rho_ic]
    type = RhoFromPressureTemperatureIC
    variable = rho
    p = ${P_out}
    T = T
    fp = sodium
  []

  [h_ic]
    type = SpecificEnthalpyFromPressureTemperatureIC
    variable = h
    p = ${P_out}
    T = T
    fp = sodium
  []

  [mdot_ic]
    type = ConstantIC
    variable = mdot
    value = 0.0
  []
[]

[AuxKernels]
  [P_out_bc]
    type = ConstantAux
    variable = P
    boundary = outlet
    value = ${P_out}
    execute_on = 'timestep_begin'
    block = subchannel
  []
  [T_in_bc]
    type = ConstantAux
    variable = T
    boundary = inlet
    value = ${T_in}
    execute_on = 'timestep_begin'
    block = subchannel
  []
  [mdot_in_bc]
    type = SCMMassFlowRateAux
    variable = mdot
    boundary = inlet
    area = S
    mass_flux = ${mass_flux_in}
    execute_on = 'timestep_begin'
    block = subchannel
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

!include SCM_output.i

[Executioner]
  type = Steady
[]

################################################################################
# A multiapp that projects data to a detailed mesh
################################################################################
[MultiApps]
  [viz]
    type = FullSolveMultiApp
    input_files = "3d.i"
    execute_on = "final"
    output_in_position = true
  []
[]

[Transfers]
  [subchannel_transfer]
      type = SCMSolutionTransfer
      to_multi_app = viz
      variable = 'mdot SumWij P DP h T rho mu S w_perim'
  []

  [pin_transfer]
      type = SCMPinSolutionTransfer
      to_multi_app = viz
      variable = 'Tpin Dpin q_prime'
  []
[]

###################################################
# Thermal-hydraulics parameters
###################################################
T_in = 653.15
P_out = 758423 # Pa
mdot_core = 5024
n_fuel_assemblies = 180
flow_area = 0.00650277
mass_flux_in = '${fparse mdot_core/n_fuel_assemblies/flow_area}' # kg/(m2.s)
###################################################
# Geometric parameters
###################################################
#f = ${fparse sqrt(3) / 2}

# units are cm - do not forget to convert to meter
scale_factor = 0.01
fuel_element_pitch = '${fparse 16.2471*scale_factor}'
inter_assembly_gap = '${fparse 0.4348*scale_factor}' # check
duct_thickness = '${fparse 0.3966*scale_factor}'
fuel_pin_pitch = '${fparse 0.8966*scale_factor}'
fuel_pin_diameter = '${fparse 0.7714*scale_factor}'
wire_z_spacing = '${fparse 20.43*scale_factor}' ###
wire_diameter = '${fparse 0.1307*scale_factor}' #check
n_rings = 10
length_entry_fuel = '${fparse 160.92*scale_factor}'
length_heated_fuel = '${fparse 85.82*scale_factor}'
length_outlet_fuel = '${fparse 233.46*scale_factor}'
orifice_plate_height = '${fparse 5*scale_factor}'
duct_outside = '${fparse fuel_element_pitch - inter_assembly_gap}'
duct_inside = '${fparse duct_outside - 2 * duct_thickness}'
###################################################

[TriSubChannelMesh]
  [subchannel]
    type = SCMTriSubChannelMeshGenerator
    nrings = '${fparse n_rings}'
    n_cells = 50
    flat_to_flat = '${fparse duct_inside}'
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pin_diameter = '${fparse fuel_pin_diameter}'
    pitch = '${fparse fuel_pin_pitch}'
    dwire = '${fparse wire_diameter}'
    hwire = '${fparse wire_z_spacing}'
    spacer_z = '${fparse orifice_plate_height} ${fparse length_entry_fuel}'
    spacer_k = '0.5 0.5'
  []

  [fuel_pins]
    type = SCMTriPinMeshGenerator
    input = subchannel
    nrings = '${fparse n_rings}'
    n_cells = 50
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pitch = '${fparse fuel_pin_pitch}'
  []

  [duct]
    type = SCMTriDuctMeshGenerator
    input = fuel_pins
    nrings = '${fparse n_rings}'
    n_cells = 50
    flat_to_flat = '${fparse duct_inside}'
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pitch = '${fparse fuel_pin_pitch}'
  []
[]
# The function is evaluated with the origin z at unheated_length_entry
[Functions]
  [axial_heat_rate]
    type = ParsedFunction
    expression = '(pi/2)*sin(pi*z/L)'
    symbol_names = 'L'
    symbol_values = '${length_heated_fuel}'
  []
[]

[AuxVariables]
  [Tduct]
    block = duct
    initial_condition = 653.15
  []
[]

[FluidProperties]
  [sodium]
    type = PBSodiumFluidProperties
  []
[]

[SubChannel]
  type = TriSubChannel1PhaseProblem
  fp = sodium
  n_blocks = 1
  P_out = ${P_out}
  CT = 1.0
  compute_density = true
  compute_viscosity = true
  compute_power = true
  P_tol = 1.0e-4
  T_tol = 1.0e-3
  implicit = true
  segregated = false
  verbose_multiapps = true
  verbose_subchannel = true
[]

[ICs]
  [S_IC]
    type = SCMTriFlowAreaIC
    variable = S
  []

  [w_perim_IC]
    type = SCMTriWettedPerimIC
    variable = w_perim
  []

  [q_prime_IC]
    type = SCMTriPowerIC
    variable = q_prime
    power = 2000000.0 #W
    filename = "pin_power_profile_uni.txt"
    axial_heat_rate = axial_heat_rate
  []

  [T_ic]
    type = ConstantIC
    variable = T
    value = ${T_in}
  []

  [P_ic]
    type = ConstantIC
    variable = P
    value = 0.0 #always set to zero
  []

  [DP_ic]
    type = ConstantIC
    variable = DP
    value = 0.0
  []

  [Dpin_ic]
    type = ConstantIC
    variable = Dpin
    value = ${fuel_pin_diameter}
  []

  [Viscosity_ic]
    type = ViscosityIC
    variable = mu
    p = ${P_out}
    T = T
    fp = sodium
  []

  [rho_ic]
    type = RhoFromPressureTemperatureIC
    variable = rho
    p = ${P_out}
    T = T
    fp = sodium
  []

  [h_ic]
    type = SpecificEnthalpyFromPressureTemperatureIC
    variable = h
    p = ${P_out}
    T = T
    fp = sodium
  []

  [mdot_ic]
    type = ConstantIC
    variable = mdot
    value = 0.0
  []
[]

[AuxKernels]
  [P_out_bc]
    type = ConstantAux
    variable = P
    boundary = outlet
    value = ${P_out}
    execute_on = 'timestep_begin'
    block = subchannel
  []
  [T_in_bc]
    type = ConstantAux
    variable = T
    boundary = inlet
    value = ${T_in}
    execute_on = 'timestep_begin'
    block = subchannel
  []
  [mdot_in_bc]
    type = SCMMassFlowRateAux
    variable = mdot
    boundary = inlet
    area = S
    mass_flux = ${mass_flux_in}
    execute_on = 'timestep_begin'
    block = subchannel
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

!include SCM_output.i

[Executioner]
  type = Steady
[]

################################################################################
# A multiapp that projects data to a detailed mesh
################################################################################
[MultiApps]
  [viz]
    type = FullSolveMultiApp
    input_files = "3d.i"
    execute_on = "final"
    output_in_position = true
  []
[]

[Transfers]
  [subchannel_transfer]
    type = SCMSolutionTransfer
    to_multi_app = viz
    variable = 'mdot SumWij P DP h T rho mu S w_perim'
  []

  [pin_transfer]
    type = SCMPinSolutionTransfer
    to_multi_app = viz
    variable = 'Tpin Dpin q_prime'
  []
[]

File that creates a 3D visualization model of the SCM solution. This is executed by each SCM simulation (grandchild app).

###################################################
# Geometric parameters
###################################################
#f = ${fparse sqrt(3) / 2}

# units are cm - do not forget to convert to meter
scale_factor = 0.01
fuel_element_pitch = '${fparse 16.2471*scale_factor}'
inter_assembly_gap = '${fparse 0.4348*scale_factor}'
duct_thickness = '${fparse 0.3966*scale_factor}'
fuel_pin_pitch = '${fparse 0.8966*scale_factor}'
fuel_pin_diameter = '${fparse 0.7714*scale_factor}'
length_entry_fuel = '${fparse 160.92*scale_factor}'
length_heated_fuel = '${fparse 85.82*scale_factor}'
length_outlet_fuel = '${fparse 233.46*scale_factor}'
duct_outside = ${fparse fuel_element_pitch - inter_assembly_gap}
duct_inside = ${fparse duct_outside -  2 * duct_thickness}
n_rings = 10

[Mesh]
  [subchannel]
    type = SCMDetailedTriSubChannelMeshGenerator
    nrings = '${fparse n_rings}'
    n_cells = 50 #100
    flat_to_flat = '${fparse duct_inside}'
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pin_diameter = '${fparse fuel_pin_diameter}'
    pitch = '${fparse fuel_pin_pitch}'
  []

  [fuel_pins]
    type = SCMDetailedTriPinMeshGenerator
    input = subchannel
    nrings = '${fparse n_rings}'
    n_cells = 50
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pitch = '${fparse fuel_pin_pitch}'
    pin_diameter = '${fparse fuel_pin_diameter}'
  []
[]

[AuxVariables]
  [mdot]
    block = subchannel
  []
  [SumWij]
    block = subchannel
  []
  [P]
    block = subchannel
  []
  [DP]
    block = subchannel
  []
  [h]
    block = subchannel
  []
  [T]
    block = subchannel
  []
  [rho]
    block = subchannel
  []
  [mu]
    block = subchannel
  []
  [S]
    block = subchannel
  []
  [w_perim]
    block = subchannel
  []
  [q_prime]
    block = fuel_pins
  []
  [Tpin]
    block = fuel_pins
  []
  [Dpin]
    block = fuel_pins
  []
[]

[Problem]
  type = NoSolveProblem
[]

[Outputs]
  exodus = true
[]

[Executioner]
  type = Steady
[]

19 assemblies

Input file to creates the wrapper/interwrapper model geometry. This only needs to be run once to generate the mesh.

# a MOOSE mesh for 19 ABR assemblies
# sqrt(3) / 2 is by how much flat to flat is smaller than corner to corner
f = ${fparse sqrt(3) / 2}

# units are cm - do not forget to convert to meter
outer_duct_out = 15.8123 # outer size of the hexagonal duct (side to side)
outer_duct_in = 15.0191 # flat to flat
inter_wrapper_width = 0.4348
height = 480.2

# discretization
n_ax = 50
ns = 10
duct_intervals_perishperic = '4 4'

[Mesh]
  [XX00]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '12'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall inter_wrapper'
    outward_interface_boundary_names = 'wall_in_00 wall_out_00'
    interface_boundary_id_shift = 100
  []

  [XX01]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '13'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_01 wall_out_01'
    interface_boundary_id_shift = 200
  []

  [XX02]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '14'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_02 wall_out_02'
    interface_boundary_id_shift = 300
  []

  [XX03]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '15'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_03 wall_out_03'
    interface_boundary_id_shift = 400
  []

  [XX04]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '16'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_04 wall_out_04'
    interface_boundary_id_shift = 500
  []

  [XX05]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '17'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_05 wall_out_05'
    interface_boundary_id_shift = 600
  []

  [XX06]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '18'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_06 wall_out_06'
    interface_boundary_id_shift = 700
  []

  [XX07]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '19'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_07 wall_out_07'
    interface_boundary_id_shift = 800
  []

  [XX08]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '20'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_08 wall_out_08'
    interface_boundary_id_shift = 900
  []

  [XX09]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '21'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_09 wall_out_09'
    interface_boundary_id_shift = 1000
  []

  [XX10]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '22'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_10 wall_out_10'
    interface_boundary_id_shift = 1100
  []

  [XX11]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '23'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_11 wall_out_11'
    interface_boundary_id_shift = 1200
  []

  [XX12]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '24'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_12 wall_out_12'
    interface_boundary_id_shift = 1300
  []

  [XX13]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '25'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_13 wall_out_13'
    interface_boundary_id_shift = 1400
  []

  [XX14]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '26'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_14 wall_out_14'
    interface_boundary_id_shift = 1500
  []

  [XX15]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '27'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_15 wall_out_15'
    interface_boundary_id_shift = 1600
  []

  [XX16]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '28'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_16 wall_out_16'
    interface_boundary_id_shift = 1700
  []

  [XX17]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '29'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_17 wall_out_17'
    interface_boundary_id_shift = 1800
  []

  [XX18]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '30'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_18 wall_out_18'
    interface_boundary_id_shift = 1900
  []

  [pattern]
    type = PatternedHexMeshGenerator
    inputs = 'XX00 XX01 XX02 XX03 XX04 XX05 XX06 XX07 XX08 XX09
              XX10 XX11 XX12 XX13 XX14 XX15 XX16 XX17 XX18'
    pattern =
            ' 15 14 13;
             16 5 4 12;
            17 6 0 3 11;
             18 1 2 10;
               7 8 9'
    pattern_boundary = none
  []

  [extrude]
    type = AdvancedExtruderGenerator
    direction = '0 0 1'
    input = pattern
    heights = '${height}'
    num_layers = '${n_ax}'
  []

  [rename]
    type = RenameBlockGenerator
    input = extrude
    old_block = '12             13              14            15             16             17             18
                                19              20            21             22             23             24
                                25              26            27             28             29             30'
    new_block = 'porous_flow_00 porous_flow_01 porous_flow_02 porous_flow_03 porous_flow_04 porous_flow_05 porous_flow_06
                                porous_flow_07 porous_flow_08 porous_flow_09 porous_flow_10 porous_flow_11 porous_flow_12
                                porous_flow_13 porous_flow_14 porous_flow_15 porous_flow_16 porous_flow_17 porous_flow_18'
  []

  [inlet_interwrapper]
    type = ParsedGenerateSideset
    input = rename
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'inter_wrapper duct_wall'
    normal = '0 0 -1'
    new_sideset_name = inlet_interwrapper
  []

  [inlet_porous_flow_hfd]
    type = ParsedGenerateSideset
    input = inlet_interwrapper
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_01
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_01
  []

  [inlet_porous_flow_p]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_hfd
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_02
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_02
  []

  [inlet_porous_flow_d1]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_p
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_03
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_03
  []

  [inlet_porous_flow_d2]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_d1
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_04
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_04
  []

  [inlet_porous_flow_k011]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_d2
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_05
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_05
  []

  [inlet_porous_flow_x402]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_k011
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_06
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_06
  []

  [inlet_central_assembly]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_x402
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_00'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_00
  []

  [inlet_porous_flow_07]
    type = ParsedGenerateSideset
    input = inlet_central_assembly
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_07'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_07
  []

  [inlet_porous_flow_08]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_07
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_08'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_08
  []

  [inlet_porous_flow_09]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_08
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_09'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_09
  []

  [inlet_porous_flow_10]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_09
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_10'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_10
  []

  [inlet_porous_flow_11]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_10
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_11'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_11
  []

  [inlet_porous_flow_12]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_11
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_12'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_12
  []

  [inlet_porous_flow_13]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_12
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_13'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_13
  []

  [inlet_porous_flow_14]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_13
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_14'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_14
  []

  [inlet_porous_flow_15]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_14
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_15'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_15
  []

  [inlet_porous_flow_16]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_15
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_16'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_16
  []

  [inlet_porous_flow_17]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_16
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_17'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_17
  []

  [inlet_porous_flow_18]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_17
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_18'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_18
  []

  [outlet_interwrapper]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_18 # outlet_interwall
    included_subdomains = 'inter_wrapper duct_wall'
    combinatorial_geometry = 'abs(z - ${fparse height}) < 1e-6'
    normal = '0 0 1'
    new_sideset_name = outlet_interwrapper
  []

  [outlet_porous_flow]
    type = ParsedGenerateSideset
    input = outlet_interwrapper
    included_subdomains = 'porous_flow_01 porous_flow_02 porous_flow_03 porous_flow_04 porous_flow_05 porous_flow_06 porous_flow_07 porous_flow_08 porous_flow_09 porous_flow_10 porous_flow_11 porous_flow_12 porous_flow_13 porous_flow_14 porous_flow_15 porous_flow_16 porous_flow_17 porous_flow_18'
    combinatorial_geometry = 'abs(z - ${fparse height}) < 1e-6'
    normal = '0 0 1'
    new_sideset_name = outlet_porous_flow
  []

  [outlet_central_assembly]
    type = ParsedGenerateSideset
    input = outlet_porous_flow
    included_subdomains = 'porous_flow_00'
    combinatorial_geometry = 'abs(z - ${fparse height}) < 1e-6'
    normal = '0 0 1'
    new_sideset_name = outlet_porous_flow_00
  []

  [rotate]
    type = TransformGenerator
    input = outlet_central_assembly
    transform = ROTATE
    vector_value = '0 0 0'
  []

  # turn into meters
  [scale]
    type = TransformGenerator
    vector_value = '0.01 0.01 0.01'
    transform = SCALE
    input = rotate
  []

  [new_wall_boundary_00]
    type = SideSetsBetweenSubdomainsGenerator
    input = scale
    new_boundary = 'prsb_interface_00'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_00'
  []

  [new_wall_boundary_01]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_00
    new_boundary = 'prsb_interface_01'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_01'
  []

  [new_wall_boundary_02]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_01
    new_boundary = 'prsb_interface_02'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_02'
  []

  [new_wall_boundary_03]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_02
    new_boundary = 'prsb_interface_03'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_03'
  []

  [new_wall_boundary_04]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_03
    new_boundary = 'prsb_interface_04'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_04'
  []

  [new_wall_boundary_05]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_04
    new_boundary = 'prsb_interface_05'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_05'
  []

  [new_wall_boundary_06]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_05
    new_boundary = 'prsb_interface_06'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_06'
  []

  [new_wall_boundary_07]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_06
    new_boundary = 'prsb_interface_07'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_07'
  []

  [new_wall_boundary_08]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_07
    new_boundary = 'prsb_interface_08'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_08'
  []

  [new_wall_boundary_09]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_08
    new_boundary = 'prsb_interface_09'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_09'
  []

  [new_wall_boundary_10]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_09
    new_boundary = 'prsb_interface_10'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_10'
  []

  [new_wall_boundary_11]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_10
    new_boundary = 'prsb_interface_11'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_11'
  []

  [new_wall_boundary_12]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_11
    new_boundary = 'prsb_interface_12'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_12'
  []

  [new_wall_boundary_13]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_12
    new_boundary = 'prsb_interface_13'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_13'
  []

  [new_wall_boundary_14]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_13
    new_boundary = 'prsb_interface_14'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_14'
  []

  [new_wall_boundary_15]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_14
    new_boundary = 'prsb_interface_15'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_15'
  []

  [new_wall_boundary_16]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_15
    new_boundary = 'prsb_interface_16'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_16'
  []

  [new_wall_boundary_17]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_16
    new_boundary = 'prsb_interface_17'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_17'
  []

  [new_wall_boundary_18]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_17
    new_boundary = 'prsb_interface_18'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_18'
  []

  [delete_assemblies]
    type = BlockDeletionGenerator
    block = 'porous_flow_00 porous_flow_01 porous_flow_02 porous_flow_03 porous_flow_04 porous_flow_05 porous_flow_06 porous_flow_07 porous_flow_08 porous_flow_09
             porous_flow_10 porous_flow_11 porous_flow_12 porous_flow_13 porous_flow_14 porous_flow_15 porous_flow_16 porous_flow_17 porous_flow_18'
    input = 'new_wall_boundary_18'
  []
[]

File that runs the main App, the heat conduction solve in the wrapper/interwrapper.

# This heat conduction model (main app) drives a multiapp setup where subchannel hexagonal ducts (solid wrappers)
# and an interwrapper (liquid sodium interwrapper between the ducts) are modeled and the temperature field is retrieved.
# This model treats the liquid sodium interwrapper as a solid and only heat conduction is considered.
# The calculated inner duct_wall axial heat flux [W/m] is transfered to 19 SCM models that run in parallel as multiapps.
# The method of coupling is domain decomposition:
# As mentioned above the heat conduction solution sends linear heat rate calculated on the inner duct surfaces to SCM
# and retrieves inner duct surface temperature from SCM. All other surfaces of the model has a Neuman BC.
##################### Instructions ###########################
# Create mesh file: abr_19assemblies_in.e
# run with : mpiexec -n N location-of-executable -i HC_master.i --keep-cout -w
# Units are SI

# physics parameters
inlet_temperature = 653.15

## fluid properties
#####  Density #####
# A12 = 1.00423e3
# A13 = -0.21390
# A14 = -1.1046e-5
# rho = '${fparse A12 + A13 * inlet_temperature + A14 * inlet_temperature * inlet_temperature}'
# #### Viscosity
# A52 = 3.6522e-5
# A53 = 0.16626
# A54 = -4.56877e1
# A55 = 2.8733e4
# mu = '${fparse A52 + A53 / inlet_temperature + A54 / inlet_temperature / inlet_temperature +
#         A55 / (inlet_temperature * inlet_temperature * inlet_temperature)} '
# #### Specific heat at constant pressure
# A28 = 7.3898e5
# A29 = 3.154e5
# A30 = 1.1340e3
# A31 = -2.2153e-1
# A32 = 1.1156e-4
# dt = '${fparse 2503.3 - inlet_temperature}'
# cp = '${fparse A28 / dt / dt + A29 / dt + A30 + A31 * dt + A32 * dt * dt}'
#### Heat conduction coefficient
A48 = 1.1045e2
A49 = -6.5112e-2
A50 = 1.5430e-5
A51 = -2.4617e-9
k_sodium = '${fparse A48 + A49 * inlet_temperature + A50 * inlet_temperature * inlet_temperature +
        A51 * inlet_temperature * inlet_temperature * inlet_temperature} '

# wrapper properties
k_wrapper = 15

wrapper_blocks = 'duct_wall'
inter_wrapper_blocks = 'inter_wrapper'

[Mesh]
  [fmesh]
    type = FileMeshGenerator
    file = 'abr_19assemblies_in.e'
  []
  coord_type = XYZ
[]

[Materials]
  [wall_heat_conductor]
    type = HeatConductionMaterial
    thermal_conductivity = ${k_wrapper}
    block = ${wrapper_blocks}
  []
  [sodium_heat_conductor]
    type = HeatConductionMaterial
    thermal_conductivity = ${k_sodium}
    block = ${inter_wrapper_blocks}
  []
[]

[Variables]
  [T_wrapper]
    order = FIRST
    family = LAGRANGE
  []
[]

[Kernels]
  [heat_conduction]
    type = HeatConduction
    variable = T_wrapper
    block = '${wrapper_blocks} ${inter_wrapper_blocks}'
  []
[]

[BCs]
  [T_wrapper_inside_wall]
    type = MatchedValueBC
    variable = T_wrapper
    boundary = 'prsb_interface_00 prsb_interface_01 prsb_interface_02 prsb_interface_03
                prsb_interface_04 prsb_interface_05 prsb_interface_06 prsb_interface_07
                prsb_interface_08 prsb_interface_09 prsb_interface_10 prsb_interface_11
                prsb_interface_12 prsb_interface_13 prsb_interface_14 prsb_interface_15
                prsb_interface_16 prsb_interface_17 prsb_interface_18'
    v = duct_surface_temperature
  []
  [outside_bc]
    type = NeumannBC
    variable = T_wrapper
    boundary = '10000 inlet_interwrapper outlet_interwrapper'
  []
[]

[AuxVariables]
  [duct_surface_temperature]
    initial_condition = ${inlet_temperature}
  []

  [duct_heat_flux]
    order = CONSTANT
    family = MONOMIAL
    block = ${wrapper_blocks}
    initial_condition = 0
  []
[]

[AuxKernels]
  [HEAT_FLUX]
    type = DiffusionFluxAux
    diffusivity = ${k_wrapper}
    variable = duct_heat_flux
    diffusion_variable = T_wrapper
    component = normal
    boundary = 'prsb_interface_00 prsb_interface_01 prsb_interface_02 prsb_interface_03
                prsb_interface_04 prsb_interface_05 prsb_interface_06 prsb_interface_07
                prsb_interface_08 prsb_interface_09 prsb_interface_10 prsb_interface_11
                prsb_interface_12 prsb_interface_13 prsb_interface_14 prsb_interface_15
                prsb_interface_16 prsb_interface_17 prsb_interface_18'
    execute_on = 'initial timestep_end'
  []
[]

[Executioner]
  type = Steady
  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  fixed_point_max_its = 5
  fixed_point_min_its = 4
  fixed_point_rel_tol = 1e-3
  fixed_point_abs_tol = 1e-3

  [Quadrature]
    order = THIRD
    side_order = FOURTH
  []
[]

[Postprocessors]
  # Total diffusive heat loss through central duct and the rest
  [center_duct_heat_loss_00]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_00
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_01]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_01
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_02]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_02
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_03]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_03
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_04]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_04
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_05]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_05
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_06]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_06
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
    [center_duct_heat_loss_07]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_07
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_08]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_08
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_09]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_09
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_10]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_10
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_11]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_11
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_12]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_12
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_13]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_13
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
    [center_duct_heat_loss_14]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_14
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_15]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_15
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_16]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_16
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_17]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_17
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_18]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_18
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
[]

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

################################################################################
######    A multiapp that couples heat conduction to subchannel via duct interface
################################################################################

[MultiApps]
  [subchannel]
    type = FullSolveMultiApp
    input_files = 'fuel_assembly_center.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i
    fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i
    fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i
    fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i
    fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i'
    execute_on = 'timestep_end'
    positions = '0          0            0
                 0.140704  -0.0812355    0
                 0.140704   0.0812355    0
                 0          0.162471     0
                -0.140704   0.0812355    0
                -0.140704  -0.0812355    0
                 0         -0.162471     0
                 0.281408  -0.162471     0
                 0.281408   0            0
                 0.281408   0.162471     0
                 0.140704   0.2437065    0
                 0          0.324942     0
                -0.140704   0.2437065    0
                -0.281408   0.162471     0
                -0.281408   0            0
                -0.281408  -0.162471     0
                -0.140704  -0.2437065    0
                 0         -0.324942     0
                 0.140704  -0.2437065    0'
    max_procs_per_app = 1
    output_in_position = true
    bounding_box_padding = '0 0 0.1'
  []
[]

[Transfers]
  [q_duct] # send duct heat flux from heat conduction solve, to SCM
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = subchannel
    source_variable = duct_heat_flux
    variable = duct_heat_flux
    greedy_search = true
    from_boundaries = 'prsb_interface_00 prsb_interface_01 prsb_interface_02 prsb_interface_03
                prsb_interface_04 prsb_interface_05 prsb_interface_06 prsb_interface_07
                prsb_interface_08 prsb_interface_09 prsb_interface_10 prsb_interface_11
                prsb_interface_12 prsb_interface_13 prsb_interface_14 prsb_interface_15
                prsb_interface_16 prsb_interface_17 prsb_interface_18'
    execute_on = 'timestep_end'
    from_postprocessors_to_be_preserved = 'center_duct_heat_loss_00 center_duct_heat_loss_01 center_duct_heat_loss_02 center_duct_heat_loss_03 center_duct_heat_loss_04 center_duct_heat_loss_05 center_duct_heat_loss_06
    center_duct_heat_loss_07 center_duct_heat_loss_08 center_duct_heat_loss_09 center_duct_heat_loss_10 center_duct_heat_loss_11
    center_duct_heat_loss_12 center_duct_heat_loss_13 center_duct_heat_loss_14 center_duct_heat_loss_15 center_duct_heat_loss_16
    center_duct_heat_loss_17 center_duct_heat_loss_18'
    to_postprocessors_to_be_preserved = 'Total_Net_Power_Through_Duct'
    allow_skipped_adjustment = true
  []

  [T_duct] # retrieve Tduct from SCM solve
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = subchannel
    source_variable = Tduct
    variable = duct_surface_temperature
    greedy_search = true
    to_boundaries = 'prsb_interface_00 prsb_interface_01 prsb_interface_02 prsb_interface_03
                prsb_interface_04 prsb_interface_05 prsb_interface_06 prsb_interface_07
                prsb_interface_08 prsb_interface_09 prsb_interface_10 prsb_interface_11
                prsb_interface_12 prsb_interface_13 prsb_interface_14 prsb_interface_15
                prsb_interface_16 prsb_interface_17 prsb_interface_18'
    execute_on = 'timestep_end'
  []
[]

File for the subchannel simulation, executed as a MultiApp.

###################################################
# Thermal-hydraulics parameters
###################################################
T_in = 653.15
P_out = 758423 # Pa
mdot_core = 5024
n_fuel_assemblies = 180
flow_area = 0.00650277
mass_flux_in = '${fparse mdot_core/n_fuel_assemblies/flow_area}' # kg/(m2.s)

###################################################
# Geometric parameters
###################################################
#f = ${fparse sqrt(3) / 2}

# units are cm - do not forget to convert to meter
scale_factor = 0.01
fuel_element_pitch = '${fparse 16.2471*scale_factor}'
inter_assembly_gap = '${fparse 0.4348*scale_factor}' # check
duct_thickness = '${fparse 0.3966*scale_factor}'
fuel_pin_pitch = '${fparse 0.8966*scale_factor}'
fuel_pin_diameter = '${fparse 0.7714*scale_factor}'
wire_z_spacing = '${fparse 20.43*scale_factor}' ###
wire_diameter = '${fparse 0.1307*scale_factor}' #check
n_rings = 10
length_entry_fuel = '${fparse 160.92*scale_factor}'
length_heated_fuel = '${fparse 85.82*scale_factor}'
length_outlet_fuel = '${fparse 233.46*scale_factor}'
orifice_plate_height = '${fparse 5*scale_factor}'
duct_outside = '${fparse fuel_element_pitch - inter_assembly_gap}'
duct_inside = '${fparse duct_outside - 2 * duct_thickness}'
###################################################

[TriSubChannelMesh]
  [subchannel]
    type = SCMTriSubChannelMeshGenerator
    nrings = '${fparse n_rings}'
    n_cells = 50
    flat_to_flat = '${fparse duct_inside}'
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pin_diameter = '${fparse fuel_pin_diameter}'
    pitch = '${fparse fuel_pin_pitch}'
    dwire = '${fparse wire_diameter}'
    hwire = '${fparse wire_z_spacing}'
    spacer_z = '${fparse orifice_plate_height} ${fparse length_entry_fuel}'
    spacer_k = '0.5 0.5'
  []

  [fuel_pins]
    type = SCMTriPinMeshGenerator
    input = subchannel
    nrings = '${fparse n_rings}'
    n_cells = 50 #100
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pitch = '${fparse fuel_pin_pitch}'
  []

  [duct]
    type = SCMTriDuctMeshGenerator
    input = fuel_pins
    nrings = '${fparse n_rings}'
    n_cells = 50 #100
    flat_to_flat = '${fparse duct_inside}'
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pitch = '${fparse fuel_pin_pitch}'
  []
[]
# The function is evaluated with the origin z at unheated_length_entry
[Functions]
  [axial_heat_rate]
    type = ParsedFunction
    expression = '(pi/2)*sin(pi*z/L)'
    symbol_names = 'L'
    symbol_values = '${length_heated_fuel}'
  []
[]

[AuxVariables]
  [Tduct]
    block = duct
    initial_condition = 653.15
  []
[]

[FluidProperties]
  [sodium]
    type = PBSodiumFluidProperties
  []
[]

[SubChannel]
  type = TriSubChannel1PhaseProblem
  fp = sodium
  n_blocks = 1
  P_out = ${P_out}
  CT = 1.0
  compute_density = true
  compute_viscosity = true
  compute_power = true
  P_tol = 1.0e-4
  T_tol = 1.0e-3
  implicit = true
  segregated = false
  verbose_multiapps = true
  verbose_subchannel = true
[]

[ICs]
  [S_IC]
    type = SCMTriFlowAreaIC
    variable = S
  []

  [w_perim_IC]
    type = SCMTriWettedPerimIC
    variable = w_perim
  []

  [q_prime_IC]
    type = SCMTriPowerIC
    variable = q_prime
    power = 2000000.0 #W
    filename = "pin_power_profile_uni.txt"
    axial_heat_rate = axial_heat_rate
  []

  [T_ic]
    type = ConstantIC
    variable = T
    value = ${T_in}
  []

  [P_ic]
    type = ConstantIC
    variable = P
    value = 0.0 #always set to zero
  []

  [DP_ic]
    type = ConstantIC
    variable = DP
    value = 0.0
  []

  [Dpin_ic]
    type = ConstantIC
    variable = Dpin
    value = ${fuel_pin_diameter}
  []

  [Viscosity_ic]
    type = ViscosityIC
    variable = mu
    p = ${P_out}
    T = T
    fp = sodium
  []

  [rho_ic]
    type = RhoFromPressureTemperatureIC
    variable = rho
    p = ${P_out}
    T = T
    fp = sodium
  []

  [h_ic]
    type = SpecificEnthalpyFromPressureTemperatureIC
    variable = h
    p = ${P_out}
    T = T
    fp = sodium
  []

  [mdot_ic]
    type = ConstantIC
    variable = mdot
    value = 0.0
  []
[]

[AuxKernels]
  [P_out_bc]
    type = ConstantAux
    variable = P
    boundary = outlet
    value = ${P_out}
    execute_on = 'timestep_begin'
    block = subchannel
  []
  [T_in_bc]
    type = ConstantAux
    variable = T
    boundary = inlet
    value = ${T_in}
    execute_on = 'timestep_begin'
    block = subchannel
  []
  [mdot_in_bc]
    type = SCMMassFlowRateAux
    variable = mdot
    boundary = inlet
    area = S
    mass_flux = ${mass_flux_in}
    execute_on = 'timestep_begin'
    block = subchannel
  []
[]

[Outputs]
  csv = true
[]

!include SCM_output.i

[Executioner]
  type = Steady
[]

################################################################################
# A multiapp that projects data to a detailed mesh
################################################################################
[MultiApps]
  [viz]
    type = FullSolveMultiApp
    input_files = "3d.i"
    execute_on = "final"
    output_in_position = true
  []
[]

[Transfers]
  [subchannel_transfer]
    type = SCMSolutionTransfer
    to_multi_app = viz
    variable = 'mdot SumWij P DP h T rho mu S w_perim'
  []

  [pin_transfer]
    type = SCMPinSolutionTransfer
    to_multi_app = viz
    variable = 'Tpin Dpin q_prime'
  []
[]

File that creates a 3D visualization model of the SCM solution. This is executed by each SCM simulation (grandchild app).

###################################################
# Geometric parameters
###################################################
#f = ${fparse sqrt(3) / 2}

# units are cm - do not forget to convert to meter
scale_factor = 0.01
fuel_element_pitch = '${fparse 16.2471*scale_factor}'
inter_assembly_gap = '${fparse 0.4348*scale_factor}'
duct_thickness = '${fparse 0.3966*scale_factor}'
fuel_pin_pitch = '${fparse 0.8966*scale_factor}'
fuel_pin_diameter = '${fparse 0.7714*scale_factor}'
length_entry_fuel = '${fparse 160.92*scale_factor}'
length_heated_fuel = '${fparse 85.82*scale_factor}'
length_outlet_fuel = '${fparse 233.46*scale_factor}'
duct_outside = ${fparse fuel_element_pitch - inter_assembly_gap}
duct_inside = ${fparse duct_outside -  2 * duct_thickness}
n_rings = 10

[Mesh]
  [subchannel]
    type = SCMDetailedTriSubChannelMeshGenerator
    nrings = '${fparse n_rings}'
    n_cells = 50 #100
    flat_to_flat = '${fparse duct_inside}'
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pin_diameter = '${fparse fuel_pin_diameter}'
    pitch = '${fparse fuel_pin_pitch}'
  []

  [fuel_pins]
    type = SCMDetailedTriPinMeshGenerator
    input = subchannel
    nrings = '${fparse n_rings}'
    n_cells = 50
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pitch = '${fparse fuel_pin_pitch}'
    pin_diameter = '${fparse fuel_pin_diameter}'
  []
[]

[AuxVariables]
  [mdot]
    block = subchannel
  []
  [SumWij]
    block = subchannel
  []
  [P]
    block = subchannel
  []
  [DP]
    block = subchannel
  []
  [h]
    block = subchannel
  []
  [T]
    block = subchannel
  []
  [rho]
    block = subchannel
  []
  [mu]
    block = subchannel
  []
  [S]
    block = subchannel
  []
  [w_perim]
    block = subchannel
  []
  [q_prime]
    block = fuel_pins
  []
  [Tpin]
    block = fuel_pins
  []
  [Dpin]
    block = fuel_pins
  []
[]

[Problem]
  type = NoSolveProblem
[]

[Outputs]
  exodus = true
[]

[Executioner]
  type = Steady
[]

Results

The temperature field for the coupled simulation is shown in Figure 7 and Figure 8

Figure 7: Temperature field for the 7 assemblies example.

Figure 8: Temperature field for the 19 assemblies example.

In both examples the center assembly has twice the power of its neighbors. The wrapper/inter-wrapper region is hotter around the center assembly. In the 19 assembly case all the assemblies have the same power. The wrapper/inter-wrapper region has a uniform temperature profile far away from the central duct.

(sfr/subchannel/multiple_SCM_assemblies/7assemblies/HC_master.i)

# This heat conduction model (main app) drives a multiapp setup where subchannel hexagonal ducts (solid wrappers)
# and an interwrapper (liquid sodium interwrapper between the ducts) are modeled and the temperature field is retrieved.
# This model treats the liquid sodium interwrapper as a solid and only heat conduction is considered.
# The calculated inner duct_wall axial heat flux [W/m] is transfered to 7 SCM models that run in parallel as multiapps.
# The method of coupling is domain decomposition:
# As mentioned above the heat conduction solution sends linear heat rate calculated on the inner duct surfaces to SCM
# and retrieves inner duct surface temperature from SCM. All other surfaces of the model have a Neuman BC.
##################### Instructions ###########################
# Create mesh file: abr_7assemblies_in.e
# run with : mpiexec -n N location-of-executable -i HC_master.i --keep-cout -w
# Units are SI

#### physics parameters
inlet_temperature = 653.15

####  Fluid properties of Sodium
####  Density #####
# A12 = 1.00423e3
# A13 = -0.21390
# A14 = -1.1046e-5
# rho = ${fparse A12 + A13 * inlet_temperature + A14 * inlet_temperature * inlet_temperature}
# #### Viscosity
# A52 = 3.6522e-5
# A53 = 0.16626
# A54 = -4.56877e1
# A55 = 2.8733e4
# mu = ${fparse A52 + A53 / inlet_temperature + A54 / inlet_temperature / inlet_temperature +
#         A55 / (inlet_temperature * inlet_temperature * inlet_temperature)}
#### Specific heat at constant pressure
# A28 = 7.3898e5
# A29 = 3.154e5
# A30 = 1.1340e3
# A31 = -2.2153e-1
# A32 = 1.1156e-4
# dt = ${fparse 2503.3 - inlet_temperature}
# cp = ${fparse A28 / dt / dt + A29 / dt + A30 + A31 * dt + A32 * dt * dt}
#### Heat conduction coefficient
A48 = 1.1045e2
A49 = -6.5112e-2
A50 = 1.5430e-5
A51 = -2.4617e-9
k_sodium = '${fparse A48 + A49 * inlet_temperature + A50 * inlet_temperature * inlet_temperature +
        A51 * inlet_temperature * inlet_temperature * inlet_temperature}'

# wrapper properties
k_wrapper = 15

wrapper_blocks = 'duct_wall'
inter_wrapper_blocks = 'inter_wrapper'

[Mesh]
  [fmesh]
    type = FileMeshGenerator
    file = 'abr_7assemblies_in.e'
  []
  coord_type = XYZ
[]

[Materials]
  [wall_heat_conductor]
    type = HeatConductionMaterial
    thermal_conductivity = ${k_wrapper}
    block = ${wrapper_blocks}
  []
  [sodium_heat_conductor]
    type = HeatConductionMaterial
    thermal_conductivity = ${k_sodium}
    block = ${inter_wrapper_blocks}
  []
[]

[Variables]
  [T_wrapper]
    order = FIRST
    family = LAGRANGE
  []
[]


[Kernels]
  [heat_conduction]
    type = HeatConduction
    variable = T_wrapper
    block = '${wrapper_blocks} ${inter_wrapper_blocks}'
  []
[]

[BCs]
  [T_wrapper_inside_wall]
    type = MatchedValueBC
    variable = T_wrapper
    boundary = 'prsb_interface_00 prsb_interface_01 prsb_interface_02 prsb_interface_03
                prsb_interface_04 prsb_interface_05 prsb_interface_06'
    v = duct_surface_temperature
  []
  [outside_bc]
   type = NeumannBC
   variable = T_wrapper
   boundary = '10000 inlet_interwrapper outlet_interwrapper'
  []
[]

[AuxVariables]
  [duct_surface_temperature]
    initial_condition = ${inlet_temperature}
  []

  [duct_heat_flux]
    order = CONSTANT
    family = MONOMIAL
    block = ${wrapper_blocks}
    initial_condition = 0
  []
[]

[AuxKernels]
  [HEAT_FLUX]
    type = DiffusionFluxAux
    diffusivity = ${k_wrapper}
    variable = duct_heat_flux
    diffusion_variable = T_wrapper
    component = normal
    boundary = 'prsb_interface_00 prsb_interface_01 prsb_interface_02 prsb_interface_03
                prsb_interface_04 prsb_interface_05 prsb_interface_06'
    execute_on = 'initial timestep_end'
  []
[]

[Executioner]
  type = Steady
  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  fixed_point_max_its = 5
  fixed_point_min_its = 4
  fixed_point_rel_tol = 1e-3
  fixed_point_abs_tol = 1e-3

  [Quadrature]
    order = THIRD
    side_order = FOURTH
  []
[]

[Postprocessors]
  # Total diffusive heat loss through central duct and the rest
  [center_duct_heat_loss_00]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_00
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_01]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_01
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_02]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_02
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_03]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_03
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_04]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_04
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_05]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_05
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_06]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_06
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
[]

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

################################################################################
######    A multiapp that couples heat conduction to subchannel via duct interface
################################################################################
[MultiApps]
  [subchannel]
    type = FullSolveMultiApp
    input_files = 'fuel_assembly_center.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i'
    execute_on = 'timestep_end'
    positions = '0         0           0
                 0.140704 -0.0812355   0
                 0.140704 0.0812355    0
                 0  0.162471           0
                 -0.140704 0.0812355   0
                 -0.140704 -0.0812355  0
                 0 -0.162471           0'
    max_procs_per_app = 1
    output_in_position = true
    bounding_box_padding = '0 0 0.1'
  []
[]

[Transfers]
  [q_duct] # send duct heat flux from heat conduction solve, to SCM
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = subchannel
    source_variable = duct_heat_flux
    variable = duct_heat_flux
    greedy_search = true
    from_boundaries = 'prsb_interface_00 prsb_interface_01  prsb_interface_02 prsb_interface_03
                       prsb_interface_04 prsb_interface_05 prsb_interface_06'
    execute_on = 'timestep_end'
    from_postprocessors_to_be_preserved = 'center_duct_heat_loss_00 center_duct_heat_loss_01 center_duct_heat_loss_02 center_duct_heat_loss_03 center_duct_heat_loss_04 center_duct_heat_loss_05 center_duct_heat_loss_06'
    to_postprocessors_to_be_preserved = 'Total_Net_Power_Through_Duct'
    allow_skipped_adjustment = true
  []

  [T_duct] # retrieve Tduct from SCM solve
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = subchannel
    source_variable = Tduct
    variable = duct_surface_temperature
    greedy_search = true
    to_boundaries = 'prsb_interface_00 prsb_interface_01  prsb_interface_02 prsb_interface_03
                     prsb_interface_04 prsb_interface_05 prsb_interface_06'
    execute_on = 'timestep_end'
  []
[]

(sfr/subchannel/multiple_SCM_assemblies/7assemblies/SCM_output.i)

[Postprocessors]
  [temp-1]
    type = SubChannelPointValue
    variable = T
    index = 527
    execute_on = 'TIMESTEP_END'
    height = 4.8
  []
  [temp-2]
    type = SubChannelPointValue
    variable = T
    index = 240
    execute_on = 'TIMESTEP_END'
    height = 4.8
  []
  [temp-3]
    type = SubChannelPointValue
    variable = T
    index = 112
    execute_on = 'TIMESTEP_END'
    height = 4.8
  []
  [temp-4]
    type = SubChannelPointValue
    variable = T
    index = 4
    execute_on = 'TIMESTEP_END'
    height = 4.8
  []
  [temp-5]
    type = SubChannelPointValue
    variable = T
    index = 1
    execute_on = 'TIMESTEP_END'
    height = 4.8
  []
  [temp-6]
    type = SubChannelPointValue
    variable = T
    index = 75
    execute_on = 'TIMESTEP_END'
    height = 4.8
  []
  [temp-7]
    type = SubChannelPointValue
    variable = T
    index = 258
    execute_on = 'TIMESTEP_END'
    height = 4.8
  []
  [temp-8]
    type = SubChannelPointValue
    variable = T
    index = 497
    execute_on = 'TIMESTEP_END'
    height = 4.8
  []
  ####### Assembly pressure drop
  [DP_SubchannelDelta]
    type = SubChannelDelta
    variable = P
    execute_on = 'TIMESTEP_END'
  []
  ##### Average temperature on plane
  [Mean_Temp]
    type = SCMPlanarMean
    variable = T
    height = 4.8
    execute_on = 'TIMESTEP_END'
  []
  #### Duct contribution to total power that goes into/out of the fluid
  [Total_Net_Power_Through_Duct]
    type = SCMDuctHeatRatePostprocessor
    execute_on ='transfer'
  []
  #### Total net power that goes into/out of the fluid
  [Total_coolant_power]
    type = SCMTHPowerPostprocessor
  []
[]

(sfr/subchannel/multiple_SCM_assemblies/7assemblies/HC_master.i)

# This heat conduction model (main app) drives a multiapp setup where subchannel hexagonal ducts (solid wrappers)
# and an interwrapper (liquid sodium interwrapper between the ducts) are modeled and the temperature field is retrieved.
# This model treats the liquid sodium interwrapper as a solid and only heat conduction is considered.
# The calculated inner duct_wall axial heat flux [W/m] is transfered to 7 SCM models that run in parallel as multiapps.
# The method of coupling is domain decomposition:
# As mentioned above the heat conduction solution sends linear heat rate calculated on the inner duct surfaces to SCM
# and retrieves inner duct surface temperature from SCM. All other surfaces of the model have a Neuman BC.
##################### Instructions ###########################
# Create mesh file: abr_7assemblies_in.e
# run with : mpiexec -n N location-of-executable -i HC_master.i --keep-cout -w
# Units are SI

#### physics parameters
inlet_temperature = 653.15

####  Fluid properties of Sodium
####  Density #####
# A12 = 1.00423e3
# A13 = -0.21390
# A14 = -1.1046e-5
# rho = ${fparse A12 + A13 * inlet_temperature + A14 * inlet_temperature * inlet_temperature}
# #### Viscosity
# A52 = 3.6522e-5
# A53 = 0.16626
# A54 = -4.56877e1
# A55 = 2.8733e4
# mu = ${fparse A52 + A53 / inlet_temperature + A54 / inlet_temperature / inlet_temperature +
#         A55 / (inlet_temperature * inlet_temperature * inlet_temperature)}
#### Specific heat at constant pressure
# A28 = 7.3898e5
# A29 = 3.154e5
# A30 = 1.1340e3
# A31 = -2.2153e-1
# A32 = 1.1156e-4
# dt = ${fparse 2503.3 - inlet_temperature}
# cp = ${fparse A28 / dt / dt + A29 / dt + A30 + A31 * dt + A32 * dt * dt}
#### Heat conduction coefficient
A48 = 1.1045e2
A49 = -6.5112e-2
A50 = 1.5430e-5
A51 = -2.4617e-9
k_sodium = '${fparse A48 + A49 * inlet_temperature + A50 * inlet_temperature * inlet_temperature +
        A51 * inlet_temperature * inlet_temperature * inlet_temperature}'

# wrapper properties
k_wrapper = 15

wrapper_blocks = 'duct_wall'
inter_wrapper_blocks = 'inter_wrapper'

[Mesh]
  [fmesh]
    type = FileMeshGenerator
    file = 'abr_7assemblies_in.e'
  []
  coord_type = XYZ
[]

[Materials]
  [wall_heat_conductor]
    type = HeatConductionMaterial
    thermal_conductivity = ${k_wrapper}
    block = ${wrapper_blocks}
  []
  [sodium_heat_conductor]
    type = HeatConductionMaterial
    thermal_conductivity = ${k_sodium}
    block = ${inter_wrapper_blocks}
  []
[]

[Variables]
  [T_wrapper]
    order = FIRST
    family = LAGRANGE
  []
[]


[Kernels]
  [heat_conduction]
    type = HeatConduction
    variable = T_wrapper
    block = '${wrapper_blocks} ${inter_wrapper_blocks}'
  []
[]

[BCs]
  [T_wrapper_inside_wall]
    type = MatchedValueBC
    variable = T_wrapper
    boundary = 'prsb_interface_00 prsb_interface_01 prsb_interface_02 prsb_interface_03
                prsb_interface_04 prsb_interface_05 prsb_interface_06'
    v = duct_surface_temperature
  []
  [outside_bc]
   type = NeumannBC
   variable = T_wrapper
   boundary = '10000 inlet_interwrapper outlet_interwrapper'
  []
[]

[AuxVariables]
  [duct_surface_temperature]
    initial_condition = ${inlet_temperature}
  []

  [duct_heat_flux]
    order = CONSTANT
    family = MONOMIAL
    block = ${wrapper_blocks}
    initial_condition = 0
  []
[]

[AuxKernels]
  [HEAT_FLUX]
    type = DiffusionFluxAux
    diffusivity = ${k_wrapper}
    variable = duct_heat_flux
    diffusion_variable = T_wrapper
    component = normal
    boundary = 'prsb_interface_00 prsb_interface_01 prsb_interface_02 prsb_interface_03
                prsb_interface_04 prsb_interface_05 prsb_interface_06'
    execute_on = 'initial timestep_end'
  []
[]

[Executioner]
  type = Steady
  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  fixed_point_max_its = 5
  fixed_point_min_its = 4
  fixed_point_rel_tol = 1e-3
  fixed_point_abs_tol = 1e-3

  [Quadrature]
    order = THIRD
    side_order = FOURTH
  []
[]

[Postprocessors]
  # Total diffusive heat loss through central duct and the rest
  [center_duct_heat_loss_00]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_00
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_01]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_01
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_02]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_02
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_03]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_03
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_04]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_04
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_05]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_05
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_06]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_06
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
[]

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

################################################################################
######    A multiapp that couples heat conduction to subchannel via duct interface
################################################################################
[MultiApps]
  [subchannel]
    type = FullSolveMultiApp
    input_files = 'fuel_assembly_center.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i'
    execute_on = 'timestep_end'
    positions = '0         0           0
                 0.140704 -0.0812355   0
                 0.140704 0.0812355    0
                 0  0.162471           0
                 -0.140704 0.0812355   0
                 -0.140704 -0.0812355  0
                 0 -0.162471           0'
    max_procs_per_app = 1
    output_in_position = true
    bounding_box_padding = '0 0 0.1'
  []
[]

[Transfers]
  [q_duct] # send duct heat flux from heat conduction solve, to SCM
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = subchannel
    source_variable = duct_heat_flux
    variable = duct_heat_flux
    greedy_search = true
    from_boundaries = 'prsb_interface_00 prsb_interface_01  prsb_interface_02 prsb_interface_03
                       prsb_interface_04 prsb_interface_05 prsb_interface_06'
    execute_on = 'timestep_end'
    from_postprocessors_to_be_preserved = 'center_duct_heat_loss_00 center_duct_heat_loss_01 center_duct_heat_loss_02 center_duct_heat_loss_03 center_duct_heat_loss_04 center_duct_heat_loss_05 center_duct_heat_loss_06'
    to_postprocessors_to_be_preserved = 'Total_Net_Power_Through_Duct'
    allow_skipped_adjustment = true
  []

  [T_duct] # retrieve Tduct from SCM solve
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = subchannel
    source_variable = Tduct
    variable = duct_surface_temperature
    greedy_search = true
    to_boundaries = 'prsb_interface_00 prsb_interface_01  prsb_interface_02 prsb_interface_03
                     prsb_interface_04 prsb_interface_05 prsb_interface_06'
    execute_on = 'timestep_end'
  []
[]

(sfr/subchannel/multiple_SCM_assemblies/7assemblies/abr_7assemblies.i)

# a MOOSE mesh for 7 ABR assemblies
# sqrt(3) / 2 is by how much flat to flat is smaller than corner to corner
f = '${fparse sqrt(3) / 2}'

# units are cm - do not forget to convert to meter
outer_duct_out = 15.8123 # outer size of the hexagonal duct (side to side)
outer_duct_in = 15.0191 # flat to flat
inter_wrapper_width = 0.4348
height = 480.2

# discretization
n_ax = 50
ns = 10
duct_intervals_perishperic = '4 4'

[Mesh]
  [XX00]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '12'
    polygon_size = '${fparse outer_duct_out / 2 + inter_wrapper_width / 2}'
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall inter_wrapper'
    outward_interface_boundary_names = 'wall_in_00 wall_out_00'
    interface_boundary_id_shift = 100
  []

  [XX01]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '13'
    polygon_size = '${fparse outer_duct_out / 2 + inter_wrapper_width / 2}'
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_01 wall_out_01'
    interface_boundary_id_shift = 200
  []

  [XX02]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '14'
    polygon_size = '${fparse outer_duct_out / 2 + inter_wrapper_width / 2}'
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_02 wall_out_02'
    interface_boundary_id_shift = 300
  []

  [XX03]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '15'
    polygon_size = '${fparse outer_duct_out / 2 + inter_wrapper_width / 2}'
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_03 wall_out_03'
    interface_boundary_id_shift = 400
  []

  [XX04]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '16'
    polygon_size = '${fparse outer_duct_out / 2 + inter_wrapper_width / 2}'
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_04 wall_out_04'
    interface_boundary_id_shift = 500
  []

  [XX05]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '17'
    polygon_size = '${fparse outer_duct_out / 2 + inter_wrapper_width / 2}'
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_05 wall_out_05'
    interface_boundary_id_shift = 600
  []

  [XX06]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '18'
    polygon_size = '${fparse outer_duct_out / 2 + inter_wrapper_width / 2}'
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_06 wall_out_06'
    interface_boundary_id_shift = 700
  []

  [pattern]
    type = PatternedHexMeshGenerator
    inputs = 'XX00 XX01 XX02 XX03 XX04 XX05 XX06'
    pattern = '5 4 ;
              6 0 3 ;
               1 2'
    pattern_boundary = none
  []

  [extrude]
    type = AdvancedExtruderGenerator
    direction = '0 0 1'
    input = pattern
    heights = '${height}'
    num_layers = '${n_ax}'
  []

  [rename]
    type = RenameBlockGenerator
    input = extrude
    old_block = '12             13             14             15             16             17             18'
    new_block = 'porous_flow_00 porous_flow_01 porous_flow_02 porous_flow_03 porous_flow_04 porous_flow_05 porous_flow_06'
  []

  [inlet_interwrapper]
    type = ParsedGenerateSideset
    input = rename
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'inter_wrapper duct_wall'
    normal = '0 0 -1'
    new_sideset_name = inlet_interwrapper
  []

  [inlet_porous_flow_hfd]
    type = ParsedGenerateSideset
    input = inlet_interwrapper
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_01
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_01
  []

  [inlet_porous_flow_p]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_hfd
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_02
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_02
  []

  [inlet_porous_flow_d1]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_p
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_03
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_03
  []

  [inlet_porous_flow_d2]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_d1
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_04
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_04
  []

  [inlet_porous_flow_k011]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_d2
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_05
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_05
  []

  [inlet_porous_flow_x402]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_k011
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_06
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_06
  []

  [inlet_central_assembly]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_x402
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_00'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_00
  []

  [outlet_interwrapper]
    type = ParsedGenerateSideset
    input = inlet_central_assembly
    included_subdomains = 'inter_wrapper duct_wall'
    combinatorial_geometry = 'abs(z - ${fparse height}) < 1e-6'
    normal = '0 0 1'
    new_sideset_name = outlet_interwrapper
  []

  [outlet_porous_flow]
    type = ParsedGenerateSideset
    input = outlet_interwrapper
    included_subdomains = 'porous_flow_01 porous_flow_02 porous_flow_03 porous_flow_04 porous_flow_05 porous_flow_06'
    combinatorial_geometry = 'abs(z - ${fparse height}) < 1e-6'
    normal = '0 0 1'
    new_sideset_name = outlet_porous_flow
  []

  [outlet_central_assembly]
    type = ParsedGenerateSideset
    input = outlet_porous_flow
    included_subdomains = 'porous_flow_00'
    combinatorial_geometry = 'abs(z - ${fparse height}) < 1e-6'
    normal = '0 0 1'
    new_sideset_name = outlet_porous_flow_00
  []

  [rotate]
    type = TransformGenerator
    input = outlet_central_assembly
    transform = ROTATE
    vector_value = '0 0 0'
  []

  # turn into meters
  [scale]
    type = TransformGenerator
    vector_value = '0.01 0.01 0.01'
    transform = SCALE
    input = rotate
  []

  [new_wall_boundary_00]
    type = SideSetsBetweenSubdomainsGenerator
    input = scale
    new_boundary = 'prsb_interface_00'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_00'
  []

  [new_wall_boundary_01]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_00
    new_boundary = 'prsb_interface_01'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_01'
  []

  [new_wall_boundary_02]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_01
    new_boundary = 'prsb_interface_02'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_02'
  []

  [new_wall_boundary_03]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_02
    new_boundary = 'prsb_interface_03'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_03'
  []

  [new_wall_boundary_04]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_03
    new_boundary = 'prsb_interface_04'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_04'
  []

  [new_wall_boundary_05]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_04
    new_boundary = 'prsb_interface_05'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_05'
  []

  [new_wall_boundary_06]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_05
    new_boundary = 'prsb_interface_06'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_06'
  []

  [delete_assemblies]
    type = BlockDeletionGenerator
    block = 'porous_flow_00 porous_flow_01 porous_flow_02 porous_flow_03 porous_flow_04 porous_flow_05 porous_flow_06'
    input = 'new_wall_boundary_06'
  []
[]

(sfr/subchannel/multiple_SCM_assemblies/7assemblies/HC_master.i)

# This heat conduction model (main app) drives a multiapp setup where subchannel hexagonal ducts (solid wrappers)
# and an interwrapper (liquid sodium interwrapper between the ducts) are modeled and the temperature field is retrieved.
# This model treats the liquid sodium interwrapper as a solid and only heat conduction is considered.
# The calculated inner duct_wall axial heat flux [W/m] is transfered to 7 SCM models that run in parallel as multiapps.
# The method of coupling is domain decomposition:
# As mentioned above the heat conduction solution sends linear heat rate calculated on the inner duct surfaces to SCM
# and retrieves inner duct surface temperature from SCM. All other surfaces of the model have a Neuman BC.
##################### Instructions ###########################
# Create mesh file: abr_7assemblies_in.e
# run with : mpiexec -n N location-of-executable -i HC_master.i --keep-cout -w
# Units are SI

#### physics parameters
inlet_temperature = 653.15

####  Fluid properties of Sodium
####  Density #####
# A12 = 1.00423e3
# A13 = -0.21390
# A14 = -1.1046e-5
# rho = ${fparse A12 + A13 * inlet_temperature + A14 * inlet_temperature * inlet_temperature}
# #### Viscosity
# A52 = 3.6522e-5
# A53 = 0.16626
# A54 = -4.56877e1
# A55 = 2.8733e4
# mu = ${fparse A52 + A53 / inlet_temperature + A54 / inlet_temperature / inlet_temperature +
#         A55 / (inlet_temperature * inlet_temperature * inlet_temperature)}
#### Specific heat at constant pressure
# A28 = 7.3898e5
# A29 = 3.154e5
# A30 = 1.1340e3
# A31 = -2.2153e-1
# A32 = 1.1156e-4
# dt = ${fparse 2503.3 - inlet_temperature}
# cp = ${fparse A28 / dt / dt + A29 / dt + A30 + A31 * dt + A32 * dt * dt}
#### Heat conduction coefficient
A48 = 1.1045e2
A49 = -6.5112e-2
A50 = 1.5430e-5
A51 = -2.4617e-9
k_sodium = '${fparse A48 + A49 * inlet_temperature + A50 * inlet_temperature * inlet_temperature +
        A51 * inlet_temperature * inlet_temperature * inlet_temperature}'

# wrapper properties
k_wrapper = 15

wrapper_blocks = 'duct_wall'
inter_wrapper_blocks = 'inter_wrapper'

[Mesh]
  [fmesh]
    type = FileMeshGenerator
    file = 'abr_7assemblies_in.e'
  []
  coord_type = XYZ
[]

[Materials]
  [wall_heat_conductor]
    type = HeatConductionMaterial
    thermal_conductivity = ${k_wrapper}
    block = ${wrapper_blocks}
  []
  [sodium_heat_conductor]
    type = HeatConductionMaterial
    thermal_conductivity = ${k_sodium}
    block = ${inter_wrapper_blocks}
  []
[]

[Variables]
  [T_wrapper]
    order = FIRST
    family = LAGRANGE
  []
[]


[Kernels]
  [heat_conduction]
    type = HeatConduction
    variable = T_wrapper
    block = '${wrapper_blocks} ${inter_wrapper_blocks}'
  []
[]

[BCs]
  [T_wrapper_inside_wall]
    type = MatchedValueBC
    variable = T_wrapper
    boundary = 'prsb_interface_00 prsb_interface_01 prsb_interface_02 prsb_interface_03
                prsb_interface_04 prsb_interface_05 prsb_interface_06'
    v = duct_surface_temperature
  []
  [outside_bc]
   type = NeumannBC
   variable = T_wrapper
   boundary = '10000 inlet_interwrapper outlet_interwrapper'
  []
[]

[AuxVariables]
  [duct_surface_temperature]
    initial_condition = ${inlet_temperature}
  []

  [duct_heat_flux]
    order = CONSTANT
    family = MONOMIAL
    block = ${wrapper_blocks}
    initial_condition = 0
  []
[]

[AuxKernels]
  [HEAT_FLUX]
    type = DiffusionFluxAux
    diffusivity = ${k_wrapper}
    variable = duct_heat_flux
    diffusion_variable = T_wrapper
    component = normal
    boundary = 'prsb_interface_00 prsb_interface_01 prsb_interface_02 prsb_interface_03
                prsb_interface_04 prsb_interface_05 prsb_interface_06'
    execute_on = 'initial timestep_end'
  []
[]

[Executioner]
  type = Steady
  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  fixed_point_max_its = 5
  fixed_point_min_its = 4
  fixed_point_rel_tol = 1e-3
  fixed_point_abs_tol = 1e-3

  [Quadrature]
    order = THIRD
    side_order = FOURTH
  []
[]

[Postprocessors]
  # Total diffusive heat loss through central duct and the rest
  [center_duct_heat_loss_00]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_00
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_01]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_01
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_02]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_02
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_03]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_03
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_04]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_04
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_05]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_05
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_06]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_06
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
[]

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

################################################################################
######    A multiapp that couples heat conduction to subchannel via duct interface
################################################################################
[MultiApps]
  [subchannel]
    type = FullSolveMultiApp
    input_files = 'fuel_assembly_center.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i'
    execute_on = 'timestep_end'
    positions = '0         0           0
                 0.140704 -0.0812355   0
                 0.140704 0.0812355    0
                 0  0.162471           0
                 -0.140704 0.0812355   0
                 -0.140704 -0.0812355  0
                 0 -0.162471           0'
    max_procs_per_app = 1
    output_in_position = true
    bounding_box_padding = '0 0 0.1'
  []
[]

[Transfers]
  [q_duct] # send duct heat flux from heat conduction solve, to SCM
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = subchannel
    source_variable = duct_heat_flux
    variable = duct_heat_flux
    greedy_search = true
    from_boundaries = 'prsb_interface_00 prsb_interface_01  prsb_interface_02 prsb_interface_03
                       prsb_interface_04 prsb_interface_05 prsb_interface_06'
    execute_on = 'timestep_end'
    from_postprocessors_to_be_preserved = 'center_duct_heat_loss_00 center_duct_heat_loss_01 center_duct_heat_loss_02 center_duct_heat_loss_03 center_duct_heat_loss_04 center_duct_heat_loss_05 center_duct_heat_loss_06'
    to_postprocessors_to_be_preserved = 'Total_Net_Power_Through_Duct'
    allow_skipped_adjustment = true
  []

  [T_duct] # retrieve Tduct from SCM solve
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = subchannel
    source_variable = Tduct
    variable = duct_surface_temperature
    greedy_search = true
    to_boundaries = 'prsb_interface_00 prsb_interface_01  prsb_interface_02 prsb_interface_03
                     prsb_interface_04 prsb_interface_05 prsb_interface_06'
    execute_on = 'timestep_end'
  []
[]

(sfr/subchannel/multiple_SCM_assemblies/7assemblies/fuel_assembly_1.i)

###################################################
# Thermal-hydraulics parameters
###################################################
T_in = 653.15
P_out = 758423 # Pa
mdot_core = 5024
n_fuel_assemblies = 180
flow_area = 0.00650277
mass_flux_in = '${fparse mdot_core/n_fuel_assemblies/flow_area}' # kg/(m2.s)

###################################################
# Geometric parameters
###################################################
#f = ${fparse sqrt(3) / 2}

# units are cm - do not forget to convert to meter
scale_factor = 0.01
fuel_element_pitch = '${fparse 16.2471*scale_factor}'
inter_assembly_gap = '${fparse 0.4348*scale_factor}'
duct_thickness = '${fparse 0.3966*scale_factor}'
fuel_pin_pitch = '${fparse 0.8966*scale_factor}'
fuel_pin_diameter = '${fparse 0.7714*scale_factor}'
wire_z_spacing = '${fparse 20.43*scale_factor}' ###
wire_diameter = '${fparse 0.1307*scale_factor}' #check
n_rings = 10
length_entry_fuel = '${fparse 160.92*scale_factor}'
length_heated_fuel = '${fparse 85.82*scale_factor}'
length_outlet_fuel =  '${fparse 233.46*scale_factor}'
orifice_plate_height = '${fparse 5*scale_factor}'
duct_outside = '${fparse fuel_element_pitch - inter_assembly_gap}'
duct_inside = '${fparse duct_outside - 2 * duct_thickness}'
###################################################

[TriSubChannelMesh]
  [subchannel]
    type = SCMTriSubChannelMeshGenerator
    nrings = '${fparse n_rings}'
    n_cells = 50
    flat_to_flat = '${fparse duct_inside}'
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pin_diameter = '${fparse fuel_pin_diameter}'
    pitch = '${fparse fuel_pin_pitch}'
    dwire = '${fparse wire_diameter}'
    hwire = '${fparse wire_z_spacing}'
    spacer_z = '${fparse orifice_plate_height} ${fparse length_entry_fuel}'
    spacer_k = '0.5 0.5'
  []

  [fuel_pins]
    type = SCMTriPinMeshGenerator
    input = subchannel
    nrings = '${fparse n_rings}'
    n_cells = 50 #100
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pitch = '${fparse fuel_pin_pitch}'
  []

  [duct]
    type = SCMTriDuctMeshGenerator
    input = fuel_pins
    nrings = '${fparse n_rings}'
    n_cells = 50 #100
    flat_to_flat = '${fparse duct_inside}'
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pitch = '${fparse fuel_pin_pitch}'
  []
[]
# The function is evaluated with the origin z at unheated_length_entry
[Functions]
  [axial_heat_rate]
    type = ParsedFunction
    expression = '(pi/2)*sin(pi*z/L)'
    symbol_names = 'L'
    symbol_values = '${length_heated_fuel}'
  []
[]

[AuxVariables]
  [Tduct]
    block = duct
    initial_condition = 653.15
  []
[]

[FluidProperties]
  [sodium]
    type = PBSodiumFluidProperties
  []
[]

[SubChannel]
  type = TriSubChannel1PhaseProblem
  fp = sodium
  n_blocks = 1
  P_out = ${P_out}
  CT = 1.0
  compute_density = true
  compute_viscosity = true
  compute_power = true
  P_tol = 1.0e-4
  T_tol = 1.0e-3
  implicit = true
  segregated = false
  verbose_multiapps = true
  verbose_subchannel = true
[]

[ICs]
  [S_IC]
    type = SCMTriFlowAreaIC
    variable = S
  []

  [w_perim_IC]
    type = SCMTriWettedPerimIC
    variable = w_perim
  []

  [q_prime_IC]
    type = SCMTriPowerIC
    variable = q_prime
    power = 1000000.0 #W
    filename = "pin_power_profile_uni.txt"
    axial_heat_rate = axial_heat_rate
  []

  [T_ic]
    type = ConstantIC
    variable = T
    value = ${T_in}
  []

  [P_ic]
    type = ConstantIC
    variable = P
    value = 0.0 #always set to zero
  []

  [DP_ic]
    type = ConstantIC
    variable = DP
    value = 0.0
  []

  [Dpin_ic]
    type = ConstantIC
    variable = Dpin
    value = ${fuel_pin_diameter}
  []

  [Viscosity_ic]
    type = ViscosityIC
    variable = mu
    p = ${P_out}
    T = T
    fp = sodium
  []

  [rho_ic]
    type = RhoFromPressureTemperatureIC
    variable = rho
    p = ${P_out}
    T = T
    fp = sodium
  []

  [h_ic]
    type = SpecificEnthalpyFromPressureTemperatureIC
    variable = h
    p = ${P_out}
    T = T
    fp = sodium
  []

  [mdot_ic]
    type = ConstantIC
    variable = mdot
    value = 0.0
  []
[]

[AuxKernels]
  [P_out_bc]
    type = ConstantAux
    variable = P
    boundary = outlet
    value = ${P_out}
    execute_on = 'timestep_begin'
    block = subchannel
  []
  [T_in_bc]
    type = ConstantAux
    variable = T
    boundary = inlet
    value = ${T_in}
    execute_on = 'timestep_begin'
    block = subchannel
  []
  [mdot_in_bc]
    type = SCMMassFlowRateAux
    variable = mdot
    boundary = inlet
    area = S
    mass_flux = ${mass_flux_in}
    execute_on = 'timestep_begin'
    block = subchannel
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

!include SCM_output.i

[Executioner]
  type = Steady
[]

################################################################################
# A multiapp that projects data to a detailed mesh
################################################################################
[MultiApps]
  [viz]
    type = FullSolveMultiApp
    input_files = "3d.i"
    execute_on = "final"
    output_in_position = true
  []
[]

[Transfers]
  [subchannel_transfer]
      type = SCMSolutionTransfer
      to_multi_app = viz
      variable = 'mdot SumWij P DP h T rho mu S w_perim'
  []

  [pin_transfer]
      type = SCMPinSolutionTransfer
      to_multi_app = viz
      variable = 'Tpin Dpin q_prime'
  []
[]

(sfr/subchannel/multiple_SCM_assemblies/7assemblies/fuel_assembly_center.i)

###################################################
# Thermal-hydraulics parameters
###################################################
T_in = 653.15
P_out = 758423 # Pa
mdot_core = 5024
n_fuel_assemblies = 180
flow_area = 0.00650277
mass_flux_in = '${fparse mdot_core/n_fuel_assemblies/flow_area}' # kg/(m2.s)
###################################################
# Geometric parameters
###################################################
#f = ${fparse sqrt(3) / 2}

# units are cm - do not forget to convert to meter
scale_factor = 0.01
fuel_element_pitch = '${fparse 16.2471*scale_factor}'
inter_assembly_gap = '${fparse 0.4348*scale_factor}' # check
duct_thickness = '${fparse 0.3966*scale_factor}'
fuel_pin_pitch = '${fparse 0.8966*scale_factor}'
fuel_pin_diameter = '${fparse 0.7714*scale_factor}'
wire_z_spacing = '${fparse 20.43*scale_factor}' ###
wire_diameter = '${fparse 0.1307*scale_factor}' #check
n_rings = 10
length_entry_fuel = '${fparse 160.92*scale_factor}'
length_heated_fuel = '${fparse 85.82*scale_factor}'
length_outlet_fuel = '${fparse 233.46*scale_factor}'
orifice_plate_height = '${fparse 5*scale_factor}'
duct_outside = '${fparse fuel_element_pitch - inter_assembly_gap}'
duct_inside = '${fparse duct_outside - 2 * duct_thickness}'
###################################################

[TriSubChannelMesh]
  [subchannel]
    type = SCMTriSubChannelMeshGenerator
    nrings = '${fparse n_rings}'
    n_cells = 50
    flat_to_flat = '${fparse duct_inside}'
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pin_diameter = '${fparse fuel_pin_diameter}'
    pitch = '${fparse fuel_pin_pitch}'
    dwire = '${fparse wire_diameter}'
    hwire = '${fparse wire_z_spacing}'
    spacer_z = '${fparse orifice_plate_height} ${fparse length_entry_fuel}'
    spacer_k = '0.5 0.5'
  []

  [fuel_pins]
    type = SCMTriPinMeshGenerator
    input = subchannel
    nrings = '${fparse n_rings}'
    n_cells = 50
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pitch = '${fparse fuel_pin_pitch}'
  []

  [duct]
    type = SCMTriDuctMeshGenerator
    input = fuel_pins
    nrings = '${fparse n_rings}'
    n_cells = 50
    flat_to_flat = '${fparse duct_inside}'
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pitch = '${fparse fuel_pin_pitch}'
  []
[]
# The function is evaluated with the origin z at unheated_length_entry
[Functions]
  [axial_heat_rate]
    type = ParsedFunction
    expression = '(pi/2)*sin(pi*z/L)'
    symbol_names = 'L'
    symbol_values = '${length_heated_fuel}'
  []
[]

[AuxVariables]
  [Tduct]
    block = duct
    initial_condition = 653.15
  []
[]

[FluidProperties]
  [sodium]
    type = PBSodiumFluidProperties
  []
[]

[SubChannel]
  type = TriSubChannel1PhaseProblem
  fp = sodium
  n_blocks = 1
  P_out = ${P_out}
  CT = 1.0
  compute_density = true
  compute_viscosity = true
  compute_power = true
  P_tol = 1.0e-4
  T_tol = 1.0e-3
  implicit = true
  segregated = false
  verbose_multiapps = true
  verbose_subchannel = true
[]

[ICs]
  [S_IC]
    type = SCMTriFlowAreaIC
    variable = S
  []

  [w_perim_IC]
    type = SCMTriWettedPerimIC
    variable = w_perim
  []

  [q_prime_IC]
    type = SCMTriPowerIC
    variable = q_prime
    power = 2000000.0 #W
    filename = "pin_power_profile_uni.txt"
    axial_heat_rate = axial_heat_rate
  []

  [T_ic]
    type = ConstantIC
    variable = T
    value = ${T_in}
  []

  [P_ic]
    type = ConstantIC
    variable = P
    value = 0.0 #always set to zero
  []

  [DP_ic]
    type = ConstantIC
    variable = DP
    value = 0.0
  []

  [Dpin_ic]
    type = ConstantIC
    variable = Dpin
    value = ${fuel_pin_diameter}
  []

  [Viscosity_ic]
    type = ViscosityIC
    variable = mu
    p = ${P_out}
    T = T
    fp = sodium
  []

  [rho_ic]
    type = RhoFromPressureTemperatureIC
    variable = rho
    p = ${P_out}
    T = T
    fp = sodium
  []

  [h_ic]
    type = SpecificEnthalpyFromPressureTemperatureIC
    variable = h
    p = ${P_out}
    T = T
    fp = sodium
  []

  [mdot_ic]
    type = ConstantIC
    variable = mdot
    value = 0.0
  []
[]

[AuxKernels]
  [P_out_bc]
    type = ConstantAux
    variable = P
    boundary = outlet
    value = ${P_out}
    execute_on = 'timestep_begin'
    block = subchannel
  []
  [T_in_bc]
    type = ConstantAux
    variable = T
    boundary = inlet
    value = ${T_in}
    execute_on = 'timestep_begin'
    block = subchannel
  []
  [mdot_in_bc]
    type = SCMMassFlowRateAux
    variable = mdot
    boundary = inlet
    area = S
    mass_flux = ${mass_flux_in}
    execute_on = 'timestep_begin'
    block = subchannel
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

!include SCM_output.i

[Executioner]
  type = Steady
[]

################################################################################
# A multiapp that projects data to a detailed mesh
################################################################################
[MultiApps]
  [viz]
    type = FullSolveMultiApp
    input_files = "3d.i"
    execute_on = "final"
    output_in_position = true
  []
[]

[Transfers]
  [subchannel_transfer]
    type = SCMSolutionTransfer
    to_multi_app = viz
    variable = 'mdot SumWij P DP h T rho mu S w_perim'
  []

  [pin_transfer]
    type = SCMPinSolutionTransfer
    to_multi_app = viz
    variable = 'Tpin Dpin q_prime'
  []
[]

(sfr/subchannel/multiple_SCM_assemblies/7assemblies/3d.i)

###################################################
# Geometric parameters
###################################################
#f = ${fparse sqrt(3) / 2}

# units are cm - do not forget to convert to meter
scale_factor = 0.01
fuel_element_pitch = '${fparse 16.2471*scale_factor}'
inter_assembly_gap = '${fparse 0.4348*scale_factor}'
duct_thickness = '${fparse 0.3966*scale_factor}'
fuel_pin_pitch = '${fparse 0.8966*scale_factor}'
fuel_pin_diameter = '${fparse 0.7714*scale_factor}'
length_entry_fuel = '${fparse 160.92*scale_factor}'
length_heated_fuel = '${fparse 85.82*scale_factor}'
length_outlet_fuel = '${fparse 233.46*scale_factor}'
duct_outside = ${fparse fuel_element_pitch - inter_assembly_gap}
duct_inside = ${fparse duct_outside -  2 * duct_thickness}
n_rings = 10

[Mesh]
  [subchannel]
    type = SCMDetailedTriSubChannelMeshGenerator
    nrings = '${fparse n_rings}'
    n_cells = 50 #100
    flat_to_flat = '${fparse duct_inside}'
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pin_diameter = '${fparse fuel_pin_diameter}'
    pitch = '${fparse fuel_pin_pitch}'
  []

  [fuel_pins]
    type = SCMDetailedTriPinMeshGenerator
    input = subchannel
    nrings = '${fparse n_rings}'
    n_cells = 50
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pitch = '${fparse fuel_pin_pitch}'
    pin_diameter = '${fparse fuel_pin_diameter}'
  []
[]

[AuxVariables]
  [mdot]
    block = subchannel
  []
  [SumWij]
    block = subchannel
  []
  [P]
    block = subchannel
  []
  [DP]
    block = subchannel
  []
  [h]
    block = subchannel
  []
  [T]
    block = subchannel
  []
  [rho]
    block = subchannel
  []
  [mu]
    block = subchannel
  []
  [S]
    block = subchannel
  []
  [w_perim]
    block = subchannel
  []
  [q_prime]
    block = fuel_pins
  []
  [Tpin]
    block = fuel_pins
  []
  [Dpin]
    block = fuel_pins
  []
[]

[Problem]
  type = NoSolveProblem
[]

[Outputs]
  exodus = true
[]

[Executioner]
  type = Steady
[]

(sfr/subchannel/multiple_SCM_assemblies/19assemblies/abr_19assemblies.i)

# a MOOSE mesh for 19 ABR assemblies
# sqrt(3) / 2 is by how much flat to flat is smaller than corner to corner
f = ${fparse sqrt(3) / 2}

# units are cm - do not forget to convert to meter
outer_duct_out = 15.8123 # outer size of the hexagonal duct (side to side)
outer_duct_in = 15.0191 # flat to flat
inter_wrapper_width = 0.4348
height = 480.2

# discretization
n_ax = 50
ns = 10
duct_intervals_perishperic = '4 4'

[Mesh]
  [XX00]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '12'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall inter_wrapper'
    outward_interface_boundary_names = 'wall_in_00 wall_out_00'
    interface_boundary_id_shift = 100
  []

  [XX01]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '13'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_01 wall_out_01'
    interface_boundary_id_shift = 200
  []

  [XX02]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '14'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_02 wall_out_02'
    interface_boundary_id_shift = 300
  []

  [XX03]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '15'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_03 wall_out_03'
    interface_boundary_id_shift = 400
  []

  [XX04]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '16'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_04 wall_out_04'
    interface_boundary_id_shift = 500
  []

  [XX05]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '17'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_05 wall_out_05'
    interface_boundary_id_shift = 600
  []

  [XX06]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '18'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_06 wall_out_06'
    interface_boundary_id_shift = 700
  []

  [XX07]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '19'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_07 wall_out_07'
    interface_boundary_id_shift = 800
  []

  [XX08]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '20'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_08 wall_out_08'
    interface_boundary_id_shift = 900
  []

  [XX09]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '21'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_09 wall_out_09'
    interface_boundary_id_shift = 1000
  []

  [XX10]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '22'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_10 wall_out_10'
    interface_boundary_id_shift = 1100
  []

  [XX11]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '23'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_11 wall_out_11'
    interface_boundary_id_shift = 1200
  []

  [XX12]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '24'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_12 wall_out_12'
    interface_boundary_id_shift = 1300
  []

  [XX13]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '25'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_13 wall_out_13'
    interface_boundary_id_shift = 1400
  []

  [XX14]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '26'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_14 wall_out_14'
    interface_boundary_id_shift = 1500
  []

  [XX15]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '27'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_15 wall_out_15'
    interface_boundary_id_shift = 1600
  []

  [XX16]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '28'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_16 wall_out_16'
    interface_boundary_id_shift = 1700
  []

  [XX17]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '29'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_17 wall_out_17'
    interface_boundary_id_shift = 1800
  []

  [XX18]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6
    num_sectors_per_side = '${ns} ${ns} ${ns} ${ns} ${ns} ${ns}'
    background_intervals = 1
    background_block_ids = '30'
    polygon_size = ${fparse outer_duct_out / 2 + inter_wrapper_width / 2}
    duct_sizes = '${fparse outer_duct_in / f /2} ${fparse outer_duct_out / f / 2}'
    duct_intervals = ${duct_intervals_perishperic}
    duct_block_names = 'duct_wall  inter_wrapper'
    outward_interface_boundary_names = 'wall_in_18 wall_out_18'
    interface_boundary_id_shift = 1900
  []

  [pattern]
    type = PatternedHexMeshGenerator
    inputs = 'XX00 XX01 XX02 XX03 XX04 XX05 XX06 XX07 XX08 XX09
              XX10 XX11 XX12 XX13 XX14 XX15 XX16 XX17 XX18'
    pattern =
            ' 15 14 13;
             16 5 4 12;
            17 6 0 3 11;
             18 1 2 10;
               7 8 9'
    pattern_boundary = none
  []

  [extrude]
    type = AdvancedExtruderGenerator
    direction = '0 0 1'
    input = pattern
    heights = '${height}'
    num_layers = '${n_ax}'
  []

  [rename]
    type = RenameBlockGenerator
    input = extrude
    old_block = '12             13              14            15             16             17             18
                                19              20            21             22             23             24
                                25              26            27             28             29             30'
    new_block = 'porous_flow_00 porous_flow_01 porous_flow_02 porous_flow_03 porous_flow_04 porous_flow_05 porous_flow_06
                                porous_flow_07 porous_flow_08 porous_flow_09 porous_flow_10 porous_flow_11 porous_flow_12
                                porous_flow_13 porous_flow_14 porous_flow_15 porous_flow_16 porous_flow_17 porous_flow_18'
  []

  [inlet_interwrapper]
    type = ParsedGenerateSideset
    input = rename
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'inter_wrapper duct_wall'
    normal = '0 0 -1'
    new_sideset_name = inlet_interwrapper
  []

  [inlet_porous_flow_hfd]
    type = ParsedGenerateSideset
    input = inlet_interwrapper
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_01
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_01
  []

  [inlet_porous_flow_p]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_hfd
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_02
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_02
  []

  [inlet_porous_flow_d1]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_p
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_03
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_03
  []

  [inlet_porous_flow_d2]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_d1
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_04
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_04
  []

  [inlet_porous_flow_k011]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_d2
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_05
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_05
  []

  [inlet_porous_flow_x402]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_k011
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = porous_flow_06
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_06
  []

  [inlet_central_assembly]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_x402
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_00'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_00
  []

  [inlet_porous_flow_07]
    type = ParsedGenerateSideset
    input = inlet_central_assembly
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_07'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_07
  []

  [inlet_porous_flow_08]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_07
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_08'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_08
  []

  [inlet_porous_flow_09]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_08
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_09'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_09
  []

  [inlet_porous_flow_10]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_09
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_10'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_10
  []

  [inlet_porous_flow_11]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_10
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_11'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_11
  []

  [inlet_porous_flow_12]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_11
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_12'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_12
  []

  [inlet_porous_flow_13]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_12
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_13'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_13
  []

  [inlet_porous_flow_14]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_13
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_14'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_14
  []

  [inlet_porous_flow_15]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_14
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_15'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_15
  []

  [inlet_porous_flow_16]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_15
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_16'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_16
  []

  [inlet_porous_flow_17]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_16
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_17'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_17
  []

  [inlet_porous_flow_18]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_17
    combinatorial_geometry = 'abs(z) < 1e-6'
    included_subdomains = 'porous_flow_18'
    normal = '0 0 -1'
    new_sideset_name = inlet_porous_flow_18
  []

  [outlet_interwrapper]
    type = ParsedGenerateSideset
    input = inlet_porous_flow_18 # outlet_interwall
    included_subdomains = 'inter_wrapper duct_wall'
    combinatorial_geometry = 'abs(z - ${fparse height}) < 1e-6'
    normal = '0 0 1'
    new_sideset_name = outlet_interwrapper
  []

  [outlet_porous_flow]
    type = ParsedGenerateSideset
    input = outlet_interwrapper
    included_subdomains = 'porous_flow_01 porous_flow_02 porous_flow_03 porous_flow_04 porous_flow_05 porous_flow_06 porous_flow_07 porous_flow_08 porous_flow_09 porous_flow_10 porous_flow_11 porous_flow_12 porous_flow_13 porous_flow_14 porous_flow_15 porous_flow_16 porous_flow_17 porous_flow_18'
    combinatorial_geometry = 'abs(z - ${fparse height}) < 1e-6'
    normal = '0 0 1'
    new_sideset_name = outlet_porous_flow
  []

  [outlet_central_assembly]
    type = ParsedGenerateSideset
    input = outlet_porous_flow
    included_subdomains = 'porous_flow_00'
    combinatorial_geometry = 'abs(z - ${fparse height}) < 1e-6'
    normal = '0 0 1'
    new_sideset_name = outlet_porous_flow_00
  []

  [rotate]
    type = TransformGenerator
    input = outlet_central_assembly
    transform = ROTATE
    vector_value = '0 0 0'
  []

  # turn into meters
  [scale]
    type = TransformGenerator
    vector_value = '0.01 0.01 0.01'
    transform = SCALE
    input = rotate
  []

  [new_wall_boundary_00]
    type = SideSetsBetweenSubdomainsGenerator
    input = scale
    new_boundary = 'prsb_interface_00'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_00'
  []

  [new_wall_boundary_01]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_00
    new_boundary = 'prsb_interface_01'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_01'
  []

  [new_wall_boundary_02]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_01
    new_boundary = 'prsb_interface_02'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_02'
  []

  [new_wall_boundary_03]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_02
    new_boundary = 'prsb_interface_03'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_03'
  []

  [new_wall_boundary_04]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_03
    new_boundary = 'prsb_interface_04'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_04'
  []

  [new_wall_boundary_05]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_04
    new_boundary = 'prsb_interface_05'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_05'
  []

  [new_wall_boundary_06]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_05
    new_boundary = 'prsb_interface_06'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_06'
  []

  [new_wall_boundary_07]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_06
    new_boundary = 'prsb_interface_07'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_07'
  []

  [new_wall_boundary_08]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_07
    new_boundary = 'prsb_interface_08'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_08'
  []

  [new_wall_boundary_09]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_08
    new_boundary = 'prsb_interface_09'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_09'
  []

  [new_wall_boundary_10]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_09
    new_boundary = 'prsb_interface_10'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_10'
  []

  [new_wall_boundary_11]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_10
    new_boundary = 'prsb_interface_11'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_11'
  []

  [new_wall_boundary_12]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_11
    new_boundary = 'prsb_interface_12'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_12'
  []

  [new_wall_boundary_13]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_12
    new_boundary = 'prsb_interface_13'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_13'
  []

  [new_wall_boundary_14]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_13
    new_boundary = 'prsb_interface_14'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_14'
  []

  [new_wall_boundary_15]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_14
    new_boundary = 'prsb_interface_15'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_15'
  []

  [new_wall_boundary_16]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_15
    new_boundary = 'prsb_interface_16'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_16'
  []

  [new_wall_boundary_17]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_16
    new_boundary = 'prsb_interface_17'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_17'
  []

  [new_wall_boundary_18]
    type = SideSetsBetweenSubdomainsGenerator
    input = new_wall_boundary_17
    new_boundary = 'prsb_interface_18'
    primary_block = 'duct_wall'
    paired_block = 'porous_flow_18'
  []

  [delete_assemblies]
    type = BlockDeletionGenerator
    block = 'porous_flow_00 porous_flow_01 porous_flow_02 porous_flow_03 porous_flow_04 porous_flow_05 porous_flow_06 porous_flow_07 porous_flow_08 porous_flow_09
             porous_flow_10 porous_flow_11 porous_flow_12 porous_flow_13 porous_flow_14 porous_flow_15 porous_flow_16 porous_flow_17 porous_flow_18'
    input = 'new_wall_boundary_18'
  []
[]

(sfr/subchannel/multiple_SCM_assemblies/19assemblies/HC_master.i)

# This heat conduction model (main app) drives a multiapp setup where subchannel hexagonal ducts (solid wrappers)
# and an interwrapper (liquid sodium interwrapper between the ducts) are modeled and the temperature field is retrieved.
# This model treats the liquid sodium interwrapper as a solid and only heat conduction is considered.
# The calculated inner duct_wall axial heat flux [W/m] is transfered to 19 SCM models that run in parallel as multiapps.
# The method of coupling is domain decomposition:
# As mentioned above the heat conduction solution sends linear heat rate calculated on the inner duct surfaces to SCM
# and retrieves inner duct surface temperature from SCM. All other surfaces of the model has a Neuman BC.
##################### Instructions ###########################
# Create mesh file: abr_19assemblies_in.e
# run with : mpiexec -n N location-of-executable -i HC_master.i --keep-cout -w
# Units are SI

# physics parameters
inlet_temperature = 653.15

## fluid properties
#####  Density #####
# A12 = 1.00423e3
# A13 = -0.21390
# A14 = -1.1046e-5
# rho = '${fparse A12 + A13 * inlet_temperature + A14 * inlet_temperature * inlet_temperature}'
# #### Viscosity
# A52 = 3.6522e-5
# A53 = 0.16626
# A54 = -4.56877e1
# A55 = 2.8733e4
# mu = '${fparse A52 + A53 / inlet_temperature + A54 / inlet_temperature / inlet_temperature +
#         A55 / (inlet_temperature * inlet_temperature * inlet_temperature)} '
# #### Specific heat at constant pressure
# A28 = 7.3898e5
# A29 = 3.154e5
# A30 = 1.1340e3
# A31 = -2.2153e-1
# A32 = 1.1156e-4
# dt = '${fparse 2503.3 - inlet_temperature}'
# cp = '${fparse A28 / dt / dt + A29 / dt + A30 + A31 * dt + A32 * dt * dt}'
#### Heat conduction coefficient
A48 = 1.1045e2
A49 = -6.5112e-2
A50 = 1.5430e-5
A51 = -2.4617e-9
k_sodium = '${fparse A48 + A49 * inlet_temperature + A50 * inlet_temperature * inlet_temperature +
        A51 * inlet_temperature * inlet_temperature * inlet_temperature} '

# wrapper properties
k_wrapper = 15

wrapper_blocks = 'duct_wall'
inter_wrapper_blocks = 'inter_wrapper'

[Mesh]
  [fmesh]
    type = FileMeshGenerator
    file = 'abr_19assemblies_in.e'
  []
  coord_type = XYZ
[]

[Materials]
  [wall_heat_conductor]
    type = HeatConductionMaterial
    thermal_conductivity = ${k_wrapper}
    block = ${wrapper_blocks}
  []
  [sodium_heat_conductor]
    type = HeatConductionMaterial
    thermal_conductivity = ${k_sodium}
    block = ${inter_wrapper_blocks}
  []
[]

[Variables]
  [T_wrapper]
    order = FIRST
    family = LAGRANGE
  []
[]

[Kernels]
  [heat_conduction]
    type = HeatConduction
    variable = T_wrapper
    block = '${wrapper_blocks} ${inter_wrapper_blocks}'
  []
[]

[BCs]
  [T_wrapper_inside_wall]
    type = MatchedValueBC
    variable = T_wrapper
    boundary = 'prsb_interface_00 prsb_interface_01 prsb_interface_02 prsb_interface_03
                prsb_interface_04 prsb_interface_05 prsb_interface_06 prsb_interface_07
                prsb_interface_08 prsb_interface_09 prsb_interface_10 prsb_interface_11
                prsb_interface_12 prsb_interface_13 prsb_interface_14 prsb_interface_15
                prsb_interface_16 prsb_interface_17 prsb_interface_18'
    v = duct_surface_temperature
  []
  [outside_bc]
    type = NeumannBC
    variable = T_wrapper
    boundary = '10000 inlet_interwrapper outlet_interwrapper'
  []
[]

[AuxVariables]
  [duct_surface_temperature]
    initial_condition = ${inlet_temperature}
  []

  [duct_heat_flux]
    order = CONSTANT
    family = MONOMIAL
    block = ${wrapper_blocks}
    initial_condition = 0
  []
[]

[AuxKernels]
  [HEAT_FLUX]
    type = DiffusionFluxAux
    diffusivity = ${k_wrapper}
    variable = duct_heat_flux
    diffusion_variable = T_wrapper
    component = normal
    boundary = 'prsb_interface_00 prsb_interface_01 prsb_interface_02 prsb_interface_03
                prsb_interface_04 prsb_interface_05 prsb_interface_06 prsb_interface_07
                prsb_interface_08 prsb_interface_09 prsb_interface_10 prsb_interface_11
                prsb_interface_12 prsb_interface_13 prsb_interface_14 prsb_interface_15
                prsb_interface_16 prsb_interface_17 prsb_interface_18'
    execute_on = 'initial timestep_end'
  []
[]

[Executioner]
  type = Steady
  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  fixed_point_max_its = 5
  fixed_point_min_its = 4
  fixed_point_rel_tol = 1e-3
  fixed_point_abs_tol = 1e-3

  [Quadrature]
    order = THIRD
    side_order = FOURTH
  []
[]

[Postprocessors]
  # Total diffusive heat loss through central duct and the rest
  [center_duct_heat_loss_00]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_00
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_01]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_01
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_02]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_02
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_03]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_03
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_04]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_04
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_05]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_05
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_06]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_06
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
    [center_duct_heat_loss_07]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_07
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_08]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_08
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_09]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_09
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_10]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_10
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_11]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_11
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_12]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_12
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_13]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_13
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
    [center_duct_heat_loss_14]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_14
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_15]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_15
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_16]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_16
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_17]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_17
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
  [center_duct_heat_loss_18]
    type        = ADSideDiffusiveFluxIntegral
    variable    = T_wrapper
    boundary    = prsb_interface_18
    diffusivity = ${k_wrapper}
    execute_on='timestep_end'
  []
[]

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

################################################################################
######    A multiapp that couples heat conduction to subchannel via duct interface
################################################################################

[MultiApps]
  [subchannel]
    type = FullSolveMultiApp
    input_files = 'fuel_assembly_center.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i
    fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i
    fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i
    fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i
    fuel_assembly_1.i fuel_assembly_1.i fuel_assembly_1.i'
    execute_on = 'timestep_end'
    positions = '0          0            0
                 0.140704  -0.0812355    0
                 0.140704   0.0812355    0
                 0          0.162471     0
                -0.140704   0.0812355    0
                -0.140704  -0.0812355    0
                 0         -0.162471     0
                 0.281408  -0.162471     0
                 0.281408   0            0
                 0.281408   0.162471     0
                 0.140704   0.2437065    0
                 0          0.324942     0
                -0.140704   0.2437065    0
                -0.281408   0.162471     0
                -0.281408   0            0
                -0.281408  -0.162471     0
                -0.140704  -0.2437065    0
                 0         -0.324942     0
                 0.140704  -0.2437065    0'
    max_procs_per_app = 1
    output_in_position = true
    bounding_box_padding = '0 0 0.1'
  []
[]

[Transfers]
  [q_duct] # send duct heat flux from heat conduction solve, to SCM
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = subchannel
    source_variable = duct_heat_flux
    variable = duct_heat_flux
    greedy_search = true
    from_boundaries = 'prsb_interface_00 prsb_interface_01 prsb_interface_02 prsb_interface_03
                prsb_interface_04 prsb_interface_05 prsb_interface_06 prsb_interface_07
                prsb_interface_08 prsb_interface_09 prsb_interface_10 prsb_interface_11
                prsb_interface_12 prsb_interface_13 prsb_interface_14 prsb_interface_15
                prsb_interface_16 prsb_interface_17 prsb_interface_18'
    execute_on = 'timestep_end'
    from_postprocessors_to_be_preserved = 'center_duct_heat_loss_00 center_duct_heat_loss_01 center_duct_heat_loss_02 center_duct_heat_loss_03 center_duct_heat_loss_04 center_duct_heat_loss_05 center_duct_heat_loss_06
    center_duct_heat_loss_07 center_duct_heat_loss_08 center_duct_heat_loss_09 center_duct_heat_loss_10 center_duct_heat_loss_11
    center_duct_heat_loss_12 center_duct_heat_loss_13 center_duct_heat_loss_14 center_duct_heat_loss_15 center_duct_heat_loss_16
    center_duct_heat_loss_17 center_duct_heat_loss_18'
    to_postprocessors_to_be_preserved = 'Total_Net_Power_Through_Duct'
    allow_skipped_adjustment = true
  []

  [T_duct] # retrieve Tduct from SCM solve
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = subchannel
    source_variable = Tduct
    variable = duct_surface_temperature
    greedy_search = true
    to_boundaries = 'prsb_interface_00 prsb_interface_01 prsb_interface_02 prsb_interface_03
                prsb_interface_04 prsb_interface_05 prsb_interface_06 prsb_interface_07
                prsb_interface_08 prsb_interface_09 prsb_interface_10 prsb_interface_11
                prsb_interface_12 prsb_interface_13 prsb_interface_14 prsb_interface_15
                prsb_interface_16 prsb_interface_17 prsb_interface_18'
    execute_on = 'timestep_end'
  []
[]

(sfr/subchannel/multiple_SCM_assemblies/19assemblies/fuel_assembly_center.i)

###################################################
# Thermal-hydraulics parameters
###################################################
T_in = 653.15
P_out = 758423 # Pa
mdot_core = 5024
n_fuel_assemblies = 180
flow_area = 0.00650277
mass_flux_in = '${fparse mdot_core/n_fuel_assemblies/flow_area}' # kg/(m2.s)

###################################################
# Geometric parameters
###################################################
#f = ${fparse sqrt(3) / 2}

# units are cm - do not forget to convert to meter
scale_factor = 0.01
fuel_element_pitch = '${fparse 16.2471*scale_factor}'
inter_assembly_gap = '${fparse 0.4348*scale_factor}' # check
duct_thickness = '${fparse 0.3966*scale_factor}'
fuel_pin_pitch = '${fparse 0.8966*scale_factor}'
fuel_pin_diameter = '${fparse 0.7714*scale_factor}'
wire_z_spacing = '${fparse 20.43*scale_factor}' ###
wire_diameter = '${fparse 0.1307*scale_factor}' #check
n_rings = 10
length_entry_fuel = '${fparse 160.92*scale_factor}'
length_heated_fuel = '${fparse 85.82*scale_factor}'
length_outlet_fuel = '${fparse 233.46*scale_factor}'
orifice_plate_height = '${fparse 5*scale_factor}'
duct_outside = '${fparse fuel_element_pitch - inter_assembly_gap}'
duct_inside = '${fparse duct_outside - 2 * duct_thickness}'
###################################################

[TriSubChannelMesh]
  [subchannel]
    type = SCMTriSubChannelMeshGenerator
    nrings = '${fparse n_rings}'
    n_cells = 50
    flat_to_flat = '${fparse duct_inside}'
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pin_diameter = '${fparse fuel_pin_diameter}'
    pitch = '${fparse fuel_pin_pitch}'
    dwire = '${fparse wire_diameter}'
    hwire = '${fparse wire_z_spacing}'
    spacer_z = '${fparse orifice_plate_height} ${fparse length_entry_fuel}'
    spacer_k = '0.5 0.5'
  []

  [fuel_pins]
    type = SCMTriPinMeshGenerator
    input = subchannel
    nrings = '${fparse n_rings}'
    n_cells = 50 #100
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pitch = '${fparse fuel_pin_pitch}'
  []

  [duct]
    type = SCMTriDuctMeshGenerator
    input = fuel_pins
    nrings = '${fparse n_rings}'
    n_cells = 50 #100
    flat_to_flat = '${fparse duct_inside}'
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pitch = '${fparse fuel_pin_pitch}'
  []
[]
# The function is evaluated with the origin z at unheated_length_entry
[Functions]
  [axial_heat_rate]
    type = ParsedFunction
    expression = '(pi/2)*sin(pi*z/L)'
    symbol_names = 'L'
    symbol_values = '${length_heated_fuel}'
  []
[]

[AuxVariables]
  [Tduct]
    block = duct
    initial_condition = 653.15
  []
[]

[FluidProperties]
  [sodium]
    type = PBSodiumFluidProperties
  []
[]

[SubChannel]
  type = TriSubChannel1PhaseProblem
  fp = sodium
  n_blocks = 1
  P_out = ${P_out}
  CT = 1.0
  compute_density = true
  compute_viscosity = true
  compute_power = true
  P_tol = 1.0e-4
  T_tol = 1.0e-3
  implicit = true
  segregated = false
  verbose_multiapps = true
  verbose_subchannel = true
[]

[ICs]
  [S_IC]
    type = SCMTriFlowAreaIC
    variable = S
  []

  [w_perim_IC]
    type = SCMTriWettedPerimIC
    variable = w_perim
  []

  [q_prime_IC]
    type = SCMTriPowerIC
    variable = q_prime
    power = 2000000.0 #W
    filename = "pin_power_profile_uni.txt"
    axial_heat_rate = axial_heat_rate
  []

  [T_ic]
    type = ConstantIC
    variable = T
    value = ${T_in}
  []

  [P_ic]
    type = ConstantIC
    variable = P
    value = 0.0 #always set to zero
  []

  [DP_ic]
    type = ConstantIC
    variable = DP
    value = 0.0
  []

  [Dpin_ic]
    type = ConstantIC
    variable = Dpin
    value = ${fuel_pin_diameter}
  []

  [Viscosity_ic]
    type = ViscosityIC
    variable = mu
    p = ${P_out}
    T = T
    fp = sodium
  []

  [rho_ic]
    type = RhoFromPressureTemperatureIC
    variable = rho
    p = ${P_out}
    T = T
    fp = sodium
  []

  [h_ic]
    type = SpecificEnthalpyFromPressureTemperatureIC
    variable = h
    p = ${P_out}
    T = T
    fp = sodium
  []

  [mdot_ic]
    type = ConstantIC
    variable = mdot
    value = 0.0
  []
[]

[AuxKernels]
  [P_out_bc]
    type = ConstantAux
    variable = P
    boundary = outlet
    value = ${P_out}
    execute_on = 'timestep_begin'
    block = subchannel
  []
  [T_in_bc]
    type = ConstantAux
    variable = T
    boundary = inlet
    value = ${T_in}
    execute_on = 'timestep_begin'
    block = subchannel
  []
  [mdot_in_bc]
    type = SCMMassFlowRateAux
    variable = mdot
    boundary = inlet
    area = S
    mass_flux = ${mass_flux_in}
    execute_on = 'timestep_begin'
    block = subchannel
  []
[]

[Outputs]
  csv = true
[]

!include SCM_output.i

[Executioner]
  type = Steady
[]

################################################################################
# A multiapp that projects data to a detailed mesh
################################################################################
[MultiApps]
  [viz]
    type = FullSolveMultiApp
    input_files = "3d.i"
    execute_on = "final"
    output_in_position = true
  []
[]

[Transfers]
  [subchannel_transfer]
    type = SCMSolutionTransfer
    to_multi_app = viz
    variable = 'mdot SumWij P DP h T rho mu S w_perim'
  []

  [pin_transfer]
    type = SCMPinSolutionTransfer
    to_multi_app = viz
    variable = 'Tpin Dpin q_prime'
  []
[]

(sfr/subchannel/multiple_SCM_assemblies/19assemblies/3d.i)

###################################################
# Geometric parameters
###################################################
#f = ${fparse sqrt(3) / 2}

# units are cm - do not forget to convert to meter
scale_factor = 0.01
fuel_element_pitch = '${fparse 16.2471*scale_factor}'
inter_assembly_gap = '${fparse 0.4348*scale_factor}'
duct_thickness = '${fparse 0.3966*scale_factor}'
fuel_pin_pitch = '${fparse 0.8966*scale_factor}'
fuel_pin_diameter = '${fparse 0.7714*scale_factor}'
length_entry_fuel = '${fparse 160.92*scale_factor}'
length_heated_fuel = '${fparse 85.82*scale_factor}'
length_outlet_fuel = '${fparse 233.46*scale_factor}'
duct_outside = ${fparse fuel_element_pitch - inter_assembly_gap}
duct_inside = ${fparse duct_outside -  2 * duct_thickness}
n_rings = 10

[Mesh]
  [subchannel]
    type = SCMDetailedTriSubChannelMeshGenerator
    nrings = '${fparse n_rings}'
    n_cells = 50 #100
    flat_to_flat = '${fparse duct_inside}'
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pin_diameter = '${fparse fuel_pin_diameter}'
    pitch = '${fparse fuel_pin_pitch}'
  []

  [fuel_pins]
    type = SCMDetailedTriPinMeshGenerator
    input = subchannel
    nrings = '${fparse n_rings}'
    n_cells = 50
    unheated_length_entry = '${fparse length_entry_fuel}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse length_outlet_fuel}'
    pitch = '${fparse fuel_pin_pitch}'
    pin_diameter = '${fparse fuel_pin_diameter}'
  []
[]

[AuxVariables]
  [mdot]
    block = subchannel
  []
  [SumWij]
    block = subchannel
  []
  [P]
    block = subchannel
  []
  [DP]
    block = subchannel
  []
  [h]
    block = subchannel
  []
  [T]
    block = subchannel
  []
  [rho]
    block = subchannel
  []
  [mu]
    block = subchannel
  []
  [S]
    block = subchannel
  []
  [w_perim]
    block = subchannel
  []
  [q_prime]
    block = fuel_pins
  []
  [Tpin]
    block = fuel_pins
  []
  [Dpin]
    block = fuel_pins
  []
[]

[Problem]
  type = NoSolveProblem
[]

[Outputs]
  exodus = true
[]

[Executioner]
  type = Steady
[]

Example Description
Coupling scheme
Energy conservation
Input files
Results