Skip to main content

feos_core/ad/dataset/
pure.rs

1use std::{io, path::Path};
2
3use nalgebra::U1;
4use ndarray::{Array1, Array2, ArrayView1, ArrayView2};
5use serde::{Deserialize, Serialize};
6
7use crate::Residual;
8use crate::ad::Gradient;
9use crate::ad::properties::*;
10
11use super::{Dataset, DatasetAD, DatasetRecord, DatasetStorage, ParametersAD};
12
13#[derive(Deserialize, Serialize)]
14pub struct VaporPressureRecord {
15    pub temperature_k: f64,
16    pub vapor_pressure_pa: f64,
17}
18
19impl DatasetRecord for VaporPressureRecord {
20    const N_INPUTS: usize = 1;
21
22    fn input(&self, _column: usize) -> f64 {
23        self.temperature_k
24    }
25
26    fn target(&self) -> f64 {
27        self.vapor_pressure_pa
28    }
29}
30
31#[derive(Deserialize, Serialize)]
32pub struct LiquidDensityRecord {
33    pub temperature_k: f64,
34    pub pressure_pa: f64,
35    pub liquid_density_kmol_m3: f64,
36}
37
38impl DatasetRecord for LiquidDensityRecord {
39    const N_INPUTS: usize = 2;
40
41    fn input(&self, column: usize) -> f64 {
42        match column {
43            0 => self.temperature_k,
44            1 => self.pressure_pa,
45            _ => unreachable!("invalid liquid density input column"),
46        }
47    }
48
49    fn target(&self) -> f64 {
50        self.liquid_density_kmol_m3
51    }
52}
53
54#[derive(Deserialize, Serialize)]
55pub struct EquilibriumLiquidDensityRecord {
56    pub temperature_k: f64,
57    pub liquid_density_kmol_m3: f64,
58}
59
60impl DatasetRecord for EquilibriumLiquidDensityRecord {
61    const N_INPUTS: usize = 1;
62
63    fn input(&self, _column: usize) -> f64 {
64        self.temperature_k
65    }
66
67    fn target(&self) -> f64 {
68        self.liquid_density_kmol_m3
69    }
70}
71
72#[derive(Deserialize, Serialize)]
73pub struct EnthalpyOfVaporizationRecord {
74    pub temperature_k: f64,
75    pub dh_vap_j_mol: f64,
76}
77
78impl DatasetRecord for EnthalpyOfVaporizationRecord {
79    const N_INPUTS: usize = 1;
80
81    fn input(&self, _column: usize) -> f64 {
82        self.temperature_k
83    }
84
85    fn target(&self) -> f64 {
86        self.dh_vap_j_mol
87    }
88}
89
90#[derive(Deserialize, Serialize)]
91pub struct ResidualIsobaricHeatCapacityRecord {
92    pub temperature_k: f64,
93    pub pressure_pa: f64,
94    pub cp_res_j_molk: f64,
95}
96
97impl DatasetRecord for ResidualIsobaricHeatCapacityRecord {
98    const N_INPUTS: usize = 2;
99
100    fn input(&self, column: usize) -> f64 {
101        match column {
102            0 => self.temperature_k,
103            1 => self.pressure_pa,
104            _ => unreachable!("invalid residual isobaric heat capacity input column"),
105        }
106    }
107
108    fn target(&self) -> f64 {
109        self.cp_res_j_molk
110    }
111}
112
113/// Expand a list of pure-component property entries into:
114/// - the [`PureProperty`] enum, metadata and dispatch methods,
115/// - constructors,
116/// - [`PureDataset::from_csv`] and [`PureDataset::from_reader`].
117///
118/// Adding a new property:
119/// - write the property file (record, `*_ad`, `*_parallel`, `*_parallel_ad`)
120/// - add entry here.
121macro_rules! pure_properties {
122    ($(
123        $variant:ident {
124            record:       $record:ty,
125            property:     $prop:ty,
126            default_name: $default:expr,
127            input_names:  $inputs:expr,
128            target_name:  $target:expr,
129            constructor:  $ctor:ident,
130        }
131    ),* $(,)?) => {
132        /// Pure-component properties.
133        #[derive(Debug, Clone, Copy, PartialEq, Eq)]
134        pub enum PureProperty {
135            $($variant,)*
136        }
137
138        impl PureProperty {
139            pub fn default_name(self) -> &'static str {
140                match self { $(Self::$variant => $default,)* }
141            }
142
143            pub fn input_names(self) -> &'static [&'static str] {
144                match self { $(Self::$variant => $inputs,)* }
145            }
146
147            pub fn target_name(self) -> &'static str {
148                match self { $(Self::$variant => $target,)* }
149            }
150
151            fn evaluate_ad<T: ParametersAD<U1>, const P: usize>(
152                self,
153                names: [String; P],
154                parameters: &[f64],
155                inputs: ArrayView2<f64>,
156            ) -> (Array1<f64>, Array2<f64>, Array1<bool>)
157            where T::Lifted<Gradient<P>>: Sync
158            {
159                match self {
160                    $(Self::$variant => <$prop>::evaluate_parallel_derivatives::<T, P>(names, parameters, inputs),)*
161                }
162            }
163
164            fn evaluate<E: Residual + Sync>(self, eos: &E, inputs: ArrayView2<f64>) -> (Array1<f64>, Array1<bool>)
165            {
166                match self {
167                    $(Self::$variant => <$prop>::evaluate_parallel(eos, inputs.view()),)*
168                }
169            }
170        }
171
172        impl PureDataset {
173            $(
174                pub fn $ctor(records: Vec<$record>) -> Self {
175                    Self {
176                        property: PureProperty::$variant,
177                        storage: DatasetStorage::from_records(records),
178                    }
179                }
180            )*
181
182            pub fn from_csv(property: PureProperty, path: &Path) -> Result<Self, csv::Error> {
183                let storage = match property {
184                    $(PureProperty::$variant => DatasetStorage::from_csv::<$record>(path)?,)*
185                };
186                Ok(Self { property, storage })
187            }
188
189            pub fn from_reader(
190                property: PureProperty,
191                reader: impl io::Read,
192            ) -> Result<Self, csv::Error> {
193                let storage = match property {
194                    $(PureProperty::$variant => DatasetStorage::from_reader::<$record>(reader)?,)*
195                };
196                Ok(Self { property, storage })
197            }
198        }
199    };
200}
201
202pure_properties! {
203    VaporPressure {
204        record:       VaporPressureRecord,
205        property:     VaporPressure,
206        default_name: "vapor pressure",
207        input_names:  &["temperature_k"],
208        target_name:  "vapor_pressure_pa",
209        constructor:  vapor_pressure,
210    },
211    LiquidDensity {
212        record:       LiquidDensityRecord,
213        property:     LiquidDensity,
214        default_name: "liquid density",
215        input_names:  &["temperature_k", "pressure_pa"],
216        target_name:  "liquid_density_kmol_m3",
217        constructor:  liquid_density,
218    },
219    EquilibriumLiquidDensity {
220        record:       EquilibriumLiquidDensityRecord,
221        property:     EquilibriumLiquidDensity,
222        default_name: "equilibrium liquid density",
223        input_names:  &["temperature_k"],
224        target_name:  "liquid_density_kmol_m3",
225        constructor:  equilibrium_liquid_density,
226    },
227    EnthalpyOfVaporization {
228        record:       EnthalpyOfVaporizationRecord,
229        property:     EnthalpyOfVaporization,
230        default_name: "enthalpy of vaporization",
231        input_names:  &["temperature_k"],
232        target_name:  "dh_vap_j_mol",
233        constructor:  enthalpy_of_vaporization,
234    },
235    ResidualIsobaricHeatCapacity {
236        record:       ResidualIsobaricHeatCapacityRecord,
237        property:     ResidualIsobaricHeatCapacity,
238        default_name: "residual isobaric heat capacity",
239        input_names:  &["temperature_k", "pressure_pa"],
240        target_name:  "cp_res_j_molk",
241        constructor:  residual_isobaric_heat_capacity,
242    },
243}
244
245/// Pure-component dataset: shared data storage plus a property tag.
246// #[derive(Clone)]
247pub struct PureDataset {
248    property: PureProperty,
249    storage: DatasetStorage,
250}
251
252impl PureDataset {
253    pub fn with_name(mut self, name: impl Into<String>) -> Self {
254        self.storage.set_name(name.into());
255        self
256    }
257
258    pub fn property(&self) -> PureProperty {
259        self.property
260    }
261
262    pub fn inputs(&self) -> ArrayView2<'_, f64> {
263        self.storage.inputs()
264    }
265
266    pub fn target(&self) -> ArrayView1<'_, f64> {
267        self.storage.target()
268    }
269
270    pub fn name(&self) -> &str {
271        self.storage.name().unwrap_or(self.property.default_name())
272    }
273
274    pub fn input_names(&self) -> &'static [&'static str] {
275        self.property.input_names()
276    }
277
278    pub fn target_name(&self) -> &'static str {
279        self.property.target_name()
280    }
281}
282
283impl Dataset for PureDataset {
284    fn inputs(&self) -> ArrayView2<'_, f64> {
285        self.inputs()
286    }
287
288    fn target(&self) -> ArrayView1<'_, f64> {
289        self.target()
290    }
291
292    fn name(&self) -> &str {
293        self.name()
294    }
295
296    fn input_names(&self) -> &'static [&'static str] {
297        self.input_names()
298    }
299
300    fn target_name(&self) -> &'static str {
301        self.target_name()
302    }
303
304    fn evaluate<E: Residual + Sync>(&self, eos: &E) -> (Array1<f64>, Array1<bool>) {
305        self.property.evaluate(eos, self.inputs().view())
306    }
307}
308
309impl DatasetAD<1> for PureDataset {
310    fn evaluate_ad_const<T: ParametersAD<U1>, const P: usize>(
311        &self,
312        names: [String; P],
313        parameters: &[f64],
314        inputs: ArrayView2<f64>,
315    ) -> (Array1<f64>, Array2<f64>, Array1<bool>)
316    where
317        T::Lifted<Gradient<P>>: Sync,
318    {
319        self.property.evaluate_ad::<T, P>(names, parameters, inputs)
320    }
321}
322
323#[cfg(test)]
324mod tests {
325    use super::*;
326    use std::io::Cursor;
327
328    fn csv(s: &str) -> Cursor<&[u8]> {
329        Cursor::new(s.as_bytes())
330    }
331
332    #[test]
333    fn vapor_pressure_from_reader() {
334        let data = "\
335temperature_k,vapor_pressure_pa
336300.0,3540.0
337350.0,41682.0
338400.0,245600.0
339";
340        let ds = PureDataset::from_reader(PureProperty::VaporPressure, csv(data)).unwrap();
341
342        assert_eq!(ds.target().len(), 3);
343        assert_eq!(ds.inputs().nrows(), 3);
344        assert_eq!(ds.inputs().ncols(), 1);
345        assert_eq!(ds.inputs()[[0, 0]], 300.0);
346        assert_eq!(ds.inputs()[[1, 0]], 350.0);
347        assert_eq!(ds.inputs()[[2, 0]], 400.0);
348        assert_eq!(ds.target()[0], 3540.0);
349        assert_eq!(ds.target()[1], 41682.0);
350        assert_eq!(ds.target()[2], 245600.0);
351        assert_eq!(ds.name(), "vapor pressure");
352    }
353
354    #[test]
355    fn vapor_pressure_from_records() {
356        let records = vec![
357            VaporPressureRecord {
358                temperature_k: 300.0,
359                vapor_pressure_pa: 3540.0,
360            },
361            VaporPressureRecord {
362                temperature_k: 350.0,
363                vapor_pressure_pa: 41682.0,
364            },
365        ];
366        let ds = PureDataset::vapor_pressure(records);
367        assert_eq!(ds.inputs()[[0, 0]], 300.0);
368        assert_eq!(ds.target()[1], 41682.0);
369    }
370
371    #[test]
372    fn liquid_density_from_reader() {
373        let data = "\
374temperature_k,pressure_pa,liquid_density_kmol_m3
375300.0,101325.0,15.2
376320.0,200000.0,14.8
377";
378        let ds = PureDataset::from_reader(PureProperty::LiquidDensity, csv(data)).unwrap();
379
380        assert_eq!(ds.inputs().nrows(), 2);
381        assert_eq!(ds.inputs().ncols(), 2);
382        assert_eq!(ds.inputs()[[0, 0]], 300.0);
383        assert_eq!(ds.inputs()[[0, 1]], 101325.0);
384        assert_eq!(ds.inputs()[[1, 0]], 320.0);
385        assert_eq!(ds.inputs()[[1, 1]], 200000.0);
386        assert_eq!(ds.target()[0], 15.2);
387        assert_eq!(ds.target()[1], 14.8);
388        assert_eq!(ds.name(), "liquid density");
389    }
390
391    #[test]
392    fn liquid_density_from_records() {
393        let records = vec![LiquidDensityRecord {
394            temperature_k: 300.0,
395            pressure_pa: 101325.0,
396            liquid_density_kmol_m3: 15.2,
397        }];
398        let ds = PureDataset::liquid_density(records);
399        assert_eq!(ds.inputs()[[0, 1]], 101325.0);
400        assert_eq!(ds.target()[0], 15.2);
401    }
402
403    #[test]
404    fn equilibrium_liquid_density_from_reader() {
405        let data = "\
406temperature_k,liquid_density_kmol_m3
407290.0,15.5
408310.0,14.9
409330.0,14.1
410";
411        let ds =
412            PureDataset::from_reader(PureProperty::EquilibriumLiquidDensity, csv(data)).unwrap();
413
414        assert_eq!(ds.inputs().ncols(), 1);
415        assert_eq!(ds.inputs().nrows(), 3);
416        assert_eq!(ds.inputs()[[2, 0]], 330.0);
417        assert_eq!(ds.target()[2], 14.1);
418        assert_eq!(ds.name(), "equilibrium liquid density");
419    }
420
421    #[test]
422    fn missing_column_returns_error() {
423        let data = "temperature_k\n300.0\n";
424        let result = PureDataset::from_reader(PureProperty::VaporPressure, csv(data));
425        assert!(result.is_err());
426    }
427
428    #[test]
429    fn wrong_type_returns_error() {
430        let data = "temperature_k,vapor_pressure_pa\n300.0,not_a_number\n";
431        let result = PureDataset::from_reader(PureProperty::VaporPressure, csv(data));
432        assert!(result.is_err());
433    }
434}