include phase fraction in PhaseEquilibrium

feos-org · prehner · Feb 14, 2026 · Jan 27, 2026 · Jan 27, 2026 · Feb 14, 2026
commit 657730565abc6a2ea92625b242b5b3c8ec5307db
diff --git a/crates/feos-core/src/ad/mod.rs b/crates/feos-core/src/ad/mod.rs
@@ -1,6 +1,6 @@
 use crate::DensityInitialization::Liquid;
 use crate::density_iteration::density_iteration;
-use crate::{FeosResult, PhaseEquilibrium, ReferenceSystem, Residual};
+use crate::{Composition, FeosResult, PhaseEquilibrium, ReferenceSystem, Residual};
 use nalgebra::{Const, SVector, U1, U2};
 #[cfg(feature = "rayon")]
 use ndarray::{Array1, Array2, ArrayView2, Zip};
@@ -215,20 +215,21 @@ pub trait PropertiesAD {
         )
     }
 
-    fn bubble_point_pressure<const P: usize>(
+    fn bubble_point_pressure<const P: usize, X: Composition<f64, U2>>(
         &self,
         temperature: Temperature,
         pressure: Option<Pressure>,
-        liquid_molefracs: SVector<f64, 2>,
+        liquid_molefracs: X,
     ) -> FeosResult<Pressure<Gradient<P>>>
     where
         Self: Residual<U2, Gradient<P>>,
     {
         let eos_f64 = self.re();
+        let (liquid_molefracs, _) = liquid_molefracs.into_molefracs(&eos_f64)?;
         let vle = PhaseEquilibrium::bubble_point(
             &eos_f64,
             temperature,
-            &liquid_molefracs,
+            liquid_molefracs,
             pressure,
             None,
             Default::default(),
@@ -265,20 +266,21 @@ pub trait PropertiesAD {
         Ok(Pressure::from_reduced(p))
     }
 
-    fn dew_point_pressure<const P: usize>(
+    fn dew_point_pressure<const P: usize, X: Composition<f64, U2>>(
         &self,
         temperature: Temperature,
         pressure: Option<Pressure>,
-        vapor_molefracs: SVector<f64, 2>,
+        vapor_molefracs: X,
     ) -> FeosResult<Pressure<Gradient<P>>>
     where
         Self: Residual<U2, Gradient<P>>,
     {
         let eos_f64 = self.re();
+        let (vapor_molefracs, _) = vapor_molefracs.into_molefracs(&eos_f64)?;
         let vle = PhaseEquilibrium::dew_point(
             &eos_f64,
             temperature,
-            &vapor_molefracs,
+            vapor_molefracs,
             pressure,
             None,
             Default::default(),
@@ -329,12 +331,8 @@ pub trait PropertiesAD {
             parameters,
             input,
             |eos: &Self::Lifted<Gradient<P>>, inp| {
-                eos.bubble_point_pressure(
-                    inp[0] * KELVIN,
-                    Some(inp[2] * PASCAL),
-                    SVector::from([inp[1], 1.0 - inp[1]]),
-                )
-                .map(|p| p.convert_into(PASCAL))
+                eos.bubble_point_pressure(inp[0] * KELVIN, Some(inp[2] * PASCAL), inp[1])
+                    .map(|p| p.convert_into(PASCAL))
             },
         )
     }
@@ -353,12 +351,8 @@ pub trait PropertiesAD {
             parameters,
             input,
             |eos: &Self::Lifted<Gradient<P>>, inp| {
-                eos.dew_point_pressure(
-                    inp[0] * KELVIN,
-                    Some(inp[2] * PASCAL),
-                    SVector::from([inp[1], 1.0 - inp[1]]),
-                )
-                .map(|p| p.convert_into(PASCAL))
+                eos.dew_point_pressure(inp[0] * KELVIN, Some(inp[2] * PASCAL), inp[1])
+                    .map(|p| p.convert_into(PASCAL))
             },
         )
     }

diff --git a/crates/feos-core/src/errors.rs b/crates/feos-core/src/errors.rs
@@ -24,6 +24,10 @@ pub enum FeosError {
     InvalidState(String, String, f64),
     #[error("Undetermined state: {0}")]
     UndeterminedState(String),
+    #[error(
+        "Extensive properties can only be evaluated for states that are initialized with extensive properties."
+    )]
+    IntensiveState,
     #[error("System is supercritical.")]
     SuperCritical,
     #[error("No phase split according to stability analysis.")]

diff --git a/crates/feos-core/src/lib.rs b/crates/feos-core/src/lib.rs
@@ -299,8 +299,8 @@ mod tests {
 
         // residual properties
         assert_relative_eq!(
-            s.helmholtz_energy(Contributions::Residual),
-            sr.residual_helmholtz_energy(),
+            s.helmholtz_energy(Contributions::Residual)?,
+            sr.residual_helmholtz_energy()?,
             max_relative = 1e-15
         );
         assert_relative_eq!(
@@ -309,8 +309,8 @@ mod tests {
             max_relative = 1e-15
         );
         assert_relative_eq!(
-            s.entropy(Contributions::Residual),
-            sr.residual_entropy(),
+            s.entropy(Contributions::Residual)?,
+            sr.residual_entropy()?,
             max_relative = 1e-15
         );
         assert_relative_eq!(
@@ -319,8 +319,8 @@ mod tests {
             max_relative = 1e-15
         );
         assert_relative_eq!(
-            s.enthalpy(Contributions::Residual),
-            sr.residual_enthalpy(),
+            s.enthalpy(Contributions::Residual)?,
+            sr.residual_enthalpy()?,
             max_relative = 1e-15
         );
         assert_relative_eq!(
@@ -329,8 +329,8 @@ mod tests {
             max_relative = 1e-15
         );
         assert_relative_eq!(
-            s.internal_energy(Contributions::Residual),
-            sr.residual_internal_energy(),
+            s.internal_energy(Contributions::Residual)?,
+            sr.residual_internal_energy()?,
             max_relative = 1e-15
         );
         assert_relative_eq!(
@@ -339,12 +339,12 @@ mod tests {
             max_relative = 1e-15
         );
         assert_relative_eq!(
-            s.gibbs_energy(Contributions::Residual)
-                - s.total_moles()
+            s.gibbs_energy(Contributions::Residual)?
+                - s.total_moles()?
                     * RGAS
                     * s.temperature
                     * s.compressibility(Contributions::Total).ln(),
-            sr.residual_gibbs_energy(),
+            sr.residual_gibbs_energy()?,
             max_relative = 1e-15
         );
         assert_relative_eq!(

diff --git a/crates/feos-core/src/phase_equilibria/bubble_dew.rs b/crates/feos-core/src/phase_equilibria/bubble_dew.rs
@@ -4,7 +4,7 @@ use crate::state::{
     Contributions,
     DensityInitialization::{InitialDensity, Liquid, Vapor},
 };
-use crate::{ReferenceSystem, Residual, SolverOptions, State, Verbosity};
+use crate::{Composition, ReferenceSystem, Residual, SolverOptions, State, Verbosity};
 use nalgebra::allocator::Allocator;
 use nalgebra::{DMatrix, DVector, DefaultAllocator, Dim, Dyn, OVector, U1};
 #[cfg(feature = "ndarray")]
@@ -141,10 +141,10 @@ where
 {
     /// Calculate a phase equilibrium for a given temperature
     /// or pressure and composition of the liquid phase.
-    pub fn bubble_point<TP: TemperatureOrPressure<D>>(
+    pub fn bubble_point<TP: TemperatureOrPressure<D>, X: Composition<D, N>>(
         eos: &E,
         temperature_or_pressure: TP,
-        liquid_molefracs: &OVector<D, N>,
+        liquid_molefracs: X,
         tp_init: Option<TP::Other>,
         vapor_molefracs: Option<&OVector<f64, N>>,
         options: (SolverOptions, SolverOptions),
@@ -162,10 +162,10 @@ where
 
     /// Calculate a phase equilibrium for a given temperature
     /// or pressure and composition of the vapor phase.
-    pub fn dew_point<TP: TemperatureOrPressure<D>>(
+    pub fn dew_point<TP: TemperatureOrPressure<D>, X: Composition<D, N>>(
         eos: &E,
         temperature_or_pressure: TP,
-        vapor_molefracs: &OVector<D, N>,
+        vapor_molefracs: X,
         tp_init: Option<TP::Other>,
         liquid_molefracs: Option<&OVector<f64, N>>,
         options: (SolverOptions, SolverOptions),
@@ -181,35 +181,43 @@ where
         )
     }
 
-    pub(super) fn bubble_dew_point<TP: TemperatureOrPressure<D>>(
+    pub(super) fn bubble_dew_point<TP: TemperatureOrPressure<D>, X: Composition<D, N>>(
         eos: &E,
         temperature_or_pressure: TP,
-        vapor_molefracs: &OVector<D, N>,
+        vapor_molefracs: X,
         tp_init: Option<TP::Other>,
         liquid_molefracs: Option<&OVector<f64, N>>,
         bubble: bool,
         options: (SolverOptions, SolverOptions),
     ) -> FeosResult<Self> {
-        let (temperature, pressure, iterate_p) =
-            temperature_or_pressure.temperature_pressure(tp_init);
-        Self::bubble_dew_point_tp(
-            eos,
-            temperature,
-            pressure,
-            vapor_molefracs,
-            liquid_molefracs,
-            bubble,
-            iterate_p,
-            options,
-        )
+        if eos.components() == 1 {
+            let mut vle = Self::pure(eos, temperature_or_pressure, None, options.1)?;
+            if bubble {
+                vle.phase_fractions = [D::from(0.0), D::from(1.0)];
+            }
+            Ok(vle)
+        } else {
+            let (temperature, pressure, iterate_p) =
+                temperature_or_pressure.temperature_pressure(tp_init);
+            Self::bubble_dew_point_tp(
+                eos,
+                temperature,
+                pressure,
+                vapor_molefracs,
+                liquid_molefracs,
+                bubble,
+                iterate_p,
+                options,
+            )
+        }
     }
 
     #[expect(clippy::too_many_arguments)]
-    fn bubble_dew_point_tp(
+    fn bubble_dew_point_tp<X: Composition<D, N>>(
         eos: &E,
         temperature: Option<Temperature<D>>,
         pressure: Option<Pressure<D>>,
-        molefracs_spec: &OVector<D, N>,
+        composition: X,
         molefracs_init: Option<&OVector<f64, N>>,
         bubble: bool,
         iterate_p: bool,
@@ -218,6 +226,7 @@ where
         let eos_re = eos.re();
         let mut temperature_re = temperature.map(|t| t.re());
         let mut pressure_re = pressure.map(|p| p.re());
+        let (molefracs_spec, total_moles) = composition.into_molefracs(eos)?;
         let molefracs_spec_re = molefracs_spec.map(|x| x.re());
         let (v1, rho2) = if iterate_p {
             // temperature is specified
@@ -306,7 +315,7 @@ where
                 Self::newton_step_t(
                     eos,
                     t,
-                    molefracs_spec,
+                    &molefracs_spec,
                     &mut p,
                     &mut molar_volume,
                     &mut rho2,
@@ -316,7 +325,7 @@ where
                 Self::newton_step_p(
                     eos,
                     &mut t,
-                    molefracs_spec,
+                    &molefracs_spec,
                     p,
                     &mut molar_volume,
                     &mut rho2,
@@ -339,11 +348,11 @@ where
             x2,
         )?;
 
-        Ok(PhaseEquilibrium(if bubble {
-            [state2, state1]
+        Ok(if bubble {
+            PhaseEquilibrium::with_vapor_phase_fraction(state2, state1, D::from(0.0), total_moles)
         } else {
-            [state1, state2]
-        }))
+            PhaseEquilibrium::with_vapor_phase_fraction(state1, state2, D::from(1.0), total_moles)
+        })
     }
 
     fn newton_step_t(