search

ClapeyronThermo / Clapeyron.jl

Clapeyron provides a framework for the development and use of fluid-thermodynamic models, including SAFT, cubic, activity, multi-parameter, and COSMO-SAC.
MIT License
188 stars 47 forks source link

LinearAlgebra.SingularException in examples/multi-component_vle_vlle_lle_crit notebook #172

Open tkeskita opened 1 year ago

tkeskita commented 1 year ago

Hi,

running this notebook stops at 3rd input cell with error:

LinearAlgebra.SingularException(5)

Stacktrace: [1] checknonsingular @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/factorization.jl:19 [inlined] [2] checknonsingular @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/factorization.jl:21 [inlined] [3] #lu!#136 @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:85 [inlined] [4] #lu#140 @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:273 [inlined] [5] lu (repeats 2 times) @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:272 [inlined] [6] (A::Matrix{Float64}, B::Vector{Float64}) @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/generic.jl:1136 [7] default_newton_linsolve(d::Vector{Float64}, B::Matrix{Float64}, g::Vector{Float64}) @ NLSolvers ~/.julia/packages/NLSolvers/lRxce/src/quasinewton/approximations/newton.jl:15 [8] solve(problem::NLSolvers.NEqProblem{NLSolvers.VectorObjective{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Clapeyron.Solvers.var"#j!#1"{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, ForwardDiff.JacobianConfig{ForwardDiff.Tag{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Float64}, Float64, 2, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Float64}, Float64, 2}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Float64}, Float64, 2}}}}, Vector{Float64}}, Clapeyron.Solvers.var"#fj!#2"{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, ForwardDiff.JacobianConfig{ForwardDiff.Tag{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Float64}, Float64, 2, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Float64}, Float64, 2}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Float64}, Float64, 2}}}}}, Clapeyron.Solvers.var"#jv!#3"}, Nothing, NLSolvers.Euclidean{Tuple{0}}, NLSolvers.InPlace}, x::Vector{Float64}, method::NLSolvers.LineSearch{NLSolvers.Newton{NLSolvers.Direct, typeof(NLSolvers.default_newton_linsolve), Nothing, Nothing}, NLSolvers.Backtracking{Float64, Int64, NLSolvers.FixedInterp, NamedTuple{(:lower, :upper), Tuple{Int64, Float64}}}, NLSolvers.InitialScaling{NLSolvers.ShannoPhua}}, options::NLSolvers.NEqOptions{Float64, Int64, Nothing}, state::NamedTuple{(:z, :d, :Fx, :Jx), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Matrix{Float64}}}) @ NLSolvers ~/.julia/packages/NLSolvers/lRxce/src/nlsolve/linesearch/newton.jl:71 [9] solve(problem::NLSolvers.NEqProblem{NLSolvers.VectorObjective{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Clapeyron.Solvers.var"#j!#1"{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, ForwardDiff.JacobianConfig{ForwardDiff.Tag{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Float64}, Float64, 2, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Float64}, Float64, 2}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Float64}, Float64, 2}}}}, Vector{Float64}}, Clapeyron.Solvers.var"#fj!#2"{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, ForwardDiff.JacobianConfig{ForwardDiff.Tag{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Float64}, Float64, 2, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Float64}, Float64, 2}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Float64}, Float64, 2}}}}}, Clapeyron.Solvers.var"#jv!#3"}, Nothing, NLSolvers.Euclidean{Tuple{0}}, NLSolvers.InPlace}, x::Vector{Float64}, method::NLSolvers.LineSearch{NLSolvers.Newton{NLSolvers.Direct, typeof(NLSolvers.default_newton_linsolve), Nothing, Nothing}, NLSolvers.Backtracking{Float64, Int64, NLSolvers.FixedInterp, NamedTuple{(:lower, :upper), Tuple{Int64, Float64}}}, NLSolvers.InitialScaling{NLSolvers.ShannoPhua}}, options::NLSolvers.NEqOptions{Float64, Int64, Nothing}) @ NLSolvers ~/.julia/packages/NLSolvers/lRxce/src/nlsolve/linesearch/newton.jl:11 [10] nlsolve(nl_problem::NLSolvers.NEqProblem{NLSolvers.VectorObjective{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Clapeyron.Solvers.var"#j!#1"{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, ForwardDiff.JacobianConfig{ForwardDiff.Tag{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Float64}, Float64, 2, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Float64}, Float64, 2}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Float64}, Float64, 2}}}}, Vector{Float64}}, Clapeyron.Solvers.var"#fj!#2"{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, ForwardDiff.JacobianConfig{ForwardDiff.Tag{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Float64}, Float64, 2, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Float64}, Float64, 2}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Clapeyron.var"#f!#261"{PCSAFT{BasicIdeal}, Float64, Vector{Float64}, Float64}, Float64}, Float64, 2}}}}}, Clapeyron.Solvers.var"#jv!#3"}, Nothing, NLSolvers.Euclidean{Tuple{0}}, NLSolvers.InPlace}, x0::Vector{Float64}, method::NLSolvers.LineSearch{NLSolvers.Newton{NLSolvers.Direct, typeof(NLSolvers.default_newton_linsolve), Nothing, Nothing}, NLSolvers.Backtracking{Float64, Int64, NLSolvers.FixedInterp, NamedTuple{(:lower, :upper), Tuple{Int64, Float64}}}, NLSolvers.InitialScaling{NLSolvers.ShannoPhua}}, options::NLSolvers.NEqOptions{Float64, Int64, Nothing}) @ Clapeyron.Solvers ~/.julia/packages/Clapeyron/kGA2x/src/solvers/nlsolve.jl:22 [11] nlsolve(f!::Function, x0::Vector{Float64}, method::NLSolvers.LineSearch{NLSolvers.Newton{NLSolvers.Direct, typeof(NLSolvers.default_newton_linsolve), Nothing, Nothing}, NLSolvers.Backtracking{Float64, Int64, NLSolvers.FixedInterp, NamedTuple{(:lower, :upper), Tuple{Int64, Float64}}}, NLSolvers.InitialScaling{NLSolvers.ShannoPhua}}, options::NLSolvers.NEqOptions{Float64, Int64, Nothing}, chunk::ForwardDiff.Chunk{2}) @ Clapeyron.Solvers ~/.julia/packages/Clapeyron/kGA2x/src/solvers/nlsolve.jl:18 [12] nlsolve(f!::Function, x0::Vector{Float64}, method::NLSolvers.LineSearch{NLSolvers.Newton{NLSolvers.Direct, typeof(NLSolvers.default_newton_linsolve), Nothing, Nothing}, NLSolvers.Backtracking{Float64, Int64, NLSolvers.FixedInterp, NamedTuple{(:lower, :upper), Tuple{Int64, Float64}}}, NLSolvers.InitialScaling{NLSolvers.ShannoPhua}}, options::NLSolvers.NEqOptions{Float64, Int64, Nothing}) @ Clapeyron.Solvers ~/.julia/packages/Clapeyron/kGA2x/src/solvers/nlsolve.jl:16 [13] bubble_pressure_impl(model::PCSAFT{BasicIdeal}, T::Float64, x::Vector{Float64}, method::ChemPotBubblePressure{Float64}) @ Clapeyron ~/.julia/packages/Clapeyron/kGA2x/src/methods/property_solvers/multicomponent/bubble_point/bubble_chempot.jl:83 [14] bubble_pressure(model::PCSAFT{BasicIdeal}, T::Float64, x::Vector{Float64}, method::ChemPotBubblePressure{Float64}) @ Clapeyron ~/.julia/packages/Clapeyron/kGA2x/src/methods/property_solvers/multicomponent/bubble_point.jl:138 [15] #bubble_pressure#274 @ ~/.julia/packages/Clapeyron/kGA2x/src/methods/property_solvers/multicomponent/bubble_point.jl:321 [inlined] [16] top-level scope @ In[3]:21

longemen3000 commented 1 year ago

reproducer:

function error_172_reproducer()
    x = [0.96611,0.01475,0.01527,0.00385]
    T = 202.694
    v0 = [-4.136285855713797, -4.131888756537859, 0.9673991775701574, 0.014192499147585259, 0.014746430039492817, 0.003661893242764558]
    model = PCSAFT(["methane","butane","isobutane","pentane"])
    bubble_pressure(model,T,x;v0 = v0)
end

it is important to note that the point that fail is the last one (the critical temp). it is inherently hard to solver at those conditions. what happened probably is that the phases merged into one, i suppose?

if we use crit_mix:

julia> crit_mix(model3,[0.96611,0.01475,0.01527,0.00385])
(202.66602773754164, 5.908687689953082e6, 7.316560337045216e-5)

seems that the critical point for that particular model is lower than the one anotated

longemen3000 commented 1 year ago

this particular case was added as a test. by a refactoring on how bubble and dew pressure points are calculated, this specific error does not appear anymore with this combination of inputs

longemen3000 commented 1 year ago

the real reproducer:

B = [-1.4519728743587657 1.451736284830506 -5.934368273445117 -1.0154003713150401 -1.4228901646890824; 5.066334307448492 -5.0664684779036735 -12.664926271615439 -64.5322835331547 -3.47142940880429; 4.537903653187883 -4.538068926327418 -12.56861224579578 -2.823162310529271 -64.73486778751243; 6.566654980630976 -6.566792578100933 212.82034111192326 224.37059305703244 223.668115325478; -0.47703263656109995 0.4769349928972182 -2.0794364795445763 -0.46597769049243537 -0.6415589206816118]

g = [1.9713002873173745e-9, -6.777883544599869e-9, -5.328619868058899e-9, -1.1566711761221415e-8, 6.496464364696421e-10]
B \ g #errors