Test Problems

The following test problems can be used to assess various models in THM:

Name	File(s)	Description
Area constriction	[input]	Constricted cross-sectional area through channel
Double rarefaction	[input]	Riemann problem that has a double-rarefaction solution
Free fall	[input]	Acceleration of a fluid due to gravity
Lax shock tube	[input]	Lax shock tube test problem
MMS	[input]	Method of Manufactured Solutions convergence rate verification
Natural circlulation	[input]	Natural circulation loop
Pressure drop	[input]	Pressure drop over a channel
Sedov blast wave	[input]	Sedov blast wave test problem
Sod shock tube	[input]	Classic Sod shock tube test problem
Square wave	[input]	Square wave problem
Super sonic tube	[input]	Super sonic tube
Water hammer	[input]	Water hammer test problem
Woodward-Colella blast wave	[input]	Woodward-Colella blast wave problem

(modules/thermal_hydraulics/test/tests/problems/area_constriction/area_constriction.i)

# This test features air flowing through a channel whose cross-sectional area
# shrinks to half its value in the right half. Assuming incompressible flow
# conditions, such as having a low Mach number, the velocity should approximately
# double from inlet to outlet.

p_outlet = 1e5

[GlobalParams]
  gravity_vector = '0 0 0'

  initial_T = 300
  initial_p = ${p_outlet}
  initial_vel = initial_vel_fn

  fp = fp
  closures = simple_closures
  f = 0

  scaling_factor_1phase = '1 1 1e-5'
[]

[FluidProperties]
  [fp]
    type = IdealGasFluidProperties
    gamma = 1.4
  []
[]

[Closures]
  [simple_closures]
    type = Closures1PhaseSimple
  []
[]

[Functions]
  [A_fn]
    type = PiecewiseConstant
    axis = x
    direction = right
    x = '0.5 1.0'
    y = '1.0 0.5'
  []
  [initial_vel_fn]
    type = PiecewiseConstant
    axis = x
    direction = right
    x = '0.5 1.0'
    y = '1.0 2'
  []
[]

[Components]
  [inlet]
    type = InletDensityVelocity1Phase
    input = 'pipe:in'
    rho = 1.16263315948279 # rho @ (p = 1e5 Pa, T = 300 K)
    vel = 1
  []

  [pipe]
    type = FlowChannel1Phase
    position = '0 0 0'
    orientation = '1 0 0'
    length = 1
    n_elems = 100

    A = A_fn
  []

  [outlet]
    type = Outlet1Phase
    input = 'pipe:out'
    p = ${p_outlet}
  []
[]

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

[Executioner]
  type = Transient
  scheme = 'bdf2'

  end_time = 10
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.001
    optimal_iterations = 5
    iteration_window = 1
    growth_factor = 1.2
  []

  steady_state_detection = true

  solve_type = PJFNK
  nl_rel_tol = 1e-10
  nl_abs_tol = 1e-8
  nl_max_its = 15

  l_tol = 1e-3
  l_max_its = 10
[]

[Outputs]
  exodus = true
  velocity_as_vector = false
  show = 'A rho vel p'
[]

(modules/thermal_hydraulics/test/tests/problems/double_rarefaction/1phase.i)

# Riemann problem that has a double-rarefaction solution

[GlobalParams]
  gravity_vector = '0 0 0'

  rdg_slope_reconstruction = minmod

  closures = simple_closures
[]

[Functions]
  [vel_ic_fn]
    type = PiecewiseConstant
    axis = x
    direction = right
    x = ' 0.0 0.1'
    y = '-1.0 1.0'
  []
[]

[FluidProperties]
  [fp]
    type = IdealGasFluidProperties
    gamma = 1.4
    molar_mass = 11.64024372
  []
[]

[Closures]
  [simple_closures]
    type = Closures1PhaseSimple
  []
[]

[Components]
  [pipe]
    type = FlowChannel1Phase

    fp = fp

    # geometry
    position = '-1 0 0'
    orientation = '1 0 0'
    length = 2.0
    n_elems = 100
    A = 1.0

    # IC
    initial_T = 0.04
    initial_p = 0.2
    initial_vel = vel_ic_fn

    f = 0
  []

  [left_boundary]
    type = FreeBoundary1Phase
    input = 'pipe:in'
  []

  [right_boundary]
    type = FreeBoundary1Phase
    input = 'pipe:out'
  []
[]

[Executioner]
  type = Transient
  [TimeIntegrator]
    type = ExplicitSSPRungeKutta
    order = 2
  []
  solve_type = LINEAR

  l_tol = 1e-4

  nl_rel_tol = 1e-20
  nl_abs_tol = 1e-8
  nl_max_its = 60

  # run to t = 0.6
  start_time = 0.0
  dt = 1e-3
  num_steps = 600
  abort_on_solve_fail = true
[]

[Outputs]
  file_base = '1phase'
  velocity_as_vector = false
  execute_on = 'initial timestep_end'
  [out]
    type = Exodus
    show = 'p T vel'
  []
[]

(modules/thermal_hydraulics/test/tests/problems/freefall/freefall.i)

# Tests acceleration of a fluid due to gravity. The flow exiting the bottom
# of the flow channel enters the top, so the flow should uniformly accelerate
# at the rate of acceleration due to gravity.

acceleration = -10.0
dt = 0.1
num_steps = 5
time = ${fparse num_steps * dt}

# The expected velocity is the following:
#   u = a * t
#     = -10 * 0.5
#     = -5

[GlobalParams]
  gravity_vector = '0 0 ${acceleration}'

  initial_p = 1e5
  initial_T = 300
  initial_vel = 0

  scaling_factor_1phase = '1 1 1e-5'

  closures = simple_closures
[]

[FluidProperties]
  [fp]
    type = StiffenedGasFluidProperties
    gamma = 2.35
    cv = 1816
    q = -1.167e6
    q_prime = 0
    p_inf = 1e9
  []
[]

[Closures]
  [simple_closures]
    type = Closures1PhaseSimple
  []
[]

[Components]
  [pipe]
    type = FlowChannel1Phase
    position = '0 0 0'
    orientation = '0 0 1'
    length = 1
    n_elems = 100
    A = 1
    f = 0
    fp = fp
  []

  [junction]
    type = JunctionOneToOne1Phase
    connections = 'pipe:in pipe:out'
  []
[]

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

[Executioner]
  type = Transient

  scheme = bdf2
  end_time = ${time}
  dt = ${dt}
  num_steps = ${num_steps}
  abort_on_solve_fail = true

  solve_type = NEWTON
  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-8
  nl_max_its = 10
  l_tol = 1e-3
  l_max_its = 10

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
[]

[Postprocessors]
  [vel_avg]
    type = ElementAverageValue
    variable = 'vel'
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[Outputs]
  velocity_as_vector = false
  [out]
    type = CSV
    execute_on = 'FINAL'
  []
[]

(modules/thermal_hydraulics/test/tests/problems/lax_shock_tube/lax_shock_tube.i)

# This test problem is the Lax shock tube test problem,
# which is a Riemann problem with the following parameters:
#   * domain = (0,1)
#   * gravity = 0
#   * EoS: Ideal gas EoS with gamma = 1.4, R = 0.71428571428571428571
#   * interface: x = 0.5
#   * typical end time: 0.15
# Left initial values:
#   * rho = 0.445
#   * vel = 0.692
#   * p = 3.52874226
# Right initial values:
#   * rho = 0.5
#   * vel = 0
#   * p = 0.571

[GlobalParams]
  gravity_vector = '0 0 0'

  rdg_slope_reconstruction = minmod

  closures = simple_closures
[]

[Functions]
  [p_ic_fn]
    type = PiecewiseConstant
    axis = x
    direction = right
    x = '0.5        1.0'
    y = '3.52874226 0.571'
  []

  [T_ic_fn]
    type = PiecewiseConstant
    axis = x
    direction = right
    x = '0.5              1.0'
    y = '11.1016610426966 1.5988'
  []

  [vel_ic_fn]
    type = PiecewiseConstant
    axis = x
    direction = right
    x = '0.5   1.0'
    y = '0.692 0.0'
  []
[]

[FluidProperties]
  [fp]
    type = IdealGasFluidProperties
    gamma = 1.4
    molar_mass = 11.64024372
  []
[]

[Closures]
  [simple_closures]
    type = Closures1PhaseSimple
  []
[]

[Components]
  [pipe]
    type = FlowChannel1Phase

    fp = fp

    # geometry
    position = '0 0 0'
    orientation = '1 0 0'
    length = 1.0
    n_elems = 100
    A = 1.0

    # IC
    initial_T = T_ic_fn
    initial_p = p_ic_fn
    initial_vel = vel_ic_fn

    f = 0
  []

  [left_boundary]
    type = FreeBoundary1Phase
    input = 'pipe:in'
  []

  [right_boundary]
    type = FreeBoundary1Phase
    input = 'pipe:out'
  []
[]

[Executioner]
  type = Transient
  [TimeIntegrator]
    type = ExplicitSSPRungeKutta
    # add order via 'cli_args' in 'tests'
  []
  solve_type = LINEAR

  l_tol = 1e-4

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

  # run to t = 0.15
  start_time = 0.0
  dt = 1e-3
  num_steps = 150
  abort_on_solve_fail = true
[]

[Outputs]
  file_base = 'lax_shock_tube'
  velocity_as_vector = false
  execute_on = 'initial timestep_end'
  [out]
    type = Exodus
    show = 'rho p vel'
  []
[]

(modules/thermal_hydraulics/test/tests/problems/mms/mms_1phase.i)

# Method of manufactured solutions (MMS) problem for 1-phase flow model.
#
# The python script mms_derivation.py derives the MMS sources used in this
# input file.
#
# To perform a convergence study, run this input file with different values of
# 'refinement_level', starting with 0. Manually create a CSV file (call it the
# "convergence CSV file") to store the error vs. mesh size data. It should have
# the columns specified in the plot script plot_convergence_1phase.py. Copy the
# CSV output from each run into the convergence CSV file. After all of the runs,
# run the plot script using python.

refinement_level = 0 # 0 is initial
n_elems_coarse = 10
n_elems = ${fparse int(n_elems_coarse * 2^refinement_level)}

dt = 1e-6
t_end = ${fparse dt * 10}

area = 1.0
gamma = 2.0
M = 0.05
A = 1
B = 1
C = 1

aA = ${fparse area}

R_univ = 8.3144598
R = ${fparse R_univ / M}
cp = ${fparse gamma * R / (gamma - 1.0)}
cv = ${fparse cp / gamma}

[GlobalParams]
  gravity_vector = '0 0 0'

  closures = simple_closures
[]

[Functions]
  # solutions
  [rho_fn]
    type = ParsedFunction
    expression = 'A * (sin(B*x + C*t) + 2)'
    symbol_names = 'A B C'
    symbol_values = '${A} ${B} ${C}'
  []
  [vel_fn]
    type = ParsedFunction
    expression = 'A * t * sin(pi * x)'
    symbol_names = 'A'
    symbol_values = '${A}'
  []
  [p_fn]
    type = ParsedFunction
    expression = 'A * (cos(B*x + C*t) + 2)'
    symbol_names = 'A B C'
    symbol_values = '${A} ${B} ${C}'
  []
  [T_fn]
    type = ParsedFunction
    expression = '(cos(B*x + C*t) + 2)/(cv*(gamma - 1)*(sin(B*x + C*t) + 2))'
    symbol_names = 'B C gamma cv'
    symbol_values = '${B} ${C} ${gamma} ${cv}'
  []

  # MMS sources
  [rho_src_fn]
    type = ParsedFunction
    expression = 'A^2*B*t*sin(pi*x)*cos(B*x + C*t) + pi*A^2*t*(sin(B*x + C*t) + 2)*cos(pi*x) + A*C*cos(B*x + C*t)'
    symbol_names = 'A B C'
    symbol_values = '${A} ${B} ${C}'
  []
  [rhou_src_fn]
    type = ParsedFunction
    expression = 'A^3*B*t^2*sin(pi*x)^2*cos(B*x + C*t) + 2*pi*A^3*t^2*(sin(B*x + C*t) + 2)*sin(pi*x)*cos(pi*x) + A^2*C*t*sin(pi*x)*cos(B*x + C*t) + A^2*(sin(B*x + C*t) + 2)*sin(pi*x) - A*B*sin(B*x + C*t)'
    symbol_names = 'A B C'
    symbol_values = '${A} ${B} ${C}'
  []
  [rhoE_src_fn]
    type = ParsedFunction
    expression = 'A*C*(A^2*t^2*sin(pi*x)^2/2 + (cos(B*x + C*t) + 2)/((gamma - 1)*(sin(B*x + C*t) + 2)))*cos(B*x + C*t) + pi*A*t*(A*(A^2*t^2*sin(pi*x)^2/2 + (cos(B*x + C*t) + 2)/((gamma - 1)*(sin(B*x + C*t) + 2)))*(sin(B*x + C*t) + 2) + A*(cos(B*x + C*t) + 2))*cos(pi*x) + A*t*(A*B*(A^2*t^2*sin(pi*x)^2/2 + (cos(B*x + C*t) + 2)/((gamma - 1)*(sin(B*x + C*t) + 2)))*cos(B*x + C*t) - A*B*sin(B*x + C*t) + A*(sin(B*x + C*t) + 2)*(pi*A^2*t^2*sin(pi*x)*cos(pi*x) - B*sin(B*x + C*t)/((gamma - 1)*(sin(B*x + C*t) + 2)) - B*(cos(B*x + C*t) + 2)*cos(B*x + C*t)/((gamma - 1)*(sin(B*x + C*t) + 2)^2)))*sin(pi*x) + A*(sin(B*x + C*t) + 2)*(A^2*t*sin(pi*x)^2 - C*sin(B*x + C*t)/((gamma - 1)*(sin(B*x + C*t) + 2)) - C*(cos(B*x + C*t) + 2)*cos(B*x + C*t)/((gamma - 1)*(sin(B*x + C*t) + 2)^2))'
    symbol_names = 'A B C gamma'
    symbol_values = '${A} ${B} ${C} ${gamma}'
  []
[]

[FluidProperties]
  [fp]
    type = IdealGasFluidProperties
    gamma = ${gamma}
    molar_mass = ${M}
  []
[]

[Closures]
  [simple_closures]
    type = Closures1PhaseSimple
  []
[]

[Components]
  [pipe]
    type = FlowChannel1Phase

    fp = fp

    # geometry
    position = '0 0 0'
    orientation = '1 0 0'
    length = 1.0
    n_elems = ${n_elems}
    A = ${area}

    # IC
    initial_p = p_fn
    initial_T = T_fn
    initial_vel = 0

    f = 0
  []

  [left_boundary]
    type = InletFunction1Phase
    input = 'pipe:in'
    p = p_fn
    rho = rho_fn
    vel = vel_fn
  []

  [right_boundary]
    type = InletFunction1Phase
    input = 'pipe:out'
    p = p_fn
    rho = rho_fn
    vel = vel_fn
  []
[]

[Kernels]
  [rho_src]
    type = BodyForce
    variable = rhoA
    function = rho_src_fn
    value = ${aA}
  []
  [rhou_src]
    type = BodyForce
    variable = rhouA
    function = rhou_src_fn
    value = ${aA}
  []
  [rhoE_src]
    type = BodyForce
    variable = rhoEA
    function = rhoE_src_fn
    value = ${aA}
  []
[]

[Postprocessors]
  [rho_err]
    type = ElementL1Error
    variable = rho
    function = rho_fn
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [vel_err]
    type = ElementL1Error
    variable = vel
    function = vel_fn
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p_err]
    type = ElementL1Error
    variable = p
    function = p_fn
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[Executioner]
  type = Transient

  [TimeIntegrator]
    type = ExplicitSSPRungeKutta
    order = 3
  []

  start_time = 0
  dt = ${dt}
  end_time = ${t_end}
  abort_on_solve_fail = true

  [Quadrature]
    type = GAUSS
    order = FIRST
  []
[]

[Outputs]
  csv = true
  execute_on = 'FINAL'
  velocity_as_vector = false
[]

(modules/thermal_hydraulics/test/tests/problems/natural_circulation/natural_circulation.i)

# Natural circulation loop
#
# The setup consists of 4 connected 1-m pipes, forming a square:
#
#                  top_pipe
#              *--------------* (1,1)
#              |              |
#              | <-        <- |                | g
#  heated_pipe | <-        <- | cooled_pipe    V
#              | <-        <- |
#              |              |
#        (0,0) *--------------*
#                 bottom_pipe
#
# Heating and cooling occurs in the range z = (0.2 m, 0.8 m) with uniform heat fluxes.

n_elems = 50

diam = 0.1
length = 1.0
heated_length = 0.6
power = 1e3

p_initial = 100e3
T_ambient = 300
htc = 25.0

area = ${fparse 0.25 * pi * diam^2}
S_heated = ${fparse pi * diam * heated_length}
S_cooled = ${fparse pi * diam * heated_length}

output_variables = 'rho p T vel rhouA'

[GlobalParams]
  gravity_vector = '0 0 -9.81'
  length = ${length}
  n_elems = ${n_elems}

  A = ${area}

  initial_T = ${T_ambient}
  initial_p = ${p_initial}
  initial_vel = 0

  fp = fp
  closures = closures
  f = 0
  Hw = ${htc}

  rdg_slope_reconstruction = full
  scaling_factor_1phase = '1 1 1e-5'
[]

[FluidProperties]
  [fp]
    type = IdealGasFluidProperties
    emit_on_nan = none
  []
[]

[Closures]
  [closures]
    type = Closures1PhaseSimple
  []
[]

[Functions]
  [heating_flux_fn]
    type = PiecewiseConstant
    axis = z
    x = '0 0.2 0.8'
    y = '0 ${fparse power / (S_heated)} 0'
  []
  [cooling_flux_fn]
    type = PiecewiseConstant
    axis = z
    x = '0 0.2 0.8'
    y = '0 ${fparse -power / (S_cooled)} 0'
  []
[]

[Components]
  [heated_pipe]
    type = FlowChannel1Phase
    position = '0 0 0'
    orientation = '0 0 1'
  []
  [top_pipe]
    type = FlowChannel1Phase
    position = '0 0 1'
    orientation = '1 0 0'
  []
  [cooled_pipe]
    type = FlowChannel1Phase
    position = '1 0 1'
    orientation = '0 0 -1'
  []
  [bottom_pipe]
    type = FlowChannel1Phase
    position = '1 0 0'
    orientation = '-1 0 0'
  []

  [junction_heated_top]
    type = JunctionOneToOne1Phase
    connections = 'heated_pipe:out top_pipe:in'
  []
  [junction_top_cooled]
    type = JunctionOneToOne1Phase
    connections = 'top_pipe:out cooled_pipe:in'
  []
  [junction_cooled_bottom]
    type = JunctionOneToOne1Phase
    connections = 'cooled_pipe:out bottom_pipe:in'
  []
  [junction_bottom_heated]
    type = JunctionOneToOne1Phase
    connections = 'bottom_pipe:out heated_pipe:in'
  []

  [heating]
    type = HeatTransferFromHeatFlux1Phase
    flow_channel = 'heated_pipe'
    q_wall = heating_flux_fn
  []
  [cooling]
    type = HeatTransferFromHeatFlux1Phase
    flow_channel = 'cooled_pipe'
    q_wall = cooling_flux_fn
  []
[]

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

[Functions]
  # Time step size function. NaN gets into residual for debug mode on first step,
  # so that step gets cut to 0.5 s, but the rest of the transient can take 1 s time steps.
  [dt_fn]
    type = ParsedFunction
    expression = 'if(t < 1e-10, 0.5, 1.0)'
  []
[]

[Executioner]
  type = Transient

  scheme = bdf2
  start_time = 0
  end_time = 50
  [TimeStepper]
    type = FunctionDT
    function = dt_fn
  []

  steady_state_detection = true

  petsc_options_iname = '-pc_type'
  petsc_options_value = ' lu     '

  nl_rel_tol = 1e-10
  nl_abs_tol = 1e-10
  nl_max_its = 15
  l_tol = 1e-4
  l_max_its = 10
[]

[VectorPostprocessors]
  [heated_pipe_vpp]
    type = ElementValueSampler
    block = 'heated_pipe'
    variable = ${output_variables}
    sort_by = z
    execute_on = 'FINAL'
  []
  [top_pipe_vpp]
    type = ElementValueSampler
    block = 'top_pipe'
    variable = ${output_variables}
    sort_by = x
    execute_on = 'FINAL'
  []
  [cooled_pipe_vpp]
    type = ElementValueSampler
    block = 'cooled_pipe'
    variable = ${output_variables}
    sort_by = z
    execute_on = 'FINAL'
  []
  [bottom_pipe_vpp]
    type = ElementValueSampler
    block = 'bottom_pipe'
    variable = ${output_variables}
    sort_by = x
    execute_on = 'FINAL'
  []
[]

[Outputs]
  xml = true
  velocity_as_vector = false
  execute_on = 'FINAL'
[]

(modules/thermal_hydraulics/test/tests/problems/pressure_drop/pressure_drop.i)

nelem = 100
friction_factor = 1e4

area = 0.176752
mfr_final = 1.0
p_out = 7e6
T_in = 300
ramp_time = 5.0

[GlobalParams]
  gravity_vector = '0 0 0'
  initial_T = ${T_in}
  initial_p = ${p_out}
  initial_vel = 0
  closures = closures
  rdg_slope_reconstruction = full
  scaling_factor_1phase = '1 1 1e-5'
[]

[FluidProperties]
  [h2]
    type = IdealGasFluidProperties
    gamma = 1.3066
    molar_mass = 2.016e-3
    k = 0.437
    mu = 3e-5
  []
[]

[Closures]
  [closures]
    type = Closures1PhaseSimple
  []
[]

[Components]
  [bc_inlet]
    type = InletMassFlowRateTemperature1Phase
    input = 'ch_1:in'
    m_dot = 0 # This value is controlled by 'mfr_ctrl'
    T = ${T_in}
  []
  [ch_1]
    type = FlowChannel1Phase
    position = '0 0 0'
    orientation = '1 0 0'
    length = 1.0
    n_elems = ${nelem}
    A = ${area}
    f = ${friction_factor}
    fp = h2
  []
  [bc_outlet]
    type = Outlet1Phase
    input = 'ch_1:out'
    p = ${p_out}
  []
[]

[Functions]
  [mfr_fn]
    type = PiecewiseLinear
    x = '0 ${ramp_time}'
    y = '0 ${mfr_final}'
  []
[]

[ControlLogic]
  [mfr_ctrl]
    type = TimeFunctionComponentControl
    component = bc_inlet
    parameter = m_dot
    function = mfr_fn
  []
[]

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

[VectorPostprocessors]
  [pressure_vpp]
    type = ADSampler1DReal
    block = 'ch_1'
    property = 'p'
    sort_by = x
    execute_on = 'FINAL'
  []
[]

[Executioner]
  type = Transient

  scheme = bdf2
  start_time = 0
  end_time = 50
  dt = 1

  steady_state_detection = true
  steady_state_start_time = ${ramp_time}

  petsc_options_iname = '-pc_type'
  petsc_options_value = ' lu     '
  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-6
  nl_max_its = 15
  l_tol = 1e-4
  l_max_its = 10
[]

[Outputs]
  [csv]
    type = CSV
    create_final_symlink = true
    execute_on = 'FINAL'
  []
[]

(modules/thermal_hydraulics/test/tests/problems/sedov_blast_wave/sedov_blast_wave.i)

# This test problem is the Sedov blast wave test problem,
# which is a Riemann problem with the following parameters:
#   * domain = (0,1)
#   * gravity = 0
#   * EoS: Ideal gas EoS with gamma = 1.4, R = 0.71428571428571428571
#   * interface: x = 0.5
#   * typical end time: 0.15
# Left initial values:
#   * rho = 0.445
#   * vel = 0.692
#   * p = 3.52874226
# Right initial values:
#   * rho = 0.5
#   * vel = 0
#   * p = 0.571

[GlobalParams]
  gravity_vector = '0 0 0'

  closures = simple_closures
[]

[Functions]
  [p_ic_fn]
    type = PiecewiseConstant
    axis = x
    direction = right
    x = '0.0025                1'
    y = '1.591549333333333e+06 6.666666666666668e-09'
  []
  [T_ic_fn]
    type = PiecewiseConstant
    axis = x
    direction = right
    x = '0.0025                1'
    y = '2.228169066666667e+06 9.333333333333334e-09'
  []
[]

[FluidProperties]
  [fp]
    type = IdealGasFluidProperties
    gamma = 1.66666666666666666667
    molar_mass = 11.64024372
  []
[]

[Closures]
  [simple_closures]
    type = Closures1PhaseSimple
  []
[]

[Components]
  [pipe]
    type = FlowChannel1Phase

    fp = fp

    # geometry
    position = '0 0 0'
    orientation = '1 0 0'
    length = 1.0
    n_elems = 400
    A = 1.0

    # IC
    initial_T = T_ic_fn
    initial_p = p_ic_fn
    initial_vel = 0

    f = 0
  []

  [left_boundary]
    type = SolidWall1Phase
    input = 'pipe:in'
  []

  [right_boundary]
    type = FreeBoundary1Phase
    input = 'pipe:out'
  []
[]

[Executioner]
  type = Transient
  [TimeIntegrator]
    type = ExplicitSSPRungeKutta
    order = 2
  []
  solve_type = LINEAR

  l_tol = 1e-4

  nl_rel_tol = 1e-20
  nl_abs_tol = 1e-8
  nl_max_its = 60

  # run to t = 0.005
  start_time = 0.0
  dt = 1e-6
  num_steps = 5000
  abort_on_solve_fail = true
[]

[Outputs]
  file_base = 'sedov_blast_wave'
  velocity_as_vector = false
  execute_on = 'initial timestep_end'
  [out]
    type = Exodus
    show = 'p T vel'
  []
[]

(modules/thermal_hydraulics/test/tests/problems/sod_shock_tube/sod_shock_tube.i)

# This test problem is the classic Sod shock tube test problem,
# which is a Riemann problem with the following parameters:
#   * domain = (0,1)
#   * gravity = 0
#   * EoS: Ideal gas EoS with gamma = 1.4, R = 0.71428571428571428571
#   * interface: x = 0.5
#   * typical end time: 0.2
# Left initial values:
#   * rho = 1
#   * vel = 0
#   * p = 1
# Right initial values:
#   * rho = 0.125
#   * vel = 0
#   * p = 0.1
#
# The output can be viewed by opening Paraview with the state file `plot.pvsm`:
#   paraview --state=plot.pvsm
# This will plot the numerical solution against the analytical solution

[GlobalParams]
  gravity_vector = '0 0 0'

  rdg_slope_reconstruction = minmod

  closures = simple_closures
[]

[Functions]
  [p_ic_fn]
    type = PiecewiseConstant
    axis = x
    direction = right
    x = '0.5 1.0'
    y = '1.0 0.1'
  []

  [T_ic_fn]
    type = PiecewiseConstant
    axis = x
    direction = right
    x = '0.5 1.0'
    y = '1.4 1.12'
  []
[]

[FluidProperties]
  [fp]
    type = IdealGasFluidProperties
    gamma = 1.4
    molar_mass = 11.64024372
  []
[]

[Closures]
  [simple_closures]
    type = Closures1PhaseSimple
  []
[]

[Components]
  [pipe]
    type = FlowChannel1Phase

    fp = fp

    # geometry
    position = '0 0 0'
    orientation = '1 0 0'
    length = 1.0
    n_elems = 100
    A = 1.0

    # IC
    initial_T = T_ic_fn
    initial_p = p_ic_fn
    initial_vel = 0

    f = 0
  []

  [left_boundary]
    type = FreeBoundary1Phase
    input = 'pipe:in'
  []

  [right_boundary]
    type = FreeBoundary1Phase
    input = 'pipe:out'
  []
[]

[Executioner]
  type = Transient
  [TimeIntegrator]
    type = ExplicitSSPRungeKutta
  []
  solve_type = LINEAR

  l_tol = 1e-4

  nl_rel_tol = 1e-20
  nl_abs_tol = 1e-8
  nl_max_its = 60

  # run to t = 0.2
  start_time = 0.0
  dt = 1e-3
  num_steps = 200
  abort_on_solve_fail = true
[]

[Outputs]
  file_base = 'sod_shock_tube'
  velocity_as_vector = false
  execute_on = 'initial timestep_end'
  [out]
    type = Exodus
    show = 'rho p vel'
  []
[]

(modules/thermal_hydraulics/test/tests/problems/square_wave/square_wave.i)

# Square wave problem

[GlobalParams]
  gravity_vector = '0 0 0'

  rdg_slope_reconstruction = minmod

  closures = simple_closures
[]

[Functions]
  [T_ic_fn]
    type = PiecewiseConstant
    axis = x
    direction = right
    x = '0.1 0.6 1.0'
    y = '2.8 1.4 2.8'
  []
[]

[FluidProperties]
  [fp]
    type = IdealGasFluidProperties
    gamma = 1.4
    molar_mass = 11.64024372
  []
[]

[Closures]
  [simple_closures]
    type = Closures1PhaseSimple
  []
[]

[Components]
  [pipe]
    type = FlowChannel1Phase

    fp = fp

    # geometry
    position = '0 0 0'
    orientation = '1 0 0'
    length = 1.0
    n_elems = 400
    A = 1.0

    # IC
    initial_T = T_ic_fn
    initial_p = 1
    initial_vel = 1

    f = 0
  []

  [left_boundary]
    type = FreeBoundary1Phase
    input = 'pipe:in'
  []

  [right_boundary]
    type = FreeBoundary1Phase
    input = 'pipe:out'
  []
[]

[Executioner]
  type = Transient
  [TimeIntegrator]
    type = ExplicitSSPRungeKutta
    order = 2
  []
  solve_type = LINEAR

  l_tol = 1e-4

  nl_rel_tol = 1e-20
  nl_abs_tol = 1e-8
  nl_max_its = 60

  # run to t = 0.3
  start_time = 0.0
  dt = 2e-4
  num_steps = 1500
  abort_on_solve_fail = true
[]

[Outputs]
  file_base = 'square_wave'
  velocity_as_vector = false
  execute_on = 'initial timestep_end'
  [out]
    type = Exodus
    show = 'p T vel'
  []
[]

(modules/thermal_hydraulics/test/tests/problems/super_sonic_tube/test.i)

[GlobalParams]
  gravity_vector = '0 0 0'

  scaling_factor_1phase = '1 1e-2 1e-4'

  initial_p = 101325
  initial_T = 300
  initial_vel = 522.676

  closures = simple_closures

  spatial_discretization = cg
[]

[FluidProperties]
  [ig]
    type = IdealGasFluidProperties
    gamma = 1.41
    molar_mass = 0.028966206103678928
  []
[]

[Closures]
  [simple_closures]
    type = Closures1PhaseSimple
  []
[]

[Components]
  [pipe]
    type = FlowChannel1Phase
    fp = ig
    # geometry
    position = '0 0 0'
    orientation = '1 0 0'
    A = 1.
    D_h = 1.12837916709551
    f = 0.0

    length = 1
    n_elems = 100
  []

  [inlet]
    type = SupersonicInlet
    input = 'pipe:in'
    p = 101325
    T = 300.0
    vel = 522.676
  []

  [outlet]
    type = FreeBoundary1Phase
    input = 'pipe:out'
  []
[]

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

[Executioner]
  type = Transient
  scheme = 'bdf2'

  start_time = 0
  dt = 1e-5
  num_steps = 10
  abort_on_solve_fail = true

  solve_type = 'PJFNK'
  nl_rel_tol = 1e-9
  nl_abs_tol = 1e-8
  nl_max_its = 30

  l_tol = 1e-3
  l_max_its = 100

  [Quadrature]
    type = TRAP
    order = FIRST
  []
[]


[Outputs]
  [out]
    type = Exodus
  []
[]

(modules/thermal_hydraulics/test/tests/problems/water_hammer/3eqn.i)

[GlobalParams]
  gravity_vector = '0 0 0'

  initial_T = 517.252072255516
  initial_vel = 0

  scaling_factor_1phase = '1.e0 1.e0 1.e-2'

  closures = simple_closures
[]

[FluidProperties]
  [eos]
    type = StiffenedGasFluidProperties
    gamma = 2.35
    q = -1167e3
    q_prime = 0
    p_inf = 1.e9
    cv = 1816
  []
[]

[Closures]
  [simple_closures]
    type = Closures1PhaseSimple
  []
[]

[Functions]
  [p_fn]
    type = PiecewiseConstant
    axis = x
    x = '0      0.5    1'
    y = '7.5e6  6.5e6  6.5e6'
  []
[]

[Components]
  [pipe1]
    type = FlowChannel1Phase
    fp = eos
    # geometry
    position = '0 0 0'
    orientation = '1 0 0'
    length = 1
    n_elems = 200

    A = 1.907720E-04
    D_h = 1.698566E-02
    f = 0.

    initial_p = p_fn
  []

  # BCs
  [left]
    type = SolidWall1Phase
    input = 'pipe1:in'
  []

  [right]
    type = SolidWall1Phase
    input = 'pipe1:out'
  []
[]

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

[Executioner]
  type = Transient
  scheme = 'bdf2'

  start_time = 0
  dt = 1e-5
  num_steps = 10
  abort_on_solve_fail = true

  solve_type = 'PJFNK'
  line_search = 'basic'
  nl_rel_tol = 1e-9
  nl_abs_tol = 1e-9
  nl_max_its = 10

  l_tol = 1e-3
  l_max_its = 100
[]

[Outputs]
  velocity_as_vector = false
  [out]
    type = Exodus
  []
[]

(modules/thermal_hydraulics/test/tests/problems/woodward_colella_blast_wave/woodward_colella_blast_wave.i)

# Woodward-Colella blast wave problem

[GlobalParams]
  gravity_vector = '0 0 0'

  closures = simple_closures
[]

[Functions]
  [p_ic_fn]
    type = PiecewiseConstant
    axis = x
    direction = right
    x = '0.1  0.9  1.0'
    y = '1000 0.01 100'
  []
  [T_ic_fn]
    type = PiecewiseConstant
    axis = x
    direction = right
    x = '0.1  0.9   1.0'
    y = '1400 0.014 140'
  []
[]

[FluidProperties]
  [fp]
    type = IdealGasFluidProperties
    gamma = 1.4
    molar_mass = 11.64024372
  []
[]

[Closures]
  [simple_closures]
    type = Closures1PhaseSimple
  []
[]

[Components]
  [pipe]
    type = FlowChannel1Phase

    fp = fp

    # geometry
    position = '0 0 0'
    orientation = '1 0 0'
    length = 1.0
    n_elems = 500
    A = 1.0

    # IC
    initial_T = T_ic_fn
    initial_p = p_ic_fn
    initial_vel = 0

    f = 0
  []

  [left_wall]
    type = SolidWall1Phase
    input = 'pipe:in'
  []

  [right_wall]
    type = SolidWall1Phase
    input = 'pipe:out'
  []
[]

[Executioner]
  type = Transient
  [TimeIntegrator]
    type = ExplicitSSPRungeKutta
    order = 2
  []
  solve_type = LINEAR

  l_tol = 1e-4

  nl_rel_tol = 1e-20
  nl_abs_tol = 1e-8
  nl_max_its = 60

  # run to t = 0.038
  start_time = 0.0
  dt = 1e-5
  num_steps = 3800
  abort_on_solve_fail = true
[]

[Outputs]
  file_base = 'woodward_colella_blast_wave'
  velocity_as_vector = false
  execute_on = 'initial timestep_end'
  [out]
    type = Exodus
    show = 'p T vel'
  []
[]