feos-org · g-bauer · Jul 17, 2023 · Jul 14, 2023 · Jul 17, 2023
diff --git a/Cargo.toml b/Cargo.toml
@@ -38,6 +38,7 @@ serde_json = "1.0"
 lazy_static = { version = "1.4", optional = true }
 indexmap = "1.8"
 rayon = { version = "1.5", optional = true }
+itertools = "0.10"
 
 [dependencies.pyo3]
 version = "0.18"

diff --git a/feos-core/CHANGELOG.md b/feos-core/CHANGELOG.md
@@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added `Components`, `Residual`, `IdealGas` and `DeBroglieWavelength` traits to decouple ideal gas models from residual models. [#158](https://github.com/feos-org/feos/pull/158)
 - Added `JobackParameters` struct that implements `Parameters` including Python bindings. [#158](https://github.com/feos-org/feos/pull/158)
 - Added `Parameter::from_model_records` as a simpler interface to generate parameters.  [#169](https://github.com/feos-org/feos/pull/169)
+- Added optional `Phase` argument for constructors of dynamic properties of `DataSet`s.  [#174](https://github.com/feos-org/feos/pull/174)
 
 ### Changed
 - Changed constructors of `Parameter` trait to return `Result`s. [#161](https://github.com/feos-org/feos/pull/161)

diff --git a/src/estimator/binary_vle.rs b/src/estimator/binary_vle.rs
@@ -1,20 +1,13 @@
-use super::{DataSet, EstimatorError};
+use super::{DataSet, EstimatorError, Phase};
 use feos_core::{
     Contributions, DensityInitialization, EosUnit, PhaseDiagram, PhaseEquilibrium, Residual, State,
 };
+use itertools::izip;
 use ndarray::{arr1, s, Array1, ArrayView1, Axis};
 use quantity::si::{SIArray1, SINumber, SIUnit};
 use std::collections::HashMap;
 use std::sync::Arc;
 
-/// Different phases of experimental data points in the `BinaryVlePressure` data set.
-#[derive(Clone, Copy)]
-#[cfg_attr(feature = "python", pyo3::pyclass)]
-pub enum Phase {
-    Vapor,
-    Liquid,
-}
-
 /// Store experimental binary VLE data for the calculation of chemical potential residuals.
 #[derive(Clone)]
 pub struct BinaryVleChemicalPotential {
@@ -63,13 +56,12 @@ impl<E: Residual> DataSet<E> for BinaryVleChemicalPotential {
 
     fn predict(&self, eos: &Arc<E>) -> Result<SIArray1, EstimatorError> {
         let mut prediction = Vec::new();
-        for (((&xi, &yi), t), p) in self
-            .liquid_molefracs
-            .iter()
-            .zip(self.vapor_molefracs.iter())
-            .zip(self.temperature.into_iter())
-            .zip(self.pressure.into_iter())
-        {
+        for (&xi, &yi, t, p) in izip!(
+            &self.liquid_molefracs,
+            &self.vapor_molefracs,
+            &self.temperature,
+            &self.pressure
+        ) {
             let liquid_moles = arr1(&[xi, 1.0 - xi]) * SIUnit::reference_moles();
             let liquid = State::new_npt(eos, t, p, &liquid_moles, DensityInitialization::Liquid)?;
             let mu_res_liquid = liquid.residual_chemical_potential();

diff --git a/src/estimator/diffusion.rs b/src/estimator/diffusion.rs
@@ -1,5 +1,6 @@
-use super::{DataSet, EstimatorError};
-use feos_core::{DensityInitialization, EntropyScaling, EosUnit, Residual, State};
+use super::{DataSet, EstimatorError, Phase};
+use feos_core::{DensityInitialization, EntropyScaling, EosUnit, State, Residual};
+use itertools::izip;
 use ndarray::{arr1, Array1};
 use quantity::si::{SIArray1, SIUnit};
 use std::collections::HashMap;
@@ -11,6 +12,7 @@ pub struct Diffusion {
     pub target: SIArray1,
     temperature: SIArray1,
     pressure: SIArray1,
+    initial_density: Vec<DensityInitialization>,
 }
 
 impl Diffusion {
@@ -19,11 +21,16 @@ impl Diffusion {
         target: SIArray1,
         temperature: SIArray1,
         pressure: SIArray1,
+        phase: Option<&Vec<Phase>>,
     ) -> Result<Self, EstimatorError> {
+        let n = temperature.len();
         Ok(Self {
             target,
             temperature,
             pressure,
+            initial_density: phase.map_or(vec![DensityInitialization::None; n], |phase| {
+                phase.iter().map(|&p| p.into()).collect()
+            }),
         })
     }
 
@@ -62,16 +69,14 @@ impl<E: Residual + EntropyScaling> DataSet<E> for Diffusion {
             .to_reduced(SIUnit::reference_pressure())
             .unwrap();
 
-        let res = ts
-            .iter()
-            .zip(ps.iter())
-            .map(|(&t, &p)| {
+        let res = izip!(&ts, &ps, &self.initial_density)
+            .map(|(&t, &p, &initial_density)| {
                 State::new_npt(
                     eos,
                     t * SIUnit::reference_temperature(),
                     p * SIUnit::reference_pressure(),
                     &moles,
-                    DensityInitialization::None,
+                    initial_density,
                 )?
                 .diffusion()?
                 .to_reduced(SIUnit::reference_diffusion())

diff --git a/src/estimator/mod.rs b/src/estimator/mod.rs
@@ -1,5 +1,5 @@
 //! Utilities for working with experimental data.
-use feos_core::EosError;
+use feos_core::{DensityInitialization, EosError};
 use quantity::QuantityError;
 use std::num::ParseFloatError;
 use thiserror::Error;
@@ -17,7 +17,7 @@ pub use vapor_pressure::VaporPressure;
 mod liquid_density;
 pub use liquid_density::{EquilibriumLiquidDensity, LiquidDensity};
 mod binary_vle;
-pub use binary_vle::{BinaryPhaseDiagram, BinaryVleChemicalPotential, BinaryVlePressure, Phase};
+pub use binary_vle::{BinaryPhaseDiagram, BinaryVleChemicalPotential, BinaryVlePressure};
 mod viscosity;
 pub use viscosity::Viscosity;
 mod thermal_conductivity;
@@ -28,6 +28,23 @@ pub use diffusion::Diffusion;
 #[cfg(feature = "python")]
 pub mod python;
 
+/// Different phases of experimental data points.
+#[derive(Clone, Copy)]
+#[cfg_attr(feature = "python", pyo3::pyclass)]
+pub enum Phase {
+    Vapor,
+    Liquid,
+}
+
+impl From<Phase> for DensityInitialization {
+    fn from(value: Phase) -> Self {
+        match value {
+            Phase::Liquid => DensityInitialization::Liquid,
+            Phase::Vapor => DensityInitialization::Vapor,
+        }
+    }
+}
+
 #[derive(Debug, Error)]
 pub enum EstimatorError {
     #[error("Input has not the same amount of data as the target.")]

diff --git a/src/estimator/python.rs b/src/estimator/python.rs
@@ -639,21 +639,26 @@ macro_rules! impl_estimator_entropy_scaling {
             ///     Temperature for experimental data points.
             /// pressure : SIArray1
             ///     Pressure for experimental data points.
+            /// phase : List[Phase], optional
+            ///     Phase of data. Used to determine the starting
+            ///     density for the density iteration. If provided,
+            ///     resulting states may not be stable.
             ///
             /// Returns
             /// -------
             /// DataSet
             #[staticmethod]
-            #[pyo3(text_signature = "(target, temperature, pressure)")]
             fn viscosity(
                 target: &PySIArray1,
                 temperature: &PySIArray1,
                 pressure: &PySIArray1,
+                phase: Option<Vec<Phase>>,
             ) -> PyResult<Self> {
                 Ok(Self(Arc::new(Viscosity::new(
                     target.clone().into(),
                     temperature.clone().into(),
                     pressure.clone().into(),
+                    phase.as_ref(),
                 )?)))
             }
 
@@ -667,21 +672,26 @@ macro_rules! impl_estimator_entropy_scaling {
             ///     Temperature for experimental data points.
             /// pressure : SIArray1
             ///     Pressure for experimental data points.
+            /// phase : List[Phase], optional
+            ///     Phase of data. Used to determine the starting
+            ///     density for the density iteration. If provided,
+            ///     resulting states may not be stable.
             ///
             /// Returns
             /// -------
             /// DataSet
             #[staticmethod]
-            #[pyo3(text_signature = "(target, temperature, pressure)")]
             fn thermal_conductivity(
                 target: &PySIArray1,
                 temperature: &PySIArray1,
                 pressure: &PySIArray1,
+                phase: Option<Vec<Phase>>,
             ) -> PyResult<Self> {
                 Ok(Self(Arc::new(ThermalConductivity::new(
                     target.clone().into(),
                     temperature.clone().into(),
                     pressure.clone().into(),
+                    phase.as_ref(),
                 )?)))
             }
 
@@ -695,21 +705,26 @@ macro_rules! impl_estimator_entropy_scaling {
             ///     Temperature for experimental data points.
             /// pressure : SIArray1
             ///     Pressure for experimental data points.
+            /// phase : List[Phase], optional
+            ///     Phase of data. Used to determine the starting
+            ///     density for the density iteration. If provided,
+            ///     resulting states may not be stable.
             ///
             /// Returns
             /// -------
             /// DataSet
             #[staticmethod]
-            #[pyo3(text_signature = "(target, temperature, pressure)")]
             fn diffusion(
                 target: &PySIArray1,
                 temperature: &PySIArray1,
                 pressure: &PySIArray1,
+                phase: Option<Vec<Phase>>,
             ) -> PyResult<Self> {
                 Ok(Self(Arc::new(Diffusion::new(
                     target.clone().into(),
                     temperature.clone().into(),
                     pressure.clone().into(),
+                    phase.as_ref(),
                 )?)))
             }
         }

diff --git a/src/estimator/thermal_conductivity.rs b/src/estimator/thermal_conductivity.rs
@@ -1,5 +1,6 @@
-use super::{DataSet, EstimatorError};
-use feos_core::{DensityInitialization, EntropyScaling, EosUnit, Residual, State};
+use super::{DataSet, EstimatorError, Phase};
+use feos_core::{DensityInitialization, EntropyScaling, EosUnit, State, Residual};
+use itertools::izip;
 use ndarray::{arr1, Array1};
 use quantity::si::{SIArray1, SIUnit};
 use std::collections::HashMap;
@@ -11,6 +12,7 @@ pub struct ThermalConductivity {
     pub target: SIArray1,
     temperature: SIArray1,
     pressure: SIArray1,
+    initial_density: Vec<DensityInitialization>,
 }
 
 impl ThermalConductivity {
@@ -19,11 +21,16 @@ impl ThermalConductivity {
         target: SIArray1,
         temperature: SIArray1,
         pressure: SIArray1,
+        phase: Option<&Vec<Phase>>,
     ) -> Result<Self, EstimatorError> {
+        let n = temperature.len();
         Ok(Self {
             target,
             temperature,
             pressure,
+            initial_density: phase.map_or(vec![DensityInitialization::None; n], |phase| {
+                phase.iter().map(|&p| p.into()).collect()
+            }),
         })
     }
 
@@ -66,16 +73,14 @@ impl<E: Residual + EntropyScaling> DataSet<E> for ThermalConductivity {
             / SIUnit::reference_temperature()
             / SIUnit::reference_length();
 
-        let res = ts
-            .iter()
-            .zip(ps.iter())
-            .map(|(&t, &p)| {
+        let res = izip!(&ts, &ps, &self.initial_density)
+            .map(|(&t, &p, &initial_density)| {
                 State::new_npt(
                     eos,
                     t * SIUnit::reference_temperature(),
                     p * SIUnit::reference_pressure(),
                     &moles,
-                    DensityInitialization::None,
+                    initial_density,
                 )?
                 .thermal_conductivity()?
                 .to_reduced(unit)

diff --git a/src/estimator/viscosity.rs b/src/estimator/viscosity.rs
@@ -1,5 +1,6 @@
-use super::{DataSet, EstimatorError};
-use feos_core::{DensityInitialization, EntropyScaling, EosUnit, Residual, State};
+use super::{DataSet, EstimatorError, Phase};
+use feos_core::{DensityInitialization, EntropyScaling, EosUnit, State, Residual};
+use itertools::izip;
 use ndarray::arr1;
 use quantity::si::{SIArray1, SIUnit};
 use std::collections::HashMap;
@@ -11,6 +12,7 @@ pub struct Viscosity {
     pub target: SIArray1,
     temperature: SIArray1,
     pressure: SIArray1,
+    initial_density: Vec<DensityInitialization>,
 }
 
 impl Viscosity {
@@ -19,11 +21,16 @@ impl Viscosity {
         target: SIArray1,
         temperature: SIArray1,
         pressure: SIArray1,
+        phase: Option<&Vec<Phase>>,
     ) -> Result<Self, EstimatorError> {
+        let n = temperature.len();
         Ok(Self {
             target,
             temperature,
             pressure,
+            initial_density: phase.map_or(vec![DensityInitialization::None; n], |phase| {
+                phase.iter().map(|&p| p.into()).collect()
+            }),
         })
     }
 
@@ -53,11 +60,9 @@ impl<E: Residual + EntropyScaling> DataSet<E> for Viscosity {
 
     fn predict(&self, eos: &Arc<E>) -> Result<SIArray1, EstimatorError> {
         let moles = arr1(&[1.0]) * SIUnit::reference_moles();
-        self.temperature
-            .into_iter()
-            .zip(self.pressure.into_iter())
-            .map(|(t, p)| {
-                State::new_npt(eos, t, p, &moles, DensityInitialization::None)?
+        izip!(&self.temperature, &self.pressure, &self.initial_density)
+            .map(|(t, p, &initial_density)| {
+                State::new_npt(eos, t, p, &moles, initial_density)?
                     .viscosity()
                     .map_err(EstimatorError::from)
             })