Added EntropyScaling in core, adjusted derive macros, started adjusti…

…ng PC-SAFT in feos, moved SolverOptions and Verbosity to state/mod.rs
feos-org · prehner · Jul 7, 2023 · Jun 1, 2023 · Jun 1, 2023 · Jun 1, 2023
commit 4476475948480c22ce482146da63ecd643b75556
diff --git a/feos-core/src/equation_of_state/__ideal_gas.rs b/feos-core/src/equation_of_state/__ideal_gas.rs
@@ -0,0 +1,100 @@
+use crate::StateHD;
+use ndarray::Array1;
+use num_dual::DualNum;
+use num_dual::*;
+use std::fmt;
+
+/// Ideal gas Helmholtz energy contribution that can
+/// be evaluated using generalized (hyper) dual numbers.
+///
+/// This trait needs to be implemented generically or for
+/// the specific types in the supertraits of [IdealGasContribution]
+/// so that the implementor can be used as an ideal gas
+/// contribution in the equation of state.
+pub trait IdealGasDual<D: DualNum<f64>> {
+    /// Return the number of components
+    fn components(&self) -> usize;
+
+    /// Return an equation of state consisting of the components
+    /// contained in component_list.
+    fn subset(&self, component_list: &[usize]) -> Self;
+
+    fn de_broglie_wavelength(&self, temperature: D) -> Array1<D>;
+
+    /// Evaluate the ideal gas contribution for a given state.
+    ///
+    /// In some cases it could be advantageous to overwrite this
+    /// implementation instead of implementing the de Broglie
+    /// wavelength.
+    fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
+        let lambda = self.de_broglie_wavelength(state.temperature);
+        ((lambda
+            + state.partial_density.mapv(|x| {
+                if x.re() == 0.0 {
+                    D::from(0.0)
+                } else {
+                    x.ln() - 1.0
+                }
+            }))
+            * &state.moles)
+            .sum()
+    }
+}
+
+pub struct DefaultIdealGas(pub usize);
+
+impl<D: DualNum<f64>> IdealGasDual<D> for DefaultIdealGas {
+    fn components(&self) -> usize {
+        self.0
+    }
+    fn subset(&self, component_list: &[usize]) -> Self {
+        Self(component_list.len())
+    }
+    fn de_broglie_wavelength(&self, temperature: D) -> Array1<D> {
+        Array1::zeros(self.0)
+    }
+}
+
+impl fmt::Display for DefaultIdealGas {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(f, "Ideal gas (default)")
+    }
+}
+
+pub trait IdealGas:
+    IdealGasDual<f64>
+    + IdealGasDual<Dual64>
+    + IdealGasDual<Dual<DualVec64<3>, f64>>
+    + IdealGasDual<HyperDual64>
+    + IdealGasDual<Dual2_64>
+    + IdealGasDual<Dual3_64>
+    + IdealGasDual<HyperDual<Dual64, f64>>
+    + IdealGasDual<HyperDual<DualVec64<2>, f64>>
+    + IdealGasDual<HyperDual<DualVec64<3>, f64>>
+    + IdealGasDual<Dual3<Dual64, f64>>
+    + IdealGasDual<Dual3<DualVec64<2>, f64>>
+    + IdealGasDual<Dual3<DualVec64<3>, f64>>
+    + fmt::Display
+    + Send
+    + Sync
+{
+}
+
+impl<T> IdealGas for T where
+    T: IdealGasDual<f64>
+        + IdealGasDual<Dual64>
+        + IdealGasDual<Dual<DualVec64<3>, f64>>
+        + IdealGasDual<HyperDual64>
+        + IdealGasDual<Dual2_64>
+        + IdealGasDual<Dual3_64>
+        + IdealGasDual<HyperDual<Dual64, f64>>
+        + IdealGasDual<HyperDual<DualVec64<2>, f64>>
+        + IdealGasDual<HyperDual<DualVec64<3>, f64>>
+        + IdealGasDual<Dual3<Dual64, f64>>
+        + IdealGasDual<Dual3<DualVec64<2>, f64>>
+        + IdealGasDual<Dual3<DualVec64<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
@@ -3,14 +3,14 @@ use num_dual::DualNum;
 use quantity::si::{SIArray1, MOL};
 use std::{fmt::Display, sync::Arc};
 
-pub use ideal_gas::IdealGas;
-pub use residual::Residual;
 pub mod debroglie;
 pub mod helmholtz_energy;
 pub mod ideal_gas;
 pub mod residual;
 use crate::StateHD;
 
+pub use ideal_gas::IdealGas;
+pub use residual::{EntropyScaling, Residual};
 pub use self::debroglie::{DeBroglieWavelength, DeBroglieWavelengthDual};
 pub use helmholtz_energy::{HelmholtzEnergy, HelmholtzEnergyDual};
 

diff --git a/feos-core/src/equation_of_state/residual.rs b/feos-core/src/equation_of_state/residual.rs
@@ -153,3 +153,28 @@ pub trait Residual: Send + Sync + fmt::Display {
             / (SIUnit::reference_density().powi(2) * SIUnit::reference_temperature()))
     }
 }
+
+/// Reference values and residual entropy correlations for entropy scaling.
+pub trait EntropyScaling {
+    fn viscosity_reference(
+        &self,
+        temperature: SINumber,
+        volume: SINumber,
+        moles: &SIArray1,
+    ) -> EosResult<SINumber>;
+    fn viscosity_correlation(&self, s_res: f64, x: &Array1<f64>) -> EosResult<f64>;
+    fn diffusion_reference(
+        &self,
+        temperature: SINumber,
+        volume: SINumber,
+        moles: &SIArray1,
+    ) -> EosResult<SINumber>;
+    fn diffusion_correlation(&self, s_res: f64, x: &Array1<f64>) -> EosResult<f64>;
+    fn thermal_conductivity_reference(
+        &self,
+        temperature: SINumber,
+        volume: SINumber,
+        moles: &SIArray1,
+    ) -> EosResult<SINumber>;
+    fn thermal_conductivity_correlation(&self, s_res: f64, x: &Array1<f64>) -> EosResult<f64>;
+}
diff --git a/feos-core/src/lib.rs b/feos-core/src/lib.rs
@@ -36,7 +36,8 @@ pub mod parameter;
 // mod phase_equilibria;
 mod state;
 pub use equation_of_state::{
-    EquationOfState, HelmholtzEnergy, HelmholtzEnergyDual, IdealGas, MolarWeight, Residual,
+    DeBroglieWavelength, EntropyScaling, EquationOfState, HelmholtzEnergy, HelmholtzEnergyDual,
+    IdealGas, MolarWeight, Residual,
 };
 pub use errors::{EosError, EosResult};
 // pub use phase_equilibria::{
@@ -199,53 +200,177 @@ mod tests {
             .build()?;
 
         // pressure
-        assert_relative_eq!(s.pressure(Contributions::Total), sr.pressure(Contributions::Total), max_relative = 1e-15);
-        assert_relative_eq!(s.pressure(Contributions::Residual), sr.pressure(Contributions::Residual), max_relative = 1e-15);
-        assert_relative_eq!(s.compressibility(Contributions::Total), sr.compressibility(Contributions::Total), max_relative = 1e-15);
-        assert_relative_eq!(s.compressibility(Contributions::Residual), sr.compressibility(Contributions::Residual), max_relative = 1e-15);
-
+        assert_relative_eq!(
+            s.pressure(Contributions::Total),
+            sr.pressure(Contributions::Total),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.pressure(Contributions::Residual),
+            sr.pressure(Contributions::Residual),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.compressibility(Contributions::Total),
+            sr.compressibility(Contributions::Total),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.compressibility(Contributions::Residual),
+            sr.compressibility(Contributions::Residual),
+            max_relative = 1e-15
+        );
+
         // residual properties
-        assert_relative_eq!(s.helmholtz_energy(Contributions::Residual), sr.residual_helmholtz_energy(), max_relative = 1e-15);
-        assert_relative_eq!(s.entropy(Contributions::Residual), sr.residual_entropy(), max_relative = 1e-15);
-        assert_relative_eq!(s.enthalpy(Contributions::Residual), sr.residual_enthalpy(), max_relative = 1e-15);
-        assert_relative_eq!(s.internal_energy(Contributions::Residual), sr.residual_internal_energy(), max_relative = 1e-15);
-        assert_relative_eq!(s.gibbs_energy(Contributions::Residual), sr.residual_gibbs_energy(), max_relative = 1e-15);
-        assert_relative_eq!(s.chemical_potential(Contributions::Residual), sr.residual_chemical_potential(), max_relative = 1e-15);
+        assert_relative_eq!(
+            s.helmholtz_energy(Contributions::Residual),
+            sr.residual_helmholtz_energy(),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.entropy(Contributions::Residual),
+            sr.residual_entropy(),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.enthalpy(Contributions::Residual),
+            sr.residual_enthalpy(),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.internal_energy(Contributions::Residual),
+            sr.residual_internal_energy(),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.gibbs_energy(Contributions::Residual),
+            sr.residual_gibbs_energy(),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.chemical_potential(Contributions::Residual),
+            sr.residual_chemical_potential(),
+            max_relative = 1e-15
+        );
 
         // pressure derivatives
-        assert_relative_eq!(s.structure_factor(), sr.structure_factor(), max_relative = 1e-15);
-        assert_relative_eq!(s.dp_dt(Contributions::Total), sr.dp_dt(Contributions::Total), max_relative = 1e-15);
-        assert_relative_eq!(s.dp_dt(Contributions::Residual), sr.dp_dt(Contributions::Residual), max_relative = 1e-15);
-        assert_relative_eq!(s.dp_dv(Contributions::Total), sr.dp_dv(Contributions::Total), max_relative = 1e-15);
-        assert_relative_eq!(s.dp_dv(Contributions::Residual), sr.dp_dv(Contributions::Residual), max_relative = 1e-15);
-        assert_relative_eq!(s.dp_drho(Contributions::Total), sr.dp_drho(Contributions::Total), max_relative = 1e-15);
-        assert_relative_eq!(s.dp_drho(Contributions::Residual), sr.dp_drho(Contributions::Residual), max_relative = 1e-15);
-        assert_relative_eq!(s.d2p_dv2(Contributions::Total), sr.d2p_dv2(Contributions::Total), max_relative = 1e-15);
-        assert_relative_eq!(s.d2p_dv2(Contributions::Residual), sr.d2p_dv2(Contributions::Residual), max_relative = 1e-15);
-        assert_relative_eq!(s.d2p_drho2(Contributions::Total), sr.d2p_drho2(Contributions::Total), max_relative = 1e-15);
-        assert_relative_eq!(s.d2p_drho2(Contributions::Residual), sr.d2p_drho2(Contributions::Residual), max_relative = 1e-15);
-        assert_relative_eq!(s.dp_dni(Contributions::Total), sr.dp_dni(Contributions::Total), max_relative = 1e-15);
-        assert_relative_eq!(s.dp_dni(Contributions::Residual), sr.dp_dni(Contributions::Residual), max_relative = 1e-15);
-
+        assert_relative_eq!(
+            s.structure_factor(),
+            sr.structure_factor(),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.dp_dt(Contributions::Total),
+            sr.dp_dt(Contributions::Total),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.dp_dt(Contributions::Residual),
+            sr.dp_dt(Contributions::Residual),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.dp_dv(Contributions::Total),
+            sr.dp_dv(Contributions::Total),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.dp_dv(Contributions::Residual),
+            sr.dp_dv(Contributions::Residual),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.dp_drho(Contributions::Total),
+            sr.dp_drho(Contributions::Total),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.dp_drho(Contributions::Residual),
+            sr.dp_drho(Contributions::Residual),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.d2p_dv2(Contributions::Total),
+            sr.d2p_dv2(Contributions::Total),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.d2p_dv2(Contributions::Residual),
+            sr.d2p_dv2(Contributions::Residual),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.d2p_drho2(Contributions::Total),
+            sr.d2p_drho2(Contributions::Total),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.d2p_drho2(Contributions::Residual),
+            sr.d2p_drho2(Contributions::Residual),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.dp_dni(Contributions::Total),
+            sr.dp_dni(Contributions::Total),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.dp_dni(Contributions::Residual),
+            sr.dp_dni(Contributions::Residual),
+            max_relative = 1e-15
+        );
+
         // entropy
-        assert_relative_eq!(s.ds_dt(Contributions::Residual), sr.ds_res_dt(), max_relative = 1e-15);
+        assert_relative_eq!(
+            s.ds_dt(Contributions::Residual),
+            sr.ds_res_dt(),
+            max_relative = 1e-15
+        );
 
         // chemical potential
-        assert_relative_eq!(s.dmu_dt(Contributions::Residual), sr.dmu_res_dt(), max_relative = 1e-15);
-        assert_relative_eq!(s.dmu_dni(Contributions::Residual), sr.dmu_res_dni(), max_relative = 1e-15);
-        assert_relative_eq!(s.dmu_dt(Contributions::Residual), sr.dmu_res_dt(), max_relative = 1e-15);
+        assert_relative_eq!(
+            s.dmu_dt(Contributions::Residual),
+            sr.dmu_res_dt(),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.dmu_dni(Contributions::Residual),
+            sr.dmu_res_dni(),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.dmu_dt(Contributions::Residual),
+            sr.dmu_res_dt(),
+            max_relative = 1e-15
+        );
 
         // fugacity
         assert_relative_eq!(s.ln_phi(), sr.ln_phi(), max_relative = 1e-15);
         assert_relative_eq!(s.dln_phi_dt(), sr.dln_phi_dt(), max_relative = 1e-15);
         assert_relative_eq!(s.dln_phi_dp(), sr.dln_phi_dp(), max_relative = 1e-15);
         assert_relative_eq!(s.dln_phi_dnj(), sr.dln_phi_dnj(), max_relative = 1e-15);
-        assert_relative_eq!(s.thermodynamic_factor(), sr.thermodynamic_factor(), max_relative = 1e-15);
+        assert_relative_eq!(
+            s.thermodynamic_factor(),
+            sr.thermodynamic_factor(),
+            max_relative = 1e-15
+        );
 
         // residual properties using multiple derivatives
-        assert_relative_eq!(s.c_v(Contributions::Residual), sr.c_v_res(), max_relative = 1e-15);
-        assert_relative_eq!(s.dc_v_dt(Contributions::Residual), sr.dc_v_res_dt(), max_relative = 1e-15);
-        assert_relative_eq!(s.c_p(Contributions::Residual), sr.c_p_res(), max_relative = 1e-15);
+        assert_relative_eq!(
+            s.c_v(Contributions::Residual),
+            sr.c_v_res(),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.dc_v_dt(Contributions::Residual),
+            sr.dc_v_res_dt(),
+            max_relative = 1e-15
+        );
+        assert_relative_eq!(
+            s.c_p(Contributions::Residual),
+            sr.c_p_res(),
+            max_relative = 1e-15
+        );
         Ok(())
     }
 }