Test Problems

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

Name	File(s)	Description
Abrupt area change	[without_junction.i], [with_junction.i]	Liquid flow through an abrupt area change. A comparison is provided to compare with and without using a junction.
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
Three-pipe shock	[input]	Water flow through a junction splitting between two pipes in the same orientation.
Water hammer	[input]	Water hammer test problem
William-Louis	[3pipes_open.i], [4pipes_closed.i]	Test problems in a multi-pipe shock tube, compared against experimental data.
Woodward-Colella blast wave	[input]	Woodward-Colella blast wave problem

(contrib/moose/modules/thermal_hydraulics/test/tests/problems/abrupt_area_change_liquid/without_junction.i)

# Version without junction

!include base_params.i
!include base.i

[Functions]
  [A_fn]
    type = PiecewiseConstant
    axis = x
    x = '0 ${xR}'
    y = '${AL} ${AR}'
  []
[]

[Components]
  [pipe]
    type = FlowChannel1Phase
    position = '0 0 0'
    orientation = '1 0 0'
    length = ${fparse lengthL + lengthR}
    n_elems = ${fparse NL + NR}
    A = A_fn
  []

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

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

[VectorPostprocessors]
  [vpp]
    type = ADSampler1DReal
    block = 'pipe'
    property = 'p vel'
    sort_by = x
    execute_on = 'FINAL'
  []
[]

[Outputs]
  file_base = 'without_junction'
[]

(contrib/moose/modules/thermal_hydraulics/test/tests/problems/abrupt_area_change_liquid/with_junction.i)

# Version with junction

!include base_params.i

# For equivalent comparison to the unsplit case, we reduce the length of the
# right pipe by one element
dx = ${fparse lengthL / NL}
NR_minus_junction = ${fparse NR - 1}
lengthR_minus_junction = ${fparse lengthR - dx}
xR_minus_junction = ${fparse xR + dx}

xJ = ${fparse lengthL + 0.5 * dx}
AJ = ${fparse AL + AR}
RJ = ${fparse sqrt(AJ / (4 * pi))} # A = 4 pi R^2
VJ = ${fparse 4/3 * pi * RJ^3}

!include base.i

[Components]
  [pipeL]
    type = FlowChannel1Phase
    position = '0 0 0'
    orientation = '1 0 0'
    length = ${lengthL}
    n_elems = ${NL}
    A = ${AL}
  []

  [pipeR]
    type = FlowChannel1Phase
    position = '${xR_minus_junction} 0 0'
    orientation = '1 0 0'
    length = ${lengthR_minus_junction}
    n_elems = ${NR_minus_junction}
    A = ${AR}
  []

  [junction]
    type = VolumeJunction1Phase
    connections = 'pipeL:out pipeR:in'
    position = '${xJ} 0 0'
    volume = ${VJ}
    initial_vel_x = 0
    initial_vel_y = 0
    initial_vel_z = 0
    scaling_factor_rhoEV = 1e-5
    apply_velocity_scaling = true
  []

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

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

[VectorPostprocessors]
  [vpp]
    type = ADSampler1DReal
    block = 'pipeL pipeR'
    property = 'p vel'
    sort_by = x
    execute_on = 'FINAL'
  []
[]

[Outputs]
  file_base = 'with_junction'
[]

(contrib/moose/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'
[]

(contrib/moose/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'
  []
[]

(contrib/moose/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'
  []
[]

(contrib/moose/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'
  []
[]

(contrib/moose/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
[]

(contrib/moose/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'
  []
[]

(contrib/moose/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'
  []
[]

(contrib/moose/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'
  []
[]

(contrib/moose/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'
  []
[]

(contrib/moose/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
  []
[]

(contrib/moose/modules/thermal_hydraulics/test/tests/problems/three_pipe_shock/three_pipe_shock.i)

# Test 8 from the following reference:
#
#   F. Daude, P. Galon. A Finite-Volume approach for compressible single- and
#   two-phase flows in flexible pipelines with fluid-structure interaction.
#   Journal of Computational Physics 362 (2018) 375-408.

L1 = 10
L2 = 3
L3 = 5

xJ = ${L1}
x_p1 = ${fparse xJ - 1.05}
x_p2 = ${fparse xJ + 0.15}
x_p3 = ${fparse xJ + 0.95}

N1 = 1000
N2 = 300
N3 = 500

D1 = 0.35682482
D2 = 0.19544100
D3 = 0.35682482

A1 = ${fparse 0.25 * pi * D1^2}
A2 = ${fparse 0.25 * pi * D2^2}
A3 = ${fparse 0.25 * pi * D3^2}
AJ = ${fparse A1 + A2 + A3}
RJ = ${fparse sqrt(AJ / (4 * pi))} # A = 4 pi R^2
VJ = ${fparse 4/3 * pi * RJ^3}

y2 = 1
y3 = -1

gamma = 2.23
p_inf = 1e9 # denoted by "pi" in reference
q = 0
cv = 2500 # arbitrary value; not given in reference

CFL = 0.8
t_end = 0.01

p_out = 80e5

initial_p = ${p_out}
initial_T = 327.1864956 # reference has rho = 1001.89 kg/m^3
initial_vel1 = 1
initial_vel2 = 0.769
initial_vel3 = 0.769

[GlobalParams]
  gravity_vector = '0 0 0'

  initial_T = ${initial_T}
  initial_p = ${initial_p}

  fp = fp
  closures = closures
  f = 0

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

[FluidProperties]
  [fp]
    type = StiffenedGasFluidProperties
    gamma = ${gamma}
    p_inf = ${p_inf}
    q = ${q}
    cv = ${cv}
  []
[]

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

[Components]
  [pipe1]
    type = FlowChannel1Phase
    position = '0 0 0'
    orientation = '1 0 0'
    length = ${L1}
    n_elems = ${N1}
    A = ${A1}
    initial_vel = ${initial_vel1}
  []
  [pipe2]
    type = FlowChannel1Phase
    position = '${xJ} ${y2} 0'
    orientation = '1 0 0'
    length = ${L2}
    n_elems = ${N2}
    A = ${A2}
    initial_vel = ${initial_vel2}
  []
  [pipe3]
    type = FlowChannel1Phase
    position = '${xJ} ${y3} 0'
    orientation = '1 0 0'
    length = ${L3}
    n_elems = ${N3}
    A = ${A3}
    initial_vel = ${initial_vel3}
  []

  [junction]
    type = VolumeJunction1Phase
    connections = 'pipe1:out pipe2:in pipe3:in'
    position = '${xJ} 0 0'
    volume = ${VJ}
    initial_vel_x = ${initial_vel2} # ?
    initial_vel_y = 0
    initial_vel_z = 0
    scaling_factor_rhoEV = 1e-5
    apply_velocity_scaling = true
  []

  [outlet1]
    type = Outlet1Phase
    input = 'pipe1:in'
    p = ${p_out}
  []
  [outlet2]
    type = Outlet1Phase
    input = 'pipe2:out'
    p = ${p_out}
  []
  [wall3]
    type = SolidWall1Phase
    input = 'pipe3:out'
  []
[]

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

[Postprocessors]
  [dt_cfl]
    type = ADCFLTimeStepSize
    CFL = ${CFL}
    vel_names = 'vel'
    c_names = 'c'
  []
  [p1]
    type = PointValue
    variable = p
    point = '${x_p1} 0 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p2]
    type = PointValue
    variable = p
    point = '${x_p2} ${y2} 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p3]
    type = PointValue
    variable = p
    point = '${x_p3} ${y3} 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[Executioner]
  type = Transient

  start_time = 0
  end_time = ${t_end}
  [TimeStepper]
    type = PostprocessorDT
    postprocessor = dt_cfl
  []
  [TimeIntegrator]
    type = ActuallyExplicitEuler
  []

  solve_type = LINEAR

  petsc_options_iname = '-pc_type'
  petsc_options_value = ' lu     '

  l_tol = 1e-4
  l_max_its = 10
[]

[Times]
  [output_times]
    type = TimeIntervalTimes
    time_interval = 1e-4
  []
[]

[Outputs]
  file_base = 'three_pipe_shock'
  [csv]
    type = CSV
    show = 'p1 p2 p3'
    sync_only = true
    sync_times_object = output_times
  []
[]

(contrib/moose/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
  []
[]

(contrib/moose/modules/thermal_hydraulics/test/tests/problems/william_louis/3pipes_open.i)

# Junction of 3 pipes:
#
#   1     3
# -----*-----
#      | 2
#
# The left end of Pipe 1 is a high-pressure region, and the rest of the system
# is at a low pressure.
#
# Pipe 1 is closed, while Pipes 2 and 3 are open.

end_time = 0.07

D_pipe = 0.01
A_pipe = ${fparse 0.25 * pi * D_pipe^2}
length_pipe1_HP = 0.53
length_pipe1_LP = 3.10
length_pipe2 = 2.595
length_pipe3 = 1.725

x_junction = ${fparse length_pipe1_HP + length_pipe1_LP}

# Numbers of elements correspond to dx ~ 1/3 cm
n_elems_pipe1_HP = 159
n_elems_pipe1_LP = 930
n_elems_pipe2 = 779
n_elems_pipe3 = 518

S_junction = ${fparse 3 * A_pipe}
r_junction = ${fparse sqrt(S_junction / (4 * pi))}
V_junction = ${fparse 4/3 * pi * r_junction^3}

p_low = 1e5
p_high = 1.15e5

T_low  = 283.5690633 # at p = 1e5 Pa,    rho = 1.23 kg/m^3
T_high = 283.5690633 # at p = 1.15e5 Pa, rho = 1.4145 kg/m^3

cfl = 0.95

[GlobalParams]
  # common FlowChannel1Phase parameters
  A = ${A_pipe}
  initial_vel = 0
  fp = fp_air
  closures = closures
  f = 0
  gravity_vector = '0 0 0'
  scaling_factor_1phase = '1 1 1e-5'
[]

[FluidProperties]
  [fp_air]
    type = IdealGasFluidProperties
    gamma = 1.4
    molar_mass = 0.029
  []
[]

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

[Functions]
  [initial_T_pipe1_fn]
    type = PiecewiseConstant
    axis = x
    x = '0 ${length_pipe1_HP}'
    y = '${T_high} ${T_low}'
  []
  [initial_p_pipe1_fn]
    type = PiecewiseConstant
    axis = x
    x = '0 ${length_pipe1_HP}'
    y = '${p_high} ${p_low}'
  []
[]

[Components]
  [pipe1_wall]
    type = SolidWall1Phase
    input = 'pipe1:in'
  []
  [pipe1]
    type = FlowChannel1Phase
    position = '0 0 0'
    orientation = '1 0 0'
    length = '${length_pipe1_HP} ${length_pipe1_LP}'
    n_elems = '${n_elems_pipe1_HP} ${n_elems_pipe1_LP}'
    initial_p = initial_p_pipe1_fn
    initial_T = initial_T_pipe1_fn
  []

  [junction]
    type = VolumeJunction1Phase
    position = '${x_junction} 0 0'
    connections = 'pipe1:out pipe2:in pipe3:in'
    initial_p = ${p_low}
    initial_T = ${T_low}
    initial_vel_x = 0
    initial_vel_y = 0
    initial_vel_z = 0
    volume = ${V_junction}
    scaling_factor_rhoEV = 1e-5
    apply_velocity_scaling = true
  []

  [pipe2]
    type = FlowChannel1Phase
    position = '${x_junction} 0 0'
    orientation = '0 -1 0'
    length = ${length_pipe2}
    n_elems = ${n_elems_pipe2}
    initial_p = ${p_low}
    initial_T = ${T_low}
  []
  [pipe2_outlet]
    type = Outlet1Phase
    input = 'pipe2:out'
    p = ${p_low}
  []

  [pipe3]
    type = FlowChannel1Phase
    position = '${x_junction} 0 0'
    orientation = '1 0 0'
    length = ${length_pipe3}
    n_elems = ${n_elems_pipe3}
    initial_p = ${p_low}
    initial_T = ${T_low}
  []
  [pipe3_outlet]
    type = Outlet1Phase
    input = 'pipe3:out'
    p = ${p_low}
  []
[]

[Postprocessors]
  [cfl_dt]
    type = ADCFLTimeStepSize
    CFL = ${cfl}
    c_names = 'c'
    vel_names = 'vel'
  []
  [p_pipe1_048]
    type = PointValue
    variable = p
    point = '${fparse x_junction - 0.48} 0 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p_pipe2_052]
    type = PointValue
    variable = p
    point = '${fparse x_junction} -0.52 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p_pipe3_048]
    type = PointValue
    variable = p
    point = '${fparse x_junction + 0.48} 0 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

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

[Executioner]
  type = Transient
  end_time = ${end_time}

  [TimeIntegrator]
    type = ExplicitSSPRungeKutta
    order = 1
  []

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = cfl_dt
  []
  abort_on_solve_fail = true

  solve_type = LINEAR
[]

[Times]
  [output_times]
    type = TimeIntervalTimes
    time_interval = 7e-4
  []
[]

[Outputs]
  file_base = '3pipes_open'
  [csv]
    type = CSV
    show = 'p_pipe1_048 p_pipe2_052 p_pipe3_048'
    sync_only = true
    sync_times_object = output_times
  []
  [console]
    type = Console
    execute_postprocessors_on = 'NONE'
  []
[]

(contrib/moose/modules/thermal_hydraulics/test/tests/problems/william_louis/4pipes_closed.i)

# Junction of 4 pipes:
#
#        4
#        |
# 1 -----*----- 3
#        |
#        2
#
# The left end of Pipe 1 is a high-pressure region, and the rest of the system
# is at a low pressure.
#
# All pipes are closed.

end_time = 0.07

D_pipe = 0.01
A_pipe = ${fparse 0.25 * pi * D_pipe^2}
length_pipe1_HP = 0.53
length_pipe1_LP = 3.10
length_pipe2 = 2.595
length_pipe3 = 1.725
length_pipe4 = 0.845

x_junction = ${fparse length_pipe1_HP + length_pipe1_LP}

# Numbers of elements correspond to dx ~ 1/3 cm
n_elems_pipe1_HP = 159
n_elems_pipe1_LP = 930
n_elems_pipe2 = 779
n_elems_pipe3 = 518
n_elems_pipe4 = 254

S_junction = ${fparse 4 * A_pipe}
r_junction = ${fparse sqrt(S_junction / (4 * pi))}
V_junction = ${fparse 4/3 * pi * r_junction^3}

p_low = 1e5
p_high = 1.15e5

T_initial = 283.5690633 # at p = 1e5 Pa, rho = 1.23 kg/m^3

cfl = 0.95

[GlobalParams]
  # common FlowChannel1Phase parameters
  A = ${A_pipe}
  initial_T = ${T_initial}
  initial_vel = 0
  fp = fp_air
  closures = closures
  f = 0
  gravity_vector = '0 0 0'
  scaling_factor_1phase = '1 1 1e-5'
[]

[FluidProperties]
  [fp_air]
    type = IdealGasFluidProperties
    gamma = 1.4
    molar_mass = 0.029
  []
[]

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

[Functions]
  [initial_p_pipe1_fn]
    type = PiecewiseConstant
    axis = x
    x = '0 ${length_pipe1_HP}'
    y = '${p_high} ${p_low}'
  []
[]

[Components]
  [pipe1_wall]
    type = SolidWall1Phase
    input = 'pipe1:in'
  []
  [pipe1]
    type = FlowChannel1Phase
    position = '0 0 0'
    orientation = '1 0 0'
    length = '${length_pipe1_HP} ${length_pipe1_LP}'
    n_elems = '${n_elems_pipe1_HP} ${n_elems_pipe1_LP}'
    initial_p = initial_p_pipe1_fn
  []

  [junction]
    type = VolumeJunction1Phase
    position = '${x_junction} 0 0'
    connections = 'pipe1:out pipe2:in pipe3:in pipe4:in'
    initial_p = ${p_low}
    initial_T = ${T_initial}
    initial_vel_x = 0
    initial_vel_y = 0
    initial_vel_z = 0
    volume = ${V_junction}
    scaling_factor_rhoEV = 1e-5
    apply_velocity_scaling = true
  []

  [pipe2]
    type = FlowChannel1Phase
    position = '${x_junction} 0 0'
    orientation = '0 -1 0'
    length = ${length_pipe2}
    n_elems = ${n_elems_pipe2}
    initial_p = ${p_low}
  []
  [pipe2_wall]
    type = SolidWall1Phase
    input = 'pipe2:out'
  []

  [pipe3]
    type = FlowChannel1Phase
    position = '${x_junction} 0 0'
    orientation = '1 0 0'
    length = ${length_pipe3}
    n_elems = ${n_elems_pipe3}
    initial_p = ${p_low}
  []
  [pipe3_wall]
    type = SolidWall1Phase
    input = 'pipe3:out'
  []

  [pipe4]
    type = FlowChannel1Phase
    position = '${x_junction} 0 0'
    orientation = '0 1 0'
    length = ${length_pipe4}
    n_elems = ${n_elems_pipe4}
    initial_p = ${p_low}
  []
  [pipe4_wall]
    type = SolidWall1Phase
    input = 'pipe4:out'
  []
[]

[Postprocessors]
  [cfl_dt]
    type = ADCFLTimeStepSize
    CFL = ${cfl}
    c_names = 'c'
    vel_names = 'vel'
  []
  [p_pipe1_048]
    type = PointValue
    variable = p
    point = '${fparse x_junction - 0.48} 0 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p_pipe2_052]
    type = PointValue
    variable = p
    point = '${fparse x_junction} -0.52 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p_pipe3_048]
    type = PointValue
    variable = p
    point = '${fparse x_junction + 0.48} 0 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p_pipe4_043]
    type = PointValue
    variable = p
    point = '${fparse x_junction} 0.43 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

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

[Executioner]
  type = Transient
  end_time = ${end_time}

  [TimeIntegrator]
    type = ExplicitSSPRungeKutta
    order = 1
  []

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = cfl_dt
  []
  abort_on_solve_fail = true

  solve_type = LINEAR
[]

[Times]
  [output_times]
    type = TimeIntervalTimes
    time_interval = 7e-4
  []
[]

[Outputs]
  file_base = '4pipes_closed'
  [csv]
    type = CSV
    show = 'p_pipe1_048 p_pipe2_052 p_pipe3_048 p_pipe4_043'
    sync_only = true
    sync_times_object = output_times
  []
  [console]
    type = Console
    execute_postprocessors_on = 'NONE'
  []
[]

(contrib/moose/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'
  []
[]