braamvandyk / SteamTables.jl

Julia implementation of steam tables according to IAPWS Industrial Formulation (1997)
Other
12 stars 4 forks source link

SpecificS_Ph failure with negative log #27

Open tkoenig1 opened 6 months ago

tkoenig1 commented 6 months ago

Carried over from a comment in #24 .

julia> SpecificS_Ph(4.5,2900)
ERROR: DomainError with -201.09208450864833:
sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
Stacktrace:
 [1] throw_complex_domainerror(f::Symbol, x::Float64)
   @ Base.Math ./math.jl:33
 [2] sqrt
   @ ./math.jl:686 [inlined]
 [3] B2bc
   @ ~/.julia/packages/SteamTables/hNd1l/src/SteamTables.jl:186 [inlined]
 [4] RegionID_Ph(P::Float64, h::Int64)
   @ SteamTables ~/.julia/packages/SteamTables/hNd1l/src/SteamTables.jl:2310
 [5] SpecificS_Ph(P::Float64, h::Int64)
   @ SteamTables ~/.julia/packages/SteamTables/hNd1l/src/SteamTables.jl:3119
 [6] top-level scope
   @ REPL[2]:1

julia> 
tkoenig1 commented 6 months ago

Here's a short script to identify dubious regions, plus its output.

Recognized two-phase flow is black. There are several areas which are white (indicating missing values), but also some areas which go into the two-phase region, which they should not.

gaps

using Plots,Colors
using SteamTables
using FlexiMaps

function my_H_ps(p,s)
  tmp = missing
  try
    tmp = SpecificH_Ps(p,s)
  catch
    # Set residual two-phase region to -1000, so it is visible
    if (p > Pc)
      return missing
    end
  end
  if !ismissing(tmp)
    return tmp
  end
  T_sat = Tsat(p)
  S_L = SatSL(T_sat)
  S_V = SatSV(T_sat)
  if (s < S_L || s > S_V)
#    println("p = ",p, " Pc = ", Pc, " T_sat = ",T_sat)
    return missing
  end
  return -4000.0

  return tmp
end

n_pix = 512

p_low = P3*1.1
p_high = 50.0
s_low = 0.01
s_high = 12.0

pv = maprange(log,p_low,p_high,n_pix)
sv = range(s_low,s_high,n_pix)

H = [my_H_ps(p,s) for p in pv, s in sv]

p = contourf(sv,pv,H,yscale=:log10)

savefig(p,"gaps.png")