temperature out of range in Water97FluidProperties

[wangzhaohao]

I want to do a thermal hydraulics problem which is in follwing:

in the THM part, there is input file.

inlet_T = 308.15 # K
m_dot_in = 12.13 # kg/s
press = '${fparse 101325 + 1e3 * 9.8 * 7.4}' # Pa

# pipe
A_pipe = '${fparse 2 * 0.0215^2 * 1.7312 - pi * 0.0181^2}'
D_h = '${fparse 4 * A_pipe / (pi * 0.0352)}'
length = 0.39 #m

[GlobalParams]
  initial_p = ${press}
  initial_vel = 0.204#m/s
  initial_T = ${inlet_T}
  gravity_vector = '0 0 0'

  rdg_slope_reconstruction = full
  closures = thm_closures
  fp = water
[]

[Closures]
  [thm_closures]
    type = Closures1PhaseNone
  []
[]

[FluidProperties]
  [water]
    type = Water97FluidProperties
  []
[]


[AuxVariables]
  [T_wall]
    family = MONOMIAL
    order = CONSTANT
  []
  [T]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [Tw_aux]
    type = ADMaterialRealAux
    block = pipe
    variable = T_wall
    property = T_wall
  []
[]

[Materials]
  # wall friction closure
  [f_mat]
    type = ADWallFrictionChurchillMaterial
    block = pipe
    D_h = D_h
    f_D = f_D
    vel = vel
    rho = rho
    mu = mu
  []
  
  # Re = \rho * vel * D_h / mu
  [Re_mat]
    type = ADReynoldsNumberMaterial
    Re = Re 
    rho = rho 
    D_h = D_h
    mu = mu
    vel = vel
    block = pipe
  []
  # Pr = C_p  * mu / k 
  [Pr_mat]
    type = ADPrandtlNumberMaterial
    Pr = Pr 
    cp = cp
    mu = mu
    k = k
    block = pipe
  []
  [Nu_mat]
    type = ADParsedMaterial
    block = pipe
    # Dittus-Boelter
    expression = '0.021 * pow(Re, 0.8) * pow(Pr, 0.4)'
    property_name = 'Nu'
    material_property_names = 'Re Pr'
  []
  [Hw_mat]
    type = ADConvectiveHeatTransferCoefficientMaterial
    D_h = D_h
    k = k
    Nu = Nu
    Hw = Hw
    block = pipe
  []
  [T_wall]
    type = ADTemperatureWall3EqnMaterial
    Hw = Hw
    T = T
    q_wall = q_wall
  []
[]


[Components]
  [inlet]
    type = InletMassFlowRateTemperature1Phase
    input = 'pipe:in'
    m_dot = ${m_dot_in}
    T = ${inlet_T}
  []
  [pipe]
    type = FlowChannel1Phase
    position = '0 0 0'
    orientation = '0 0 1'
    length = ${length}
    n_elems = 39 
    A = ${A_pipe} 
    D_h = ${D_h} 
  []
  [hxconn]
    type = HeatTransferFromExternalAppHeatFlux1Phase
    flow_channel = pipe
    P_hf = ${fparse pi * 0.0181 * 2}
  []
  [outlet]
    type = Outlet1Phase
    input = 'pipe:out'
    p = ${press}
  []
[]

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

[Executioner]
  type = Transient
  line_search = basic

  start_time = 0
  end_time = 31536000

  steady_state_detection = true
  steady_state_tolerance = 1e-8

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-8
  nl_max_its = 25
  [TimeStepper]
    type = FunctionDT
    function = dts
  []
[]

[Functions]
  [dts]  
    type = PiecewiseLinear
    x = '0 60 3600 86400 31536000'
    y = '10 10 360 86400 86400'
  []
[]

[Outputs]
  exodus = true

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

the wrong is

libMesh terminating:
Temperature -376998 is out of range in water: inRegion()
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 0
[unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=1

and I try to use gdb follow code, I find the temperature is right in RhoEAFromPressureTemperatureVelocityIC, why temperature wrong in Water97FluidProperties? Thanks for you help.

[GiudGiud]

Hello

Do you have a backtrace from gdb? This would help us look at the routines in Water97FP

Guillaume

[wangzhaohao]

Sure, I am not very understand the GDB, so i write the steps which I have done clearly.
I reed the debug part in MOOSE web,

b MPI_Abort
bt

and part of output is

 0x00007fffef624092 in Water97FluidProperties::inRegion (this=0x555555d0def0, pressure=173845, temperature=-376997.58328186366)
    at /home/wang/projects/moose/modules/fluid_properties/src/fluidproperties/Water97FluidProperties.C:1157
#6  0x00007fffef680d22 in Water97FluidProperties::rho_from_p_T_template<MetaPhysicL::DualNumber<double, MetaPhysicL::NumberArray<2ul, double>, false> > (this=0x555555d0def0, pressure=..., temperature=...)
    at /home/wang/projects/moose/modules/fluid_properties/build/header_symlinks/Water97FluidProperties.h:2557
#7  0x00007fffef6636af in Water97FluidProperties::rho_from_p_T_template<double>(double const&, double const&, double&, double&, double&) const::{lambda(auto:1 const&, auto:2 const&)#1}::operator()<MetaPhysicL::DualNumber<double, MetaPhysicL::NumberArray<2ul, double>, false>, MetaPhysicL::DualNumber<double, MetaPhysicL::NumberArray<2ul, double>, false> >(MetaPhysicL::DualNumber<double, MetaPhysicL::NumberArray<2ul, double>, false> const&, MetaPhysicL::DualNumber<double, MetaPhysicL::NumberArray<2ul, double>, false> const&) const (this=0x555555d0def0, pressure=..., temperature=...)
    at /home/wang/projects/moose/modules/fluid_properties/build/header_symlinks/Water97FluidProperties.h:2003
#15 0x00007ffff105de05 in RhoEAFromPressureTemperatureFunctionVelocityIC::value (this=0x555555f28320, p=...)
    at /home/wang/projects/moose/modules/thermal_hydraulics/src/ics/RhoEAFromPressureTemperatureFunctionVelocityIC.C:46

so I set breakpoints in RhoEAFromPressureTemperatureFunctionVelocityI and Water97FluidProperties

Thread 1 "pulse-dbg" hit Breakpoint 3, RhoEAFromPressureTemperatureFunctionVelocityIC::value (this=0x555555f28320, p=...)
    at /home/wang/projects/moose/modules/thermal_hydraulics/src/ics/RhoEAFromPressureTemperatureFunctionVelocityIC.C:46
46        const Real e = _fp.e_from_p_rho(_p[_qp], rho);
(gdb) info locals
vel = 0.20399999999999999
rho = 994.07047349434276
e = 6.9533462538108351e-310
(gdb) print _T[_qp]
$1 = (const double &) @0x5555562eda50: 308.14999999999992
(gdb) n
Currently Executing
  Performing Initial Setup
    Projecting Initial Solutions
Thread 1 "pulse-dbg" hit Breakpoint 4, Water97FluidProperties::inRegion (this=0x555555d0def0, pressure=173845, temperature=-376997.58328186366)
    at /home/wang/projects/moose/modules/fluid_properties/src/fluidproperties/Water97FluidProperties.C:1157
1157        mooseException("Temperature ", temperature, " is out of range in ", name(), ": inRegion()");
(gdb) info locals
region = 32767
(gdb) print temperature
$2 = -376997.58328186366
(gdb) 

I don't konw the Water97 how to get the temperature, because I not find such as coupledVariable function in Water97FluidProperties.C.

[wangzhaohao]

I find where set the temperature, I resee the GDB and find up method.

    at /home/wang/projects/moose/modules/fluid_properties/build/header_symlinks/Water97FluidProperties.h:2436
2436      const auto temperature = T_drhodT_from_p_rho(p, rho).first;

and use step to see the T_drhodT_from_p_rho(p, rho).first
where the code is

template <typename T>
std::pair<T, T>
Water97FluidProperties::T_drhodT_from_p_rho(const T & p, const T & rho) const
{
  // For NewtonSolve
  // y = rho
  // x = p
  // z = T
  auto lambda =
      [this](const T & p, const T & temperature, T & rho, T & drho_dp, T & drho_dtemperature)
  { rho_from_p_T_template(p, temperature, rho, drho_dp, drho_dtemperature); };
  return FluidPropertiesUtils::NewtonSolve(
      p, rho, _T_initial_guess, _tolerance, lambda, name() + "::T_drhodT_from_p_rho");
}

so something wrong in FluidPropertiesUtils::NewtonSolve function, try to output the parameters variable p rho _T_intial_guess tolerance, all is right. SO, I need to have a look at NewtonSolve?

[GiudGiud]

you can step in NewtonSolve in NewtonInversion to see how the newton solve is doing.
This is a 1D newton solve so it should be very simple.
I dont imagine something is wrong there though.

p rho _T_intial_guess tolerance, all is right.

what are p and rho ? If you compute rho(p, T=-30000) do you get close to the rho there?

[wangzhaohao]

what are p and rho ? If you compute rho(p, T=-30000) do you get close to the rho there?

i output the p rho, seems in reasonable range

(gdb) print p
$12 = (const double &) @0x7fffffff9e80: 173845
(gdb) print rho
$13 = (const double &) @0x7fffffff9e78: 994.07047349434276
(gdb) print _T_initial_gusss
No symbol "_T_initial_gusss" in current context.
(gdb) print _T_initial_guess
$14 = 400

p is the press = '${fparse 101325 + 1e3 * 9.8 * 7.4}' # Pa

[wangzhaohao]

I find a useful method
if i change the T_initial_guess, code will continue, so do you have something about how to set it?

[GiudGiud]

T_initial_guess can be set in the input file.
No there is no guarantee.

I would want to know if there are two solutions maybe for the Newton inversion.
One with T > 0 and one with T = -30000