Skip to content

Use single function evaluation in the calculation of critical points #96

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
prehner merged 3 commits into from
Dec 21, 2022
Merged

Use single function evaluation in the calculation of critical points #96

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
3294cf1
Use single function evaluation in the calculation of critical points
prehner Dec 20, 2022
d137524
update changelog
prehner Dec 21, 2022
7adf2ef
finetune changelog entry
prehner Dec 21, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions feos-core/CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
- 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)
- 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)

## [0.3.1] - 2022-08-25
### Added
Expand Down
4 changes: 3 additions & 1 deletion feos-core/src/equation_of_state.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,9 @@ 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, Dual2_64};
use num_dual::{
Dual, Dual2_64, Dual3, Dual3_64, Dual64, DualNum, DualVec64, HyperDual, HyperDual64,
};
use num_traits::{One, Zero};
use quantity::{QuantityArray1, QuantityScalar};
use std::fmt;
Expand Down
6 changes: 1 addition & 5 deletions feos-core/src/phase_equilibria/bubble_dew.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -488,11 +488,7 @@ where
let mut lnpstep = -f / df;

// catch too big p-steps
if lnpstep < -MAX_LNPSTEP {
lnpstep = -MAX_LNPSTEP;
} else if lnpstep > MAX_LNPSTEP {
lnpstep = MAX_LNPSTEP;
}
lnpstep = lnpstep.clamp(-MAX_LNPSTEP, MAX_LNPSTEP);

// Update p
*p = *p * lnpstep.exp();
Expand Down
34 changes: 16 additions & 18 deletions feos-core/src/state/critical_point.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -123,17 +123,15 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {

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

// calculate Newton step
let h = arr2(&[
[res_t[0].eps[0], res_r[0].eps[0]],
[res_t[1].eps[0], res_r[1].eps[0]],
]);
let mut delta = LU::new(h)?.solve(&res);

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

fn critical_point_objective<E: EquationOfState>(
eos: &Arc<E>,
temperature: Dual64,
density: Dual64,
temperature: DualVec64<2>,
density: DualVec64<2>,
moles: &Array1<f64>,
) -> EosResult<Array1<Dual64>> {
) -> EosResult<StaticVec<DualVec64<2>, 2>> {
// calculate second partial derivatives w.r.t. moles
let t = HyperDual::from_re(temperature);
let v = HyperDual::from_re(density.recip() * moles.sum());
let qij = Array2::from_shape_fn((eos.components(), eos.components()), |(i, j)| {
let mut m = moles.mapv(HyperDual::from);
m[i].eps1[0] = Dual64::one();
m[j].eps2[0] = Dual64::one();
m[i].eps1[0] = DualVec64::one();
m[j].eps2[0] = DualVec64::one();
let state = StateHD::new(t, v, m);
(eos.evaluate_residual(&state).eps1eps2[(0, 0)]
+ eos.ideal_gas().evaluate(&state).eps1eps2[(0, 0)])
Expand All @@ -493,10 +491,10 @@ fn critical_point_objective<E: EquationOfState>(
// evaluate third partial derivative w.r.t. s
let moles_hd = Array1::from_shape_fn(eos.components(), |i| {
Dual3::new(
Dual64::from(moles[i]),
DualVec64::from(moles[i]),
evec[i] * moles[i].sqrt(),
Dual64::zero(),
Dual64::zero(),
DualVec64::zero(),
DualVec64::zero(),
)
});
let state_s = StateHD::new(
Expand All @@ -505,7 +503,7 @@ fn critical_point_objective<E: EquationOfState>(
moles_hd,
);
let res = eos.evaluate_residual(&state_s) + eos.ideal_gas().evaluate(&state_s);
Ok(arr1(&[eval, res.v3]))
Ok(StaticVec::new_vec([eval, res.v3]))
}

fn critical_point_objective_t<E: EquationOfState>(
Expand Down
10 changes: 5 additions & 5 deletions feos-core/src/state/properties.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,34 +78,34 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
let new_state = self.derive0();
let computation =
|| self.eos.evaluate_residual(&new_state) * new_state.temperature;
cache.get_or_insert_with_f64(&computation) * U::reference_energy()
cache.get_or_insert_with_f64(computation) * U::reference_energy()
}
PartialDerivative::First(v) => {
let new_state = self.derive1(v);
let computation =
|| self.eos.evaluate_residual(&new_state) * new_state.temperature;
cache.get_or_insert_with_d64(v, &computation) * U::reference_energy()
cache.get_or_insert_with_d64(v, computation) * U::reference_energy()
/ v.reference()
}
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()
cache.get_or_insert_with_d2_64(v, computation) * U::reference_energy()
/ (v.reference() * v.reference())
}
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()
cache.get_or_insert_with_hd64(v1, v2, computation) * U::reference_energy()
/ (v1.reference() * v2.reference())
}
PartialDerivative::Third(v) => {
let new_state = self.derive3(v);
let computation =
|| self.eos.evaluate_residual(&new_state) * new_state.temperature;
cache.get_or_insert_with_hd364(v, &computation) * U::reference_energy()
cache.get_or_insert_with_hd364(v, computation) * U::reference_energy()
/ (v.reference() * v.reference() * v.reference())
}
}),
Expand Down