MultiApps problems with Subchannel Solver

[garcs2]

Hello,

I am attempting to create a multiphysics model for a hexagonal lattice liquid-metal cooled microeactor, using BISON for heat conduction and Subchannel for the thermal hydraulics. I have been having issues with Subchannel not solving and throwing out errors when attempting to couple the two. For debugging purposes I have: 1) ensured that the linear heat rate (LHR) in BISON is identical to what is being sent to in Subchannel 2) successfully ran standalone Subchannel cases with constant and non-constant LHR and 3) ensured that the solutions in the standalone Subchannel make sense. At this point I believed it might have been a transfer issue but then attempted to run a multiapp but have no transfers done. This is to say that I had the working standalone subchannel model run as the child app of the BISON model. From what I understand this should effectively amount to having two physics modules run one after the other but with no transfers they are not connected. Since the standalone SC model works when not in MultiApp I figured it should work here as well but I run into the same error being thrown. This to me points to their being an issue with the timing of execution of some process within subchannel. I would appreciate any insight on this! Below are my inputs for the BISON model and the SC model. Note that for the SC model, I've also attempted to not use AuxKernels and set the initial conditions manually but it yielded the same error. Below is the error I get when attempting to run MultiApps:

sc0: *** ERROR ***
sc0: /home/garcsamu/Serpent/SNAP-REACTORS-PRIVATE/snapReactors/reactor_models/Wet_Experiment_Models/standard_conditions/sc_test/pow_den/sc_core.i:168.1:
sc0: The following error occurred in the Problem 'MOOSE Problem' of type TriSubChannel1PhaseProblem.
sc0: 
sc0: MOOSE Problem : Calculation of negative Enthalpy h_out = : -5.14372e+06 Axial Level= : 1

Below are my models for Bison and SC.

Bison Model:

inlet_T_fluid             = 860.0  # (K) 
ht_coeff                  = 4539.6
fuel_blocks = 'Fuel'
ceram_blocks = 'Ceramic'
clad_blocks = 'Clad'
intref_blocks = 'Reflector'
# ==============================================================================
# GEOMETRY AND MESH
# ==============================================================================
[Mesh]
    [core_unextruded]
        type = FileMeshGenerator
        file = /home/garcsamu/Serpent/SNAP-REACTORS-PRIVATE/snapReactors/reactor_models/Wet_Experiment_Models/meshes/SNAP_mesh_1_in.e
    []
[]
# ==============================================================================
# VARIABLES AND KERNELS
# ==============================================================================
[Variables]
    [bison_temp]
        initial_condition = '${inlet_T_fluid}'
    []
[]

[Kernels]
    # Equation 0 (Heat conduction equation with heat generation).
    [heat_conduction]
        type = ADHeatConduction
        variable = bison_temp
    []
    [heat_source]
        type = ADCoupledForce
        variable = bison_temp
        v = bison_norm_power_density
        block = 'Fuel'
        #coef = .93
    []
    # [heat_source_extref]
    #     type = CoupledForce
    #     variable = bison_temp
    #     v = bison_norm_power_density
    #     block = 'Reflector'
    #     coef = .07
    # []
[]

# ==============================================================================
# AUXVARIABLES AND AUXKERNELS
# ==============================================================================
[AuxVariables]
    [bison_power_density]
        family = L2_LAGRANGE 
        order = FIRST 
        block = 'Fuel'
        #initial_condition = '${pow_dens}' #
    []
    [bison_T_inf]
        initial_condition = '${inlet_T_fluid}'
    []
    [bison_Tfuel]
        block = ${fuel_blocks}
    []

    # [bison_Tref]
    #     block = ${extref_blocks}
    # []
    [bison_norm_power_density]
        family = L2_LAGRANGE
        order = FIRST
        block = 'Fuel'
        #initial_condition = 2.01E+03
    []
    [HTC]
        family = MONOMIAL
        order = CONSTANT
        initial_condition = '${ht_coeff}'
    []
    [bison_pow_lin]
        family = L2_LAGRANGE
        order = FIRST
        block = 'Fuel'
    []
    [bison_dummy_lin]
        family = L2_LAGRANGE
        order = FIRST
        block = 'Fuel'
    []
[]

[AuxKernels]
    # [norm_Tfuel]
    #     type = NormalizationAux
    #     variable = bison_Tfuel
    #     source_variable = bison_temp
    #     normal_factor = 1
    #     execute_on = 'timestep_end' #check
    # []
    # [norm_Tref]
    #     type = NormalizationAux
    #     variable = bison_Tref
    #     source_variable = bison_temp
    #     normal_factor = 1
    #     execute_on = 'timestep_begin' #check
    # []
    [aux_power]
        type = FunctionAux
        function = axial_heat_rate
        variable = bison_power_density
        block = 'Fuel'
        execute_on = 'timestep_begin'
    []
    [norm_power_density]
        type = NormalizationAux
        variable = bison_norm_power_density
        source_variable = bison_power_density
        normal_factor = 1#1.2658064684
        execute_on = 'timestep_begin' #check
    []  
    [make_powdens_linear]
        type = NormalizationAux
        variable = bison_pow_lin
        source_variable = bison_power_density
        normal_factor = 0.000155981906099#1.2658064684
        execute_on = 'timestep_end' #check
    []    
    # [pow_dens_to_node]
    #     type = ProjectionAux
    #     v = pow_lin
    #     variable = dummy_lin
    #     execute_on = 'timestep_end'
    # []
[]

# ==============================================================================
# MULTIAPPS AND TRANSFERS
# ==============================================================================
[MultiApps]
    [sc]
      type = FullSolveMultiApp
      #app_type = ThermalHydraulicsApp
      input_files = '/home/garcsamu/Serpent/SNAP-REACTORS-PRIVATE/snapReactors/reactor_models/Wet_Experiment_Models/standard_conditions/sc_test/pow_den/sc_core.i'
      execute_on =  timestep_end
      bounding_box_padding = '0.1 0.1 0'
      positions = '0 0 0'
      output_in_position = true
    []
[]
  
#   [Transfers]
#     [flux_to_SC]
#         type = MultiAppGeneralFieldNearestLocationTransfer
#         to_multi_app = sc
#         source_variable = bison_pow_lin
#         # variable = q_dens
#         # variable = q_prime_element
#         variable = q_prime
#         from_blocks = 'Fuel'
#         to_blocks = fuel_pins
#         greedy_search = true
#     []
#   []

# ==============================================================================
# INITIAL CONDITIONS AND FUNCTIONS
# ==============================================================================
[Functions]
    # [tk_f]
    #     type = ParsedFunction
    #     symbol_names = 'bison_temp'
    #     symbol_values = 'temp_av'
    #     expression = '27.73992+bison_temp*0.027428444'#Models/SNAP10A_dimensions
    # []
    # [cp_f]
    #     type = ParsedFunction
    #     symbol_names = 'bison_temp'
    #     symbol_values = 'temp_av'
    #     expression = '472.27104+bison_temp*0.7275728'#Models/SNAP10A_dimensions
    # []
    [axial_heat_rate]
        type = ParsedFunction
        expression = '55760474.8606'#"if(z > l1 & z < l2, 55760474.8606, 0)"
        symbol_names = 'l1 l2'
        symbol_values = '0 35.56'
    []
[]

# ==============================================================================
# FLUID PROPERTIES, MATERIALS, AND USER OBJECTS
# ==============================================================================

[Materials]
    [Reflector]
      type = ADHeatConductionMaterial
      thermal_conductivity = 216.0
      specific_heat = 1925.0
      temp = bison_temp
      block = 'Reflector'
    []
    [Fuel]
      type = ADHeatConductionMaterial
      thermal_conductivity = 38.31
      specific_heat = 763.3
      temp = bison_temp
      block = 'Fuel'
    []
    [Ceramic]
      type = ADHeatConductionMaterial
      thermal_conductivity = 1.7307
      specific_heat = 837.36
      temp = bison_temp
      block = 'Ceramic'
    []
    [Clad]
      type = ADHeatConductionMaterial
      thermal_conductivity = 20.3
      specific_heat = 565.0
      temp = bison_temp
      block = 'Clad'
    []
    [GAPHE]
      type = ADHeatConductionMaterial
      thermal_conductivity = 0.3271
      specific_heat = 5193.0
      temp = bison_temp
      block = 'GAPHE'
    []
    [barrel]
      type = ADHeatConductionMaterial
      thermal_conductivity = 15.0
      specific_heat = 460.0
      temp = bison_temp
      block = 'barrel'
    []
  []

# ==============================================================================
# BOUNDARY CONDITIONS
# ==============================================================================
[BCs]
    # Convective BC outer surface fuel pin
    [convective_boundary]
        type = CoupledConvectiveHeatFluxBC
        variable = bison_temp
        boundary = 'fluid_solid_interface'
        T_infinity = bison_T_inf
        htc = HTC
    []
    # Convective BC outer surface fuel pin
[]
# ==============================================================================
# EXECUTION PARAMETERS
# ==============================================================================
[Executioner]
    type = Steady
    nl_rel_tol = 1e-8
    nl_abs_tol = 1e-8
    nl_abs_step_tol = 1e-8
    l_tol = 1e-8
    solve_type = NEWTON
    petsc_options_iname = '-pc_type -pc_hypre_type'
    petsc_options_value = 'hypre boomeramg'
[]
[VectorPostprocessors]
    [pow_dens]
        type = LineValueSampler
        start_point = '0 0 0.05'
        end_point = '0 0 0.355'
        num_points = 10
        variable = bison_pow_lin
        sort_by = 'z'
    []
[]

[Outputs]
    [exodus]
        type = Exodus
        execute_on = 'timestep_end'
    []
    [csv]
        type = CSV
        execute_on = 'initial timestep_end'
    []
    # [nemesis]
    #     type = Nemesis
    # []
[]

SC Model

# Following Advanced Burner Test Reactor Preconceptual Design Report
# Vailable at: https://www.ne.anl.gov/eda/ABTR_1cv2_ws.pdf
###################################################
# Thermal-hydraulics parameters
###################################################
T_in = 866.0 
P_out = 253727.1   # Pa
reactor_power = 671337.24 #WTh
#fuel_assemblies_per_power_unit = '${fparse 1}'
#fuel_pins_per_assembly = 211
#pin_power = '${fparse reactor_power/(fuel_assemblies_per_power_unit*fuel_pins_per_assembly)}' # Approx.
mass_flow = '${fparse 6.15}' # kg/(s)

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

# units are cm - do not forget to convert to meter
scale_factor = 0.01
#duct_thickness = '${fparse 0.3*scale_factor}'
fuel_pin_pitch = '${fparse 1.4478*scale_factor}'
fuel_pin_diameter = '${fparse 1.4268*scale_factor}'
wire_z_spacing = '${fparse 0*scale_factor}'
wire_diameter = '${fparse 0*scale_factor}'
n_rings = 9
#length_entry_fuel = '${fparse 0*scale_factor}'
length_heated_fuel = '${fparse 35.56*scale_factor}'
#length_outlet_fuel = '${fparse 0*scale_factor}'
#height = '${fparse length_entry_fuel+length_heated_fuel+length_outlet_fuel}'
#orifice_plate_height = '${fparse 0*scale_factor}'
duct_inside = '${fparse 11.43*2*scale_factor}'

# entry1 = '${fparse 0.79502/100}'
# entry2 = '${fparse 0.9652/100}'
# entry3 = '${fparse 2.1717/100}'
entry1 = '${fparse 0/100}'
entry2 = '${fparse 0/100}'
entry3 = '${fparse 0/100}'
entry_length = '${fparse entry1 + entry2 + entry3}'
# exit1 = '${fparse 2.9083/100}'
# exit2 = '${fparse 0.2286/100}'
# exit3 = '${fparse 0.87376/100}'
exit1 = '${fparse 0/100}'
exit2 = '${fparse 0/100}'
exit3 = '${fparse 0/100}'
exit_length = '${fparse exit1 + exit2 + exit3}'#'${fparse exit2 + exit3}'#
###################################################

[TriSubChannelMesh]
  [subchannel]
    type = SCMTriSubChannelMeshGenerator
    nrings = '${fparse n_rings}'
    n_cells = 10
    flat_to_flat = '${fparse duct_inside}'
    unheated_length_entry = '${fparse entry_length}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse exit_length}'
    pin_diameter = '${fparse fuel_pin_diameter}'
    pitch = '${fparse fuel_pin_pitch}'
    dwire = '${fparse wire_diameter}'
    hwire = '${fparse wire_z_spacing}'
    spacer_z = '0'
    spacer_k = '0'
  []

  [fuel_pins]
    type = SCMTriPinMeshGenerator
    input = subchannel
    nrings = '${fparse n_rings}'
    n_cells = 10
    unheated_length_entry = '${fparse entry_length}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse exit_length}'
    pitch = '${fparse fuel_pin_pitch}'
  []

  [duct]
    type = SCMTriDuctMeshGenerator
    input = fuel_pins
    nrings = '${fparse n_rings}'
    n_cells = 10
    flat_to_flat = '${fparse duct_inside}'
    unheated_length_entry = '${fparse entry_length}'
    heated_length = '${fparse length_heated_fuel}'
    unheated_length_exit = '${fparse exit_length}'
    pitch = '${fparse fuel_pin_pitch}'
  []
[]


[Functions]
  [axial_heat_rate]
    type = ParsedFunction
    expression = '1.0'
    symbol_names = 'l1 l2'
    symbol_values = '${entry_length} ${fparse length_heated_fuel}'
  []
[]

[AuxVariables]
  [mdot]
    block = subchannel
  []
  [SumWij]
    block = subchannel
  []
  [P]
    block = subchannel
  []
  [DP]
    block = subchannel
  []
  [h]
    block = subchannel
  []
  [T]
    block = subchannel
  []
  [Tpin]
    block = fuel_pins
  []
  [Dpin]
    block = fuel_pins
  []
  [rho]
    block = subchannel
  []
  [S]
    block = subchannel
  []
  [w_perim]
    block = subchannel
  []
  [q_prime]
    block = fuel_pins
    #initial_condition = 7000
  []
  [mu]
    block = subchannel
  []
  [q_prime_duct]
    block = duct
    initial_condition = 0
  []
  [Tduct]
    block = duct
  []
  [displacement]
    block = subchannel
    initial_condition = 0
  []
[]

[FluidProperties]
  [sodium]
    type = SimpleFluidProperties
    molar_mass = 0.0355
    cp = 873.0
    cv = 873.0
    specific_entropy = 1055
    viscosity = 0.0001582
    thermal_conductivity = 25.9
    thermal_expansion = 2.77e-4
  []
[]

[Problem]
  type = TriSubChannel1PhaseProblem
  fp = sodium
  n_blocks = 1
  P_out = ${P_out}
  CT = 1.0
  compute_density = false #true
  compute_viscosity = true #true
  compute_power = true #true
  P_tol = 1.0e-2
  T_tol = 1.0e-2
  implicit = true
  segregated = false
  staggered_pressure = false
  monolithic_thermal = false
  verbose_multiapps = true
  verbose_subchannel = true
  # type = NoSolveProblem
[]


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

  [w_perim_IC]
    type = SCMTriWettedPerimIC
    variable = w_perim
  []
  
  [q_prime_IC]
    type = SCMTriPowerIC
    variable = q_prime
    power = ${reactor_power} # W
    filename = "/home/garcsamu/Serpent/SNAP-REACTORS-PRIVATE/snapReactors/reactor_models/Wet_Experiment_Models/standard_conditions/sc_test/sc_standalone/S8ER_pin.txt"
    axial_heat_rate = axial_heat_rate
  []

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

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

  [P_ic]
    type = ConstantIC
    variable = P
    value = 0.0
  []

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

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

  [T_duct_ic]
    type = ConstantIC
    variable = Tduct
    value = ${T_in}
  []
  [mdot_ic]
    type = ConstantIC
    variable = mdot
    value = 0.0
  []
[]

[AuxKernels]
  [T_in_bc]
    type = ConstantAux
    variable = T
    boundary = inlet
    value = ${T_in}
    execute_on = 'timestep_begin'
    block = subchannel
  []
  [mdot_in_bc]
    type = SCMFlatMassFlowRateAux
    variable = mdot
    boundary = inlet
    #area = S
    mass_flow = ${mass_flow}
    execute_on = 'timestep_begin'
    block = subchannel
  []
[]

[UserObjects]
  active = ''
  [Tpin_avg_uo]
    type = NearestPointLayeredAverage
    direction = z
    num_layers = 10
    variable = Tpin
    block = fuel_pins
    points = '${fparse 0.012} 0.0 0.0'
    execute_on = timestep_end
  []
[]

[VectorPostprocessors]
  [pow_dens]
      type = LineValueSampler
      start_point = '0 0 0.05'
      end_point = '0 0 0.355'
      num_points = 10
      variable = q_prime
      sort_by = 'z'
      execute_on = 'initial timestep_begin'
  []
[]
[Outputs]
  [exodus]
    type = Exodus
    execute_on = 'timestep_end'
[]
[csv]
    type = CSV
    execute_on = 'initial timestep_end'
[]
[]

[Executioner]
  type = Steady
  petsc_options_value = 'hypre boomeramg'
  petsc_options_iname = '-pc_type -pc_hypre_type'
  # steady_state_detection = true
  # steady_state_tolerance = 1e-4

[]

################################################################################
# A multiapp that projects data to a detailed mesh
################################################################################
[MultiApps]
  #active = ''
  [viz]
    type = FullSolveMultiApp
    input_files = "/home/garcsamu/Serpent/SNAP-REACTORS-PRIVATE/snapReactors/reactor_models/Wet_Experiment_Models/standard_conditions/sc_test/pow_den/sc_core_viz.i"
    execute_on = "final"
  []
[]

[Transfers]
  #active = ''
  [subchannel_transfer]
    type = SCMSolutionTransfer
    to_multi_app = viz
    variable = 'mdot SumWij P DP h T rho mu S'
  []
  [pin_transfer]
    type = SCMPinSolutionTransfer
    to_multi_app = viz
    variable = 'Tpin q_prime'
  []
[]

[GiudGiud]

@kyriv1980

sc0: MOOSE Problem : Calculation of negative Enthalpy h_out = : -5.14372e+06 Axial Level= : 1

We should let the solver survive this in hopes it recovers. Depending on the fluid properties and the enthalpy value at the reference temperature, negative does not always mean completely wrong

[kyriv1980]

some initial comments:

1. You don't have to supply a power profile to the IC q_prime_IC if it's gonna be a uniform profile. That is the default value.

2. You don't have to supply a q_prime_IC IC at all, since you are transferring directly to q_prime. Remember the ICs aim is to calculate initial q_prime given total power and power shape (radial, axial).

3)Your solution tolerances (1.0e-2) are quite high. Try lowering them to 1e-4 or less.

4. Try using monolithic_thermal = true.

Errors such as the one you are getting usually happen when we have very high power and relatively low mass-flow rate.
If you say the stand alone SC runs with specific values for power and mass flow... and that those same values get transferred from BISON that means that they are not transferred correctly. I suspect that the IC you are using are not the same with the q_prime from BISON.

Regarding Guillaume's suggestions. You have access to the source code, you can track these error messages and deactivate them, but I suspect that you will get some seg fault error when the solver tries to calculate FPs with negative enthalpy.

[garcs2]

Thank you for the initial suggestions I will attempt to lower them and run again. My main point of confusion is that I'm running a BISON parent app with the working standalone SC input as the child app, without any transfers. So there is no q_prime being transferred from BISON to SC, and SC already has everything it needs to solve on its own; I'm effectively just running a working BISON followed by a working SC input.

The issue I'm trying to point out is that the SC input works when it's standalone, but fails when as a child in MultiApp when there's no transfers occuring. I'm effectively removing the transfers to help debug the issue, so now it can't be a transfers issue since there are no transfers. I'm not sure if this is more clear on what the problem that I'm experiencing is but please let me know if more clarification would be helpful!

[kyriv1980]

Can you check that the q_prime is what you would expect?

[garcs2]

Yes so q_prime is 8700 which is what I would expect for this reactor, but m_dot is showing up to be zero after the first node when I run the MultiApp with no transfers. When I run the exact same SC input but as a standalone, q_prime is 8700 and m_dot is 0.014 (roughly 6.15 kg/s divided by 438 subchannels) at every point in the subchannels. What I am concluding from this is that there is something happening with the execution timing of the mass flow rate initialization that is specifically an issue when SC is the child app of a MultiApp, see below for images.

I am able to circumvent this by setting the initial conditions of the mass flow via the IC block (note that I need to comment out the boundary), and commenting the mass flow aux kernel SCMFlatMassFlowRateAux:

  [mdot_ic]
    type = ConstantIC
    variable = mdot
    value = ${fparse mass_flow*10000}
    block = subchannel
    #boundary = inlet
  []

However, this still doesn't solve and leaves me with the following error:

sc0: *** ERROR ***
sc0: /home/garcsamu/Serpent/SNAP-REACTORS-PRIVATE/snapReactors/reactor_models/Wet_Experiment_Models/standard_conditions/sc_test/pow_den/sc_core.i:168.1:
sc0: The following error occurred in the Problem 'MOOSE Problem' of type TriSubChannel1PhaseProblem.
sc0: 
sc0: MOOSE Problem : Calculation of negative Enthalpy h_out = : -12.4469 Axial Level= : 1

So in summary:

1. When I run a standalone SC input, it runs fine
2. When I run a MultiApp with BISON as the parent, and the same standalone SC as the child, with no transfers m_dot is not populated
3. When I populate m_dot manually, and use an unphysically high mass flow rate, I get a lower negative enthalpy but still no solution

From what I understand m_dot populates at each axial level as it solves, but there is nothing different between the standalone and coupled MultiApp input at this point and there's no transfers occuring. It seems there's an issue with the Initial Condition initialization when MultiApps is used.

[kyriv1980]

"m_dot is showing up to be zero after the first node when I run the MultiApp with no transfers": Thats correct. You set the boundary condition of massflow correctly and the initial massflow to zero. What's not good is that the solver doesn't solve so these values don't get updated by the solution. So when you are looking at the exodus file you see the initial values.

Maybe you can compare the two versions of the problem (multiapp, standalone) to see if the BCs an ICs are the same?
Since the standalone solves without issue the multiapp should solve without issue as long as the BCs and ICs and power is the same.

[GiudGiud]

do you have the S8ER_pin.txt file?

[kyriv1980]

BTW. I run the standalone SCM input with a uniform radial profile. It runs fine.

[garcs2]

The S8ER_pin.txt just provides a uniform radial profile as @kyriv1980 has done in the comment above this!

Is there a working example of SCM being used in MultiApps as a child? I'd be more than happy to look through examples but didn't see any when I checked, the only MultiApps example I see is to viz which transfers on 'Final'

[kyriv1980]

This one is with SCM as main app:

https://mooseframework.inl.gov/modules/subchannel/examples/quad_thermo_mech.html

We have more cases were SCM is the sub app. They are going to be uploaded in the future.

[GiudGiud]

interestingly, this passes with the multiapp executed on timestep_begin

[GiudGiud]

ok so the petsc options are contaminating from the main app to the subapp. We'll fix this very soon, I'll have the patch ready this morning. For now you can only solve subchannel problems either:

• before any other solve
• OR using the PC_JACOBI in any other solve

[kyriv1980]

Thank you Guillaume.

[garcs2]

Thank you @kyriv1980 and @GiudGiud!