diff --git a/feos-core/CHANGELOG.md b/feos-core/CHANGELOG.md index d428a2d69..bc26a00a2 100644 --- a/feos-core/CHANGELOG.md +++ b/feos-core/CHANGELOG.md @@ -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 diff --git a/feos-core/src/equation_of_state.rs b/feos-core/src/equation_of_state.rs index 406da7d0a..8677be69e 100644 --- a/feos-core/src/equation_of_state.rs +++ b/feos-core/src/equation_of_state.rs @@ -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; diff --git a/feos-core/src/phase_equilibria/bubble_dew.rs b/feos-core/src/phase_equilibria/bubble_dew.rs index c2217793d..0c2bbdff9 100644 --- a/feos-core/src/phase_equilibria/bubble_dew.rs +++ b/feos-core/src/phase_equilibria/bubble_dew.rs @@ -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(); diff --git a/feos-core/src/state/critical_point.rs b/feos-core/src/state/critical_point.rs index 2976b6f59..3007b11fc 100644 --- a/feos-core/src/state/critical_point.rs +++ b/feos-core/src/state/critical_point.rs @@ -123,17 +123,15 @@ impl State { 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 @@ -470,17 +468,17 @@ impl State { fn critical_point_objective( eos: &Arc, - temperature: Dual64, - density: Dual64, + temperature: DualVec64<2>, + density: DualVec64<2>, moles: &Array1, -) -> EosResult> { +) -> EosResult, 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)]) @@ -493,10 +491,10 @@ fn critical_point_objective( // 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( @@ -505,7 +503,7 @@ fn critical_point_objective( 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( diff --git a/feos-core/src/state/properties.rs b/feos-core/src/state/properties.rs index 9a580926b..ecc404d88 100644 --- a/feos-core/src/state/properties.rs +++ b/feos-core/src/state/properties.rs @@ -78,34 +78,34 @@ impl State { 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()) } }),