Skip to content

Commit 6769a63

Browse files
prehnerprehner
authored
Use single function evaluation in the calculation of critical points (#96)
* Use single function evaluation in the calculation of critical points * update changelog * finetune changelog entry
1 parent a1785cb commit 6769a63

5 files changed

Lines changed: 26 additions & 29 deletions

File tree

‎feos-core/CHANGELOG.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
1515
- Added `PartialDerivative::Second` to enable non-mixed second order partial derivatives. [#94](https://github.com/feos-org/feos/pull/94)
1616
- Changed `dp_dv_` and `ds_dt_` to use `Dual2_64` instead of `HyperDual64`. [#94](https://github.com/feos-org/feos/pull/94)
1717
- Added `get_or_insert_with_d2_64` to `Cache`. [#94](https://github.com/feos-org/feos/pull/94)
18+
- The critical point algorithm now uses vector dual numbers to reduce the number of model evaluations and computation times. [#96](https://github.com/feos-org/feos/pull/96)
1819

1920
## [0.3.1] - 2022-08-25
2021
### Added

‎feos-core/src/equation_of_state.rs‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,9 @@ 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, Dual2_64};
5+
use num_dual::{
6+
Dual, Dual2_64, Dual3, Dual3_64, Dual64, DualNum, DualVec64, HyperDual, HyperDual64,
7+
};
68
use num_traits::{One, Zero};
79
use quantity::{QuantityArray1, QuantityScalar};
810
use std::fmt;

‎feos-core/src/phase_equilibria/bubble_dew.rs‎

Lines changed: 1 addition & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -488,11 +488,7 @@ where
488488
let mut lnpstep = -f / df;
489489

490490
// catch too big p-steps
491-
if lnpstep < -MAX_LNPSTEP {
492-
lnpstep = -MAX_LNPSTEP;
493-
} else if lnpstep > MAX_LNPSTEP {
494-
lnpstep = MAX_LNPSTEP;
495-
}
491+
lnpstep = lnpstep.clamp(-MAX_LNPSTEP, MAX_LNPSTEP);
496492

497493
// Update p
498494
*p = *p * lnpstep.exp();

‎feos-core/src/state/critical_point.rs‎

Lines changed: 16 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -123,17 +123,15 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
123123

124124
for i in 1..=max_iter {
125125
// calculate residuals and derivatives w.r.t. temperature and density
126-
let res_t =
127-
critical_point_objective(eos, Dual64::from(t).derive(), Dual64::from(rho), &n)?;
128-
let res_r =
129-
critical_point_objective(eos, Dual64::from(t), Dual64::from(rho).derive(), &n)?;
130-
let res = res_t.map(Dual64::re);
126+
let [t_dual, rho_dual] = *StaticVec::new_vec([t, rho])
127+
.map(DualVec64::<2>::from_re)
128+
.derive()
129+
.raw_array();
130+
let res = critical_point_objective(eos, t_dual, rho_dual, &n)?;
131+
let h = arr2(res.jacobian().raw_data());
132+
let res = arr1(res.map(|r| r.re()).raw_array());
131133

132134
// calculate Newton step
133-
let h = arr2(&[
134-
[res_t[0].eps[0], res_r[0].eps[0]],
135-
[res_t[1].eps[0], res_r[1].eps[0]],
136-
]);
137135
let mut delta = LU::new(h)?.solve(&res);
138136

139137
// reduce step if necessary
@@ -470,17 +468,17 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
470468

471469
fn critical_point_objective<E: EquationOfState>(
472470
eos: &Arc<E>,
473-
temperature: Dual64,
474-
density: Dual64,
471+
temperature: DualVec64<2>,
472+
density: DualVec64<2>,
475473
moles: &Array1<f64>,
476-
) -> EosResult<Array1<Dual64>> {
474+
) -> EosResult<StaticVec<DualVec64<2>, 2>> {
477475
// calculate second partial derivatives w.r.t. moles
478476
let t = HyperDual::from_re(temperature);
479477
let v = HyperDual::from_re(density.recip() * moles.sum());
480478
let qij = Array2::from_shape_fn((eos.components(), eos.components()), |(i, j)| {
481479
let mut m = moles.mapv(HyperDual::from);
482-
m[i].eps1[0] = Dual64::one();
483-
m[j].eps2[0] = Dual64::one();
480+
m[i].eps1[0] = DualVec64::one();
481+
m[j].eps2[0] = DualVec64::one();
484482
let state = StateHD::new(t, v, m);
485483
(eos.evaluate_residual(&state).eps1eps2[(0, 0)]
486484
+ eos.ideal_gas().evaluate(&state).eps1eps2[(0, 0)])
@@ -493,10 +491,10 @@ fn critical_point_objective<E: EquationOfState>(
493491
// evaluate third partial derivative w.r.t. s
494492
let moles_hd = Array1::from_shape_fn(eos.components(), |i| {
495493
Dual3::new(
496-
Dual64::from(moles[i]),
494+
DualVec64::from(moles[i]),
497495
evec[i] * moles[i].sqrt(),
498-
Dual64::zero(),
499-
Dual64::zero(),
496+
DualVec64::zero(),
497+
DualVec64::zero(),
500498
)
501499
});
502500
let state_s = StateHD::new(
@@ -505,7 +503,7 @@ fn critical_point_objective<E: EquationOfState>(
505503
moles_hd,
506504
);
507505
let res = eos.evaluate_residual(&state_s) + eos.ideal_gas().evaluate(&state_s);
508-
Ok(arr1(&[eval, res.v3]))
506+
Ok(StaticVec::new_vec([eval, res.v3]))
509507
}
510508

511509
fn critical_point_objective_t<E: EquationOfState>(

‎feos-core/src/state/properties.rs‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -78,34 +78,34 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
7878
let new_state = self.derive0();
7979
let computation =
8080
|| self.eos.evaluate_residual(&new_state) * new_state.temperature;
81-
cache.get_or_insert_with_f64(&computation) * U::reference_energy()
81+
cache.get_or_insert_with_f64(computation) * U::reference_energy()
8282
}
8383
PartialDerivative::First(v) => {
8484
let new_state = self.derive1(v);
8585
let computation =
8686
|| self.eos.evaluate_residual(&new_state) * new_state.temperature;
87-
cache.get_or_insert_with_d64(v, &computation) * U::reference_energy()
87+
cache.get_or_insert_with_d64(v, computation) * U::reference_energy()
8888
/ v.reference()
8989
}
9090
PartialDerivative::Second(v) => {
9191
let new_state = self.derive2(v);
9292
let computation =
9393
|| self.eos.evaluate_residual(&new_state) * new_state.temperature;
94-
cache.get_or_insert_with_d2_64(v, &computation) * U::reference_energy()
94+
cache.get_or_insert_with_d2_64(v, computation) * U::reference_energy()
9595
/ (v.reference() * v.reference())
9696
}
9797
PartialDerivative::SecondMixed(v1, v2) => {
9898
let new_state = self.derive2_mixed(v1, v2);
9999
let computation =
100100
|| self.eos.evaluate_residual(&new_state) * new_state.temperature;
101-
cache.get_or_insert_with_hd64(v1, v2, &computation) * U::reference_energy()
101+
cache.get_or_insert_with_hd64(v1, v2, computation) * U::reference_energy()
102102
/ (v1.reference() * v2.reference())
103103
}
104104
PartialDerivative::Third(v) => {
105105
let new_state = self.derive3(v);
106106
let computation =
107107
|| self.eos.evaluate_residual(&new_state) * new_state.temperature;
108-
cache.get_or_insert_with_hd364(v, &computation) * U::reference_energy()
108+
cache.get_or_insert_with_hd364(v, computation) * U::reference_energy()
109109
/ (v.reference() * v.reference() * v.reference())
110110
}
111111
}),

0 commit comments

Comments
 (0)