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.

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

Fix detection of real cubic roots #185

Merged
merged 2 commits into from
Jun 22, 2023

Conversation

Liozou
Copy link
Contributor

@Liozou Liozou commented Jun 22, 2023

Thank you for this great package!

I encountered the following behaviour when looking at CO2 with a cubic equation of state:

julia> system = EPPR78(["carbon dioxide"]);

julia> Clapeyron.volume(system , 3311.0u"bar", 400.0u"K")
ERROR: UndefVarError: `vl` not defined
Stacktrace:
 [1] volume_impl(model::PR{…}, p::Float64, T::Float64, z::SVector{…}, phase::Symbol, threaded::Bool, vol0::Nothing)
   @ Clapeyron ~/.julia/packages/Clapeyron/JJGLA/src/models/cubic/equations.jl:269
 [2] volume(model::PR{…}, p::Float64, T::Float64, z::SVector{…}; phase::Symbol, threaded::Bool, vol0::Nothing)
   @ Clapeyron ~/.julia/packages/Clapeyron/JJGLA/src/methods/property_solvers/volume.jl:118 [inlined]
 [3] volume(model::PR{…}, p::Quantity{…}, T::Quantity{…}, z::SVector{…}; phase::Symbol, output::Unitful.FreeUnits{…})
   @ Clapeyron ~/.julia/packages/Clapeyron/JJGLA/src/methods/unitful_methods.jl:58 [inlined]
 [4] volume(model::PR{…}, p::Quantity{…}, T::Quantity{…}, z::SVector{…})
   @ Clapeyron ~/.julia/packages/Clapeyron/JJGLA/src/methods/unitful_methods.jl:55 [inlined]
 [5] top-level scope
   @ REPL[3]:1
Some type information was truncated. Use `show(err)` to see complete types.

Arguably worse, for a slightly different pressure, there is no error but the answer is obviously wrong (the volume is negative):

julia> Clapeyron.volume(system , 3363.1u"bar", 400.0u"K")
-5.070066466695018e-5 m^3

After investigating, it turns out that this is due to a faulty heuristic used to determine which roots of a cubic polynomial are real. This PR implements a cleaner heuristic that checks the value of the polynomial between the real parts of the roots, and compares the sign of these values to determine more precisely whether the roots are real. This yields more satisfying answers:

julia> Clapeyron.volume(system , 3311.0u"bar", 400.0u"K")
3.363139140634349e-5 m^3

julia> Clapeyron.volume(system , 3363.1u"bar", 400.0u"K")
3.3545657072357283e-5 m^3

And as a bonus, it seems that this bugfix incidentally makes the test for #172 pass (whereas it was marked as @test_broken).

@Liozou
Copy link
Contributor Author

Liozou commented Jun 22, 2023

Nevermind the "bonus", it seems that CI is a bit flaky considering the last runs. Should I erase the @test_broken into @test conversion?

@longemen3000
Copy link
Member

longemen3000 commented Jun 22, 2023

Nice!
Just comment that test out for now. Can you update the version number?, so we can release a new version as soon this is merged

@Liozou
Copy link
Contributor Author

Liozou commented Jun 22, 2023

Sure, I just did that (and added a small comment by the test to explain what the issue is)!

@codecov
Copy link

codecov bot commented Jun 22, 2023

Codecov Report

Merging #185 (76f4769) into master (45724c3) will decrease coverage by 0.01%.
The diff coverage is 100.00%.

❗ Current head 76f4769 differs from pull request most recent head 3b9eda7. Consider uploading reports for the commit 3b9eda7 to get more accurate results

@@            Coverage Diff             @@
##           master     #185      +/-   ##
==========================================
- Coverage   85.76%   85.76%   -0.01%     
==========================================
  Files         222      222              
  Lines       15099    15097       -2     
==========================================
- Hits        12950    12948       -2     
  Misses       2149     2149              
Impacted Files Coverage Δ
src/models/cubic/equations.jl 85.27% <100.00%> (+0.59%) ⬆️
src/modules/solvers/poly.jl 100.00% <100.00%> (ø)

... and 11 files with indirect coverage changes

@longemen3000 longemen3000 merged commit 86b8464 into ClapeyronThermo:master Jun 22, 2023
10 checks passed
@Liozou Liozou deleted the cubicrealroots branch June 22, 2023 18:57
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants