Merge branch 'master' of https://github.com/ClapeyronThermo/Clapeyron.jl into implicit-solvent

longemen3000 · longemen3000 · commit dabc55b3daf1 · 2026-02-19T20:44:23.000-03:00
diff --git a/NEWS.md b/NEWS.md
@@ -2,6 +2,11 @@
 
 ## New Features
 
+- Implicit differentiation: the core implicit differentiation routines were replaced by `IFTDuals.ift`.
+
 ## package deprecations
 
 ## Bug fixes
+
+- Fixes in Multifluid initial volume
+- Fixes in `MultiPhaseTPFlash`
diff --git a/Project.toml b/Project.toml
@@ -62,7 +62,7 @@ Downloads = "1"
 EoSSuperancillaries = "1.2.3"
 FillArrays = "^0.12, 0.13, 1"
 ForwardDiff = "0.10.36,1"
-IFTDuals = "0.1.2"
+IFTDuals = "0.1.4"
 JSON = "1"
 JutulDarcy = "0.2, 0.3"
 LinearAlgebra = "1"
diff --git a/README.md b/README.md
@@ -24,7 +24,7 @@ We support many equations of state and properties. Some examples of figures you
 
 - Isobaric heat capacity of carbon dioxide at 20 MPa:
 
-  ![CO2_cp](docs/CO2_cp.svg) 
+  ![CO2_cp](docs/CO2_cp.svg)
 
 - Water VLE envelope:
 
@@ -38,7 +38,7 @@ We support many equations of state and properties. Some examples of figures you
 
 ![CH3OH_CyHx](docs/CH3OH_CyHex.svg)
 
-If you find any issues, feel free to contact us directly on the [Zulip Channel](https://julialang.zulipchat.com/#narrow/stream/265161-Clapeyron.2Ejl), or open a Github issue. 
+If you find any issues, feel free to contact us directly on the [Zulip Channel](https://julialang.zulipchat.com/#narrow/stream/265161-Clapeyron.2Ejl), or open a Github issue.
 
 # Installing Clapeyron
 
@@ -93,7 +93,7 @@ uv add pyclapeyron
 
 If you are using Clapeyron for your research work, please cite the following:
 
-```
+```json
 @article{Clapeyron-2022,
     title={Clapeyron.jl: An Extensible, Open-Source Fluid Thermodynamics Toolkit},
     author={Pierre J. Walker, Hon-Wa Yew, and Andrés Riedemann},
@@ -112,4 +112,4 @@ If you are using Clapeyron for your research work, please cite the following:
 
 Note that at its current stage, Clapeyron is still in the early stages of development, and things may be moving around or changing rapidly, but we are very excited to see where this project may go!
 
-We are open to contributions. New models, improved methods and more databases are always appreciated!
+We are open to contributions. New models, improved methods and more databases are always appreciated!
diff --git a/docs/src/installation.md b/docs/src/installation.md
@@ -56,12 +56,12 @@ Here is a list of packages the developers of Clapeyron recommend using:
 Clapeyron.jl is also available in python via the [pyclapeyron](github.com/ClapeyronThermo/pyclapeyron) package.
 You can install `pyclapeyron` from PyPI via `pip` or `uv`:
 
-```
+```sh
 pip install pyclapeyron
 ```
 
 or
 
-```
+```sh
 uv add pyclapeyron
-```
+```
diff --git a/src/methods/pT.jl b/src/methods/pT.jl
@@ -1052,7 +1052,7 @@ export gibbs_energy,helmholtz_energy,gibbs_energy_res,helmholtz_energy_res
 export mass_enthalpy,mass_entropy,mass_internal_energy,mass_gibbs_energy,mass_gibbs_free_energy,mass_helmholtz_energy,mass_helmholtz_free_energy
 #second derivative order properties
 export isochoric_heat_capacity, isobaric_heat_capacity,adiabatic_index
-export mass_isobaric_heat_capacity,mass_isobaric_heat_capacity
+export mass_isobaric_heat_capacity,mass_isochoric_heat_capacity
 export isothermal_compressibility, isentropic_compressibility, speed_of_sound
 export isobaric_expansivity, joule_thomson_coefficient, inversion_temperature
 #higher derivative order properties
diff --git a/src/methods/property_solvers/multicomponent/bubble_point.jl b/src/methods/property_solvers/multicomponent/bubble_point.jl
@@ -525,7 +525,7 @@ function bubble_temperature_init(model,p,x,vol0,T0,y0,volatiles)
 end
 
 """
-    bubble_temperature(model::EoSModel, p, x,method::BubblePointMethod = ChemPotBubbleTemperature())
+    bubble_temperature(model::EoSModel, p, x,method::ThermodynamicMethod = ChemPotBubbleTemperature())
 
 Calculates the bubble temperature and properties at a given pressure `p`.
 Returns a tuple, containing:
diff --git a/src/methods/property_solvers/multicomponent/tp_flash/multiphase.jl b/src/methods/property_solvers/multicomponent/tp_flash/multiphase.jl
@@ -1035,7 +1035,7 @@ function set_idx_vapour!(idx_vapour,model,result)
     v = result.volumes
     β = result.fractions
     x = result.compositions
-    
+
     for i in 1:np
         if is_vapour(identify_phase(model,p,T,x[i],vol = v[i]))
             idx_vapour[] = i
@@ -1171,15 +1171,25 @@ function multi_g_obj(model,p,T,z,_result,ss_cache)
         xnp .= z ./ t
         vnp = vols[np]
         phase_np = __mpflash_phase(idx_vapour[],np)
-        g = modified_gibbs(model,p,T,xnp,phase_np,vnp)[1]
+
         #g = βnp*(eos(model,vnp,T,xnp) + p*vnp)
+        if has_a_res(model)
+            g = βnp*(eos(model,vnp,T,xnp) + p*vnp)
+        else
+            g = βnp*modified_gibbs(model,p,T,xnp,phase_np,vnp)[1]
+        end
+
         for i in 1:np-1
             Ki = viewn(𝕏,nc,i)
             xi .= xnp .* exp.(Ki)
             vi = exp(vols[i])
-            #g += β[i]*(eos(model,vi,T,xi) + p*vi)
             phase_i = __mpflash_phase(idx_vapour[],i)
-            g += modified_gibbs(model,p,T,xi,phase_i,vi)[1]
+            #g += β[i]*(eos(model,vi,T,xi) + p*vi)
+            if has_a_res(model)
+                g += β[i]*(eos(model,vi,T,xi) + p*vi)
+            else
+                g += β[i]*modified_gibbs(model,p,T,xi,phase_i,vi)[1]
+            end
         end
         return g/(Rgas(model)*T)
     end
@@ -1197,7 +1207,7 @@ function initial_beta!(comps,z)
         reconstitute_x!(comps,z,βi,imin)
     end
     βi ./= sum(βi)
-    return βi 
+    return βi
 end
 
 function reconstitute_x!(comps,z,bi,i0)
diff --git a/src/models/CompositeModel/GammaPhi.jl b/src/models/CompositeModel/GammaPhi.jl
@@ -485,7 +485,7 @@ function modified_∂lnϕ∂n(wrapper::PTFlashWrapper{<:GammaPhi}, p, T, z, cach
     end
 end
 
-function ∂lnϕ∂n∂P∂T(wrapper::PTFlashWrapper, p, T, z=SA[1.],cache = ∂lnϕ_cache(model,p,T,z,Val{true}());
+function ∂lnϕ∂n∂P∂T(wrapper::PTFlashWrapper, p, T, z=SA[1.],cache = ∂lnϕ_cache(wrapper,p,T,z,Val{true}());
             phase=:unknown,
             vol0=nothing,
             threaded = true,
@@ -511,7 +511,7 @@ function ∂lnϕ∂n∂P∂T(wrapper::PTFlashWrapper, p, T, z=SA[1.],cache = ∂
     end
 end
 
-function ∂lnϕ∂n∂P(wrapper::PTFlashWrapper, p, T, z=SA[1.],cache = ∂lnϕ_cache(model,p,T,z,Val{false}());
+function ∂lnϕ∂n∂P(wrapper::PTFlashWrapper, p, T, z=SA[1.],cache = ∂lnϕ_cache(wrapper,p,T,z,Val{false}());
             phase=:unknown,
             vol0=nothing,
             threaded = true,
@@ -537,7 +537,7 @@ function ∂lnϕ∂n∂P(wrapper::PTFlashWrapper, p, T, z=SA[1.],cache = ∂lnϕ
     end
 end
 
-function ∂lnϕ∂P(wrapper::PTFlashWrapper, p, T, z=SA[1.], cache = ∂lnϕ_cache(model,p,T,z,Val{false}());
+function ∂lnϕ∂P(wrapper::PTFlashWrapper, p, T, z=SA[1.], cache = ∂lnϕ_cache(wrapper,p,T,z,Val{false}());
             phase=:unknown,
             vol0=nothing,
             threaded = true,
@@ -560,7 +560,7 @@ function ∂lnϕ∂P(wrapper::PTFlashWrapper, p, T, z=SA[1.], cache = ∂lnϕ_ca
     end
 end
 
-function ∂lnϕ∂T(wrapper::PTFlashWrapper, p, T, z=SA[1.], cache = ∂lnϕ_cache(model,p,T,z,Val{false}());
+function ∂lnϕ∂T(wrapper::PTFlashWrapper, p, T, z=SA[1.], cache = ∂lnϕ_cache(wrapper,p,T,z,Val{true}());
             phase=:unknown,
             vol0=nothing,
             threaded = true,
diff --git a/src/models/CompositeModel/LiquidVolumeModel/LiquidVolumeModel.jl b/src/models/CompositeModel/LiquidVolumeModel/LiquidVolumeModel.jl
@@ -1,6 +1,6 @@
 abstract type LiquidVolumeModel <: GibbsBasedModel end
 
-include("NaNLiquid/NaNLiquid.jl")
+include("ZeroLiquid/ZeroLiquid.jl")
 include("Rackett/RackettLiquid.jl")
 include("COSTALD/costald.jl")
 include("DIPPR105Liquid/DIPPR105Liquid.jl")
diff --git a/src/models/CompositeModel/LiquidVolumeModel/ZeroLiquid/ZeroLiquid.jl b/src/models/CompositeModel/LiquidVolumeModel/ZeroLiquid/ZeroLiquid.jl
@@ -25,7 +25,8 @@ liquid = ZeroLiquid(["water","carbon dioxide"])
 ZeroLiquid
 
 export ZeroLiquid
-const NaNLiquid = ZeroLiquid
+
+@deprecate NaNLiquid(args...;kwargs...) ZeroLiquid(args...;kwargs...)
 
 function volume_impl(model::ZeroLiquidModel,p,T,z,phase,threaded,vol0)
     _0 = zero(first(z))
diff --git a/src/models/CompositeModel/SaturationModel/SaturationModel.jl b/src/models/CompositeModel/SaturationModel/SaturationModel.jl
@@ -35,10 +35,10 @@ function saturation_temperature_impl(model::SaturationModel,p,method::Saturation
     return sol,nan,nan
 end
 
-#this method allows to use a Saturation Model as a whole fluid model. it supposes ideal gas and no info about the liquid phase (NaNLiquid)
+#this method allows to use a Saturation Model as a whole fluid model. it supposes ideal gas and no info about the liquid phase (ZeroLiquid)
 function init_puremodel(model::SaturationModel,components,userlocations,verbose)
     _components = format_components(components)
-    fluid = CompositeModel(_components,gas=BasicIdeal(),liquid=NaNLiquid(),saturation = model)
+    fluid = CompositeModel(_components,gas=BasicIdeal(),liquid=ZeroLiquid(),saturation = model)
     return EoSVectorParam(fluid,_components)
 end 
 
diff --git a/src/models/EmpiricHelmholtz/MultiFluid/multifluid.jl b/src/models/EmpiricHelmholtz/MultiFluid/multifluid.jl
@@ -229,9 +229,14 @@ function x0_volume_liquid(model::MultiFluid,p,T,z)
     lb_v = lb_volume(model,T,z)
     for (i,pure) in pairs(model.pures)
         if T > pure.properties.Tc
+            lb_i = lb_volume(pure,T,SA[1.0])
             v0 += z[i]*1.01*lb_volume(pure,T,SA[1.0])
         else
-            v0 += z[i]*x0_volume_liquid(pure,p,T,SA[1.0])
+            if p > x0_psat(pure,T)
+                v0 += z[i]*1.01*lb_volume(pure,T,SA[1.0])
+            else
+                v0 += z[i]*x0_volume_liquid(pure,p,T,SA[1.0])
+            end
         end
     end
     p0 = pressure(model,v0,T,z)
diff --git a/src/models/cubic/alphas/CPAAlpha.jl b/src/models/cubic/alphas/CPAAlpha.jl
@@ -53,7 +53,7 @@ end
 const sCPAAlphaModel = CPAAlphaModel
 
 @newmodelsimple sCPAAlpha CPAAlphaModel CPAAlphaParam
-
+export sCPAAlpha
 """
     sCPAAlpha <: sCPAAlphaModel
     
@@ -91,6 +91,7 @@ alpha = sCPAAlpha(["water","carbon dioxide"];userlocations = (;c1 = [0.67,0.76])
 sCPAAlpha
 default_locations(::Type{sCPAAlpha}) = ["SAFT/CPA/sCPA/sCPA_like.csv"]
 
+#=
 @newmodelsimple eCPAAlpha CPAAlphaModel CPAAlphaParam
 
 """
@@ -128,4 +129,5 @@ alpha = eCPAAlpha(["water","carbon dioxide"];userlocations = (;c1 = [0.67,0.76])
 
 """
 eCPAAlpha
-default_locations(::Type{eCPAAlpha}) = ["SAFT/CPA/eCPA/eCPA_like.csv"]
+default_locations(::Type{eCPAAlpha}) = ["SAFT/CPA/eCPA/eCPA_like.csv"]
+=#
diff --git a/src/models/cubic/alphas/RKAlpha.jl b/src/models/cubic/alphas/RKAlpha.jl
@@ -1,6 +1,7 @@
 abstract type RKAlphaModel <: AlphaModel end
 
 @newmodelsingleton RKAlpha RKAlphaModel
+export RKAlpha
 
 """
     RKAlpha <: RKAlphaModel
@@ -45,4 +46,4 @@ function α_function(model::CubicModel,V,T,z::SingleComp,alpha_model::RKAlphaMod
     Tc = model.params.Tc.values[1]
     Tr = T/Tc
     α = 1 /sqrt(Tr)
-end
+end
diff --git a/src/utils/macros.jl b/src/utils/macros.jl
@@ -1,19 +1,4 @@
 #const IDEALTYPE = Union{T,Type{T}} where T<:EoSModel
-
-"""
-    arbitraryparam(params)
-
-Returns the first field in the struct that is a subtype of `ClapeyronParam`. Errors if it finds none.
-"""
-function arbitraryparam(params)
-    paramstype = typeof(params)
-    idx = findfirst(z->z <: ClapeyronParam,fieldtypes(paramstype))
-    if isnothing(idx)
-        error("The parameter struct ", paramstype, " must contain at least one ClapeyronParam")
-    end
-     return fieldnames(paramstype)[idx] |> z->getfield(params,z)
-end
-
 """
     @groups
 
@@ -622,6 +607,7 @@ the necessary traits to make the model compatible with Clapeyron routines.
     This macro is a no-op from Clapeyron 0.5 onwards.
 """
 macro registermodel(model)
+    Base.depwarn("`@registermodel` is deprecated, it is a no-op.")
     esc(model)
 end
 
diff --git a/test/test_methods_api_flash.jl b/test/test_methods_api_flash.jl
@@ -97,6 +97,11 @@
         res3 = Clapeyron.tp_flash2(system3, p, T, z, MultiPhaseTPFlash())
         @test Clapeyron.numphases(res3) == 3
         @test res3.fractions ≈ [0.3126977407489071, 0.3221079660567595, 0.3651942931943334] rtol = 1e-6
+
+        #issue #546
+        system4 = EOS_CG(["carbon dioxide","water"])
+        res4 = Clapeyron.tp_flash2(system4, 1e6, 293.15, [0.5,0.5], MultiPhaseTPFlash())
+        @test res4.fractions ≈ [0.5025679516689612, 0.4974320483310388] rtol = 1e-6
     end
 
     GC.gc()
diff --git a/test/test_models_others.jl b/test/test_models_others.jl
@@ -260,6 +260,10 @@ end
         dep = departure_functions(system)
         @test count(!iszero,dep) == 1
         @test Clapeyron.eos(system,Vx,Tx,zx) ≈ -6020.0044 rtol = 5e-6
+
+        #546
+        system_546 = EOS_CG(["carbon dioxide","water"])
+        @test volume(system_546,1e6,293.15,[0.1742127426126829, 0.8257872573873171]) ≈ 2.080669615020994e-5 rtol = 1e-6
     end
 
     @testset "LKP" begin
