Reimplement Molarweight trait (#258)

feos-org · Oct 30, 2024 · 950a96e · 950a96e
1 parent 576b56b
commit 950a96e
Show file tree

Hide file tree

Showing 28 changed files with 238 additions and 164 deletions.
diff --git a/feos-core/src/cubic.rs b/feos-core/src/cubic.rs
@@ -4,7 +4,7 @@
 //! of state - with a single contribution to the Helmholtz energy - can be implemented.
 //! The implementation closely follows the form of the equations given in
 //! [this wikipedia article](https://en.wikipedia.org/wiki/Cubic_equations_of_state#Peng%E2%80%93Robinson_equation_of_state).
-use crate::equation_of_state::{Components, Residual};
+use crate::equation_of_state::{Components, Molarweight, Residual};
 use crate::parameter::{Identifier, Parameter, ParameterError, PureRecord};
 use crate::state::StateHD;
 use ndarray::{Array1, Array2, ScalarOperand};
@@ -214,7 +214,9 @@ impl Residual for PengRobinson {
             self.residual_helmholtz_energy(state),
         )]
     }
+}
 
+impl Molarweight for PengRobinson {
     fn molar_weight(&self) -> MolarWeight<Array1<f64>> {
         &self.parameters.molarweight * (GRAM / MOL)
     }

diff --git a/feos-core/src/equation_of_state/mod.rs b/feos-core/src/equation_of_state/mod.rs
@@ -9,7 +9,7 @@ mod ideal_gas;
 mod residual;
 
 pub use ideal_gas::IdealGas;
-pub use residual::{EntropyScaling, NoResidual, Residual};
+pub use residual::{EntropyScaling, Molarweight, NoResidual, Residual};
 
 /// The number of components that the model is initialized for.
 pub trait Components {
@@ -91,7 +91,9 @@ impl<I: IdealGas, R: Residual> Residual for EquationOfState<I, R> {
     ) -> Vec<(String, D)> {
         self.residual.residual_helmholtz_energy_contributions(state)
     }
+}
 
+impl<I, R: Molarweight> Molarweight for EquationOfState<I, R> {
     fn molar_weight(&self) -> MolarWeight<Array1<f64>> {
         self.residual.molar_weight()
     }

diff --git a/feos-core/src/equation_of_state/residual.rs b/feos-core/src/equation_of_state/residual.rs
@@ -8,7 +8,14 @@ use quantity::*;
 use std::ops::Div;
 use typenum::Quot;
 
-/// A reisdual Helmholtz energy model.
+/// Molar weight of all components.
+///
+/// Enables calculation of (mass) specific properties.
+pub trait Molarweight {
+    fn molar_weight(&self) -> MolarWeight<Array1<f64>>;
+}
+
+/// A residual Helmholtz energy model.
 pub trait Residual: Components + Send + Sync {
     /// Return the maximum density in Angstrom^-3.
     ///
@@ -18,11 +25,6 @@ pub trait Residual: Components + Send + Sync {
     /// equation of state anyways).
     fn compute_max_density(&self, moles: &Array1<f64>) -> f64;
 
-    /// Molar weight of all components.
-    ///
-    /// Enables calculation of (mass) specific properties.
-    fn molar_weight(&self) -> MolarWeight<Array1<f64>>;
-
     /// Evaluate the reduced Helmholtz energy of each individual contribution
     /// and return them together with a string representation of the contribution.
     fn residual_helmholtz_energy_contributions<D: DualNum<f64> + Copy + ScalarOperand>(
@@ -193,8 +195,4 @@ impl Residual for NoResidual {
     ) -> Vec<(String, D)> {
         vec![]
     }
-
-    fn molar_weight(&self) -> MolarWeight<Array1<f64>> {
-        panic!("No mass specific properties are available for this model!")
-    }
 }
diff --git a/feos-core/src/lib.rs b/feos-core/src/lib.rs
@@ -33,7 +33,7 @@ pub mod parameter;
 mod phase_equilibria;
 mod state;
 pub use equation_of_state::{
-    Components, EntropyScaling, EquationOfState, IdealGas, NoResidual, Residual,
+    Components, EntropyScaling, EquationOfState, IdealGas, Molarweight, NoResidual, Residual,
 };
 pub use errors::{EosError, EosResult};
 pub use phase_equilibria::{

diff --git a/feos-core/src/python/phase_equilibria.rs b/feos-core/src/python/phase_equilibria.rs
@@ -943,21 +943,23 @@ macro_rules! impl_phase_equilibrium {
                 if different_pressures {
                     dict.insert(String::from("pressure vapor"), p_v);
                     dict.insert(String::from("pressure liquid"), p_l);
-                }else {
+                } else {
                     dict.insert(String::from("pressure"), p_v);
                 }
                 dict.insert(String::from("density liquid"), self.0.liquid().density().convert_to(MOL / METER.powi::<P3>()).into_raw_vec_and_offset().0);
                 dict.insert(String::from("density vapor"), self.0.vapor().density().convert_to(MOL / METER.powi::<P3>()).into_raw_vec_and_offset().0);
-                dict.insert(String::from("mass density liquid"), self.0.liquid().mass_density().convert_to(KILOGRAM / METER.powi::<P3>()).into_raw_vec_and_offset().0);
-                dict.insert(String::from("mass density vapor"), self.0.vapor().mass_density().convert_to(KILOGRAM / METER.powi::<P3>()).into_raw_vec_and_offset().0);
                 dict.insert(String::from("molar enthalpy liquid"), self.0.liquid().molar_enthalpy(contributions).convert_to(KILO * JOULE / MOL).into_raw_vec_and_offset().0);
                 dict.insert(String::from("molar enthalpy vapor"), self.0.vapor().molar_enthalpy(contributions).convert_to(KILO * JOULE / MOL).into_raw_vec_and_offset().0);
                 dict.insert(String::from("molar entropy liquid"), self.0.liquid().molar_entropy(contributions).convert_to(KILO * JOULE / KELVIN / MOL).into_raw_vec_and_offset().0);
                 dict.insert(String::from("molar entropy vapor"), self.0.vapor().molar_entropy(contributions).convert_to(KILO * JOULE / KELVIN / MOL).into_raw_vec_and_offset().0);
-                dict.insert(String::from("specific enthalpy liquid"), self.0.liquid().specific_enthalpy(contributions).convert_to(KILO * JOULE / KILOGRAM).into_raw_vec_and_offset().0);
-                dict.insert(String::from("specific enthalpy vapor"), self.0.vapor().specific_enthalpy(contributions).convert_to(KILO * JOULE / KILOGRAM).into_raw_vec_and_offset().0);
-                dict.insert(String::from("specific entropy liquid"), self.0.liquid().specific_entropy(contributions).convert_to(KILO * JOULE / KELVIN / KILOGRAM).into_raw_vec_and_offset().0);
-                dict.insert(String::from("specific entropy vapor"), self.0.vapor().specific_entropy(contributions).convert_to(KILO * JOULE / KELVIN / KILOGRAM).into_raw_vec_and_offset().0);
+                if self.0.states[0].liquid().eos.residual.has_molar_weight() {
+                    dict.insert(String::from("mass density liquid"), self.0.liquid().mass_density().convert_to(KILOGRAM / METER.powi::<P3>()).into_raw_vec_and_offset().0);
+                    dict.insert(String::from("mass density vapor"), self.0.vapor().mass_density().convert_to(KILOGRAM / METER.powi::<P3>()).into_raw_vec_and_offset().0);
+                    dict.insert(String::from("specific enthalpy liquid"), self.0.liquid().specific_enthalpy(contributions).convert_to(KILO * JOULE / KILOGRAM).into_raw_vec_and_offset().0);
+                    dict.insert(String::from("specific enthalpy vapor"), self.0.vapor().specific_enthalpy(contributions).convert_to(KILO * JOULE / KILOGRAM).into_raw_vec_and_offset().0);
+                    dict.insert(String::from("specific entropy liquid"), self.0.liquid().specific_entropy(contributions).convert_to(KILO * JOULE / KELVIN / KILOGRAM).into_raw_vec_and_offset().0);
+                    dict.insert(String::from("specific entropy vapor"), self.0.vapor().specific_entropy(contributions).convert_to(KILO * JOULE / KELVIN / KILOGRAM).into_raw_vec_and_offset().0);
+                }
                 dict
             }
 

diff --git a/feos-core/src/python/state.rs b/feos-core/src/python/state.rs
@@ -1254,7 +1254,7 @@ macro_rules! impl_state {
             /// SIArray1
             #[pyo3(signature = (contributions=Contributions::Total), text_signature = "($self, contributions)")]
             fn molar_entropy(&self, contributions: Contributions) -> MolarEntropy<Array1<f64>> {
-                StateVec::from(self).molar_entropy(contributions).into()
+                StateVec::from(self).molar_entropy(contributions)
             }
 
             /// Return mass specific entropy.
@@ -1270,7 +1270,7 @@ macro_rules! impl_state {
             /// SIArray1
             #[pyo3(signature = (contributions=Contributions::Total), text_signature = "($self, contributions)")]
             fn specific_entropy(&self, contributions: Contributions) -> SpecificEntropy<Array1<f64>> {
-                StateVec::from(self).specific_entropy(contributions).into()
+                StateVec::from(self).specific_entropy(contributions)
             }
 
             /// Return molar enthalpy.
@@ -1286,7 +1286,7 @@ macro_rules! impl_state {
             /// SIArray1
             #[pyo3(signature = (contributions=Contributions::Total), text_signature = "($self, contributions)")]
             fn molar_enthalpy(&self, contributions: Contributions) -> MolarEnergy<Array1<f64>> {
-                StateVec::from(self).molar_enthalpy(contributions).into()
+                StateVec::from(self).molar_enthalpy(contributions)
             }
 
             /// Return mass specific enthalpy.
@@ -1302,18 +1302,18 @@ macro_rules! impl_state {
             /// SIArray1
             #[pyo3(signature = (contributions=Contributions::Total), text_signature = "($self, contributions)")]
             fn specific_enthalpy(&self, contributions: Contributions) -> SpecificEnergy<Array1<f64>> {
-                StateVec::from(self).specific_enthalpy(contributions).into()
+                StateVec::from(self).specific_enthalpy(contributions)
             }
 
 
             #[getter]
             fn get_temperature(&self) -> Temperature<Array1<f64>> {
-                StateVec::from(self).temperature().into()
+                StateVec::from(self).temperature()
             }
 
             #[getter]
             fn get_pressure(&self) -> Pressure<Array1<f64>> {
-                StateVec::from(self).pressure().into()
+                StateVec::from(self).pressure()
             }
 
             #[getter]
@@ -1323,12 +1323,12 @@ macro_rules! impl_state {
 
             #[getter]
             fn get_density(&self) -> Density<Array1<f64>> {
-                StateVec::from(self).density().into()
+                StateVec::from(self).density()
             }
 
             #[getter]
             fn get_moles<'py>(&self, py: Python<'py>) -> Moles<Array2<f64>> {
-                StateVec::from(self).moles().into()
+                StateVec::from(self).moles()
             }
 
             #[getter]
@@ -1337,13 +1337,13 @@ macro_rules! impl_state {
             }
 
             #[getter]
-            fn get_mass_density(&self) -> MassDensity<Array1<f64>> {
-                StateVec::from(self).mass_density().into()
+            fn get_mass_density(&self) -> Option<MassDensity<Array1<f64>>> {
+                self.0[0].eos.residual.has_molar_weight().then(|| StateVec::from(self).mass_density())
             }
 
             #[getter]
-            fn get_massfracs<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray2<f64>> {
-                StateVec::from(self).massfracs().into_pyarray_bound(py)
+            fn get_massfracs<'py>(&self, py: Python<'py>) -> Option<Bound<'py, PyArray2<f64>>> {
+                self.0[0].eos.residual.has_molar_weight().then(|| StateVec::from(self).massfracs().into_pyarray_bound(py))
             }
 
             /// Returns selected properties of a StateVec as dictionary.
@@ -1385,11 +1385,13 @@ macro_rules! impl_state {
                 dict.insert(String::from("temperature"), states.temperature().convert_to(KELVIN).into_raw_vec_and_offset().0);
                 dict.insert(String::from("pressure"), states.pressure().convert_to(PASCAL).into_raw_vec_and_offset().0);
                 dict.insert(String::from("density"), states.density().convert_to(MOL / METER.powi::<P3>()).into_raw_vec_and_offset().0);
-                dict.insert(String::from("mass density"), states.mass_density().convert_to(KILOGRAM / METER.powi::<P3>()).into_raw_vec_and_offset().0);
                 dict.insert(String::from("molar enthalpy"), states.molar_enthalpy(contributions).convert_to(KILO * JOULE / MOL).into_raw_vec_and_offset().0);
                 dict.insert(String::from("molar entropy"), states.molar_entropy(contributions).convert_to(KILO * JOULE / KELVIN / MOL).into_raw_vec_and_offset().0);
-                dict.insert(String::from("specific enthalpy"), states.specific_enthalpy(contributions).convert_to(KILO * JOULE / KILOGRAM).into_raw_vec_and_offset().0);
-                dict.insert(String::from("specific entropy"), states.specific_entropy(contributions).convert_to(KILO * JOULE / KELVIN / KILOGRAM).into_raw_vec_and_offset().0);
+                if states.0[0].eos.residual.has_molar_weight() {
+                    dict.insert(String::from("mass density"), states.mass_density().convert_to(KILOGRAM / METER.powi::<P3>()).into_raw_vec_and_offset().0);
+                    dict.insert(String::from("specific enthalpy"), states.specific_enthalpy(contributions).convert_to(KILO * JOULE / KILOGRAM).into_raw_vec_and_offset().0);
+                    dict.insert(String::from("specific entropy"), states.specific_entropy(contributions).convert_to(KILO * JOULE / KELVIN / KILOGRAM).into_raw_vec_and_offset().0);
+                }
                 dict
             }
         }

diff --git a/feos-core/src/python/user_defined.rs b/feos-core/src/python/user_defined.rs
@@ -1,5 +1,5 @@
 #![allow(non_snake_case)]
-use crate::{Components, IdealGas, Residual, StateHD};
+use crate::{Components, IdealGas, Molarweight, Residual, StateHD};
 use ndarray::{Array1, ScalarOperand};
 use num_dual::*;
 use numpy::convert::IntoPyArray;
@@ -203,7 +203,9 @@ macro_rules! impl_residual {
                 ) -> Vec<(String, D)> {
                 vec![("Python".to_string(), self.residual_helmholtz_energy(state))]
             }
+        }
 
+        impl Molarweight for PyResidual {
             fn molar_weight(&self) -> MolarWeight<Array1<f64>> {
                 Python::with_gil(|py| {
                     let py_result = self.0.bind(py).call_method0("molar_weight").unwrap();

diff --git a/feos-core/src/state/properties.rs b/feos-core/src/state/properties.rs
@@ -1,5 +1,5 @@
 use super::{Contributions, Derivative::*, PartialDerivative, State};
-use crate::equation_of_state::{IdealGas, Residual};
+use crate::equation_of_state::{IdealGas, Molarweight, Residual};
 use crate::ReferenceSystem;
 use ndarray::Array1;
 use quantity::*;
@@ -74,14 +74,6 @@ impl<E: Residual + IdealGas> State<E> {
         self.temperature * self.ds_dt(contributions) / self.total_moles
     }
 
-    /// Specific isochoric heat capacity: $c_v^{(m)}=\frac{C_v}{m}$
-    pub fn specific_isochoric_heat_capacity(
-        &self,
-        contributions: Contributions,
-    ) -> SpecificEntropy {
-        self.molar_isochoric_heat_capacity(contributions) / self.total_molar_weight()
-    }
-
     /// Partial derivative of the molar isochoric heat capacity w.r.t. temperature: $\left(\frac{\partial c_V}{\partial T}\right)_{V,N_i}$
     pub fn dc_v_dt(
         &self,
@@ -103,11 +95,6 @@ impl<E: Residual + IdealGas> State<E> {
         }
     }
 
-    /// Specific isobaric heat capacity: $c_p^{(m)}=\frac{C_p}{m}$
-    pub fn specific_isobaric_heat_capacity(&self, contributions: Contributions) -> SpecificEntropy {
-        self.molar_isobaric_heat_capacity(contributions) / self.total_molar_weight()
-    }
-
     /// Entropy: $S=-\left(\frac{\partial A}{\partial T}\right)_{V,N_i}$
     pub fn entropy(&self, contributions: Contributions) -> Entropy {
         Entropy::from_reduced(
@@ -120,11 +107,6 @@ impl<E: Residual + IdealGas> State<E> {
         self.entropy(contributions) / self.total_moles
     }
 
-    /// Specific entropy: $s^{(m)}=\frac{S}{m}$
-    pub fn specific_entropy(&self, contributions: Contributions) -> SpecificEntropy {
-        self.molar_entropy(contributions) / self.total_molar_weight()
-    }
-
     /// Partial molar entropy: $s_i=\left(\frac{\partial S}{\partial N_i}\right)_{T,p,N_j}$
     pub fn partial_molar_entropy(&self) -> MolarEntropy<Array1<f64>> {
         let c = Contributions::Total;
@@ -160,11 +142,6 @@ impl<E: Residual + IdealGas> State<E> {
         self.enthalpy(contributions) / self.total_moles
     }
 
-    /// Specific enthalpy: $h^{(m)}=\frac{H}{m}$
-    pub fn specific_enthalpy(&self, contributions: Contributions) -> SpecificEnergy {
-        self.molar_enthalpy(contributions) / self.total_molar_weight()
-    }
-
     /// Partial molar enthalpy: $h_i=\left(\frac{\partial H}{\partial N_i}\right)_{T,p,N_j}$
     pub fn partial_molar_enthalpy(&self) -> MolarEnergy<Array1<f64>> {
         let s = self.partial_molar_entropy();
@@ -184,11 +161,6 @@ impl<E: Residual + IdealGas> State<E> {
         self.helmholtz_energy(contributions) / self.total_moles
     }
 
-    /// Specific Helmholtz energy: $a^{(m)}=\frac{A}{m}$
-    pub fn specific_helmholtz_energy(&self, contributions: Contributions) -> SpecificEnergy {
-        self.molar_helmholtz_energy(contributions) / self.total_molar_weight()
-    }
-
     /// Internal energy: $U=A+TS$
     pub fn internal_energy(&self, contributions: Contributions) -> Energy {
         self.temperature * self.entropy(contributions) + self.helmholtz_energy(contributions)
@@ -199,11 +171,6 @@ impl<E: Residual + IdealGas> State<E> {
         self.internal_energy(contributions) / self.total_moles
     }
 
-    /// Specific internal energy: $u^{(m)}=\frac{U}{m}$
-    pub fn specific_internal_energy(&self, contributions: Contributions) -> SpecificEnergy {
-        self.molar_internal_energy(contributions) / self.total_molar_weight()
-    }
-
     /// Gibbs energy: $G=A+pV$
     pub fn gibbs_energy(&self, contributions: Contributions) -> Energy {
         self.pressure(contributions) * self.volume + self.helmholtz_energy(contributions)
@@ -214,11 +181,6 @@ impl<E: Residual + IdealGas> State<E> {
         self.gibbs_energy(contributions) / self.total_moles
     }
 
-    /// Specific Gibbs energy: $g^{(m)}=\frac{G}{m}$
-    pub fn specific_gibbs_energy(&self, contributions: Contributions) -> SpecificEnergy {
-        self.molar_gibbs_energy(contributions) / self.total_molar_weight()
-    }
-
     /// Joule Thomson coefficient: $\mu_{JT}=\left(\frac{\partial T}{\partial p}\right)_{H,N_i}$
     pub fn joule_thomson(&self) -> <Temperature as Div<Pressure>>::Output {
         let c = Contributions::Total;
@@ -278,6 +240,46 @@ impl<E: Residual + IdealGas> State<E> {
         }
         res
     }
+}
+
+impl<E: Residual + Molarweight + IdealGas> State<E> {
+    /// Specific isochoric heat capacity: $c_v^{(m)}=\frac{C_v}{m}$
+    pub fn specific_isochoric_heat_capacity(
+        &self,
+        contributions: Contributions,
+    ) -> SpecificEntropy {
+        self.molar_isochoric_heat_capacity(contributions) / self.total_molar_weight()
+    }
+
+    /// Specific isobaric heat capacity: $c_p^{(m)}=\frac{C_p}{m}$
+    pub fn specific_isobaric_heat_capacity(&self, contributions: Contributions) -> SpecificEntropy {
+        self.molar_isobaric_heat_capacity(contributions) / self.total_molar_weight()
+    }
+
+    /// Specific entropy: $s^{(m)}=\frac{S}{m}$
+    pub fn specific_entropy(&self, contributions: Contributions) -> SpecificEntropy {
+        self.molar_entropy(contributions) / self.total_molar_weight()
+    }
+
+    /// Specific enthalpy: $h^{(m)}=\frac{H}{m}$
+    pub fn specific_enthalpy(&self, contributions: Contributions) -> SpecificEnergy {
+        self.molar_enthalpy(contributions) / self.total_molar_weight()
+    }
+
+    /// Specific Helmholtz energy: $a^{(m)}=\frac{A}{m}$
+    pub fn specific_helmholtz_energy(&self, contributions: Contributions) -> SpecificEnergy {
+        self.molar_helmholtz_energy(contributions) / self.total_molar_weight()
+    }
+
+    /// Specific internal energy: $u^{(m)}=\frac{U}{m}$
+    pub fn specific_internal_energy(&self, contributions: Contributions) -> SpecificEnergy {
+        self.molar_internal_energy(contributions) / self.total_molar_weight()
+    }
+
+    /// Specific Gibbs energy: $g^{(m)}=\frac{G}{m}$
+    pub fn specific_gibbs_energy(&self, contributions: Contributions) -> SpecificEnergy {
+        self.molar_gibbs_energy(contributions) / self.total_molar_weight()
+    }
 
     /// Speed of sound: $c=\sqrt{\left(\frac{\partial p}{\partial\rho^{(m)}}\right)_{S,N_i}}$
     pub fn speed_of_sound(&self) -> Velocity {

diff --git a/feos-core/src/state/residual_properties.rs b/feos-core/src/state/residual_properties.rs
@@ -1,5 +1,5 @@
 use super::{Contributions, Derivative::*, PartialDerivative, State};
-use crate::equation_of_state::{EntropyScaling, Residual};
+use crate::equation_of_state::{EntropyScaling, Molarweight, Residual};
 use crate::errors::EosResult;
 use crate::phase_equilibria::PhaseEquilibrium;
 use crate::ReferenceSystem;
@@ -484,7 +484,9 @@ impl<E: Residual> State<E> {
     pub fn residual_molar_gibbs_energy(&self) -> MolarEnergy {
         self.residual_gibbs_energy() / self.total_moles
     }
+}
 
+impl<E: Residual + Molarweight> State<E> {
     /// Total molar weight: $MW=\sum_ix_iMW_i$
     pub fn total_molar_weight(&self) -> MolarWeight {
         (self.eos.molar_weight() * Dimensionless::new(&self.molefracs)).sum()