T_gas_in = 800.0 # K
m_dot_in = 2.0e-2 # kg/s
press_gas = 1.0e6 # Pa
press_coolant = 1.01325e5 # Pa

up_pipe_dia = '${units 30.0 mm -> m}'
up_pipe_length = 2.0 # m
up_pipe_n_elems =50

fir_heater_length = 2.4 # m
fir_heater_n_elems = 50
fir_heater_dia = '${units 20.0 mm -> m}'
fir_shell_dia = 0.45 # m
fir_rods_num = 10
fir_heater_A = '${fparse pi * (fir_shell_dia * fir_shell_dia - fir_rods_num * fir_heater_dia * fir_heater_dia) / 4. }'
fir_heater_dh = '${fparse 4 * fir_heater_A / (fir_rods_num * fir_heater_dia * pi + fir_shell_dia * pi)}'


# cooling section
col_length1 = 2.0 # m horizontal
col_length2 = 0.2 # m vertical
col_connect = 0.4 # m
col_n_elems_v = 10 # vertical direction elements number
col_n_elems = 50  # horizontal direction elements number
col_gas_dia = '${units 20.0 mm -> m}'
col_wall_thickness = '${units 2.0 mm -> m}'
col_coolant_dia = '${units 50.0 mm -> m}'
col_radius_wall = '${fparse col_gas_dia / 2. + col_wall_thickness}'
m_dot_col_in = 0.66 # kg/s

[GlobalParams]
  initial_p = ${press_gas}
  initial_vel = 3
  initial_T = ${T_gas_in}
  initial_vel_x = 0.02
  initial_vel_y = 0.0
  initial_vel_z = 0.02
  gravity_vector = '0 0 0' 

  rdg_slope_reconstruction = minmod
  scaling_factor_1phase = '1 1e-2 1e-4'
  scaling_factor_rhoV = 1
  scaling_factor_rhouV = 1e-2
  scaling_factor_rhovV = 1e-2
  scaling_factor_rhowV = 1e-2
  scaling_factor_rhoEV = 1e-4
  closures = thm_closures
[]

[FluidProperties]
  [ar]
    type = IdealGasFluidProperties
    molar_mass = 0.039948 # kg/mol
    gamma = 1.669 
    k = 0.04071
    mu = 4.7872e-5
  []

  [water]
    type = StiffenedGasFluidProperties
    gamma = 2.35
    cv = 1816.0
    q = -1.167e6
    p_inf = 1.0e9
    q_prime = 0
  []
[]

[Closures]
  [thm_closures]
    type = Closures1PhaseTHM
  []
  [none_closures]
    type = Closures1PhaseNone
  []
[]

[Materials]
  [Re_mat]
    type = ADReynoldsNumberMaterial
    Re = Re
    rho = rho
    vel = vel
    D_h = D_h
    mu = mu
    block = 'col_sys/gcol_chan1 col_sys/gcol_chan2 col_sys/gcol_chan3 col_sys/gcol_chan4 col_sys/gcol_chan5'
  []

  [f_mat]
    type = ADParsedMaterial
    property_name = f_D
    constant_names = 'a b c'
    constant_expressions = '1 0.1 -0.5'
    material_property_names = 'Re'
    expression = 'a + b * Re^c'
    block = 'col_sys/gcol_chan1 col_sys/gcol_chan2 col_sys/gcol_chan3 col_sys/gcol_chan4 col_sys/gcol_chan5'
  []

  [Pr_mat]
    type = ADPrandtlNumberMaterial
    Pr = Pr
    cp = cp
    mu = mu
    k = k
    block = 'col_sys/gcol_chan1 col_sys/gcol_chan2 col_sys/gcol_chan3 col_sys/gcol_chan4 col_sys/gcol_chan5'
  []

  [Nu_mat]
    type = ADParsedMaterial
    property_name = 'Nu'
    constant_names = 'a b c'
    constant_expressions = '0.03 0.9 0.5'
    material_property_names = 'Re Pr'
    expression = 'a * Re ^b * Pr^c'
    block = 'col_sys/gcol_chan1 col_sys/gcol_chan2 col_sys/gcol_chan3 col_sys/gcol_chan4 col_sys/gcol_chan5'
  []

  [Hw_mat]
    type = ADConvectiveHeatTransferCoefficientMaterial
    D_h = D_h
    k = k
    Nu = Nu
    Hw = Hw
    block = 'col_sys/gcol_chan1 col_sys/gcol_chan2 col_sys/gcol_chan3 col_sys/gcol_chan4 col_sys/gcol_chan5'
  []

[]

[SolidProperties]
  [steel]
    type = ThermalFunctionSolidProperties
    rho = 8050.
    k = 45.
    cp = 466.
  []
[]

[Components]

  [first_heater]

    [fir_ht_chan]
      type = FlowChannel1Phase
      position = '-${col_connect} 0 ${fparse up_pipe_length - col_length2 * 2}'
      orientation = '1 0 0'
      length = ${fir_heater_length}
      n_elems = ${fir_heater_n_elems}
      roughness = 6.3e-6
      A = '${ fir_heater_A }'
      D_h = ${ fir_heater_dh }
      initial_p = 1.0e6
      initial_T = 500.0
      initial_vel = 3.0
      fp = ar
      # closures = none_closures
    []
    
    [fir_heating]
      type = HeatSourceVolumetric1Phase
      flow_channel = first_heater/fir_ht_chan
      q = 12000.0 
    []

  []

  [jct1]
    type = JunctionOneToOne1Phase
    connections = 'first_heater/fir_ht_chan:out bypass_pipe:in'
  []

  [bypass_pipe]
    type = FlowChannel1Phase
    position = '${fparse fir_heater_length - col_connect} 0 ${fparse up_pipe_length - col_length2 * 2}'
    orientation = '0 0 -1'
    length = ${fparse up_pipe_length - col_length2 * 2}
    n_elems = ${fir_heater_n_elems}
    roughness = 6.3e-6
    A = '${fparse pi * col_gas_dia * col_gas_dia / 4.}'
    D_h = ${col_gas_dia}
    initial_p = 1.0e6
    initial_T = 600.0
    initial_vel = 3.0
    fp = ar
  []

  [jct2]
    type = JunctionOneToOne1Phase
    connections = 'bypass_pipe:out col_sys/gcol_chan1:in'
  []
  
  [col_sys]
 
    [gcol_chan1]
      type = FlowChannel1Phase
      position = '${col_length1} 0 0'
      orientation = '-1 0 0'
      length = ${col_length1}
      n_elems = ${col_n_elems}
      roughness = 6.3e-6
      A = '${fparse pi * col_gas_dia * col_gas_dia / 4.}'
      D_h = ${col_gas_dia}
      initial_p = 1.0e6
      initial_T = 600.0
      initial_vel = 3.0
      fp = ar
      closures = none_closures
    []

    [gct1]
      type = JunctionOneToOne1Phase
      connections = 'col_sys/gcol_chan1:out col_sys/gcol_chan2:in'
    []

    [gcol_chan2]
      type = FlowChannel1Phase
      position = '0 0 0'
      orientation = '0 0 -1'
      length = ${col_length2}
      n_elems = ${col_n_elems_v}
      roughness = 6.3e-6
      A = '${fparse pi * col_gas_dia * col_gas_dia / 4.}'
      D_h = ${col_gas_dia}
      initial_p = 1.0e6
      initial_T = 500.0
      initial_vel = 3.0
      fp = ar
      closures = none_closures
    []

    [gct2]
      type = JunctionOneToOne1Phase
      connections = 'col_sys/gcol_chan2:out col_sys/gcol_chan3:in'
    []

    [gcol_chan3]
      type = FlowChannel1Phase
      position = '0 0 -${col_length2}'
      orientation = '1 0 0'
      length = ${col_length1}
      n_elems = ${col_n_elems}
      roughness = 6.3e-6
      A = '${fparse pi * col_gas_dia * col_gas_dia / 4.}'
      D_h = ${col_gas_dia}
      initial_p = 1.0e6
      initial_T = 400.0
      initial_vel = 3.0
      fp = ar
      closures = none_closures
    []

    [gct3]
      type = JunctionOneToOne1Phase
      connections = 'col_sys/gcol_chan3:out col_sys/gcol_chan4:in'
    []

    [gcol_chan4]
      type = FlowChannel1Phase
      position = '${col_length1} 0 -${col_length2}'
      orientation = '0 0 -1'
      length = ${col_length2}
      n_elems = ${col_n_elems_v}
      roughness = 6.3e-6
      A = '${fparse pi * col_gas_dia * col_gas_dia / 4.}'
      D_h = ${col_gas_dia}
      initial_p = 1.0e6
      initial_T = 400.0
      initial_vel = 3.0
      fp = ar
      closures = none_closures
    []

    [gct4]
      type = JunctionOneToOne1Phase
      connections = 'col_sys/gcol_chan4:out col_sys/gcol_chan5:in'
    []

    [gcol_chan5]
      type = FlowChannel1Phase
      position = '${col_length1} 0 -${fparse 2 * col_length2}'
      orientation = '-1 0 0'
      length = ${col_length1}
      n_elems = ${col_n_elems}
      roughness = 6.3e-6
      A = '${fparse pi * col_gas_dia * col_gas_dia / 4.}'
      D_h = ${col_gas_dia}
      initial_p = 1.0e6
      initial_T = 400.0
      initial_vel = 3.0
      fp = ar
      closures = none_closures
    []

    [gct5]
      type = JunctionOneToOne1Phase
      connections = 'col_sys/gcol_chan5:out col_sys/gcol_connect:in'
    []

    [gcol_connect]
      type = FlowChannel1Phase
      position = '0 0 -${fparse 2 * col_length2}'
      orientation = '-1 0 0'
      length = ${col_connect}
      n_elems = ${col_n_elems_v}
      roughness = 6.3e-6
      A = '${fparse pi * col_gas_dia * col_gas_dia / 4.}'
      D_h = ${col_gas_dia}
      initial_p = 1.0e6
      initial_T = 320.0
      initial_vel = 3.0
      fp = ar
    []

    [col_gas1]
      type = HeatTransferFromHeatStructure1Phase
      hs = col_sys/wall1
      hs_side = inner
      flow_channel = col_sys/gcol_chan1
      P_hf = '${fparse pi * col_gas_dia}'
    []

    [col_coolant1]
      type = HeatTransferFromHeatStructure1Phase
      hs = col_sys/wall1
      hs_side = outer
      flow_channel = col_sys/wcol_chan5
      P_hf = '${fparse 2 * pi * col_radius_wall}'
    []

    [col_gas2]
      type = HeatTransferFromHeatStructure1Phase
      hs = col_sys/wall2
      hs_side = inner
      flow_channel = col_sys/gcol_chan2
      P_hf = '${fparse pi * col_gas_dia}'
    []

    [col_coolant2]
      type = HeatTransferFromHeatStructure1Phase
      hs = col_sys/wall2
      hs_side = outer
      flow_channel = col_sys/wcol_chan4
      P_hf = '${fparse 2 * pi * col_radius_wall}'
    []

    [col_gas3]
      type = HeatTransferFromHeatStructure1Phase
      hs = col_sys/wall3
      hs_side = inner
      flow_channel = col_sys/gcol_chan3
      P_hf = '${fparse pi * col_gas_dia}'
    []

    [col_coolant3]
      type = HeatTransferFromHeatStructure1Phase
      hs = col_sys/wall3
      hs_side = outer
      flow_channel = col_sys/wcol_chan3
      P_hf = '${fparse 2 * pi * col_radius_wall}'
    []

    [col_gas4]
      type = HeatTransferFromHeatStructure1Phase
      hs = col_sys/wall4
      hs_side = inner
      flow_channel = col_sys/gcol_chan4
      P_hf = '${fparse pi * col_gas_dia}'
    []

    [col_coolant4]
      type = HeatTransferFromHeatStructure1Phase
      hs = col_sys/wall4
      hs_side = outer
      flow_channel = col_sys/wcol_chan2
      P_hf = '${fparse 2 * pi * col_radius_wall}'
    []

    [col_gas5]
      type = HeatTransferFromHeatStructure1Phase
      hs = col_sys/wall5
      hs_side = inner
      flow_channel = col_sys/gcol_chan5
      P_hf = '${fparse pi * col_gas_dia}'
    []

    [col_coolant5]
      type = HeatTransferFromHeatStructure1Phase
      hs = col_sys/wall5
      hs_side = outer
      flow_channel = col_sys/wcol_chan1
      P_hf = '${fparse 2 * pi * col_radius_wall}'
    []

    [wall1]
      type = HeatStructureCylindrical
      position = '${col_length1} 0 0'
      orientation = '-1 0 0'
      length = ${col_length1}
      n_elems = ${col_n_elems}
      widths = '${col_wall_thickness}'
      n_part_elems = '3'
      solid_properties = 'steel'
      solid_properties_T_ref = '300'
      names = '1'
      inner_radius = '${fparse col_gas_dia / 2.}'
    []

    [wall2]
      type = HeatStructureCylindrical
      position = '0 0 0'
      orientation = '0 0 -1'
      length = ${col_length2}
      n_elems = ${col_n_elems_v}
      widths = '${col_wall_thickness}'
      n_part_elems = '3'
      solid_properties = 'steel'
      solid_properties_T_ref = '300'
      names = '2'
      inner_radius = '${fparse col_gas_dia / 2.}'
    []

    [wall3]
      type = HeatStructureCylindrical
      position = '0 0 -${col_length2}'
      orientation = '1 0 0'
      length = ${col_length1}
      n_elems = ${col_n_elems}
      widths = '${col_wall_thickness}'
      n_part_elems = '3'
      solid_properties = 'steel'
      solid_properties_T_ref = '300'
      names = '3'
      inner_radius = '${fparse col_gas_dia / 2.}'
    []

    [wall4]
      type = HeatStructureCylindrical
      position = '${col_length1} 0 -${fparse col_length2}'
      orientation = '0 0 -1'
      length = ${col_length2}
      n_elems = ${col_n_elems_v}
      widths = '${col_wall_thickness}'
      n_part_elems = '3'
      solid_properties = 'steel'
      solid_properties_T_ref = '300'
      names = '4'
      inner_radius = '${fparse col_gas_dia / 2.}'
    []

    [wall5]
      type = HeatStructureCylindrical
      position = '${col_length1} 0 -${fparse 2 * col_length2}'
      orientation = '-1 0 0'
      length = ${col_length1}
      n_elems = ${col_n_elems}
      widths = '${col_wall_thickness}'
      n_part_elems = '3'
      solid_properties = 'steel'
      solid_properties_T_ref = '300'
      names = '5'
      inner_radius = '${fparse col_gas_dia / 2.}'
    []

    [col_coolant_inlet]
      type = InletMassFlowRateTemperature1Phase
      input = 'col_sys/wcol_chan1:in'
      m_dot = ${m_dot_col_in}
      T = 300.0
    []

    [wcol_chan1]
      type = FlowChannel1Phase
      position = '0 0 -${fparse 2 * col_length2}'
      orientation = '1 0 0'
      length = ${col_length1}
      n_elems = ${col_n_elems}
      A = '${fparse pi * (col_coolant_dia * col_coolant_dia / 4. - col_radius_wall * col_radius_wall)}'
      D_h = '${fparse col_coolant_dia - (2 * col_radius_wall)}'
      fp = water
      initial_T = 300.0
    []

    [wct1]
      type = JunctionOneToOne1Phase
      connections = 'col_sys/wcol_chan1:out col_sys/wcol_chan2:in'
    []

    [wcol_chan2]
      type = FlowChannel1Phase
      position = '${col_length1} 0 -${fparse 2 * col_length2}'
      orientation = '0 0 1'
      length = ${col_length2}
      n_elems = ${col_n_elems_v}
      A = '${fparse pi * (col_coolant_dia * col_coolant_dia / 4. - col_radius_wall * col_radius_wall)}'
      D_h = '${fparse col_coolant_dia - (2 * col_radius_wall)}'
      fp = water
      initial_T = 300.0
    []

    [wct2]
      type = JunctionOneToOne1Phase
      connections = 'col_sys/wcol_chan2:out col_sys/wcol_chan3:in'
    []

    [wcol_chan3]
      type = FlowChannel1Phase
      position = '${col_length1} 0 -${col_length2}'
      orientation = '-1 0 0'
      length = ${col_length1}
      n_elems = ${col_n_elems}
      A = '${fparse pi * (col_coolant_dia * col_coolant_dia / 4. - col_radius_wall * col_radius_wall)}'
      D_h = '${fparse col_coolant_dia - (2 * col_radius_wall)}'
      fp = water
      initial_T = 300.0
    []

    [wct3]
      type = JunctionOneToOne1Phase
      connections = 'col_sys/wcol_chan3:out col_sys/wcol_chan4:in'
    []

    [wcol_chan4]
      type = FlowChannel1Phase
      position = '0 0 -${col_length2}'
      orientation = '0 0 1'
      length = ${col_length2}
      n_elems = ${col_n_elems_v}
      A = '${fparse pi * (col_coolant_dia * col_coolant_dia / 4. - col_radius_wall * col_radius_wall)}'
      D_h = '${fparse col_coolant_dia - (2 * col_radius_wall)}'
      fp = water
      initial_T = 300.0
    []

    [wct4]
      type = JunctionOneToOne1Phase
      connections = 'col_sys/wcol_chan4:out col_sys/wcol_chan5:in'
    []

    [wcol_chan5]
      type = FlowChannel1Phase
      position = '0 0 0'
      orientation = '1 0 0'
      length = ${col_length1}
      n_elems = ${col_n_elems}
      A = '${fparse pi * (col_coolant_dia * col_coolant_dia / 4. - col_radius_wall * col_radius_wall)}'
      D_h = '${fparse col_coolant_dia - (2 * col_radius_wall)}'
      fp = water
      initial_T = 300.0
    []
  
    [col_coolant_outlet]
      type = Outlet1Phase
      input = 'col_sys/wcol_chan5:out'
      p = ${press_coolant}    
    []
  
  []

  [pump]
    type = Pump1Phase
    position = '-${col_connect} 0 -${fparse 2 * col_length2}'
    connections = 'col_sys/gcol_connect:out up_pipe:in'
    volume = 1e-4
    A_ref = '${fparse pi * col_gas_dia * col_gas_dia / 4.}'
    head = 5.
  []

  [up_pipe]
    type = FlowChannel1Phase
    position = '-${col_connect} 0 -${fparse 2 * col_length2}'
    orientation = '0 0 1'
    length = ${up_pipe_length}
    n_elems = ${up_pipe_n_elems}
    A = '${fparse pi * up_pipe_dia * up_pipe_dia / 4.}'
    D_h = ${up_pipe_dia}
    initial_p = 1.0e6
    initial_T = 300.0
    initial_vel = 3.0
    fp = ar
  []

  [jct3]
    type = JunctionOneToOne1Phase
    connections = 'up_pipe:out first_heater/fir_ht_chan:in'
  []
  
[]

[ControlLogic]
  [set_point]
    type = GetFunctionValueControl
    function = ${m_dot_in}
  []

  [pid]
    type = PIDControl
    initial_value = 10.0
    set_point = set_point:value
    input = m_dot_pump
    K_p = 150.
    K_i = 450.
    K_d = 150.
  []

  [set_pump_head]
    type = SetComponentRealValueControl
    component = pump
    parameter = head
    value = pid:output
  []
[]

[Postprocessors]
  
  [m_dot_pump]
    type = ADFlowJunctionFlux1Phase
    boundary = col_sys/gcol_chan1:in
    connection_index = 1
    equation = mass
    junction = col_sys/gct1
  []

  [col_sys_T_1_in] # temprature of gas from pipe1_out in col_sys 
    type = SideAverageValue
    boundary = col_sys/gcol_chan1:in
    variable = T
  []
  
  [col_sys_T_1_out] # temprature of gas from pipe1_out in col_sys 
    type = SideAverageValue
    boundary = col_sys/gcol_chan1:out
    variable = T
  []

  [col_sys_T_2_out] # temprature of gas from pipe2_out in col_sys 
    type = SideAverageValue
    boundary = col_sys/gcol_chan2:out
    variable = T
  []

  [col_sys_T_3_out] # temprature of gas from pipe3_out in col_sys 
    type = SideAverageValue
    boundary = col_sys/gcol_chan3:out
    variable = T
  []

  [col_sys_T_4_out] # temprature of gas from pipe4_out in col_sys 
    type = SideAverageValue
    boundary = col_sys/gcol_chan4:out
    variable = T
  []

  [col_sys_T_5_out] # temprature of gas from pipe5_out in col_sys 
    type = SideAverageValue
    boundary = col_sys/gcol_chan5:out
    variable = T
  []

  [col_sys_colant_Tin]
    type = SideAverageValue
    boundary = col_sys/wcol_chan1:in
    variable = T
  []

  [col_sys_colant_T1_out]
    type = SideAverageValue
    boundary = col_sys/wcol_chan1:out
    variable = T
  []

  [col_sys_colant_T2_out]
    type = SideAverageValue
    boundary = col_sys/wcol_chan2:out
    variable = T
  []

  [col_sys_colant_T3_out]
    type = SideAverageValue
    boundary = col_sys/wcol_chan3:out
    variable = T
  []

  [col_sys_colant_T4_out]
    type = SideAverageValue
    boundary = col_sys/wcol_chan4:out
    variable = T
  []

  [col_sys_colant_T5out]
    type = SideAverageValue
    boundary = col_sys/wcol_chan5:out
    variable = T
  []
  
  [pump_out_temp]
    type = SideAverageValue
    boundary = up_pipe:in
    variable = T
  []
  
  [m_dot_coolant]
    type = ADFlowBoundaryFlux1Phase
    boundary = col_sys/col_coolant_outlet
    equation = mass
  []

  [m_dot_gas]
    type = ADFlowJunctionFlux1Phase
    boundary = col_sys/gcol_chan1:in
    connection_index = 1
    equation = mass
    junction = col_sys/gct1
  []

[]

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

[Executioner]
  type = Transient
  start_time = 0

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.000001
  []
  dtmax = 5
  end_time = 5000.0

  line_search = basic
  solve_type = NEWTON

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-8
  nl_max_its = 25

[]

[Outputs]
  exodus = true
  csv = true

  [console]
    type = Console
    max_rows = 1
    outlier_variable_norms = false
  []
  print_linear_residuals = false
[]