From d3b4c2305b70239e9e3ba080a76c42160eed58d1 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Wed, 7 Feb 2024 18:01:48 +0100 Subject: [PATCH 1/2] Reuse temperature-dependent properties across contributions --- feos-core/src/cubic.rs | 18 +- .../src/equation_of_state/helmholtz_energy.rs | 65 +++---- feos-core/src/equation_of_state/mod.rs | 9 +- feos-core/src/equation_of_state/residual.rs | 169 ++++++++++-------- feos-core/src/lib.rs | 2 +- 5 files changed, 148 insertions(+), 115 deletions(-) diff --git a/feos-core/src/cubic.rs b/feos-core/src/cubic.rs index ccd6fff7d..88f84c0d0 100644 --- a/feos-core/src/cubic.rs +++ b/feos-core/src/cubic.rs @@ -8,11 +8,13 @@ use crate::equation_of_state::{Components, HelmholtzEnergy, HelmholtzEnergyDual, use crate::parameter::{Identifier, Parameter, ParameterError, PureRecord}; use crate::si::{MolarWeight, GRAM, MOL}; use crate::state::StateHD; +use crate::Properties; use ndarray::{Array1, Array2}; use num_dual::DualNum; use serde::{Deserialize, Serialize}; use std::f64::consts::SQRT_2; use std::fmt; +use std::marker::PhantomData; use std::sync::Arc; const KB_A3: f64 = 13806490.0; @@ -150,8 +152,8 @@ struct PengRobinsonContribution { parameters: Arc, } -impl + Copy> HelmholtzEnergyDual for PengRobinsonContribution { - fn helmholtz_energy(&self, state: &StateHD) -> D { +impl + Copy> HelmholtzEnergyDual, D> for PengRobinsonContribution { + fn helmholtz_energy(&self, state: &StateHD, properties: &PhantomData) -> D { // temperature dependent a parameter let p = &self.parameters; let x = &state.molefracs; @@ -188,13 +190,13 @@ pub struct PengRobinson { /// Parameters parameters: Arc, /// Non-ideal contributions to the Helmholtz energy - contributions: Vec>, + contributions: Vec>>, } impl PengRobinson { /// Create a new equation of state from a set of parameters. pub fn new(parameters: Arc) -> Self { - let contributions: Vec> = + let contributions: Vec>> = vec![Box::new(PengRobinsonContribution { parameters: parameters.clone(), })]; @@ -222,12 +224,18 @@ impl Components for PengRobinson { } impl Residual for PengRobinson { + type Properties = PhantomData; + + fn properties>(&self, _: D) -> PhantomData { + PhantomData + } + fn compute_max_density(&self, moles: &Array1) -> f64 { let b = (moles * &self.parameters.b).sum() / moles.sum(); 0.9 / b } - fn contributions(&self) -> &[Box] { + fn contributions(&self) -> &[Box>] { &self.contributions } diff --git a/feos-core/src/equation_of_state/helmholtz_energy.rs b/feos-core/src/equation_of_state/helmholtz_energy.rs index 11d2c2703..af438b4c1 100644 --- a/feos-core/src/equation_of_state/helmholtz_energy.rs +++ b/feos-core/src/equation_of_state/helmholtz_energy.rs @@ -1,4 +1,4 @@ -use crate::StateHD; +use crate::{Residual, StateHD}; use num_dual::*; use std::fmt; @@ -9,49 +9,52 @@ use std::fmt; /// the specific types in the supertraits of [HelmholtzEnergy] /// so that the implementor can be used as a Helmholtz energy /// contribution in the equation of state. -pub trait HelmholtzEnergyDual> { +pub trait HelmholtzEnergyDual> { /// The Helmholtz energy contribution $\beta A$ of a given state in reduced units. - fn helmholtz_energy(&self, state: &StateHD) -> D; + fn helmholtz_energy(&self, state: &StateHD, properties: &P) -> D; } /// Object safe version of the [HelmholtzEnergyDual] trait. /// /// The trait is implemented automatically for every struct that implements /// the supertraits. -pub trait HelmholtzEnergy: - HelmholtzEnergyDual - + HelmholtzEnergyDual - + HelmholtzEnergyDual, f64>> - + HelmholtzEnergyDual - + HelmholtzEnergyDual - + HelmholtzEnergyDual - + HelmholtzEnergyDual> - + HelmholtzEnergyDual, f64>> - + HelmholtzEnergyDual, f64>> - + HelmholtzEnergyDual> - + HelmholtzEnergyDual> - + HelmholtzEnergyDual, f64>> - + HelmholtzEnergyDual, f64>> +pub trait HelmholtzEnergy: HelmholtzEnergyDual, f64> + + HelmholtzEnergyDual, Dual64> + + HelmholtzEnergyDual, f64>>, Dual, f64>> + + HelmholtzEnergyDual, HyperDual64> + + HelmholtzEnergyDual, Dual2_64> + + HelmholtzEnergyDual, Dual3_64> + + HelmholtzEnergyDual>, HyperDual> + + HelmholtzEnergyDual, f64>>, HyperDual, f64>> + + HelmholtzEnergyDual, f64>>, HyperDual, f64>> + + HelmholtzEnergyDual>, Dual2> + + HelmholtzEnergyDual>, Dual3> + + HelmholtzEnergyDual, f64>>, Dual3, f64>> + + HelmholtzEnergyDual, f64>>, Dual3, f64>> + fmt::Display + Send + Sync { } -impl HelmholtzEnergy for T where - T: HelmholtzEnergyDual - + HelmholtzEnergyDual - + HelmholtzEnergyDual, f64>> - + HelmholtzEnergyDual - + HelmholtzEnergyDual - + HelmholtzEnergyDual - + HelmholtzEnergyDual> - + HelmholtzEnergyDual, f64>> - + HelmholtzEnergyDual, f64>> - + HelmholtzEnergyDual> - + HelmholtzEnergyDual> - + HelmholtzEnergyDual, f64>> - + HelmholtzEnergyDual, f64>> +impl HelmholtzEnergy

for T where + T: HelmholtzEnergyDual, f64> + + HelmholtzEnergyDual, Dual64> + + HelmholtzEnergyDual, f64>>, Dual, f64>> + + HelmholtzEnergyDual, HyperDual64> + + HelmholtzEnergyDual, Dual2_64> + + HelmholtzEnergyDual, Dual3_64> + + HelmholtzEnergyDual>, HyperDual> + + HelmholtzEnergyDual< + P::Properties, f64>>, + HyperDual, f64>, + > + HelmholtzEnergyDual< + P::Properties, f64>>, + HyperDual, f64>, + > + HelmholtzEnergyDual>, Dual2> + + HelmholtzEnergyDual>, Dual3> + + HelmholtzEnergyDual, f64>>, Dual3, f64>> + + HelmholtzEnergyDual, f64>>, Dual3, 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 index a1ac42a64..30f272bea 100644 --- a/feos-core/src/equation_of_state/mod.rs +++ b/feos-core/src/equation_of_state/mod.rs @@ -11,7 +11,7 @@ mod residual; pub use helmholtz_energy::{HelmholtzEnergy, HelmholtzEnergyDual}; pub use ideal_gas::{DeBroglieWavelength, DeBroglieWavelengthDual, IdealGas}; -pub use residual::{EntropyScaling, NoResidual, Residual}; +pub use residual::{EntropyScaling, NoResidual, Properties, Residual}; /// The number of components that the model is initialized for. pub trait Components { @@ -79,11 +79,16 @@ impl IdealGas for EquationOfState Residual for EquationOfState { + type Properties = R::Properties; + + fn properties>(&self, temperature: D) -> Self::Properties { + self.residual.properties(temperature) + } fn compute_max_density(&self, moles: &Array1) -> f64 { self.residual.compute_max_density(moles) } - fn contributions(&self) -> &[Box] { + fn contributions(&self) -> &[Box>] { self.residual.contributions() } diff --git a/feos-core/src/equation_of_state/residual.rs b/feos-core/src/equation_of_state/residual.rs index e8e2819e2..1d2ec5375 100644 --- a/feos-core/src/equation_of_state/residual.rs +++ b/feos-core/src/equation_of_state/residual.rs @@ -5,10 +5,19 @@ use crate::{EosError, EosResult}; use ndarray::prelude::*; use num_dual::*; use num_traits::{One, Zero}; +use std::marker::PhantomData; use std::ops::Div; +pub trait Properties { + type Values; +} + /// A reisdual Helmholtz energy model. pub trait Residual: Components + Send + Sync { + type Properties; + + fn properties>(&self, temperature: D) -> Self::Properties; + /// Return the maximum density in Angstrom^-3. /// /// This value is used as an estimate for a liquid phase for phase @@ -19,7 +28,7 @@ pub trait Residual: Components + Send + Sync { /// Return a slice of the individual contributions (excluding the ideal gas) /// of the equation of state. - fn contributions(&self) -> &[Box]; + fn contributions(&self) -> &[Box>]; /// Molar weight of all components. /// @@ -29,11 +38,12 @@ pub trait Residual: Components + Send + Sync { /// Evaluate the residual reduced Helmholtz energy $\beta A^\mathrm{res}$. fn evaluate_residual + Copy>(&self, state: &StateHD) -> D where - dyn HelmholtzEnergy: HelmholtzEnergyDual, + dyn HelmholtzEnergy: HelmholtzEnergyDual, D>, { + let properties = self.properties(state.temperature); self.contributions() .iter() - .map(|c| c.helmholtz_energy(state)) + .map(|c| c.helmholtz_energy(state, &properties)) .sum() } @@ -44,11 +54,12 @@ pub trait Residual: Components + Send + Sync { state: &StateHD, ) -> Vec<(String, D)> where - dyn HelmholtzEnergy: HelmholtzEnergyDual, + dyn HelmholtzEnergy: HelmholtzEnergyDual, D>, { + let properties = self.properties(state.temperature); self.contributions() .iter() - .map(|c| (c.to_string(), c.helmholtz_energy(state))) + .map(|c| (c.to_string(), c.helmholtz_energy(state, &properties))) .collect() } @@ -82,76 +93,76 @@ pub trait Residual: Components + Send + Sync { Ok(Density::from_reduced(self.compute_max_density(&mr))) } - /// Calculate the second virial coefficient $B(T)$ - fn second_virial_coefficient( - &self, - temperature: Temperature, - moles: Option<&Moles>>, - ) -> EosResult<>::Output> { - let mr = self.validate_moles(moles)?; - let x = (&mr / mr.sum()).into_value(); - let mut rho = HyperDual64::zero(); - rho.eps1 = 1.0; - rho.eps2 = 1.0; - let t = HyperDual64::from(temperature.to_reduced()); - let s = StateHD::new_virial(t, rho, x); - Ok(Quantity::from_reduced( - self.evaluate_residual(&s).eps1eps2 * 0.5, - )) - } - - /// Calculate the third virial coefficient $C(T)$ - #[allow(clippy::type_complexity)] - fn third_virial_coefficient( - &self, - temperature: Temperature, - moles: Option<&Moles>>, - ) -> EosResult<<>::Output as Div>::Output> { - let mr = self.validate_moles(moles)?; - let x = (&mr / mr.sum()).into_value(); - let rho = Dual3_64::zero().derivative(); - let t = Dual3_64::from(temperature.to_reduced()); - let s = StateHD::new_virial(t, rho, x); - Ok(Quantity::from_reduced(self.evaluate_residual(&s).v3 / 3.0)) - } - - /// Calculate the temperature derivative of the second virial coefficient $B'(T)$ - #[allow(clippy::type_complexity)] - fn second_virial_coefficient_temperature_derivative( - &self, - temperature: Temperature, - moles: Option<&Moles>>, - ) -> EosResult<<>::Output as Div>::Output> { - let mr = self.validate_moles(moles)?; - let x = (&mr / mr.sum()).into_value(); - let mut rho = HyperDual::zero(); - rho.eps1 = Dual64::one(); - rho.eps2 = Dual64::one(); - let t = HyperDual::from_re(Dual64::from(temperature.to_reduced()).derivative()); - let s = StateHD::new_virial(t, rho, x); - Ok(Quantity::from_reduced( - self.evaluate_residual(&s).eps1eps2.eps * 0.5, - )) - } - - /// Calculate the temperature derivative of the third virial coefficient $C'(T)$ - #[allow(clippy::type_complexity)] - fn third_virial_coefficient_temperature_derivative( - &self, - temperature: Temperature, - moles: Option<&Moles>>, - ) -> EosResult< - <<>::Output as Div>::Output as Div>::Output, - > { - let mr = self.validate_moles(moles)?; - let x = (&mr / mr.sum()).into_value(); - let rho = Dual3::zero().derivative(); - let t = Dual3::from_re(Dual64::from(temperature.to_reduced()).derivative()); - let s = StateHD::new_virial(t, rho, x); - Ok(Quantity::from_reduced( - self.evaluate_residual(&s).v3.eps / 3.0, - )) - } + // /// Calculate the second virial coefficient $B(T)$ + // fn second_virial_coefficient( + // &self, + // temperature: Temperature, + // moles: Option<&Moles>>, + // ) -> EosResult<>::Output> { + // let mr = self.validate_moles(moles)?; + // let x = (&mr / mr.sum()).into_value(); + // let mut rho = HyperDual64::zero(); + // rho.eps1 = 1.0; + // rho.eps2 = 1.0; + // let t = HyperDual64::from(temperature.to_reduced()); + // let s = StateHD::new_virial(t, rho, x); + // Ok(Quantity::from_reduced( + // self.evaluate_residual(&s).eps1eps2 * 0.5, + // )) + // } + + // /// Calculate the third virial coefficient $C(T)$ + // #[allow(clippy::type_complexity)] + // fn third_virial_coefficient( + // &self, + // temperature: Temperature, + // moles: Option<&Moles>>, + // ) -> EosResult<<>::Output as Div>::Output> { + // let mr = self.validate_moles(moles)?; + // let x = (&mr / mr.sum()).into_value(); + // let rho = Dual3_64::zero().derivative(); + // let t = Dual3_64::from(temperature.to_reduced()); + // let s = StateHD::new_virial(t, rho, x); + // Ok(Quantity::from_reduced(self.evaluate_residual(&s).v3 / 3.0)) + // } + + // /// Calculate the temperature derivative of the second virial coefficient $B'(T)$ + // #[allow(clippy::type_complexity)] + // fn second_virial_coefficient_temperature_derivative( + // &self, + // temperature: Temperature, + // moles: Option<&Moles>>, + // ) -> EosResult<<>::Output as Div>::Output> { + // let mr = self.validate_moles(moles)?; + // let x = (&mr / mr.sum()).into_value(); + // let mut rho = HyperDual::zero(); + // rho.eps1 = Dual64::one(); + // rho.eps2 = Dual64::one(); + // let t = HyperDual::from_re(Dual64::from(temperature.to_reduced()).derivative()); + // let s = StateHD::new_virial(t, rho, x); + // Ok(Quantity::from_reduced( + // self.evaluate_residual(&s).eps1eps2.eps * 0.5, + // )) + // } + + // /// Calculate the temperature derivative of the third virial coefficient $C'(T)$ + // #[allow(clippy::type_complexity)] + // fn third_virial_coefficient_temperature_derivative( + // &self, + // temperature: Temperature, + // moles: Option<&Moles>>, + // ) -> EosResult< + // <<>::Output as Div>::Output as Div>::Output, + // > { + // let mr = self.validate_moles(moles)?; + // let x = (&mr / mr.sum()).into_value(); + // let rho = Dual3::zero().derivative(); + // let t = Dual3::from_re(Dual64::from(temperature.to_reduced()).derivative()); + // let s = StateHD::new_virial(t, rho, x); + // Ok(Quantity::from_reduced( + // self.evaluate_residual(&s).v3.eps / 3.0, + // )) + // } } /// Reference values and residual entropy correlations for entropy scaling. @@ -193,11 +204,17 @@ impl Components for NoResidual { } impl Residual for NoResidual { + type Properties = PhantomData; + + fn properties>(&self, temperature: D) -> PhantomData { + PhantomData + } + fn compute_max_density(&self, _: &Array1) -> f64 { 1.0 } - fn contributions(&self) -> &[Box] { + fn contributions(&self) -> &[Box>] { &[] } diff --git a/feos-core/src/lib.rs b/feos-core/src/lib.rs index 5baeb764c..90ddeee49 100644 --- a/feos-core/src/lib.rs +++ b/feos-core/src/lib.rs @@ -33,7 +33,7 @@ pub mod si; mod state; pub use equation_of_state::{ Components, DeBroglieWavelength, DeBroglieWavelengthDual, EntropyScaling, EquationOfState, - HelmholtzEnergy, HelmholtzEnergyDual, IdealGas, NoResidual, Residual, + HelmholtzEnergy, HelmholtzEnergyDual, IdealGas, NoResidual, Properties, Residual, }; pub use errors::{EosError, EosResult}; pub use phase_equilibria::{ From 54b287d83ed1060d055e693a2798c4d2c07d5368 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Wed, 7 Feb 2024 22:10:50 +0100 Subject: [PATCH 2/2] Remove extra trait --- feos-core/src/cubic.rs | 4 +- feos-core/src/equation_of_state/mod.rs | 10 +- feos-core/src/equation_of_state/residual.rs | 162 ++++++++++---------- feos-core/src/lib.rs | 2 +- feos-dft/src/functional.rs | 25 ++- feos-dft/src/functional_contribution.rs | 11 +- feos-dft/src/ideal_chain_contribution.rs | 5 +- src/association/mod.rs | 2 +- src/hard_sphere/mod.rs | 4 +- src/pcsaft/eos/dispersion.rs | 4 +- src/pcsaft/eos/hard_chain.rs | 4 +- src/pcsaft/eos/mod.rs | 12 +- src/pcsaft/eos/polar.rs | 8 +- 13 files changed, 145 insertions(+), 108 deletions(-) diff --git a/feos-core/src/cubic.rs b/feos-core/src/cubic.rs index 88f84c0d0..0c615714a 100644 --- a/feos-core/src/cubic.rs +++ b/feos-core/src/cubic.rs @@ -8,7 +8,6 @@ use crate::equation_of_state::{Components, HelmholtzEnergy, HelmholtzEnergyDual, use crate::parameter::{Identifier, Parameter, ParameterError, PureRecord}; use crate::si::{MolarWeight, GRAM, MOL}; use crate::state::StateHD; -use crate::Properties; use ndarray::{Array1, Array2}; use num_dual::DualNum; use serde::{Deserialize, Serialize}; @@ -153,7 +152,7 @@ struct PengRobinsonContribution { } impl + Copy> HelmholtzEnergyDual, D> for PengRobinsonContribution { - fn helmholtz_energy(&self, state: &StateHD, properties: &PhantomData) -> D { + fn helmholtz_energy(&self, state: &StateHD, _: &PhantomData) -> D { // temperature dependent a parameter let p = &self.parameters; let x = &state.molefracs; @@ -225,6 +224,7 @@ impl Components for PengRobinson { impl Residual for PengRobinson { type Properties = PhantomData; + type Inner = Self; fn properties>(&self, _: D) -> PhantomData { PhantomData diff --git a/feos-core/src/equation_of_state/mod.rs b/feos-core/src/equation_of_state/mod.rs index 30f272bea..bca430600 100644 --- a/feos-core/src/equation_of_state/mod.rs +++ b/feos-core/src/equation_of_state/mod.rs @@ -11,7 +11,7 @@ mod residual; pub use helmholtz_energy::{HelmholtzEnergy, HelmholtzEnergyDual}; pub use ideal_gas::{DeBroglieWavelength, DeBroglieWavelengthDual, IdealGas}; -pub use residual::{EntropyScaling, NoResidual, Properties, Residual}; +pub use residual::{EntropyScaling, NoResidual, Residual}; /// The number of components that the model is initialized for. pub trait Components { @@ -80,15 +80,19 @@ impl IdealGas for EquationOfState Residual for EquationOfState { type Properties = R::Properties; + type Inner = R::Inner; - fn properties>(&self, temperature: D) -> Self::Properties { + fn properties>( + &self, + temperature: D, + ) -> ::Properties { self.residual.properties(temperature) } fn compute_max_density(&self, moles: &Array1) -> f64 { self.residual.compute_max_density(moles) } - fn contributions(&self) -> &[Box>] { + fn contributions(&self) -> &[Box>] { self.residual.contributions() } diff --git a/feos-core/src/equation_of_state/residual.rs b/feos-core/src/equation_of_state/residual.rs index 1d2ec5375..82b8d64b4 100644 --- a/feos-core/src/equation_of_state/residual.rs +++ b/feos-core/src/equation_of_state/residual.rs @@ -8,15 +8,16 @@ use num_traits::{One, Zero}; use std::marker::PhantomData; use std::ops::Div; -pub trait Properties { - type Values; -} - /// A reisdual Helmholtz energy model. pub trait Residual: Components + Send + Sync { type Properties; + /// For wrapper structs, the inner type. In basic cases `Self`. + type Inner: Residual; - fn properties>(&self, temperature: D) -> Self::Properties; + fn properties>( + &self, + temperature: D, + ) -> ::Properties; /// Return the maximum density in Angstrom^-3. /// @@ -28,7 +29,7 @@ pub trait Residual: Components + Send + Sync { /// Return a slice of the individual contributions (excluding the ideal gas) /// of the equation of state. - fn contributions(&self) -> &[Box>]; + fn contributions(&self) -> &[Box>]; /// Molar weight of all components. /// @@ -38,7 +39,8 @@ pub trait Residual: Components + Send + Sync { /// Evaluate the residual reduced Helmholtz energy $\beta A^\mathrm{res}$. fn evaluate_residual + Copy>(&self, state: &StateHD) -> D where - dyn HelmholtzEnergy: HelmholtzEnergyDual, D>, + dyn HelmholtzEnergy: + HelmholtzEnergyDual<::Properties, D>, { let properties = self.properties(state.temperature); self.contributions() @@ -54,7 +56,8 @@ pub trait Residual: Components + Send + Sync { state: &StateHD, ) -> Vec<(String, D)> where - dyn HelmholtzEnergy: HelmholtzEnergyDual, D>, + dyn HelmholtzEnergy: + HelmholtzEnergyDual<::Properties, D>, { let properties = self.properties(state.temperature); self.contributions() @@ -93,76 +96,76 @@ pub trait Residual: Components + Send + Sync { Ok(Density::from_reduced(self.compute_max_density(&mr))) } - // /// Calculate the second virial coefficient $B(T)$ - // fn second_virial_coefficient( - // &self, - // temperature: Temperature, - // moles: Option<&Moles>>, - // ) -> EosResult<>::Output> { - // let mr = self.validate_moles(moles)?; - // let x = (&mr / mr.sum()).into_value(); - // let mut rho = HyperDual64::zero(); - // rho.eps1 = 1.0; - // rho.eps2 = 1.0; - // let t = HyperDual64::from(temperature.to_reduced()); - // let s = StateHD::new_virial(t, rho, x); - // Ok(Quantity::from_reduced( - // self.evaluate_residual(&s).eps1eps2 * 0.5, - // )) - // } - - // /// Calculate the third virial coefficient $C(T)$ - // #[allow(clippy::type_complexity)] - // fn third_virial_coefficient( - // &self, - // temperature: Temperature, - // moles: Option<&Moles>>, - // ) -> EosResult<<>::Output as Div>::Output> { - // let mr = self.validate_moles(moles)?; - // let x = (&mr / mr.sum()).into_value(); - // let rho = Dual3_64::zero().derivative(); - // let t = Dual3_64::from(temperature.to_reduced()); - // let s = StateHD::new_virial(t, rho, x); - // Ok(Quantity::from_reduced(self.evaluate_residual(&s).v3 / 3.0)) - // } - - // /// Calculate the temperature derivative of the second virial coefficient $B'(T)$ - // #[allow(clippy::type_complexity)] - // fn second_virial_coefficient_temperature_derivative( - // &self, - // temperature: Temperature, - // moles: Option<&Moles>>, - // ) -> EosResult<<>::Output as Div>::Output> { - // let mr = self.validate_moles(moles)?; - // let x = (&mr / mr.sum()).into_value(); - // let mut rho = HyperDual::zero(); - // rho.eps1 = Dual64::one(); - // rho.eps2 = Dual64::one(); - // let t = HyperDual::from_re(Dual64::from(temperature.to_reduced()).derivative()); - // let s = StateHD::new_virial(t, rho, x); - // Ok(Quantity::from_reduced( - // self.evaluate_residual(&s).eps1eps2.eps * 0.5, - // )) - // } - - // /// Calculate the temperature derivative of the third virial coefficient $C'(T)$ - // #[allow(clippy::type_complexity)] - // fn third_virial_coefficient_temperature_derivative( - // &self, - // temperature: Temperature, - // moles: Option<&Moles>>, - // ) -> EosResult< - // <<>::Output as Div>::Output as Div>::Output, - // > { - // let mr = self.validate_moles(moles)?; - // let x = (&mr / mr.sum()).into_value(); - // let rho = Dual3::zero().derivative(); - // let t = Dual3::from_re(Dual64::from(temperature.to_reduced()).derivative()); - // let s = StateHD::new_virial(t, rho, x); - // Ok(Quantity::from_reduced( - // self.evaluate_residual(&s).v3.eps / 3.0, - // )) - // } + /// Calculate the second virial coefficient $B(T)$ + fn second_virial_coefficient( + &self, + temperature: Temperature, + moles: Option<&Moles>>, + ) -> EosResult<>::Output> { + let mr = self.validate_moles(moles)?; + let x = (&mr / mr.sum()).into_value(); + let mut rho = HyperDual64::zero(); + rho.eps1 = 1.0; + rho.eps2 = 1.0; + let t = HyperDual64::from(temperature.to_reduced()); + let s = StateHD::new_virial(t, rho, x); + Ok(Quantity::from_reduced( + self.evaluate_residual(&s).eps1eps2 * 0.5, + )) + } + + /// Calculate the third virial coefficient $C(T)$ + #[allow(clippy::type_complexity)] + fn third_virial_coefficient( + &self, + temperature: Temperature, + moles: Option<&Moles>>, + ) -> EosResult<<>::Output as Div>::Output> { + let mr = self.validate_moles(moles)?; + let x = (&mr / mr.sum()).into_value(); + let rho = Dual3_64::zero().derivative(); + let t = Dual3_64::from(temperature.to_reduced()); + let s = StateHD::new_virial(t, rho, x); + Ok(Quantity::from_reduced(self.evaluate_residual(&s).v3 / 3.0)) + } + + /// Calculate the temperature derivative of the second virial coefficient $B'(T)$ + #[allow(clippy::type_complexity)] + fn second_virial_coefficient_temperature_derivative( + &self, + temperature: Temperature, + moles: Option<&Moles>>, + ) -> EosResult<<>::Output as Div>::Output> { + let mr = self.validate_moles(moles)?; + let x = (&mr / mr.sum()).into_value(); + let mut rho = HyperDual::zero(); + rho.eps1 = Dual64::one(); + rho.eps2 = Dual64::one(); + let t = HyperDual::from_re(Dual64::from(temperature.to_reduced()).derivative()); + let s = StateHD::new_virial(t, rho, x); + Ok(Quantity::from_reduced( + self.evaluate_residual(&s).eps1eps2.eps * 0.5, + )) + } + + /// Calculate the temperature derivative of the third virial coefficient $C'(T)$ + #[allow(clippy::type_complexity)] + fn third_virial_coefficient_temperature_derivative( + &self, + temperature: Temperature, + moles: Option<&Moles>>, + ) -> EosResult< + <<>::Output as Div>::Output as Div>::Output, + > { + let mr = self.validate_moles(moles)?; + let x = (&mr / mr.sum()).into_value(); + let rho = Dual3::zero().derivative(); + let t = Dual3::from_re(Dual64::from(temperature.to_reduced()).derivative()); + let s = StateHD::new_virial(t, rho, x); + Ok(Quantity::from_reduced( + self.evaluate_residual(&s).v3.eps / 3.0, + )) + } } /// Reference values and residual entropy correlations for entropy scaling. @@ -205,8 +208,9 @@ impl Components for NoResidual { impl Residual for NoResidual { type Properties = PhantomData; + type Inner = Self; - fn properties>(&self, temperature: D) -> PhantomData { + fn properties>(&self, _: D) -> PhantomData { PhantomData } diff --git a/feos-core/src/lib.rs b/feos-core/src/lib.rs index 90ddeee49..5baeb764c 100644 --- a/feos-core/src/lib.rs +++ b/feos-core/src/lib.rs @@ -33,7 +33,7 @@ pub mod si; mod state; pub use equation_of_state::{ Components, DeBroglieWavelength, DeBroglieWavelengthDual, EntropyScaling, EquationOfState, - HelmholtzEnergy, HelmholtzEnergyDual, IdealGas, NoResidual, Properties, Residual, + HelmholtzEnergy, HelmholtzEnergyDual, IdealGas, NoResidual, Residual, }; pub use errors::{EosError, EosResult}; pub use phase_equilibria::{ diff --git a/feos-dft/src/functional.rs b/feos-dft/src/functional.rs index 913406f10..d95e7f2cd 100644 --- a/feos-dft/src/functional.rs +++ b/feos-dft/src/functional.rs @@ -15,6 +15,7 @@ use petgraph::graph::{Graph, UnGraph}; use petgraph::visit::EdgeRef; use petgraph::Directed; use std::borrow::Cow; +use std::marker::PhantomData; use std::ops::{Deref, MulAssign}; use std::sync::Arc; @@ -95,11 +96,18 @@ impl Components for DFT { } impl Residual for DFT { + type Properties = PhantomData; + type Inner = Self; + + fn properties>(&self, _: D) -> PhantomData { + PhantomData + } + fn compute_max_density(&self, moles: &Array1) -> f64 { self.0.compute_max_density(moles) } - fn contributions(&self) -> &[Box] { + fn contributions(&self) -> &[Box>] { unreachable!() } @@ -109,14 +117,16 @@ impl Residual for DFT { fn evaluate_residual + Copy>(&self, state: &StateHD) -> D where - dyn HelmholtzEnergy: HelmholtzEnergyDual, + dyn HelmholtzEnergy: HelmholtzEnergyDual, D>, { self.0 .contributions() .iter() - .map(|c| (c as &dyn HelmholtzEnergy).helmholtz_energy(state)) + .map(|c| (c as &dyn HelmholtzEnergy).helmholtz_energy(state, &PhantomData)) .sum::() - + self.ideal_chain_contribution().helmholtz_energy(state) + + self + .ideal_chain_contribution() + .helmholtz_energy(state, &PhantomData) } fn evaluate_residual_contributions + Copy>( @@ -124,7 +134,7 @@ impl Residual for DFT { state: &StateHD, ) -> Vec<(String, D)> where - dyn HelmholtzEnergy: HelmholtzEnergyDual, + dyn HelmholtzEnergy: HelmholtzEnergyDual, D>, { let mut res: Vec<(String, D)> = self .0 @@ -133,13 +143,14 @@ impl Residual for DFT { .map(|c| { ( c.to_string(), - (c as &dyn HelmholtzEnergy).helmholtz_energy(state), + (c as &dyn HelmholtzEnergy).helmholtz_energy(state, &PhantomData), ) }) .collect(); res.push(( self.ideal_chain_contribution().to_string(), - self.ideal_chain_contribution().helmholtz_energy(state), + self.ideal_chain_contribution() + .helmholtz_energy(state, &PhantomData), )); res } diff --git a/feos-dft/src/functional_contribution.rs b/feos-dft/src/functional_contribution.rs index 89bbf8499..02f6916b2 100644 --- a/feos-dft/src/functional_contribution.rs +++ b/feos-dft/src/functional_contribution.rs @@ -5,11 +5,18 @@ use ndarray::RemoveAxis; use num_dual::*; use num_traits::{One, Zero}; use std::fmt::Display; +use std::marker::PhantomData; macro_rules! impl_helmholtz_energy { ($number:ty) => { - impl HelmholtzEnergyDual<$number> for Box { - fn helmholtz_energy(&self, state: &StateHD<$number>) -> $number { + impl HelmholtzEnergyDual, $number> + for Box + { + fn helmholtz_energy( + &self, + state: &StateHD<$number>, + _: &PhantomData<$number>, + ) -> $number { // calculate weight functions let weight_functions = self.weight_functions(state.temperature); diff --git a/feos-dft/src/ideal_chain_contribution.rs b/feos-dft/src/ideal_chain_contribution.rs index d42e68a00..a38e6e544 100644 --- a/feos-dft/src/ideal_chain_contribution.rs +++ b/feos-dft/src/ideal_chain_contribution.rs @@ -3,6 +3,7 @@ use feos_core::{EosResult, HelmholtzEnergyDual, StateHD}; use ndarray::*; use num_dual::DualNum; use std::fmt; +use std::marker::PhantomData; #[derive(Clone)] pub struct IdealChainContribution { @@ -19,8 +20,8 @@ impl IdealChainContribution { } } -impl + Copy> HelmholtzEnergyDual for IdealChainContribution { - fn helmholtz_energy(&self, state: &StateHD) -> D { +impl + Copy> HelmholtzEnergyDual, D> for IdealChainContribution { + fn helmholtz_energy(&self, state: &StateHD, _: &PhantomData) -> D { let segments = self.component_index.len(); if self.component_index[segments - 1] + 1 != segments { return D::zero(); diff --git a/src/association/mod.rs b/src/association/mod.rs index c3c74dff0..4c8468e90 100644 --- a/src/association/mod.rs +++ b/src/association/mod.rs @@ -332,7 +332,7 @@ impl Association

{ } } -impl + Copy + ScalarOperand, P: HardSphereProperties> HelmholtzEnergyDual +impl + Copy + ScalarOperand, P: HardSphereProperties> HelmholtzEnergyDual for Association

{ fn helmholtz_energy(&self, state: &StateHD) -> D { diff --git a/src/hard_sphere/mod.rs b/src/hard_sphere/mod.rs index 62b494795..786b2c5f8 100644 --- a/src/hard_sphere/mod.rs +++ b/src/hard_sphere/mod.rs @@ -116,7 +116,9 @@ impl

HardSphere

{ } } -impl + Copy, P: HardSphereProperties> HelmholtzEnergyDual for HardSphere

{ +impl + Copy, P: HardSphereProperties> HelmholtzEnergyDual + for HardSphere

+{ fn helmholtz_energy(&self, state: &StateHD) -> D { let p = &self.parameters; let zeta = p.zeta(state.temperature, &state.partial_density, [0, 1, 2, 3]); diff --git a/src/pcsaft/eos/dispersion.rs b/src/pcsaft/eos/dispersion.rs index 416778685..75d67ce9a 100644 --- a/src/pcsaft/eos/dispersion.rs +++ b/src/pcsaft/eos/dispersion.rs @@ -1,4 +1,4 @@ -use super::PcSaftParameters; +use super::{PcSaft, PcSaftParameters}; use crate::hard_sphere::HardSphereProperties; use feos_core::{HelmholtzEnergyDual, StateHD}; use num_dual::DualNum; @@ -65,7 +65,7 @@ pub struct Dispersion { pub parameters: Arc, } -impl + Copy> HelmholtzEnergyDual for Dispersion { +impl + Copy> HelmholtzEnergyDual for Dispersion { fn helmholtz_energy(&self, state: &StateHD) -> D { // auxiliary variables let n = self.parameters.m.len(); diff --git a/src/pcsaft/eos/hard_chain.rs b/src/pcsaft/eos/hard_chain.rs index a046c893c..87adb1167 100644 --- a/src/pcsaft/eos/hard_chain.rs +++ b/src/pcsaft/eos/hard_chain.rs @@ -1,4 +1,4 @@ -use super::PcSaftParameters; +use super::{PcSaft, PcSaftParameters}; use crate::hard_sphere::HardSphereProperties; use feos_core::{HelmholtzEnergyDual, StateHD}; use ndarray::Array; @@ -10,7 +10,7 @@ pub struct HardChain { pub parameters: Arc, } -impl + Copy> HelmholtzEnergyDual for HardChain { +impl + Copy> HelmholtzEnergyDual for HardChain { fn helmholtz_energy(&self, state: &StateHD) -> D { let p = &self.parameters; let d = self.parameters.hs_diameter(state.temperature); diff --git a/src/pcsaft/eos/mod.rs b/src/pcsaft/eos/mod.rs index af2ad08d7..358477b28 100644 --- a/src/pcsaft/eos/mod.rs +++ b/src/pcsaft/eos/mod.rs @@ -9,6 +9,7 @@ use feos_core::{ use ndarray::Array1; use std::f64::consts::{FRAC_PI_6, PI}; use std::fmt; +use std::marker::PhantomData; use std::sync::Arc; use typenum::P2; @@ -44,7 +45,7 @@ impl Default for PcSaftOptions { pub struct PcSaft { parameters: Arc, options: PcSaftOptions, - contributions: Vec>, + contributions: Vec>>, } impl PcSaft { @@ -108,13 +109,20 @@ impl Components for PcSaft { } impl Residual for PcSaft { + type Properties = PhantomData; + type Inner = Self; + + fn properties>(&self, temperature: D) -> PhantomData { + todo!() + } + fn compute_max_density(&self, moles: &Array1) -> f64 { self.options.max_eta * moles.sum() / (FRAC_PI_6 * &self.parameters.m * self.parameters.sigma.mapv(|v| v.powi(3)) * moles) .sum() } - fn contributions(&self) -> &[Box] { + fn contributions(&self) -> &[Box>] { &self.contributions } diff --git a/src/pcsaft/eos/polar.rs b/src/pcsaft/eos/polar.rs index 9efe1c428..b6129afc0 100644 --- a/src/pcsaft/eos/polar.rs +++ b/src/pcsaft/eos/polar.rs @@ -1,4 +1,4 @@ -use super::PcSaftParameters; +use super::{PcSaft, PcSaftParameters}; use crate::hard_sphere::HardSphereProperties; use feos_core::{HelmholtzEnergyDual, StateHD}; use ndarray::prelude::*; @@ -165,7 +165,7 @@ pub struct Dipole { pub parameters: Arc, } -impl + Copy> HelmholtzEnergyDual for Dipole { +impl + Copy> HelmholtzEnergyDual for Dipole { fn helmholtz_energy(&self, state: &StateHD) -> D { let m = MeanSegmentNumbers::new(&self.parameters, Multipole::Dipole); let p = &self.parameters; @@ -242,7 +242,7 @@ pub struct Quadrupole { pub parameters: Arc, } -impl + Copy> HelmholtzEnergyDual for Quadrupole { +impl + Copy> HelmholtzEnergyDual for Quadrupole { fn helmholtz_energy(&self, state: &StateHD) -> D { let m = MeanSegmentNumbers::new(&self.parameters, Multipole::Quadrupole); let p = &self.parameters; @@ -329,7 +329,7 @@ pub struct DipoleQuadrupole { pub variant: DQVariants, } -impl + Copy> HelmholtzEnergyDual for DipoleQuadrupole { +impl + Copy> HelmholtzEnergyDual for DipoleQuadrupole { fn helmholtz_energy(&self, state: &StateHD) -> D { let p = &self.parameters;