From 8b24b9a1c61b17d84711805d84a20df8ecfa6f77 Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Mon, 19 Dec 2022 18:37:58 +0100 Subject: [PATCH 1/9] add second (non-mixed) derivative to cache and property evaluation --- Cargo.toml | 4 +++ benches/dual_numbers.rs | 2 +- feos-core/src/equation_of_state.rs | 6 +++- feos-core/src/python/statehd.rs | 16 ++++++++-- feos-core/src/python/user_defined.rs | 7 +++-- feos-core/src/state/cache.rs | 34 +++++++++++++++++--- feos-core/src/state/mod.rs | 18 +++++++++-- feos-core/src/state/properties.rs | 42 ++++++++++++++++++------- feos-dft/src/functional_contribution.rs | 3 ++ 9 files changed, 108 insertions(+), 24 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 58f415ba0..046042203 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -52,6 +52,10 @@ criterion = "0.4" name = "state_properties" harness = false +[[bench]] +name = "state_creation" +harness = false + [[bench]] name = "dual_numbers" harness = false diff --git a/benches/dual_numbers.rs b/benches/dual_numbers.rs index bf9f37f38..a0d67a586 100644 --- a/benches/dual_numbers.rs +++ b/benches/dual_numbers.rs @@ -49,7 +49,7 @@ fn bench_dual_numbers( b.iter(|| a_res((&state.eos, &state.derive1(Derivative::DV)))) }); group.bench_function("a_hyperdual", |b| { - b.iter(|| a_res((&state.eos, &state.derive2(Derivative::DV, Derivative::DV)))) + b.iter(|| a_res((&state.eos, &state.derive2partial(Derivative::DV, Derivative::DV)))) }); group.bench_function("a_dual3", |b| { b.iter(|| a_res((&state.eos, &state.derive3(Derivative::DV)))) diff --git a/feos-core/src/equation_of_state.rs b/feos-core/src/equation_of_state.rs index e2efca5c7..406da7d0a 100644 --- a/feos-core/src/equation_of_state.rs +++ b/feos-core/src/equation_of_state.rs @@ -2,7 +2,7 @@ use crate::errors::{EosError, EosResult}; use crate::state::StateHD; use crate::EosUnit; use ndarray::prelude::*; -use num_dual::{Dual, Dual3, Dual3_64, Dual64, DualNum, DualVec64, HyperDual, HyperDual64}; +use num_dual::{Dual, Dual3, Dual3_64, Dual64, DualNum, DualVec64, HyperDual, HyperDual64, Dual2_64}; use num_traits::{One, Zero}; use quantity::{QuantityArray1, QuantityScalar}; use std::fmt; @@ -28,6 +28,7 @@ pub trait HelmholtzEnergy: + HelmholtzEnergyDual + HelmholtzEnergyDual, f64>> + HelmholtzEnergyDual + + HelmholtzEnergyDual + HelmholtzEnergyDual + HelmholtzEnergyDual> + HelmholtzEnergyDual, f64>> @@ -46,6 +47,7 @@ impl HelmholtzEnergy for T where + HelmholtzEnergyDual + HelmholtzEnergyDual, f64>> + HelmholtzEnergyDual + + HelmholtzEnergyDual + HelmholtzEnergyDual + HelmholtzEnergyDual> + HelmholtzEnergyDual, f64>> @@ -99,6 +101,7 @@ pub trait IdealGasContribution: + IdealGasContributionDual + IdealGasContributionDual, f64>> + IdealGasContributionDual + + IdealGasContributionDual + IdealGasContributionDual + IdealGasContributionDual> + IdealGasContributionDual, f64>> @@ -115,6 +118,7 @@ impl IdealGasContribution for T where + IdealGasContributionDual + IdealGasContributionDual, f64>> + IdealGasContributionDual + + IdealGasContributionDual + IdealGasContributionDual + IdealGasContributionDual> + IdealGasContributionDual, f64>> diff --git a/feos-core/src/python/statehd.rs b/feos-core/src/python/statehd.rs index ab4dc2ae6..25af15a2b 100644 --- a/feos-core/src/python/statehd.rs +++ b/feos-core/src/python/statehd.rs @@ -1,7 +1,9 @@ use crate::StateHD; use ndarray::Array1; -use num_dual::python::{PyDual3Dual64, PyDual3_64, PyDual64, PyHyperDual64, PyHyperDualDual64}; -use num_dual::{Dual, Dual3, Dual3_64, Dual64, DualVec64, HyperDual, HyperDual64}; +use num_dual::python::{ + PyDual2_64, PyDual3Dual64, PyDual3_64, PyDual64, PyHyperDual64, PyHyperDualDual64, +}; +use num_dual::{Dual, Dual2_64, Dual3, Dual3_64, Dual64, DualVec64, HyperDual, HyperDual64}; use pyo3::prelude::*; use std::convert::From; @@ -75,6 +77,16 @@ impl From, f64>>> for PyStateDDV3 { } } +#[pyclass(name = "StateD2")] +#[derive(Clone)] +pub struct PyStateD2(StateHD); + +impl From> for PyStateD2 { + fn from(s: StateHD) -> Self { + Self(s) + } +} + #[pyclass(name = "StateD3")] #[derive(Clone)] pub struct PyStateD3(StateHD); diff --git a/feos-core/src/python/user_defined.rs b/feos-core/src/python/user_defined.rs index e218ee60e..939e73fd8 100644 --- a/feos-core/src/python/user_defined.rs +++ b/feos-core/src/python/user_defined.rs @@ -1,8 +1,10 @@ use crate::python::statehd::*; use crate::*; use ndarray::prelude::*; -use num_dual::python::{PyDual3Dual64, PyDual3_64, PyDual64, PyHyperDual64, PyHyperDualDual64}; -use num_dual::{Dual, Dual3, Dual3_64, Dual64, DualVec64, HyperDual, HyperDual64}; +use num_dual::python::{ + PyDual2_64, PyDual3Dual64, PyDual3_64, PyDual64, PyHyperDual64, PyHyperDualDual64, +}; +use num_dual::{Dual, Dual2_64, Dual3, Dual3_64, Dual64, DualVec64, HyperDual, HyperDual64}; use numpy::convert::IntoPyArray; use pyo3::prelude::*; use quantity::python::PySIArray1; @@ -215,6 +217,7 @@ impl_helmholtz_energy!( PyHyperDualDualVec64_3, HyperDual, f64> ); +impl_helmholtz_energy!(PyStateD2, PyDual2_64, Dual2_64); impl_helmholtz_energy!(PyStateD3, PyDual3_64, Dual3_64); impl_helmholtz_energy!(PyStateD3D, PyDual3Dual64, Dual3); impl_helmholtz_energy!(PyStateD3DV2, PyDual3DualVec64_2, Dual3, f64>); diff --git a/feos-core/src/state/cache.rs b/feos-core/src/state/cache.rs index 3723ed7e6..505a70627 100644 --- a/feos-core/src/state/cache.rs +++ b/feos-core/src/state/cache.rs @@ -50,6 +50,26 @@ impl Cache { } } + pub fn get_or_insert_with_d2_64 Dual2_64>( + &mut self, + derivative: Derivative, + f: F, + ) -> f64 { + if let Some(&value) = self.map.get(&PartialDerivative::Second(derivative)) { + self.hit += 1; + value + } else { + self.miss += 1; + let value = f(); + self.map.insert(PartialDerivative::Zeroth, value.re); + self.map + .insert(PartialDerivative::First(derivative), value.v1[0]); + self.map + .insert(PartialDerivative::Second(derivative), value.v2[0]); + value.v2[0] + } + } + pub fn get_or_insert_with_hd64 HyperDual64>( &mut self, derivative1: Derivative, @@ -58,7 +78,7 @@ impl Cache { ) -> f64 { let d1 = min(derivative1, derivative2); let d2 = max(derivative1, derivative2); - if let Some(&value) = self.map.get(&PartialDerivative::Second(d1, d2)) { + if let Some(&value) = self.map.get(&PartialDerivative::SecondPartial(d1, d2)) { self.hit += 1; value } else { @@ -69,8 +89,10 @@ impl Cache { .insert(PartialDerivative::First(derivative1), value.eps1[0]); self.map .insert(PartialDerivative::First(derivative2), value.eps2[0]); - self.map - .insert(PartialDerivative::Second(d1, d2), value.eps1eps2[(0, 0)]); + self.map.insert( + PartialDerivative::SecondPartial(d1, d2), + value.eps1eps2[(0, 0)], + ); value.eps1eps2[(0, 0)] } } @@ -90,7 +112,11 @@ impl Cache { self.map .insert(PartialDerivative::First(derivative), value.v1); self.map - .insert(PartialDerivative::Second(derivative, derivative), value.v2); + .insert(PartialDerivative::Second(derivative), value.v2); + self.map.insert( + PartialDerivative::SecondPartial(derivative, derivative), + value.v2, + ); self.map .insert(PartialDerivative::Third(derivative), value.v3); value.v3 diff --git a/feos-core/src/state/mod.rs b/feos-core/src/state/mod.rs index 13ff9792a..fcb28c985 100644 --- a/feos-core/src/state/mod.rs +++ b/feos-core/src/state/mod.rs @@ -211,7 +211,8 @@ impl Derivative { pub(crate) enum PartialDerivative { Zeroth, First(Derivative), - Second(Derivative, Derivative), + Second(Derivative), + SecondPartial(Derivative, Derivative), Third(Derivative), } @@ -693,7 +694,20 @@ impl State { } /// Creates a [StateHD] taking the first and second (partial) derivatives. - pub fn derive2( + pub fn derive2(&self, derivative: Derivative) -> StateHD { + let mut t = Dual2_64::from(self.reduced_temperature); + let mut v = Dual2_64::from(self.reduced_volume); + let mut n = self.reduced_moles.mapv(Dual2_64::from); + match derivative { + Derivative::DT => t = t.derive(), + Derivative::DV => v = v.derive(), + Derivative::DN(i) => n[i] = n[i].derive(), + } + StateHD::new(t, v, n) + } + + /// Creates a [StateHD] taking the first and second (partial) derivatives. + pub fn derive2partial( &self, derivative1: Derivative, derivative2: Derivative, diff --git a/feos-core/src/state/properties.rs b/feos-core/src/state/properties.rs index e84be6645..09a5bd1ee 100644 --- a/feos-core/src/state/properties.rs +++ b/feos-core/src/state/properties.rs @@ -50,8 +50,13 @@ impl State { -(new_state.moles.sum() * new_state.temperature * new_state.volume.ln()).eps[0] * (U::reference_energy() / v.reference()) } - PartialDerivative::Second(v1, v2) => { - let new_state = self.derive2(v1, v2); + PartialDerivative::Second(v) => { + let new_state = self.derive2(v); + -(new_state.moles.sum() * new_state.temperature * new_state.volume.ln()).v2[0] + * (U::reference_energy() / (v.reference() * v.reference())) + } + PartialDerivative::SecondPartial(v1, v2) => { + let new_state = self.derive2partial(v1, v2); -(new_state.moles.sum() * new_state.temperature * new_state.volume.ln()) .eps1eps2[(0, 0)] * (U::reference_energy() / (v1.reference() * v2.reference())) @@ -82,8 +87,15 @@ impl State { cache.get_or_insert_with_d64(v, &computation) * U::reference_energy() / v.reference() } - PartialDerivative::Second(v1, v2) => { - let new_state = self.derive2(v1, v2); + PartialDerivative::Second(v) => { + let new_state = self.derive2(v); + let computation = + || self.eos.evaluate_residual(&new_state) * new_state.temperature; + cache.get_or_insert_with_d2_64(v, &computation) * U::reference_energy() + / (v.reference() * v.reference()) + } + PartialDerivative::SecondPartial(v1, v2) => { + let new_state = self.derive2partial(v1, v2); let computation = || self.eos.evaluate_residual(&new_state) * new_state.temperature; cache.get_or_insert_with_hd64(v1, v2, &computation) * U::reference_energy() @@ -114,8 +126,14 @@ impl State { * U::reference_energy() / v.reference() } - PartialDerivative::Second(v1, v2) => { - let new_state = self.derive2(v1, v2); + PartialDerivative::Second(v) => { + let new_state = self.derive2(v); + (self.eos.ideal_gas().evaluate(&new_state) * new_state.temperature).v2[0] + * U::reference_energy() + / (v.reference() * v.reference()) + } + PartialDerivative::SecondPartial(v1, v2) => { + let new_state = self.derive2partial(v1, v2); (self.eos.ideal_gas().evaluate(&new_state) * new_state.temperature).eps1eps2 [(0, 0)] * U::reference_energy() @@ -190,16 +208,16 @@ impl State { } fn dp_dv_(&self, evaluate: Evaluate) -> QuantityScalar { - -self.get_or_compute_derivative(PartialDerivative::Second(DV, DV), evaluate) + -self.get_or_compute_derivative(PartialDerivative::Second(DV), evaluate) } fn dp_dt_(&self, evaluate: Evaluate) -> QuantityScalar { - -self.get_or_compute_derivative(PartialDerivative::Second(DV, DT), evaluate) + -self.get_or_compute_derivative(PartialDerivative::SecondPartial(DV, DT), evaluate) } fn dp_dni_(&self, evaluate: Evaluate) -> QuantityArray1 { QuantityArray::from_shape_fn(self.eos.components(), |i| { - -self.get_or_compute_derivative(PartialDerivative::Second(DV, DN(i)), evaluate) + -self.get_or_compute_derivative(PartialDerivative::SecondPartial(DV, DN(i)), evaluate) }) } @@ -209,19 +227,19 @@ impl State { fn dmu_dt_(&self, evaluate: Evaluate) -> QuantityArray1 { QuantityArray::from_shape_fn(self.eos.components(), |i| { - self.get_or_compute_derivative(PartialDerivative::Second(DT, DN(i)), evaluate) + self.get_or_compute_derivative(PartialDerivative::SecondPartial(DT, DN(i)), evaluate) }) } fn dmu_dni_(&self, evaluate: Evaluate) -> QuantityArray2 { let n = self.eos.components(); QuantityArray::from_shape_fn((n, n), |(i, j)| { - self.get_or_compute_derivative(PartialDerivative::Second(DN(i), DN(j)), evaluate) + self.get_or_compute_derivative(PartialDerivative::SecondPartial(DN(i), DN(j)), evaluate) }) } fn ds_dt_(&self, evaluate: Evaluate) -> QuantityScalar { - -self.get_or_compute_derivative(PartialDerivative::Second(DT, DT), evaluate) + -self.get_or_compute_derivative(PartialDerivative::Second(DT), evaluate) } fn d2s_dt2_(&self, evaluate: Evaluate) -> QuantityScalar { diff --git a/feos-dft/src/functional_contribution.rs b/feos-dft/src/functional_contribution.rs index cc5dc099f..b458d386b 100644 --- a/feos-dft/src/functional_contribution.rs +++ b/feos-dft/src/functional_contribution.rs @@ -36,6 +36,7 @@ impl_helmholtz_energy!(f64); impl_helmholtz_energy!(Dual64); impl_helmholtz_energy!(Dual, f64>); impl_helmholtz_energy!(HyperDual64); +impl_helmholtz_energy!(Dual2_64); impl_helmholtz_energy!(Dual3_64); impl_helmholtz_energy!(HyperDual); impl_helmholtz_energy!(HyperDual, f64>); @@ -76,6 +77,7 @@ pub trait FunctionalContribution: + FunctionalContributionDual + FunctionalContributionDual, f64>> + FunctionalContributionDual + + FunctionalContributionDual + FunctionalContributionDual + FunctionalContributionDual> + FunctionalContributionDual, f64>> @@ -161,6 +163,7 @@ impl FunctionalContribution for T where + FunctionalContributionDual + FunctionalContributionDual, f64>> + FunctionalContributionDual + + FunctionalContributionDual + FunctionalContributionDual + FunctionalContributionDual> + FunctionalContributionDual, f64>> From 7db1a5d9fd3d1136640fff7259aca91fbd3157ec Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Mon, 19 Dec 2022 20:02:27 +0100 Subject: [PATCH 2/9] renamed derivative to SecondMixed, renamed some tests --- benches/dual_numbers.rs | 2 +- benches/state_creation.rs | 58 +++++++++++++++++++++++++++++++ feos-core/src/state/cache.rs | 6 ++-- feos-core/src/state/mod.rs | 4 +-- feos-core/src/state/properties.rs | 20 +++++------ 5 files changed, 74 insertions(+), 16 deletions(-) create mode 100644 benches/state_creation.rs diff --git a/benches/dual_numbers.rs b/benches/dual_numbers.rs index a0d67a586..a96709ac7 100644 --- a/benches/dual_numbers.rs +++ b/benches/dual_numbers.rs @@ -49,7 +49,7 @@ fn bench_dual_numbers( b.iter(|| a_res((&state.eos, &state.derive1(Derivative::DV)))) }); group.bench_function("a_hyperdual", |b| { - b.iter(|| a_res((&state.eos, &state.derive2partial(Derivative::DV, Derivative::DV)))) + b.iter(|| a_res((&state.eos, &state.derive2_mixed(Derivative::DV, Derivative::DV)))) }); group.bench_function("a_dual3", |b| { b.iter(|| a_res((&state.eos, &state.derive3(Derivative::DV)))) diff --git a/benches/state_creation.rs b/benches/state_creation.rs new file mode 100644 index 000000000..fee9060b7 --- /dev/null +++ b/benches/state_creation.rs @@ -0,0 +1,58 @@ +use criterion::{criterion_group, criterion_main, Criterion}; +use feos::pcsaft::{PcSaft, PcSaftParameters}; +use feos_core::{ + parameter::{IdentifierOption, Parameter}, + DensityInitialization, PhaseEquilibrium, State, +}; +use ndarray::arr1; +use quantity::si::*; +use std::sync::Arc; + +/// Evaluate NPT constructor +fn npt( + (eos, t, p, n, rho0): ( + &Arc, + SINumber, + SINumber, + &SIArray1, + DensityInitialization, + ), +) { + State::new_npt(eos, t, p, n, rho0).unwrap(); +} + +/// Evaluate critical point constructor +fn critical_point((eos, n): (&Arc, Option<&SIArray1>)) { + State::critical_point(eos, n, None, Default::default()).unwrap(); +} + +fn tp_flash((eos, t, p, feed): (&Arc, SINumber, SINumber, &SIArray1)) { + PhaseEquilibrium::tp_flash(eos, t, p, feed, None, Default::default(), None).unwrap(); +} + +fn states_pcsaft(c: &mut Criterion) { + let parameters = PcSaftParameters::from_json( + vec!["methane", "ethane", "propane"], + "./parameters/pcsaft/gross2001.json", + None, + IdentifierOption::Name, + ) + .unwrap(); + let eos = Arc::new(PcSaft::new(Arc::new(parameters))); + let t = 300.0 * KELVIN; + let x = arr1(&[1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0]); + let n = &x * 100.0 * MOL; + let p = BAR; + + let mut group = c.benchmark_group("methane_ethane_propane"); + group.bench_function("NPT", |b| { + b.iter(|| npt((&eos, t, p, &n, DensityInitialization::None))) + }); + group.bench_function("critical_point", |b| { + b.iter(|| critical_point((&eos, Some(&n)))) + }); + group.bench_function("tp_flash", |b| b.iter(|| tp_flash((&eos, 315.0 * KELVIN, 72.0 * BAR, &n)))); +} + +criterion_group!(bench, states_pcsaft); +criterion_main!(bench); diff --git a/feos-core/src/state/cache.rs b/feos-core/src/state/cache.rs index 505a70627..5ef34a5e8 100644 --- a/feos-core/src/state/cache.rs +++ b/feos-core/src/state/cache.rs @@ -78,7 +78,7 @@ impl Cache { ) -> f64 { let d1 = min(derivative1, derivative2); let d2 = max(derivative1, derivative2); - if let Some(&value) = self.map.get(&PartialDerivative::SecondPartial(d1, d2)) { + if let Some(&value) = self.map.get(&PartialDerivative::SecondMixed(d1, d2)) { self.hit += 1; value } else { @@ -90,7 +90,7 @@ impl Cache { self.map .insert(PartialDerivative::First(derivative2), value.eps2[0]); self.map.insert( - PartialDerivative::SecondPartial(d1, d2), + PartialDerivative::SecondMixed(d1, d2), value.eps1eps2[(0, 0)], ); value.eps1eps2[(0, 0)] @@ -114,7 +114,7 @@ impl Cache { self.map .insert(PartialDerivative::Second(derivative), value.v2); self.map.insert( - PartialDerivative::SecondPartial(derivative, derivative), + PartialDerivative::SecondMixed(derivative, derivative), value.v2, ); self.map diff --git a/feos-core/src/state/mod.rs b/feos-core/src/state/mod.rs index fcb28c985..159972784 100644 --- a/feos-core/src/state/mod.rs +++ b/feos-core/src/state/mod.rs @@ -212,7 +212,7 @@ pub(crate) enum PartialDerivative { Zeroth, First(Derivative), Second(Derivative), - SecondPartial(Derivative, Derivative), + SecondMixed(Derivative, Derivative), Third(Derivative), } @@ -707,7 +707,7 @@ impl State { } /// Creates a [StateHD] taking the first and second (partial) derivatives. - pub fn derive2partial( + pub fn derive2_mixed( &self, derivative1: Derivative, derivative2: Derivative, diff --git a/feos-core/src/state/properties.rs b/feos-core/src/state/properties.rs index 09a5bd1ee..9a580926b 100644 --- a/feos-core/src/state/properties.rs +++ b/feos-core/src/state/properties.rs @@ -55,8 +55,8 @@ impl State { -(new_state.moles.sum() * new_state.temperature * new_state.volume.ln()).v2[0] * (U::reference_energy() / (v.reference() * v.reference())) } - PartialDerivative::SecondPartial(v1, v2) => { - let new_state = self.derive2partial(v1, v2); + PartialDerivative::SecondMixed(v1, v2) => { + let new_state = self.derive2_mixed(v1, v2); -(new_state.moles.sum() * new_state.temperature * new_state.volume.ln()) .eps1eps2[(0, 0)] * (U::reference_energy() / (v1.reference() * v2.reference())) @@ -94,8 +94,8 @@ impl State { cache.get_or_insert_with_d2_64(v, &computation) * U::reference_energy() / (v.reference() * v.reference()) } - PartialDerivative::SecondPartial(v1, v2) => { - let new_state = self.derive2partial(v1, v2); + PartialDerivative::SecondMixed(v1, v2) => { + let new_state = self.derive2_mixed(v1, v2); let computation = || self.eos.evaluate_residual(&new_state) * new_state.temperature; cache.get_or_insert_with_hd64(v1, v2, &computation) * U::reference_energy() @@ -132,8 +132,8 @@ impl State { * U::reference_energy() / (v.reference() * v.reference()) } - PartialDerivative::SecondPartial(v1, v2) => { - let new_state = self.derive2partial(v1, v2); + PartialDerivative::SecondMixed(v1, v2) => { + let new_state = self.derive2_mixed(v1, v2); (self.eos.ideal_gas().evaluate(&new_state) * new_state.temperature).eps1eps2 [(0, 0)] * U::reference_energy() @@ -212,12 +212,12 @@ impl State { } fn dp_dt_(&self, evaluate: Evaluate) -> QuantityScalar { - -self.get_or_compute_derivative(PartialDerivative::SecondPartial(DV, DT), evaluate) + -self.get_or_compute_derivative(PartialDerivative::SecondMixed(DV, DT), evaluate) } fn dp_dni_(&self, evaluate: Evaluate) -> QuantityArray1 { QuantityArray::from_shape_fn(self.eos.components(), |i| { - -self.get_or_compute_derivative(PartialDerivative::SecondPartial(DV, DN(i)), evaluate) + -self.get_or_compute_derivative(PartialDerivative::SecondMixed(DV, DN(i)), evaluate) }) } @@ -227,14 +227,14 @@ impl State { fn dmu_dt_(&self, evaluate: Evaluate) -> QuantityArray1 { QuantityArray::from_shape_fn(self.eos.components(), |i| { - self.get_or_compute_derivative(PartialDerivative::SecondPartial(DT, DN(i)), evaluate) + self.get_or_compute_derivative(PartialDerivative::SecondMixed(DT, DN(i)), evaluate) }) } fn dmu_dni_(&self, evaluate: Evaluate) -> QuantityArray2 { let n = self.eos.components(); QuantityArray::from_shape_fn((n, n), |(i, j)| { - self.get_or_compute_derivative(PartialDerivative::SecondPartial(DN(i), DN(j)), evaluate) + self.get_or_compute_derivative(PartialDerivative::SecondMixed(DN(i), DN(j)), evaluate) }) } From f448a02e5a22450369ed074f812c4daba3d63e3c Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Tue, 20 Dec 2022 13:55:24 +0100 Subject: [PATCH 3/9] change cache to only store second order derivatives in SecondMixed --- benches/state_creation.rs | 1 + feos-core/src/state/cache.rs | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/benches/state_creation.rs b/benches/state_creation.rs index fee9060b7..ad7e5753f 100644 --- a/benches/state_creation.rs +++ b/benches/state_creation.rs @@ -26,6 +26,7 @@ fn critical_point((eos, n): (&Arc, Option<&SIArray1>)) { State::critical_point(eos, n, None, Default::default()).unwrap(); } +/// Evaluate temperature, pressure flash. fn tp_flash((eos, t, p, feed): (&Arc, SINumber, SINumber, &SIArray1)) { PhaseEquilibrium::tp_flash(eos, t, p, feed, None, Default::default(), None).unwrap(); } diff --git a/feos-core/src/state/cache.rs b/feos-core/src/state/cache.rs index 5ef34a5e8..d7401e9d2 100644 --- a/feos-core/src/state/cache.rs +++ b/feos-core/src/state/cache.rs @@ -55,7 +55,7 @@ impl Cache { derivative: Derivative, f: F, ) -> f64 { - if let Some(&value) = self.map.get(&PartialDerivative::Second(derivative)) { + if let Some(&value) = self.map.get(&PartialDerivative::SecondMixed(derivative, derivative)) { self.hit += 1; value } else { @@ -65,7 +65,7 @@ impl Cache { self.map .insert(PartialDerivative::First(derivative), value.v1[0]); self.map - .insert(PartialDerivative::Second(derivative), value.v2[0]); + .insert(PartialDerivative::SecondMixed(derivative, derivative), value.v2[0]); value.v2[0] } } From 3667b4bb7e81d6f7b364ff35cb4ad15774fbc2d6 Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Tue, 20 Dec 2022 14:17:55 +0100 Subject: [PATCH 4/9] remove redundant hashmap entry in cache --- feos-core/src/state/cache.rs | 2 -- 1 file changed, 2 deletions(-) diff --git a/feos-core/src/state/cache.rs b/feos-core/src/state/cache.rs index d7401e9d2..daf2a3e4c 100644 --- a/feos-core/src/state/cache.rs +++ b/feos-core/src/state/cache.rs @@ -111,8 +111,6 @@ impl Cache { self.map.insert(PartialDerivative::Zeroth, value.re); self.map .insert(PartialDerivative::First(derivative), value.v1); - self.map - .insert(PartialDerivative::Second(derivative), value.v2); self.map.insert( PartialDerivative::SecondMixed(derivative, derivative), value.v2, From e1726b8c8db31aaff35767c35a3b732155a8dff4 Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Tue, 20 Dec 2022 15:29:51 +0100 Subject: [PATCH 5/9] Add more benchmarks and systems, better names for benchmark groups --- benches/dual_numbers.rs | 11 +-- benches/state_creation.rs | 140 ++++++++++++++++++++++++++++++------ benches/state_properties.rs | 2 +- 3 files changed, 128 insertions(+), 25 deletions(-) diff --git a/benches/dual_numbers.rs b/benches/dual_numbers.rs index a96709ac7..b38607b4b 100644 --- a/benches/dual_numbers.rs +++ b/benches/dual_numbers.rs @@ -48,6 +48,9 @@ fn bench_dual_numbers( group.bench_function("a_dual", |b| { b.iter(|| a_res((&state.eos, &state.derive1(Derivative::DV)))) }); + group.bench_function("a_dual2", |b| { + b.iter(|| a_res((&state.eos, &state.derive2(Derivative::DV)))) + }); group.bench_function("a_hyperdual", |b| { b.iter(|| a_res((&state.eos, &state.derive2_mixed(Derivative::DV, Derivative::DV)))) }); @@ -66,7 +69,7 @@ fn pcsaft(c: &mut Criterion) { IdentifierOption::Name, ) .unwrap(); - bench_dual_numbers(c, "methane", state_pcsaft(parameters)); + bench_dual_numbers(c, "dual_numbers_pcsaft_methane", state_pcsaft(parameters)); // water (4C, polar) let parameters = PcSaftParameters::from_json( @@ -76,7 +79,7 @@ fn pcsaft(c: &mut Criterion) { IdentifierOption::Name, ) .unwrap(); - bench_dual_numbers(c, "water_4c_polar", state_pcsaft(parameters)); + bench_dual_numbers(c, "dual_numbers_pcsaft_water_4c_polar", state_pcsaft(parameters)); // methane, ethane, propane let parameters = PcSaftParameters::from_json( @@ -86,7 +89,7 @@ fn pcsaft(c: &mut Criterion) { IdentifierOption::Name, ) .unwrap(); - bench_dual_numbers(c, "methane_ethane_propane", state_pcsaft(parameters)); + bench_dual_numbers(c, "dual_numbers_pcsaft_methane_ethane_propane", state_pcsaft(parameters)); } /// Benchmark for the PC-SAFT equation of state. @@ -116,7 +119,7 @@ fn methane_co2_pcsaft(c: &mut Criterion) { let x = arr1(&[0.15, 0.85]); let moles = &x * 10.0 * MOL; let state = State::new_nvt(&eos, temperature, volume, &moles).unwrap(); - bench_dual_numbers(c, "methane_co2", state); + bench_dual_numbers(c, "dual_numbers_pcsaft_methane_co2", state); } criterion_group!(bench, pcsaft, methane_co2_pcsaft); diff --git a/benches/state_creation.rs b/benches/state_creation.rs index ad7e5753f..81d4aee75 100644 --- a/benches/state_creation.rs +++ b/benches/state_creation.rs @@ -2,16 +2,16 @@ use criterion::{criterion_group, criterion_main, Criterion}; use feos::pcsaft::{PcSaft, PcSaftParameters}; use feos_core::{ parameter::{IdentifierOption, Parameter}, - DensityInitialization, PhaseEquilibrium, State, + Contributions, DensityInitialization, EquationOfState, PhaseEquilibrium, State, }; -use ndarray::arr1; +use ndarray::{Array, Array1}; use quantity::si::*; use std::sync::Arc; /// Evaluate NPT constructor -fn npt( +fn npt( (eos, t, p, n, rho0): ( - &Arc, + &Arc, SINumber, SINumber, &SIArray1, @@ -22,38 +22,138 @@ fn npt( } /// Evaluate critical point constructor -fn critical_point((eos, n): (&Arc, Option<&SIArray1>)) { +fn critical_point((eos, n): (&Arc, Option<&SIArray1>)) { State::critical_point(eos, n, None, Default::default()).unwrap(); } /// Evaluate temperature, pressure flash. -fn tp_flash((eos, t, p, feed): (&Arc, SINumber, SINumber, &SIArray1)) { +fn tp_flash((eos, t, p, feed): (&Arc, SINumber, SINumber, &SIArray1)) { PhaseEquilibrium::tp_flash(eos, t, p, feed, None, Default::default(), None).unwrap(); } -fn states_pcsaft(c: &mut Criterion) { - let parameters = PcSaftParameters::from_json( - vec!["methane", "ethane", "propane"], - "./parameters/pcsaft/gross2001.json", +fn bubble_point((eos, t, x): (&Arc, SINumber, &Array1)) { + PhaseEquilibrium::bubble_point( + eos, + t, + x, None, - IdentifierOption::Name, + None, + (Default::default(), Default::default()), ) .unwrap(); - let eos = Arc::new(PcSaft::new(Arc::new(parameters))); - let t = 300.0 * KELVIN; - let x = arr1(&[1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0]); +} + +fn dew_point((eos, t, y): (&Arc, SINumber, &Array1)) { + PhaseEquilibrium::dew_point( + eos, + t, + y, + None, + None, + (Default::default(), Default::default()), + ) + .unwrap(); +} + +fn bench_states(c: &mut Criterion, group_name: &str, eos: &Arc) { + let ncomponents = eos.components(); + let x = Array::from_elem(ncomponents, 1.0 / ncomponents as f64); let n = &x * 100.0 * MOL; - let p = BAR; + let crit = State::critical_point(&eos, Some(&n), None, Default::default()).unwrap(); + let vle = if ncomponents == 1 { + PhaseEquilibrium::pure(&eos, crit.temperature * 0.95, None, Default::default()).unwrap() + } else { + PhaseEquilibrium::tp_flash( + &eos, + crit.temperature, + crit.pressure(Contributions::Total) * 0.95, + &crit.moles, + None, + Default::default(), + None, + ) + .unwrap() + }; - let mut group = c.benchmark_group("methane_ethane_propane"); - group.bench_function("NPT", |b| { - b.iter(|| npt((&eos, t, p, &n, DensityInitialization::None))) + let mut group = c.benchmark_group(group_name); + group.bench_function("new_npt_liquid", |b| { + b.iter(|| { + npt(( + &eos, + vle.liquid().temperature, + vle.liquid().pressure(Contributions::Total) * 1.01, + &n, + DensityInitialization::Liquid, + )) + }) + }); + group.bench_function("new_npt_vapor", |b| { + b.iter(|| { + npt(( + &eos, + vle.vapor().temperature, + vle.vapor().pressure(Contributions::Total) * 0.99, + &n, + DensityInitialization::Vapor, + )) + }) }); group.bench_function("critical_point", |b| { b.iter(|| critical_point((&eos, Some(&n)))) }); - group.bench_function("tp_flash", |b| b.iter(|| tp_flash((&eos, 315.0 * KELVIN, 72.0 * BAR, &n)))); + if ncomponents != 1 { + group.bench_function("tp_flash", |b| { + b.iter(|| { + tp_flash(( + &eos, + crit.temperature, + crit.pressure(Contributions::Total) * 0.99, + &n, + )) + }) + }); + + group.bench_function("bubble_point", |b| { + b.iter(|| bubble_point((&eos, vle.liquid().temperature, &vle.liquid().molefracs))) + }); + + group.bench_function("dew_point", |b| { + b.iter(|| dew_point((&eos, vle.vapor().temperature, &vle.vapor().molefracs))) + }); + } +} + +fn pcsaft(c: &mut Criterion) { + let parameters = PcSaftParameters::from_json( + vec!["methane"], + "./parameters/pcsaft/gross2001.json", + None, + IdentifierOption::Name, + ) + .unwrap(); + let eos = Arc::new(PcSaft::new(Arc::new(parameters))); + bench_states(c, "state_creation_pcsaft_methane", &eos); + + let parameters = PcSaftParameters::from_json( + vec!["methane", "ethane", "propane"], + "./parameters/pcsaft/gross2001.json", + None, + IdentifierOption::Name, + ) + .unwrap(); + let eos = Arc::new(PcSaft::new(Arc::new(parameters))); + bench_states(c, "state_creation_pcsaft_methane_ethane", &eos); + + let parameters = PcSaftParameters::from_json( + vec!["methane", "ethane", "propane"], + "./parameters/pcsaft/gross2001.json", + None, + IdentifierOption::Name, + ) + .unwrap(); + let eos = Arc::new(PcSaft::new(Arc::new(parameters))); + bench_states(c, "state_creation_pcsaft_methane_ethane_propane", &eos); } -criterion_group!(bench, states_pcsaft); +criterion_group!(bench, pcsaft); criterion_main!(bench); diff --git a/benches/state_properties.rs b/benches/state_properties.rs index e9b386c2b..db4373568 100644 --- a/benches/state_properties.rs +++ b/benches/state_properties.rs @@ -50,7 +50,7 @@ fn properties_pcsaft(c: &mut Criterion) { let x = arr1(&[1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0]); let m = &x * 100.0 * MOL; - let mut group = c.benchmark_group("methane_ethane_propane"); + let mut group = c.benchmark_group("state_properties_pcsaft_methane_ethane_propane"); group.bench_function("a", |b| { b.iter(|| property((&eos, S::helmholtz_energy, t, v, &m, Contributions::Total))) }); From 524632ebef2dc85cffe143480eaf315cfc2b85c1 Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Tue, 20 Dec 2022 15:34:29 +0100 Subject: [PATCH 6/9] Use derive methods of dual numbers in the derive method for state --- feos-core/src/state/mod.rs | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/feos-core/src/state/mod.rs b/feos-core/src/state/mod.rs index 159972784..db6778e5a 100644 --- a/feos-core/src/state/mod.rs +++ b/feos-core/src/state/mod.rs @@ -686,9 +686,9 @@ impl State { let mut v = Dual64::from(self.reduced_volume); let mut n = self.reduced_moles.mapv(Dual64::from); match derivative { - Derivative::DT => t.eps[0] = 1.0, - Derivative::DV => v.eps[0] = 1.0, - Derivative::DN(i) => n[i].eps[0] = 1.0, + Derivative::DT => t = t.derive(), + Derivative::DV => v = v.derive(), + Derivative::DN(i) => n[i] = n[i].derive(), } StateHD::new(t, v, n) } @@ -716,14 +716,14 @@ impl State { let mut v = HyperDual64::from(self.reduced_volume); let mut n = self.reduced_moles.mapv(HyperDual64::from); match derivative1 { - Derivative::DT => t.eps1[0] = 1.0, - Derivative::DV => v.eps1[0] = 1.0, - Derivative::DN(i) => n[i].eps1[0] = 1.0, + Derivative::DT => t = t.derive1(), + Derivative::DV => v = v.derive1(), + Derivative::DN(i) => n[i] = n[i].derive1(), } match derivative2 { - Derivative::DT => t.eps2[0] = 1.0, - Derivative::DV => v.eps2[0] = 1.0, - Derivative::DN(i) => n[i].eps2[0] = 1.0, + Derivative::DT => t = t.derive2(), + Derivative::DV => v = v.derive2(), + Derivative::DN(i) => n[i] = n[i].derive2(), } StateHD::new(t, v, n) } From 806629f9205e12b8a83586b8df2f621aa28cf36f Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Tue, 20 Dec 2022 16:05:18 +0100 Subject: [PATCH 7/9] Updated readme for benchmarks, added benchmarks for VLE of pure substances for given temperature and pressure --- benches/README.md | 3 ++- benches/state_creation.rs | 12 ++++++++++++ 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/benches/README.md b/benches/README.md index bc9ab59e5..0689e2cbe 100644 --- a/benches/README.md +++ b/benches/README.md @@ -13,4 +13,5 @@ cargo bench --profile=release-lto --features=pcsaft --bench=dual_numbers |Name|Description|Cargo features| |--|--|--| |`dual_numbers`|Helmholtz energy function evaluated using `StateHD` with different dual number types.|`pcsaft`| -|`state_properties`|Properties of `State`. Including state creation using the natural variables of the Helmholtz energy (no density iteration).|`pcsaft`| \ No newline at end of file +|`state_properties`|Properties of `State`. Including state creation using the natural variables of the Helmholtz energy (no density iteration).|`pcsaft`| +|`state_creation`|Different constructors of `State` and `PhaseEquilibrium` including critical point calculations. For pure substances and mixtures.|`pcsaft`| \ No newline at end of file diff --git a/benches/state_creation.rs b/benches/state_creation.rs index 81d4aee75..55f1f3a11 100644 --- a/benches/state_creation.rs +++ b/benches/state_creation.rs @@ -26,6 +26,11 @@ fn critical_point((eos, n): (&Arc, Option<&SIArray1>)) { State::critical_point(eos, n, None, Default::default()).unwrap(); } +/// VLE for pure substance for given temperature or pressure +fn pure((eos, t_or_p): (&Arc, SINumber)) { + PhaseEquilibrium::pure(eos, t_or_p, None, Default::default()).unwrap(); +} + /// Evaluate temperature, pressure flash. fn tp_flash((eos, t, p, feed): (&Arc, SINumber, SINumber, &SIArray1)) { PhaseEquilibrium::tp_flash(eos, t, p, feed, None, Default::default(), None).unwrap(); @@ -120,6 +125,13 @@ fn bench_states(c: &mut Criterion, group_name: &str, eos: &A group.bench_function("dew_point", |b| { b.iter(|| dew_point((&eos, vle.vapor().temperature, &vle.vapor().molefracs))) }); + } else { + group.bench_function("pure_t", |b| { + b.iter(|| pure((&eos, vle.vapor().temperature))) + }); + group.bench_function("pure_p", |b| { + b.iter(|| pure((&eos, vle.vapor().pressure(Contributions::Total)))) + }); } } From 51c628404d48f6ce5ae2a33de46d0c6b5b7f3128 Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Tue, 20 Dec 2022 16:14:10 +0100 Subject: [PATCH 8/9] Implemented HelmholtzEnergy trait for Python StateHD-version that uses Dual2_64 --- feos-core/src/python/statehd.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/feos-core/src/python/statehd.rs b/feos-core/src/python/statehd.rs index 25af15a2b..c1f70f99e 100644 --- a/feos-core/src/python/statehd.rs +++ b/feos-core/src/python/statehd.rs @@ -179,6 +179,7 @@ impl_state_hd!( HyperDual ); impl_state_hd!(PyStateD, PyDual64, StateHD, Dual64); +impl_state_hd!(PyStateD2, PyDual2_64, StateHD, Dual2_64); impl_state_hd!(PyStateD3, PyDual3_64, StateHD, Dual3_64); impl_state_hd!( PyStateD3D, From 200ef35ae2ff964ff2fdf1392030f50611c69de2 Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Tue, 20 Dec 2022 18:03:20 +0100 Subject: [PATCH 9/9] Updated changelog of core --- feos-core/CHANGELOG.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/feos-core/CHANGELOG.md b/feos-core/CHANGELOG.md index 475f68937..d428a2d69 100644 --- a/feos-core/CHANGELOG.md +++ b/feos-core/CHANGELOG.md @@ -10,6 +10,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed - Added `Sync` and `Send` as supertraits to `EquationOfState`. [#57](https://github.com/feos-org/feos/pull/57) +- Added `Dual2_64` dual number to `HelmholtzEnergy` trait to facilitate efficient non-mixed second order derivatives. [#94](https://github.com/feos-org/feos/pull/94) +- Renamed `PartialDerivative::Second` to `PartialDerivative::SecondMixed`. [#94](https://github.com/feos-org/feos/pull/94) +- Added `PartialDerivative::Second` to enable non-mixed second order partial derivatives. [#94](https://github.com/feos-org/feos/pull/94) +- Changed `dp_dv_` and `ds_dt_` to use `Dual2_64` instead of `HyperDual64`. [#94](https://github.com/feos-org/feos/pull/94) +- Added `get_or_insert_with_d2_64` to `Cache`. [#94](https://github.com/feos-org/feos/pull/94) ## [0.3.1] - 2022-08-25 ### Added