feos_ad/core/
total.rs

1use super::{ParametersAD, ResidualHelmholtzEnergy};
2use nalgebra::SVector;
3use num_dual::{
4    first_derivative, gradient, hessian, second_derivative, Dual, Dual2, Dual2Vec, DualNum, DualVec,
5};
6
7/// Implementation of an ideal gas Helmholtz energy contribution.
8pub trait IdealGasAD: ParametersAD {
9    /// The name of the model.
10    const IDEAL_GAS: &str;
11
12    /// The logarithmic cubed thermal de Broglie wavelength for the given temperature.
13    fn ln_lambda3_dual<D: DualNum<f64> + Copy>(
14        parameters: &Self::Parameters<D>,
15        temperature: D,
16    ) -> D;
17}
18
19/// An equation of state consisting of a residual model and an ideal gas model.
20pub struct EquationOfStateAD<I, R, const N: usize> {
21    ideal_gas: [I; N],
22    residual: R,
23}
24
25impl<I, R, const N: usize> EquationOfStateAD<I, R, N> {
26    pub fn new(ideal_gas: [I; N], residual: R) -> Self {
27        Self {
28            ideal_gas,
29            residual,
30        }
31    }
32}
33
34impl<I: ParametersAD, R: ParametersAD, const N: usize> ParametersAD for EquationOfStateAD<I, R, N> {
35    type Parameters<D: DualNum<f64> + Copy> = ([I::Parameters<D>; N], R::Parameters<D>);
36
37    fn params<D: DualNum<f64> + Copy>(&self) -> Self::Parameters<D> {
38        (
39            self.ideal_gas.each_ref().map(I::params),
40            self.residual.params(),
41        )
42    }
43
44    fn params_from_inner<D: DualNum<f64> + Copy, D2: DualNum<f64, Inner = D> + Copy>(
45        (ideal_gas, residual): &([I::Parameters<D>; N], R::Parameters<D>),
46    ) -> Self::Parameters<D2> {
47        (
48            ideal_gas.each_ref().map(|ig| I::params_from_inner(ig)),
49            R::params_from_inner(residual),
50        )
51    }
52}
53
54impl<I: ParametersAD, R: ResidualHelmholtzEnergy<N>, const N: usize> ResidualHelmholtzEnergy<N>
55    for EquationOfStateAD<I, R, N>
56{
57    const RESIDUAL: &str = R::RESIDUAL;
58
59    fn compute_max_density(&self, molefracs: &SVector<f64, N>) -> f64 {
60        self.residual.compute_max_density(molefracs)
61    }
62
63    fn residual_helmholtz_energy_density<D: DualNum<f64> + Copy>(
64        (_, residual): &([I::Parameters<D>; N], R::Parameters<D>),
65        temperature: D,
66        partial_density: &SVector<D, N>,
67    ) -> D {
68        R::residual_helmholtz_energy_density(residual, temperature, partial_density)
69    }
70}
71
72/// Methods of [EquationOfStateAD] extracted in a trait for genericness.
73pub trait TotalHelmholtzEnergy<const N: usize>: ResidualHelmholtzEnergy<N> {
74    const IDEAL_GAS: &str;
75
76    fn ln_lambda3<D: DualNum<f64> + Copy>(
77        parameters: &Self::Parameters<D>,
78        temperature: D,
79    ) -> SVector<D, N>;
80
81    fn helmholtz_energy_density<D: DualNum<f64> + Copy>(
82        parameters: &Self::Parameters<D>,
83        temperature: D,
84        partial_density: &SVector<D, N>,
85    ) -> D {
86        let ln_lambda_3 = Self::ln_lambda3(parameters, temperature);
87        let ig = partial_density
88            .component_mul(
89                &(partial_density.map(|d| d.ln()) + ln_lambda_3 - SVector::from([D::from(1.0); N])),
90            )
91            .sum()
92            * temperature;
93        Self::residual_helmholtz_energy_density(parameters, temperature, partial_density) + ig
94    }
95
96    fn molar_helmholtz_energy<D: DualNum<f64> + Copy>(
97        parameters: &Self::Parameters<D>,
98        temperature: D,
99        molar_volume: D,
100        molefracs: &SVector<D, N>,
101    ) -> D {
102        let partial_density = molefracs / molar_volume;
103        Self::helmholtz_energy_density(parameters, temperature, &partial_density) * molar_volume
104    }
105
106    fn chemical_potential<D: DualNum<f64> + Copy>(
107        parameters: &Self::Parameters<D>,
108        temperature: D,
109        molar_volume: D,
110        molefracs: &SVector<D, N>,
111    ) -> SVector<D, N> {
112        let params = Self::params_from_inner(parameters);
113        let temperature = DualVec::from_re(temperature);
114        let molar_volume = DualVec::from_re(molar_volume);
115        let (_, mu) = gradient(
116            |molefracs| {
117                Self::molar_helmholtz_energy(&params, temperature, molar_volume, &molefracs)
118            },
119            *molefracs,
120        );
121        mu
122    }
123
124    fn molar_entropy<D: DualNum<f64> + Copy>(
125        parameters: &Self::Parameters<D>,
126        temperature: D,
127        molar_volume: D,
128        molefracs: &SVector<D, N>,
129    ) -> D {
130        let params = Self::params_from_inner(parameters);
131        let molar_volume = Dual::from_re(molar_volume);
132        let molefracs = molefracs.map(Dual::from_re);
133        let (_, da_dt) = first_derivative(
134            |temperature| {
135                Self::molar_helmholtz_energy(&params, temperature, molar_volume, &molefracs)
136            },
137            temperature,
138        );
139        -da_dt
140    }
141
142    fn molar_enthalpy<D: DualNum<f64> + Copy>(
143        parameters: &Self::Parameters<D>,
144        temperature: D,
145        molar_volume: D,
146        molefracs: &SVector<D, N>,
147    ) -> D {
148        let params = Self::params_from_inner(parameters);
149        let molefracs = molefracs.map(DualVec::from_re);
150        let (a, da) = gradient(
151            |x| {
152                let [temperature, molar_volume] = x.data.0[0];
153                Self::molar_helmholtz_energy(&params, temperature, molar_volume, &molefracs)
154            },
155            SVector::from([temperature, molar_volume]),
156        );
157        let [da_dt, da_dv] = da.data.0[0];
158        a - temperature * da_dt - molar_volume * da_dv
159    }
160
161    fn molar_isochoric_heat_capacity<D: DualNum<f64> + Copy>(
162        parameters: &Self::Parameters<D>,
163        temperature: D,
164        molar_volume: D,
165        molefracs: &SVector<D, N>,
166    ) -> D {
167        let params = Self::params_from_inner(parameters);
168        let molar_volume = Dual2::from_re(molar_volume);
169        let molefracs = molefracs.map(Dual2::from_re);
170        let (_, _, d2a) = second_derivative(
171            |temperature| {
172                Self::molar_helmholtz_energy(&params, temperature, molar_volume, &molefracs)
173            },
174            temperature,
175        );
176        -temperature * d2a
177    }
178
179    fn molar_isobaric_heat_capacity<D: DualNum<f64> + Copy>(
180        parameters: &Self::Parameters<D>,
181        temperature: D,
182        molar_volume: D,
183        molefracs: &SVector<D, N>,
184    ) -> D {
185        let params = Self::params_from_inner(parameters);
186        let molefracs = molefracs.map(Dual2Vec::from_re);
187        let (_, _, d2a) = hessian(
188            |x| {
189                let [temperature, molar_volume] = x.data.0[0];
190                Self::molar_helmholtz_energy(&params, temperature, molar_volume, &molefracs)
191            },
192            SVector::from([temperature, molar_volume]),
193        );
194        let [[a_tt, a_tv], [_, a_vv]] = d2a.data.0;
195        temperature * (a_tv * a_tv / a_vv - a_tt)
196    }
197
198    fn pressure_entropy<D: DualNum<f64> + Copy>(
199        parameters: &Self::Parameters<D>,
200        temperature: D,
201        molar_volume: D,
202        molefracs: &SVector<D, N>,
203    ) -> SVector<D, 2> {
204        let params = Self::params_from_inner(parameters);
205        let molefracs = molefracs.map(DualVec::from_re);
206        gradient(
207            |x| {
208                let [molar_volume, temperature] = x.data.0[0];
209                -Self::molar_helmholtz_energy(&params, temperature, molar_volume, &molefracs)
210            },
211            SVector::from([molar_volume, temperature]),
212        )
213        .1
214    }
215
216    fn pressure_enthalpy<D: DualNum<f64> + Copy>(
217        parameters: &Self::Parameters<D>,
218        temperature: D,
219        molar_volume: D,
220        molefracs: &SVector<D, N>,
221    ) -> SVector<D, 2> {
222        let params = Self::params_from_inner(parameters);
223        let molefracs = molefracs.map(DualVec::from_re);
224        let (a, da) = gradient(
225            |x| {
226                let [temperature, molar_volume] = x.data.0[0];
227                Self::molar_helmholtz_energy(&params, temperature, molar_volume, &molefracs)
228            },
229            SVector::from([temperature, molar_volume]),
230        );
231        let [da_dt, da_dv] = da.data.0[0];
232        let h = a - temperature * da_dt - molar_volume * da_dv;
233        let p = -da_dv;
234        SVector::from([p, h])
235    }
236}
237
238impl<I: IdealGasAD, R: ResidualHelmholtzEnergy<N>, const N: usize> TotalHelmholtzEnergy<N>
239    for EquationOfStateAD<I, R, N>
240{
241    const IDEAL_GAS: &str = I::IDEAL_GAS;
242
243    fn ln_lambda3<D: DualNum<f64> + Copy>(
244        (ideal_gas, _): &([I::Parameters<D>; N], R::Parameters<D>),
245        temperature: D,
246    ) -> SVector<D, N> {
247        SVector::from(
248            ideal_gas
249                .each_ref()
250                .map(|ig| I::ln_lambda3_dual(ig, temperature)),
251        )
252    }
253}
feos_ad/core/total.rs

feos_ad/core/
total.rs
