Proper implementation of cubic equations of state

Cargo.toml

@@ -40,6 +40,7 @@ indexmap = "2.0"
 rayon = { version = "1.7", optional = true }
 itertools = "0.13"
 typenum = "1.16"
+enum_dispatch = "0.3.13"
 
 [dependencies.pyo3]
 version = "0.21"
@@ -73,9 +74,10 @@ uvtheory = ["lazy_static"]
 pets = []
 saftvrqmie = []
 saftvrmie = []
+cubic = []
 rayon = ["dep:rayon", "ndarray/rayon", "feos-core/rayon", "feos-dft?/rayon"]
 python = ["pyo3", "numpy", "quantity/python", "feos-core/python", "feos-dft?/python", "rayon"]
-all_models = ["dft", "estimator", "pcsaft", "epcsaft", "gc_pcsaft", "uvtheory", "pets", "saftvrqmie", "saftvrmie"]
+all_models = ["dft", "estimator", "pcsaft", "epcsaft", "gc_pcsaft", "uvtheory", "pets", "saftvrqmie", "saftvrmie", "cubic"]
 
 [[bench]]
 name = "state_properties"

src/cubic/alpha/mathias_copeman.rs

@@ -0,0 +1,52 @@
+use std::sync::Arc;
+
+use feos_core::{parameter::ParameterError, EosResult};
+use ndarray::{Array1, Zip};
+use num_dual::DualNum;
+use serde::{Deserialize, Serialize};
+
+use crate::cubic::parameters::CubicParameters;
+
+use super::AlphaFunction;
+
+#[derive(Serialize, Deserialize, Clone, Debug)]
+pub struct MathiasCopeman(pub Vec<[f64; 3]>);
+
+impl AlphaFunction for MathiasCopeman {
+    #[inline]
+    fn alpha<D: DualNum<f64> + Copy>(
+        &self,
+        _: &Array1<f64>,
+        reduced_temperature: &Array1<D>,
+    ) -> Array1<D> {
+        Zip::from(reduced_temperature)
+            .and(&self.0)
+            .map_collect(|tr, c| {
+                let trsq = -tr.sqrt() + 1.0;
+                let a1 = trsq + c[0] + 1.0;
+                let a2 = match tr {
+                    tr if tr.re() < 1.0 => trsq * (c[1] + c[2]),
+                    _ => D::zero(),
+                };
+                (a1 + a2).powi(2)
+            })
+    }
+
+    fn validate(&self, parameters: &Arc<CubicParameters>) -> EosResult<()> {
+        if self.0.len() == parameters.tc.len() {
+            Ok(())
+        } else {
+            Err(ParameterError::IncompatibleParameters(
+                format!(
+                    "Mathias Copeman alpha function was initialized for {} components, but the equation of state contains {}.",
+                    self.0.len(), parameters.tc.len()
+                )
+            ).into())
+        }
+    }
+
+    fn subset(&self, component_list: &[usize]) -> Self {
+        let mi = component_list.iter().map(|&i| self.0[i]).collect();
+        Self(mi)
+    }
+}

src/cubic/alpha/mod.rs

@@ -0,0 +1,47 @@
+use std::sync::Arc;
+
+use enum_dispatch::enum_dispatch;
+use feos_core::EosResult;
+use ndarray::{Array1, ScalarOperand};
+use num_dual::DualNum;
+
+use super::parameters::CubicParameters;
+
+mod soave;
+pub use soave::{
+    PengRobinson1976, PengRobinson1978, PengRobinson2019, RedlichKwong1972, RedlichKwong2019, Soave,
+};
+mod mathias_copeman;
+pub use mathias_copeman::MathiasCopeman;
+mod twu;
+pub use twu::{GeneralizedTwu, Twu};
+
+#[enum_dispatch]
+pub trait AlphaFunction {
+    fn alpha<D: DualNum<f64> + Copy + ScalarOperand>(
+        &self,
+        acentric_factor: &Array1<f64>,
+        reduced_temperature: &Array1<D>,
+    ) -> Array1<D>;
+
+    /// Check for validity of alpha function against parameters, e.g.
+    /// to assert that the number of components match.
+    fn validate(&self, parameters: &Arc<CubicParameters>) -> EosResult<()>;
+
+    /// Generate the alpha function for a subset of components.
+    fn subset(&self, component_list: &[usize]) -> Self;
+}
+
+#[enum_dispatch(AlphaFunction)]
+#[derive(Debug, Clone)]
+pub enum Alpha {
+    Soave,
+    PengRobinson1976,
+    PengRobinson1978,
+    PengRobinson2019,
+    RedlichKwong1972,
+    RedlichKwong2019,
+    MathiasCopeman,
+    GeneralizedTwu,
+    Twu,
+}

src/cubic/alpha/soave.rs

@@ -0,0 +1,173 @@
+use std::sync::Arc;
+
+use feos_core::EosResult;
+use ndarray::{Array1, Zip};
+use num_dual::DualNum;
+use serde::{Deserialize, Serialize};
+
+use crate::cubic::parameters::CubicParameters;
+
+use super::AlphaFunction;
+
+#[derive(Serialize, Deserialize, Clone, Debug)]
+/// Generic version of Soave's function using 3rd order polynomial.
+pub struct Soave {
+    /// coefficients for m-polynomial
+    mi: Array1<f64>,
+}
+
+impl Soave {
+    pub fn new(mi: Array1<f64>) -> Self {
+        Soave { mi }
+    }
+}
+
+impl AlphaFunction for Soave {
+    #[inline]
+    fn alpha<D: DualNum<f64>>(
+        &self,
+        acentric_factor: &Array1<f64>,
+        reduced_temperature: &Array1<D>,
+    ) -> Array1<D> {
+        let m = self
+            .mi
+            .iter()
+            .enumerate()
+            .fold(Array1::zeros(acentric_factor.len()), |m, (i, &mi)| {
+                &m + acentric_factor.mapv(|w| w.powi(i as i32)) * mi
+            });
+        ((-reduced_temperature.mapv(|t| t.sqrt()) + 1.0) * m + 1.0).mapv(|a| a.powi(2))
+    }
+
+    fn validate(&self, _: &Arc<CubicParameters>) -> EosResult<()> {
+        Ok(())
+    }
+
+    fn subset(&self, _: &[usize]) -> Self {
+        self.clone()
+    }
+}
+
+#[derive(Serialize, Deserialize, Clone, Debug)]
+pub struct RedlichKwong1972;
+
+impl AlphaFunction for RedlichKwong1972 {
+    #[inline]
+    fn alpha<D: DualNum<f64>>(
+        &self,
+        acentric_factor: &Array1<f64>,
+        reduced_temperature: &Array1<D>,
+    ) -> Array1<D> {
+        let m = acentric_factor.mapv(|w| 0.48 + w * (1.574 - w * 0.176));
+        ((-reduced_temperature.mapv(|t| t.sqrt()) + 1.0) * m + 1.0).mapv(|a| a.powi(2))
+    }
+
+    fn validate(&self, _: &Arc<CubicParameters>) -> EosResult<()> {
+        Ok(())
+    }
+
+    fn subset(&self, _: &[usize]) -> Self {
+        Self
+    }
+}
+
+#[derive(Serialize, Deserialize, Clone, Debug)]
+pub struct PengRobinson1976;
+
+impl AlphaFunction for PengRobinson1976 {
+    #[inline]
+    fn alpha<D: DualNum<f64>>(
+        &self,
+        acentric_factor: &Array1<f64>,
+        reduced_temperature: &Array1<D>,
+    ) -> Array1<D> {
+        let m = acentric_factor.mapv(|w| 0.37464 + w * (1.54226 - w * 0.26992));
+        ((-reduced_temperature.mapv(|t| t.sqrt()) + 1.0) * m + 1.0).mapv(|a| a.powi(2))
+    }
+
+    fn validate(&self, _: &Arc<CubicParameters>) -> EosResult<()> {
+        Ok(())
+    }
+
+    fn subset(&self, _: &[usize]) -> Self {
+        Self
+    }
+}
+
+#[derive(Serialize, Deserialize, Clone, Debug)]
+pub struct PengRobinson1978;
+
+impl AlphaFunction for PengRobinson1978 {
+    #[inline]
+    fn alpha<D: DualNum<f64> + Copy>(
+        &self,
+        acentric_factor: &Array1<f64>,
+        reduced_temperature: &Array1<D>,
+    ) -> Array1<D> {
+        Zip::from(acentric_factor)
+            .and(reduced_temperature)
+            .map_collect(|&w, &tr| {
+                let m = if w <= 0.491 {
+                    0.37464 + w * (1.54226 - w * 0.26992)
+                } else {
+                    // use higher-order polynomial if w > w(n-decane)
+                    0.379642 + w * (1.48503 + w * (-0.164423 + w * 0.016666))
+                };
+                ((-tr.sqrt() + 1.0) * m + 1.0).powi(2)
+            })
+    }
+
+    fn validate(&self, _: &Arc<CubicParameters>) -> EosResult<()> {
+        Ok(())
+    }
+
+    fn subset(&self, _: &[usize]) -> Self {
+        Self
+    }
+}
+
+#[derive(Serialize, Deserialize, Clone, Debug)]
+pub struct RedlichKwong2019;
+
+impl AlphaFunction for RedlichKwong2019 {
+    #[inline]
+    fn alpha<D: DualNum<f64>>(
+        &self,
+        acentric_factor: &Array1<f64>,
+        reduced_temperature: &Array1<D>,
+    ) -> Array1<D> {
+        let m = acentric_factor.mapv(|w| 0.481 + w * (1.5963 + w * (-0.2963 + w * 0.1223)));
+        ((-reduced_temperature.mapv(|t| t.sqrt()) + 1.0) * m + 1.0).mapv(|a| a.powi(2))
+    }
+
+    fn validate(&self, _: &Arc<CubicParameters>) -> EosResult<()> {
+        Ok(())
+    }
+
+    fn subset(&self, _: &[usize]) -> Self {
+        Self
+    }
+}
+
+#[derive(Serialize, Deserialize, Clone, Debug)]
+pub struct PengRobinson2019;
+
+impl AlphaFunction for PengRobinson2019 {
+    #[inline]
+    fn alpha<D: DualNum<f64>>(
+        &self,
+        acentric_factor: &Array1<f64>,
+        reduced_temperature: &Array1<D>,
+    ) -> Array1<D> {
+        let m = acentric_factor.mapv(|w| 0.3919 + w * (1.4996 + w * (-0.2721 + w * 0.1063)));
+        ((-reduced_temperature.mapv(|t| t.sqrt()) + 1.0) * m + 1.0).mapv(|a| a.powi(2))
+    }
+
+    fn validate(&self, _: &Arc<CubicParameters>) -> EosResult<()> {
+        Ok(())
+    }
+
+    fn subset(&self, _: &[usize]) -> Self {
+        Self
+    }
+}