PC SAFT done. Checking benches

feos-org · g-bauer · Feb 26, 2024 · Feb 12, 2024 · Feb 12, 2024 · Feb 13, 2024
commit 4ab68a99d3bc3cfb31bed59548ce1d21bfd1b040
diff --git a/benches/contributions.rs b/benches/contributions.rs
@@ -78,7 +78,7 @@ fn pcsaft(c: &mut Criterion) {
             let name1 = comp1.identifier.name.as_deref().unwrap();
             let name2 = comp2.identifier.name.as_deref().unwrap();
             let mix = format!("{name1}_{name2}");
-            group.bench_function(mix, |b| b.iter(|| eos.evaluate_residual(&state_hd)));
+            group.bench_function(mix, |b| b.iter(|| eos.residual_helmholtz_energy(&state_hd)));
         }
     }
 }

diff --git a/benches/dual_numbers.rs b/benches/dual_numbers.rs
@@ -8,7 +8,7 @@ use feos::pcsaft::{PcSaft, PcSaftParameters};
 use feos_core::si::*;
 use feos_core::{
     parameter::{IdentifierOption, Parameter},
-    Derivative, HelmholtzEnergy, HelmholtzEnergyDual, Residual, State, StateHD,
+    Derivative, Residual, State, StateHD,
 };
 use ndarray::{arr1, Array};
 use num_dual::DualNum;
@@ -29,11 +29,8 @@ fn state_pcsaft(parameters: PcSaftParameters) -> State<PcSaft> {
 }
 
 /// Residual Helmholtz energy given an equation of state and a StateHD.
-fn a_res<D: DualNum<f64> + Copy, E: Residual>(inp: (&Arc<E>, &StateHD<D>)) -> D
-where
-    (dyn HelmholtzEnergy + 'static): HelmholtzEnergyDual<D>,
-{
-    inp.0.evaluate_residual(inp.1)
+fn a_res<D: DualNum<f64> + Copy, E: Residual>(inp: (&Arc<E>, &StateHD<D>)) -> D {
+    inp.0.residual_helmholtz_energy(inp.1)
 }
 
 /// Benchmark for evaluation of the Helmholtz energy for different dual number types.

diff --git a/examples/gc_pcsaft_functional.ipynb b/examples/gc_pcsaft_functional.ipynb
diff --git a/examples/pcsaft_surface_tension.ipynb b/examples/pcsaft_surface_tension.ipynb
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, HelmholtzEnergy, HelmholtzEnergyDual, Residual};
+use crate::equation_of_state::{Components, Residual};
 use crate::parameter::{Identifier, Parameter, ParameterError, PureRecord};
 use crate::si::{MolarWeight, GRAM, MOL};
 use crate::state::StateHD;
@@ -150,33 +150,6 @@ struct PengRobinsonContribution {
     parameters: Arc<PengRobinsonParameters>,
 }
 
-impl<D: DualNum<f64> + Copy> HelmholtzEnergyDual<D> for PengRobinsonContribution {
-    fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
-        // temperature dependent a parameter
-        let p = &self.parameters;
-        let x = &state.molefracs;
-        let ak = (&p.tc.mapv(|tc| (D::one() - (state.temperature / tc).sqrt())) * &p.kappa + 1.0)
-            .mapv(|x| x.powi(2))
-            * &p.a;
-
-        // Mixing rules
-        let mut ak_mix = D::zero();
-        for i in 0..ak.len() {
-            for j in 0..ak.len() {
-                ak_mix += (ak[i] * ak[j]).sqrt() * (x[i] * x[j] * (1.0 - p.k_ij[(i, j)]));
-            }
-        }
-        let b = (x * &p.b).sum();
-
-        // Helmholtz energy
-        let n = state.moles.sum();
-        let v = state.volume;
-        n * ((v / (v - b * n)).ln()
-            - ak_mix / (b * SQRT_2 * 2.0 * state.temperature)
-                * ((v + b * n * (1.0 + SQRT_2)) / (v + b * n * (1.0 - SQRT_2))).ln())
-    }
-}
-
 impl fmt::Display for PengRobinsonContribution {
     fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
         write!(f, "Peng Robinson")
@@ -187,21 +160,12 @@ impl fmt::Display for PengRobinsonContribution {
 pub struct PengRobinson {
     /// Parameters
     parameters: Arc<PengRobinsonParameters>,
-    /// Non-ideal contributions to the Helmholtz energy
-    contributions: Vec<Box<dyn HelmholtzEnergy>>,
 }
 
 impl PengRobinson {
     /// Create a new equation of state from a set of parameters.
     pub fn new(parameters: Arc<PengRobinsonParameters>) -> Self {
-        let contributions: Vec<Box<dyn HelmholtzEnergy>> =
-            vec![Box::new(PengRobinsonContribution {
-                parameters: parameters.clone(),
-            })];
-        Self {
-            parameters,
-            contributions,
-        }
+        Self { parameters }
     }
 }
 
@@ -227,8 +191,38 @@ impl Residual for PengRobinson {
         0.9 / b
     }
 
-    fn contributions(&self) -> &[Box<dyn HelmholtzEnergy>] {
-        &self.contributions
+    fn residual_helmholtz_energy<D: DualNum<f64> + Copy>(&self, state: &StateHD<D>) -> D {
+        let p = &self.parameters;
+        let x = &state.molefracs;
+        let ak = (&p.tc.mapv(|tc| (D::one() - (state.temperature / tc).sqrt())) * &p.kappa + 1.0)
+            .mapv(|x| x.powi(2))
+            * &p.a;
+
+        // Mixing rules
+        let mut ak_mix = D::zero();
+        for i in 0..ak.len() {
+            for j in 0..ak.len() {
+                ak_mix += (ak[i] * ak[j]).sqrt() * (x[i] * x[j] * (1.0 - p.k_ij[(i, j)]));
+            }
+        }
+        let b = (x * &p.b).sum();
+
+        // Helmholtz energy
+        let n = state.moles.sum();
+        let v = state.volume;
+        n * ((v / (v - b * n)).ln()
+            - ak_mix / (b * SQRT_2 * 2.0 * state.temperature)
+                * ((v + b * n * (1.0 + SQRT_2)) / (v + b * n * (1.0 - SQRT_2))).ln())
+    }
+
+    fn residual_helmholtz_energy_contributions<D: DualNum<f64> + Copy>(
+        &self,
+        state: &StateHD<D>,
+    ) -> Vec<(String, D)> {
+        vec![(
+            "Peng Robinson".to_string(),
+            self.residual_helmholtz_energy(state),
+        )]
     }
 
     fn molar_weight(&self) -> MolarWeight<Array1<f64>> {

diff --git a/feos-core/src/density_iteration.rs b/feos-core/src/density_iteration.rs
@@ -135,7 +135,7 @@ pub fn density_iteration<E: Residual>(
     }
 }
 
-fn pressure_spinodal<E: Residual>(
+pub(crate) fn pressure_spinodal<E: Residual>(
     eos: &Arc<E>,
     temperature: Temperature,
     rho_init: Density,

diff --git a/feos-core/src/equation_of_state/helmholtz_energy.rs b/feos-core/src/equation_of_state/helmholtz_energy.rs
diff --git a/feos-core/src/equation_of_state/ideal_gas.rs b/feos-core/src/equation_of_state/ideal_gas.rs
@@ -2,24 +2,21 @@ use super::Components;
 use crate::StateHD;
 use ndarray::Array1;
 use num_dual::DualNum;
-use num_dual::*;
-use std::fmt;
 
 /// Ideal gas Helmholtz energy contribution.
 pub trait IdealGas: Components + Sync + Send {
     // Return a reference to the implementation of the de Broglie wavelength.
-    fn ideal_gas_model(&self) -> &dyn DeBroglieWavelength;
+    fn ideal_gas_name(&self) -> String;
+
+    fn ln_lambda3<D: DualNum<f64> + Copy>(&self, temperature: D) -> Array1<D>;
 
     /// Evaluate the ideal gas Helmholtz energy contribution for a given state.
     ///
     /// In some cases it could be advantageous to overwrite this
     /// implementation instead of implementing the de Broglie
     /// wavelength.
-    fn evaluate_ideal_gas<D: DualNum<f64> + Copy>(&self, state: &StateHD<D>) -> D
-    where
-        for<'a> dyn DeBroglieWavelength + 'a: DeBroglieWavelengthDual<D>,
-    {
-        let ln_lambda3 = self.ideal_gas_model().ln_lambda3(state.temperature);
+    fn ideal_gas_helmholtz_energy<D: DualNum<f64> + Copy>(&self, state: &StateHD<D>) -> D {
+        let ln_lambda3 = self.ln_lambda3(state.temperature);
         ((ln_lambda3
             + state.partial_density.mapv(|x| {
                 if x.re() == 0.0 {
@@ -32,59 +29,3 @@ pub trait IdealGas: Components + Sync + Send {
             .sum()
     }
 }
-
-/// Implementation of an ideal gas model in terms of the
-/// logarithm of the cubic thermal de Broglie wavelength
-/// in units ln(A³).
-///
-/// This trait needs to be implemented generically or for
-/// the specific types in the supertraits of [DeBroglieWavelength]
-/// so that the implementor can be used as an ideal gas
-/// contribution in the equation of state.
-pub trait DeBroglieWavelengthDual<D: DualNum<f64>> {
-    fn ln_lambda3(&self, temperature: D) -> Array1<D>;
-}
-
-/// Object safe version of the [DeBroglieWavelengthDual] trait.
-///
-/// The trait is implemented automatically for every struct that implements
-/// the supertraits.
-pub trait DeBroglieWavelength:
-    DeBroglieWavelengthDual<f64>
-    + DeBroglieWavelengthDual<Dual64>
-    + DeBroglieWavelengthDual<Dual<DualSVec64<3>, f64>>
-    + DeBroglieWavelengthDual<HyperDual64>
-    + DeBroglieWavelengthDual<Dual2_64>
-    + DeBroglieWavelengthDual<Dual3_64>
-    + DeBroglieWavelengthDual<HyperDual<Dual64, f64>>
-    + DeBroglieWavelengthDual<HyperDual<DualSVec64<2>, f64>>
-    + DeBroglieWavelengthDual<HyperDual<DualSVec64<3>, f64>>
-    + DeBroglieWavelengthDual<Dual2<Dual64, f64>>
-    + DeBroglieWavelengthDual<Dual3<Dual64, f64>>
-    + DeBroglieWavelengthDual<Dual3<DualSVec64<2>, f64>>
-    + DeBroglieWavelengthDual<Dual3<DualSVec64<3>, f64>>
-    + fmt::Display
-    + Send
-    + Sync
-{
-}
-
-impl<T> DeBroglieWavelength for T where
-    T: DeBroglieWavelengthDual<f64>
-        + DeBroglieWavelengthDual<Dual64>
-        + DeBroglieWavelengthDual<Dual<DualSVec64<3>, f64>>
-        + DeBroglieWavelengthDual<HyperDual64>
-        + DeBroglieWavelengthDual<Dual2_64>
-        + DeBroglieWavelengthDual<Dual3_64>
-        + DeBroglieWavelengthDual<HyperDual<Dual64, f64>>
-        + DeBroglieWavelengthDual<HyperDual<DualSVec64<2>, f64>>
-        + DeBroglieWavelengthDual<HyperDual<DualSVec64<3>, f64>>
-        + DeBroglieWavelengthDual<Dual2<Dual64, f64>>
-        + DeBroglieWavelengthDual<Dual3<Dual64, f64>>
-        + DeBroglieWavelengthDual<Dual3<DualSVec64<2>, f64>>
-        + DeBroglieWavelengthDual<Dual3<DualSVec64<3>, f64>>
-        + fmt::Display
-        + Send
-        + Sync
-{
-}
diff --git a/feos-core/src/equation_of_state/mod.rs b/feos-core/src/equation_of_state/mod.rs
@@ -5,12 +5,10 @@ use crate::{
 use ndarray::Array1;
 use std::sync::Arc;
 
-mod helmholtz_energy;
 mod ideal_gas;
 mod residual;
 
-pub use helmholtz_energy::{HelmholtzEnergy, HelmholtzEnergyDual};
-pub use ideal_gas::{DeBroglieWavelength, DeBroglieWavelengthDual, IdealGas};
+pub use ideal_gas::IdealGas;
 pub use residual::{EntropyScaling, NoResidual, Residual};
 
 /// The number of components that the model is initialized for.
@@ -73,8 +71,12 @@ impl<I: Components, R: Components> Components for EquationOfState<I, R> {
 }
 
 impl<I: IdealGas, R: Components + Sync + Send> IdealGas for EquationOfState<I, R> {
-    fn ideal_gas_model(&self) -> &dyn DeBroglieWavelength {
-        self.ideal_gas.ideal_gas_model()
+    fn ideal_gas_name(&self) -> String {
+        self.ideal_gas.ideal_gas_name()
+    }
+
+    fn ln_lambda3<D: num_dual::DualNum<f64> + Copy>(&self, temperature: D) -> Array1<D> {
+        self.ideal_gas.ln_lambda3(temperature)
     }
 }
 
@@ -83,8 +85,18 @@ impl<I: IdealGas, R: Residual> Residual for EquationOfState<I, R> {
         self.residual.compute_max_density(moles)
     }
 
-    fn contributions(&self) -> &[Box<dyn HelmholtzEnergy>] {
-        self.residual.contributions()
+    fn residual_helmholtz_energy<D: num_dual::DualNum<f64> + Copy>(
+        &self,
+        state: &crate::StateHD<D>,
+    ) -> D {
+        self.residual.residual_helmholtz_energy(state)
+    }
+
+    fn residual_helmholtz_energy_contributions<D: num_dual::DualNum<f64> + Copy>(
+        &self,
+        state: &crate::StateHD<D>,
+    ) -> Vec<(String, D)> {
+        self.residual.residual_helmholtz_energy_contributions(state)
     }
 
     fn molar_weight(&self) -> MolarWeight<Array1<f64>> {