braamvandyk / SteamTables.jl

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

failures in SpecificH(P, T) for various points in the domain #14

Closed stevengj closed 2 years ago

stevengj commented 2 years ago

I tried plotting SpecificH for a grid of pressure and temperature points, and noticed that it throws seemingly spurious exceptions at a bunch of points in the domain. In the plot below, the points where it threw errors are colored white: image

Sample errors:

julia> SpecificH(30, 628)
ERROR: Roots.ConvergenceFailed("Algorithm failed to converge")
Stacktrace:
 [1] #find_zero#15
   @ ~/.julia/packages/Roots/bSIXy/src/find_zero.jl:711 [inlined]
 [2] find_zero
   @ ~/.julia/packages/Roots/bSIXy/src/find_zero.jl:709 [inlined]
 [3] #find_zero#16
   @ ~/.julia/packages/Roots/bSIXy/src/find_zero.jl:719 [inlined]
 [4] find_zero
   @ ~/.julia/packages/Roots/bSIXy/src/find_zero.jl:719 [inlined]
 [5] Region3(Output::Symbol, P::Int64, T::Int64)
   @ SteamTables ~/.julia/packages/SteamTables/pMvX5/src/SteamTables.jl:1757
 [6] SpecificH(P::Int64, T::Int64)
   @ SteamTables ~/.julia/packages/SteamTables/pMvX5/src/SteamTables.jl:3018
 [7] top-level scope
   @ REPL[4]:1

julia> SpecificH(60, 636)
ERROR: DomainError with -150.77268715489907:
log will only return a complex result if called with a complex argument. Try log(Complex(x)).
Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float64)
    @ Base.Math ./math.jl:33
  [2] _log(x::Float64, base::Val{:ℯ}, func::Symbol)
    @ Base.Math ./special/log.jl:304
  [3] log
    @ ./special/log.jl:269 [inlined]
  [4] Region3_ρ(Output::Symbol, ρ::Float64, T::Int64)
    @ SteamTables ~/.julia/packages/SteamTables/pMvX5/src/SteamTables.jl:1718
  [5] f
    @ ~/.julia/packages/SteamTables/pMvX5/src/SteamTables.jl:1756 [inlined]
  [6] Callable_Function
    @ ~/.julia/packages/Roots/bSIXy/src/find_zero.jl:299 [inlined]
  [7] update_state
    @ ~/.julia/packages/Roots/bSIXy/src/derivative_free.jl:117 [inlined]
  [8] update_state
    @ ~/.julia/packages/Roots/bSIXy/src/derivative_free.jl:107 [inlined]
  [9] solve!(𝐙::Roots.ZeroProblemIterator{Roots.Secant, Roots.AlefeldPotraShi, Roots.Callable_Function{Val{1}, Val{false}, SteamTables.var"#f#77"{Int64, Int64}, Nothing}, Roots.UnivariateZeroState{Float64, Float64}, Roots.UnivariateZeroOptions{Float64, Float64, Float64, Float64}, Roots.NullTracks}; verbose::Bool)
    @ Roots ~/.julia/packages/Roots/bSIXy/src/order0.jl:86
 [10] #solve#21
    @ ~/.julia/packages/Roots/bSIXy/src/find_zero.jl:996 [inlined]
 [11] #find_zero#15
    @ ~/.julia/packages/Roots/bSIXy/src/find_zero.jl:709 [inlined]
 [12] find_zero
    @ ~/.julia/packages/Roots/bSIXy/src/find_zero.jl:709 [inlined]
 [13] #find_zero#16
    @ ~/.julia/packages/Roots/bSIXy/src/find_zero.jl:719 [inlined]
 [14] find_zero
    @ ~/.julia/packages/Roots/bSIXy/src/find_zero.jl:719 [inlined]
 [15] Region3(Output::Symbol, P::Int64, T::Int64)
    @ SteamTables ~/.julia/packages/SteamTables/pMvX5/src/SteamTables.jl:1757
 [16] SpecificH(P::Int64, T::Int64)
    @ SteamTables ~/.julia/packages/SteamTables/pMvX5/src/SteamTables.jl:3018
 [17] top-level scope
    @ REPL[5]:1
braamvandyk commented 2 years ago

I shall dig into this. Thanks for logging the issue.

braamvandyk commented 2 years ago

From the first example, the code is looking for a root of f(ρ) = Region3_ρ(:Pressure, ρ, T) - P, but in the range 0..500 for ρ, there are no roots. Looks like a region identification error. Will have to dig a little deeper.

braamvandyk commented 2 years ago

julia> SpecificH(30, 628) 1640.4031042944146

julia> SpecificH(60, 636) 1002.5901301042508

I am admittedly embarrased to report on the source of the error ... In a (no doubt late night) fit of laziness, I initiated the function using the ideal gas law, then did root finding for the actual density. Starting from ideal gas at 30MPa. is where the problem originated, with some help from a density curve that "wobbles" near the region boundary, which is where my search started. Silliness!

I'll update.

braamvandyk commented 2 years ago

Fix (I hope) included in v1.3.0 which is now available.