Skip to main content

feos_core/ad/properties/
mod.rs

1use super::Gradient;
2use crate::{FeosResult, Residual};
3use nalgebra::{DefaultAllocator, Dim, allocator::Allocator};
4#[cfg(feature = "ndarray")]
5use ndarray::{Array1, Array2, ArrayView2};
6use num_dual::DualNum;
7use quantity::Quantity;
8
9mod boiling_temperature;
10mod bubble_point_pressure;
11mod dew_point_pressure;
12mod enthalpy_of_vaporization;
13mod equilibrium_liquid_density;
14mod liquid_density;
15mod residual_isobaric_heat_capacity;
16mod vapor_pressure;
17
18pub use boiling_temperature::BoilingTemperature;
19pub use bubble_point_pressure::BubblePointPressure;
20pub use dew_point_pressure::DewPointPressure;
21pub use enthalpy_of_vaporization::EnthalpyOfVaporization;
22pub use equilibrium_liquid_density::EquilibriumLiquidDensity;
23pub use liquid_density::LiquidDensity;
24pub use residual_isobaric_heat_capacity::ResidualIsobaricHeatCapacity;
25pub use vapor_pressure::VaporPressure;
26
27/// Properties that can be rapidly evaluated in parallel together with
28/// their gradients with respect to model parameters
29pub trait PropertyAD<N: Dim>: for<'a> From<&'a [f64]>
30where
31    DefaultAllocator: Allocator<N>,
32{
33    type Unit;
34    const REFERENCE: Quantity<f64, Self::Unit>;
35
36    /// Evaluate the property for an arbitrary derivative.
37    fn evaluate<E: Residual<N, D>, D: DualNum<f64, Inner = f64> + Copy>(
38        &self,
39        eos: &E,
40    ) -> FeosResult<Quantity<D, Self::Unit>>;
41
42    /// Evaluate the property for the first derivative w.r.t. model parameters.
43    ///
44    /// This can be overridden if there is a more performant implementation than
45    /// the general implementation in `evaluate`.
46    fn evaluate_gradient<E: Residual<N, Gradient<P>>, const P: usize>(
47        &self,
48        eos: &E,
49    ) -> FeosResult<Quantity<Gradient<P>, Self::Unit>> {
50        self.evaluate(eos)
51    }
52
53    /// Evaluate the property for all inputs in parallel.
54    ///
55    /// Return the property values and the success of the calculations.
56    #[cfg(feature = "ndarray")]
57    fn evaluate_parallel<E: Residual<N> + Sync>(
58        eos: &E,
59        input: ArrayView2<f64>,
60    ) -> (Array1<f64>, Array1<bool>) {
61        #[cfg(feature = "rayon")]
62        let values = ndarray::Zip::from(input.rows()).par_map_collect(|inp| {
63            let inp = inp.as_slice().expect("Input array is not contiguous!");
64            Self::from(inp)
65                .evaluate(eos)
66                .map(|d| d.convert_into(Self::REFERENCE))
67        });
68
69        #[cfg(not(feature = "rayon"))]
70        let values = ndarray::Zip::from(input.rows()).map_collect(|inp| {
71            let inp = inp.as_slice().expect("Input array is not contiguous!");
72            Self::from(inp)
73                .evaluate(eos)
74                .map(|d| d.convert_into(Self::REFERENCE))
75        });
76
77        let n = input.nrows();
78        let status: Array1<bool> = values.iter().map(|r| r.is_ok()).collect();
79        let mut value = Array1::from_elem(n, f64::NAN);
80        for (i, result) in values.into_iter().enumerate() {
81            if let Ok(v) = result {
82                value[i] = v;
83            }
84        }
85        (value, status)
86    }
87
88    /// Evaluate the property and its gradients for all inputs in parallel.
89    ///
90    /// Return the property values, the gradients, and the success of the calculations.  
91    #[cfg(feature = "ndarray")]
92    fn evaluate_parallel_derivatives<E: super::ParametersAD<N>, const P: usize>(
93        parameter_names: [String; P],
94        parameters: &[f64],
95        input: ArrayView2<f64>,
96    ) -> (Array1<f64>, Array2<f64>, Array1<bool>)
97    where
98        E::Lifted<Gradient<P>>: Sync,
99    {
100        let parameter_names = parameter_names.each_ref().map(|s| s as &str);
101        let eos = E::seed_derivatives(parameters, parameter_names);
102
103        #[cfg(feature = "rayon")]
104        let value_dual = ndarray::Zip::from(input.rows()).par_map_collect(|inp| {
105            let inp = inp.as_slice().expect("Input array is not contiguous!");
106            Self::from(inp)
107                .evaluate_gradient(&eos)
108                .map(|d| d.convert_into(Self::REFERENCE))
109        });
110
111        #[cfg(not(feature = "rayon"))]
112        let value_dual = ndarray::Zip::from(input.rows()).map_collect(|inp| {
113            let inp = inp.as_slice().expect("Input array is not contiguous!");
114            Self::from(inp)
115                .evaluate_gradient(&eos)
116                .map(|d| d.convert_into(Self::REFERENCE))
117        });
118
119        let n = input.nrows();
120        let status = value_dual.iter().map(|p| p.is_ok()).collect();
121        let mut value = Array1::from_elem(n, f64::NAN);
122        let mut grad = Array2::zeros([n, P]);
123        for (i, result) in value_dual.into_iter().enumerate() {
124            if let Ok(p_dual) = result {
125                value[i] = p_dual.re;
126                let eps = p_dual
127                    .eps
128                    .unwrap_generic(nalgebra::Const::<P>, nalgebra::U1);
129                for (g, &e) in grad.row_mut(i).iter_mut().zip(eps.data.0[0].iter()) {
130                    *g = e;
131                }
132            }
133        }
134        (value, grad, status)
135    }
136
137    /// Evaluate the property and its gradients for all inputs and parameters in parallel.
138    ///
139    /// Return the property values, the gradients, and the success of the calculations.  
140    #[cfg(feature = "ndarray")]
141    fn evaluate_parallel_derivatives_params<E: super::ParametersAD<N>, const P: usize>(
142        parameter_names: [String; P],
143        parameters: ArrayView2<f64>,
144        input: ArrayView2<f64>,
145    ) -> (Array1<f64>, Array2<f64>, Array1<bool>) {
146        let parameter_names = parameter_names.each_ref().map(|s| s as &str);
147
148        #[cfg(feature = "rayon")]
149        let value_dual = ndarray::Zip::from(parameters.rows())
150            .and(input.rows())
151            .par_map_collect(|par, inp| {
152                let par = par.as_slice().expect("Parameter array is not contiguous!");
153                let inp = inp.as_slice().expect("Input array is not contiguous!");
154                let eos = E::seed_derivatives(par, parameter_names);
155                Self::from(inp)
156                    .evaluate_gradient(&eos)
157                    .map(|d| d.convert_into(Self::REFERENCE))
158            });
159
160        #[cfg(not(feature = "rayon"))]
161        let value_dual = ndarray::Zip::from(parameters.rows())
162            .and(input.rows())
163            .map_collect(|par, inp| {
164                let par = par.as_slice().expect("Parameter array is not contiguous!");
165                let inp = inp.as_slice().expect("Input array is not contiguous!");
166                let eos = E::seed_derivatives(par, parameter_names);
167                Self::from(inp)
168                    .evaluate_gradient(&eos)
169                    .map(|d| d.convert_into(Self::REFERENCE))
170            });
171
172        let n = parameters.nrows();
173        let status = value_dual.iter().map(|p| p.is_ok()).collect();
174        let mut value = Array1::from_elem(n, f64::NAN);
175        let mut grad = Array2::zeros([n, P]);
176        for (i, result) in value_dual.into_iter().enumerate() {
177            if let Ok(p_dual) = result {
178                value[i] = p_dual.re;
179                let eps = p_dual
180                    .eps
181                    .unwrap_generic(nalgebra::Const::<P>, nalgebra::U1);
182                for (g, &e) in grad.row_mut(i).iter_mut().zip(eps.data.0[0].iter()) {
183                    *g = e;
184                }
185            }
186        }
187        (value, grad, status)
188    }
189}