Fix state constructor for obscure specification (#26)

prehner · web-flow · commit 2ed2d5c41e58 · 2022-02-11T14:47:57.000+01:00
* Fix state constructor for obscure specification

* update changelog

* no latex in MD on github 😞
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -5,6 +5,8 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased]
+### Fixed
+- Fix state constructor for `T`, `p`, `V`, `x_i` specification. [#26](https://github.com/feos-org/feos-core/pull/26)
 
 ## [0.1.3] - 2022-01-21
 ### Added
diff --git a/src/cubic.rs b/src/cubic.rs
@@ -72,7 +72,7 @@ impl std::fmt::Display for PengRobinsonParameters {
     fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
         self.pure_records
             .iter()
-            .try_for_each(|pr| writeln!(f, "{}", pr.to_string()))?;
+            .try_for_each(|pr| writeln!(f, "{}", pr))?;
         writeln!(f, "\nk_ij:\n{}", self.k_ij)
     }
 }
diff --git a/src/python/cubic.rs b/src/python/cubic.rs
@@ -19,6 +19,14 @@ use std::rc::Rc;
 #[derive(Clone)]
 pub struct PyPengRobinsonRecord(PengRobinsonRecord);
 
+#[pymethods]
+impl PyPengRobinsonRecord {
+    #[new]
+    fn new(tc: f64, pc: f64, acentric_factor: f64) -> Self {
+        Self(PengRobinsonRecord::new(tc, pc, acentric_factor))
+    }
+}
+
 #[pyproto]
 impl pyo3::class::basic::PyObjectProtocol for PyPengRobinsonRecord {
     fn __repr__(&self) -> PyResult<String> {
diff --git a/src/state/mod.rs b/src/state/mod.rs
@@ -355,55 +355,60 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
             )));
         }
         let x = partial_density
-            .map(|pd| pd / pd.sum())
-            .or_else(|| moles.map(|ms| ms / ms.sum()));
+            .map(|pd| pd.to_reduced(pd.sum()))
+            .or_else(|| moles.map(|ms| ms.to_reduced(ms.sum())))
+            .transpose()?;
         let x_u = match (x, molefracs, eos.components()) {
             (Some(_), Some(_), _) => {
                 return Err(EosError::UndeterminedState(String::from(
                     "Composition is overdetermined.",
                 )))
             }
             (Some(x), None, _) => x,
-            (None, Some(x), _) => x.clone().into(),
-            (None, None, 1) => arr1(&[1.0]).into(),
+            (None, Some(x), _) => x.clone(),
+            (None, None, 1) => arr1(&[1.0]),
             _ => {
                 return Err(EosError::UndeterminedState(String::from(
                     "Missing composition.",
                 )))
             }
         };
 
-        // If no extensive property is given, moles is set to reference moles.
-        // Therefore, moles and composition are known from this point on.
-        let n_t = n.unwrap_or_else(U::reference_moles);
-        let n_i = x_u * n_t;
-        let v = volume.or_else(|| rho.map(|d| n_t / d));
+        // If no extensive property is given, moles is set to the reference value.
+        if let (None, None) = (volume, n) {
+            n = Some(U::reference_moles())
+        }
+        let n_i = n.map(|n| &x_u * n);
+        let v = volume.or_else(|| rho.and_then(|d| n.map(|n| n / d)));
 
         // check if new state can be created using default constructor
-        if let (Some(v), Some(t)) = (v, temperature) {
-            return State::new_nvt(eos, t, v, &n_i);
+        if let (Some(v), Some(t), Some(n_i)) = (v, temperature, &n_i) {
+            return State::new_nvt(eos, t, v, n_i);
         }
 
         // Check if new state can be created using density iteration
-        if let (Some(p), Some(t)) = (pressure, temperature) {
-            return State::new_npt(eos, t, p, &n_i, density_initialization);
+        if let (Some(p), Some(t), Some(n_i)) = (pressure, temperature, &n_i) {
+            return State::new_npt(eos, t, p, n_i, density_initialization);
+        }
+        if let (Some(p), Some(t), Some(v)) = (pressure, temperature, v) {
+            return State::new_npvx(eos, t, p, v, &x_u, density_initialization);
         }
 
         // Check if new state can be created using molar_enthalpy and temperature
-        if let (Some(p), Some(h)) = (pressure, molar_enthalpy) {
-            return State::new_nph(eos, p, h, &n_i, density_initialization, initial_temperature);
+        if let (Some(p), Some(h), Some(n_i)) = (pressure, molar_enthalpy, &n_i) {
+            return State::new_nph(eos, p, h, n_i, density_initialization, initial_temperature);
         }
-        if let (Some(p), Some(s)) = (pressure, molar_entropy) {
-            return State::new_nps(eos, p, s, &n_i, density_initialization, initial_temperature);
+        if let (Some(p), Some(s), Some(n_i)) = (pressure, molar_entropy, &n_i) {
+            return State::new_nps(eos, p, s, n_i, density_initialization, initial_temperature);
         }
-        if let (Some(t), Some(h)) = (temperature, molar_enthalpy) {
-            return State::new_nth(eos, t, h, &n_i, density_initialization);
+        if let (Some(t), Some(h), Some(n_i)) = (temperature, molar_enthalpy, &n_i) {
+            return State::new_nth(eos, t, h, n_i, density_initialization);
         }
-        if let (Some(t), Some(s)) = (temperature, molar_entropy) {
-            return State::new_nts(eos, t, s, &n_i, density_initialization);
+        if let (Some(t), Some(s), Some(n_i)) = (temperature, molar_entropy, &n_i) {
+            return State::new_nts(eos, t, s, n_i, density_initialization);
         }
-        if let (Some(u), Some(v)) = (molar_internal_energy, volume) {
-            return State::new_nvu(eos, v, u, &n_i, initial_temperature);
+        if let (Some(u), Some(v), Some(n_i)) = (molar_internal_energy, volume, &n_i) {
+            return State::new_nvu(eos, v, u, n_i, initial_temperature);
         }
         Err(EosError::UndeterminedState(String::from(
             "Missing input parameters.",
@@ -478,6 +483,21 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
         }
     }
 
+    /// Return a new `State` for given pressure $p$, volume $V$, temperature $T$ and composition $x_i$.
+    pub fn new_npvx(
+        eos: &Rc<E>,
+        temperature: QuantityScalar<U>,
+        pressure: QuantityScalar<U>,
+        volume: QuantityScalar<U>,
+        molefracs: &Array1<f64>,
+        density_initialization: DensityInitialization<U>,
+    ) -> EosResult<Self> {
+        let moles = molefracs * U::reference_moles();
+        let state = Self::new_npt(eos, temperature, pressure, &moles, density_initialization)?;
+        let moles = state.partial_density * volume;
+        Self::new_nvt(eos, temperature, volume, &moles)
+    }
+
     /// Return a new `State` for given pressure $p$ and molar enthalpy $h$.
     pub fn new_nph(
         eos: &Rc<E>,
diff --git a/src/utils/estimator.rs b/src/utils/estimator.rs
@@ -115,7 +115,7 @@ where
 impl<U: EosUnit, E: EquationOfState> fmt::Display for Estimator<U, E> {
     fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
         for d in self.data.iter() {
-            writeln!(f, "{}", d.to_string())?;
+            writeln!(f, "{}", d)?;
         }
         Ok(())
     }

Original file line number	Diff line number	Diff line change
`@@ -72,7 +72,7 @@ impl std::fmt::Display for PengRobinsonParameters {`
`72`	`72`	`fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {`
`73`	`73`	`self.pure_records`
`74`	`74`	`.iter()`
`75`		`- .try_for_each(\|pr\| writeln!(f, "{}", pr.to_string()))?;`
	`75`	`+ .try_for_each(\|pr\| writeln!(f, "{}", pr))?;`
`76`	`76`	`writeln!(f, "\nk_ij:\n{}", self.k_ij)`
`77`	`77`	`}`
`78`	`78`	`}`
Original file line number	Diff line number	Diff line change
`@@ -115,7 +115,7 @@ where`
`115`	`115`	`impl<U: EosUnit, E: EquationOfState> fmt::Display for Estimator<U, E> {`
`116`	`116`	`fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {`
`117`	`117`	`for d in self.data.iter() {`
`118`		`- writeln!(f, "{}", d.to_string())?;`
	`118`	`+ writeln!(f, "{}", d)?;`
`119`	`119`	`}`
`120`	`120`	`Ok(())`
`121`	`121`	`}`