MSFR Griffin-Pronghorn Model

Contact: Mauricio Tano, mauricio.tanoretamales.at.inl.gov

Model link: Griffin-Pronghorn Steady-State Model

The MultiApp system is used to separate the neutronics and the fluid dynamics problems. The fluid system is solved by the subapp and it uses the run_ns.i input files. (Here "ns" is an abbreviation for Navier-Stokes.)

The fluid system includes conservation equations for fluid mass, momentum, and energy as well as the conservation of delayed neutron precursors.

Conservation of fluid mass

The conservation of mass is,

$var element = document.getElementById("moose-equation-6357506e-0ea8-49b2-8870-9df3fc57ecb6");katex.render(" \\nabla \\cdot \\rho \\vec{u} = 0,", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-f7309fee-a422-4b3d-9780-a9451f6df1d5");katex.render("\\rho", element, {displayMode:false,throwOnError:false});$ is the fluid density and $var element = document.getElementById("moose-equation-3dd4b1f4-980b-4eb3-92f3-f9d652309ad3");katex.render("\\vec{u}", element, {displayMode:false,throwOnError:false});$ is the velocity vector.

Here the system will be simplified by modeling the flow as incompressible. (The effect of Buoyancy will be re-introduced later with the Boussinesq approximation.) The simplified conservation of mass is then given by,

$var element = document.getElementById("moose-equation-141a7ef7-b9a0-42a9-b6b5-5aa60549ed5e");katex.render(" \\nabla \\cdot \\vec{u} = 0.", element, {displayMode:true,throwOnError:false});$

This conservation equation is automatically added by the NavierStokesFV action when selecting the incompressible model.


[Modules]
  [NavierStokesFV]
    # General parameters
    compressibility = 'incompressible'
    ...

Legacy syntax for the mass equation is included here.

Conservation of fluid momentum

This system also includes the conservation of momentum in the $var element = document.getElementById("moose-equation-90edd814-85b1-4f80-9165-0dba8999800b");katex.render("x", element, {displayMode:false,throwOnError:false});$ -direction. A fairly general form of the steady-state condition is,

$var element = document.getElementById("moose-equation-1325a403-976f-4867-9b00-dad8c4dbb187");katex.render(" \\nabla \\cdot \\rho \\vec{u} u = -\\frac{\\partial}{\\partial x} P + \\mu \\nabla^2 u + f_{\\text{fric},x} + \\rho \\vec{g} \\cdot \\hat{x},", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-6cce0b36-573a-4181-beb0-f1db8bb1aaf9");katex.render("u", element, {displayMode:false,throwOnError:false});$ is the $var element = document.getElementById("moose-equation-541e6446-6d36-4644-beeb-f9fa2cf7dab4");katex.render("x", element, {displayMode:false,throwOnError:false});$ component of the velocity, $var element = document.getElementById("moose-equation-04689b1f-02c8-47f1-8abd-a248af4cf696");katex.render("P", element, {displayMode:false,throwOnError:false});$ is the pressure, $var element = document.getElementById("moose-equation-75e53d8b-664b-473d-a4e7-af7d54e7fcf0");katex.render("f_{\\text{fric},x}", element, {displayMode:false,throwOnError:false});$ is the $var element = document.getElementById("moose-equation-36f223db-af42-40a9-acc6-e38377987d77");katex.render("x", element, {displayMode:false,throwOnError:false});$ component of the viscous friction force, and $var element = document.getElementById("moose-equation-e37e9441-a147-4390-bae1-139f6be79b83");katex.render("\\vec{g}", element, {displayMode:false,throwOnError:false});$ is the gravity vector.

In this model, gravity will point in the negative $var element = document.getElementById("moose-equation-7daeea46-014a-41ce-b732-7345d2b4b383");katex.render("y", element, {displayMode:false,throwOnError:false});$ -direction so the quantity $var element = document.getElementById("moose-equation-ad8bccef-bed3-4f8d-ace3-492721bd5fc4");katex.render("\\vec{g} \\cdot \\hat{x}", element, {displayMode:false,throwOnError:false});$ is zero,

$var element = document.getElementById("moose-equation-41c21a16-f859-4117-9fa4-8686808c9def");katex.render(" \\nabla \\cdot \\rho \\vec{u} u = -\\frac{\\partial}{\\partial x} P + \\mu \\nabla^2 u + f_{\\text{fric},x}.", element, {displayMode:true,throwOnError:false});$

Practical simulations require modifications to the momentum equations in order to model the effects of turbulence without explicitly resolving the turbulent structures. Here, we will apply the Reynolds-averaging procedure and the Boussinesq hypothesis so that the effect of turbulent momentum transfer is modeled with a term analogous to viscous shear,

$var element = document.getElementById("moose-equation-fe15014d-b088-4467-b368-0b1bac379bdd");katex.render(" \\nabla \\cdot \\rho \\vec{u} u = -\\frac{\\partial}{\\partial x} P + \\left( \\mu + \\mu_t \\right) \\nabla^2 u + f_{\\text{fric},x},", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-85ed55d9-d536-4a6f-928f-4ba51d0d805d");katex.render("\\mu_t", element, {displayMode:false,throwOnError:false});$ is the turbulent viscosity.

Here, a zero-equation model based on the mixing length model is used. In this model the turbulent viscosity is defined as: $var element = document.getElementById("moose-equation-e232c4ce-705d-44db-911d-f1d304c56693");katex.render(" \\mu_t = \\rho \\cdot {l_m}^2 \\cdot |2\\overline{\\overline S} : \\overline{\\overline S}|", element, {displayMode:true,throwOnError:false});$ and $var element = document.getElementById("moose-equation-7d1a1427-9e4f-459c-9c10-2356470ce961");katex.render(" \\overline{\\overline S} = 0.5 \\cdot \\left( \\nabla u + \\nabla u^t \\right).", element, {displayMode:true,throwOnError:false});$

The standard Prandtl's mixing length model dictates that $var element = document.getElementById("moose-equation-53586e2f-6948-4e94-81b3-45b7203d4491");katex.render("l_m", element, {displayMode:false,throwOnError:false});$ has a linear dependence on the distance to the nearest wall. However, for this simulation we implement a capped mixing length model (Escudier, 1966) that defines the mixing length as

$var element = document.getElementById("moose-equation-510d36e7-da86-4fc5-9afe-bc6f2db1a0b3");katex.render(" l_m = \\kappa y_d, \\quad \\text{if} \\: \\kappa y_d < \\kappa_0 \\delta, \\\\ l_m = \\kappa_0 \\delta, \\quad \\text{if} \\: \\kappa y_d \\geq \\kappa_0 \\delta,", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-ec66cb19-5451-45b6-8df3-6dd3ea02ee8f");katex.render("\\kappa =0.41", element, {displayMode:false,throwOnError:false});$ is the Von Karman constant, $var element = document.getElementById("moose-equation-9c1bde1d-d5b2-4249-a77d-c96e6d58740a");katex.render("\\kappa_0 = 0.09", element, {displayMode:false,throwOnError:false});$ as in Escudier's model and $var element = document.getElementById("moose-equation-3ab94f17-8b6f-4bc2-83ba-494bb130b1dc");katex.render("\\delta", element, {displayMode:false,throwOnError:false});$ has length units and represents the thickness of the velocity boundary layer.

The viscous friction is treated with two different models for different regions of the reactor. The bulk of the friction is expected to occur in the heat exchanger region where there is a large surface area in contact with the salt to facilitate heat transfer. Rather than explicitly modeling the heat exchanger geometry, a simple, tunable model will be used, $var element = document.getElementById("moose-equation-faef7d3c-fe6f-4183-aa68-b9f08ecd956d");katex.render(" f_{\\text{fric},x} = -C_q u | u |,", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-e3049221-6736-4115-a9eb-12deecd19aaf");katex.render("C_q", element, {displayMode:false,throwOnError:false});$ is a tunable volumetric friction coefficient, which is selected to match the desired pressure drop, in conjunction with the pump head value.

A typical viscous friction model could be used for the regions outside of the heat exchanger. However, the viscous effects will be negligible due to the large, uniform eddy viscosity model used here. Consequently, the viscous force will be neglected outside of the heat exchanger, $var element = document.getElementById("moose-equation-f9bbdbcb-f62d-4476-b707-c2506aa21c06");katex.render(" f_{\\text{fric},x} = 0.", element, {displayMode:true,throwOnError:false});$

By convention, we must collect all of the terms on one side of the equation. This gives the form that is implemented for the MSFR model,

$(1)var element = document.getElementById("moose-equation-1bdc211c-0907-45ff-a98a-a6b9d11f52a9");katex.render(" \\nabla \\cdot \\rho \\vec{u} u + \\frac{\\partial}{\\partial x} P - \\left( \\mu + \\mu_t \\right) \\nabla^2 u - f_{\\text{fric},x} = 0. ", element, {displayMode:true,throwOnError:false});$

The first few terms in Eq. (1), including the advection of momentum, the pressure gradient and the laminar diffusion are automatically handled by the NavierStokesFV action.

The mixing length turbulence model is added by specifying additional parameters to the action. This basic model is currently the most limiting issue in terms of accuracy.


[Modules]
  [NavierStokesFV]
    ...
    # Turbulence parameters
    turbulence_handling = 'mixing-length'
    turbulent_prandtl = ${Pr_t}
    von_karman_const = ${von_karman_const}
    mixing_length_delta = 0.1
    mixing_length_walls = 'shield_wall reflector_wall'
    mixing_length_aux_execute_on = 'initial'
    ...

The friction in the heat exchanger is specified by passing the following parameters to the NavierStokesFV action. A quadratic friction model with a coefficient tuned to obtain the desired mass flow rate is selected.


[Modules]
  [NavierStokesFV]
    ...
    # Heat exchanger
    friction_blocks = 'hx'
    friction_types = 'FORCHHEIMER'
    friction_coeffs = ${friction}
    ...

This model does not include a friction force for the other blocks so no kernel needs to be specified for those blocks.

The conservation of momentum in the $var element = document.getElementById("moose-equation-aaefe751-82b2-44ce-94a2-fdd48889eb0a");katex.render("y", element, {displayMode:false,throwOnError:false});$ -direction is analogous, but it also includes the Boussinesq approximation in order to capture the effect of buoyancy and a body force in the pump region. Note that the Boussinesq term is needed because of the approximation that the fluid density is uniform and constant,

$var element = document.getElementById("moose-equation-656b4430-d52b-47a2-9544-80548c420d68");katex.render(" \\nabla \\cdot \\rho \\vec{u} v + \\frac{\\partial}{\\partial y} P - \\left( \\mu + \\mu_t \\right) \\nabla^2 v - f_{\\text{fric}, y} - f_{\\text{pump}} - \\rho \\alpha \\vec{g} \\left( T - T_0 \\right) = 0,", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-abeab497-cd85-490a-819a-f97cc43f2b3c");katex.render("f_{\\text{pump}}", element, {displayMode:false,throwOnError:false});$ is the pump head driving the flow, $var element = document.getElementById("moose-equation-1e8c521b-5036-4618-a2aa-35bbac6f2e5c");katex.render("\\alpha", element, {displayMode:false,throwOnError:false});$ is the expansion coefficient, $var element = document.getElementById("moose-equation-beae63b3-f26e-47b5-892a-686188f901b4");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the fluid temperature, and $var element = document.getElementById("moose-equation-17379e57-28c8-4486-80b6-78b030205cd8");katex.render("T_0", element, {displayMode:false,throwOnError:false});$ is a reference temperature value. The pump head is tuned such that the imposed mass flow rate is ~18500 kg/s.

For each kernel describing the $var element = document.getElementById("moose-equation-21466ab2-ada6-4730-9aaf-9a3a02b2b0b0");katex.render("x", element, {displayMode:false,throwOnError:false});$ -momentum equation, there is a corresponding kernel for the $var element = document.getElementById("moose-equation-117eca6a-a974-4dad-aa3c-c7a95a992b49");katex.render("y", element, {displayMode:false,throwOnError:false});$ -momentum equation. They are similarly added by the NavierStokesFV. An additional kernel is added to add a volumetric pump component to the system:

[FVKernels]
  [pump]
    type = INSFVBodyForce
    variable = vel_y
    functor = ${pump_force}
    block = 'pump'
    momentum_component = 'y'
  []
[]

The Boussinesq buoyancy approximation is added to the model using this parameter:


[Modules]
  [NavierStokesFV]
    ...
    boussinesq_approximation = true
    ...

Boundary conditions include standard velocity wall functions at the walls to account for the non-linearity of the velocity in the boundary layer given the coarse mesh and symmetry at the center axis of the MSFR. They are added by the action based on user-input parameters:


[Modules]
  [NavierStokesFV]
    ...
    # Boundary conditions
    wall_boundaries = 'shield_wall reflector_wall fluid_symmetry'
    momentum_wall_types = 'wallfunction wallfunction symmetry'
    ...

Auxkernels are used to compute the wall shear stress obtained by the standard wall function model, the dimensionless wall distance $var element = document.getElementById("moose-equation-10548181-b6b2-4749-a8d4-f17e4fd728cb");katex.render("y^+", element, {displayMode:false,throwOnError:false});$ and the value for the eddy viscosity. These are used for analysis purposes.

[AuxKernels]
  [mixing_len]
    type = WallDistanceMixingLengthAux
    walls = 'shield_wall reflector_wall'
    variable = mixing_len
    execute_on = 'initial'
    block = 'fuel'
    von_karman_const = ${von_karman_const}
    delta = 0.1 # capping parameter (m)
  []
  [wall_shear_stress]
    type = WallFunctionWallShearStressAux
    variable = wall_shear_stress
    walls = 'shield_wall reflector_wall'
    block = 'fuel'
    mu = 'total_viscosity'
  []
  [wall_yplus]
    type = WallFunctionYPlusAux
    variable = wall_yplus
    walls = 'shield_wall reflector_wall'
    block = 'fuel'
    mu = 'total_viscosity'
  []
  [turbulent_viscosity]
    type = INSFVMixingLengthTurbulentViscosityAux
    variable = eddy_viscosity
    block = 'fuel pump hx'
    execute_on = 'initial'
  []
[]

Legacy syntax for the momentum equation is included here.

Conservation of fluid energy

The steady-state conservation of energy can be expressed as, $var element = document.getElementById("moose-equation-2618dccc-be97-4749-be13-e1dde9e54eac");katex.render(" \\nabla \\cdot \\rho h \\vec{u} - \\nabla \\cdot \\lambda \\nabla T = Q_q,", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-0b77552a-4cdc-4c4f-8426-66b4a1fcb2b1");katex.render("h", element, {displayMode:false,throwOnError:false});$ is the fluid specific enthalpy, $var element = document.getElementById("moose-equation-e5703118-063e-4935-b51d-9bfd01d19112");katex.render("\\lambda", element, {displayMode:false,throwOnError:false});$ is the thermal conductivity, and $var element = document.getElementById("moose-equation-7a85f2d8-d762-493a-9810-a5278fb40954");katex.render("Q_q", element, {displayMode:false,throwOnError:false});$ is the volumetric heat generation rate.

Here it is expected that the energy released from nuclear reactions will be very large compared to pressure work terms. Consequently, we will use the simplified form, $var element = document.getElementById("moose-equation-912aa50c-9602-4eeb-8a2d-b369735f096e");katex.render(" \\nabla \\cdot \\rho c_p T \\vec{u} - \\nabla \\cdot \\lambda \\nabla T = Q_q.", element, {displayMode:true,throwOnError:false});$

As with the momentum equations, a practical solver requires a model for the turbulent transport of energy, $var element = document.getElementById("moose-equation-408d8a69-e87b-42c9-b0cf-951b1583e6aa");katex.render(" \\nabla \\cdot \\rho c_p T \\vec{u} - \\nabla \\cdot \\lambda \\nabla T - \\nabla \\cdot \\rho c_p \\epsilon_q \\nabla T = Q_q", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-8222a305-2e55-4047-b12d-07cc777adb2f");katex.render("\\epsilon_q", element, {displayMode:false,throwOnError:false});$ is the eddy diffusivity for heat. It is related to the eddy diffusivity for momentum (i.e. the kinematic eddy viscosity) as $var element = document.getElementById("moose-equation-698d3b25-cb1c-4b8d-98fb-444247a37857");katex.render("\\text{Pr}_t = \\nu_t / \\epsilon_q", element, {displayMode:false,throwOnError:false});$ where $var element = document.getElementById("moose-equation-90de2b9f-e292-47c6-9c29-14f0b4d7a16d");katex.render("\\text{Pr}_t", element, {displayMode:false,throwOnError:false});$ is the turbulent Prandtl number.

As with the molecular viscosity in the momentum equations, the simple turbulence model used here will overwhelm the thermal conductivity. Neglecting that term and moving the heat generation to the left-hand-side of the equation gives, $var element = document.getElementById("moose-equation-3e76500c-8de8-4076-851b-a212d0278852");katex.render(" \\nabla \\cdot \\rho c_p T \\vec{u} - \\nabla \\cdot \\rho c_p \\epsilon_q \\nabla T - Q_q = 0.", element, {displayMode:true,throwOnError:false});$

The heat generation due to nuclear reactions is computed by the neutronics solver, and this distribution will be used directly in the $var element = document.getElementById("moose-equation-e7b6bc8b-14f7-4ad2-9c9a-1fa1016688c7");katex.render("Q_q", element, {displayMode:false,throwOnError:false});$ term. The effect of the heat exchanger will also be included in the $var element = document.getElementById("moose-equation-87d89a23-4254-4e52-86a0-5e54e3f50365");katex.render("Q_q", element, {displayMode:false,throwOnError:false});$ term as a volumetric heat loss per the model, $var element = document.getElementById("moose-equation-6457668d-2cd6-463c-a311-e5d66e1e0b8c");katex.render(" Q_q = -\\alpha (T - T_\\text{ambient}),", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-fa508db0-ba21-4ce7-9686-dd0b163ed3e9");katex.render("\\alpha", element, {displayMode:false,throwOnError:false});$ is a coefficient (equal to the surface area density times the heat transfer coefficient) and $var element = document.getElementById("moose-equation-1501dd4d-0337-4a68-8b6e-7f84fb226167");katex.render("T_\\text{ambient}", element, {displayMode:false,throwOnError:false});$ is the temperature of the coolant on the secondary side of the heat exchanger.

Note that all material properties, including $var element = document.getElementById("moose-equation-5fdfa226-76f9-4488-948d-602af564eb0d");katex.render("\\rho", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-b34b7a62-27a9-4e7e-828f-07010e9befa8");katex.render("c_p", element, {displayMode:false,throwOnError:false});$ , are assumed constant. It will therefore be convenient to move the $var element = document.getElementById("moose-equation-327dd504-70d4-48ea-b372-b9e3703544d3");katex.render("\\rho c_p", element, {displayMode:false,throwOnError:false});$ factors outside of the divergence operators and divide the entire equation by that factor, $(2)var element = document.getElementById("moose-equation-b13997a6-d0c0-4782-b381-969838e4d149");katex.render(" \\nabla \\cdot T \\vec{u} - \\nabla \\cdot \\epsilon_q \\nabla T - \\frac{1}{\\rho c_p} Q_q = 0. ", element, {displayMode:true,throwOnError:false});$

Eq. (2) is the final equation that is implemented in this model. It is added to the equation system using a parameter in the NavierStokesFV action:


[Modules]
  [NavierStokesFV]
    ...
    add_energy_equation = true
    ...

The first term—energy advection— is automatically added by the action syntax. The second term—turbulent diffusion of heat— is added based on the momentum turbulent term, with the Prandtl number also to be specified to adjust the diffusion coefficient.

The third term—heat source and loss—is covered by two kernels. First, the nuclear heating computed by the neutronics solver is included with,


[Modules]
  [NavierStokesFV]
    ...
    # Heat source
    external_heat_source = power_density
    ...

And second, the heat loss through the heat exchanger is implemented with the kernel,


[Modules]
  [NavierStokesFV]
    ...
    # Heat exchanger
    ambient_convection_blocks = 'hx'
    ambient_convection_alpha = ${fparse 600 * 20e3} # HX specifications
    ambient_temperature = ${T_HX}
    ...

The boundary conditions are added similarly, by simply specifying 0 heat flux for symmetry and adiabatic boundaries.


[Modules]
  [NavierStokesFV]
    ...
    energy_wall_types = 'heatflux heatflux heatflux'
    energy_wall_function = '0 0 0'

Legacy syntax for the energy equation is included here.

Conservation of delayed neutron precursors

The steady-state conservation of delayed neutron precursors (DNP) can be expressed as, $(3)var element = document.getElementById("moose-equation-fac2458c-a13e-4dec-bf07-f5788c5e5620");katex.render(" \\nabla \\cdot (\\vec{u}c_i) - \\nabla \\cdot \\epsilon_c \\nabla c_i - \\lambda_i c_i = \\beta_i f,\\quad i=1,\\cdots,I, ", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-2956aa0a-dda5-4b2b-8f45-bfec86676157");katex.render("c_i", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-03b50016-7b5a-4c32-a3fc-b93e14f5d09a");katex.render("\\lambda_i", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-ea051820-0697-471c-8a8f-1b3afadf0743");katex.render("\\beta_i", element, {displayMode:false,throwOnError:false});$ are the concentration, the decay constant and the production fraction per neutron generated from fission of the $var element = document.getElementById("moose-equation-43c6c493-fc52-46a4-a26b-7d7bda015f36");katex.render("i", element, {displayMode:false,throwOnError:false});$ -th group of DNP respectively. $var element = document.getElementById("moose-equation-929b7ba2-7794-4554-8d0b-2643e368afc0");katex.render("I", element, {displayMode:false,throwOnError:false});$ is the total number of groups of DNP. $var element = document.getElementById("moose-equation-7c55b273-3129-4149-9182-bcfdb0becf51");katex.render("f", element, {displayMode:false,throwOnError:false});$ is the fission neutron production rate. $var element = document.getElementById("moose-equation-a8bd26ac-58b1-4323-950c-207c8d4c5406");katex.render("\\epsilon_c = \\frac{\\nu_t}{Sc_t}", element, {displayMode:false,throwOnError:false});$ is the eddy diffusivity for DNP, where $var element = document.getElementById("moose-equation-ced0fe93-beb0-4602-9feb-a62162dc753f");katex.render("\\nu_t", element, {displayMode:false,throwOnError:false});$ is the turbulnt viscosity and $var element = document.getElementById("moose-equation-b30ee941-5510-4ba8-ba3d-67a06b1d0c1f");katex.render("Sc_t", element, {displayMode:false,throwOnError:false});$ is the turbulent Schmidt number. It is noted the DNP term in the neutron transport equation is scaled with $var element = document.getElementById("moose-equation-9ebabe54-255a-46b7-81b2-0399c38646a4");katex.render("k", element, {displayMode:false,throwOnError:false});$ -effective as the prompt fission term for steady-state eigenvalue calculations, where $var element = document.getElementById("moose-equation-b897c131-3c85-4ca7-9fb0-7ed9c6a86673");katex.render("k", element, {displayMode:false,throwOnError:false});$ -effective (the eigenvalue) is part of the solution.

It is added to the equation system using a parameter in the NavierStokesFV action:


[Modules]
  [NavierStokesFV]
    ...
    add_scalar_equation = true
    ...

The terms of Eq. (2) are added with


[Modules]
  [NavierStokesFV]
    ...
    # Precursor advection, diffusion and source term
    passive_scalar_names = 'c1 c2 c3 c4 c5 c6'
    passive_scalar_schmidt_number = '${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t}'
    passive_scalar_coupled_source = 'fission_source c1; fission_source c2; fission_source c3;
                                     fission_source c4; fission_source c5; fission_source c6;'
    passive_scalar_coupled_source_coeff = '${beta1} ${lambda1_m}; ${beta2} ${lambda2_m};
                                           ${beta3} ${lambda3_m}; ${beta4} ${lambda4_m};
                                           ${beta5} ${lambda5_m}; ${beta6} ${lambda6_m}'

Neutronics

With Griffin, the process of converting the basic conservation equations into MOOSE variables and kernels is automated with the TransportSystems block:

[TransportSystems]
  particle = neutron
  equation_type = eigenvalue

  G = 6

  VacuumBoundary = 'outer'

  [diff]
    scheme = CFEM-Diffusion
    n_delay_groups = 6
    external_dnp_variable = 'dnp'
    family = LAGRANGE
    order = FIRST
    fission_source_aux = true

    # For PJFNKMO
    assemble_scattering_jacobian = true
    assemble_fission_jacobian = true
  []
[]

Details about neutron transport equations can be found in Griffin theory manual.

Here we are specifying an eigenvalue neutronics problem using 6 energy groups (G = 6) solved via the diffusion approximation with a continuous finite element discretization scheme (scheme = CFEM-Diffusion).

Note the external_dnp_variable = 'dnp' parameter. This is a special option needed for liquid-fueled MSRs which signals that the conservation equations for DNP will be handled "externally" from the default Griffin implementation which assumes that the precursors do not have the turbulent treatment. This parameter is referencing the dnp array auxiliary variable defined by the MSRNeutronicsFluidCoupling class, which in this case is an array auxiliary variable with 6 components, corresponding to the 6 DNP groups used here. Support within the Framework for array variables is somewhat limited. For example, not all of the multiapp transfers work with array variables, and the Navier-Stokes module does not include the kernels that are needed to advect an array variable. For this reason, there is also a separate auxiliary variable for each DNP group, which are also created within MSRNeutronicsFluidCoupling.

The MSRNeutronicsFluidCoupling block,

[MSRNeutronicsFluidCoupling]
  fluid_blocks = 'fuel pump hx'
  solid_blocks = 'reflector shield'
  n_dnp = 6
  use_transient_multi_app = false
  fluid_app_name = ns
  fluid_input_file = run_ns.i
  initial_fluid_temperature = 873.15
  fluid_temperature_name_neutronics_app = tfuel
  fluid_temperature_name_solid_suffix = _constant
  fluid_temperature_name_fluid_app = T_fluid
  dnp_name_prefix = c
  dnp_array_name = dnp
[]

creates a variety of objects needed on the neutronics side of the neutronics-fluid coupling for MSRs:

aux variables for the aforementioned DNP concentration variables and additionally, fluid temperature,
the aux kernel to form the DNP array from its components,
initial conditions for the fluid temperature if not restarting the simulation,
a MultiApp (TransientMultiApp for transients, FullSolveMultiApp otherwise) to run the fluid application, and
transfers for the DNP concentrations, power density, fission source, and fluid temperature.

The run_ns.i subapp is responsible for computing the precursor distributions, and the distributions are transferred from the subapp to the main app by MSRNeutronicsFluidCoupling and then internally copied from the individual DNP group variables into the DNP array variable using an aux kernel.

Also note that solving the neutronics problem requires a set of multigroup cross sections. Generating cross sections is a topic that is left outside the scope of this example. A set has been generated for the MSFR problem and stored in the repository using Griffin's ISOXML format. These cross sections are included by the blocks,

[Materials]
  [fuel]
    type = CoupledFeedbackNeutronicsMaterial
    library_name = 'msfr_xs'
    library_file = '../mgxs/msfr_xs.xml'
    grid_names = 'tfuel'
    grid_variables = 'tfuel'
    plus = true
    isotopes = 'pseudo'
    densities = '1.0'
    is_meter = true
    material_id = 2
    block = 'fuel pump hx'
  []
  [shield]
    type = CoupledFeedbackNeutronicsMaterial
    library_name = 'msfr_xs'
    library_file = '../mgxs/msfr_xs.xml'
    grid_names = 'tfuel'
    grid_variables = 'tfuel_constant'
    isotopes = 'pseudo'
    densities = '1.0'
    is_meter = true
    material_id = 5
    block = 'shield'
  []
  [reflector]
    type = CoupledFeedbackNeutronicsMaterial
    library_name = 'msfr_xs'
    library_file = '../mgxs/msfr_xs.xml'
    grid_names = 'tfuel'
    grid_variables = 'tfuel_constant'
    isotopes = 'pseudo'
    densities = '1.0'
    is_meter = true
    material_id = 1
    block = 'reflector'
  []
[]

CoupledFeedbackNeutronicsMaterial is able to use the temperature transferred from the fluid system for evaluating multigroup cross sections based on a table lookup on element quadrature points to bring in the feedback effect. It also has the capability of adjusting cross sections based on fluid density. In this model, the fluid density change is negligible.

The neutronics solution is normalized to the rated power 3000 MW with the PowerDensity input block:

[PowerDensity]
  power = 3e9
  power_density_variable = power_density
  power_scaling_postprocessor = power_scaling
  family = MONOMIAL
  order = CONSTANT
[]

The power density is evaluated with the normalized neutronics solution. It provides the source of the fluid energy equation. Because the fluid energy equation is discretized with FV, we evaluate the power density variable with constant monomial.

The calculation is driven by Eigenvalue, an executioner available in the MOOSE framework. The PJFNKMO option for the solve_type parameter is able to drive the eigenvalue calculation with the contribution of DNP to the neutron transport equation as an external source scaled with $var element = document.getElementById("moose-equation-cc00c523-3262-4766-9fe6-83b0387461c4");katex.render("k", element, {displayMode:false,throwOnError:false});$ -effective.

[Executioner]
  type = Eigenvalue
  solve_type = PJFNKMO

  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart '
  petsc_options_value = 'hypre boomeramg 50'
  l_max_its = 50

  free_power_iterations = 4 # important to obtain fundamental mode eigenvalue

  nl_abs_tol = 1e-9
  fixed_point_abs_tol = 1e-9
  fixed_point_max_its = 20
[]

References

M.P Escudier. The Distribution of Mixing-Length in Turbulent Flows Near Walls. PhD thesis, Imperial College, 1966.

@phdthesis{escudier1966,
    author = "Escudier, M.P",
    publisher = "Imperial College",
    school = "Imperial College",
    title = "{The Distribution of Mixing-Length in Turbulent Flows Near Walls}",
    year = "1966"
}

(msr/msfr/steady/run_ns.i)

################################################################################
## Molten Salt Fast Reactor - Euratom EVOL + Rosatom MARS Design              ##
## Pronghorn Sub-Application input file                                       ##
## Relaxation towards Steady state 3D thermal hydraulics model                ##
################################################################################
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html

# The flow in this simulation should be initialized with a previous flow
# solution (isothermal, heated or multiphysics) OR using a viscosity rampdown
# An isothermal solution can be generated using 'run_ns_initial.i', and is
# saved in 'restart/run_ns_initial_restart.e'
# A heated solution, with a flat power distribution, can be generated with
# this script 'run_ns.i', and is saved in 'restart/run_ns_restart.e'
# A coupled neutronics-coarse TH solution can be generated with
# 'run_neutronics.i', saved in 'restart/run_neutronics_ns_restart.e'

# Material properties
rho = 4284 # density [kg / m^3]  (@1000K)
cp = 1594 # specific heat capacity [J / kg / K]
drho_dT = 0.882 # derivative of density w.r.t. temperature [kg / m^3 / K]
mu = 0.0166 # viscosity [Pa s]
k = 1.7 # thermal conductivity [W / m / K]
# https://www.researchgate.net/publication/337161399_Development_of_a_control-\
# oriented_power_plant_simulator_for_the_molten_salt_fast_reactor/fulltext/5dc9\
# 5c9da6fdcc57503eec39/Development-of-a-control-oriented-power-plant-simulator-\
# for-the-molten-salt-fast-reactor.pdf
von_karman_const = 0.41

# Turbulent properties
Pr_t = 0.9 # turbulent Prandtl number
Sc_t = 1 # turbulent Schmidt number

# Derived material properties
alpha = '${fparse drho_dT / rho}' # thermal expansion coefficient

# Operating parameters
T_HX = 873.15 # heat exchanger temperature [K]

# Mass flow rate tuning, for heat exchanger pressure and temperature drop
friction = 4e3 # [kg / m^4]
pump_force = -20000. # [N / m^3]

# Delayed neutron precursor parameters. Lambda values are decay constants in
# [1 / s]. Beta values are production fractions.
lambda1_m = -0.0133104
lambda2_m = -0.0305427
lambda3_m = -0.115179
lambda4_m = -0.301152
lambda5_m = -0.879376
lambda6_m = -2.91303
beta1 = 8.42817e-05
beta2 = 0.000684616
beta3 = 0.000479796
beta4 = 0.00103883
beta5 = 0.000549185
beta6 = 0.000184087

[GlobalParams]
  rhie_chow_user_object = 'ins_rhie_chow_interpolator'
[]

################################################################################
# GEOMETRY
################################################################################

[Mesh]
  coord_type = 'RZ'
  rz_coord_axis = Y
  allow_renumbering = false

  [restart]
    type = FileMeshGenerator
    use_for_exodus_restart = true
    # Depending on the file chosen, the initialization of variables should be
    # adjusted. The following variables can be initalized:
    # - vel_x, vel_y, p from isothermal simulation
    # file = 'restart/run_ns_initial_restart.e'
    # Below are initialization points created from this input file
    # The variable IC should be set from_file_var for temperature and precursors
    # - vel_x, vel_y, p, T_fluid, c_i from cosine heated simulation
    # file = 'restart/run_ns_restart.e'
    # - vel_x, vel_y, p, T_fluid, c_i from coupled multiphysics simulation
    file = 'restart/run_neutronics_out_ns0_restart.e'
  []
  [hx_top]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'y > 0'
    included_subdomains = '3'
    included_neighbors = '1'
    fixed_normal = true
    normal = '0 1 0'
    new_sideset_name = 'hx_top'
    input = 'restart'
  []
  [hx_bot]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'y <-0.6'
    included_subdomains = '3'
    included_neighbors = '1'
    fixed_normal = true
    normal = '0 -1 0'
    new_sideset_name = 'hx_bot'
    input = 'hx_top'
  []
[]

################################################################################
# EQUATIONS: VARIABLES, KERNELS, BOUNDARY CONDITIONS
################################################################################

scalar_systems = 'prec1 prec2 prec3 prec4 prec5 prec6'

[Problem]
  # velocity, pressure restarted, precursors are not
  allow_initial_conditions_with_restart = true
  nl_sys_names = 'nl0 ${scalar_systems}'
[]

[Physics]
  [NavierStokes]
    [Flow/salt]
      # General parameters
      compressibility = 'incompressible'
      block = 'fuel pump hx'
      initialize_variables_from_mesh_file = true

      # Variables, defined below for the Exodus restart
      velocity_variable = 'vel_x vel_y'
      pressure_variable = 'pressure'
      solve_for_dynamic_pressure = true

      # Material properties
      density = ${rho}
      dynamic_viscosity = ${mu}
      thermal_expansion = ${alpha}

      # Boussinesq parameters
      boussinesq_approximation = true
      gravity = '0 -9.81 0'
      ref_temperature = ${T_HX}

      # Boundary conditions
      wall_boundaries = 'shield_wall reflector_wall fluid_symmetry'
      momentum_wall_types = 'wallfunction wallfunction symmetry'

      # Pressure pin for incompressible flow
      pin_pressure = true
      pinned_pressure_type = average
      # pressure pin for dynamic pressure: 0
      # pressure pin for total pressure: 1e5
      pinned_pressure_value = 0

      # Numerical scheme
      # time_derivative_contributes_to_RC_coefficients = false
      momentum_advection_interpolation = 'upwind'
      mass_advection_interpolation = 'upwind'

      # Heat exchanger friction
      standard_friction_formulation = false
      friction_blocks = 'hx'
      friction_types = 'FORCHHEIMER'
      friction_coeffs = 'friction_coeff_vector'
    []
    [FluidHeatTransfer/salt]
      block = 'fuel pump hx'
      fluid_temperature_variable = 'T_fluid'
      initial_temperature = ${T_HX}
      energy_advection_interpolation = 'upwind'

      # Material properties
      thermal_conductivity = ${k}
      specific_heat = 'cp'

      # See Flow for wall boundaries
      energy_wall_types = 'heatflux heatflux heatflux'
      energy_wall_functors = '0 0 0'

      # Heat source
      external_heat_source = power_density

      ambient_convection_blocks = 'hx'
      ambient_convection_alpha = '${fparse 600 * 20e3}' # HX specifications
      ambient_temperature = ${T_HX}

      # Numerical parameters
      energy_scaling = 100
    []
    [ScalarTransport/salt]
      block = 'fuel pump hx'

      # Precursor advection, diffusion and source term
      passive_scalar_names = 'c1 c2 c3 c4 c5 c6'
      passive_scalar_coupled_source = 'fission_source c1; fission_source c2; fission_source c3;
                                       fission_source c4; fission_source c5; fission_source c6'
      passive_scalar_coupled_source_coeff = '${beta1} ${lambda1_m}; ${beta2} ${lambda2_m};
                                             ${beta3} ${lambda3_m}; ${beta4} ${lambda4_m};
                                             ${beta5} ${lambda5_m}; ${beta6} ${lambda6_m}'

      # Numerical parameters
      passive_scalar_advection_interpolation = 'upwind'
      system_names = '${scalar_systems}'
    []
    [Turbulence/salt]
      block = 'fuel pump hx'

      fluid_heat_transfer_physics = salt
      scalar_transport_physics = salt

      # Turbulence parameters
      turbulence_handling = 'mixing-length'
      turbulent_prandtl = ${Pr_t}
      Sc_t = '${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t}'
      von_karman_const = ${von_karman_const}
      mixing_length_delta = 0.1
      mixing_length_walls = 'shield_wall reflector_wall'
      mixing_length_aux_execute_on = 'initial'
    []
  []
[]

[FVICs]
  [c1]
    type = FVFunctionIC
    variable = c1
    function = 'cosine_guess'
    scaling_factor = 0.02
  []
  [c2]
    type = FVFunctionIC
    variable = c2
    function = 'cosine_guess'
    scaling_factor = 0.1
  []
  [c3]
    type = FVFunctionIC
    variable = c3
    function = 'cosine_guess'
    scaling_factor = 0.03
  []
  [c4]
    type = FVFunctionIC
    variable = c4
    function = 'cosine_guess'
    scaling_factor = 0.04
  []
  [c5]
    type = FVFunctionIC
    variable = c5
    function = 'cosine_guess'
    scaling_factor = 0.01
  []
  [c6]
    type = FVFunctionIC
    variable = c6
    function = 'cosine_guess'
    scaling_factor = 0.001
  []
  # Power density is re-initalized by a transfer from neutronics
  [power]
    type = FVFunctionIC
    variable = power_density
    function = 'cosine_guess'
    scaling_factor = '${fparse 3e9/2.81543}'
  []
  # Fission source is re-initalized by a transfer from neutronics
  [fission_source]
    type = FVFunctionIC
    variable = fission_source
    function = 'cosine_guess'
    scaling_factor = '${fparse 6.303329e+01/2.81543}'
  []
[]

[AuxVariables]
  [power_density]
    type = MooseVariableFVReal
    block = 'fuel pump hx'
  []
  [fission_source]
    type = MooseVariableFVReal
    block = 'fuel pump hx'
  []
[]

[Functions]
  # Guess to have a 3D power distribution
  [cosine_guess]
    type = ParsedFunction
    expression = 'max(0, cos(x*pi/2/1.2))*max(0, cos(y*pi/2/1.1))'
  []
[]

[FVKernels]
  [pump]
    type = INSFVBodyForce
    variable = vel_y
    functor = ${pump_force}
    block = 'pump'
    momentum_component = 'y'
  []
[]

################################################################################
# MATERIALS
################################################################################

[FunctorMaterials]
  # Most of these constants could be specified directly to the action
  [mu]
    type = ADGenericFunctorMaterial
    prop_names = 'mu'
    prop_values = '${mu}'
    block = 'fuel pump hx'
  []
  # [not_used]
  #   type = ADGenericFunctorMaterial
  #   prop_names = 'not_used'
  #   prop_values = 0
  #   block = 'shield reflector'
  # []
  [cp]
    type = ADGenericFunctorMaterial
    prop_names = 'cp'
    prop_values = '${cp}'
    block = 'fuel pump hx'
  []
  [friction_coeff_mat]
    type = ADGenericVectorFunctorMaterial
    prop_names = 'friction_coeff_vector'
    prop_values = '${friction} ${friction} ${friction}'
  []
[]

################################################################################
# EXECUTION / SOLVE
################################################################################

[Functions]
  [dts]
    type = PiecewiseConstant
    x = '0 10  100'
    y = '2 4 5'
  []
[]

[Executioner]
  type = Transient

  # Time stepping parameters
  start_time = 0.0
  end_time = 50
  # end_time will depend on the restart file chosen
  # though steady state detection can also be used
  # from _initial/no heating : 150 - 200s enough
  # from _ns/_ns_coupled/heated: 10s enough

  # [TimeStepper]
  #   # This time stepper makes the time step grow exponentially
  #   # It can only be used with proper initialization
  #   type = IterationAdaptiveDT
  #   dt = 1 # chosen to obtain convergence with first coupled iteration
  #   growth_factor = 2
  # []
  [TimeStepper]
    type = FunctionDT
    function = dts
  []
  steady_state_detection = true
  steady_state_tolerance = 1e-8
  steady_state_start_time = 10

  # Time integration scheme
  scheme = 'implicit-euler'

  # Solver parameters
  solve_type = 'NEWTON'
  line_search = 'none'

  # nonlinear solver parameters
  nl_rel_tol = 2e-8
  nl_abs_tol = 2e-8
  nl_abs_div_tol = 1e11
  nl_max_its = 15

  # linear solver parameters
  l_max_its = 50

  automatic_scaling = true
[]

# Try a direct solve. The precursor advection problem should be linear
petsc_options_iname_prec = '-pc_type -pc_factor_shift_type'
petsc_options_value_prec = 'lu NONZERO'

[Preconditioning]
  [flow]
    type = SMP
    full = true
    nl_sys = "nl0"
    petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_mat_solver_package'
    petsc_options_value = 'lu NONZERO superlu_dist'
  []
  [scalar1]
    type = SMP
    nl_sys = "prec1"
    petsc_options_iname = ${petsc_options_iname_prec}
    petsc_options_value = ${petsc_options_value_prec}
  []
  [scalar2]
    type = SMP
    nl_sys = "prec2"
    petsc_options_iname = ${petsc_options_iname_prec}
    petsc_options_value = ${petsc_options_value_prec}
  []
  [scalar3]
    type = SMP
    nl_sys = "prec3"
    petsc_options_iname = ${petsc_options_iname_prec}
    petsc_options_value = ${petsc_options_value_prec}
  []
  [scalar4]
    type = SMP
    nl_sys = "prec4"
    petsc_options_iname = ${petsc_options_iname_prec}
    petsc_options_value = ${petsc_options_value_prec}
  []
  [scalar5]
    type = SMP
    nl_sys = "prec5"
    petsc_options_iname = ${petsc_options_iname_prec}
    petsc_options_value = ${petsc_options_value_prec}
  []
  [scalar6]
    type = SMP
    nl_sys = "prec6"
    petsc_options_iname = ${petsc_options_iname_prec}
    petsc_options_value = ${petsc_options_value_prec}
  []
[]

################################################################################
# SIMULATION OUTPUTS
################################################################################

[Outputs]
  csv = true
  hide = 'flow_hx_bot flow_hx_top min_flow_T max_flow_T'
  [restart]
    type = Exodus
    overwrite = true
  []
  # Reduce base output
  print_linear_converged_reason = false
  print_linear_residuals = false
  print_nonlinear_converged_reason = false
[]

[Postprocessors]
  [max_v]
    type = ElementExtremeValue
    variable = vel_x
    value_type = max
    block = 'fuel pump hx'
  []
  # TODO: weakly compressible, switch to mass flow rate
  [flow_hx_bot]
    type = VolumetricFlowRate
    boundary = 'hx_bot'
    vel_x = vel_x
    vel_y = vel_y
    advected_quantity = 1
  []
  [flow_hx_top]
    type = VolumetricFlowRate
    boundary = 'hx_top'
    vel_x = vel_x
    vel_y = vel_y
    advected_quantity = 1
  []
  [max_flow_T]
    type = VolumetricFlowRate
    boundary = 'hx_top'
    vel_x = vel_x
    vel_y = vel_y
    advected_quantity = 'T_fluid'
  []
  [min_flow_T]
    type = VolumetricFlowRate
    boundary = 'hx_bot'
    vel_x = vel_x
    vel_y = vel_y
    advected_quantity = 'T_fluid'
  []
  [dT]
    type = ParsedPostprocessor
    expression = '-max_flow_T / flow_hx_bot + min_flow_T / flow_hx_top'
    pp_names = 'max_flow_T min_flow_T flow_hx_bot flow_hx_top'
  []
  [total_power]
    type = ElementIntegralVariablePostprocessor
    variable = power_density
    block = 'fuel pump hx'
  []
  [total_fission_source]
    type = ElementIntegralVariablePostprocessor
    variable = fission_source
    block = 'fuel pump hx'
  []
[]

(msr/msfr/steady/legacy/run_ns.i)

################################################################################
## Molten Salt Fast Reactor - Euratom EVOL + Rosatom MARS Design              ##
## Pronghorn Sub-Application input file                                       ##
## Relaxation towards Steady state 3D thermal hydraulics model                ##
################################################################################
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html

## WARNING

# This input is present for documentation purposes. It is not actively
# maintained and should not be used under any circumstance.

## WARNING

# The flow in this simulation should be initialized with a previous flow
# solution (isothermal, heated or multiphysics) OR using a viscosity rampdown
# An isothermal solution can be generated using 'run_ns_initial.i', and is
# saved in 'restart/run_ns_initial_restart.e'
# A heated solution, with a flat power distribution, can be generated with
# this script 'run_ns.i', and is saved in 'restart/run_ns_restart.e'
# A coupled neutronics-coarse TH solution can be generated with
# 'run_neutronics.i', saved in 'restart/run_neutronics_ns_restart.e'

# Material properties
rho = 4284 # density [kg / m^3]  (@1000K)
cp = 1594 # specific heat capacity [J / kg / K]
drho_dT = 0.882 # derivative of density w.r.t. temperature [kg / m^3 / K]
mu = 0.0166 # viscosity [Pa s]
k = 1.7 # thermal conductivity [W / m / K]
# https://www.researchgate.net/publication/337161399_Development_of_a_control-\
# oriented_power_plant_simulator_for_the_molten_salt_fast_reactor/fulltext/5dc9\
# 5c9da6fdcc57503eec39/Development-of-a-control-oriented-power-plant-simulator-\
# for-the-molten-salt-fast-reactor.pdf

# Turbulent properties
Pr_t = 10 # turbulent Prandtl number
Sc_t = 1 # turbulent Schmidt number

# Derived material properties
alpha = '${fparse drho_dT / rho}' # thermal expansion coefficient

# Operating parameters
T_HX = 873.15 # heat exchanger temperature [K]

# Mass flow rate tuning, for heat exchanger pressure and temperature drop
friction = 4e3 # [kg / m^4]
pump_force = -20000. # [N / m^3]

# Delayed neutron precursor parameters. Lambda values are decay constants in
# [1 / s]. Beta values are production fractions.
lambda1 = 0.0133104
lambda2 = 0.0305427
lambda3 = 0.115179
lambda4 = 0.301152
lambda5 = 0.879376
lambda6 = 2.91303
beta1 = 8.42817e-05
beta2 = 0.000684616
beta3 = 0.000479796
beta4 = 0.00103883
beta5 = 0.000549185
beta6 = 0.000184087

[GlobalParams]
  u = v_x
  v = v_y
  pressure = pressure
  temperature = T

  rhie_chow_user_object = 'rc'
  advected_interp_method = 'upwind'
  velocity_interp_method = 'rc'
  mu = 'mu'
  rho = ${rho}
  mixing_length = 'mixing_len'
[]

################################################################################
# GEOMETRY
################################################################################

[Mesh]
  [restart]
    type = FileMeshGenerator
    use_for_exodus_restart = true
    # Depending on the file chosen, the initialization of variables should be
    # adjusted. The following variables can be initalized:
    # - v_x, v_y, p, lambda from isothermal simulation
    # file = 'restart/run_ns_initial_restart.e'
    # - v_x, v_y, p, T, lambda, c_ifrom cosine heated simulation
    file = 'restart/run_ns_restart.e'
    # - v_x, v_y, p, T, lambda, c_i from coupled multiphysics simulation
    # file = 'restart/run_ns_coupled_restart.e'
  []
  [min_radius]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'abs(y-0.) < 1e-10 & x < 1.8'
    input = restart
    new_sideset_name = min_core_radius
    normal = '0 1 0'
  []
  [hx_top]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'y > 0'
    included_subdomains = '3'
    included_neighbors = '1'
    fixed_normal = true
    normal = '0 1 0'
    new_sideset_name = 'hx_top'
    input = 'min_radius'
  []
  [hx_bot]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'y <-0.6'
    included_subdomains = '3'
    included_neighbors = '1'
    fixed_normal = true
    normal = '0 -1 0'
    new_sideset_name = 'hx_bot'
    input = 'hx_top'
  []
[]

[Problem]
  coord_type = 'RZ'
  rz_coord_axis = Y
[]

################################################################################
# EQUATIONS: VARIABLES, KERNELS
################################################################################

[Variables]
  [v_x]
    type = INSFVVelocityVariable
    block = 'fuel pump hx'
    initial_from_file_var = v_x
  []
  [v_y]
    type = INSFVVelocityVariable
    block = 'fuel pump hx'
    initial_from_file_var = v_y
  []
  [pressure]
    type = INSFVPressureVariable
    block = 'fuel pump hx'
    initial_from_file_var = pressure
  []
  [lambda]
    family = SCALAR
    order = FIRST
    block = 'fuel pump hx'
    initial_from_file_var = lambda
  []
  [T]
    type = MooseVariableFVReal
    block = 'fuel pump hx'
    initial_condition = ${T_HX}
    # initial_from_file_var = T
  []
  [c1]
    type = MooseVariableFVReal
    block = 'fuel pump hx'
    # initial_from_file_var = c1
  []
  [c2]
    type = MooseVariableFVReal
    block = 'fuel pump hx'
    # initial_from_file_var = c2
  []
  [c3]
    type = MooseVariableFVReal
    block = 'fuel pump hx'
    # initial_from_file_var = c3
  []
  [c4]
    type = MooseVariableFVReal
    block = 'fuel pump hx'
    # initial_from_file_var = c4
  []
  [c5]
    type = MooseVariableFVReal
    block = 'fuel pump hx'
    # initial_from_file_var = c5
  []
  [c6]
    type = MooseVariableFVReal
    block = 'fuel pump hx'
    # initial_from_file_var = c6
  []
[]

[FVICs]
  [c1]
    type = FVFunctionIC
    variable = c1
    function = 'cosine_guess'
    scaling_factor = 0.02
  []
  [c2]
    type = FVFunctionIC
    variable = c2
    function = 'cosine_guess'
    scaling_factor = 0.1
  []
  [c3]
    type = FVFunctionIC
    variable = c3
    function = 'cosine_guess'
    scaling_factor = 0.03
  []
  [c4]
    type = FVFunctionIC
    variable = c4
    function = 'cosine_guess'
    scaling_factor = 0.04
  []
  [c5]
    type = FVFunctionIC
    variable = c5
    function = 'cosine_guess'
    scaling_factor = 0.01
  []
  [c6]
    type = FVFunctionIC
    variable = c6
    function = 'cosine_guess'
    scaling_factor = 0.001
  []
  # Power density is re-initalized by a transfer from neutronics
  [power]
    type = FVFunctionIC
    variable = power_density
    function = 'cosine_guess'
    scaling_factor = '${fparse 3e9/2.81543}'
  []
  # Fission source is re-initalized by a transfer from neutronics
  [fission_source]
    type = FVFunctionIC
    variable = fission_source
    function = 'cosine_guess'
    scaling_factor = '${fparse 6.303329e+01/2.81543}'
  []
[]

[AuxVariables]
  [mixing_len]
    type = MooseVariableFVReal
    initial_from_file_var = mixing_len
    block = 'fuel pump hx'
  []
  [power_density]
    type = MooseVariableFVReal
    block = 'fuel pump hx'
  []
  [fission_source]
    type = MooseVariableFVReal
    block = 'fuel pump hx'
  []

  # For visualization purposes
  [eddy_viscosity]
    type = MooseVariableFVReal
    block = 'fuel pump hx'
  []
  [wall_shear_stress]
    type = MooseVariableFVReal
    block = 'fuel pump hx'
  []
  [wall_yplus]
    type = MooseVariableFVReal
    block = 'fuel pump hx'
  []
[]

[Functions]
  # Guess to have a 3D power distribution
  [cosine_guess]
    type = ParsedFunction
    value = 'max(0, cos(x*pi/2/1.2))*max(0, cos(y*pi/2/1.1))'
  []
[]

[UserObjects]
  [rc]
    type = INSFVRhieChowInterpolator
    block = 'fuel pump hx'
  []
[]

[FVKernels]
  [mass]
    type = INSFVMassAdvection
    variable = pressure
    block = 'fuel pump hx'
  []
  [mean_zero_pressure]
    type = FVScalarLagrangeMultiplier
    variable = pressure
    lambda = lambda
    block = 'fuel pump hx'
  []

  [u_time]
    type = INSFVMomentumTimeDerivative
    variable = v_x
    momentum_component = 'x'
  []
  [u_advection]
    type = INSFVMomentumAdvection
    variable = v_x
    block = 'fuel pump hx'
    momentum_component = 'x'
  []
  [u_turbulent_diffusion_rans]
    type = INSFVMixingLengthReynoldsStress
    variable = v_x
    momentum_component = 'x'
  []
  [u_molecular_diffusion]
    type = INSFVMomentumDiffusion
    variable = v_x
    block = 'fuel pump hx'
    momentum_component = 'x'
  []
  [u_pressure]
    type = INSFVMomentumPressure
    variable = v_x
    momentum_component = 'x'
    block = 'fuel pump hx'
  []

  [v_time]
    type = INSFVMomentumTimeDerivative
    variable = v_y
    momentum_component = 'y'
  []
  [v_advection]
    type = INSFVMomentumAdvection
    variable = v_y
    block = 'fuel pump hx'
    momentum_component = 'y'
  []
  [v_turbulent_diffusion_rans]
    type = INSFVMixingLengthReynoldsStress
    variable = v_y
    momentum_component = 'y'
  []
  [v_molecular_diffusion]
    type = INSFVMomentumDiffusion
    variable = v_y
    block = 'fuel pump hx'
    momentum_component = 'y'
  []
  [v_pressure]
    type = INSFVMomentumPressure
    variable = v_y
    momentum_component = 'y'
    block = 'fuel pump hx'
  []
  [v_buoyancy]
    type = INSFVMomentumBoussinesq
    variable = v_y
    T_fluid = T
    gravity = '0 -9.81 0'
    ref_temperature = 1000
    momentum_component = 'y'
    block = 'fuel'
    alpha_name = 'alpha_b'
  []
  [v_gravity]
    type = INSFVMomentumGravity
    variable = v_y
    gravity = '0 -9.81 0'
    block = 'fuel pump hx'
    momentum_component = 'y'
  []

  [pump]
    type = INSFVBodyForce
    variable = v_y
    functor = ${pump_force}
    block = 'pump'
    momentum_component = 'y'
  []

  [friction_hx_x]
    type = INSFVMomentumFriction
    variable = v_x
    quadratic_coef_name = 'friction_coef'
    block = 'hx'
    momentum_component = 'x'
  []
  [friction_hx_y]
    type = INSFVMomentumFriction
    variable = v_y
    quadratic_coef_name = 'friction_coef'
    block = 'hx'
    momentum_component = 'y'
  []

  [heat_time]
    type = INSFVEnergyTimeDerivative
    variable = T
    cp_name = 'cp'
  []
  [heat_advection]
    type = INSFVEnergyAdvection
    variable = T
    block = 'fuel pump hx'
  []
  [heat_diffusion]
    type = FVDiffusion
    coeff = '${k}'
    variable = T
    block = 'fuel pump hx'
  []
  [heat_turb_diffusion]
    type = WCNSFVMixingLengthEnergyDiffusion
    schmidt_number = ${Pr_t}
    variable = T
    block = 'fuel pump hx'
    cp = 'cp'
  []
  [heat_src]
    type = FVCoupledForce
    variable = T
    v = power_density
    block = 'fuel pump hx'
  []
  [heat_sink]
    type = NSFVEnergyAmbientConvection
    variable = T
    # Compute the coefficient as 600 m^2 / m^3 surface area density times a heat
    # transfer coefficient of 20 kW / m^2 / K
    alpha = 'alpha'
    block = 'hx'
    T_ambient = ${T_HX}
  []

  [c1_advection]
    type = INSFVScalarFieldAdvection
    variable = c1
    block = 'fuel pump hx'
  []
  [c2_advection]
    type = INSFVScalarFieldAdvection
    variable = c2
    block = 'fuel pump hx'
  []
  [c3_advection]
    type = INSFVScalarFieldAdvection
    variable = c3
    block = 'fuel pump hx'
  []
  [c4_advection]
    type = INSFVScalarFieldAdvection
    variable = c4
    block = 'fuel pump hx'
  []
  [c5_advection]
    type = INSFVScalarFieldAdvection
    variable = c5
    block = 'fuel pump hx'
  []
  [c6_advection]
    type = INSFVScalarFieldAdvection
    variable = c6
    block = 'fuel pump hx'
  []
  [c1_turb_diffusion]
    type = INSFVMixingLengthScalarDiffusion
    schmidt_number = ${Sc_t}
    variable = c1
    block = 'fuel pump hx'
  []
  [c2_turb_diffusion]
    type = INSFVMixingLengthScalarDiffusion
    schmidt_number = ${Sc_t}
    variable = c2
    block = 'fuel pump hx'
  []
  [c3_turb_diffusion]
    type = INSFVMixingLengthScalarDiffusion
    schmidt_number = ${Sc_t}
    variable = c3
    block = 'fuel pump hx'
  []
  [c4_turb_diffusion]
    type = INSFVMixingLengthScalarDiffusion
    schmidt_number = ${Sc_t}
    variable = c4
    block = 'fuel pump hx'
  []
  [c5_turb_diffusion]
    type = INSFVMixingLengthScalarDiffusion
    schmidt_number = ${Sc_t}
    variable = c5
    block = 'fuel pump hx'
  []
  [c6_turb_diffusion]
    type = INSFVMixingLengthScalarDiffusion
    schmidt_number = ${Sc_t}
    variable = c6
    block = 'fuel pump hx'
  []
  [c1_src]
    type = FVCoupledForce
    variable = c1
    v = fission_source
    coef = ${beta1}
    block = 'fuel pump hx'
  []
  [c2_src]
    type = FVCoupledForce
    variable = c2
    v = fission_source
    coef = ${beta2}
    block = 'fuel pump hx'
  []
  [c3_src]
    type = FVCoupledForce
    variable = c3
    v = fission_source
    coef = ${beta3}
    block = 'fuel pump hx'
  []
  [c4_src]
    type = FVCoupledForce
    variable = c4
    v = fission_source
    coef = ${beta4}
    block = 'fuel pump hx'
  []
  [c5_src]
    type = FVCoupledForce
    variable = c5
    v = fission_source
    coef = ${beta5}
    block = 'fuel pump hx'
  []
  [c6_src]
    type = FVCoupledForce
    variable = c6
    v = fission_source
    coef = ${beta6}
    block = 'fuel pump hx'
  []
  [c1_decay]
    type = FVReaction
    variable = c1
    rate = ${lambda1}
    block = 'fuel pump hx'
  []
  [c2_decay]
    type = FVReaction
    variable = c2
    rate = ${lambda2}
    block = 'fuel pump hx'
  []
  [c3_decay]
    type = FVReaction
    variable = c3
    rate = ${lambda3}
    block = 'fuel pump hx'
  []
  [c4_decay]
    type = FVReaction
    variable = c4
    rate = ${lambda4}
    block = 'fuel pump hx'
  []
  [c5_decay]
    type = FVReaction
    variable = c5
    rate = ${lambda5}
    block = 'fuel pump hx'
  []
  [c6_decay]
    type = FVReaction
    variable = c6
    rate = ${lambda6}
    block = 'fuel pump hx'
  []
[]

[AuxKernels]
  [mixing_len]
    type = WallDistanceMixingLengthAux
    walls = 'shield_wall reflector_wall'
    variable = mixing_len
    execute_on = 'initial'
    block = 'fuel'
    von_karman_const = ${von_karman_const}
    delta = 0.1 # capping parameter (m)
  []
  [wall_shear_stress]
    type = WallFunctionWallShearStressAux
    variable = wall_shear_stress
    walls = 'shield_wall reflector_wall'
    block = 'fuel'
    mu = 'total_viscosity'
  []
  [wall_yplus]
    type = WallFunctionYPlusAux
    variable = wall_yplus
    walls = 'shield_wall reflector_wall'
    block = 'fuel'
    mu = 'total_viscosity'
  []
  [turbulent_viscosity]
    type = INSFVMixingLengthTurbulentViscosityAux
    variable = eddy_viscosity
    block = 'fuel pump hx'
    execute_on = 'initial'
  []
[]

################################################################################
# BOUNDARY CONDITIONS
################################################################################

[FVBCs]
  [walls_u]
    type = INSFVWallFunctionBC
    variable = v_x
    boundary = 'shield_wall reflector_wall'
    momentum_component = 'x'
  []
  [walls_v]
    type = INSFVWallFunctionBC
    variable = v_y
    boundary = 'shield_wall reflector_wall'
    momentum_component = 'y'
  []

  [symmetry_u]
    type = INSFVSymmetryVelocityBC
    variable = v_x
    boundary = 'fluid_symmetry'
    momentum_component = 'x'
    mu = 'total_viscosity'
  []
  [symmetry_v]
    type = INSFVSymmetryVelocityBC
    variable = v_y
    boundary = 'fluid_symmetry'
    momentum_component = 'y'
    mu = 'total_viscosity'
  []
  [symmetry_pressure]
    type = INSFVSymmetryPressureBC
    boundary = 'fluid_symmetry'
    variable = pressure
  []
[]

################################################################################
# MATERIALS
################################################################################

[Materials]
  [heat_exchanger_coefficient]
    type = ADGenericFunctionMaterial
    prop_names = 'alpha'
    prop_values = '${fparse 600 * 20e3}'
    block = 'fuel pump hx'
  []
[]

[FunctorMaterials]
  [boussinesq]
    type = ADGenericFunctorMaterial
    prop_names = 'alpha_b'
    prop_values = '${alpha}'
  []
  [mu]
    type = ADGenericFunctorMaterial
    prop_names = 'mu'
    prop_values = '${mu}'
    block = 'fuel pump hx'
  []
  [total_viscosity]
    type = MixingLengthTurbulentViscosityFunctorMaterial
    mu = 'mu'
    block = 'fuel pump hx'
  []
  [ins_fv]
    type = INSFVEnthalpyFunctorMaterial
    block = 'fuel pump hx'
  []
  # [not_used]
  #   type = ADGenericFunctorMaterial
  #   prop_names = 'not_used'
  #   prop_values = 0
  #   block = 'shield reflector'
  # []
  [friction]
    type = ADGenericFunctorMaterial
    prop_names = 'friction_coef'
    prop_values = '${friction} '
    block = 'hx'
  []
  [cp]
    type = ADGenericFunctorMaterial
    prop_names = 'cp'
    prop_values = '${cp}'
    block = 'fuel pump hx'
  []
[]

################################################################################
# EXECUTION / SOLVE
################################################################################

[Functions]
  [dts]
    type = PiecewiseConstant
    x = '0    100'
    y = '0.75 2.5'
  []
[]

[Executioner]
  type = Transient

  # Time stepping parameters
  start_time = 0.0
  end_time = 20
  # [TimeStepper]
  #   # This time stepper makes the time step grow exponentially
  #   # It can only be used with proper initialization
  #   type = IterationAdaptiveDT
  #   dt = 10  # chosen to obtain convergence with first coupled iteration
  #   growth_factor = 2
  # []
  [TimeStepper]
    type = FunctionDT
    function = dts
  []
  steady_state_detection = true
  steady_state_tolerance = 1e-8
  steady_state_start_time = 10

  # Solver parameters
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -pc_factor_shift_type -ksp_gmres_restart'
  petsc_options_value = 'lu NONZERO 50'
  line_search = 'none'

  nl_rel_tol = 1e-9
  nl_abs_tol = 2e-8
  nl_max_its = 15
  l_max_its = 50

  automatic_scaling = true
  resid_vs_jac_scaling_param = 1
[]

################################################################################
# SIMULATION OUTPUTS
################################################################################

[Outputs]
  exodus = true
  csv = true
  hide = 'flow_hx_bot flow_hx_top min_flow_T max_flow_T'
  [restart]
    type = Exodus
    overwrite = true
  []
  # Reduce base output
  print_linear_converged_reason = false
  print_linear_residuals = false
  print_nonlinear_converged_reason = false
[]

[Postprocessors]
  [max_v]
    type = ElementExtremeValue
    variable = v_x
    value_type = max
    block = 'fuel pump hx'
  []
  [mdot]
    type = InternalVolumetricFlowRate
    boundary = 'min_core_radius'
    vel_x = v_x
    vel_y = v_y
    advected_mat_prop = ${rho}
  []
  # TODO: weakly compressible, switch to mass flow rate
  [flow_hx_bot]
    type = InternalVolumetricFlowRate
    boundary = 'hx_bot'
    vel_x = v_x
    vel_y = v_y
  []
  [flow_hx_top]
    type = InternalVolumetricFlowRate
    boundary = 'hx_top'
    vel_x = v_x
    vel_y = v_y
  []
  [max_flow_T]
    type = InternalVolumetricFlowRate
    boundary = 'hx_top'
    vel_x = v_x
    vel_y = v_y
    advected_variable = 'T'
  []
  [min_flow_T]
    type = InternalVolumetricFlowRate
    boundary = 'hx_bot'
    vel_x = v_x
    vel_y = v_y
    advected_variable = 'T'
  []
  [dT]
    type = ParsedPostprocessor
    function = '-max_flow_T / flow_hx_bot + min_flow_T / flow_hx_top'
    pp_names = 'max_flow_T min_flow_T flow_hx_bot flow_hx_top'
  []
  [total_power]
    type = ElementIntegralVariablePostprocessor
    variable = power_density
    block = 'fuel pump hx'
  []
  [total_fission_source]
    type = ElementIntegralVariablePostprocessor
    variable = fission_source
    block = 'fuel pump hx'
  []
[]

(msr/msfr/steady/run_neutronics.i)

################################################################################
## Molten Salt Fast Reactor - Euratom EVOL + Rosatom MARS Design              ##
## Griffin Main Application input file                                        ##
## Steady state neutronics model                                              ##
## Neutron diffusion with delayed precursor source, no equivalence            ##
################################################################################
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html

!include run_neutronics_base.i

[MSRNeutronicsFluidCoupling]
  fluid_blocks = 'fuel pump hx'
  solid_blocks = 'reflector shield'
  n_dnp = 6
  use_transient_multi_app = false
  fluid_app_name = ns
  fluid_input_file = run_ns.i
  initial_fluid_temperature = 873.15
  fluid_temperature_name_neutronics_app = tfuel
  fluid_temperature_name_solid_suffix = _constant
  fluid_temperature_name_fluid_app = T_fluid
  dnp_name_prefix = c
  dnp_array_name = dnp
[]

(msr/msfr/steady/run_neutronics.i)

################################################################################
## Molten Salt Fast Reactor - Euratom EVOL + Rosatom MARS Design              ##
## Griffin Main Application input file                                        ##
## Steady state neutronics model                                              ##
## Neutron diffusion with delayed precursor source, no equivalence            ##
################################################################################
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html

!include run_neutronics_base.i

[MSRNeutronicsFluidCoupling]
  fluid_blocks = 'fuel pump hx'
  solid_blocks = 'reflector shield'
  n_dnp = 6
  use_transient_multi_app = false
  fluid_app_name = ns
  fluid_input_file = run_ns.i
  initial_fluid_temperature = 873.15
  fluid_temperature_name_neutronics_app = tfuel
  fluid_temperature_name_solid_suffix = _constant
  fluid_temperature_name_fluid_app = T_fluid
  dnp_name_prefix = c
  dnp_array_name = dnp
[]

(msr/msfr/steady/run_neutronics.i)

################################################################################
## Molten Salt Fast Reactor - Euratom EVOL + Rosatom MARS Design              ##
## Griffin Main Application input file                                        ##
## Steady state neutronics model                                              ##
## Neutron diffusion with delayed precursor source, no equivalence            ##
################################################################################
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html

!include run_neutronics_base.i

[MSRNeutronicsFluidCoupling]
  fluid_blocks = 'fuel pump hx'
  solid_blocks = 'reflector shield'
  n_dnp = 6
  use_transient_multi_app = false
  fluid_app_name = ns
  fluid_input_file = run_ns.i
  initial_fluid_temperature = 873.15
  fluid_temperature_name_neutronics_app = tfuel
  fluid_temperature_name_solid_suffix = _constant
  fluid_temperature_name_fluid_app = T_fluid
  dnp_name_prefix = c
  dnp_array_name = dnp
[]

(msr/msfr/steady/run_neutronics.i)

################################################################################
## Molten Salt Fast Reactor - Euratom EVOL + Rosatom MARS Design              ##
## Griffin Main Application input file                                        ##
## Steady state neutronics model                                              ##
## Neutron diffusion with delayed precursor source, no equivalence            ##
################################################################################
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html

!include run_neutronics_base.i

[MSRNeutronicsFluidCoupling]
  fluid_blocks = 'fuel pump hx'
  solid_blocks = 'reflector shield'
  n_dnp = 6
  use_transient_multi_app = false
  fluid_app_name = ns
  fluid_input_file = run_ns.i
  initial_fluid_temperature = 873.15
  fluid_temperature_name_neutronics_app = tfuel
  fluid_temperature_name_solid_suffix = _constant
  fluid_temperature_name_fluid_app = T_fluid
  dnp_name_prefix = c
  dnp_array_name = dnp
[]

(msr/msfr/steady/run_neutronics_base.i)

################################################################################
## Molten Salt Fast Reactor - Euratom EVOL + Rosatom MARS Design              ##
## Griffin Main Application base input file                                   ##
## Steady state neutronics model                                              ##
## Neutron diffusion with delayed precursor source, no equivalence            ##
################################################################################
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html

[Mesh]
  uniform_refine = 1
  coord_type = 'RZ'
  [fmg]
    type = FileMeshGenerator
    file = '../mesh/msfr_rz_mesh.e'
  []
[]

[TransportSystems]
  particle = neutron
  equation_type = eigenvalue

  G = 6

  VacuumBoundary = 'outer'

  [diff]
    scheme = CFEM-Diffusion
    n_delay_groups = 6
    external_dnp_variable = 'dnp'
    family = LAGRANGE
    order = FIRST
    fission_source_aux = true

    # For PJFNKMO
    assemble_scattering_jacobian = true
    assemble_fission_jacobian = true
  []
[]

[PowerDensity]
  power = 3e9
  power_density_variable = power_density
  power_scaling_postprocessor = power_scaling
  family = MONOMIAL
  order = CONSTANT
[]

################################################################################
# CROSS SECTIONS
################################################################################

[Materials]
  [fuel]
    type = CoupledFeedbackNeutronicsMaterial
    library_name = 'msfr_xs'
    library_file = '../mgxs/msfr_xs.xml'
    grid_names = 'tfuel'
    grid_variables = 'tfuel'
    plus = true
    isotopes = 'pseudo'
    densities = '1.0'
    is_meter = true
    material_id = 2
    block = 'fuel pump hx'
  []
  [shield]
    type = CoupledFeedbackNeutronicsMaterial
    library_name = 'msfr_xs'
    library_file = '../mgxs/msfr_xs.xml'
    grid_names = 'tfuel'
    grid_variables = 'tfuel_constant'
    isotopes = 'pseudo'
    densities = '1.0'
    is_meter = true
    material_id = 5
    block = 'shield'
  []
  [reflector]
    type = CoupledFeedbackNeutronicsMaterial
    library_name = 'msfr_xs'
    library_file = '../mgxs/msfr_xs.xml'
    grid_names = 'tfuel'
    grid_variables = 'tfuel_constant'
    isotopes = 'pseudo'
    densities = '1.0'
    is_meter = true
    material_id = 1
    block = 'reflector'
  []
[]

################################################################################
# EXECUTION / SOLVE
################################################################################

[Executioner]
  type = Eigenvalue
  solve_type = PJFNKMO

  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart '
  petsc_options_value = 'hypre boomeramg 50'
  l_max_its = 50

  free_power_iterations = 4  # important to obtain fundamental mode eigenvalue

  nl_abs_tol = 1e-9
  fixed_point_abs_tol = 1e-9
  fixed_point_max_its = 20
[]

################################################################################
# SIMULATION OUTPUTS
################################################################################

[Outputs]
  exodus = true
  csv = true
  checkpoint = true
  [restart]
    type = Exodus
    execute_on = 'final'
    file_base = 'run_neutronics_restart'
  []
  print_linear_converged_reason = false
  print_linear_residuals = false
  print_nonlinear_converged_reason = false
  hide = 'power_scaling'
[]

Conservation of fluid mass
Conservation of fluid momentum
Conservation of fluid energy
Conservation of delayed neutron precursors
Neutronics
References