diff --git a/CHANGELOG.md b/CHANGELOG.md index e22425c..ffaef91 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## Unreleased ### Changed - Removed `State` from `EntropyScaling` trait and adjusted associated methods to use temperature, volume and moles instead of state. [#36](https://github.com/feos-org/feos-core/pull/36) +- Replaced the outer loop iterations for the critical point of binary systems with dedicated algorithms. [#34](https://github.com/feos-org/feos-core/pull/34) ## [0.1.5] - 2022-02-21 ### Fixed diff --git a/Cargo.toml b/Cargo.toml index b5aefb2..94faf5c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -23,7 +23,6 @@ num-dual = { version = "0.4", features = ["ndarray"] } ndarray = { version = "0.15", features = ["serde"] } num-traits = "0.2" thiserror = "1.0" -argmin = "0.4" serde = { version = "1.0", features = ["derive"] } serde_json = "1.0" indexmap = "1.7" diff --git a/src/equation_of_state.rs b/src/equation_of_state.rs index 8590ba5..9d3e86f 100644 --- a/src/equation_of_state.rs +++ b/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::{Dual3, Dual3_64, Dual64, DualNum, HyperDual, HyperDual64}; +use num_dual::{Dual, Dual3, Dual3_64, Dual64, DualNum, DualVec64, HyperDual, HyperDual64}; use num_traits::{One, Zero}; use quantity::{QuantityArray1, QuantityScalar}; use std::fmt; @@ -26,10 +26,15 @@ pub trait HelmholtzEnergyDual> { pub trait HelmholtzEnergy: HelmholtzEnergyDual + HelmholtzEnergyDual + + HelmholtzEnergyDual, f64>> + HelmholtzEnergyDual + HelmholtzEnergyDual + HelmholtzEnergyDual> + + HelmholtzEnergyDual, f64>> + + HelmholtzEnergyDual, f64>> + HelmholtzEnergyDual> + + HelmholtzEnergyDual, f64>> + + HelmholtzEnergyDual, f64>> + fmt::Display { } @@ -37,10 +42,15 @@ pub trait HelmholtzEnergy: impl HelmholtzEnergy for T where T: HelmholtzEnergyDual + HelmholtzEnergyDual + + HelmholtzEnergyDual, f64>> + HelmholtzEnergyDual + HelmholtzEnergyDual + HelmholtzEnergyDual> + + HelmholtzEnergyDual, f64>> + + HelmholtzEnergyDual, f64>> + HelmholtzEnergyDual> + + HelmholtzEnergyDual, f64>> + + HelmholtzEnergyDual, f64>> + fmt::Display { } @@ -83,10 +93,15 @@ pub trait IdealGasContributionDual> { pub trait IdealGasContribution: IdealGasContributionDual + IdealGasContributionDual + + IdealGasContributionDual, f64>> + IdealGasContributionDual + IdealGasContributionDual + IdealGasContributionDual> + + IdealGasContributionDual, f64>> + + IdealGasContributionDual, f64>> + IdealGasContributionDual> + + IdealGasContributionDual, f64>> + + IdealGasContributionDual, f64>> + fmt::Display { } @@ -94,15 +109,20 @@ pub trait IdealGasContribution: impl IdealGasContribution for T where T: IdealGasContributionDual + IdealGasContributionDual + + IdealGasContributionDual, f64>> + IdealGasContributionDual + IdealGasContributionDual + IdealGasContributionDual> + + IdealGasContributionDual, f64>> + + IdealGasContributionDual, f64>> + IdealGasContributionDual> + + IdealGasContributionDual, f64>> + + IdealGasContributionDual, f64>> + fmt::Display { } -struct DefaultIdealGasContribution(); +struct DefaultIdealGasContribution; impl> IdealGasContributionDual for DefaultIdealGasContribution { fn de_broglie_wavelength(&self, _: D, components: usize) -> Array1 { Array1::zeros(components) @@ -185,7 +205,7 @@ pub trait EquationOfState { /// required (e.g. for the calculation of enthalpies) this function /// has to be overwritten. fn ideal_gas(&self) -> &dyn IdealGasContribution { - &DefaultIdealGasContribution() + &DefaultIdealGasContribution } /// Check if the provided optional mole number is consistent with the diff --git a/src/errors.rs b/src/errors.rs index c746d34..d087ae3 100644 --- a/src/errors.rs +++ b/src/errors.rs @@ -1,5 +1,4 @@ use crate::parameter::ParameterError; -use argmin::core::Error as ArgminError; use num_dual::linalg::LinAlgError; use quantity::QuantityError; use thiserror::Error; @@ -28,8 +27,6 @@ pub enum EosError { #[error(transparent)] ParameterError(#[from] ParameterError), #[error(transparent)] - ArgminError(#[from] ArgminError), - #[error(transparent)] LinAlgError(#[from] LinAlgError), } diff --git a/src/parameter/mod.rs b/src/parameter/mod.rs index 6e9b0fd..5731d3d 100644 --- a/src/parameter/mod.rs +++ b/src/parameter/mod.rs @@ -200,8 +200,7 @@ where // empty, if no binary segment records are provided let binary_map: HashMap<_, _> = binary_segment_records .into_iter() - .map(|seg| seg.into_iter()) - .flatten() + .flat_map(|seg| seg.into_iter()) .map(|br| ((br.id1, br.id2), br.model_record)) .collect(); diff --git a/src/phase_equilibria/phase_diagram_binary.rs b/src/phase_equilibria/phase_diagram_binary.rs index c779093..bd7d280 100644 --- a/src/phase_equilibria/phase_diagram_binary.rs +++ b/src/phase_equilibria/phase_diagram_binary.rs @@ -1,10 +1,10 @@ use super::{PhaseEquilibrium, SolverOptions}; use crate::equation_of_state::EquationOfState; use crate::errors::{EosError, EosResult}; -use num_dual::linalg::{norm, LU}; use crate::state::{Contributions, DensityInitialization, State, StateBuilder, TPSpec}; use crate::EosUnit; use ndarray::{arr1, arr2, concatenate, s, Array1, Array2, Axis}; +use num_dual::linalg::{norm, LU}; use quantity::{QuantityArray1, QuantityScalar}; use std::rc::Rc; @@ -110,12 +110,14 @@ impl PhaseDiagramBinary { let (x_lim, vle_lim, bubble) = match vle_sat { [None, None] => return Err(EosError::SuperCritical()), [Some(vle2), None] => { - let cp = State::critical_point_binary(eos, tp, SolverOptions::default())?; + let cp = + State::critical_point_binary(eos, tp, None, None, SolverOptions::default())?; let cp_vle = PhaseEquilibrium::from_states(cp.clone(), cp.clone()); ([0.0, cp.molefracs[0]], (vle2, cp_vle), bubble) } [None, Some(vle1)] => { - let cp = State::critical_point_binary(eos, tp, SolverOptions::default())?; + let cp = + State::critical_point_binary(eos, tp, None, None, SolverOptions::default())?; let cp_vle = PhaseEquilibrium::from_states(cp.clone(), cp.clone()); ([1.0, cp.molefracs[0]], (vle1, cp_vle), bubble) } diff --git a/src/python/state.rs b/src/python/state.rs index 6473a7b..3e16bc7 100644 --- a/src/python/state.rs +++ b/src/python/state.rs @@ -191,6 +191,8 @@ macro_rules! impl_state { /// The equation of state to use. /// temperature: SINumber /// temperature. + /// initial_molefracs: [float], optional + /// An initial guess for the composition. /// max_iter : int, optional /// The maximum number of iterations. /// tol: float, optional @@ -202,10 +204,11 @@ macro_rules! impl_state { /// ------- /// State : State at critical conditions. #[staticmethod] - #[pyo3(text_signature = "(eos, temperature, max_iter=None, tol=None, verbosity=None)")] + #[pyo3(text_signature = "(eos, temperature, initial_molefracs=None, max_iter=None, tol=None, verbosity=None)")] fn critical_point_binary_t( eos: $py_eos, temperature: PySINumber, + initial_molefracs: Option<[f64; 2]>, max_iter: Option, tol: Option, verbosity: Option, @@ -213,6 +216,7 @@ macro_rules! impl_state { Ok(PyState(State::critical_point_binary_t( &eos.0, temperature.into(), + initial_molefracs, (max_iter, tol, verbosity.map(|v| v.0)).into(), )?)) } @@ -226,6 +230,10 @@ macro_rules! impl_state { /// The equation of state to use. /// pressure: SINumber /// pressure. + /// initial_temperature: SINumber, optional + /// The initial temperature. + /// initial_molefracs: [float], optional + /// An initial guess for the composition. /// max_iter : int, optional /// The maximum number of iterations. /// tol: float, optional @@ -237,10 +245,12 @@ macro_rules! impl_state { /// ------- /// State : State at critical conditions. #[staticmethod] - #[pyo3(text_signature = "(eos, temperature, max_iter=None, tol=None, verbosity=None)")] + #[pyo3(text_signature = "(eos, pressure, initial_temperature=None, initial_molefracs=None, max_iter=None, tol=None, verbosity=None)")] fn critical_point_binary_p( eos: $py_eos, pressure: PySINumber, + initial_temperature: Option, + initial_molefracs: Option<[f64; 2]>, max_iter: Option, tol: Option, verbosity: Option, @@ -248,6 +258,8 @@ macro_rules! impl_state { Ok(PyState(State::critical_point_binary_p( &eos.0, pressure.into(), + initial_temperature.map(|t| t.into()), + initial_molefracs, (max_iter, tol, verbosity.map(|v| v.0)).into(), )?)) } diff --git a/src/python/statehd.rs b/src/python/statehd.rs index 51f5172..ab4dc2a 100644 --- a/src/python/statehd.rs +++ b/src/python/statehd.rs @@ -1,7 +1,7 @@ use crate::StateHD; use ndarray::Array1; use num_dual::python::{PyDual3Dual64, PyDual3_64, PyDual64, PyHyperDual64, PyHyperDualDual64}; -use num_dual::{Dual3, Dual3_64, Dual64, HyperDual, HyperDual64}; +use num_dual::{Dual, Dual3, Dual3_64, Dual64, DualVec64, HyperDual, HyperDual64}; use pyo3::prelude::*; use std::convert::From; @@ -35,6 +35,26 @@ impl From>> for PyStateHDD { } } +#[pyclass(name = "StateHDDV2")] +#[derive(Clone)] +pub struct PyStateHDDV2(StateHD, f64>>); + +impl From, f64>>> for PyStateHDDV2 { + fn from(s: StateHD, f64>>) -> Self { + Self(s) + } +} + +#[pyclass(name = "StateHDDV3")] +#[derive(Clone)] +pub struct PyStateHDDV3(StateHD, f64>>); + +impl From, f64>>> for PyStateHDDV3 { + fn from(s: StateHD, f64>>) -> Self { + Self(s) + } +} + #[pyclass(name = "StateD")] #[derive(Clone)] pub struct PyStateD(StateHD); @@ -45,7 +65,17 @@ impl From> for PyStateD { } } -#[pyclass(name = "StateHD3")] +#[pyclass(name = "StateDDV3")] +#[derive(Clone)] +pub struct PyStateDDV3(StateHD, f64>>); + +impl From, f64>>> for PyStateDDV3 { + fn from(s: StateHD, f64>>) -> Self { + Self(s) + } +} + +#[pyclass(name = "StateD3")] #[derive(Clone)] pub struct PyStateD3(StateHD); @@ -55,7 +85,7 @@ impl From> for PyStateD3 { } } -#[pyclass(name = "StateHD3D")] +#[pyclass(name = "StateD3D")] #[derive(Clone)] pub struct PyStateD3D(StateHD>); @@ -65,6 +95,26 @@ impl From>> for PyStateD3D { } } +#[pyclass(name = "StateD3DV2")] +#[derive(Clone)] +pub struct PyStateD3DV2(StateHD, f64>>); + +impl From, f64>>> for PyStateD3DV2 { + fn from(s: StateHD, f64>>) -> Self { + Self(s) + } +} + +#[pyclass(name = "StateD3DV3")] +#[derive(Clone)] +pub struct PyStateD3DV3(StateHD, f64>>); + +impl From, f64>>> for PyStateD3DV3 { + fn from(s: StateHD, f64>>) -> Self { + Self(s) + } +} + macro_rules! impl_state_hd { ($pystate:ty, $pyhd:ty, $state:ty, $hd:ty) => { #[pymethods] diff --git a/src/python/user_defined.rs b/src/python/user_defined.rs index 5f59806..d773635 100644 --- a/src/python/user_defined.rs +++ b/src/python/user_defined.rs @@ -2,7 +2,7 @@ use crate::python::{statehd::*, PyContributions, PyVerbosity}; use crate::*; use ndarray::prelude::*; use num_dual::python::{PyDual3Dual64, PyDual3_64, PyDual64, PyHyperDual64, PyHyperDualDual64}; -use num_dual::{Dual3, Dual3_64, Dual64, HyperDual, HyperDual64}; +use num_dual::{Dual, Dual3, Dual3_64, Dual64, DualVec64, HyperDual, HyperDual64}; use numpy::convert::{IntoPyArray, ToPyArray}; use numpy::{PyArray1, PyArray2}; use pyo3::exceptions::PyValueError; @@ -183,11 +183,74 @@ impl fmt::Display for PyHelmholtzEnergy { } } +#[pyclass(name = "DualDualVec64_2")] +#[derive(Clone)] +pub struct PyDualDualVec64_3(Dual, f64>); + +impl From for Dual, f64> { + fn from(value: PyDualDualVec64_3) -> Self { + value.0 + } +} + +#[pyclass(name = "HyperDualDualVec64_2")] +#[derive(Clone)] +pub struct PyHyperDualDualVec64_2(HyperDual, f64>); + +impl From for HyperDual, f64> { + fn from(value: PyHyperDualDualVec64_2) -> Self { + value.0 + } +} + +#[pyclass(name = "HyperDualDualVec64_3")] +#[derive(Clone)] +pub struct PyHyperDualDualVec64_3(HyperDual, f64>); + +impl From for HyperDual, f64> { + fn from(value: PyHyperDualDualVec64_3) -> Self { + value.0 + } +} + +#[pyclass(name = "Dual3DualVec64_2")] +#[derive(Clone)] +pub struct PyDual3DualVec64_2(Dual3, f64>); + +impl From for Dual3, f64> { + fn from(value: PyDual3DualVec64_2) -> Self { + value.0 + } +} + +#[pyclass(name = "Dual3DualVec64_3")] +#[derive(Clone)] +pub struct PyDual3DualVec64_3(Dual3, f64>); + +impl From for Dual3, f64> { + fn from(value: PyDual3DualVec64_3) -> Self { + value.0 + } +} + impl_helmholtz_energy!(PyStateD, PyDual64, Dual64); +impl_helmholtz_energy!(PyStateDDV3, PyDualDualVec64_3, Dual, f64>); impl_helmholtz_energy!(PyStateHD, PyHyperDual64, HyperDual64); impl_helmholtz_energy!(PyStateHDD, PyHyperDualDual64, HyperDual); +impl_helmholtz_energy!( + PyStateHDDV2, + PyHyperDualDualVec64_2, + HyperDual, f64> +); +impl_helmholtz_energy!( + PyStateHDDV3, + PyHyperDualDualVec64_3, + HyperDual, f64> +); impl_helmholtz_energy!(PyStateD3, PyDual3_64, Dual3_64); impl_helmholtz_energy!(PyStateD3D, PyDual3Dual64, Dual3); +impl_helmholtz_energy!(PyStateD3DV2, PyDual3DualVec64_2, Dual3, f64>); +impl_helmholtz_energy!(PyStateD3DV3, PyDual3DualVec64_3, Dual3, f64>); impl_helmholtz_energy!(PyStateF, f64, f64); impl_equation_of_state!(PyUserDefinedEos); @@ -200,10 +263,15 @@ impl_vle_state!(PyEoSObj, PyUserDefinedEos); pub fn user_defined(_py: Python, m: &PyModule) -> PyResult<()> { m.add_class::()?; m.add_class::()?; + m.add_class::()?; m.add_class::()?; m.add_class::()?; m.add_class::()?; + m.add_class::()?; + m.add_class::()?; m.add_class::()?; + m.add_class::()?; + m.add_class::()?; m.add_class::()?; m.add_class::()?; m.add_class::()?; diff --git a/src/state/critical_point.rs b/src/state/critical_point.rs index 25abcf1..347d4e3 100644 --- a/src/state/critical_point.rs +++ b/src/state/critical_point.rs @@ -1,13 +1,11 @@ -use super::{Contributions, State, StateHD, TPSpec}; +use super::{State, StateHD, TPSpec}; use crate::equation_of_state::EquationOfState; use crate::errors::{EosError, EosResult}; use crate::phase_equilibria::{SolverOptions, Verbosity}; use crate::EosUnit; -use argmin::prelude::{ArgminOp, Error, Executor}; -use argmin::solver::brent::Brent; use ndarray::{arr1, arr2, Array1, Array2}; use num_dual::linalg::{norm, smallest_ev, LU}; -use num_dual::{Dual3, Dual64, DualNum, HyperDual}; +use num_dual::{Dual, Dual3, Dual64, DualNum, DualVec64, HyperDual, StaticVec}; use num_traits::{One, Zero}; use quantity::{QuantityArray1, QuantityScalar}; use std::rc::Rc; @@ -38,47 +36,28 @@ impl State { .collect() } - /// Calculate the critical point of a binary system for given temperature. - pub fn critical_point_binary_t( - eos: &Rc, - temperature: QuantityScalar, - options: SolverOptions, - ) -> EosResult - where - QuantityScalar: std::fmt::Display, - { - Self::critical_point_binary(eos, TPSpec::Temperature(temperature), options) - } - - /// Calculate the critical point of a binary system for given pressure. - pub fn critical_point_binary_p( - eos: &Rc, - pressure: QuantityScalar, - options: SolverOptions, - ) -> EosResult - where - QuantityScalar: std::fmt::Display, - { - Self::critical_point_binary(eos, TPSpec::Pressure(pressure), options) - } - pub(crate) fn critical_point_binary( eos: &Rc, tp: TPSpec, + initial_temperature: Option>, + initial_molefracs: Option<[f64; 2]>, options: SolverOptions, ) -> EosResult where QuantityScalar: std::fmt::Display, { - let solver = Brent::new(1e-10, 1.0 - 1e-10, options.tol.unwrap_or(TOL_CRIT_POINT)); - let cost = CritOp::new(eos, tp); - let x = Executor::new(cost, solver, 0.5) - .max_iters(options.max_iter.unwrap_or(MAX_ITER_CRIT_POINT) as u64) - .run()? - .state - .best_param; - let moles = arr1(&[x, 1.0 - x]) * U::reference_moles(); - State::critical_point(eos, Some(&moles), None, SolverOptions::default()) + match tp { + TPSpec::Temperature(t) => { + Self::critical_point_binary_t(eos, t, initial_molefracs, options) + } + TPSpec::Pressure(p) => Self::critical_point_binary_p( + eos, + p, + initial_temperature, + initial_molefracs, + options, + ), + } } /// Calculate the critical point of a system for given moles. @@ -194,6 +173,191 @@ impl State { } Err(EosError::NotConverged(String::from("Critical point"))) } + + /// Calculate the critical point of a binary system for given temperature. + pub fn critical_point_binary_t( + eos: &Rc, + temperature: QuantityScalar, + initial_molefracs: Option<[f64; 2]>, + options: SolverOptions, + ) -> EosResult + where + QuantityScalar: std::fmt::Display, + { + let (max_iter, tol, verbosity) = options.unwrap_or(MAX_ITER_CRIT_POINT, TOL_CRIT_POINT); + + let t = temperature.to_reduced(U::reference_temperature())?; + let x = StaticVec::new_vec(initial_molefracs.unwrap_or([0.5, 0.5])); + let max_density = eos + .max_density(Some(&(arr1(x.raw_array()) * U::reference_moles())))? + .to_reduced(U::reference_density())?; + let mut rho = x * 0.3 * max_density; + + log_iter!( + verbosity, + " iter | residual | density 1 | density 2 " + ); + log_iter!(verbosity, "{:-<69}", ""); + log_iter!( + verbosity, + " {:4} | | {:12.8} | {:12.8}", + 0, + rho[0] * U::reference_density(), + rho[1] * U::reference_density(), + ); + + for i in 1..=max_iter { + // calculate residuals and derivatives w.r.t. partial densities + let r = StaticVec::new_vec([DualVec64::from_re(rho[0]), DualVec64::from_re(rho[1])]) + .derive(); + let res = critical_point_objective_t(eos, t, r)?; + + // calculate Newton step + let h = res.jacobian(); + let res = res.map(|r| r.re); + let mut delta = StaticVec::new_vec([ + h[(1, 1)] * res[0] - h[(0, 1)] * res[1], + h[(0, 0)] * res[1] - h[(1, 0)] * res[0], + ]) / (h[(0, 0)] * h[(1, 1)] - h[(0, 1)] * h[(1, 0)]); + + // reduce step if necessary + for i in 0..2 { + if delta[i].abs() > 0.03 * max_density { + delta *= 0.03 * max_density / delta[i].abs() + } + } + + // apply step + rho -= delta; + rho[0] = f64::max(rho[0], 1e-4 * max_density); + rho[1] = f64::max(rho[1], 1e-4 * max_density); + + log_iter!( + verbosity, + " {:4} | {:14.8e} | {:12.8} | {:12.8}", + i, + res.norm(), + rho[0] * U::reference_density(), + rho[1] * U::reference_density(), + ); + + // check convergence + if res.norm() < tol { + log_result!( + verbosity, + "Critical point calculation converged in {} step(s)\n", + i + ); + return State::new_nvt( + eos, + t * U::reference_temperature(), + U::reference_volume(), + &(arr1(rho.raw_array()) * U::reference_moles()), + ); + } + } + Err(EosError::NotConverged(String::from("Critical point"))) + } + + /// Calculate the critical point of a binary system for given pressure. + pub fn critical_point_binary_p( + eos: &Rc, + pressure: QuantityScalar, + initial_temperature: Option>, + initial_molefracs: Option<[f64; 2]>, + options: SolverOptions, + ) -> EosResult + where + QuantityScalar: std::fmt::Display, + { + let (max_iter, tol, verbosity) = options.unwrap_or(MAX_ITER_CRIT_POINT, TOL_CRIT_POINT); + + let p = pressure.to_reduced(U::reference_pressure())?; + let mut t = initial_temperature + .map(|t| t.to_reduced(U::reference_temperature())) + .transpose()? + .unwrap_or(300.0); + let x = StaticVec::new_vec(initial_molefracs.unwrap_or([0.5, 0.5])); + let max_density = eos + .max_density(Some(&(arr1(x.raw_array()) * U::reference_moles())))? + .to_reduced(U::reference_density())?; + let mut rho = x * 0.3 * max_density; + + log_iter!( + verbosity, + " iter | residual | temperature | density 1 | density 2 " + ); + log_iter!(verbosity, "{:-<87}", ""); + log_iter!( + verbosity, + " {:4} | | {:13.8} | {:12.8} | {:12.8}", + 0, + t * U::reference_temperature(), + rho[0] * U::reference_density(), + rho[1] * U::reference_density(), + ); + + for i in 1..=max_iter { + // calculate residuals and derivatives w.r.t. temperature and partial densities + let x = StaticVec::new_vec([ + DualVec64::from_re(t), + DualVec64::from_re(rho[0]), + DualVec64::from_re(rho[1]), + ]) + .derive(); + let r = StaticVec::new_vec([x[1], x[2]]); + let res = critical_point_objective_p(eos, p, x[0], r)?; + + // calculate Newton step + let h = arr2(res.jacobian().raw_data()); + let res = arr1(res.map(|r| r.re).raw_array()); + let mut delta = LU::new(h)?.solve(&res); + + // reduce step if necessary + if delta[0].abs() > 0.25 * t { + delta *= 0.25 * t / delta[0].abs() + } + if delta[1].abs() > 0.03 * max_density { + delta *= 0.03 * max_density / delta[1].abs() + } + if delta[2].abs() > 0.03 * max_density { + delta *= 0.03 * max_density / delta[2].abs() + } + + // apply step + t -= delta[0]; + rho[0] -= delta[1]; + rho[1] -= delta[2]; + rho[0] = f64::max(rho[0], 1e-4 * max_density); + rho[1] = f64::max(rho[1], 1e-4 * max_density); + + log_iter!( + verbosity, + " {:4} | {:14.8e} | {:13.8} | {:12.8} | {:12.8}", + i, + norm(&res), + t * U::reference_temperature(), + rho[0] * U::reference_density(), + rho[1] * U::reference_density(), + ); + + // check convergence + if norm(&res) < tol { + log_result!( + verbosity, + "Critical point calculation converged in {} step(s)\n", + i + ); + return State::new_nvt( + eos, + t * U::reference_temperature(), + U::reference_volume(), + &(arr1(rho.raw_array()) * U::reference_moles()), + ); + } + } + Err(EosError::NotConverged(String::from("Critical point"))) + } } pub fn critical_point_objective( @@ -236,40 +400,84 @@ pub fn critical_point_objective( Ok(arr1(&[eval, res.v3])) } -struct CritOp { - eos: Rc, - tp: TPSpec, -} +fn critical_point_objective_t( + eos: &Rc, + temperature: f64, + density: StaticVec, 2>, +) -> EosResult, 2>> { + // calculate second partial derivatives w.r.t. moles + let t = HyperDual::from(temperature); + let v = HyperDual::from(1.0); + let qij = Array2::from_shape_fn((eos.components(), eos.components()), |(i, j)| { + let mut m = density.map(HyperDual::from_re); + m[i].eps1[0] = DualVec64::one(); + m[j].eps2[0] = DualVec64::one(); + let state = StateHD::new(t, v, arr1(&[m[0], m[1]])); + (eos.evaluate_residual(&state).eps1eps2[(0, 0)] + + eos.ideal_gas().evaluate(&state).eps1eps2[(0, 0)]) + * (density[i] * density[j]).sqrt() + }); -impl CritOp { - fn new(eos: &Rc, tp: TPSpec) -> Self { - Self { - eos: eos.clone(), - tp, - } - } + // calculate smallest eigenvalue and corresponding eigenvector of q + let (eval, evec) = smallest_ev(qij); + + // evaluate third partial derivative w.r.t. s + let moles_hd = Array1::from_shape_fn(eos.components(), |i| { + Dual3::new( + density[i], + evec[i] * density[i].sqrt(), + DualVec64::zero(), + DualVec64::zero(), + ) + }); + let state_s = StateHD::new(Dual3::from(temperature), Dual3::from(1.0), moles_hd); + let res = eos.evaluate_residual(&state_s) + eos.ideal_gas().evaluate(&state_s); + Ok(StaticVec::new_vec([eval, res.v3])) } -impl ArgminOp for CritOp -where - QuantityScalar: std::fmt::Display, -{ - type Param = f64; - type Output = f64; - type Jacobian = (); - type Hessian = (); - type Float = f64; - - fn apply(&self, p: &Self::Param) -> Result { - let moles = arr1(&[*p, 1.0 - *p]) * U::reference_moles(); - let state = State::critical_point(&self.eos, Some(&moles), None, SolverOptions::default())?; - match self.tp { - TPSpec::Pressure(p) => Ok( - (state.pressure(Contributions::Total) - p).to_reduced(U::reference_pressure())? - ), - TPSpec::Temperature(t) => { - Ok((state.temperature - t).to_reduced(U::reference_temperature())?) - } - } - } +fn critical_point_objective_p( + eos: &Rc, + pressure: f64, + temperature: DualVec64<3>, + density: StaticVec, 2>, +) -> EosResult, 3>> { + // calculate second partial derivatives w.r.t. moles + let t = HyperDual::from_re(temperature); + let v = HyperDual::from(1.0); + let qij = Array2::from_shape_fn((eos.components(), eos.components()), |(i, j)| { + let mut m = density.map(HyperDual::from_re); + m[i].eps1[0] = DualVec64::one(); + m[j].eps2[0] = DualVec64::one(); + let state = StateHD::new(t, v, arr1(&[m[0], m[1]])); + (eos.evaluate_residual(&state).eps1eps2[(0, 0)] + + eos.ideal_gas().evaluate(&state).eps1eps2[(0, 0)]) + * (density[i] * density[j]).sqrt() + }); + + // calculate smallest eigenvalue and corresponding eigenvector of q + let (eval, evec) = smallest_ev(qij); + + // evaluate third partial derivative w.r.t. s + let moles_hd = Array1::from_shape_fn(eos.components(), |i| { + Dual3::new( + density[i], + evec[i] * density[i].sqrt(), + DualVec64::zero(), + DualVec64::zero(), + ) + }); + let state_s = StateHD::new(Dual3::from_re(temperature), Dual3::from(1.0), moles_hd); + let res = eos.evaluate_residual(&state_s) + eos.ideal_gas().evaluate(&state_s); + + // calculate pressure + let v = Dual::from(1.0).derive(); + let m = arr1(&[Dual::from_re(density[0]), Dual::from_re(density[1])]); + let state_p = StateHD::new(Dual::from_re(temperature), v, m); + let p = eos.evaluate_residual(&state_p) + eos.ideal_gas().evaluate(&state_p); + + Ok(StaticVec::new_vec([ + eval, + res.v3, + p.eps[0] * temperature + pressure, + ])) }