diff --git a/Cargo.toml b/Cargo.toml index ff44725de..2d1af650f 100644 --- 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 index 8f77aa789..9277492fd 100644 --- 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 index 63c668677..f6a992265 100644 --- 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 DataSet for BinaryVleChemicalPotential { fn predict(&self, eos: &Arc) -> Result { 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 index d921ead42..c9bbd883e 100644 --- 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, } impl Diffusion { @@ -19,11 +21,16 @@ impl Diffusion { target: SIArray1, temperature: SIArray1, pressure: SIArray1, + phase: Option<&Vec>, ) -> Result { + 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 DataSet 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 index 4f5894df9..271008a2e 100644 --- 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 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 index c7b77a72e..49542a050 100644 --- 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>, ) -> PyResult { 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>, ) -> PyResult { 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>, ) -> PyResult { 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 index ec7003edd..a53ab8f40 100644 --- 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, } impl ThermalConductivity { @@ -19,11 +21,16 @@ impl ThermalConductivity { target: SIArray1, temperature: SIArray1, pressure: SIArray1, + phase: Option<&Vec>, ) -> Result { + 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 DataSet 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 index ad05b168f..eab7edc0c 100644 --- 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, } impl Viscosity { @@ -19,11 +21,16 @@ impl Viscosity { target: SIArray1, temperature: SIArray1, pressure: SIArray1, + phase: Option<&Vec>, ) -> Result { + 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 DataSet for Viscosity { fn predict(&self, eos: &Arc) -> Result { 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) })