Integrated impl_cubics features into main

.github/workflows/test.yml

@@ -28,7 +28,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        model: [pcsaft, epcsaft, gc_pcsaft, pets, uvtheory, saftvrqmie, saftvrmie]
+        model: [pcsaft, epcsaft, gc_pcsaft, pets, uvtheory, saftvrqmie, saftvrmie, cubic]
 
     steps:
       - uses: actions/checkout@v4

Cargo.toml

@@ -34,6 +34,7 @@ serde_json = "1.0"
 indexmap = "2.0"
 itertools = "0.14"
 typenum = "1.16"
+enum_dispatch = "0.3.13"
 rayon = "1.5"
 petgraph = "0.8"
 rustdct = "0.7"

crates/feos/Cargo.toml

@@ -13,7 +13,7 @@ description = "FeOs - A framework for equations of state and classical density f
 [dependencies]
 quantity = { workspace = true }
 num-dual = { workspace = true }
-ndarray = { workspace = true }
+ndarray = { workspace = true } # updated without this feature?, features = ["serde"] }
 petgraph = { workspace = true, optional = true }
 thiserror = { workspace = true }
 num-traits = { workspace = true }
@@ -23,6 +23,7 @@ rayon = { workspace = true, optional = true }
 itertools = { workspace = true }
 typenum = { workspace = true }
 arrayvec = { workspace = true, features = ["serde"] }
+enum_dispatch = {workspace = true, optional = true}
 
 feos-core = { workspace = true }
 feos-derive = { workspace = true }
@@ -46,6 +47,7 @@ pets = []
 saftvrqmie = []
 saftvrmie = ["association"]
 rayon = ["dep:rayon", "ndarray/rayon", "feos-core/rayon", "feos-dft?/rayon"]
+cubic = ["enum_dispatch"]
 all_models = [
     "dft",
     "pcsaft",
@@ -55,6 +57,7 @@ all_models = [
     "pets",
     "saftvrqmie",
     "saftvrmie",
+    "cubic"
 ]
 
 [[bench]]

crates/feos/src/cubic/alpha/mathias_copeman.rs

@@ -0,0 +1,49 @@
+use super::AlphaFunction;
+use crate::cubic::parameters::CubicParameters;
+use feos_core::{FeosError, FeosResult};
+use ndarray::{Array1, Zip};
+use num_dual::DualNum;
+use serde::{Deserialize, Serialize};
+use std::sync::Arc;
+
+#[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>) -> FeosResult<()> {
+        if self.0.len() == parameters.tc.len() {
+            Ok(())
+        } else {
+            Err(FeosError::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)
+    }
+}

crates/feos/src/cubic/alpha/mod.rs

@@ -0,0 +1,45 @@
+use super::parameters::CubicParameters;
+use enum_dispatch::enum_dispatch;
+use feos_core::FeosResult;
+pub use mathias_copeman::MathiasCopeman;
+use ndarray::{Array1, ScalarOperand};
+use num_dual::DualNum;
+pub use soave::{
+    PengRobinson1976, PengRobinson1978, PengRobinson2019, RedlichKwong1972, RedlichKwong2019, Soave,
+};
+use std::sync::Arc;
+pub use twu::{GeneralizedTwu, Twu};
+
+mod mathias_copeman;
+mod soave;
+mod 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>) -> FeosResult<()>;
+
+    /// 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,
+}

crates/feos/src/cubic/alpha/soave.rs

@@ -0,0 +1,170 @@
+use super::AlphaFunction;
+use crate::cubic::parameters::CubicParameters;
+use feos_core::FeosResult;
+use ndarray::{Array1, Zip};
+use num_dual::DualNum;
+use serde::{Deserialize, Serialize};
+use std::sync::Arc;
+
+/// Generic version of Soave's function using 3rd order polynomial.
+#[derive(Serialize, Deserialize, Clone, Debug)]
+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>) -> FeosResult<()> {
+        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>) -> FeosResult<()> {
+        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>) -> FeosResult<()> {
+        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>) -> FeosResult<()> {
+        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>) -> FeosResult<()> {
+        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>) -> FeosResult<()> {
+        Ok(())
+    }
+
+    fn subset(&self, _: &[usize]) -> Self {
+        Self
+    }
+}