feos_ad/core/
residual.rs

1use super::{FeOsWrapper, HelmholtzEnergyWrapper};
2use nalgebra::{SMatrix, SVector};
3use num_dual::{
4    first_derivative, gradient, hessian, partial_hessian, second_derivative, Dual, Dual2, Dual2Vec,
5    DualNum, DualVec, HyperDualVec,
6};
7use std::sync::Arc;
8
9/// A model that can be evaluated with derivatives of its parameters.
10pub trait ParametersAD: Send + Sync + Sized {
11    /// The type of the structure that stores the parameters internally.
12    type Parameters<D: DualNum<f64> + Copy>: Clone;
13
14    /// Return the parameters in the given data type.
15    fn params<D: DualNum<f64> + Copy>(&self) -> Self::Parameters<D>;
16
17    /// Lift the parameters to the given type of dual number.
18    fn params_from_inner<D: DualNum<f64> + Copy, D2: DualNum<f64, Inner = D> + Copy>(
19        parameters: &Self::Parameters<D>,
20    ) -> Self::Parameters<D2>;
21
22    /// Wraps the model in the [HelmholtzEnergyWrapper] struct, so that it can be used
23    /// as an argument to [StateAD](crate::StateAD) and [PhaseEquilibriumAD](crate::PhaseEquilibriumAD) constructors.
24    fn wrap<const N: usize>(self) -> HelmholtzEnergyWrapper<Self, f64, N> {
25        let parameters = self.params();
26        HelmholtzEnergyWrapper {
27            eos: Arc::new(FeOsWrapper(self)),
28            parameters,
29        }
30    }
31}
32
33/// Implementation of a residual Helmholtz energy model.
34pub trait ResidualHelmholtzEnergy<const N: usize>: ParametersAD {
35    /// The name of the model.
36    const RESIDUAL: &str;
37
38    /// Return a density (in reduced units) that corresponds to a dense liquid phase.
39    fn compute_max_density(&self, molefracs: &SVector<f64, N>) -> f64;
40
41    fn residual_helmholtz_energy_density<D: DualNum<f64> + Copy>(
42        parameters: &Self::Parameters<D>,
43        temperature: D,
44        partial_density: &SVector<D, N>,
45    ) -> D;
46
47    fn residual_molar_helmholtz_energy<D: DualNum<f64> + Copy>(
48        parameters: &Self::Parameters<D>,
49        temperature: D,
50        molar_volume: D,
51        molefracs: &SVector<D, N>,
52    ) -> D {
53        let partial_density = molefracs / molar_volume;
54        Self::residual_helmholtz_energy_density(parameters, temperature, &partial_density)
55            * molar_volume
56    }
57
58    fn residual_chemical_potential<D: DualNum<f64> + Copy>(
59        parameters: &Self::Parameters<D>,
60        temperature: D,
61        molar_volume: D,
62        molefracs: &SVector<D, N>,
63    ) -> SVector<D, N> {
64        let params = Self::params_from_inner(parameters);
65        let temperature = DualVec::from_re(temperature);
66        let molar_volume = DualVec::from_re(molar_volume);
67        let (_, mu) = gradient(
68            |molefracs| {
69                Self::residual_molar_helmholtz_energy(
70                    &params,
71                    temperature,
72                    molar_volume,
73                    &molefracs,
74                )
75            },
76            *molefracs,
77        );
78        mu
79    }
80
81    fn pressure<D: DualNum<f64> + Copy>(
82        parameters: &Self::Parameters<D>,
83        temperature: D,
84        molar_volume: D,
85        molefracs: &SVector<D, N>,
86    ) -> D {
87        let params = Self::params_from_inner(parameters);
88        let t = Dual::from_re(temperature);
89        let molefracs = molefracs.map(Dual::from_re);
90        let (_, dadv) = first_derivative(
91            |v| Self::residual_molar_helmholtz_energy(&params, t, v, &molefracs),
92            molar_volume,
93        );
94        -dadv + temperature / molar_volume
95    }
96
97    fn residual_molar_entropy<D: DualNum<f64> + Copy>(
98        parameters: &Self::Parameters<D>,
99        temperature: D,
100        molar_volume: D,
101        molefracs: &SVector<D, N>,
102    ) -> D {
103        let params = Self::params_from_inner(parameters);
104        let molar_volume = Dual::from_re(molar_volume);
105        let molefracs = molefracs.map(Dual::from_re);
106        let (_, da_dt) = first_derivative(
107            |temperature| {
108                Self::residual_molar_helmholtz_energy(
109                    &params,
110                    temperature,
111                    molar_volume,
112                    &molefracs,
113                )
114            },
115            temperature,
116        );
117        -da_dt
118    }
119
120    fn residual_molar_enthalpy<D: DualNum<f64> + Copy>(
121        parameters: &Self::Parameters<D>,
122        temperature: D,
123        molar_volume: D,
124        molefracs: &SVector<D, N>,
125    ) -> D {
126        let params = Self::params_from_inner(parameters);
127        let molefracs = molefracs.map(DualVec::from_re);
128        let (a, da) = gradient(
129            |x| {
130                let [temperature, molar_volume] = x.data.0[0];
131                Self::residual_molar_helmholtz_energy(
132                    &params,
133                    temperature,
134                    molar_volume,
135                    &molefracs,
136                )
137            },
138            SVector::from([temperature, molar_volume]),
139        );
140        let [da_dt, da_dv] = da.data.0[0];
141        a - temperature * da_dt - molar_volume * da_dv
142    }
143
144    fn dp_drho<D: DualNum<f64> + Copy>(
145        parameters: &Self::Parameters<D>,
146        temperature: D,
147        molar_volume: D,
148        molefracs: &SVector<D, N>,
149    ) -> (D, D, D) {
150        let params = Self::params_from_inner(parameters);
151        let t = Dual2::from_re(temperature);
152        let x = molefracs.map(Dual2::from_re);
153        let (a, da, d2a) = second_derivative(
154            |molar_volume| Self::residual_molar_helmholtz_energy(&params, t, molar_volume, &x),
155            molar_volume,
156        );
157        let density = molar_volume.recip();
158        (
159            a * density,
160            -da + temperature * density,
161            molar_volume * molar_volume * d2a + temperature,
162        )
163    }
164
165    /// calculates p, mu_res, dp_drho, dmu_drho
166    fn dmu_drho<D: DualNum<f64> + Copy>(
167        parameters: &Self::Parameters<D>,
168        temperature: D,
169        partial_density: &SVector<D, N>,
170    ) -> (D, SVector<D, N>, SVector<D, N>, SMatrix<D, N, N>) {
171        let params = Self::params_from_inner(parameters);
172        let t = Dual2Vec::from_re(temperature);
173        let (f_res, mu_res, dmu_res) = hessian(
174            |rho| Self::residual_helmholtz_energy_density(&params, t, &rho),
175            *partial_density,
176        );
177        let p = mu_res.dot(partial_density) - f_res + temperature * partial_density.sum();
178        let dmu = dmu_res + SMatrix::from_diagonal(&partial_density.map(|d| temperature / d));
179        let dp = dmu * partial_density;
180        (p, mu_res, dp, dmu)
181    }
182
183    /// calculates p, mu_res, dp_dv, dmu_dv
184    fn dmu_dv<D: DualNum<f64> + Copy>(
185        parameters: &Self::Parameters<D>,
186        temperature: D,
187        molar_volume: D,
188        molefracs: &SVector<D, N>,
189    ) -> (D, SVector<D, N>, D, SVector<D, N>) {
190        let params = Self::params_from_inner(parameters);
191        let t = HyperDualVec::from_re(temperature);
192        let (_, mu_res, a_res_v, mu_res_v) = partial_hessian(
193            |x, v| Self::residual_molar_helmholtz_energy(&params, t, v[0], &x),
194            *molefracs,
195            SVector::from([molar_volume]),
196        );
197        let p = (-a_res_v)[0] + temperature / molar_volume;
198        let mu_v = mu_res_v.map(|m| m - temperature / molar_volume);
199        let p_v = mu_v.dot(molefracs) / molar_volume;
200        (p, mu_res, p_v, mu_v)
201    }
202}
feos_ad/core/residual.rs

feos_ad/core/
residual.rs
