feos-org · g-bauer · Mar 10, 2022 · Feb 28, 2022 · Mar 4, 2022 · Mar 4, 2022
diff --git 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
@@ -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
@@ -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,21 +26,31 @@ pub trait HelmholtzEnergyDual<D: DualNum<f64>> {
 pub trait HelmholtzEnergy:
     HelmholtzEnergyDual<f64>
     + HelmholtzEnergyDual<Dual64>
+    + HelmholtzEnergyDual<Dual<DualVec64<3>, f64>>
     + HelmholtzEnergyDual<HyperDual64>
     + HelmholtzEnergyDual<Dual3_64>
     + HelmholtzEnergyDual<HyperDual<Dual64, f64>>
+    + HelmholtzEnergyDual<HyperDual<DualVec64<2>, f64>>
+    + HelmholtzEnergyDual<HyperDual<DualVec64<3>, f64>>
     + HelmholtzEnergyDual<Dual3<Dual64, f64>>
+    + HelmholtzEnergyDual<Dual3<DualVec64<2>, f64>>
+    + HelmholtzEnergyDual<Dual3<DualVec64<3>, f64>>
     + fmt::Display
 {
 }
 
 impl<T> HelmholtzEnergy for T where
     T: HelmholtzEnergyDual<f64>
         + HelmholtzEnergyDual<Dual64>
+        + HelmholtzEnergyDual<Dual<DualVec64<3>, f64>>
         + HelmholtzEnergyDual<HyperDual64>
         + HelmholtzEnergyDual<Dual3_64>
         + HelmholtzEnergyDual<HyperDual<Dual64, f64>>
+        + HelmholtzEnergyDual<HyperDual<DualVec64<2>, f64>>
+        + HelmholtzEnergyDual<HyperDual<DualVec64<3>, f64>>
         + HelmholtzEnergyDual<Dual3<Dual64, f64>>
+        + HelmholtzEnergyDual<Dual3<DualVec64<2>, f64>>
+        + HelmholtzEnergyDual<Dual3<DualVec64<3>, f64>>
         + fmt::Display
 {
 }
@@ -83,26 +93,36 @@ pub trait IdealGasContributionDual<D: DualNum<f64>> {
 pub trait IdealGasContribution:
     IdealGasContributionDual<f64>
     + IdealGasContributionDual<Dual64>
+    + IdealGasContributionDual<Dual<DualVec64<3>, f64>>
     + IdealGasContributionDual<HyperDual64>
     + IdealGasContributionDual<Dual3_64>
     + IdealGasContributionDual<HyperDual<Dual64, f64>>
+    + IdealGasContributionDual<HyperDual<DualVec64<2>, f64>>
+    + IdealGasContributionDual<HyperDual<DualVec64<3>, f64>>
     + IdealGasContributionDual<Dual3<Dual64, f64>>
+    + IdealGasContributionDual<Dual3<DualVec64<2>, f64>>
+    + IdealGasContributionDual<Dual3<DualVec64<3>, f64>>
     + fmt::Display
 {
 }
 
 impl<T> IdealGasContribution for T where
     T: IdealGasContributionDual<f64>
         + IdealGasContributionDual<Dual64>
+        + IdealGasContributionDual<Dual<DualVec64<3>, f64>>
         + IdealGasContributionDual<HyperDual64>
         + IdealGasContributionDual<Dual3_64>
         + IdealGasContributionDual<HyperDual<Dual64, f64>>
+        + IdealGasContributionDual<HyperDual<DualVec64<2>, f64>>
+        + IdealGasContributionDual<HyperDual<DualVec64<3>, f64>>
         + IdealGasContributionDual<Dual3<Dual64, f64>>
+        + IdealGasContributionDual<Dual3<DualVec64<2>, f64>>
+        + IdealGasContributionDual<Dual3<DualVec64<3>, f64>>
         + fmt::Display
 {
 }
 
-struct DefaultIdealGasContribution();
+struct DefaultIdealGasContribution;
 impl<D: DualNum<f64>> IdealGasContributionDual<D> for DefaultIdealGasContribution {
     fn de_broglie_wavelength(&self, _: D, components: usize) -> Array1<D> {
         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
@@ -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
@@ -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
@@ -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<U: EosUnit, E: EquationOfState> PhaseDiagramBinary<U, E> {
         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
@@ -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,17 +204,19 @@ 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<usize>,
                 tol: Option<f64>,
                 verbosity: Option<PyVerbosity>,
             ) -> PyResult<Self> {
                 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,17 +245,21 @@ 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<PySINumber>,
+                initial_molefracs: Option<[f64; 2]>,
                 max_iter: Option<usize>,
                 tol: Option<f64>,
                 verbosity: Option<PyVerbosity>,
             ) -> PyResult<Self> {
                 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
@@ -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<StateHD<HyperDual<Dual64, f64>>> for PyStateHDD {
     }
 }
 
+#[pyclass(name = "StateHDDV2")]
+#[derive(Clone)]
+pub struct PyStateHDDV2(StateHD<HyperDual<DualVec64<2>, f64>>);
+
+impl From<StateHD<HyperDual<DualVec64<2>, f64>>> for PyStateHDDV2 {
+    fn from(s: StateHD<HyperDual<DualVec64<2>, f64>>) -> Self {
+        Self(s)
+    }
+}
+
+#[pyclass(name = "StateHDDV3")]
+#[derive(Clone)]
+pub struct PyStateHDDV3(StateHD<HyperDual<DualVec64<3>, f64>>);
+
+impl From<StateHD<HyperDual<DualVec64<3>, f64>>> for PyStateHDDV3 {
+    fn from(s: StateHD<HyperDual<DualVec64<3>, f64>>) -> Self {
+        Self(s)
+    }
+}
+
 #[pyclass(name = "StateD")]
 #[derive(Clone)]
 pub struct PyStateD(StateHD<Dual64>);
@@ -45,7 +65,17 @@ impl From<StateHD<Dual64>> for PyStateD {
     }
 }
 
-#[pyclass(name = "StateHD3")]
+#[pyclass(name = "StateDDV3")]
+#[derive(Clone)]
+pub struct PyStateDDV3(StateHD<Dual<DualVec64<3>, f64>>);
+
+impl From<StateHD<Dual<DualVec64<3>, f64>>> for PyStateDDV3 {
+    fn from(s: StateHD<Dual<DualVec64<3>, f64>>) -> Self {
+        Self(s)
+    }
+}
+
+#[pyclass(name = "StateD3")]
 #[derive(Clone)]
 pub struct PyStateD3(StateHD<Dual3_64>);
 
@@ -55,7 +85,7 @@ impl From<StateHD<Dual3_64>> for PyStateD3 {
     }
 }
 
-#[pyclass(name = "StateHD3D")]
+#[pyclass(name = "StateD3D")]
 #[derive(Clone)]
 pub struct PyStateD3D(StateHD<Dual3<Dual64, f64>>);
 
@@ -65,6 +95,26 @@ impl From<StateHD<Dual3<Dual64, f64>>> for PyStateD3D {
     }
 }
 
+#[pyclass(name = "StateD3DV2")]
+#[derive(Clone)]
+pub struct PyStateD3DV2(StateHD<Dual3<DualVec64<2>, f64>>);
+
+impl From<StateHD<Dual3<DualVec64<2>, f64>>> for PyStateD3DV2 {
+    fn from(s: StateHD<Dual3<DualVec64<2>, f64>>) -> Self {
+        Self(s)
+    }
+}
+
+#[pyclass(name = "StateD3DV3")]
+#[derive(Clone)]
+pub struct PyStateD3DV3(StateHD<Dual3<DualVec64<3>, f64>>);
+
+impl From<StateHD<Dual3<DualVec64<3>, f64>>> for PyStateD3DV3 {
+    fn from(s: StateHD<Dual3<DualVec64<3>, f64>>) -> Self {
+        Self(s)
+    }
+}
+
 macro_rules! impl_state_hd {
     ($pystate:ty, $pyhd:ty, $state:ty, $hd:ty) => {
         #[pymethods]