New struct StateVec and simpler interfaces for phase equilibria

feos-org · g-bauer · Mar 26, 2022 · Mar 21, 2022 · Mar 21, 2022 · Mar 21, 2022
commit aa56a74f63f6ebc54053495371a2355cf66f9150
diff --git a/build_wheel/src/cubic.rs b/build_wheel/src/cubic.rs
@@ -38,16 +38,14 @@ impl_virial_coefficients!(PyPengRobinson);
 
 impl_state!(PengRobinson, PyPengRobinson);
 impl_state_molarweight!(PengRobinson, PyPengRobinson);
-impl_vle_state!(PengRobinson, PyPengRobinson);
+impl_phase_equilibrium!(PengRobinson, PyPengRobinson);
 
 #[pymodule]
 pub fn cubic(_py: Python<'_>, m: &PyModule) -> PyResult<()> {
     m.add_class::<PyPengRobinson>()?;
     m.add_class::<PyPengRobinsonParameters>()?;
     m.add_class::<PyState>()?;
-    m.add_class::<PyPhaseDiagramPure>()?;
-    m.add_class::<PyPhaseDiagramBinary>()?;
-    m.add_class::<PyPhaseDiagramHetero>()?;
+    m.add_class::<PyPhaseDiagram>()?;
     m.add_class::<PyPhaseEquilibrium>()?;
     Ok(())
 }
diff --git a/build_wheel/src/user_defined.rs b/build_wheel/src/user_defined.rs
@@ -1,5 +1,5 @@
-use feos_core::python::user_defined::*;
 use feos_core::python::statehd::*;
+use feos_core::python::user_defined::*;
 use feos_core::*;
 use numpy::convert::ToPyArray;
 use numpy::{PyArray1, PyArray2};
@@ -10,7 +10,6 @@ use quantity::si::*;
 use std::collections::HashMap;
 use std::rc::Rc;
 
-
 /// Equation of state implemented as python class.
 ///
 /// Parameters
@@ -43,7 +42,7 @@ impl_equation_of_state!(PyUserDefinedEos);
 impl_virial_coefficients!(PyUserDefinedEos);
 impl_state!(PyEoSObj, PyUserDefinedEos);
 impl_state_molarweight!(PyEoSObj, PyUserDefinedEos);
-impl_vle_state!(PyEoSObj, PyUserDefinedEos);
+impl_phase_equilibrium!(PyEoSObj, PyUserDefinedEos);
 
 #[pymodule]
 pub fn user_defined(_py: Python, m: &PyModule) -> PyResult<()> {
@@ -61,9 +60,7 @@ pub fn user_defined(_py: Python, m: &PyModule) -> PyResult<()> {
     m.add_class::<PyUserDefinedEos>()?;
     m.add_class::<PyState>()?;
     m.add_class::<PyPhaseEquilibrium>()?;
-    m.add_class::<PyPhaseDiagramPure>()?;
-    m.add_class::<PyPhaseDiagramBinary>()?;
-    m.add_class::<PyPhaseDiagramHetero>()?;
+    m.add_class::<PyPhaseDiagram>()?;
     m.add_class::<PyPhaseEquilibrium>()?;
     Ok(())
 }
diff --git a/example/user_defined_eos.ipynb b/example/user_defined_eos.ipynb
diff --git a/src/errors.rs b/src/errors.rs
@@ -19,9 +19,11 @@ pub enum EosError {
     #[error("Undetermined state: {0}.")]
     UndeterminedState(String),
     #[error("System is supercritical.")]
-    SuperCritical(),
+    SuperCritical,
     #[error("No phase split according to stability analysis.")]
     NoPhaseSplit,
+    #[error("Wrong input units. Expected: {0}, got {1}")]
+    WrongUnits(String, String),
     #[error(transparent)]
     QuantityError(#[from] QuantityError),
     #[error(transparent)]

diff --git a/src/lib.rs b/src/lib.rs
@@ -6,6 +6,7 @@
 use quantity::si::*;
 use quantity::*;
 
+/// Print messages with level `Verbosity::Iter` or higher.
 #[macro_export]
 macro_rules! log_iter {
     ($verbosity:expr, $($arg:tt)*) => {
@@ -15,6 +16,7 @@ macro_rules! log_iter {
     }
 }
 
+/// Print messages with level `Verbosity::Result` or higher.
 #[macro_export]
 macro_rules! log_result {
     ($verbosity:expr, $($arg:tt)*) => {
@@ -38,8 +40,7 @@ pub use equation_of_state::{
 };
 pub use errors::{EosError, EosResult};
 pub use phase_equilibria::{
-    PhaseDiagramBinary, PhaseDiagramHetero, PhaseDiagramPure, PhaseEquilibrium, SolverOptions,
-    Verbosity,
+    PhaseDiagram, PhaseDiagramHetero, PhaseEquilibrium, SolverOptions, StateVec, Verbosity,
 };
 pub use state::{Contributions, DensityInitialization, State, StateBuilder, StateHD};
 

diff --git a/src/phase_equilibria/bubble_dew.rs b/src/phase_equilibria/bubble_dew.rs
@@ -9,6 +9,7 @@ use crate::{equation_of_state::EquationOfState, EosUnit};
 use ndarray::*;
 use num_dual::linalg::{norm, LU};
 use quantity::{QuantityArray1, QuantityScalar};
+use std::convert::TryFrom;
 use std::rc::Rc;
 
 const MAX_ITER_INNER: usize = 5;
@@ -68,36 +69,12 @@ where
 /// # Bubble and dew point calculations
 impl<U: EosUnit, E: EquationOfState> PhaseEquilibrium<U, E, 2> {
     /// Calculate a phase equilibrium for a given temperature
-    /// and composition of the liquid phase.
-    pub fn bubble_point_tx(
+    /// or pressure and composition of the liquid phase.
+    pub fn bubble_point(
         eos: &Rc<E>,
-        temperature: QuantityScalar<U>,
-        pressure: Option<QuantityScalar<U>>,
-        liquid_molefracs: &Array1<f64>,
-        vapor_molefracs: Option<&Array1<f64>>,
-        options: (SolverOptions, SolverOptions),
-    ) -> EosResult<Self>
-    where
-        QuantityScalar<U>: std::fmt::Display,
-    {
-        Self::bubble_dew_point_with_options(
-            eos,
-            TPSpec::Temperature(temperature),
-            pressure,
-            liquid_molefracs,
-            vapor_molefracs,
-            true,
-            options,
-        )
-    }
-
-    /// Calculate a phase equilibrium for a given pressure
-    /// and composition of the liquid phase.
-    pub fn bubble_point_px(
-        eos: &Rc<E>,
-        pressure: QuantityScalar<U>,
-        temperature: Option<QuantityScalar<U>>,
+        temperature_or_pressure: QuantityScalar<U>,
         liquid_molefracs: &Array1<f64>,
+        tp_init: Option<QuantityScalar<U>>,
         vapor_molefracs: Option<&Array1<f64>>,
         options: (SolverOptions, SolverOptions),
     ) -> EosResult<Self>
@@ -106,8 +83,8 @@ impl<U: EosUnit, E: EquationOfState> PhaseEquilibrium<U, E, 2> {
     {
         Self::bubble_dew_point_with_options(
             eos,
-            TPSpec::Pressure(pressure),
-            temperature,
+            TPSpec::try_from(temperature_or_pressure)?,
+            tp_init,
             liquid_molefracs,
             vapor_molefracs,
             true,
@@ -116,36 +93,12 @@ impl<U: EosUnit, E: EquationOfState> PhaseEquilibrium<U, E, 2> {
     }
 
     /// Calculate a phase equilibrium for a given temperature
-    /// and composition of the vapor phase.
-    pub fn dew_point_tx(
+    /// or pressure and composition of the vapor phase.
+    pub fn dew_point(
         eos: &Rc<E>,
-        temperature: QuantityScalar<U>,
-        pressure: Option<QuantityScalar<U>>,
-        vapor_molefracs: &Array1<f64>,
-        liquid_molefracs: Option<&Array1<f64>>,
-        options: (SolverOptions, SolverOptions),
-    ) -> EosResult<Self>
-    where
-        QuantityScalar<U>: std::fmt::Display,
-    {
-        Self::bubble_dew_point_with_options(
-            eos,
-            TPSpec::Temperature(temperature),
-            pressure,
-            vapor_molefracs,
-            liquid_molefracs,
-            false,
-            options,
-        )
-    }
-
-    /// Calculate a phase equilibrium for a given pressure
-    /// and composition of the vapor phase.
-    pub fn dew_point_px(
-        eos: &Rc<E>,
-        pressure: QuantityScalar<U>,
-        temperature: Option<QuantityScalar<U>>,
+        temperature_or_pressure: QuantityScalar<U>,
         vapor_molefracs: &Array1<f64>,
+        tp_init: Option<QuantityScalar<U>>,
         liquid_molefracs: Option<&Array1<f64>>,
         options: (SolverOptions, SolverOptions),
     ) -> EosResult<Self>
@@ -154,8 +107,8 @@ impl<U: EosUnit, E: EquationOfState> PhaseEquilibrium<U, E, 2> {
     {
         Self::bubble_dew_point_with_options(
             eos,
-            TPSpec::Pressure(pressure),
-            temperature,
+            TPSpec::try_from(temperature_or_pressure)?,
+            tp_init,
             vapor_molefracs,
             liquid_molefracs,
             false,

diff --git a/src/phase_equilibria/mod.rs b/src/phase_equilibria/mod.rs
@@ -13,8 +13,8 @@ mod phase_diagram_pure;
 mod stability_analysis;
 mod tp_flash;
 mod vle_pure;
-pub use phase_diagram_binary::{PhaseDiagramBinary, PhaseDiagramHetero};
-pub use phase_diagram_pure::PhaseDiagramPure;
+pub use phase_diagram_binary::PhaseDiagramHetero;
+pub use phase_diagram_pure::{PhaseDiagram, StateVec};
 
 /// Level of detail in the iteration output.
 #[derive(Copy, Clone, PartialOrd, PartialEq)]