diff --git a/CHANGELOG.md b/CHANGELOG.md index 3d69100a7..85400f656 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,16 +12,23 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Added the `Composition` trait to allow more flexibility in the creation of states and phase equilibria. [#330](https://github.com/feos-org/feos/pull/330) - Added `PhaseEquilibrium::ph_flash` and `PhaseEquilibrium::ps_flash`. [#338](https://github.com/feos-org/feos/pull/338) - Added getters for `vapor_phase_fraction`, `molar_enthalpy`, `molar_entropy`, `total_moles`, `enthalpy`, and `entropy` to `PhaseEquilibrium`. [#338](https://github.com/feos-org/feos/pull/338) +- Added `PropertyAD` trait in `feos_core::ad` with one struct per property for uniform evaluation with or without parameter derivatives, including parallel variants. [#358](https://github.com/feos-org/feos/pull/358) +- Added `feos_core::ad::dataset` module with `PureDataset` and `BinaryDataset` types, constructible from records, CSV files, or readers, for use in parameter fits. [#358](https://github.com/feos-org/feos/pull/358) +- Exposed `Property`, `PureDataset`, and `BinaryDataset` in `py-feos`. [#358](https://github.com/feos-org/feos/pull/358) ### Changed - Removed any assumptions about the total number of moles in a `State` or `PhaseEquilibrium`. Evaluating extensive properties now returns a `Result`. [#330](https://github.com/feos-org/feos/pull/330) - Redesigned the `IdealGas` trait and added `IdealGasAD` in analogy to `ResidualDyn` and `Residual`. [#330](https://github.com/feos-org/feos/pull/330) +- Replaced the `PropertiesAD` blanket-impl trait with per-property `PropertyAD` types. [#358](https://github.com/feos-org/feos/pull/358) +- Replaced `ParametersAD::named_derivatives` with a `build` constructor and `seed_derivatives(&values, names)`. [#358](https://github.com/feos-org/feos/pull/358) +- Removed the per-property `*_derivatives` free functions in Python in favour of static methods on the new `Property` class. [#358](https://github.com/feos-org/feos/pull/358) ### Removed - Removed the `StateBuilder` struct, because it is mostly obsolete with the addition of the `Composition` trait. [#330](https://github.com/feos-org/feos/pull/330) ### Packaging - Updated `quantity` dependency to 0.13 and removed the `typenum` dependency. [#328](https://github.com/feos-org/feos/pull/328) +- Added `csv` as a `feos-core` dependency for the new dataset module. [#358](https://github.com/feos-org/feos/pull/358) ## [Unreleased] ### Added diff --git a/Cargo.toml b/Cargo.toml index de98c170f..eb351252c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -42,6 +42,7 @@ gauss-quad = "0.2" approx = "0.5" criterion = "0.8" paste = "1.0" +csv = "1.0" feos-core = { version = "0.9", path = "crates/feos-core" } feos-dft = { version = "0.9", path = "crates/feos-dft" } diff --git a/crates/feos-core/Cargo.toml b/crates/feos-core/Cargo.toml index f33cd1255..20e516cc3 100644 --- a/crates/feos-core/Cargo.toml +++ b/crates/feos-core/Cargo.toml @@ -25,6 +25,7 @@ serde = { workspace = true, features = ["derive"] } serde_json = { workspace = true, features = ["preserve_order"] } indexmap = { workspace = true, features = ["serde"] } rayon = { workspace = true, optional = true } +csv = { workspace = true } itertools = { workspace = true } [dev-dependencies] diff --git a/crates/feos-core/src/ad/dataset/binary.rs b/crates/feos-core/src/ad/dataset/binary.rs new file mode 100644 index 000000000..2bf797f6e --- /dev/null +++ b/crates/feos-core/src/ad/dataset/binary.rs @@ -0,0 +1,304 @@ +use std::{io, path::Path}; + +use nalgebra::U2; +use ndarray::{Array1, Array2, ArrayView1, ArrayView2}; +use serde::{Deserialize, Serialize}; + +use crate::Residual; +use crate::ad::Gradient; +use crate::ad::properties::{BubblePointPressure, DewPointPressure, PropertyAD}; + +use super::{Dataset, DatasetAD, DatasetRecord, DatasetStorage, ParametersAD}; + +/// The pressure column doubles as the initial guess passed to the VLE solver. +#[derive(Deserialize, Serialize)] +pub struct BubblePointRecord { + pub temperature_k: f64, + pub liquid_molefrac_1: f64, + pub bubble_pressure_pa: f64, +} + +impl DatasetRecord for BubblePointRecord { + const N_INPUTS: usize = 3; + + fn input(&self, column: usize) -> f64 { + match column { + 0 => self.temperature_k, + 1 => self.liquid_molefrac_1, + 2 => self.bubble_pressure_pa, + _ => unreachable!("invalid bubble point input column"), + } + } + + fn target(&self) -> f64 { + self.bubble_pressure_pa + } +} + +/// The pressure column doubles as the initial guess passed to the VLE solver. +#[derive(Deserialize, Serialize)] +pub struct DewPointRecord { + pub temperature_k: f64, + pub vapor_molefrac_1: f64, + pub dew_pressure_pa: f64, +} + +impl DatasetRecord for DewPointRecord { + const N_INPUTS: usize = 3; + + fn input(&self, column: usize) -> f64 { + match column { + 0 => self.temperature_k, + 1 => self.vapor_molefrac_1, + 2 => self.dew_pressure_pa, + _ => unreachable!("invalid dew point input column"), + } + } + + fn target(&self) -> f64 { + self.dew_pressure_pa + } +} + +/// Expand a list of binary-mixture property entries into: +/// - the [`BinaryProperty`] enum, metadata and dispatch methods, +/// - constructors, +/// - [`BinaryDataset::from_csv`] and [`BinaryDataset::from_reader`] match arms. +macro_rules! binary_properties { + ($( + $variant:ident { + record: $record:ty, + property: $prop:ty, + default_name: $default:expr, + input_names: $inputs:expr, + target_name: $target:expr, + constructor: $ctor:ident, + } + ),* $(,)?) => { + /// Binary-mixture properties. + #[derive(Debug, Clone, Copy, PartialEq, Eq)] + pub enum BinaryProperty { + $($variant,)* + } + + impl BinaryProperty { + pub fn default_name(self) -> &'static str { + match self { $(Self::$variant => $default,)* } + } + + pub fn input_names(self) -> &'static [&'static str] { + match self { $(Self::$variant => $inputs,)* } + } + + pub fn target_name(self) -> &'static str { + match self { $(Self::$variant => $target,)* } + } + + fn evaluate_ad, const P: usize>( + self, + names: [String; P], + parameters: &[f64], + inputs: ArrayView2, + ) -> (Array1, Array2, Array1) + where + T::Lifted>: Sync, + { + match self { + $(Self::$variant => <$prop>::evaluate_parallel_derivatives::(names, parameters, inputs),)* + } + } + + fn evaluate(self, eos: &E, inputs: ArrayView2) -> (Array1, Array1) + where + E: Residual + Sync, + { + match self { + $(Self::$variant => <$prop>::evaluate_parallel(eos, inputs),)* + } + } + } + + impl BinaryDataset { + $( + pub fn $ctor(records: Vec<$record>) -> Self { + Self { + property: BinaryProperty::$variant, + storage: DatasetStorage::from_records(records), + } + } + )* + + pub fn from_csv(property: BinaryProperty, path: &Path) -> Result { + let storage = match property { + $(BinaryProperty::$variant => DatasetStorage::from_csv::<$record>(path)?,)* + }; + Ok(Self { property, storage }) + } + + pub fn from_reader( + property: BinaryProperty, + reader: impl io::Read, + ) -> Result { + let storage = match property { + $(BinaryProperty::$variant => DatasetStorage::from_reader::<$record>(reader)?,)* + }; + Ok(Self { property, storage }) + } + } + }; +} + +binary_properties! { + BubblePointPressure { + record: BubblePointRecord, + property: BubblePointPressure, + default_name: "bubble point pressure", + input_names: &["temperature_k", "liquid_molefrac_1"], + target_name: "bubble_pressure_pa", + constructor: bubble_point_pressure, + }, + DewPointPressure { + record: DewPointRecord, + property: DewPointPressure, + default_name: "dew point pressure", + input_names: &["temperature_k", "vapor_molefrac_1"], + target_name: "dew_pressure_pa", + constructor: dew_point_pressure, + }, +} + +/// Binary-mixture dataset: shared data storage plus a property tag. +#[derive(Clone)] +pub struct BinaryDataset { + property: BinaryProperty, + storage: DatasetStorage, +} + +impl BinaryDataset { + pub fn with_name(mut self, name: impl Into) -> Self { + self.storage.set_name(name.into()); + self + } + + pub fn property(&self) -> BinaryProperty { + self.property + } + + pub fn inputs(&self) -> ArrayView2<'_, f64> { + self.storage.inputs() + } + + pub fn target(&self) -> ArrayView1<'_, f64> { + self.storage.target() + } + + pub fn name(&self) -> &str { + self.storage.name().unwrap_or(self.property.default_name()) + } + + pub fn input_names(&self) -> &'static [&'static str] { + self.property.input_names() + } + + pub fn target_name(&self) -> &'static str { + self.property.target_name() + } +} + +impl Dataset for BinaryDataset { + fn inputs(&self) -> ArrayView2<'_, f64> { + self.inputs() + } + + fn target(&self) -> ArrayView1<'_, f64> { + self.target() + } + + fn name(&self) -> &str { + self.name() + } + + fn input_names(&self) -> &'static [&'static str] { + self.input_names() + } + + fn target_name(&self) -> &'static str { + self.target_name() + } + + fn evaluate(&self, eos: &E) -> (Array1, Array1) { + self.property.evaluate(eos, self.inputs()) + } +} + +impl DatasetAD<2> for BinaryDataset { + fn evaluate_ad_const, const P: usize>( + &self, + names: [String; P], + parameters: &[f64], + inputs: ArrayView2, + ) -> (Array1, Array2, Array1) + where + T::Lifted>: Sync, + { + self.property.evaluate_ad::(names, parameters, inputs) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use std::io::Cursor; + + fn csv(s: &str) -> Cursor<&[u8]> { + Cursor::new(s.as_bytes()) + } + + #[test] + fn bubble_point_from_reader() { + let data = "\ +temperature_k,liquid_molefrac_1,bubble_pressure_pa +300.0,0.3,500000.0 +320.0,0.5,800000.0 +"; + let ds = + BinaryDataset::from_reader(BinaryProperty::BubblePointPressure, csv(data)).unwrap(); + + assert_eq!(ds.inputs().ncols(), 3); + assert_eq!(ds.inputs().nrows(), 2); + assert_eq!(ds.inputs()[[0, 0]], 300.0); + assert_eq!(ds.inputs()[[0, 1]], 0.3); + assert_eq!(ds.inputs()[[0, 2]], 500000.0); + assert_eq!(ds.target()[0], 500000.0); + assert_eq!(ds.target()[1], 800000.0); + assert_eq!(ds.name(), "bubble point pressure"); + } + + #[test] + fn dew_point_from_reader() { + let data = "\ +temperature_k,vapor_molefrac_1,dew_pressure_pa +310.0,0.7,400000.0 +330.0,0.9,700000.0 +"; + let ds = BinaryDataset::from_reader(BinaryProperty::DewPointPressure, csv(data)).unwrap(); + + assert_eq!(ds.inputs().ncols(), 3); + assert_eq!(ds.inputs()[[0, 0]], 310.0); + assert_eq!(ds.inputs()[[0, 1]], 0.7); + assert_eq!(ds.inputs()[[0, 2]], 400000.0); + assert_eq!(ds.target()[1], 700000.0); + assert_eq!(ds.name(), "dew point pressure"); + } + + #[test] + fn dew_point_pressure_doubles_as_initial_guess() { + let records = vec![DewPointRecord { + temperature_k: 310.0, + vapor_molefrac_1: 0.7, + dew_pressure_pa: 400000.0, + }]; + let ds = BinaryDataset::dew_point_pressure(records); + assert_eq!(ds.inputs()[[0, 2]], ds.target()[0]); + } +} diff --git a/crates/feos-core/src/ad/dataset/mod.rs b/crates/feos-core/src/ad/dataset/mod.rs new file mode 100644 index 000000000..7be48a6fb --- /dev/null +++ b/crates/feos-core/src/ad/dataset/mod.rs @@ -0,0 +1,176 @@ +mod binary; +mod pure; + +use std::{io, path::Path, sync::Arc}; + +use nalgebra::Const; +use ndarray::{Array1, Array2, ArrayView1, ArrayView2}; +use serde::de::DeserializeOwned; + +use crate::Residual; +use crate::ad::Gradient; + +use super::ParametersAD; + +pub use binary::*; +pub use pure::*; + +/// Shared numerical data for all datasets. +struct DatasetData { + inputs: Array2, + target: Array1, +} + +/// Shared representation for all datasets. +#[derive(Clone)] +struct DatasetStorage { + data: Arc, + name: Option, +} + +impl DatasetStorage { + fn from_records(records: Vec) -> Self { + let n = records.len(); + let inputs = Array2::from_shape_fn((n, R::N_INPUTS), |(i, j)| records[i].input(j)); + let target = Array1::from_iter(records.iter().map(DatasetRecord::target)); + Self { + data: Arc::new(DatasetData { inputs, target }), + name: None, + } + } + + fn from_csv(path: &Path) -> Result { + let records = csv::Reader::from_path(path)? + .deserialize() + .collect::, _>>()?; + Ok(Self::from_records(records)) + } + + fn from_reader(reader: impl io::Read) -> Result { + let records = csv::Reader::from_reader(reader) + .deserialize() + .collect::, _>>()?; + Ok(Self::from_records(records)) + } + + fn inputs(&self) -> ArrayView2<'_, f64> { + self.data.inputs.view() + } + + fn target(&self) -> ArrayView1<'_, f64> { + self.data.target.view() + } + + fn name(&self) -> Option<&str> { + self.name.as_deref() + } + + fn set_name(&mut self, name: String) { + self.name = Some(name); + } +} + +/// A record that can be collected into a dataset. +pub trait DatasetRecord: DeserializeOwned { + /// Number of columns for inputs. + const N_INPUTS: usize; + + /// Value of EoS input column. + fn input(&self, column: usize) -> f64; + + /// Target value. + fn target(&self) -> f64; +} + +/// Dataset that can be evaluated by an equation of state. +pub trait Dataset { + /// Inputs for EoS evaluation, shape `[n_points, k]`. + fn inputs(&self) -> ArrayView2<'_, f64>; + + /// Target values, shape `[n_points]`. + fn target(&self) -> ArrayView1<'_, f64>; + + /// Property name. + /// + /// Used for logging and diagnostics. + fn name(&self) -> &str; + + /// Names of independent input columns. + fn input_names(&self) -> &'static [&'static str]; + + /// Name of the target property. + fn target_name(&self) -> &'static str; + + /// Evaluate this dataset's property with an equation of state. + /// + /// Returns `(predicted, converged)`: + /// - `predicted`: shape `[n_points]`, in SI units; `NaN` where the + /// underlying solver did not converge. + /// - `converged`: shape `[n_points]`. + fn evaluate(&self, eos: &E) -> (Array1, Array1); +} + +/// Build [`GRADIENT_SLOTS`] and [`DatasetAD`] trait from a single list of 'slots'. +/// +/// For each slot in [`GRADIENT_SLOTS`] a compile-time constant `P` variant is generated. +macro_rules! define_dataset_ad { + ($($p:literal),+ $(,)?) => { + /// Compile-time gradient slot supported by [`DatasetAD::evaluate_ad`]. + pub const GRADIENT_SLOTS: &[usize] = &[$($p),+]; + + /// Dataset that supports parameter-gradient evaluation + /// for equations of state implementing [`ParametersAD`]. + pub trait DatasetAD: Dataset { + /// Evaluate the property and its `P` parameter gradients. + fn evaluate_ad_const>, const P: usize>( + &self, + names: [String; P], + parameters: &[f64], + inputs: ArrayView2, + ) -> (Array1, Array2, Array1) + where + T::Lifted>: Sync; + + /// Evaluate the property and its parameter gradients at the given parameters. + /// + /// - `param_names`: names of the `P` parameters being differentiated. + /// - `params`: the full parameter vector. Only entries listed in `param_names` are seeded. + /// + /// This function dispatches the const-P methods at run-time. + fn evaluate_ad>>( + &self, + param_names: &[String], + parameters: &[f64], + ) -> (Array1, Array2, Array1) + where + $(T::Lifted>: Sync,)* + { + fn to_const(names: &[String]) -> [String; P] { + names.to_vec().try_into().expect("parameter count mismatch") + } + + match param_names.len() { + $( + $p => self.evaluate_ad_const::( + to_const(param_names), + parameters, + self.inputs().view(), + ), + )+ + p => unreachable!( + "parameter count {p} is not a member of GRADIENT_SLOTS={:?}", + GRADIENT_SLOTS, + ), + } + } + } + }; +} + +// We define the number of slots here. +// +// Note: might be good to investigate whether a smaller list makes sense here. +// LLVM vectorises across entries in DualSVec. We might see no perf. difference +// when using e.g. 3 vs 4 slots (even if only 3 are needed by the user). +// If the number of monomophised variants ever gets problematic, we could reduce it that way. +define_dataset_ad!(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14); diff --git a/crates/feos-core/src/ad/dataset/pure.rs b/crates/feos-core/src/ad/dataset/pure.rs new file mode 100644 index 000000000..4a7e74aa8 --- /dev/null +++ b/crates/feos-core/src/ad/dataset/pure.rs @@ -0,0 +1,434 @@ +use std::{io, path::Path}; + +use nalgebra::U1; +use ndarray::{Array1, Array2, ArrayView1, ArrayView2}; +use serde::{Deserialize, Serialize}; + +use crate::Residual; +use crate::ad::Gradient; +use crate::ad::properties::*; + +use super::{Dataset, DatasetAD, DatasetRecord, DatasetStorage, ParametersAD}; + +#[derive(Deserialize, Serialize)] +pub struct VaporPressureRecord { + pub temperature_k: f64, + pub vapor_pressure_pa: f64, +} + +impl DatasetRecord for VaporPressureRecord { + const N_INPUTS: usize = 1; + + fn input(&self, _column: usize) -> f64 { + self.temperature_k + } + + fn target(&self) -> f64 { + self.vapor_pressure_pa + } +} + +#[derive(Deserialize, Serialize)] +pub struct LiquidDensityRecord { + pub temperature_k: f64, + pub pressure_pa: f64, + pub liquid_density_kmol_m3: f64, +} + +impl DatasetRecord for LiquidDensityRecord { + const N_INPUTS: usize = 2; + + fn input(&self, column: usize) -> f64 { + match column { + 0 => self.temperature_k, + 1 => self.pressure_pa, + _ => unreachable!("invalid liquid density input column"), + } + } + + fn target(&self) -> f64 { + self.liquid_density_kmol_m3 + } +} + +#[derive(Deserialize, Serialize)] +pub struct EquilibriumLiquidDensityRecord { + pub temperature_k: f64, + pub liquid_density_kmol_m3: f64, +} + +impl DatasetRecord for EquilibriumLiquidDensityRecord { + const N_INPUTS: usize = 1; + + fn input(&self, _column: usize) -> f64 { + self.temperature_k + } + + fn target(&self) -> f64 { + self.liquid_density_kmol_m3 + } +} + +#[derive(Deserialize, Serialize)] +pub struct EnthalpyOfVaporizationRecord { + pub temperature_k: f64, + pub dh_vap_j_mol: f64, +} + +impl DatasetRecord for EnthalpyOfVaporizationRecord { + const N_INPUTS: usize = 1; + + fn input(&self, _column: usize) -> f64 { + self.temperature_k + } + + fn target(&self) -> f64 { + self.dh_vap_j_mol + } +} + +#[derive(Deserialize, Serialize)] +pub struct ResidualIsobaricHeatCapacityRecord { + pub temperature_k: f64, + pub pressure_pa: f64, + pub cp_res_j_molk: f64, +} + +impl DatasetRecord for ResidualIsobaricHeatCapacityRecord { + const N_INPUTS: usize = 2; + + fn input(&self, column: usize) -> f64 { + match column { + 0 => self.temperature_k, + 1 => self.pressure_pa, + _ => unreachable!("invalid residual isobaric heat capacity input column"), + } + } + + fn target(&self) -> f64 { + self.cp_res_j_molk + } +} + +/// Expand a list of pure-component property entries into: +/// - the [`PureProperty`] enum, metadata and dispatch methods, +/// - constructors, +/// - [`PureDataset::from_csv`] and [`PureDataset::from_reader`]. +/// +/// Adding a new property: +/// - write the property file (record, `*_ad`, `*_parallel`, `*_parallel_ad`) +/// - add entry here. +macro_rules! pure_properties { + ($( + $variant:ident { + record: $record:ty, + property: $prop:ty, + default_name: $default:expr, + input_names: $inputs:expr, + target_name: $target:expr, + constructor: $ctor:ident, + } + ),* $(,)?) => { + /// Pure-component properties. + #[derive(Debug, Clone, Copy, PartialEq, Eq)] + pub enum PureProperty { + $($variant,)* + } + + impl PureProperty { + pub fn default_name(self) -> &'static str { + match self { $(Self::$variant => $default,)* } + } + + pub fn input_names(self) -> &'static [&'static str] { + match self { $(Self::$variant => $inputs,)* } + } + + pub fn target_name(self) -> &'static str { + match self { $(Self::$variant => $target,)* } + } + + fn evaluate_ad, const P: usize>( + self, + names: [String; P], + parameters: &[f64], + inputs: ArrayView2, + ) -> (Array1, Array2, Array1) + where T::Lifted>: Sync + { + match self { + $(Self::$variant => <$prop>::evaluate_parallel_derivatives::(names, parameters, inputs),)* + } + } + + fn evaluate(self, eos: &E, inputs: ArrayView2) -> (Array1, Array1) + { + match self { + $(Self::$variant => <$prop>::evaluate_parallel(eos, inputs.view()),)* + } + } + } + + impl PureDataset { + $( + pub fn $ctor(records: Vec<$record>) -> Self { + Self { + property: PureProperty::$variant, + storage: DatasetStorage::from_records(records), + } + } + )* + + pub fn from_csv(property: PureProperty, path: &Path) -> Result { + let storage = match property { + $(PureProperty::$variant => DatasetStorage::from_csv::<$record>(path)?,)* + }; + Ok(Self { property, storage }) + } + + pub fn from_reader( + property: PureProperty, + reader: impl io::Read, + ) -> Result { + let storage = match property { + $(PureProperty::$variant => DatasetStorage::from_reader::<$record>(reader)?,)* + }; + Ok(Self { property, storage }) + } + } + }; +} + +pure_properties! { + VaporPressure { + record: VaporPressureRecord, + property: VaporPressure, + default_name: "vapor pressure", + input_names: &["temperature_k"], + target_name: "vapor_pressure_pa", + constructor: vapor_pressure, + }, + LiquidDensity { + record: LiquidDensityRecord, + property: LiquidDensity, + default_name: "liquid density", + input_names: &["temperature_k", "pressure_pa"], + target_name: "liquid_density_kmol_m3", + constructor: liquid_density, + }, + EquilibriumLiquidDensity { + record: EquilibriumLiquidDensityRecord, + property: EquilibriumLiquidDensity, + default_name: "equilibrium liquid density", + input_names: &["temperature_k"], + target_name: "liquid_density_kmol_m3", + constructor: equilibrium_liquid_density, + }, + EnthalpyOfVaporization { + record: EnthalpyOfVaporizationRecord, + property: EnthalpyOfVaporization, + default_name: "enthalpy of vaporization", + input_names: &["temperature_k"], + target_name: "dh_vap_j_mol", + constructor: enthalpy_of_vaporization, + }, + ResidualIsobaricHeatCapacity { + record: ResidualIsobaricHeatCapacityRecord, + property: ResidualIsobaricHeatCapacity, + default_name: "residual isobaric heat capacity", + input_names: &["temperature_k", "pressure_pa"], + target_name: "cp_res_j_molk", + constructor: residual_isobaric_heat_capacity, + }, +} + +/// Pure-component dataset: shared data storage plus a property tag. +// #[derive(Clone)] +pub struct PureDataset { + property: PureProperty, + storage: DatasetStorage, +} + +impl PureDataset { + pub fn with_name(mut self, name: impl Into) -> Self { + self.storage.set_name(name.into()); + self + } + + pub fn property(&self) -> PureProperty { + self.property + } + + pub fn inputs(&self) -> ArrayView2<'_, f64> { + self.storage.inputs() + } + + pub fn target(&self) -> ArrayView1<'_, f64> { + self.storage.target() + } + + pub fn name(&self) -> &str { + self.storage.name().unwrap_or(self.property.default_name()) + } + + pub fn input_names(&self) -> &'static [&'static str] { + self.property.input_names() + } + + pub fn target_name(&self) -> &'static str { + self.property.target_name() + } +} + +impl Dataset for PureDataset { + fn inputs(&self) -> ArrayView2<'_, f64> { + self.inputs() + } + + fn target(&self) -> ArrayView1<'_, f64> { + self.target() + } + + fn name(&self) -> &str { + self.name() + } + + fn input_names(&self) -> &'static [&'static str] { + self.input_names() + } + + fn target_name(&self) -> &'static str { + self.target_name() + } + + fn evaluate(&self, eos: &E) -> (Array1, Array1) { + self.property.evaluate(eos, self.inputs().view()) + } +} + +impl DatasetAD<1> for PureDataset { + fn evaluate_ad_const, const P: usize>( + &self, + names: [String; P], + parameters: &[f64], + inputs: ArrayView2, + ) -> (Array1, Array2, Array1) + where + T::Lifted>: Sync, + { + self.property.evaluate_ad::(names, parameters, inputs) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use std::io::Cursor; + + fn csv(s: &str) -> Cursor<&[u8]> { + Cursor::new(s.as_bytes()) + } + + #[test] + fn vapor_pressure_from_reader() { + let data = "\ +temperature_k,vapor_pressure_pa +300.0,3540.0 +350.0,41682.0 +400.0,245600.0 +"; + let ds = PureDataset::from_reader(PureProperty::VaporPressure, csv(data)).unwrap(); + + assert_eq!(ds.target().len(), 3); + assert_eq!(ds.inputs().nrows(), 3); + assert_eq!(ds.inputs().ncols(), 1); + assert_eq!(ds.inputs()[[0, 0]], 300.0); + assert_eq!(ds.inputs()[[1, 0]], 350.0); + assert_eq!(ds.inputs()[[2, 0]], 400.0); + assert_eq!(ds.target()[0], 3540.0); + assert_eq!(ds.target()[1], 41682.0); + assert_eq!(ds.target()[2], 245600.0); + assert_eq!(ds.name(), "vapor pressure"); + } + + #[test] + fn vapor_pressure_from_records() { + let records = vec![ + VaporPressureRecord { + temperature_k: 300.0, + vapor_pressure_pa: 3540.0, + }, + VaporPressureRecord { + temperature_k: 350.0, + vapor_pressure_pa: 41682.0, + }, + ]; + let ds = PureDataset::vapor_pressure(records); + assert_eq!(ds.inputs()[[0, 0]], 300.0); + assert_eq!(ds.target()[1], 41682.0); + } + + #[test] + fn liquid_density_from_reader() { + let data = "\ +temperature_k,pressure_pa,liquid_density_kmol_m3 +300.0,101325.0,15.2 +320.0,200000.0,14.8 +"; + let ds = PureDataset::from_reader(PureProperty::LiquidDensity, csv(data)).unwrap(); + + assert_eq!(ds.inputs().nrows(), 2); + assert_eq!(ds.inputs().ncols(), 2); + assert_eq!(ds.inputs()[[0, 0]], 300.0); + assert_eq!(ds.inputs()[[0, 1]], 101325.0); + assert_eq!(ds.inputs()[[1, 0]], 320.0); + assert_eq!(ds.inputs()[[1, 1]], 200000.0); + assert_eq!(ds.target()[0], 15.2); + assert_eq!(ds.target()[1], 14.8); + assert_eq!(ds.name(), "liquid density"); + } + + #[test] + fn liquid_density_from_records() { + let records = vec![LiquidDensityRecord { + temperature_k: 300.0, + pressure_pa: 101325.0, + liquid_density_kmol_m3: 15.2, + }]; + let ds = PureDataset::liquid_density(records); + assert_eq!(ds.inputs()[[0, 1]], 101325.0); + assert_eq!(ds.target()[0], 15.2); + } + + #[test] + fn equilibrium_liquid_density_from_reader() { + let data = "\ +temperature_k,liquid_density_kmol_m3 +290.0,15.5 +310.0,14.9 +330.0,14.1 +"; + let ds = + PureDataset::from_reader(PureProperty::EquilibriumLiquidDensity, csv(data)).unwrap(); + + assert_eq!(ds.inputs().ncols(), 1); + assert_eq!(ds.inputs().nrows(), 3); + assert_eq!(ds.inputs()[[2, 0]], 330.0); + assert_eq!(ds.target()[2], 14.1); + assert_eq!(ds.name(), "equilibrium liquid density"); + } + + #[test] + fn missing_column_returns_error() { + let data = "temperature_k\n300.0\n"; + let result = PureDataset::from_reader(PureProperty::VaporPressure, csv(data)); + assert!(result.is_err()); + } + + #[test] + fn wrong_type_returns_error() { + let data = "temperature_k,vapor_pressure_pa\n300.0,not_a_number\n"; + let result = PureDataset::from_reader(PureProperty::VaporPressure, csv(data)); + assert!(result.is_err()); + } +} diff --git a/crates/feos-core/src/ad/mod.rs b/crates/feos-core/src/ad/mod.rs index f48fb153c..b8d94ab6c 100644 --- a/crates/feos-core/src/ad/mod.rs +++ b/crates/feos-core/src/ad/mod.rs @@ -1,395 +1,74 @@ -use crate::DensityInitialization::Liquid; -use crate::density_iteration::density_iteration; -use crate::{Composition, FeosResult, PhaseEquilibrium, ReferenceSystem, Residual}; -use nalgebra::{Const, SVector, U1, U2}; -#[cfg(feature = "rayon")] -use ndarray::{Array1, Array2, ArrayView2, Zip}; -use num_dual::{Derivative, DualNum, DualSVec, DualStruct, first_derivative, partial2}; -use quantity::{Density, Pressure, Temperature}; -#[cfg(feature = "rayon")] -use quantity::{KELVIN, KILO, METER, MOL, PASCAL}; +use crate::Residual; +use nalgebra::{Const, DefaultAllocator, Dim, U1, allocator::Allocator}; +use num_dual::{Derivative, DualNum, DualSVec}; -type Gradient = DualSVec; +#[cfg(feature = "ndarray")] +mod dataset; +mod properties; +#[cfg(feature = "ndarray")] +pub use dataset::*; +pub use properties::*; -/// A model that can be evaluated with derivatives of its parameters. -pub trait ParametersAD: for<'a> From<&'a [f64]> + Residual> { - /// Return a mutable reference to the parameter named by `index` from the parameter set. - fn index_parameters_mut<'a, const P: usize>( - eos: &'a mut Self::Lifted>, - index: &str, - ) -> &'a mut Gradient

; - - /// Return the parameters with the appropriate derivatives. - fn named_derivatives( - &self, - parameter_names: [&str; P], - ) -> Self::Lifted> { - let mut eos = self.lift::>(); - for (i, p) in parameter_names.into_iter().enumerate() { - Self::index_parameters_mut(&mut eos, p).eps = - Derivative::derivative_generic(Const::

, U1, i) - } - eos - } -} - -/// Properties that can be evaluated with derivatives of model parameters. -pub trait PropertiesAD { - fn vapor_pressure( - &self, - temperature: Temperature, - ) -> FeosResult>> - where - Self: Residual>, - { - let eos_f64 = self.re(); - let (_, [vapor_density, liquid_density]) = - PhaseEquilibrium::pure_t(&eos_f64, temperature, None, Default::default())?; - - // implicit differentiation is implemented here instead of just calling pure_t with dual - // numbers, because for the first derivative, we can avoid calculating density derivatives. - let v1 = 1.0 / liquid_density.to_reduced(); - let v2 = 1.0 / vapor_density.to_reduced(); - let t = temperature.into_reduced(); - let (a1, a2) = { - let t = Gradient::from(t); - let v1 = Gradient::from(v1); - let v2 = Gradient::from(v2); - let x = Self::pure_molefracs(); - - let a1 = self.residual_helmholtz_energy(t, v1, &x); - let a2 = self.residual_helmholtz_energy(t, v2, &x); - (a1, a2) - }; - - let p = -(a1 - a2 + t * (v2 / v1).ln()) / (v1 - v2); - Ok(Pressure::from_reduced(p)) - } - - fn boiling_temperature( - &self, - pressure: Pressure, - ) -> FeosResult>> - where - Self: Residual>, - { - let eos_f64 = self.re(); - let (temperature, [vapor_density, liquid_density]) = - PhaseEquilibrium::pure_p(&eos_f64, pressure, None, Default::default())?; - - // implicit differentiation is implemented here instead of just calling pure_t with dual - // numbers, because for the first derivative, we can avoid calculating density derivatives. - let t = temperature.into_reduced(); - let v1 = 1.0 / liquid_density.to_reduced(); - let v2 = 1.0 / vapor_density.to_reduced(); - let p = pressure.into_reduced(); - let t = Gradient::from(t); - let t = t + { - let v1 = Gradient::from(v1); - let v2 = Gradient::from(v2); - let p = Gradient::from(p); - let x = Self::pure_molefracs(); - - let residual_entropy = |v| { - let (a, s) = first_derivative( - partial2( - |t, &v, x| self.lift().residual_helmholtz_energy(t, v, x), - &v, - &x, - ), - t, - ); - (a, -s) - }; - let (a1, s1) = residual_entropy(v1); - let (a2, s2) = residual_entropy(v2); - - let ln_rho = (v1 / v2).ln(); - (p * (v2 - v1) + (a2 - a1 + t * ln_rho)) / (s2 - s1 - ln_rho) - }; - Ok(Temperature::from_reduced(t)) - } - - fn equilibrium_liquid_density( - &self, - temperature: Temperature, - ) -> FeosResult<(Pressure>, Density>)> - where - Self: Residual>, - { - let t = Temperature::from_inner(&temperature); - PhaseEquilibrium::pure_t(self, t, None, Default::default()).map(|(p, [_, rho])| (p, rho)) - } - - fn liquid_density( - &self, - temperature: Temperature, - pressure: Pressure, - ) -> FeosResult>> - where - Self: Residual>, - { - let x = Self::pure_molefracs(); - let t = Temperature::from_inner(&temperature); - let p = Pressure::from_inner(&pressure); - density_iteration(self, t, p, &x, Some(Liquid)) - } - - #[cfg(feature = "rayon")] - fn vapor_pressure_parallel( - parameter_names: [String; P], - parameters: ArrayView2, - input: ArrayView2, - ) -> (Array1, Array2, Array1) - where - Self: ParametersAD<1>, - { - parallelize::<_, Self, _, _>( - parameter_names, - parameters, - input, - |eos: &Self::Lifted>, inp| { - eos.vapor_pressure(inp[0] * KELVIN) - .map(|p| p.convert_into(PASCAL)) - }, - ) - } - - #[cfg(feature = "rayon")] - fn boiling_temperature_parallel( - parameter_names: [String; P], - parameters: ArrayView2, - input: ArrayView2, - ) -> (Array1, Array2, Array1) - where - Self: ParametersAD<1>, - { - parallelize::<_, Self, _, _>( - parameter_names, - parameters, - input, - |eos: &Self::Lifted>, inp| { - eos.boiling_temperature(inp[0] * PASCAL) - .map(|p| p.convert_into(KELVIN)) - }, - ) - } +pub(crate) type Gradient = DualSVec; - #[cfg(feature = "rayon")] - fn liquid_density_parallel( - parameter_names: [String; P], - parameters: ArrayView2, - input: ArrayView2, - ) -> (Array1, Array2, Array1) - where - Self: ParametersAD<1>, - { - parallelize::<_, Self, _, _>( - parameter_names, - parameters, - input, - |eos: &Self::Lifted>, inp| { - eos.liquid_density(inp[0] * KELVIN, inp[1] * PASCAL) - .map(|d| d.convert_into(KILO * MOL / (METER * METER * METER))) - }, - ) - } - - #[cfg(feature = "rayon")] - fn equilibrium_liquid_density_parallel( - parameter_names: [String; P], - parameters: ArrayView2, - input: ArrayView2, - ) -> (Array1, Array2, Array1) - where - Self: ParametersAD<1>, - { - parallelize::<_, Self, _, _>( - parameter_names, - parameters, - input, - |eos: &Self::Lifted>, inp| { - eos.equilibrium_liquid_density(inp[0] * KELVIN) - .map(|(_, d)| d.convert_into(KILO * MOL / (METER * METER * METER))) - }, - ) - } - - fn bubble_point_pressure>( - &self, - temperature: Temperature, - pressure: Option, - liquid_molefracs: X, - ) -> FeosResult>> - where - Self: Residual>, - { - let eos_f64 = self.re(); - let (liquid_molefracs, _) = liquid_molefracs.into_molefracs(&eos_f64)?; - let vle = PhaseEquilibrium::bubble_point( - &eos_f64, - temperature, - liquid_molefracs, - pressure, - None, - Default::default(), - )?; - - // implicit differentiation is implemented here instead of just calling bubble_point with dual - // numbers, because for the first derivative, we can avoid calculating density derivatives. - let v_l = 1.0 / vle.liquid().density.to_reduced(); - let v_v = 1.0 / vle.vapor().density.to_reduced(); - let y = &vle.vapor().molefracs; - let y: SVector<_, 2> = SVector::from_fn(|i, _| y[i]); - let t = temperature.into_reduced(); - let (a_l, a_v, v_l, v_v) = { - let t = Gradient::from(t); - let v_l = Gradient::from(v_l); - let v_v = Gradient::from(v_v); - let y = y.map(Gradient::from); - let x = liquid_molefracs.map(Gradient::from); - - let a_v = self.residual_helmholtz_energy(t, v_v, &y); - let (p_l, mu_res_l, dp_l, dmu_l) = self.dmu_dv(t, v_l, &x); - let vi_l = dmu_l / dp_l; - let v_l = vi_l.dot(&y); - let a_l = (mu_res_l - vi_l * p_l).dot(&y); - (a_l, a_v, v_l, v_v) - }; - let rho_l = vle.liquid().partial_density().to_reduced(); - let rho_l = [rho_l[0], rho_l[1]]; - let rho_v = vle.vapor().partial_density().to_reduced(); - let rho_v = [rho_v[0], rho_v[1]]; - let p = -(a_v - a_l - + t * (y[0] * (rho_v[0] / rho_l[0]).ln() + y[1] * (rho_v[1] / rho_l[1]).ln() - 1.0)) - / (v_v - v_l); - Ok(Pressure::from_reduced(p)) - } - - fn dew_point_pressure>( - &self, - temperature: Temperature, - pressure: Option, - vapor_molefracs: X, - ) -> FeosResult>> - where - Self: Residual>, - { - let eos_f64 = self.re(); - let (vapor_molefracs, _) = vapor_molefracs.into_molefracs(&eos_f64)?; - let vle = PhaseEquilibrium::dew_point( - &eos_f64, - temperature, - vapor_molefracs, - pressure, - None, - Default::default(), - )?; - - // implicit differentiation is implemented here instead of just calling dew_point with dual - // numbers, because for the first derivative, we can avoid calculating density derivatives. - let v_l = 1.0 / vle.liquid().density.to_reduced(); - let v_v = 1.0 / vle.vapor().density.to_reduced(); - let x = &vle.liquid().molefracs; - let x: SVector<_, 2> = SVector::from_fn(|i, _| x[i]); - let t = temperature.into_reduced(); - let (a_l, a_v, v_l, v_v) = { - let t = Gradient::from(t); - let v_l = Gradient::from(v_l); - let v_v = Gradient::from(v_v); - let x = x.map(Gradient::from); - let y = vapor_molefracs.map(Gradient::from); - - let a_l = self.residual_helmholtz_energy(t, v_l, &x); - let (p_v, mu_res_v, dp_v, dmu_v) = self.dmu_dv(t, v_v, &y); - let vi_v = dmu_v / dp_v; - let v_v = vi_v.dot(&x); - let a_v = (mu_res_v - vi_v * p_v).dot(&x); - (a_l, a_v, v_l, v_v) - }; - let rho_l = vle.liquid().partial_density().to_reduced(); - let rho_l = [rho_l[0], rho_l[1]]; - let rho_v = vle.vapor().partial_density().to_reduced(); - let rho_v = [rho_v[0], rho_v[1]]; - let p = -(a_l - a_v - + t * (x[0] * (rho_l[0] / rho_v[0]).ln() + x[1] * (rho_l[1] / rho_v[1]).ln() - 1.0)) - / (v_l - v_v); - Ok(Pressure::from_reduced(p)) +/// A model that can be evaluated with derivatives of its parameters. +pub trait ParametersAD: Residual +where + DefaultAllocator: Allocator, +{ + /// Build the model by requesting each parameter by name. + /// + /// Call `f(name, differentiable)` for each parameter. The order of calls + /// defines the canonical parameter order. + /// + /// Set `differentiable` to `false` for fixed parameters. + fn build + Copy>( + f: impl FnMut(&'static str, bool) -> D, + ) -> Self::Lifted; + + /// Canonical parameter names in the order defined by [`build`](Self::build). + fn parameter_names() -> Vec<&'static str> { + let mut names = Vec::new(); + let _ = Self::build(|name, _| { + names.push(name); + 0.0 + }); + names } - #[cfg(feature = "rayon")] - fn bubble_point_pressure_parallel( - parameter_names: [String; P], - parameters: ArrayView2, - input: ArrayView2, - ) -> (Array1, Array2, Array1) - where - Self: ParametersAD<2>, - { - parallelize::<_, Self, _, _>( - parameter_names, - parameters, - input, - |eos: &Self::Lifted>, inp| { - eos.bubble_point_pressure(inp[0] * KELVIN, Some(inp[2] * PASCAL), inp[1]) - .map(|p| p.convert_into(PASCAL)) - }, - ) + /// Parameter names that can be differentiated, in canonical order. + fn differentiable_parameters() -> Vec<&'static str> { + let mut names = Vec::new(); + let _ = Self::build(|name, differentiable| { + if differentiable { + names.push(name); + } + 0.0 + }); + names } - #[cfg(feature = "rayon")] - fn dew_point_pressure_parallel( - parameter_names: [String; P], - parameters: ArrayView2, - input: ArrayView2, - ) -> (Array1, Array2, Array1) - where - Self: ParametersAD<2>, - { - parallelize::<_, Self, _, _>( - parameter_names, - parameters, - input, - |eos: &Self::Lifted>, inp| { - eos.dew_point_pressure(inp[0] * KELVIN, Some(inp[2] * PASCAL), inp[1]) - .map(|p| p.convert_into(PASCAL)) - }, - ) + /// Construct the model with derivative seeds for the `P` named parameters. + /// + /// - `parameter_values`: all parameter values in the canonical order + /// defined by [`build`](Self::build). + /// - `derivative_names`: names of the parameters to differentiate with + /// respect to. Gradient component `i` corresponds to + /// `derivative_names[i]`. + fn seed_derivatives( + parameter_values: &[f64], + derivative_names: [&str; P], + ) -> Self::Lifted> { + let mut idx = 0; + Self::build(|name, _differentiable| { + let i = idx; + idx += 1; + let mut d = Gradient::

::from(parameter_values[i]); + if let Some(seed_idx) = derivative_names.iter().position(|&n| n == name) { + d.eps = + Derivative::<_, _, Const

, _>::derivative_generic(Const::

, U1, seed_idx); + } + d + }) } } - -impl PropertiesAD for T {} - -#[cfg(feature = "rayon")] -fn parallelize, const N: usize, const P: usize>( - parameter_names: [String; P], - parameters: ArrayView2, - input: ArrayView2, - f: F, -) -> (Array1, Array2, Array1) -where - F: Fn(&E::Lifted>, &[f64]) -> FeosResult> + Sync, -{ - let parameter_names = parameter_names.each_ref().map(|s| s as &str); - let value_dual = Zip::from(parameters.rows()) - .and(input.rows()) - .par_map_collect(|par, inp| { - let par = par.as_slice().expect("Parameter array is not contiguous!"); - let inp = inp.as_slice().expect("Input array is not contiguous!"); - let eos = E::from(par).named_derivatives(parameter_names); - f(&eos, inp) - }); - let status = value_dual.iter().map(|p| p.is_ok()).collect(); - let value_dual: Array1<_> = value_dual.into_iter().flatten().collect(); - let mut value = Array1::zeros(value_dual.len()); - let mut grad = Array2::zeros([value_dual.len(), P]); - Zip::from(grad.rows_mut()) - .and(&mut value) - .and(&value_dual) - .for_each(|mut grad, p, p_dual| { - *p = p_dual.re; - let eps = p_dual.eps.unwrap_generic(Const::

, U1).data.0[0].to_vec(); - grad.assign(&Array1::from(eps)); - }); - (value, grad, status) -} diff --git a/crates/feos-core/src/ad/properties/boiling_temperature.rs b/crates/feos-core/src/ad/properties/boiling_temperature.rs new file mode 100644 index 000000000..c597fb30e --- /dev/null +++ b/crates/feos-core/src/ad/properties/boiling_temperature.rs @@ -0,0 +1,71 @@ +use super::PropertyAD; +use crate::ad::Gradient; +use crate::{FeosResult, PhaseEquilibrium, ReferenceSystem, Residual}; +use nalgebra::allocator::Allocator; +use nalgebra::{DefaultAllocator, U1}; +use num_dual::{DualNum, DualStruct, Gradients, first_derivative, partial2}; +use quantity::{_Temperature, KELVIN, PASCAL, Pressure, Temperature}; + +/// Boiling temperature of a pure component as function of pressure. +pub struct BoilingTemperature(pub Pressure); + +impl<'a> From<&'a [f64]> for BoilingTemperature { + fn from(value: &'a [f64]) -> Self { + Self(value[0] * PASCAL) + } +} + +impl PropertyAD for BoilingTemperature +where + DefaultAllocator: Allocator + Allocator + Allocator, +{ + type Unit = _Temperature; + const REFERENCE: Temperature = KELVIN; + + fn evaluate, D: DualNum + Copy>( + &self, + eos: &E, + ) -> FeosResult> { + let p = Pressure::from_inner(&self.0); + PhaseEquilibrium::pure_p(eos, p, None, Default::default()).map(|(t, _)| t) + } + + fn evaluate_gradient>, const P: usize>( + &self, + eos: &E, + ) -> FeosResult, Self::Unit>> { + let eos_f64 = eos.re(); + let (temperature, [vapor_density, liquid_density]) = + PhaseEquilibrium::pure_p(&eos_f64, self.0, None, Default::default())?; + + let t = temperature.into_reduced(); + let v1 = 1.0 / liquid_density.to_reduced(); + let v2 = 1.0 / vapor_density.to_reduced(); + let p = self.0.into_reduced(); + let t = Gradient::from(t); + let t = t + { + let v1 = Gradient::from(v1); + let v2 = Gradient::from(v2); + let p = Gradient::from(p); + let x = E::pure_molefracs(); + + let residual_entropy = |v| { + let (a, s) = first_derivative( + partial2( + |t, &v, x| eos.lift().residual_helmholtz_energy(t, v, x), + &v, + &x, + ), + t, + ); + (a, -s) + }; + let (a1, s1) = residual_entropy(v1); + let (a2, s2) = residual_entropy(v2); + + let ln_rho = (v1 / v2).ln(); + (p * (v2 - v1) + (a2 - a1 + t * ln_rho)) / (s2 - s1 - ln_rho) + }; + Ok(Temperature::from_reduced(t)) + } +} diff --git a/crates/feos-core/src/ad/properties/bubble_point_pressure.rs b/crates/feos-core/src/ad/properties/bubble_point_pressure.rs new file mode 100644 index 000000000..9e838b9f7 --- /dev/null +++ b/crates/feos-core/src/ad/properties/bubble_point_pressure.rs @@ -0,0 +1,91 @@ +use super::PropertyAD; +use crate::Contributions; +use crate::ad::Gradient; +use crate::{Composition, FeosResult, PhaseEquilibrium, ReferenceSystem, Residual}; +use nalgebra::allocator::Allocator; +use nalgebra::{DefaultAllocator, U1}; +use num_dual::{DualNum, DualStruct, Gradients}; +use quantity::{_Pressure, KELVIN, PASCAL, Pressure, Temperature}; + +/// Bubble point pressure of a binary mixture as function of temperature and +/// molefracs of the first component. +/// +/// An initial value for the pressure can be passed as optional argument to +/// increase robustness and speed. +pub struct BubblePointPressure(pub Temperature, pub f64, pub Option); + +impl<'a> From<&'a [f64]> for BubblePointPressure { + fn from(value: &'a [f64]) -> Self { + Self(value[0] * KELVIN, value[1], Some(value[2] * PASCAL)) + } +} + +impl PropertyAD for BubblePointPressure +where + DefaultAllocator: Allocator + Allocator + Allocator, + f64: Composition, +{ + type Unit = _Pressure; + const REFERENCE: Pressure = PASCAL; + + fn evaluate, D: DualNum + Copy>( + &self, + eos: &E, + ) -> FeosResult> + where + DefaultAllocator: Allocator + Allocator + Allocator, + { + let t = Temperature::from_inner(&self.0); + let p = Option::from_inner(&self.2); + let (x, _) = self.1.into_molefracs(&eos.re())?; + let x = x.map(D::from); + let vle = PhaseEquilibrium::bubble_point(eos, t, x, p, None, Default::default())?; + Ok(vle.vapor().pressure(Contributions::Total)) + } + + fn evaluate_gradient>, const P: usize>( + &self, + eos: &E, + ) -> FeosResult, Self::Unit>> + where + DefaultAllocator: Allocator + Allocator + Allocator, + { + let eos_f64 = eos.re(); + let (liquid_molefracs, _) = self.1.into_molefracs(&eos_f64)?; + let vle = PhaseEquilibrium::bubble_point( + &eos_f64, + self.0, + &liquid_molefracs, + self.2, + None, + Default::default(), + )?; + + let v_l = 1.0 / vle.liquid().density.to_reduced(); + let v_v = 1.0 / vle.vapor().density.to_reduced(); + let y = &vle.vapor().molefracs; + let t = self.0.into_reduced(); + let (a_l, a_v, v_l, v_v) = { + let t = Gradient::from(t); + let v_l = Gradient::from(v_l); + let v_v = Gradient::from(v_v); + let y = y.map(Gradient::from); + let x = liquid_molefracs.map(Gradient::from); + + let a_v = eos.residual_helmholtz_energy(t, v_v, &y); + let (p_l, mu_res_l, dp_l, dmu_l) = eos.dmu_dv(t, v_l, &x); + let vi_l = dmu_l / dp_l; + let v_l = vi_l.dot(&y); + let a_l = (mu_res_l - vi_l * p_l).dot(&y); + (a_l, a_v, v_l, v_v) + }; + let rho_l = vle.liquid().partial_density().to_reduced(); + let rho_l = [rho_l[0], rho_l[1]]; + let rho_v = vle.vapor().partial_density().to_reduced(); + let rho_v = [rho_v[0], rho_v[1]]; + let p = -(a_v - a_l + + t * (y[0] * (rho_v[0] / rho_l[0]).ln() + y[1] * (rho_v[1] / rho_l[1]).ln() - 1.0)) + / (v_v - v_l); + Ok(Pressure::from_reduced(p)) + } +} diff --git a/crates/feos-core/src/ad/properties/dew_point_pressure.rs b/crates/feos-core/src/ad/properties/dew_point_pressure.rs new file mode 100644 index 000000000..4ce733484 --- /dev/null +++ b/crates/feos-core/src/ad/properties/dew_point_pressure.rs @@ -0,0 +1,90 @@ +use super::PropertyAD; +use crate::ad::Gradient; +use crate::{Composition, Contributions, FeosResult, PhaseEquilibrium, ReferenceSystem, Residual}; +use nalgebra::allocator::Allocator; +use nalgebra::{DefaultAllocator, U1}; +use num_dual::{DualNum, DualStruct, Gradients}; +use quantity::{_Pressure, KELVIN, PASCAL, Pressure, Temperature}; + +/// Dew point pressure of a binary mixture as function of temperature and +/// molefracs of the first component. +/// +/// An initial value for the pressure can be passed as optional argument to +/// increase robustness and speed. +pub struct DewPointPressure(pub Temperature, pub f64, pub Option); + +impl<'a> From<&'a [f64]> for DewPointPressure { + fn from(value: &'a [f64]) -> Self { + Self(value[0] * KELVIN, value[1], Some(value[2] * PASCAL)) + } +} + +impl PropertyAD for DewPointPressure +where + DefaultAllocator: Allocator + Allocator + Allocator, + f64: Composition, +{ + type Unit = _Pressure; + const REFERENCE: Pressure = PASCAL; + + fn evaluate, D: DualNum + Copy>( + &self, + eos: &E, + ) -> FeosResult> + where + DefaultAllocator: Allocator + Allocator + Allocator, + { + let t = Temperature::from_inner(&self.0); + let p = Option::from_inner(&self.2); + let (y, _) = self.1.into_molefracs(&eos.re())?; + let y = y.map(D::from); + let vle = PhaseEquilibrium::dew_point(eos, t, y, p, None, Default::default())?; + Ok(vle.vapor().pressure(Contributions::Total)) + } + + fn evaluate_gradient>, const P: usize>( + &self, + eos: &E, + ) -> FeosResult, Self::Unit>> + where + DefaultAllocator: Allocator + Allocator + Allocator, + { + let eos_f64 = eos.re(); + let (vapor_molefracs, _) = self.1.into_molefracs(&eos_f64)?; + let vle = PhaseEquilibrium::dew_point( + &eos_f64, + self.0, + &vapor_molefracs, + self.2, + None, + Default::default(), + )?; + + let v_l = 1.0 / vle.liquid().density.to_reduced(); + let v_v = 1.0 / vle.vapor().density.to_reduced(); + let x = &vle.liquid().molefracs; + let t = self.0.into_reduced(); + let (a_l, a_v, v_l, v_v) = { + let t = Gradient::from(t); + let v_l = Gradient::from(v_l); + let v_v = Gradient::from(v_v); + let x = x.map(Gradient::from); + let y = vapor_molefracs.map(Gradient::from); + + let a_l = eos.residual_helmholtz_energy(t, v_l, &x); + let (p_v, mu_res_v, dp_v, dmu_v) = eos.dmu_dv(t, v_v, &y); + let vi_v = dmu_v / dp_v; + let v_v = vi_v.dot(&x); + let a_v = (mu_res_v - vi_v * p_v).dot(&x); + (a_l, a_v, v_l, v_v) + }; + let rho_l = vle.liquid().partial_density().to_reduced(); + let rho_l = [rho_l[0], rho_l[1]]; + let rho_v = vle.vapor().partial_density().to_reduced(); + let rho_v = [rho_v[0], rho_v[1]]; + let p = -(a_l - a_v + + t * (x[0] * (rho_l[0] / rho_v[0]).ln() + x[1] * (rho_l[1] / rho_v[1]).ln() - 1.0)) + / (v_l - v_v); + Ok(Pressure::from_reduced(p)) + } +} diff --git a/crates/feos-core/src/ad/properties/enthalpy_of_vaporization.rs b/crates/feos-core/src/ad/properties/enthalpy_of_vaporization.rs new file mode 100644 index 000000000..ee0d3d3aa --- /dev/null +++ b/crates/feos-core/src/ad/properties/enthalpy_of_vaporization.rs @@ -0,0 +1,34 @@ +use super::PropertyAD; +use crate::{FeosResult, PhaseEquilibrium, Residual}; +use nalgebra::allocator::Allocator; +use nalgebra::{DefaultAllocator, U1}; +use num_dual::{DualNum, DualStruct, Gradients}; +use quantity::{_MolarEnergy, KELVIN, MolarEnergy, Temperature}; + +/// Enthalpy of vaporization of a pure component as function of temperature. +pub struct EnthalpyOfVaporization(pub Temperature); + +impl<'a> From<&'a [f64]> for EnthalpyOfVaporization { + fn from(value: &'a [f64]) -> Self { + Self(value[0] * KELVIN) + } +} + +impl PropertyAD for EnthalpyOfVaporization +where + DefaultAllocator: Allocator + Allocator + Allocator, +{ + type Unit = _MolarEnergy; + const REFERENCE: MolarEnergy = MolarEnergy::new(1.0); + + fn evaluate, D: DualNum + Copy>( + &self, + eos: &E, + ) -> FeosResult> { + let t = Temperature::from_inner(&self.0); + let vle = PhaseEquilibrium::pure(eos, t, None, Default::default())?; + let h_v = vle.vapor().residual_molar_enthalpy(); + let h_l = vle.liquid().residual_molar_enthalpy(); + Ok(h_v - h_l) + } +} diff --git a/crates/feos-core/src/ad/properties/equilibrium_liquid_density.rs b/crates/feos-core/src/ad/properties/equilibrium_liquid_density.rs new file mode 100644 index 000000000..85ded0885 --- /dev/null +++ b/crates/feos-core/src/ad/properties/equilibrium_liquid_density.rs @@ -0,0 +1,31 @@ +use super::PropertyAD; +use crate::{FeosResult, PhaseEquilibrium, Residual}; +use nalgebra::allocator::Allocator; +use nalgebra::{DefaultAllocator, U1}; +use num_dual::{DualNum, DualStruct, Gradients}; +use quantity::{_Density, Density, KELVIN, Temperature}; + +/// Equilibrium liquid density of a pure component as function of temperature. +pub struct EquilibriumLiquidDensity(pub Temperature); + +impl<'a> From<&'a [f64]> for EquilibriumLiquidDensity { + fn from(value: &'a [f64]) -> Self { + Self(value[0] * KELVIN) + } +} + +impl PropertyAD for EquilibriumLiquidDensity +where + DefaultAllocator: Allocator + Allocator + Allocator, +{ + type Unit = _Density; + const REFERENCE: Density = Density::new(1000.0); + + fn evaluate, D: DualNum + Copy>( + &self, + eos: &E, + ) -> FeosResult> { + let t = Temperature::from_inner(&self.0); + PhaseEquilibrium::pure_t(eos, t, None, Default::default()).map(|(_, [_, r])| r) + } +} diff --git a/crates/feos-core/src/ad/properties/liquid_density.rs b/crates/feos-core/src/ad/properties/liquid_density.rs new file mode 100644 index 000000000..f8992c3fe --- /dev/null +++ b/crates/feos-core/src/ad/properties/liquid_density.rs @@ -0,0 +1,35 @@ +use super::PropertyAD; +use crate::DensityInitialization::Liquid; +use crate::density_iteration::density_iteration; +use crate::{FeosResult, Residual}; +use nalgebra::allocator::Allocator; +use nalgebra::{DefaultAllocator, Dim}; +use num_dual::{DualNum, DualStruct}; +use quantity::{_Density, Density, KELVIN, PASCAL, Pressure, Temperature}; + +/// Liquid density of a pure component as function of temperature and pressure. +pub struct LiquidDensity(pub Temperature, pub Pressure); + +impl<'a> From<&'a [f64]> for LiquidDensity { + fn from(value: &'a [f64]) -> Self { + Self(value[0] * KELVIN, value[1] * PASCAL) + } +} + +impl PropertyAD for LiquidDensity +where + DefaultAllocator: Allocator, +{ + type Unit = _Density; + const REFERENCE: Density = Density::new(1000.0); + + fn evaluate, D: DualNum + Copy>( + &self, + eos: &E, + ) -> FeosResult> { + let x = E::pure_molefracs(); + let t = Temperature::from_inner(&self.0); + let p = Pressure::from_inner(&self.1); + density_iteration(eos, t, p, &x, Some(Liquid)) + } +} diff --git a/crates/feos-core/src/ad/properties/mod.rs b/crates/feos-core/src/ad/properties/mod.rs new file mode 100644 index 000000000..660716e1d --- /dev/null +++ b/crates/feos-core/src/ad/properties/mod.rs @@ -0,0 +1,189 @@ +use super::Gradient; +use crate::{FeosResult, Residual}; +use nalgebra::{DefaultAllocator, Dim, allocator::Allocator}; +#[cfg(feature = "ndarray")] +use ndarray::{Array1, Array2, ArrayView2}; +use num_dual::DualNum; +use quantity::Quantity; + +mod boiling_temperature; +mod bubble_point_pressure; +mod dew_point_pressure; +mod enthalpy_of_vaporization; +mod equilibrium_liquid_density; +mod liquid_density; +mod residual_isobaric_heat_capacity; +mod vapor_pressure; + +pub use boiling_temperature::BoilingTemperature; +pub use bubble_point_pressure::BubblePointPressure; +pub use dew_point_pressure::DewPointPressure; +pub use enthalpy_of_vaporization::EnthalpyOfVaporization; +pub use equilibrium_liquid_density::EquilibriumLiquidDensity; +pub use liquid_density::LiquidDensity; +pub use residual_isobaric_heat_capacity::ResidualIsobaricHeatCapacity; +pub use vapor_pressure::VaporPressure; + +/// Properties that can be rapidly evaluated in parallel together with +/// their gradients with respect to model parameters +pub trait PropertyAD: for<'a> From<&'a [f64]> +where + DefaultAllocator: Allocator, +{ + type Unit; + const REFERENCE: Quantity; + + /// Evaluate the property for an arbitrary derivative. + fn evaluate, D: DualNum + Copy>( + &self, + eos: &E, + ) -> FeosResult>; + + /// Evaluate the property for the first derivative w.r.t. model parameters. + /// + /// This can be overridden if there is a more performant implementation than + /// the general implementation in `evaluate`. + fn evaluate_gradient>, const P: usize>( + &self, + eos: &E, + ) -> FeosResult, Self::Unit>> { + self.evaluate(eos) + } + + /// Evaluate the property for all inputs in parallel. + /// + /// Return the property values and the success of the calculations. + #[cfg(feature = "ndarray")] + fn evaluate_parallel + Sync>( + eos: &E, + input: ArrayView2, + ) -> (Array1, Array1) { + #[cfg(feature = "rayon")] + let values = ndarray::Zip::from(input.rows()).par_map_collect(|inp| { + let inp = inp.as_slice().expect("Input array is not contiguous!"); + Self::from(inp) + .evaluate(eos) + .map(|d| d.convert_into(Self::REFERENCE)) + }); + + #[cfg(not(feature = "rayon"))] + let values = ndarray::Zip::from(input.rows()).map_collect(|inp| { + let inp = inp.as_slice().expect("Input array is not contiguous!"); + Self::from(inp) + .evaluate(eos) + .map(|d| d.convert_into(Self::REFERENCE)) + }); + + let n = input.nrows(); + let status: Array1 = values.iter().map(|r| r.is_ok()).collect(); + let mut value = Array1::from_elem(n, f64::NAN); + for (i, result) in values.into_iter().enumerate() { + if let Ok(v) = result { + value[i] = v; + } + } + (value, status) + } + + /// Evaluate the property and its gradients for all inputs in parallel. + /// + /// Return the property values, the gradients, and the success of the calculations. + #[cfg(feature = "ndarray")] + fn evaluate_parallel_derivatives, const P: usize>( + parameter_names: [String; P], + parameters: &[f64], + input: ArrayView2, + ) -> (Array1, Array2, Array1) + where + E::Lifted>: Sync, + { + let parameter_names = parameter_names.each_ref().map(|s| s as &str); + let eos = E::seed_derivatives(parameters, parameter_names); + + #[cfg(feature = "rayon")] + let value_dual = ndarray::Zip::from(input.rows()).par_map_collect(|inp| { + let inp = inp.as_slice().expect("Input array is not contiguous!"); + Self::from(inp) + .evaluate_gradient(&eos) + .map(|d| d.convert_into(Self::REFERENCE)) + }); + + #[cfg(not(feature = "rayon"))] + let value_dual = ndarray::Zip::from(input.rows()).map_collect(|inp| { + let inp = inp.as_slice().expect("Input array is not contiguous!"); + Self::from(inp) + .evaluate_gradient(&eos) + .map(|d| d.convert_into(Self::REFERENCE)) + }); + + let n = input.nrows(); + let status = value_dual.iter().map(|p| p.is_ok()).collect(); + let mut value = Array1::from_elem(n, f64::NAN); + let mut grad = Array2::zeros([n, P]); + for (i, result) in value_dual.into_iter().enumerate() { + if let Ok(p_dual) = result { + value[i] = p_dual.re; + let eps = p_dual + .eps + .unwrap_generic(nalgebra::Const::

, nalgebra::U1); + for (g, &e) in grad.row_mut(i).iter_mut().zip(eps.data.0[0].iter()) { + *g = e; + } + } + } + (value, grad, status) + } + + /// Evaluate the property and its gradients for all inputs and parameters in parallel. + /// + /// Return the property values, the gradients, and the success of the calculations. + #[cfg(feature = "ndarray")] + fn evaluate_parallel_derivatives_params, const P: usize>( + parameter_names: [String; P], + parameters: ArrayView2, + input: ArrayView2, + ) -> (Array1, Array2, Array1) { + let parameter_names = parameter_names.each_ref().map(|s| s as &str); + + #[cfg(feature = "rayon")] + let value_dual = ndarray::Zip::from(parameters.rows()) + .and(input.rows()) + .par_map_collect(|par, inp| { + let par = par.as_slice().expect("Parameter array is not contiguous!"); + let inp = inp.as_slice().expect("Input array is not contiguous!"); + let eos = E::seed_derivatives(par, parameter_names); + Self::from(inp) + .evaluate_gradient(&eos) + .map(|d| d.convert_into(Self::REFERENCE)) + }); + + #[cfg(not(feature = "rayon"))] + let value_dual = ndarray::Zip::from(parameters.rows()) + .and(input.rows()) + .map_collect(|par, inp| { + let par = par.as_slice().expect("Parameter array is not contiguous!"); + let inp = inp.as_slice().expect("Input array is not contiguous!"); + let eos = E::seed_derivatives(par, parameter_names); + Self::from(inp) + .evaluate_gradient(&eos) + .map(|d| d.convert_into(Self::REFERENCE)) + }); + + let n = parameters.nrows(); + let status = value_dual.iter().map(|p| p.is_ok()).collect(); + let mut value = Array1::from_elem(n, f64::NAN); + let mut grad = Array2::zeros([n, P]); + for (i, result) in value_dual.into_iter().enumerate() { + if let Ok(p_dual) = result { + value[i] = p_dual.re; + let eps = p_dual + .eps + .unwrap_generic(nalgebra::Const::

, nalgebra::U1); + for (g, &e) in grad.row_mut(i).iter_mut().zip(eps.data.0[0].iter()) { + *g = e; + } + } + } + (value, grad, status) + } +} diff --git a/crates/feos-core/src/ad/properties/residual_isobaric_heat_capacity.rs b/crates/feos-core/src/ad/properties/residual_isobaric_heat_capacity.rs new file mode 100644 index 000000000..b8f962757 --- /dev/null +++ b/crates/feos-core/src/ad/properties/residual_isobaric_heat_capacity.rs @@ -0,0 +1,36 @@ +use super::PropertyAD; +use crate::DensityInitialization::Liquid; +use crate::{FeosResult, Residual, State}; +use nalgebra::DefaultAllocator; +use nalgebra::allocator::Allocator; +use num_dual::{DualNum, DualStruct, Gradients}; +use quantity::{_MolarEntropy, KELVIN, MolarEntropy, PASCAL, Pressure, Temperature}; + +/// Liquid residual isobaric heat capacity of a pure component as function of temperature +/// and pressure. +pub struct ResidualIsobaricHeatCapacity(pub Temperature, pub Pressure); + +impl<'a> From<&'a [f64]> for ResidualIsobaricHeatCapacity { + fn from(value: &'a [f64]) -> Self { + Self(value[0] * KELVIN, value[1] * PASCAL) + } +} + +impl PropertyAD for ResidualIsobaricHeatCapacity +where + DefaultAllocator: Allocator, +{ + type Unit = _MolarEntropy; + const REFERENCE: MolarEntropy = MolarEntropy::new(1.0); + + fn evaluate, D: DualNum + Copy>( + &self, + eos: &E, + ) -> FeosResult> { + let x = E::pure_molefracs(); + let t = Temperature::from_inner(&self.0); + let p = Pressure::from_inner(&self.1); + let state = State::new_npt(eos, t, p, x, Some(Liquid))?; + Ok(state.residual_molar_isobaric_heat_capacity()) + } +} diff --git a/crates/feos-core/src/ad/properties/vapor_pressure.rs b/crates/feos-core/src/ad/properties/vapor_pressure.rs new file mode 100644 index 000000000..8192d3ff0 --- /dev/null +++ b/crates/feos-core/src/ad/properties/vapor_pressure.rs @@ -0,0 +1,60 @@ +use super::PropertyAD; +use crate::ad::Gradient; +use crate::{FeosResult, PhaseEquilibrium, ReferenceSystem, Residual}; +use nalgebra::allocator::Allocator; +use nalgebra::{DefaultAllocator, U1}; +use num_dual::{DualNum, DualStruct, Gradients}; +use quantity::{_Pressure, KELVIN, PASCAL, Pressure, Temperature}; + +/// Vapor pressure of a pure component as function of temperature. +pub struct VaporPressure(pub Temperature); + +impl<'a> From<&'a [f64]> for VaporPressure { + fn from(value: &'a [f64]) -> Self { + Self(value[0] * KELVIN) + } +} + +impl PropertyAD for VaporPressure +where + DefaultAllocator: Allocator + Allocator + Allocator, +{ + type Unit = _Pressure; + const REFERENCE: Pressure = PASCAL; + + fn evaluate, D: DualNum + Copy>( + &self, + eos: &E, + ) -> FeosResult> { + let t = Temperature::from_inner(&self.0); + PhaseEquilibrium::pure_t(eos, t, None, Default::default()).map(|(p, _)| p) + } + + fn evaluate_gradient>, const P: usize>( + &self, + eos: &E, + ) -> FeosResult>> { + let eos_f64 = eos.re(); + let (_, [vapor_density, liquid_density]) = + PhaseEquilibrium::pure_t(&eos_f64, self.0, None, Default::default())?; + + // implicit differentiation is implemented here instead of just calling pure_t with dual + // numbers, because for the first derivative, we can avoid calculating density derivatives. + let v1 = 1.0 / liquid_density.to_reduced(); + let v2 = 1.0 / vapor_density.to_reduced(); + let t = self.0.into_reduced(); + let (a1, a2) = { + let t = Gradient::from(t); + let v1 = Gradient::from(v1); + let v2 = Gradient::from(v2); + let x = E::pure_molefracs(); + + let a1 = eos.residual_helmholtz_energy(t, v1, &x); + let a2 = eos.residual_helmholtz_energy(t, v2, &x); + (a1, a2) + }; + + let p = -(a1 - a2 + t * (v2 / v1).ln()) / (v1 - v2); + Ok(Pressure::from_reduced(p)) + } +} diff --git a/crates/feos-core/src/density_iteration.rs b/crates/feos-core/src/density_iteration.rs index aaee68ec1..e27ca0fba 100644 --- a/crates/feos-core/src/density_iteration.rs +++ b/crates/feos-core/src/density_iteration.rs @@ -85,10 +85,7 @@ where let t = Dual::from_re(temperature); let x = molefracs.map(Dual::from); let (a_res, da_res) = first_derivative( - |molar_volume| { - eos.lift() - .residual_helmholtz_energy(t, molar_volume, &x) - }, + |molar_volume| eos.lift().residual_helmholtz_energy(t, molar_volume, &x), molar_volume, ); a_res - da_res * molar_volume + temperature * density.ln() diff --git a/crates/feos-core/src/lib.rs b/crates/feos-core/src/lib.rs index cc280cf0a..156894648 100644 --- a/crates/feos-core/src/lib.rs +++ b/crates/feos-core/src/lib.rs @@ -24,7 +24,7 @@ macro_rules! log_result { } } -mod ad; +pub mod ad; pub mod cubic; mod density_iteration; mod equation_of_state; @@ -32,7 +32,6 @@ mod errors; pub mod parameter; mod phase_equilibria; mod state; -pub use ad::{ParametersAD, PropertiesAD}; pub use equation_of_state::{ EntropyScaling, EquationOfState, IdealGas, IdealGasAD, Molarweight, NoResidual, Residual, ResidualDyn, Subset, Total, diff --git a/crates/feos-derive/src/dft.rs b/crates/feos-derive/src/dft.rs index a9ceb4f6a..d7dcd8e63 100644 --- a/crates/feos-derive/src/dft.rs +++ b/crates/feos-derive/src/dft.rs @@ -1,4 +1,4 @@ -use crate::{implement, OPT_IMPLS}; +use crate::{OPT_IMPLS, implement}; use quote::quote; use syn::DeriveInput; diff --git a/crates/feos/benches/README.md b/crates/feos/benches/README.md index b1921cd1a..4d0f3b756 100644 --- a/crates/feos/benches/README.md +++ b/crates/feos/benches/README.md @@ -7,11 +7,18 @@ For example, to run the benchmarks in `dual_numbers`, which uses PC-SAFT, use ``` cargo bench --profile=release-lto --bench=dual_numbers -``` +``` + +The static-vs-dynamic vector dual benchmark uses mimalloc. Run it with LTO via + +``` +cargo bench --profile=release-lto --bench=dual_static_vs_dynamic --features=pcsaft +``` |Name|Description| |--|--| |`dual_numbers`|Helmholtz energy function evaluated using `StateHD` with different dual number types using the PC-SAFT equation of state.| +|`dual_static_vs_dynamic`|PC-SAFT-like pure-component Helmholtz energy expression evaluated with static (`DualSVec64

`) and dynamic (`DualDVec64`) vector dual numbers.| |`dual_numbers_saftvrmie`|Helmholtz energy function evaluated using `StateHD` with different dual number types using the SAFT-VR-Mie equation of state.| |`state_properties`|Properties of `State`. Including state creation using the natural variables of the Helmholtz energy (no density iteration).| |`state_creation`|Different constructors of `State` and `PhaseEquilibrium` including critical point calculations. For pure substances and mixtures.| diff --git a/crates/feos/src/gc_pcsaft/dft/mod.rs b/crates/feos/src/gc_pcsaft/dft/mod.rs index 5005fa315..18ba4c793 100644 --- a/crates/feos/src/gc_pcsaft/dft/mod.rs +++ b/crates/feos/src/gc_pcsaft/dft/mod.rs @@ -1,6 +1,6 @@ use super::eos::GcPcSaftOptions; use super::record::GcPcSaftAssociationRecord; -use crate::association::{Association, YuWuAssociationFunctional, AssociationStrength}; +use crate::association::{Association, AssociationStrength, YuWuAssociationFunctional}; use crate::gc_pcsaft::GcPcSaftParameters; use crate::hard_sphere::{FMTContribution, FMTVersion, HardSphereProperties, MonomerShape}; use feos_core::{FeosResult, Molarweight, ResidualDyn, StateHD, Subset}; @@ -108,7 +108,8 @@ impl HelmholtzEnergyFunctionalDyn for GcPcSaftFunctional { fn contributions<'a>(&'a self) -> impl Iterator> { let mut contributions = Vec::with_capacity(4); - let assoc = YuWuAssociationFunctional::new(&self.params, &self.parameters, self.association); + let assoc = + YuWuAssociationFunctional::new(&self.params, &self.parameters, self.association); // Hard sphere contribution let hs = FMTContribution::new(&self.params, self.fmt_version); diff --git a/crates/feos/src/pcsaft/dft/mod.rs b/crates/feos/src/pcsaft/dft/mod.rs index a7f6b4566..b3afa5435 100644 --- a/crates/feos/src/pcsaft/dft/mod.rs +++ b/crates/feos/src/pcsaft/dft/mod.rs @@ -103,7 +103,8 @@ impl HelmholtzEnergyFunctionalDyn for PcSaftFunctional { fn contributions<'a>(&'a self) -> impl Iterator> { let mut contributions = Vec::with_capacity(4); - let assoc = YuWuAssociationFunctional::new(&self.params, &self.parameters, self.association); + let assoc = + YuWuAssociationFunctional::new(&self.params, &self.parameters, self.association); if matches!( self.fmt_version, diff --git a/crates/feos/src/pcsaft/eos/mod.rs b/crates/feos/src/pcsaft/eos/mod.rs index 37a885248..6a68f67f7 100644 --- a/crates/feos/src/pcsaft/eos/mod.rs +++ b/crates/feos/src/pcsaft/eos/mod.rs @@ -598,12 +598,25 @@ mod tests_parameter_fit { use super::*; use approx::assert_relative_eq; use feos_core::DensityInitialization::Liquid; - use feos_core::{Contributions, PropertiesAD, ReferenceSystem, SolverOptions}; - use feos_core::{FeosResult, ParametersAD, PhaseEquilibrium, State}; - use nalgebra::{U1, U3, U8, vector}; + use feos_core::ad::{ + BoilingTemperature, BubblePointPressure, DewPointPressure, EquilibriumLiquidDensity, + LiquidDensity, PropertyAD, VaporPressure, + }; + use feos_core::{Contributions, ReferenceSystem, SolverOptions}; + use feos_core::{FeosResult, PhaseEquilibrium, State, ad::ParametersAD}; + use nalgebra::{U1, U3, U8}; use num_dual::{Dual64, DualStruct, DualVec, partial}; use quantity::{BAR, KELVIN, LITER, MOL, PASCAL}; + fn flat_binary_params(b: &PcSaftBinary) -> Vec { + b.0.0 + .iter() + .flatten() + .copied() + .chain(std::iter::once(b.0.1)) + .collect() + } + fn pcsaft_non_assoc() -> PcSaftPure { let m = 1.5; let sigma = 3.4; @@ -626,9 +639,9 @@ mod tests_parameter_fit { "nb", ]; let (pcsaft, _) = pcsaft()?; - let pcsaft_ad = pcsaft.named_derivatives(pcsaft_params); + let pcsaft_ad = PcSaftPure::::seed_derivatives(&pcsaft.0, pcsaft_params); let temperature = 250.0 * KELVIN; - let p = pcsaft_ad.vapor_pressure(temperature)?; + let p = VaporPressure(temperature).evaluate(&pcsaft_ad)?; let p = p.convert_into(PASCAL); let (p, grad) = (p.re, p.eps.unwrap_generic(U8, U1)); @@ -658,9 +671,10 @@ mod tests_parameter_fit { #[test] fn test_vapor_pressure_derivatives_fit() -> FeosResult<()> { let pcsaft = pcsaft_non_assoc(); - let pcsaft_ad = pcsaft.named_derivatives(["m", "sigma", "epsilon_k"]); + let pcsaft_ad = + PcSaftPure::::seed_derivatives(&pcsaft.0, ["m", "sigma", "epsilon_k"]); let temperature = 150.0 * KELVIN; - let p = pcsaft_ad.vapor_pressure(temperature)?; + let p = VaporPressure(temperature).evaluate(&pcsaft_ad)?; let p = p.convert_into(PASCAL); let (p, grad) = (p.re, p.eps.unwrap_generic(U3, U1)); @@ -690,9 +704,10 @@ mod tests_parameter_fit { #[test] fn test_boiling_temperature_derivatives_fit() -> FeosResult<()> { let pcsaft = pcsaft_non_assoc(); - let pcsaft_ad = pcsaft.named_derivatives(["m", "sigma", "epsilon_k"]); + let pcsaft_ad = + PcSaftPure::::seed_derivatives(&pcsaft.0, ["m", "sigma", "epsilon_k"]); let pressure = BAR; - let t = pcsaft_ad.boiling_temperature(pressure)?; + let t = BoilingTemperature(pressure).evaluate(&pcsaft_ad)?; let t = t.convert_into(KELVIN); let (t, grad) = (t.re, t.eps.unwrap_generic(U3, U1)); @@ -734,16 +749,14 @@ mod tests_parameter_fit { #[test] fn test_equilibrium_liquid_density_derivatives_fit() -> FeosResult<()> { let pcsaft = pcsaft_non_assoc(); - let pcsaft_ad = pcsaft.named_derivatives(["m", "sigma", "epsilon_k"]); + let pcsaft_ad = + PcSaftPure::::seed_derivatives(&pcsaft.0, ["m", "sigma", "epsilon_k"]); let temperature = 150.0 * KELVIN; - let (p, rho) = pcsaft_ad.equilibrium_liquid_density(temperature)?; - let p = p.convert_into(PASCAL); + let rho = EquilibriumLiquidDensity(temperature).evaluate(&pcsaft_ad)?; let rho = rho.convert_into(MOL / LITER); - let (p, p_grad) = (p.re, p.eps.unwrap_generic(U3, U1)); let (rho, rho_grad) = (rho.re, rho.eps.unwrap_generic(U3, U1)); - println!("{p:.5} {rho:.5}"); - println!("{p_grad:.5?}"); + println!("{rho:.5}"); println!("{rho_grad:.5?}"); for (i, par) in ["m", "sigma", "epsilon_k"].into_iter().enumerate() { @@ -751,22 +764,17 @@ mod tests_parameter_fit { let h = params[i] * 1e-7; params[i] += h; let pcsaft_h = PcSaftPure(params); - let (p_h, [_, rho_h]) = + let (_, [_, rho_h]) = PhaseEquilibrium::pure_t(&pcsaft_h, temperature, None, Default::default())?; - let dp_h = (p_h.convert_into(PASCAL) - p) / h; let drho_h = (rho_h.convert_into(MOL / LITER) - rho) / h; - let dp = p_grad[i]; let drho = rho_grad[i]; println!( - "{par:12}: {:11.5} {:11.5} {:.3e} {:11.5} {:11.5} {:.3e}", - dp_h, - dp, - ((dp_h - dp) / dp).abs(), + "{par:12}: {:11.5} {:11.5} {:.3e}", drho_h, drho, ((drho_h - drho) / drho).abs() ); - assert_relative_eq!(dp, dp_h, max_relative = 1e-6); + assert_relative_eq!(drho, drho_h, max_relative = 1e-6); } Ok(()) } @@ -774,10 +782,11 @@ mod tests_parameter_fit { #[test] fn test_liquid_density_derivatives_fit() -> FeosResult<()> { let pcsaft = pcsaft_non_assoc(); - let pcsaft_ad = pcsaft.named_derivatives(["m", "sigma", "epsilon_k"]); + let pcsaft_ad = + PcSaftPure::::seed_derivatives(&pcsaft.0, ["m", "sigma", "epsilon_k"]); let temperature = 150.0 * KELVIN; let pressure = BAR; - let rho = pcsaft_ad.liquid_density(temperature, pressure)?; + let rho = LiquidDensity(temperature, pressure).evaluate(&pcsaft_ad)?; let rho = rho.convert_into(MOL / LITER); let (rho, grad) = (rho.re, rho.eps.unwrap_generic(U3, U1)); @@ -806,10 +815,11 @@ mod tests_parameter_fit { #[test] fn test_bubble_point_pressure() -> FeosResult<()> { let (pcsaft, _) = pcsaft_binary()?; - let pcsaft_ad = pcsaft.named_derivatives(["k_ij"]); + let pcsaft_ad = + PcSaftBinary::::seed_derivatives(&flat_binary_params(&pcsaft), ["k_ij"]); let temperature = 500.0 * KELVIN; - let x = vector![0.5, 0.5]; - let p = pcsaft_ad.bubble_point_pressure(temperature, None, x)?; + let x = 0.5; + let p = BubblePointPressure(temperature, x, None).evaluate(&pcsaft_ad)?; let p = p.convert_into(BAR); let (p, [[grad]]) = (p.re, p.eps.unwrap_generic(U1, U1).data.0); @@ -844,10 +854,11 @@ mod tests_parameter_fit { #[test] fn test_dew_point_pressure() -> FeosResult<()> { let (pcsaft, _) = pcsaft_binary()?; - let pcsaft_ad = pcsaft.named_derivatives(["k_ij"]); + let pcsaft_ad = + PcSaftBinary::::seed_derivatives(&flat_binary_params(&pcsaft), ["k_ij"]); let temperature = 500.0 * KELVIN; let y = 0.5; - let p = pcsaft_ad.dew_point_pressure(temperature, None, y)?; + let p = DewPointPressure(temperature, y, None).evaluate(&pcsaft_ad)?; let p = p.convert_into(BAR); let (p, [[grad]]) = (p.re, p.eps.unwrap_generic(U1, U1).data.0); @@ -876,7 +887,8 @@ mod tests_parameter_fit { #[test] fn test_bubble_point_temperature() -> FeosResult<()> { let (pcsaft, _) = pcsaft_binary()?; - let pcsaft_ad = pcsaft.named_derivatives(["k_ij"]); + let pcsaft_ad = + PcSaftBinary::::seed_derivatives(&flat_binary_params(&pcsaft), ["k_ij"]); let pressure = Pressure::from_reduced(DualVec::from(45. * BAR.into_reduced())); let t_init = Temperature::from_reduced(DualVec::from(500.0)); let x = DualVec::from(0.5); @@ -924,7 +936,8 @@ mod tests_parameter_fit { #[test] fn test_dew_point_temperature() -> FeosResult<()> { let (pcsaft, _) = pcsaft_binary()?; - let pcsaft_ad = pcsaft.named_derivatives(["k_ij"]); + let pcsaft_ad = + PcSaftBinary::::seed_derivatives(&flat_binary_params(&pcsaft), ["k_ij"]); let pressure = Pressure::from_reduced(DualVec::from(45. * BAR.into_reduced())); let t_init = Temperature::from_reduced(DualVec::from(500.0)); let x = DualVec::from(0.5); @@ -1060,7 +1073,11 @@ mod tests_parameter_fit { #[test] fn test_tp_flash() -> FeosResult<()> { let (pcsaft, _) = pcsaft_binary()?; - let pcsaft_ad = pcsaft.named_derivatives(["k_ij"]); + let (params, mut kij) = pcsaft.0; + let mut flat_params: Vec = params[0].to_vec(); + flat_params.extend_from_slice(¶ms[1]); + flat_params.push(kij); + let pcsaft_ad = PcSaftBinary::::seed_derivatives(&flat_params, ["k_ij"]); let temperature = 500.0 * KELVIN; let pressure = 44.6 * BAR; let x = 0.5; @@ -1080,8 +1097,6 @@ mod tests_parameter_fit { println!("{beta:.5}"); println!("{grad:.5?}"); - - let (params, mut kij) = pcsaft.0; let h = 1e-7; kij += h; let pcsaft_h = PcSaftBinary::new(params, kij); diff --git a/crates/feos/src/pcsaft/eos/pcsaft_binary.rs b/crates/feos/src/pcsaft/eos/pcsaft_binary.rs index 71b15ae88..160b3bd27 100644 --- a/crates/feos/src/pcsaft/eos/pcsaft_binary.rs +++ b/crates/feos/src/pcsaft/eos/pcsaft_binary.rs @@ -1,8 +1,8 @@ use super::dispersion::{A0, A1, A2, B0, B1, B2}; use super::polar::{AD, BD, CD}; -use feos_core::{ParametersAD, Residual, StateHD}; +use feos_core::{Residual, StateHD, ad::ParametersAD}; use nalgebra::{SVector, U2}; -use num_dual::{DualNum, DualSVec64, DualVec, jacobian}; +use num_dual::{DualNum, DualVec, jacobian}; use std::f64::consts::{FRAC_PI_6, PI}; const PI_SQ_43: f64 = 4.0 / 3.0 * PI * PI; @@ -19,45 +19,59 @@ impl PcSaftBinary { } } -impl + Copy, const N: usize> From<&[f64]> for PcSaftBinary { - fn from(parameters: &[f64]) -> Self { - if parameters.len() != 2 * N + 1 { - panic!( - "This version of PC-SAFT requires exactly {} parameters!", - 2 * N + 1 - ) - } - let (Ok(p1), Ok(p2)): (Result<[f64; N], _>, Result<[f64; N], _>) = - (parameters[..N].try_into(), parameters[N..2 * N].try_into()) - else { - unreachable!() - }; - let kij = D::from(parameters[2 * N]); - Self::new([p1.map(D::from), p2.map(D::from)], kij) - } -} - -impl ParametersAD<2> for PcSaftBinary { - fn index_parameters_mut<'a, const P: usize>( - eos: &'a mut Self::Lifted>, - index: &str, - ) -> &'a mut DualSVec64

{ - match index { - "k_ij" => &mut eos.0.1, - _ => panic!("{index} is not a valid binary PC-SAFT parameter!"), - } +impl ParametersAD for PcSaftBinary { + fn build + Copy>( + mut f: impl FnMut(&'static str, bool) -> D, + ) -> PcSaftBinary { + PcSaftBinary::new( + [ + [ + f("m1", true), + f("sigma1", true), + f("epsilon_k1", true), + f("mu1", true), + ], + [ + f("m2", true), + f("sigma2", true), + f("epsilon_k2", true), + f("mu2", true), + ], + ], + f("k_ij", true), + ) } } -impl ParametersAD<2> for PcSaftBinary { - fn index_parameters_mut<'a, const P: usize>( - eos: &'a mut Self::Lifted>, - index: &str, - ) -> &'a mut DualSVec64

{ - match index { - "k_ij" => &mut eos.0.1, - _ => panic!("{index} is not a valid binary PC-SAFT parameter!"), - } +impl ParametersAD for PcSaftBinary { + fn build + Copy>( + mut f: impl FnMut(&'static str, bool) -> D, + ) -> PcSaftBinary { + PcSaftBinary::new( + [ + [ + f("m1", true), + f("sigma1", true), + f("epsilon_k1", true), + f("mu1", true), + f("kappa_ab1", true), + f("epsilon_k_ab1", true), + f("na1", false), + f("nb1", false), + ], + [ + f("m2", true), + f("sigma2", true), + f("epsilon_k2", true), + f("mu2", true), + f("kappa_ab2", true), + f("epsilon_k_ab2", true), + f("na2", false), + f("nb2", false), + ], + ], + f("k_ij", true), + ) } } diff --git a/crates/feos/src/pcsaft/eos/pcsaft_pure.rs b/crates/feos/src/pcsaft/eos/pcsaft_pure.rs index f9ea2983f..46155f4cf 100644 --- a/crates/feos/src/pcsaft/eos/pcsaft_pure.rs +++ b/crates/feos/src/pcsaft/eos/pcsaft_pure.rs @@ -1,8 +1,8 @@ use super::dispersion::{A0, A1, A2, B0, B1, B2}; use super::polar::{AD, BD, CD}; -use feos_core::{ParametersAD, Residual, StateHD}; +use feos_core::{Residual, StateHD, ad::ParametersAD}; use nalgebra::{SVector, U1}; -use num_dual::{DualNum, DualSVec64}; +use num_dual::DualNum; use std::f64::consts::{FRAC_PI_6, PI}; const PI_SQ_43: f64 = 4.0 / 3.0 * PI * PI; @@ -183,46 +183,33 @@ impl + Copy> Residual for PcSaftPure { } } -impl + Copy, const N: usize> From<&[f64]> for PcSaftPure { - fn from(parameters: &[f64]) -> Self { - let Ok(parameters): Result<[f64; N], _> = parameters.try_into() else { - panic!("This version of PC-SAFT requires exactly {N} parameters!") - }; - Self(parameters.map(D::from)) +impl ParametersAD for PcSaftPure { + fn build + Copy>( + mut f: impl FnMut(&'static str, bool) -> D, + ) -> PcSaftPure { + PcSaftPure([ + f("m", true), + f("sigma", true), + f("epsilon_k", true), + f("mu", true), + ]) } } -impl ParametersAD<1> for PcSaftPure { - fn index_parameters_mut<'a, const P: usize>( - eos: &'a mut Self::Lifted>, - index: &str, - ) -> &'a mut DualSVec64

{ - match index { - "m" => &mut eos.0[0], - "sigma" => &mut eos.0[1], - "epsilon_k" => &mut eos.0[2], - "mu" => &mut eos.0[3], - _ => panic!("{index} is not a valid PC-SAFT parameter!"), - } - } -} - -impl ParametersAD<1> for PcSaftPure { - fn index_parameters_mut<'a, const P: usize>( - eos: &'a mut Self::Lifted>, - index: &str, - ) -> &'a mut DualSVec64

{ - match index { - "m" => &mut eos.0[0], - "sigma" => &mut eos.0[1], - "epsilon_k" => &mut eos.0[2], - "mu" => &mut eos.0[3], - "kappa_ab" => &mut eos.0[4], - "epsilon_k_ab" => &mut eos.0[5], - "na" => &mut eos.0[6], - "nb" => &mut eos.0[7], - _ => panic!("{index} is not a valid PC-SAFT parameter!"), - } +impl ParametersAD for PcSaftPure { + fn build + Copy>( + mut f: impl FnMut(&'static str, bool) -> D, + ) -> PcSaftPure { + PcSaftPure([ + f("m", true), + f("sigma", true), + f("epsilon_k", true), + f("mu", true), + f("kappa_ab", true), + f("epsilon_k_ab", true), + f("na", false), + f("nb", false), + ]) } } diff --git a/crates/feos/tests/pcsaft/px_flashes.rs b/crates/feos/tests/pcsaft/px_flashes.rs index 4cfb179b6..3f102a5e5 100644 --- a/crates/feos/tests/pcsaft/px_flashes.rs +++ b/crates/feos/tests/pcsaft/px_flashes.rs @@ -2,8 +2,8 @@ use approx::assert_relative_eq; use feos::ideal_gas::Joback; use feos::pcsaft::PcSaftBinary; use feos_core::{ - Contributions, EquationOfState, FeosResult, IdealGasAD, ParametersAD, PhaseEquilibrium, - ReferenceSystem, SolverOptions, Verbosity, + Contributions, EquationOfState, FeosResult, IdealGasAD, PhaseEquilibrium, ReferenceSystem, + SolverOptions, Verbosity, ad::ParametersAD, }; use nalgebra::U1; use num_dual::{DualStruct, DualVec}; @@ -38,7 +38,10 @@ fn test_ph_flash() -> FeosResult<()> { println!("{h}\n{}", vle.molar_enthalpy()); assert_relative_eq!(h, vle.molar_enthalpy(), max_relative = 1e-10); - let pcsaft_ad = pcsaft.named_derivatives(["k_ij"]); + let mut flat_params: Vec = params[0].to_vec(); + flat_params.extend_from_slice(¶ms[1]); + flat_params.push(kij); + let pcsaft_ad = PcSaftBinary::::seed_derivatives(&flat_params, ["k_ij"]); let joback_ad = joback.each_ref().map(|j| j.lift()); let eos_ad = EquationOfState::new(joback_ad, pcsaft_ad); let vle_ad = PhaseEquilibrium::ph_flash( @@ -105,7 +108,10 @@ fn test_ps_flash() -> FeosResult<()> { println!("{s}\n{}", vle.molar_entropy()); assert_relative_eq!(s, vle.molar_entropy(), max_relative = 1e-10); - let pcsaft_ad = pcsaft.named_derivatives(["k_ij"]); + let mut flat_params: Vec = params[0].to_vec(); + flat_params.extend_from_slice(¶ms[1]); + flat_params.push(kij); + let pcsaft_ad = PcSaftBinary::::seed_derivatives(&flat_params, ["k_ij"]); let joback_ad = joback.each_ref().map(|j| j.lift()); let eos_ad = EquationOfState::new(joback_ad, pcsaft_ad); let vle_ad = PhaseEquilibrium::ps_flash( diff --git a/py-feos/src/ad/dataset.rs b/py-feos/src/ad/dataset.rs new file mode 100644 index 000000000..27022c58f --- /dev/null +++ b/py-feos/src/ad/dataset.rs @@ -0,0 +1,551 @@ +use feos_core::ad::{ + BinaryDataset, BinaryProperty, BubblePointRecord, Dataset, DewPointRecord, + EnthalpyOfVaporizationRecord, EquilibriumLiquidDensityRecord, LiquidDensityRecord, PureDataset, + PureProperty, ResidualIsobaricHeatCapacityRecord, VaporPressureRecord, +}; +use ndarray::{Array2, ArrayView1}; +use numpy::{PyArray1, PyArray2, PyReadonlyArray1, ToPyArray}; +use pyo3::exceptions::{PyTypeError, PyValueError}; +use pyo3::prelude::*; + +use crate::eos::PyEquationOfState; + +/// Run a `Dataset::evaluate` against a single equation of state or a list of them. +/// +/// If `eos` extracts as a single `PyEquationOfState`, returns the +/// `(predicted, converged)` arrays as 1D. If it extracts as a sequence, +/// returns them stacked as 2D arrays with shape `[n_points, n_eos]`. +fn evaluate_eos<'py, D: Dataset>( + py: Python<'py>, + dataset: &D, + eos: &Bound<'py, PyAny>, +) -> PyResult<(Bound<'py, PyAny>, Bound<'py, PyAny>)> { + // Evaluate a single EoS. + if let Ok(e) = eos.extract::>() { + let (pred, status) = dataset.evaluate(&e.0); + return Ok(( + pred.to_pyarray(py).into_any(), + status.to_pyarray(py).into_any(), + )); + } + + // Construct a references to multiple EoSs. + let eos_refs: Vec> = eos.extract().map_err(|_| { + PyTypeError::new_err( + "expected an EquationOfState or a sequence of EquationOfState instances", + ) + })?; + + let n_points = dataset.target().len(); + let n_eos = eos_refs.len(); + let mut pred = Array2::::from_elem((n_points, n_eos), f64::NAN); + let mut status = Array2::::from_elem((n_points, n_eos), false); + + // Iterate through models > evaluate dataset > collect results as columns. + // Column index is EoS-index. + for (j, e) in eos_refs.iter().enumerate() { + let (p, c) = dataset.evaluate(&e.0); + pred.column_mut(j).assign(&p); + status.column_mut(j).assign(&c); + } + Ok(( + pred.to_pyarray(py).into_any(), + status.to_pyarray(py).into_any(), + )) +} + +fn ensure_same_len(arrays: &[(&str, usize)]) -> PyResult { + let Some((_, n)) = arrays.first() else { + return Ok(0); + }; + if arrays.iter().any(|(_, len)| len != n) { + let names = arrays + .iter() + .map(|(name, _)| *name) + .collect::>() + .join(", "); + return Err(PyValueError::new_err(format!( + "all arrays must have the same length: {names}" + ))); + } + Ok(*n) +} + +fn collect_records_2( + a: (&str, ArrayView1<'_, f64>), + b: (&str, ArrayView1<'_, f64>), + f: impl Fn(f64, f64) -> R, +) -> PyResult> { + ensure_same_len(&[(a.0, a.1.len()), (b.0, b.1.len())])?; + Ok(a.1.iter().zip(b.1.iter()).map(|(&a, &b)| f(a, b)).collect()) +} + +fn collect_records_3( + a: (&str, ArrayView1<'_, f64>), + b: (&str, ArrayView1<'_, f64>), + c: (&str, ArrayView1<'_, f64>), + f: impl Fn(f64, f64, f64) -> R, +) -> PyResult> { + ensure_same_len(&[(a.0, a.1.len()), (b.0, b.1.len()), (c.0, c.1.len())])?; + Ok(a.1 + .iter() + .zip(b.1.iter()) + .zip(c.1.iter()) + .map(|((&a, &b), &c)| f(a, b, c)) + .collect()) +} + +// Parser methods: map from str to enum +// Preferred here to declutter enums/structs exported in feos. + +fn parse_pure_property(property: &str) -> PyResult { + match property { + "vapor_pressure" => Ok(PureProperty::VaporPressure), + "liquid_density" => Ok(PureProperty::LiquidDensity), + "equilibrium_liquid_density" => Ok(PureProperty::EquilibriumLiquidDensity), + "enthalpy_of_vaporization" => Ok(PureProperty::EnthalpyOfVaporization), + "residual_isobaric_heat_capacity" => Ok(PureProperty::ResidualIsobaricHeatCapacity), + _ => Err(PyValueError::new_err(format!( + "unknown pure property '{property}'; valid: \ + 'vapor_pressure', 'liquid_density', 'equilibrium_liquid_density', \ + 'enthalpy_of_vaporization', 'residual_isobaric_heat_capacity'" + ))), + } +} + +fn parse_binary_property(property: &str) -> PyResult { + match property { + "bubble_point_pressure" | "bubble_point" => Ok(BinaryProperty::BubblePointPressure), + "dew_point_pressure" | "dew_point" => Ok(BinaryProperty::DewPointPressure), + _ => Err(PyValueError::new_err(format!( + "unknown binary property '{property}'; valid: \ + 'bubble_point_pressure', 'dew_point_pressure'" + ))), + } +} + +#[pyclass(name = "PureDataset")] +pub struct PyPureDataset { + pub(crate) inner: PureDataset, +} + +#[pymethods] +impl PyPureDataset { + /// Load pure-component data from CSV. + /// + /// Args: + /// path (str): Path to the CSV file. + /// property (str): Property identifier. Valid values are + /// ``"vapor_pressure"``, ``"liquid_density"``, + /// ``"equilibrium_liquid_density"``, ``"enthalpy_of_vaporization"``, + /// and ``"residual_isobaric_heat_capacity"``. + /// name (str, optional): Dataset name used in regressor diagnostics. + #[staticmethod] + #[pyo3(signature = (path, property, name=None))] + pub fn from_csv(path: &str, property: &str, name: Option<&str>) -> PyResult { + Self::from_csv_for_property(parse_pure_property(property)?, path, name) + } + + /// Load vapor pressure data from CSV. + /// + /// CSV columns: ``temperature_k, vapor_pressure_pa``. + #[staticmethod] + #[pyo3(signature = (path, name=None))] + pub fn vapor_pressure_from_csv(path: &str, name: Option<&str>) -> PyResult { + Self::from_csv_for_property(PureProperty::VaporPressure, path, name) + } + + /// Construct vapor pressure data from numpy arrays. + #[staticmethod] + #[pyo3(signature = (temperature_k, vapor_pressure_pa, name=None))] + pub fn vapor_pressure( + temperature_k: PyReadonlyArray1, + vapor_pressure_pa: PyReadonlyArray1, + name: Option<&str>, + ) -> PyResult { + let temperature_k = temperature_k.as_array(); + let vapor_pressure_pa = vapor_pressure_pa.as_array(); + let records = collect_records_2( + ("temperature_k", temperature_k), + ("vapor_pressure_pa", vapor_pressure_pa), + |temperature_k, vapor_pressure_pa| VaporPressureRecord { + temperature_k, + vapor_pressure_pa, + }, + )?; + Ok(Self::with_optional_name( + PureDataset::vapor_pressure(records), + name, + )) + } + + /// Load liquid density data from CSV. + /// + /// CSV columns: ``temperature_k, pressure_pa, liquid_density_kmol_m3``. + #[staticmethod] + #[pyo3(signature = (path, name=None))] + pub fn liquid_density_from_csv(path: &str, name: Option<&str>) -> PyResult { + Self::from_csv_for_property(PureProperty::LiquidDensity, path, name) + } + + /// Construct liquid density data from numpy arrays. + #[staticmethod] + #[pyo3(signature = (temperature_k, pressure_pa, liquid_density_kmol_m3, name=None))] + pub fn liquid_density( + temperature_k: PyReadonlyArray1, + pressure_pa: PyReadonlyArray1, + liquid_density_kmol_m3: PyReadonlyArray1, + name: Option<&str>, + ) -> PyResult { + let temperature_k = temperature_k.as_array(); + let pressure_pa = pressure_pa.as_array(); + let liquid_density_kmol_m3 = liquid_density_kmol_m3.as_array(); + let records = collect_records_3( + ("temperature_k", temperature_k), + ("pressure_pa", pressure_pa), + ("liquid_density_kmol_m3", liquid_density_kmol_m3), + |temperature_k, pressure_pa, liquid_density_kmol_m3| LiquidDensityRecord { + temperature_k, + pressure_pa, + liquid_density_kmol_m3, + }, + )?; + Ok(Self::with_optional_name( + PureDataset::liquid_density(records), + name, + )) + } + + /// Load saturated liquid density data from CSV. + /// + /// CSV columns: ``temperature_k, liquid_density_kmol_m3``. + #[staticmethod] + #[pyo3(signature = (path, name=None))] + pub fn equilibrium_liquid_density_from_csv(path: &str, name: Option<&str>) -> PyResult { + Self::from_csv_for_property(PureProperty::EquilibriumLiquidDensity, path, name) + } + + /// Construct saturated liquid density data from numpy arrays. + #[staticmethod] + #[pyo3(signature = (temperature_k, liquid_density_kmol_m3, name=None))] + pub fn equilibrium_liquid_density( + temperature_k: PyReadonlyArray1, + liquid_density_kmol_m3: PyReadonlyArray1, + name: Option<&str>, + ) -> PyResult { + let temperature_k = temperature_k.as_array(); + let liquid_density_kmol_m3 = liquid_density_kmol_m3.as_array(); + let records = collect_records_2( + ("temperature_k", temperature_k), + ("liquid_density_kmol_m3", liquid_density_kmol_m3), + |temperature_k, liquid_density_kmol_m3| EquilibriumLiquidDensityRecord { + temperature_k, + liquid_density_kmol_m3, + }, + )?; + Ok(Self::with_optional_name( + PureDataset::equilibrium_liquid_density(records), + name, + )) + } + + /// Load enthalpy of vaporization data from CSV. + /// + /// CSV columns: ``temperature_k, dh_vap_j_mol``. + #[staticmethod] + #[pyo3(signature = (path, name=None))] + pub fn enthalpy_of_vaporization_from_csv(path: &str, name: Option<&str>) -> PyResult { + Self::from_csv_for_property(PureProperty::EnthalpyOfVaporization, path, name) + } + + /// Construct enthalpy of vaporization data from numpy arrays. + #[staticmethod] + #[pyo3(signature = (temperature_k, dh_vap_j_mol, name=None))] + pub fn enthalpy_of_vaporization( + temperature_k: PyReadonlyArray1, + dh_vap_j_mol: PyReadonlyArray1, + name: Option<&str>, + ) -> PyResult { + let temperature_k = temperature_k.as_array(); + let dh_vap_j_mol = dh_vap_j_mol.as_array(); + let records = collect_records_2( + ("temperature_k", temperature_k), + ("dh_vap_j_mol", dh_vap_j_mol), + |temperature_k, dh_vap_j_mol| EnthalpyOfVaporizationRecord { + temperature_k, + dh_vap_j_mol, + }, + )?; + Ok(Self::with_optional_name( + PureDataset::enthalpy_of_vaporization(records), + name, + )) + } + + /// Load residual isobaric heat capacity data from CSV. + /// + /// CSV columns: ``temperature_k, pressure_pa, cp_res_j_molk``. + #[staticmethod] + #[pyo3(signature = (path, name=None))] + pub fn residual_isobaric_heat_capacity_from_csv( + path: &str, + name: Option<&str>, + ) -> PyResult { + Self::from_csv_for_property(PureProperty::ResidualIsobaricHeatCapacity, path, name) + } + + /// Construct residual isobaric heat capacity data from numpy arrays. + #[staticmethod] + #[pyo3(signature = (temperature_k, pressure_pa, cp_res_j_molk, name=None))] + pub fn residual_isobaric_heat_capacity( + temperature_k: PyReadonlyArray1, + pressure_pa: PyReadonlyArray1, + cp_res_j_molk: PyReadonlyArray1, + name: Option<&str>, + ) -> PyResult { + let temperature_k = temperature_k.as_array(); + let pressure_pa = pressure_pa.as_array(); + let cp_res_j_molk = cp_res_j_molk.as_array(); + let records = collect_records_3( + ("temperature_k", temperature_k), + ("pressure_pa", pressure_pa), + ("cp_res_j_molk", cp_res_j_molk), + |temperature_k, pressure_pa, cp_res_j_molk| ResidualIsobaricHeatCapacityRecord { + temperature_k, + pressure_pa, + cp_res_j_molk, + }, + )?; + Ok(Self::with_optional_name( + PureDataset::residual_isobaric_heat_capacity(records), + name, + )) + } + + /// Property name. + #[getter] + pub fn name(&self) -> &str { + self.inner.name() + } + + /// Number of data points. + pub fn __len__(&self) -> usize { + self.inner.target().len() + } + + /// Target values. + pub fn target<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray1> { + self.inner.target().to_owned().to_pyarray(py) + } + + /// Input values. + pub fn inputs<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray2> { + self.inner.inputs().to_owned().to_pyarray(py) + } + + /// Evaluate the dataset's property for one or more equations of state. + /// + /// Args: + /// eos: A single ``EquationOfState`` or a list of them. + /// Each must describe a single substance. + /// + /// Returns: + /// ``(predicted, converged)``. For a single EoS, both are 1D arrays + /// of length ``n_points``. For a list of ``n_eos`` EoS, both are + /// 2D arrays of shape ``[n_points, n_eos]``; column ``k`` + /// corresponds to ``eos[k]``. Non-converged points are reported as + /// ``NaN`` in ``predicted`` and ``False`` in ``converged``. + pub fn evaluate<'py>( + &self, + py: Python<'py>, + eos: &Bound<'py, PyAny>, + ) -> PyResult<(Bound<'py, PyAny>, Bound<'py, PyAny>)> { + evaluate_eos(py, &self.inner, eos) + } + + pub fn __repr__(&self) -> String { + format!( + "PureDataset(property={}, n={})", + self.inner.name(), + self.inner.target().len() + ) + } +} + +impl PyPureDataset { + fn from_csv_for_property( + property: PureProperty, + path: &str, + name: Option<&str>, + ) -> PyResult { + PureDataset::from_csv(property, std::path::Path::new(path)) + .map(|inner| Self::with_optional_name(inner, name)) + .map_err(|e| PyValueError::new_err(e.to_string())) + } + + fn with_optional_name(mut inner: PureDataset, name: Option<&str>) -> Self { + if let Some(n) = name { + inner = inner.with_name(n); + } + Self { inner } + } +} + +#[pyclass(name = "BinaryDataset")] +pub struct PyBinaryDataset { + pub(crate) inner: BinaryDataset, +} + +#[pymethods] +impl PyBinaryDataset { + /// Load binary-mixture data from CSV. + /// + /// Args: + /// path (str): Path to the CSV file. + /// property (str): Property identifier. Valid values are + /// ``"bubble_point_pressure"`` and ``"dew_point_pressure"``. + /// Short aliases ``"bubble_point"`` and ``"dew_point"`` are also accepted. + /// name (str, optional): Dataset name used in regressor diagnostics. + #[staticmethod] + #[pyo3(signature = (path, property, name=None))] + pub fn from_csv(path: &str, property: &str, name: Option<&str>) -> PyResult { + Self::from_csv_for_property(parse_binary_property(property)?, path, name) + } + + /// Load bubble point pressure data from CSV. + /// + /// CSV columns: ``temperature_k, liquid_molefrac_1, bubble_pressure_pa``. + #[staticmethod] + #[pyo3(signature = (path, name=None))] + pub fn bubble_point_pressure_from_csv(path: &str, name: Option<&str>) -> PyResult { + Self::from_csv_for_property(BinaryProperty::BubblePointPressure, path, name) + } + + /// Construct bubble point pressure data from numpy arrays. + #[staticmethod] + #[pyo3(signature = (temperature_k, liquid_molefrac_1, bubble_pressure_pa, name=None))] + pub fn bubble_point_pressure( + temperature_k: PyReadonlyArray1, + liquid_molefrac_1: PyReadonlyArray1, + bubble_pressure_pa: PyReadonlyArray1, + name: Option<&str>, + ) -> PyResult { + let temperature_k = temperature_k.as_array(); + let liquid_molefrac_1 = liquid_molefrac_1.as_array(); + let bubble_pressure_pa = bubble_pressure_pa.as_array(); + let records = collect_records_3( + ("temperature_k", temperature_k), + ("liquid_molefrac_1", liquid_molefrac_1), + ("bubble_pressure_pa", bubble_pressure_pa), + |temperature_k, liquid_molefrac_1, bubble_pressure_pa| BubblePointRecord { + temperature_k, + liquid_molefrac_1, + bubble_pressure_pa, + }, + )?; + Ok(Self::with_optional_name( + BinaryDataset::bubble_point_pressure(records), + name, + )) + } + + /// Load dew point pressure data from CSV. + /// + /// CSV columns: ``temperature_k, vapor_molefrac_1, dew_pressure_pa``. + #[staticmethod] + #[pyo3(signature = (path, name=None))] + pub fn dew_point_pressure_from_csv(path: &str, name: Option<&str>) -> PyResult { + Self::from_csv_for_property(BinaryProperty::DewPointPressure, path, name) + } + + /// Construct dew point pressure data from numpy arrays. + #[staticmethod] + #[pyo3(signature = (temperature_k, vapor_molefrac_1, dew_pressure_pa, name=None))] + pub fn dew_point_pressure( + temperature_k: PyReadonlyArray1, + vapor_molefrac_1: PyReadonlyArray1, + dew_pressure_pa: PyReadonlyArray1, + name: Option<&str>, + ) -> PyResult { + let temperature_k = temperature_k.as_array(); + let vapor_molefrac_1 = vapor_molefrac_1.as_array(); + let dew_pressure_pa = dew_pressure_pa.as_array(); + let records = collect_records_3( + ("temperature_k", temperature_k), + ("vapor_molefrac_1", vapor_molefrac_1), + ("dew_pressure_pa", dew_pressure_pa), + |temperature_k, vapor_molefrac_1, dew_pressure_pa| DewPointRecord { + temperature_k, + vapor_molefrac_1, + dew_pressure_pa, + }, + )?; + Ok(Self::with_optional_name( + BinaryDataset::dew_point_pressure(records), + name, + )) + } + + /// Property name. + #[getter] + pub fn name(&self) -> &str { + self.inner.name() + } + + /// Number of data points. + pub fn __len__(&self) -> usize { + self.inner.target().len() + } + + /// Target values. + pub fn target<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray1> { + self.inner.target().to_owned().to_pyarray(py) + } + + /// Evaluate the dataset's property for one or more equations of state (no gradients). + /// + /// Args: + /// eos: A single ``EquationOfState`` or a list of them. Each + /// must describe a binary system (``components() == 2``). + /// + /// Returns: + /// ``(predicted, converged)``. For a single EoS, both are 1D arrays + /// of length ``n_points``. For a list of ``n_eos`` EoS, both are + /// 2D arrays of shape ``[n_points, n_eos]``; column ``k`` + /// corresponds to ``eos[k]``. Non-converged points are reported as + /// ``NaN`` in ``predicted`` and ``False`` in ``converged``. + pub fn evaluate<'py>( + &self, + py: Python<'py>, + eos: &Bound<'py, PyAny>, + ) -> PyResult<(Bound<'py, PyAny>, Bound<'py, PyAny>)> { + evaluate_eos(py, &self.inner, eos) + } + + pub fn __repr__(&self) -> String { + format!( + "BinaryDataset(property={}, n={})", + self.inner.name(), + self.inner.target().len() + ) + } +} + +impl PyBinaryDataset { + fn from_csv_for_property( + property: BinaryProperty, + path: &str, + name: Option<&str>, + ) -> PyResult { + BinaryDataset::from_csv(property, std::path::Path::new(path)) + .map(|inner| Self::with_optional_name(inner, name)) + .map_err(|e| PyValueError::new_err(e.to_string())) + } + + fn with_optional_name(mut inner: BinaryDataset, name: Option<&str>) -> Self { + if let Some(n) = name { + inner = inner.with_name(n); + } + Self { inner } + } +} diff --git a/py-feos/src/ad/mod.rs b/py-feos/src/ad/mod.rs index 3f775a2da..71eac67a0 100644 --- a/py-feos/src/ad/mod.rs +++ b/py-feos/src/ad/mod.rs @@ -1,9 +1,19 @@ +use crate::eos::PyEquationOfState; use feos::pcsaft::{PcSaftBinary, PcSaftPure}; -use feos_core::{ParametersAD, PropertiesAD}; -use numpy::{PyArray1, PyArray2, PyReadonlyArray2, ToPyArray}; +use feos_core::ad::{ + BoilingTemperature, BubblePointPressure, DewPointPressure, EnthalpyOfVaporization, + EquilibriumLiquidDensity, LiquidDensity, ParametersAD, PropertyAD, + ResidualIsobaricHeatCapacity, VaporPressure, +}; +use nalgebra::{U1, U2}; +use num_dual::DualSVec; +use numpy::{PyArray1, PyArray2, PyReadonlyArray1, PyReadonlyArray2, ToPyArray}; use paste::paste; use pyo3::prelude::*; +pub mod dataset; +pub use dataset::{PyBinaryDataset, PyPureDataset}; + #[pyclass(name = "EquationOfStateAD", eq, eq_int)] #[derive(Clone, Copy, PartialEq)] pub enum PyEquationOfStateAD { @@ -25,168 +35,411 @@ impl From for BinaryModels { } } +type EvalResult<'py> = (Bound<'py, PyArray1>, Bound<'py, PyArray1>); + type GradResult<'py> = ( Bound<'py, PyArray1>, Bound<'py, PyArray2>, Bound<'py, PyArray1>, ); -/// Calculate vapor pressures and derivatives w.r.t. model parameters. -/// -/// Parameters -/// ---------- -/// model: EquationOfStateAD -/// The equation of state to use. -/// parameter_names: List[string] -/// The name of the parameters for which derivatives are calculated. -/// parameters: np.ndarray[float] -/// The parameters for every data point. -/// input: np.ndarray[float] -/// The temperature (in K) for every data point. -/// -/// Returns -/// ------- -/// (np.ndarray[float], np.ndarray[float], np.ndarray[bool]): The vapor pressures (in Pa), gradients, and convergence status. -#[pyfunction] -pub fn vapor_pressure_derivatives<'py>( - model: PyEquationOfStateAD, - parameter_names: &Bound<'py, PyAny>, - parameters: PyReadonlyArray2, - input: PyReadonlyArray2, -) -> GradResult<'py> { - _vapor_pressure_derivatives(model, parameter_names, parameters, input) -} +#[pyclass(name = "Property")] +pub struct PyPropertyAD; -/// Calculate boiling temperatures and derivatives w.r.t. model parameters. -/// -/// Parameters -/// ---------- -/// model: EquationOfStateAD -/// The equation of state to use. -/// parameter_names: List[string] -/// The name of the parameters for which derivatives are calculated. -/// parameters: np.ndarray[float] -/// The parameters for every data point. -/// input: np.ndarray[float] -/// The pressure (in Pa) for every data point. -/// -/// Returns -/// ------- -/// (np.ndarray[float], np.ndarray[float], np.ndarray[bool]): The boiling temperature (in K), gradients, and convergence status. -#[pyfunction] -pub fn boiling_temperature_derivatives<'py>( - model: PyEquationOfStateAD, - parameter_names: &Bound<'py, PyAny>, - parameters: PyReadonlyArray2, - input: PyReadonlyArray2, -) -> GradResult<'py> { - _boiling_temperature_derivatives(model, parameter_names, parameters, input) -} +#[pymethods] +impl PyPropertyAD { + /// Calculate vapor pressures in parallel. + /// + /// Parameters + /// ---------- + /// eos: EquationOfState + /// The equation of state to use. + /// input: np.ndarray[float] + /// The temperature (in K) for every data point. + /// + /// Returns + /// ------- + /// (np.ndarray[float], np.ndarray[bool]): The vapor pressures (in Pa), and convergence status. + #[staticmethod] + pub fn vapor_pressure<'py>( + py: Python<'py>, + eos: &PyEquationOfState, + input: PyReadonlyArray2, + ) -> EvalResult<'py> { + let (value, status) = VaporPressure::evaluate_parallel(&eos.0, input.as_array()); + (value.to_pyarray(py), status.to_pyarray(py)) + } -/// Calculate liquid densities and derivatives w.r.t. model parameters. -/// -/// Parameters -/// ---------- -/// model: EquationOfStateAD -/// The equation of state to use. -/// parameter_names: List[string] -/// The name of the parameters for which derivatives are calculated. -/// parameters: np.ndarray[float] -/// The parameters for every data point. -/// input: np.ndarray[float] -/// The temperature (in K) and pressure (in Pa) for every data point. -/// -/// Returns -/// ------- -/// (np.ndarray[float], np.ndarray[float], np.ndarray[bool]): The liquid densities (in kmol/m³), gradients, and convergence status. -#[pyfunction] -pub fn liquid_density_derivatives<'py>( - model: PyEquationOfStateAD, - parameter_names: &Bound<'py, PyAny>, - parameters: PyReadonlyArray2, - input: PyReadonlyArray2, -) -> GradResult<'py> { - _liquid_density_derivatives(model, parameter_names, parameters, input) -} + /// Calculate vapor pressures and derivatives w.r.t. model parameters. + /// + /// Parameters + /// ---------- + /// model: EquationOfStateAD + /// The equation of state to use. + /// parameter_names: List[string] + /// The name of the parameters for which derivatives are calculated. + /// parameters: np.ndarray[float] + /// The parameters for every data point. + /// input: np.ndarray[float] + /// The temperature (in K) for every data point. + /// + /// Returns + /// ------- + /// (np.ndarray[float], np.ndarray[float], np.ndarray[bool]): The vapor pressures (in Pa), gradients, and convergence status. + #[staticmethod] + pub fn vapor_pressure_derivatives<'py>( + model: PyEquationOfStateAD, + parameter_names: &Bound<'py, PyAny>, + parameters: &Bound<'py, PyAny>, + input: PyReadonlyArray2, + ) -> GradResult<'py> { + _vapor_pressure_derivatives(model, parameter_names, parameters, input) + } -/// Calculate liquid densities at saturation and derivatives w.r.t. model parameters. -/// -/// Parameters -/// ---------- -/// model: EquationOfStateAD -/// The equation of state to use. -/// parameter_names: List[string] -/// The name of the parameters for which derivatives are calculated. -/// parameters: np.ndarray[float] -/// The parameters for every data point. -/// input: np.ndarray[float] -/// The temperature (in K) for every data point. -/// -/// Returns -/// ------- -/// (np.ndarray[float], np.ndarray[float], np.ndarray[bool]): The liquid densities (in kmol/m³), gradients, and convergence status. -#[pyfunction] -pub fn equilibrium_liquid_density_derivatives<'py>( - model: PyEquationOfStateAD, - parameter_names: &Bound<'py, PyAny>, - parameters: PyReadonlyArray2, - input: PyReadonlyArray2, -) -> GradResult<'py> { - _equilibrium_liquid_density_derivatives(model, parameter_names, parameters, input) -} + /// Calculate boiling temperatures in parallel. + /// + /// Parameters + /// ---------- + /// eos: EquationOfState + /// The equation of state to use. + /// input: np.ndarray[float] + /// The pressure (in Pa) for every data point. + /// + /// Returns + /// ------- + /// (np.ndarray[float], np.ndarray[bool]): The boiling temperature (in K), and convergence status. + #[staticmethod] + pub fn boiling_temperature<'py>( + py: Python<'py>, + eos: &PyEquationOfState, + input: PyReadonlyArray2, + ) -> EvalResult<'py> { + let (value, status) = BoilingTemperature::evaluate_parallel(&eos.0, input.as_array()); + (value.to_pyarray(py), status.to_pyarray(py)) + } -/// Calculate bubble point pressures of binary mixtures and derivatives w.r.t. model parameters. -/// -/// Parameters -/// ---------- -/// model: EquationOfStateAD -/// The equation of state to use. -/// parameter_names: List[string] -/// The name of the parameters for which derivatives are calculated. -/// parameters: np.ndarray[float] -/// The parameters for every data point. -/// input: np.ndarray[float] -/// The temperature (in K), composition of the first component, and an initial guess for the -/// pressure (in Pa) for every data point. -/// -/// Returns -/// ------- -/// (np.ndarray[float], np.ndarray[float], np.ndarray[bool]): The bubble point pressures (in Pa), gradients, and convergence status. -#[pyfunction] -pub fn bubble_point_pressure_derivatives<'py>( - model: PyEquationOfStateAD, - parameter_names: &Bound<'py, PyAny>, - parameters: PyReadonlyArray2, - input: PyReadonlyArray2, -) -> GradResult<'py> { - _bubble_point_pressure_derivatives(model, parameter_names, parameters, input) -} + /// Calculate boiling temperatures and derivatives w.r.t. model parameters. + /// + /// Parameters + /// ---------- + /// model: EquationOfStateAD + /// The equation of state to use. + /// parameter_names: List[string] + /// The name of the parameters for which derivatives are calculated. + /// parameters: np.ndarray[float] + /// The parameters for every data point. + /// input: np.ndarray[float] + /// The pressure (in Pa) for every data point. + /// + /// Returns + /// ------- + /// (np.ndarray[float], np.ndarray[float], np.ndarray[bool]): The boiling temperature (in K), gradients, and convergence status. + #[staticmethod] + pub fn boiling_temperature_derivatives<'py>( + model: PyEquationOfStateAD, + parameter_names: &Bound<'py, PyAny>, + parameters: &Bound<'py, PyAny>, + input: PyReadonlyArray2, + ) -> GradResult<'py> { + _boiling_temperature_derivatives(model, parameter_names, parameters, input) + } + + /// Calculate liquid densities in parallel. + /// + /// Parameters + /// ---------- + /// eos: EquationOfState + /// The equation of state to use. + /// input: np.ndarray[float] + /// The temperature (in K) and pressure (in Pa) for every data point. + /// + /// Returns + /// ------- + /// (np.ndarray[float], np.ndarray[bool]): The liquid densities (in kmol/m³), and convergence status. + #[staticmethod] + pub fn liquid_density<'py>( + py: Python<'py>, + eos: &PyEquationOfState, + input: PyReadonlyArray2, + ) -> EvalResult<'py> { + let (value, status) = LiquidDensity::evaluate_parallel(&eos.0, input.as_array()); + (value.to_pyarray(py), status.to_pyarray(py)) + } + + /// Calculate liquid densities and derivatives w.r.t. model parameters. + /// + /// Parameters + /// ---------- + /// model: EquationOfStateAD + /// The equation of state to use. + /// parameter_names: List[string] + /// The name of the parameters for which derivatives are calculated. + /// parameters: np.ndarray[float] + /// The parameters for every data point. + /// input: np.ndarray[float] + /// The temperature (in K) and pressure (in Pa) for every data point. + /// + /// Returns + /// ------- + /// (np.ndarray[float], np.ndarray[float], np.ndarray[bool]): The liquid densities (in kmol/m³), gradients, and convergence status. + #[staticmethod] + pub fn liquid_density_derivatives<'py>( + model: PyEquationOfStateAD, + parameter_names: &Bound<'py, PyAny>, + parameters: &Bound<'py, PyAny>, + input: PyReadonlyArray2, + ) -> GradResult<'py> { + _liquid_density_derivatives(model, parameter_names, parameters, input) + } -/// Calculate dew point pressures of binary mixtures and derivatives w.r.t. model parameters. -/// -/// Parameters -/// ---------- -/// model: EquationOfStateAD -/// The equation of state to use. -/// parameter_names: List[string] -/// The name of the parameters for which derivatives are calculated. -/// parameters: np.ndarray[float] -/// The parameters for every data point. -/// input: np.ndarray[float] -/// The temperature (in K), composition of the first component, and an initial guess for the -/// pressure (in Pa) for every data point. -/// -/// Returns -/// ------- -/// (np.ndarray[float], np.ndarray[float], np.ndarray[bool]): The dew point pressures (in Pa), gradients, and convergence status. -#[pyfunction] -pub fn dew_point_pressure_derivatives<'py>( - model: PyEquationOfStateAD, - parameter_names: &Bound<'py, PyAny>, - parameters: PyReadonlyArray2, - input: PyReadonlyArray2, -) -> GradResult<'py> { - _dew_point_pressure_derivatives(model, parameter_names, parameters, input) + /// Calculate liquid densities at saturation in parallel. + /// + /// Parameters + /// ---------- + /// eos: EquationOfState + /// The equation of state to use. + /// input: np.ndarray[float] + /// The temperature (in K) for every data point. + /// + /// Returns + /// ------- + /// (np.ndarray[float], np.ndarray[bool]): The liquid densities (in kmol/m³), and convergence status. + #[staticmethod] + pub fn equilibrium_liquid_density<'py>( + py: Python<'py>, + eos: &PyEquationOfState, + input: PyReadonlyArray2, + ) -> EvalResult<'py> { + let (value, status) = EquilibriumLiquidDensity::evaluate_parallel(&eos.0, input.as_array()); + (value.to_pyarray(py), status.to_pyarray(py)) + } + + /// Calculate liquid densities at saturation and derivatives w.r.t. model parameters. + /// + /// Parameters + /// ---------- + /// model: EquationOfStateAD + /// The equation of state to use. + /// parameter_names: List[string] + /// The name of the parameters for which derivatives are calculated. + /// parameters: np.ndarray[float] + /// The parameters for every data point. + /// input: np.ndarray[float] + /// The temperature (in K) for every data point. + /// + /// Returns + /// ------- + /// (np.ndarray[float], np.ndarray[float], np.ndarray[bool]): The liquid densities (in kmol/m³), gradients, and convergence status. + #[staticmethod] + pub fn equilibrium_liquid_density_derivatives<'py>( + model: PyEquationOfStateAD, + parameter_names: &Bound<'py, PyAny>, + parameters: &Bound<'py, PyAny>, + input: PyReadonlyArray2, + ) -> GradResult<'py> { + _equilibrium_liquid_density_derivatives(model, parameter_names, parameters, input) + } + + /// Calculate enthalpy of vaporization in parallel. + /// + /// Parameters + /// ---------- + /// eos: EquationOfState + /// The equation of state to use. + /// input: np.ndarray[float] + /// The temperature (in K) for every data point. + /// + /// Returns + /// ------- + /// (np.ndarray[float], np.ndarray[bool]): + /// The enthalpies of vaporization (in J/mol), and convergence status. + #[staticmethod] + pub fn enthalpy_of_vaporization<'py>( + py: Python<'py>, + eos: &PyEquationOfState, + input: PyReadonlyArray2, + ) -> EvalResult<'py> { + let (value, status) = EnthalpyOfVaporization::evaluate_parallel(&eos.0, input.as_array()); + (value.to_pyarray(py), status.to_pyarray(py)) + } + + /// Calculate enthalpy of vaporization and derivatives w.r.t. model parameters. + /// + /// Parameters + /// ---------- + /// model: EquationOfStateAD + /// The equation of state to use. + /// parameter_names: List[string] + /// The name of the parameters for which derivatives are calculated. + /// parameters: np.ndarray[float] + /// The parameters for every data point. + /// input: np.ndarray[float] + /// The temperature (in K) for every data point. + /// + /// Returns + /// ------- + /// (np.ndarray[float], np.ndarray[float], np.ndarray[bool]): + /// The enthalpies of vaporization (in J/mol), gradients, and convergence status. + #[staticmethod] + pub fn enthalpy_of_vaporization_derivatives<'py>( + model: PyEquationOfStateAD, + parameter_names: &Bound<'py, PyAny>, + parameters: &Bound<'py, PyAny>, + input: PyReadonlyArray2, + ) -> GradResult<'py> { + _enthalpy_of_vaporization_derivatives(model, parameter_names, parameters, input) + } + + /// Calculate residual isobaric molar heat capacities (liquid phase) in parallel. + /// + /// Parameters + /// ---------- + /// eos: EquationOfState + /// The equation of state to use. + /// input: np.ndarray[float] + /// The temperature (in K) and pressure (in Pa) for every data point. + /// + /// Returns + /// ------- + /// (np.ndarray[float], np.ndarray[bool]): + /// The residual isobaric heat capacities (in J/(mol·K)), and convergence status. + #[staticmethod] + pub fn residual_isobaric_heat_capacity<'py>( + py: Python<'py>, + eos: &PyEquationOfState, + input: PyReadonlyArray2, + ) -> EvalResult<'py> { + let (value, status) = + ResidualIsobaricHeatCapacity::evaluate_parallel(&eos.0, input.as_array()); + (value.to_pyarray(py), status.to_pyarray(py)) + } + + /// Calculate residual isobaric molar heat capacities (liquid phase) and derivatives w.r.t. model parameters. + /// + /// Parameters + /// ---------- + /// model: EquationOfStateAD + /// The equation of state to use. + /// parameter_names: List[string] + /// The name of the parameters for which derivatives are calculated. + /// parameters: np.ndarray[float] + /// The parameters for every data point. + /// input: np.ndarray[float] + /// The temperature (in K) and pressure (in Pa) for every data point. + /// + /// Returns + /// ------- + /// (np.ndarray[float], np.ndarray[float], np.ndarray[bool]): + /// The residual isobaric heat capacities (in J/(mol·K)), gradients, and convergence status. + #[staticmethod] + pub fn residual_isobaric_heat_capacity_derivatives<'py>( + model: PyEquationOfStateAD, + parameter_names: &Bound<'py, PyAny>, + parameters: &Bound<'py, PyAny>, + input: PyReadonlyArray2, + ) -> GradResult<'py> { + _residual_isobaric_heat_capacity_derivatives(model, parameter_names, parameters, input) + } + + /// Calculate bubble point pressures of binary mixtures in parallel. + /// + /// Parameters + /// ---------- + /// eos: EquationOfState + /// The equation of state to use. + /// input: np.ndarray[float] + /// The temperature (in K), composition of the first component, and an initial guess for the + /// pressure (in Pa) for every data point. + /// + /// Returns + /// ------- + /// (np.ndarray[float], np.ndarray[bool]): The bubble point pressures (in Pa), and convergence status. + #[staticmethod] + pub fn bubble_point_pressure<'py>( + py: Python<'py>, + eos: &PyEquationOfState, + input: PyReadonlyArray2, + ) -> EvalResult<'py> { + let (value, status) = BubblePointPressure::evaluate_parallel(&eos.0, input.as_array()); + (value.to_pyarray(py), status.to_pyarray(py)) + } + + /// Calculate bubble point pressures of binary mixtures and derivatives w.r.t. model parameters. + /// + /// Parameters + /// ---------- + /// model: EquationOfStateAD + /// The equation of state to use. + /// parameter_names: List[string] + /// The name of the parameters for which derivatives are calculated. + /// parameters: np.ndarray[float] + /// The parameters for every data point. + /// input: np.ndarray[float] + /// The temperature (in K), composition of the first component, and an initial guess for the + /// pressure (in Pa) for every data point. + /// + /// Returns + /// ------- + /// (np.ndarray[float], np.ndarray[float], np.ndarray[bool]): The bubble point pressures (in Pa), gradients, and convergence status. + #[staticmethod] + pub fn bubble_point_pressure_derivatives<'py>( + model: PyEquationOfStateAD, + parameter_names: &Bound<'py, PyAny>, + parameters: &Bound<'py, PyAny>, + input: PyReadonlyArray2, + ) -> GradResult<'py> { + _bubble_point_pressure_derivatives(model, parameter_names, parameters, input) + } + + /// Calculate dew point pressures of binary mixtures in parallel. + /// + /// Parameters + /// ---------- + /// eos: EquationOfState + /// The equation of state to use. + /// input: np.ndarray[float] + /// The temperature (in K), composition of the first component, and an initial guess for the + /// pressure (in Pa) for every data point. + /// + /// Returns + /// ------- + /// (np.ndarray[float], np.ndarray[bool]): The dew point pressures (in Pa), and convergence status. + #[staticmethod] + pub fn dew_point_pressure<'py>( + py: Python<'py>, + eos: &PyEquationOfState, + input: PyReadonlyArray2, + ) -> EvalResult<'py> { + let (value, status) = DewPointPressure::evaluate_parallel(&eos.0, input.as_array()); + (value.to_pyarray(py), status.to_pyarray(py)) + } + + /// Calculate dew point pressures of binary mixtures and derivatives w.r.t. model parameters. + /// + /// Parameters + /// ---------- + /// model: EquationOfStateAD + /// The equation of state to use. + /// parameter_names: List[string] + /// The name of the parameters for which derivatives are calculated. + /// parameters: np.ndarray[float] + /// The parameters for every data point. + /// input: np.ndarray[float] + /// The temperature (in K), composition of the first component, and an initial guess for the + /// pressure (in Pa) for every data point. + /// + /// Returns + /// ------- + /// (np.ndarray[float], np.ndarray[float], np.ndarray[bool]): The dew point pressures (in Pa), gradients, and convergence status. + #[staticmethod] + pub fn dew_point_pressure_derivatives<'py>( + model: PyEquationOfStateAD, + parameter_names: &Bound<'py, PyAny>, + parameters: &Bound<'py, PyAny>, + input: PyReadonlyArray2, + ) -> GradResult<'py> { + _dew_point_pressure_derivatives(model, parameter_names, parameters, input) + } } macro_rules! expand_models { @@ -196,7 +449,7 @@ macro_rules! expand_models { fn [<_ $prop _derivatives>]<'py>( model: PyEquationOfStateAD, parameter_names: &Bound<'py, PyAny>, - parameters: PyReadonlyArray2, + parameters: &Bound<'py, PyAny>, input: PyReadonlyArray2, ) -> GradResult<'py> { match <$enum>::from(model) { @@ -210,50 +463,67 @@ macro_rules! expand_models { } macro_rules! impl_evaluate_gradients { - (pure, [$($prop:ident),*], $models:tt) => { - $(impl_evaluate_gradients!(1,PyEquationOfStateAD,$prop,$models,0,1,2,3,4,5,max:6);)* + (pure, [$($prop:ident: $prop_type:ty),*], $models:tt) => { + $(impl_evaluate_gradients!(U1,PyEquationOfStateAD,$prop,$prop_type,$models,0,1,2,3,4,5,max:6);)* }; - (binary, [$($prop:ident),*], $models:tt) => { - $(impl_evaluate_gradients!(2,BinaryModels,$prop,$models,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,max:15);)* + (binary, [$($prop:ident: $prop_type:ty),*], $models:tt) => { + $(impl_evaluate_gradients!(U2,BinaryModels,$prop,$prop_type,$models,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,max:15);)* }; - ($n:literal, $enum:ty, $prop:ident, {$($model:ident: $type:ty),*}, $($p:literal,)* max: $max:literal) => { + ($n:ty, $enum:ty, $prop:ident, $prop_type:ty, {$($model:ident: $type:ty),*}, $($p:literal,)* max: $max:literal) => { expand_models!($enum, $prop, $($model: $type),*); - paste!( fn $prop<'py, R: ParametersAD<$n>>( parameter_names: &Bound<'py, PyAny>, - parameters: PyReadonlyArray2, + parameters: &Bound<'py, PyAny>, input: PyReadonlyArray2, ) -> ( Bound<'py, PyArray1>, Bound<'py, PyArray2>, Bound<'py, PyArray1>, - ) { + ) + where + $(R::Lifted>: Sync,)* + R::Lifted>: Sync + { let (value, grad, status) = - $( - if let Ok(p) = parameter_names.extract::<[String; $p]>() { - R::[<$prop _parallel>](p, parameters.as_array(), input.as_array()) - } else)* if let Ok(p) = parameter_names.extract::<[String; $max]>() { - R::[<$prop _parallel>](p, parameters.as_array(), input.as_array()) + if let Ok(pars) = parameters.extract::>() { + let pars = pars.as_slice().expect("Parameter array is not contiguous!"); + $( + if let Ok(p) = parameter_names.extract::<[String; $p]>() { + <$prop_type>::evaluate_parallel_derivatives::(p, pars, input.as_array()) + } else)* if let Ok(p) = parameter_names.extract::<[String; $max]>() { + <$prop_type>::evaluate_parallel_derivatives::(p, pars, input.as_array()) + } else { + panic!("Gradients can only be evaluated for up to {} parameters!", $max) + } + } else if let Ok(pars) = parameters.extract::>() { + $( + if let Ok(p) = parameter_names.extract::<[String; $p]>() { + <$prop_type>::evaluate_parallel_derivatives_params::(p, pars.as_array(), input.as_array()) + } else)* if let Ok(p) = parameter_names.extract::<[String; $max]>() { + <$prop_type>::evaluate_parallel_derivatives_params::(p, pars.as_array(), input.as_array()) + } else { + panic!("Gradients can only be evaluated for up to {} parameters!", $max) + } } else { - panic!("Gradients can only be evaluated for up to {} parameters!", $max) + panic!("Argument `parameters` needs to be a 1D or 2D array!") }; ( value.to_pyarray(parameter_names.py()), grad.to_pyarray(parameter_names.py()), status.to_pyarray(parameter_names.py()), ) - }); + } }; } impl_evaluate_gradients!( pure, - [vapor_pressure, boiling_temperature, liquid_density, equilibrium_liquid_density], + [vapor_pressure: VaporPressure, boiling_temperature: BoilingTemperature, liquid_density: LiquidDensity, equilibrium_liquid_density: EquilibriumLiquidDensity, enthalpy_of_vaporization: EnthalpyOfVaporization, residual_isobaric_heat_capacity: ResidualIsobaricHeatCapacity], {PcSaftNonAssoc: PcSaftPure, PcSaftFull: PcSaftPure} ); impl_evaluate_gradients!( binary, - [bubble_point_pressure, dew_point_pressure], + [bubble_point_pressure: BubblePointPressure, dew_point_pressure: DewPointPressure], {PcSaftNonAssoc: PcSaftBinary, PcSaftFull: PcSaftBinary} ); diff --git a/py-feos/src/lib.rs b/py-feos/src/lib.rs index 948a24d6b..017dbb4a6 100644 --- a/py-feos/src/lib.rs +++ b/py-feos/src/lib.rs @@ -179,25 +179,15 @@ fn feos(m: &Bound<'_, PyModule>) -> PyResult<()> { // Equation of state m.add_class::()?; - // // Estimator - // m.add_class::()?; - // m.add_class::()?; - // m.add_class::()?; - // m.add_class::()?; - // AD #[cfg(feature = "ad")] { - m.add_function(wrap_pyfunction!(ad::vapor_pressure_derivatives, m)?)?; - m.add_function(wrap_pyfunction!(ad::boiling_temperature_derivatives, m)?)?; - m.add_function(wrap_pyfunction!(ad::liquid_density_derivatives, m)?)?; - m.add_function(wrap_pyfunction!( - ad::equilibrium_liquid_density_derivatives, - m - )?)?; - m.add_function(wrap_pyfunction!(ad::bubble_point_pressure_derivatives, m)?)?; - m.add_function(wrap_pyfunction!(ad::dew_point_pressure_derivatives, m)?)?; m.add_class::()?; + m.add_class::()?; + + // Datasets + m.add_class::()?; + m.add_class::()?; } #[cfg(feature = "dft")] diff --git a/py-feos/src/user_defined.rs b/py-feos/src/user_defined.rs index 38bbe38ca..641238f0b 100644 --- a/py-feos/src/user_defined.rs +++ b/py-feos/src/user_defined.rs @@ -14,7 +14,10 @@ impl PyIdealGas { pub fn new(obj: Bound<'_, PyAny>) -> PyResult { let attr = obj.hasattr("ln_lambda3")?; if !attr { - panic!("{}", "Python Class has to have a method 'ln_lambda3' with signature:\n\tdef ln_lambda3(self, temperature: HD) -> HD\nwhere 'HD' has to be any (hyper-) dual number.") + panic!( + "{}", + "Python Class has to have a method 'ln_lambda3' with signature:\n\tdef ln_lambda3(self, temperature: HD) -> HD\nwhere 'HD' has to be any (hyper-) dual number." + ) } Ok(Self(obj.unbind())) } @@ -57,23 +60,34 @@ impl PyResidual { pub fn new(obj: Bound<'_, PyAny>) -> PyResult { let attr = obj.hasattr("components")?; if !attr { - panic!("Python Class has to have a method 'components' with signature:\n\tdef signature(self) -> int") + panic!( + "Python Class has to have a method 'components' with signature:\n\tdef signature(self) -> int" + ) } let attr = obj.hasattr("subset")?; if !attr { - panic!("Python Class has to have a method 'subset' with signature:\n\tdef subset(self, component_list: List[int]) -> Self") + panic!( + "Python Class has to have a method 'subset' with signature:\n\tdef subset(self, component_list: List[int]) -> Self" + ) } let attr = obj.hasattr("molar_weight")?; if !attr { - panic!("Python Class has to have a method 'molar_weight' with signature:\n\tdef molar_weight(self) -> SIArray1\nwhere the size of the returned array has to be 'components'.") + panic!( + "Python Class has to have a method 'molar_weight' with signature:\n\tdef molar_weight(self) -> SIArray1\nwhere the size of the returned array has to be 'components'." + ) } let attr = obj.hasattr("max_density")?; if !attr { - panic!("Python Class has to have a method 'max_density' with signature:\n\tdef max_density(self, moles: numpy.ndarray[float]) -> float\nwhere the size of the input array has to be 'components'.") + panic!( + "Python Class has to have a method 'max_density' with signature:\n\tdef max_density(self, moles: numpy.ndarray[float]) -> float\nwhere the size of the input array has to be 'components'." + ) } let attr = obj.hasattr("helmholtz_energy")?; if !attr { - panic!("{}", "Python Class has to have a method 'helmholtz_energy' with signature:\n\tdef helmholtz_energy(self, state: StateHD) -> HD\nwhere 'HD' has to be any of {{float, Dual64, HyperDual64, HyperDualDual64, Dual3Dual64, Dual3_64}}.") + panic!( + "{}", + "Python Class has to have a method 'helmholtz_energy' with signature:\n\tdef helmholtz_energy(self, state: StateHD) -> HD\nwhere 'HD' has to be any of {{float, Dual64, HyperDual64, HyperDualDual64, Dual3Dual64, Dual3_64}}." + ) } Ok(Self(obj.unbind())) }