feos_core/ad/
mod.rs

1use crate::DensityInitialization::Liquid;
2use crate::density_iteration::density_iteration;
3use crate::{FeosResult, PhaseEquilibrium, ReferenceSystem, Residual};
4use nalgebra::{Const, SVector, U1, U2};
5#[cfg(feature = "rayon")]
6use ndarray::{Array1, Array2, ArrayView2, Zip};
7use num_dual::{Derivative, DualSVec, DualStruct};
8use quantity::{Density, Pressure, Temperature};
9#[cfg(feature = "rayon")]
10use quantity::{KELVIN, KILO, METER, MOL, PASCAL};
11
12type Gradient<const P: usize> = DualSVec<f64, f64, P>;
13
14/// A model that can be evaluated with derivatives of its parameters.
15pub trait ParametersAD<const N: usize>: for<'a> From<&'a [f64]> + Residual<Const<N>> {
16    /// Return a mutable reference to the parameter named by `index` from the parameter set.
17    fn index_parameters_mut<'a, const P: usize>(
18        eos: &'a mut Self::Lifted<Gradient<P>>,
19        index: &str,
20    ) -> &'a mut Gradient<P>;
21
22    /// Return the parameters with the appropriate derivatives.
23    fn named_derivatives<const P: usize>(
24        &self,
25        parameter_names: [&str; P],
26    ) -> Self::Lifted<Gradient<P>> {
27        let mut eos = self.lift::<Gradient<P>>();
28        for (i, p) in parameter_names.into_iter().enumerate() {
29            Self::index_parameters_mut(&mut eos, p).eps =
30                Derivative::derivative_generic(Const::<P>, U1, i)
31        }
32        eos
33    }
34}
35
36/// Properties that can be evaluated with derivatives of model parameters.
37pub trait PropertiesAD {
38    fn vapor_pressure<const P: usize>(
39        &self,
40        temperature: Temperature,
41    ) -> FeosResult<Pressure<Gradient<P>>>
42    where
43        Self: Residual<U1, Gradient<P>>,
44    {
45        let eos_f64 = self.re();
46        let (_, [vapor_density, liquid_density]) =
47            PhaseEquilibrium::pure_t(&eos_f64, temperature, None, Default::default())?;
48
49        // implicit differentiation is implemented here instead of just calling pure_t with dual
50        // numbers, because for the first derivative, we can avoid calculating density derivatives.
51        let v1 = 1.0 / liquid_density.to_reduced();
52        let v2 = 1.0 / vapor_density.to_reduced();
53        let t = temperature.into_reduced();
54        let (a1, a2) = {
55            let t = Gradient::from(t);
56            let v1 = Gradient::from(v1);
57            let v2 = Gradient::from(v2);
58            let x = Self::pure_molefracs();
59
60            let a1 = self.residual_molar_helmholtz_energy(t, v1, &x);
61            let a2 = self.residual_molar_helmholtz_energy(t, v2, &x);
62            (a1, a2)
63        };
64
65        let p = -(a1 - a2 + t * (v2 / v1).ln()) / (v1 - v2);
66        Ok(Pressure::from_reduced(p))
67    }
68
69    fn equilibrium_liquid_density<const P: usize>(
70        &self,
71        temperature: Temperature,
72    ) -> FeosResult<(Pressure<Gradient<P>>, Density<Gradient<P>>)>
73    where
74        Self: Residual<U1, Gradient<P>>,
75    {
76        let t = Temperature::from_inner(&temperature);
77        PhaseEquilibrium::pure_t(self, t, None, Default::default()).map(|(p, [_, rho])| (p, rho))
78    }
79
80    fn liquid_density<const P: usize>(
81        &self,
82        temperature: Temperature,
83        pressure: Pressure,
84    ) -> FeosResult<Density<Gradient<P>>>
85    where
86        Self: Residual<U1, Gradient<P>>,
87    {
88        let x = Self::pure_molefracs();
89        let t = Temperature::from_inner(&temperature);
90        let p = Pressure::from_inner(&pressure);
91        density_iteration(self, t, p, &x, Some(Liquid))
92    }
93
94    #[cfg(feature = "rayon")]
95    fn vapor_pressure_parallel<const P: usize>(
96        parameter_names: [String; P],
97        parameters: ArrayView2<f64>,
98        input: ArrayView2<f64>,
99    ) -> (Array1<f64>, Array2<f64>, Array1<bool>)
100    where
101        Self: ParametersAD<1>,
102    {
103        parallelize::<_, Self, _, _>(
104            parameter_names,
105            parameters,
106            input,
107            |eos: &Self::Lifted<Gradient<P>>, inp| {
108                eos.vapor_pressure(inp[0] * KELVIN)
109                    .map(|p| p.convert_into(PASCAL))
110            },
111        )
112    }
113
114    #[cfg(feature = "rayon")]
115    fn liquid_density_parallel<const P: usize>(
116        parameter_names: [String; P],
117        parameters: ArrayView2<f64>,
118        input: ArrayView2<f64>,
119    ) -> (Array1<f64>, Array2<f64>, Array1<bool>)
120    where
121        Self: ParametersAD<1>,
122    {
123        parallelize::<_, Self, _, _>(
124            parameter_names,
125            parameters,
126            input,
127            |eos: &Self::Lifted<Gradient<P>>, inp| {
128                eos.liquid_density(inp[0] * KELVIN, inp[1] * PASCAL)
129                    .map(|d| d.convert_into(KILO * MOL / (METER * METER * METER)))
130            },
131        )
132    }
133
134    #[cfg(feature = "rayon")]
135    fn equilibrium_liquid_density_parallel<const P: usize>(
136        parameter_names: [String; P],
137        parameters: ArrayView2<f64>,
138        input: ArrayView2<f64>,
139    ) -> (Array1<f64>, Array2<f64>, Array1<bool>)
140    where
141        Self: ParametersAD<1>,
142    {
143        parallelize::<_, Self, _, _>(
144            parameter_names,
145            parameters,
146            input,
147            |eos: &Self::Lifted<Gradient<P>>, inp| {
148                eos.equilibrium_liquid_density(inp[0] * KELVIN)
149                    .map(|(_, d)| d.convert_into(KILO * MOL / (METER * METER * METER)))
150            },
151        )
152    }
153
154    fn bubble_point_pressure<const P: usize>(
155        &self,
156        temperature: Temperature,
157        pressure: Option<Pressure>,
158        liquid_molefracs: SVector<f64, 2>,
159    ) -> FeosResult<Pressure<Gradient<P>>>
160    where
161        Self: Residual<U2, Gradient<P>>,
162    {
163        let eos_f64 = self.re();
164        let vle = PhaseEquilibrium::bubble_point(
165            &eos_f64,
166            temperature,
167            &liquid_molefracs,
168            pressure,
169            None,
170            Default::default(),
171        )?;
172
173        // implicit differentiation is implemented here instead of just calling bubble_point with dual
174        // numbers, because for the first derivative, we can avoid calculating density derivatives.
175        let v_l = 1.0 / vle.liquid().density.to_reduced();
176        let v_v = 1.0 / vle.vapor().density.to_reduced();
177        let y = &vle.vapor().molefracs;
178        let y: SVector<_, 2> = SVector::from_fn(|i, _| y[i]);
179        let t = temperature.into_reduced();
180        let (a_l, a_v, v_l, v_v) = {
181            let t = Gradient::from(t);
182            let v_l = Gradient::from(v_l);
183            let v_v = Gradient::from(v_v);
184            let y = y.map(Gradient::from);
185            let x = liquid_molefracs.map(Gradient::from);
186
187            let a_v = self.residual_molar_helmholtz_energy(t, v_v, &y);
188            let (p_l, mu_res_l, dp_l, dmu_l) = self.dmu_dv(t, v_l, &x);
189            let vi_l = dmu_l / dp_l;
190            let v_l = vi_l.dot(&y);
191            let a_l = (mu_res_l - vi_l * p_l).dot(&y);
192            (a_l, a_v, v_l, v_v)
193        };
194        let rho_l = vle.liquid().partial_density.to_reduced();
195        let rho_l = [rho_l[0], rho_l[1]];
196        let rho_v = vle.vapor().partial_density.to_reduced();
197        let rho_v = [rho_v[0], rho_v[1]];
198        let p = -(a_v - a_l
199            + t * (y[0] * (rho_v[0] / rho_l[0]).ln() + y[1] * (rho_v[1] / rho_l[1]).ln() - 1.0))
200            / (v_v - v_l);
201        Ok(Pressure::from_reduced(p))
202    }
203
204    fn dew_point_pressure<const P: usize>(
205        &self,
206        temperature: Temperature,
207        pressure: Option<Pressure>,
208        vapor_molefracs: SVector<f64, 2>,
209    ) -> FeosResult<Pressure<Gradient<P>>>
210    where
211        Self: Residual<U2, Gradient<P>>,
212    {
213        let eos_f64 = self.re();
214        let vle = PhaseEquilibrium::dew_point(
215            &eos_f64,
216            temperature,
217            &vapor_molefracs,
218            pressure,
219            None,
220            Default::default(),
221        )?;
222
223        // implicit differentiation is implemented here instead of just calling dew_point with dual
224        // numbers, because for the first derivative, we can avoid calculating density derivatives.
225        let v_l = 1.0 / vle.liquid().density.to_reduced();
226        let v_v = 1.0 / vle.vapor().density.to_reduced();
227        let x = &vle.liquid().molefracs;
228        let x: SVector<_, 2> = SVector::from_fn(|i, _| x[i]);
229        let t = temperature.into_reduced();
230        let (a_l, a_v, v_l, v_v) = {
231            let t = Gradient::from(t);
232            let v_l = Gradient::from(v_l);
233            let v_v = Gradient::from(v_v);
234            let x = x.map(Gradient::from);
235            let y = vapor_molefracs.map(Gradient::from);
236
237            let a_l = self.residual_molar_helmholtz_energy(t, v_l, &x);
238            let (p_v, mu_res_v, dp_v, dmu_v) = self.dmu_dv(t, v_v, &y);
239            let vi_v = dmu_v / dp_v;
240            let v_v = vi_v.dot(&x);
241            let a_v = (mu_res_v - vi_v * p_v).dot(&x);
242            (a_l, a_v, v_l, v_v)
243        };
244        let rho_l = vle.liquid().partial_density.to_reduced();
245        let rho_l = [rho_l[0], rho_l[1]];
246        let rho_v = vle.vapor().partial_density.to_reduced();
247        let rho_v = [rho_v[0], rho_v[1]];
248        let p = -(a_l - a_v
249            + t * (x[0] * (rho_l[0] / rho_v[0]).ln() + x[1] * (rho_l[1] / rho_v[1]).ln() - 1.0))
250            / (v_l - v_v);
251        Ok(Pressure::from_reduced(p))
252    }
253
254    #[cfg(feature = "rayon")]
255    fn bubble_point_pressure_parallel<const P: usize>(
256        parameter_names: [String; P],
257        parameters: ArrayView2<f64>,
258        input: ArrayView2<f64>,
259    ) -> (Array1<f64>, Array2<f64>, Array1<bool>)
260    where
261        Self: ParametersAD<2>,
262    {
263        parallelize::<_, Self, _, _>(
264            parameter_names,
265            parameters,
266            input,
267            |eos: &Self::Lifted<Gradient<P>>, inp| {
268                eos.bubble_point_pressure(
269                    inp[0] * KELVIN,
270                    Some(inp[2] * PASCAL),
271                    SVector::from([inp[1], 1.0 - inp[1]]),
272                )
273                .map(|p| p.convert_into(PASCAL))
274            },
275        )
276    }
277
278    #[cfg(feature = "rayon")]
279    fn dew_point_pressure_parallel<const P: usize>(
280        parameter_names: [String; P],
281        parameters: ArrayView2<f64>,
282        input: ArrayView2<f64>,
283    ) -> (Array1<f64>, Array2<f64>, Array1<bool>)
284    where
285        Self: ParametersAD<2>,
286    {
287        parallelize::<_, Self, _, _>(
288            parameter_names,
289            parameters,
290            input,
291            |eos: &Self::Lifted<Gradient<P>>, inp| {
292                eos.dew_point_pressure(
293                    inp[0] * KELVIN,
294                    Some(inp[2] * PASCAL),
295                    SVector::from([inp[1], 1.0 - inp[1]]),
296                )
297                .map(|p| p.convert_into(PASCAL))
298            },
299        )
300    }
301}
302
303impl<T> PropertiesAD for T {}
304
305#[cfg(feature = "rayon")]
306fn parallelize<F, E: ParametersAD<N>, const N: usize, const P: usize>(
307    parameter_names: [String; P],
308    parameters: ArrayView2<f64>,
309    input: ArrayView2<f64>,
310    f: F,
311) -> (Array1<f64>, Array2<f64>, Array1<bool>)
312where
313    F: Fn(&E::Lifted<Gradient<P>>, &[f64]) -> FeosResult<Gradient<P>> + Sync,
314{
315    let parameter_names = parameter_names.each_ref().map(|s| s as &str);
316    let value_dual = Zip::from(parameters.rows())
317        .and(input.rows())
318        .par_map_collect(|par, inp| {
319            let par = par.as_slice().expect("Parameter array is not contiguous!");
320            let inp = inp.as_slice().expect("Input array is not contiguous!");
321            let eos = E::from(par).named_derivatives(parameter_names);
322            f(&eos, inp)
323        });
324    let status = value_dual.iter().map(|p| p.is_ok()).collect();
325    let value_dual: Array1<_> = value_dual.into_iter().flatten().collect();
326    let mut value = Array1::zeros(value_dual.len());
327    let mut grad = Array2::zeros([value_dual.len(), P]);
328    Zip::from(grad.rows_mut())
329        .and(&mut value)
330        .and(&value_dual)
331        .for_each(|mut grad, p, p_dual| {
332            *p = p_dual.re;
333            let eps = p_dual.eps.unwrap_generic(Const::<P>, U1).data.0[0].to_vec();
334            grad.assign(&Array1::from(eps));
335        });
336    (value, grad, status)
337}
feos_core/ad/mod.rs

feos_core/ad/
mod.rs
