Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

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

Open
tkeskita opened this issue Apr 17, 2023 · 3 comments
Open

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

tkeskita opened this issue Apr 17, 2023 · 3 comments

Comments

@tkeskita
Copy link
Contributor

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
Copy link
Member

longemen3000 commented Apr 17, 2023 •
edited

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
Copy link
Member

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
Copy link
Member

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

@Liozou Liozou mentioned this issue Jun 22, 2023
Merged
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@tkeskita @longemen3000