Skip to content

Commit df9c2a9

Browse files
g-bauerg-bauer
authored
Add non-mixed second order derivatives to State property evaluation (#94)
* add second (non-mixed) derivative to cache and property evaluation * renamed derivative to SecondMixed, renamed some tests * change cache to only store second order derivatives in SecondMixed * remove redundant hashmap entry in cache * Add more benchmarks and systems, better names for benchmark groups * Use derive methods of dual numbers in the derive method for state * Updated readme for benchmarks, added benchmarks for VLE of pure substances for given temperature and pressure * Implemented HelmholtzEnergy trait for Python StateHD-version that uses Dual2_64 * Updated changelog of core
1 parent 2d821f3 commit df9c2a9

File tree

13 files changed

+303
-40
lines changed

13 files changed

+303
-40
lines changed

‎Cargo.toml‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,10 @@ criterion = "0.4"
5252
name = "state_properties"
5353
harness = false
5454

55+
[[bench]]
56+
name = "state_creation"
57+
harness = false
58+
5559
[[bench]]
5660
name = "dual_numbers"
5761
harness = false

‎benches/README.md‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,4 +13,5 @@ cargo bench --profile=release-lto --features=pcsaft --bench=dual_numbers
1313
|Name|Description|Cargo features|
1414
|--|--|--|
1515
|`dual_numbers`|Helmholtz energy function evaluated using `StateHD` with different dual number types.|`pcsaft`|
16-
|`state_properties`|Properties of `State`. Including state creation using the natural variables of the Helmholtz energy (no density iteration).|`pcsaft`|
16+
|`state_properties`|Properties of `State`. Including state creation using the natural variables of the Helmholtz energy (no density iteration).|`pcsaft`|
17+
|`state_creation`|Different constructors of `State` and `PhaseEquilibrium` including critical point calculations. For pure substances and mixtures.|`pcsaft`|

‎benches/dual_numbers.rs‎

Lines changed: 8 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -48,8 +48,11 @@ fn bench_dual_numbers<E: EquationOfState>(
4848
group.bench_function("a_dual", |b| {
4949
b.iter(|| a_res((&state.eos, &state.derive1(Derivative::DV))))
5050
});
51+
group.bench_function("a_dual2", |b| {
52+
b.iter(|| a_res((&state.eos, &state.derive2(Derivative::DV))))
53+
});
5154
group.bench_function("a_hyperdual", |b| {
52-
b.iter(|| a_res((&state.eos, &state.derive2(Derivative::DV, Derivative::DV))))
55+
b.iter(|| a_res((&state.eos, &state.derive2_mixed(Derivative::DV, Derivative::DV))))
5356
});
5457
group.bench_function("a_dual3", |b| {
5558
b.iter(|| a_res((&state.eos, &state.derive3(Derivative::DV))))
@@ -66,7 +69,7 @@ fn pcsaft(c: &mut Criterion) {
6669
IdentifierOption::Name,
6770
)
6871
.unwrap();
69-
bench_dual_numbers(c, "methane", state_pcsaft(parameters));
72+
bench_dual_numbers(c, "dual_numbers_pcsaft_methane", state_pcsaft(parameters));
7073

7174
// water (4C, polar)
7275
let parameters = PcSaftParameters::from_json(
@@ -76,7 +79,7 @@ fn pcsaft(c: &mut Criterion) {
7679
IdentifierOption::Name,
7780
)
7881
.unwrap();
79-
bench_dual_numbers(c, "water_4c_polar", state_pcsaft(parameters));
82+
bench_dual_numbers(c, "dual_numbers_pcsaft_water_4c_polar", state_pcsaft(parameters));
8083

8184
// methane, ethane, propane
8285
let parameters = PcSaftParameters::from_json(
@@ -86,7 +89,7 @@ fn pcsaft(c: &mut Criterion) {
8689
IdentifierOption::Name,
8790
)
8891
.unwrap();
89-
bench_dual_numbers(c, "methane_ethane_propane", state_pcsaft(parameters));
92+
bench_dual_numbers(c, "dual_numbers_pcsaft_methane_ethane_propane", state_pcsaft(parameters));
9093
}
9194

9295
/// Benchmark for the PC-SAFT equation of state.
@@ -116,7 +119,7 @@ fn methane_co2_pcsaft(c: &mut Criterion) {
116119
let x = arr1(&[0.15, 0.85]);
117120
let moles = &x * 10.0 * MOL;
118121
let state = State::new_nvt(&eos, temperature, volume, &moles).unwrap();
119-
bench_dual_numbers(c, "methane_co2", state);
122+
bench_dual_numbers(c, "dual_numbers_pcsaft_methane_co2", state);
120123
}
121124

122125
criterion_group!(bench, pcsaft, methane_co2_pcsaft);

‎benches/state_creation.rs‎

Lines changed: 171 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,171 @@
1+
use criterion::{criterion_group, criterion_main, Criterion};
2+
use feos::pcsaft::{PcSaft, PcSaftParameters};
3+
use feos_core::{
4+
parameter::{IdentifierOption, Parameter},
5+
Contributions, DensityInitialization, EquationOfState, PhaseEquilibrium, State,
6+
};
7+
use ndarray::{Array, Array1};
8+
use quantity::si::*;
9+
use std::sync::Arc;
10+
11+
/// Evaluate NPT constructor
12+
fn npt<E: EquationOfState>(
13+
(eos, t, p, n, rho0): (
14+
&Arc<E>,
15+
SINumber,
16+
SINumber,
17+
&SIArray1,
18+
DensityInitialization<SIUnit>,
19+
),
20+
) {
21+
State::new_npt(eos, t, p, n, rho0).unwrap();
22+
}
23+
24+
/// Evaluate critical point constructor
25+
fn critical_point<E: EquationOfState>((eos, n): (&Arc<E>, Option<&SIArray1>)) {
26+
State::critical_point(eos, n, None, Default::default()).unwrap();
27+
}
28+
29+
/// VLE for pure substance for given temperature or pressure
30+
fn pure<E: EquationOfState>((eos, t_or_p): (&Arc<E>, SINumber)) {
31+
PhaseEquilibrium::pure(eos, t_or_p, None, Default::default()).unwrap();
32+
}
33+
34+
/// Evaluate temperature, pressure flash.
35+
fn tp_flash<E: EquationOfState>((eos, t, p, feed): (&Arc<E>, SINumber, SINumber, &SIArray1)) {
36+
PhaseEquilibrium::tp_flash(eos, t, p, feed, None, Default::default(), None).unwrap();
37+
}
38+
39+
fn bubble_point<E: EquationOfState>((eos, t, x): (&Arc<E>, SINumber, &Array1<f64>)) {
40+
PhaseEquilibrium::bubble_point(
41+
eos,
42+
t,
43+
x,
44+
None,
45+
None,
46+
(Default::default(), Default::default()),
47+
)
48+
.unwrap();
49+
}
50+
51+
fn dew_point<E: EquationOfState>((eos, t, y): (&Arc<E>, SINumber, &Array1<f64>)) {
52+
PhaseEquilibrium::dew_point(
53+
eos,
54+
t,
55+
y,
56+
None,
57+
None,
58+
(Default::default(), Default::default()),
59+
)
60+
.unwrap();
61+
}
62+
63+
fn bench_states<E: EquationOfState>(c: &mut Criterion, group_name: &str, eos: &Arc<E>) {
64+
let ncomponents = eos.components();
65+
let x = Array::from_elem(ncomponents, 1.0 / ncomponents as f64);
66+
let n = &x * 100.0 * MOL;
67+
let crit = State::critical_point(&eos, Some(&n), None, Default::default()).unwrap();
68+
let vle = if ncomponents == 1 {
69+
PhaseEquilibrium::pure(&eos, crit.temperature * 0.95, None, Default::default()).unwrap()
70+
} else {
71+
PhaseEquilibrium::tp_flash(
72+
&eos,
73+
crit.temperature,
74+
crit.pressure(Contributions::Total) * 0.95,
75+
&crit.moles,
76+
None,
77+
Default::default(),
78+
None,
79+
)
80+
.unwrap()
81+
};
82+
83+
let mut group = c.benchmark_group(group_name);
84+
group.bench_function("new_npt_liquid", |b| {
85+
b.iter(|| {
86+
npt((
87+
&eos,
88+
vle.liquid().temperature,
89+
vle.liquid().pressure(Contributions::Total) * 1.01,
90+
&n,
91+
DensityInitialization::Liquid,
92+
))
93+
})
94+
});
95+
group.bench_function("new_npt_vapor", |b| {
96+
b.iter(|| {
97+
npt((
98+
&eos,
99+
vle.vapor().temperature,
100+
vle.vapor().pressure(Contributions::Total) * 0.99,
101+
&n,
102+
DensityInitialization::Vapor,
103+
))
104+
})
105+
});
106+
group.bench_function("critical_point", |b| {
107+
b.iter(|| critical_point((&eos, Some(&n))))
108+
});
109+
if ncomponents != 1 {
110+
group.bench_function("tp_flash", |b| {
111+
b.iter(|| {
112+
tp_flash((
113+
&eos,
114+
crit.temperature,
115+
crit.pressure(Contributions::Total) * 0.99,
116+
&n,
117+
))
118+
})
119+
});
120+
121+
group.bench_function("bubble_point", |b| {
122+
b.iter(|| bubble_point((&eos, vle.liquid().temperature, &vle.liquid().molefracs)))
123+
});
124+
125+
group.bench_function("dew_point", |b| {
126+
b.iter(|| dew_point((&eos, vle.vapor().temperature, &vle.vapor().molefracs)))
127+
});
128+
} else {
129+
group.bench_function("pure_t", |b| {
130+
b.iter(|| pure((&eos, vle.vapor().temperature)))
131+
});
132+
group.bench_function("pure_p", |b| {
133+
b.iter(|| pure((&eos, vle.vapor().pressure(Contributions::Total))))
134+
});
135+
}
136+
}
137+
138+
fn pcsaft(c: &mut Criterion) {
139+
let parameters = PcSaftParameters::from_json(
140+
vec!["methane"],
141+
"./parameters/pcsaft/gross2001.json",
142+
None,
143+
IdentifierOption::Name,
144+
)
145+
.unwrap();
146+
let eos = Arc::new(PcSaft::new(Arc::new(parameters)));
147+
bench_states(c, "state_creation_pcsaft_methane", &eos);
148+
149+
let parameters = PcSaftParameters::from_json(
150+
vec!["methane", "ethane", "propane"],
151+
"./parameters/pcsaft/gross2001.json",
152+
None,
153+
IdentifierOption::Name,
154+
)
155+
.unwrap();
156+
let eos = Arc::new(PcSaft::new(Arc::new(parameters)));
157+
bench_states(c, "state_creation_pcsaft_methane_ethane", &eos);
158+
159+
let parameters = PcSaftParameters::from_json(
160+
vec!["methane", "ethane", "propane"],
161+
"./parameters/pcsaft/gross2001.json",
162+
None,
163+
IdentifierOption::Name,
164+
)
165+
.unwrap();
166+
let eos = Arc::new(PcSaft::new(Arc::new(parameters)));
167+
bench_states(c, "state_creation_pcsaft_methane_ethane_propane", &eos);
168+
}
169+
170+
criterion_group!(bench, pcsaft);
171+
criterion_main!(bench);

‎benches/state_properties.rs‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -50,7 +50,7 @@ fn properties_pcsaft(c: &mut Criterion) {
5050
let x = arr1(&[1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0]);
5151
let m = &x * 100.0 * MOL;
5252

53-
let mut group = c.benchmark_group("methane_ethane_propane");
53+
let mut group = c.benchmark_group("state_properties_pcsaft_methane_ethane_propane");
5454
group.bench_function("a", |b| {
5555
b.iter(|| property((&eos, S::helmholtz_energy, t, v, &m, Contributions::Total)))
5656
});

‎feos-core/CHANGELOG.md‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
1010

1111
### Changed
1212
- Added `Sync` and `Send` as supertraits to `EquationOfState`. [#57](https://github.com/feos-org/feos/pull/57)
13+
- 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)
14+
- Renamed `PartialDerivative::Second` to `PartialDerivative::SecondMixed`. [#94](https://github.com/feos-org/feos/pull/94)
15+
- Added `PartialDerivative::Second` to enable non-mixed second order partial derivatives. [#94](https://github.com/feos-org/feos/pull/94)
16+
- Changed `dp_dv_` and `ds_dt_` to use `Dual2_64` instead of `HyperDual64`. [#94](https://github.com/feos-org/feos/pull/94)
17+
- Added `get_or_insert_with_d2_64` to `Cache`. [#94](https://github.com/feos-org/feos/pull/94)
1318

1419
## [0.3.1] - 2022-08-25
1520
### Added

‎feos-core/src/equation_of_state.rs‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@ use crate::errors::{EosError, EosResult};
22
use crate::state::StateHD;
33
use crate::EosUnit;
44
use ndarray::prelude::*;
5-
use num_dual::{Dual, Dual3, Dual3_64, Dual64, DualNum, DualVec64, HyperDual, HyperDual64};
5+
use num_dual::{Dual, Dual3, Dual3_64, Dual64, DualNum, DualVec64, HyperDual, HyperDual64, Dual2_64};
66
use num_traits::{One, Zero};
77
use quantity::{QuantityArray1, QuantityScalar};
88
use std::fmt;
@@ -28,6 +28,7 @@ pub trait HelmholtzEnergy:
2828
+ HelmholtzEnergyDual<Dual64>
2929
+ HelmholtzEnergyDual<Dual<DualVec64<3>, f64>>
3030
+ HelmholtzEnergyDual<HyperDual64>
31+
+ HelmholtzEnergyDual<Dual2_64>
3132
+ HelmholtzEnergyDual<Dual3_64>
3233
+ HelmholtzEnergyDual<HyperDual<Dual64, f64>>
3334
+ HelmholtzEnergyDual<HyperDual<DualVec64<2>, f64>>
@@ -46,6 +47,7 @@ impl<T> HelmholtzEnergy for T where
4647
+ HelmholtzEnergyDual<Dual64>
4748
+ HelmholtzEnergyDual<Dual<DualVec64<3>, f64>>
4849
+ HelmholtzEnergyDual<HyperDual64>
50+
+ HelmholtzEnergyDual<Dual2_64>
4951
+ HelmholtzEnergyDual<Dual3_64>
5052
+ HelmholtzEnergyDual<HyperDual<Dual64, f64>>
5153
+ HelmholtzEnergyDual<HyperDual<DualVec64<2>, f64>>
@@ -99,6 +101,7 @@ pub trait IdealGasContribution:
99101
+ IdealGasContributionDual<Dual64>
100102
+ IdealGasContributionDual<Dual<DualVec64<3>, f64>>
101103
+ IdealGasContributionDual<HyperDual64>
104+
+ IdealGasContributionDual<Dual2_64>
102105
+ IdealGasContributionDual<Dual3_64>
103106
+ IdealGasContributionDual<HyperDual<Dual64, f64>>
104107
+ IdealGasContributionDual<HyperDual<DualVec64<2>, f64>>
@@ -115,6 +118,7 @@ impl<T> IdealGasContribution for T where
115118
+ IdealGasContributionDual<Dual64>
116119
+ IdealGasContributionDual<Dual<DualVec64<3>, f64>>
117120
+ IdealGasContributionDual<HyperDual64>
121+
+ IdealGasContributionDual<Dual2_64>
118122
+ IdealGasContributionDual<Dual3_64>
119123
+ IdealGasContributionDual<HyperDual<Dual64, f64>>
120124
+ IdealGasContributionDual<HyperDual<DualVec64<2>, f64>>

‎feos-core/src/python/statehd.rs‎

Lines changed: 15 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,9 @@
11
use crate::StateHD;
22
use ndarray::Array1;
3-
use num_dual::python::{PyDual3Dual64, PyDual3_64, PyDual64, PyHyperDual64, PyHyperDualDual64};
4-
use num_dual::{Dual, Dual3, Dual3_64, Dual64, DualVec64, HyperDual, HyperDual64};
3+
use num_dual::python::{
4+
PyDual2_64, PyDual3Dual64, PyDual3_64, PyDual64, PyHyperDual64, PyHyperDualDual64,
5+
};
6+
use num_dual::{Dual, Dual2_64, Dual3, Dual3_64, Dual64, DualVec64, HyperDual, HyperDual64};
57
use pyo3::prelude::*;
68
use std::convert::From;
79

@@ -75,6 +77,16 @@ impl From<StateHD<Dual<DualVec64<3>, f64>>> for PyStateDDV3 {
7577
}
7678
}
7779

80+
#[pyclass(name = "StateD2")]
81+
#[derive(Clone)]
82+
pub struct PyStateD2(StateHD<Dual2_64>);
83+
84+
impl From<StateHD<Dual2_64>> for PyStateD2 {
85+
fn from(s: StateHD<Dual2_64>) -> Self {
86+
Self(s)
87+
}
88+
}
89+
7890
#[pyclass(name = "StateD3")]
7991
#[derive(Clone)]
8092
pub struct PyStateD3(StateHD<Dual3_64>);
@@ -167,6 +179,7 @@ impl_state_hd!(
167179
HyperDual<Dual64, f64>
168180
);
169181
impl_state_hd!(PyStateD, PyDual64, StateHD<Dual64>, Dual64);
182+
impl_state_hd!(PyStateD2, PyDual2_64, StateHD<Dual2_64>, Dual2_64);
170183
impl_state_hd!(PyStateD3, PyDual3_64, StateHD<Dual3_64>, Dual3_64);
171184
impl_state_hd!(
172185
PyStateD3D,

‎feos-core/src/python/user_defined.rs‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,10 @@
11
use crate::python::statehd::*;
22
use crate::*;
33
use ndarray::prelude::*;
4-
use num_dual::python::{PyDual3Dual64, PyDual3_64, PyDual64, PyHyperDual64, PyHyperDualDual64};
5-
use num_dual::{Dual, Dual3, Dual3_64, Dual64, DualVec64, HyperDual, HyperDual64};
4+
use num_dual::python::{
5+
PyDual2_64, PyDual3Dual64, PyDual3_64, PyDual64, PyHyperDual64, PyHyperDualDual64,
6+
};
7+
use num_dual::{Dual, Dual2_64, Dual3, Dual3_64, Dual64, DualVec64, HyperDual, HyperDual64};
68
use numpy::convert::IntoPyArray;
79
use pyo3::prelude::*;
810
use quantity::python::PySIArray1;
@@ -215,6 +217,7 @@ impl_helmholtz_energy!(
215217
PyHyperDualDualVec64_3,
216218
HyperDual<DualVec64<3>, f64>
217219
);
220+
impl_helmholtz_energy!(PyStateD2, PyDual2_64, Dual2_64);
218221
impl_helmholtz_energy!(PyStateD3, PyDual3_64, Dual3_64);
219222
impl_helmholtz_energy!(PyStateD3D, PyDual3Dual64, Dual3<Dual64, f64>);
220223
impl_helmholtz_energy!(PyStateD3DV2, PyDual3DualVec64_2, Dual3<DualVec64<2>, f64>);

0 commit comments

Comments
 (0)