feos_core/equation_of_state/
mod.rs

1use crate::ReferenceSystem;
2use crate::state::StateHD;
3use nalgebra::{
4    Const, DVector, DefaultAllocator, Dim, Dyn, OVector, SVector, allocator::Allocator,
5};
6use num_dual::DualNum;
7use quantity::{Dimensionless, MolarEnergy, MolarVolume, Temperature};
8use std::ops::Deref;
9
10mod residual;
11pub use residual::{EntropyScaling, Molarweight, NoResidual, Residual, ResidualDyn, Subset};
12
13/// An equation of state consisting of an ideal gas model
14/// and a residual Helmholtz energy model.
15#[derive(Clone)]
16pub struct EquationOfState<I, R> {
17    pub ideal_gas: I,
18    pub residual: R,
19}
20
21impl<I, R> Deref for EquationOfState<I, R> {
22    type Target = R;
23    fn deref(&self) -> &R {
24        &self.residual
25    }
26}
27
28impl<I, R> EquationOfState<I, R> {
29    /// Return a new [EquationOfState] with the given ideal gas
30    /// and residual models.
31    pub fn new(ideal_gas: I, residual: R) -> Self {
32        Self {
33            ideal_gas,
34            residual,
35        }
36    }
37}
38
39impl<I> EquationOfState<Vec<I>, NoResidual> {
40    /// Return a new [EquationOfState] that only consists of
41    /// an ideal gas models.
42    pub fn ideal_gas(ideal_gas: Vec<I>) -> Self {
43        let residual = NoResidual(ideal_gas.len());
44        Self {
45            ideal_gas,
46            residual,
47        }
48    }
49}
50
51impl<I, R: ResidualDyn> ResidualDyn for EquationOfState<Vec<I>, R> {
52    fn components(&self) -> usize {
53        self.residual.components()
54    }
55
56    fn compute_max_density<D: DualNum<f64> + Copy>(&self, molefracs: &DVector<D>) -> D {
57        self.residual.compute_max_density(molefracs)
58    }
59
60    fn reduced_helmholtz_energy_density_contributions<D: DualNum<f64> + Copy>(
61        &self,
62        state: &StateHD<D>,
63    ) -> Vec<(&'static str, D)> {
64        self.residual
65            .reduced_helmholtz_energy_density_contributions(state)
66    }
67}
68
69impl<I: Clone, R: Subset> Subset for EquationOfState<Vec<I>, R> {
70    fn subset(&self, component_list: &[usize]) -> Self {
71        let ideal_gas = component_list
72            .iter()
73            .map(|&i| self.ideal_gas[i].clone())
74            .collect();
75        EquationOfState {
76            ideal_gas,
77            residual: self.residual.subset(component_list),
78        }
79    }
80}
81
82impl<I: Clone, R: Residual<Const<N>, D>, D: DualNum<f64> + Copy, const N: usize>
83    Residual<Const<N>, D> for EquationOfState<[I; N], R>
84{
85    fn components(&self) -> usize {
86        N
87    }
88
89    type Real = EquationOfState<[I; N], R::Real>;
90    type Lifted<D2: DualNum<f64, Inner = D> + Copy> = EquationOfState<[I; N], R::Lifted<D2>>;
91    fn re(&self) -> Self::Real {
92        EquationOfState::new(self.ideal_gas.clone(), self.residual.re())
93    }
94    fn lift<D2: DualNum<f64, Inner = D> + Copy>(&self) -> Self::Lifted<D2> {
95        EquationOfState::new(self.ideal_gas.clone(), self.residual.lift())
96    }
97
98    fn compute_max_density(&self, molefracs: &SVector<D, N>) -> D {
99        self.residual.compute_max_density(molefracs)
100    }
101
102    fn reduced_helmholtz_energy_density_contributions(
103        &self,
104        state: &StateHD<D, Const<N>>,
105    ) -> Vec<(&'static str, D)> {
106        self.residual
107            .reduced_helmholtz_energy_density_contributions(state)
108    }
109
110    fn reduced_residual_helmholtz_energy_density(&self, state: &StateHD<D, Const<N>>) -> D {
111        self.residual
112            .reduced_residual_helmholtz_energy_density(state)
113    }
114}
115
116/// Ideal gas Helmholtz energy contribution.
117pub trait IdealGas {
118    /// Implementation of an ideal gas model in terms of the
119    /// logarithm of the cubic thermal de Broglie wavelength
120    /// in units ln(A³) for each component in the system.
121    fn ln_lambda3<D: DualNum<f64> + Copy>(&self, temperature: D) -> D;
122
123    /// The name of the ideal gas model.
124    fn ideal_gas_model(&self) -> &'static str;
125}
126
127/// Ideal gas Helmholtz energy contribution with automatic differentiation with
128/// respect to parameters.
129pub trait IdealGasAD<D = f64>: Clone {
130    type Real: IdealGasAD;
131    type Lifted<D2: DualNum<f64, Inner = D> + Copy>: IdealGasAD<D2>;
132    fn re(&self) -> Self::Real;
133    fn lift<D2: DualNum<f64, Inner = D> + Copy>(&self) -> Self::Lifted<D2>;
134
135    /// Implementation of an ideal gas model in terms of the
136    /// logarithm of the cubic thermal de Broglie wavelength
137    /// in units ln(A³) for each component in the system.
138    fn ln_lambda3(&self, temperature: D) -> D;
139
140    /// The name of the ideal gas model.
141    fn ideal_gas_model(&self) -> &'static str;
142}
143
144/// A total Helmholtz energy model consisting of a [Residual] model and an [IdealGas] part.
145pub trait Total<N: Dim = Dyn, D: DualNum<f64> + Copy = f64>: Residual<N, D>
146where
147    DefaultAllocator: Allocator<N>,
148{
149    type RealTotal: Total<N, f64>;
150    type LiftedTotal<D2: DualNum<f64, Inner = D> + Copy>: Total<N, D2>;
151    fn re_total(&self) -> Self::RealTotal;
152    fn lift_total<D2: DualNum<f64, Inner = D> + Copy>(&self) -> Self::LiftedTotal<D2>;
153
154    fn ideal_gas_model(&self) -> &'static str;
155
156    fn ln_lambda3(&self, temperature: D) -> OVector<D, N>;
157
158    fn ideal_gas_molar_helmholtz_energy(
159        &self,
160        temperature: D,
161        molar_volume: D,
162        molefracs: &OVector<D, N>,
163    ) -> D {
164        let partial_density = molefracs / molar_volume;
165        let mut res = D::from(0.0);
166        for (&l, &r) in self
167            .ln_lambda3(temperature)
168            .iter()
169            .zip(partial_density.iter())
170        {
171            let ln_rho_m1 = if r.re() == 0.0 {
172                D::from(0.0)
173            } else {
174                r.ln() - 1.0
175            };
176            res += r * (l + ln_rho_m1)
177        }
178        res * molar_volume * temperature
179    }
180
181    fn ideal_gas_helmholtz_energy(
182        &self,
183        temperature: Temperature<D>,
184        volume: MolarVolume<D>,
185        moles: &OVector<D, N>,
186    ) -> MolarEnergy<D> {
187        let total_moles = moles.sum();
188        let molefracs = moles / total_moles;
189        let molar_volume = volume.into_reduced() / total_moles;
190        MolarEnergy::from_reduced(self.ideal_gas_molar_helmholtz_energy(
191            temperature.into_reduced(),
192            molar_volume,
193            &molefracs,
194        )) * Dimensionless::new(total_moles)
195    }
196}
197
198impl<
199    I: IdealGas + 'static,
200    C: Deref<Target = EquationOfState<Vec<I>, R>> + Clone,
201    R: ResidualDyn + 'static,
202    D: DualNum<f64> + Copy,
203> Total<Dyn, D> for C
204{
205    type RealTotal = Self;
206    type LiftedTotal<D2: DualNum<f64, Inner = D> + Copy> = Self;
207    fn re_total(&self) -> Self::RealTotal {
208        self.clone()
209    }
210    fn lift_total<D2: DualNum<f64, Inner = D> + Copy>(&self) -> Self::LiftedTotal<D2> {
211        self.clone()
212    }
213
214    fn ideal_gas_model(&self) -> &'static str {
215        self.ideal_gas[0].ideal_gas_model()
216    }
217
218    fn ln_lambda3(&self, temperature: D) -> DVector<D> {
219        DVector::from_vec(
220            self.ideal_gas
221                .iter()
222                .map(|i| i.ln_lambda3(temperature))
223                .collect(),
224        )
225    }
226}
227
228impl<I: IdealGasAD<D>, R: Residual<Const<N>, D>, D: DualNum<f64> + Copy, const N: usize>
229    Total<Const<N>, D> for EquationOfState<[I; N], R>
230{
231    type RealTotal = EquationOfState<[I::Real; N], R::Real>;
232    type LiftedTotal<D2: DualNum<f64, Inner = D> + Copy> =
233        EquationOfState<[I::Lifted<D2>; N], R::Lifted<D2>>;
234    fn re_total(&self) -> Self::RealTotal {
235        EquationOfState::new(
236            self.ideal_gas.each_ref().map(|i| i.re()),
237            self.residual.re(),
238        )
239    }
240    fn lift_total<D2: DualNum<f64, Inner = D> + Copy>(&self) -> Self::LiftedTotal<D2> {
241        EquationOfState::new(
242            self.ideal_gas.each_ref().map(|i| i.lift()),
243            self.residual.lift(),
244        )
245    }
246
247    fn ideal_gas_model(&self) -> &'static str {
248        self.ideal_gas[0].ideal_gas_model()
249    }
250
251    fn ln_lambda3(&self, temperature: D) -> SVector<D, N> {
252        SVector::from(self.ideal_gas.each_ref().map(|i| i.ln_lambda3(temperature)))
253    }
254}
feos_core/equation_of_state/mod.rs

feos_core/equation_of_state/
mod.rs
