Remove unused impls; updated Python interface

feos-org · g-bauer · Oct 18, 2022 · Sep 27, 2022 · Sep 28, 2022 · Sep 30, 2022
commit 51715e2fae6a7a15003f7843305993e4edaf3961
diff --git a/feos-core/src/lib.rs b/feos-core/src/lib.rs
@@ -55,24 +55,6 @@ pub fn initialize_global_thread_pool(num_threads: usize) -> Result<(), String> {
         .map_err(|e| e.to_string())
 }
 
-/// Initialize global thread pool.
-///
-/// Parameters
-/// ----------
-/// num_threads : int
-///     number of threads to use
-#[cfg(feature = "python")]
-#[cfg_attr(
-    all(feature = "python", feature = "rayon"),
-    pyo3::pyfunction,
-    pyo3(name = "initialize_global_thread_pool"),
-    pyo3(text_signature = "(num_threads)")
-)]
-pub fn initialize_global_thread_pool_py(num_threads: usize) -> pyo3::PyResult<()> {
-    initialize_global_thread_pool(num_threads)
-        .map_err(|e| pyo3::PyErr::new::<pyo3::exceptions::PyTypeError, _>(e))
-}
-
 /// Consistent conversions between quantities and reduced properties.
 pub trait EosUnit: Unit + Send + Sync {
     fn reference_temperature() -> QuantityScalar<Self>;

diff --git a/feos-core/src/phase_equilibria/phase_diagram_pure.rs b/feos-core/src/phase_equilibria/phase_diagram_pure.rs
@@ -3,11 +3,12 @@ use crate::equation_of_state::EquationOfState;
 use crate::errors::EosResult;
 use crate::state::{State, StateVec};
 use crate::EosUnit;
+#[cfg(feature = "rayon")]
+use ndarray::{Array1, ArrayView1, Axis};
 use quantity::{QuantityArray1, QuantityScalar};
 #[cfg(feature = "rayon")]
 use rayon_::prelude::*;
-#[cfg(feature = "rayon")]
-use ndarray::{Array1, ArrayView1, Axis};
+use rayon_::ThreadPool;
 use std::sync::Arc;
 
 /// Pure component and binary mixture phase diagrams.
@@ -68,17 +69,17 @@ impl<U: EosUnit, E: EquationOfState> PhaseDiagram<U, E> {
 
 #[cfg(feature = "rayon")]
 impl<U: EosUnit, E: EquationOfState> PhaseDiagram<U, E> {
-    fn solve_range(
+    fn solve_temperatures(
         eos: &Arc<E>,
-        range: ArrayView1<f64>,
+        temperatures: ArrayView1<f64>,
         options: SolverOptions,
     ) -> EosResult<Vec<PhaseEquilibrium<U, E, 2>>>
     where
         QuantityScalar<U>: std::fmt::Display + std::fmt::LowerExp,
     {
-        let mut states = Vec::with_capacity(range.len());
+        let mut states = Vec::with_capacity(temperatures.len());
         let mut vle = None;
-        for ti in range {
+        for ti in temperatures {
             vle = PhaseEquilibrium::pure(
                 eos,
                 *ti * U::reference_temperature(),
@@ -97,8 +98,9 @@ impl<U: EosUnit, E: EquationOfState> PhaseDiagram<U, E> {
         eos: &Arc<E>,
         min_temperature: QuantityScalar<U>,
         npoints: usize,
-        critical_temperature: Option<QuantityScalar<U>>,
         chunksize: usize,
+        thread_pool: ThreadPool,
+        critical_temperature: Option<QuantityScalar<U>>,
         options: SolverOptions,
     ) -> EosResult<Self>
     where
@@ -114,23 +116,16 @@ impl<U: EosUnit, E: EquationOfState> PhaseDiagram<U, E> {
             npoints - 1,
         );
 
-        let mut states: Vec<PhaseEquilibrium<U, E, 2>> = temperatures
-            .axis_chunks_iter(Axis(0), chunksize)
-            .into_par_iter()
-            .filter_map(|t| Self::solve_range(eos, t, options).ok())
-            .flatten()
-            .collect();
+        let mut states: Vec<PhaseEquilibrium<U, E, 2>> = thread_pool.install(|| {
+            temperatures
+                .axis_chunks_iter(Axis(0), chunksize)
+                .into_par_iter()
+                .filter_map(|t| Self::solve_temperatures(eos, t, options).ok())
+                .flatten()
+                .collect()
+        });
 
         states.push(PhaseEquilibrium::from_states(sc.clone(), sc));
-        // let mut vle = None;
-        // for ti in temperatures.into_raw_vec().par_iter() {
-        //     vle = PhaseEquilibrium::pure(eos, ti, vle.as_ref(), options).ok();
-        //     if let Some(vle) = vle.as_ref() {
-        //         states.push(vle.clone());
-        //     }
-        // }
-        // states.push(PhaseEquilibrium::from_states(sc.clone(), sc));
-
         Ok(PhaseDiagram { states })
     }
 }
diff --git a/feos-core/src/python/phase_equilibria.rs b/feos-core/src/python/phase_equilibria.rs
@@ -513,25 +513,60 @@ macro_rules! impl_phase_equilibrium {
                 Ok(Self(dia))
             }
 
+            /// Calculate a pure component phase diagram in parallel.
+            ///
+            /// Parameters
+            /// ----------
+            /// eos : Eos
+            ///     The equation of state.
+            /// min_temperature: SINumber
+            ///     The lower limit for the temperature.
+            /// npoints : int
+            ///     The number of points.
+            /// chunksize : int
+            ///     The number of points that are calculated in sequence
+            ///     within a thread.
+            /// nthreads : int
+            ///     Number of threads.
+            /// critical_temperature: SINumber, optional
+            ///     An estimate for the critical temperature to initialize
+            ///     the calculation if necessary. For most components not necessary.
+            ///     Defaults to `None`.
+            /// max_iter : int, optional
+            ///     The maximum number of iterations.
+            /// tol: float, optional
+            ///     The solution tolerance.
+            /// verbosity : Verbosity, optional
+            ///     The verbosity.
+            ///
+            /// Returns
+            /// -------
+            /// PhaseDiagram 
             #[cfg(feature = "rayon")]
             #[staticmethod]
-            #[pyo3(text_signature = "(eos, min_temperature, npoints, critical_temperature=None, max_iter=None, tol=None, verbosity=None)")]
+            #[pyo3(text_signature = "(eos, min_temperature, npoints, chunksize, nthreads, critical_temperature=None, max_iter=None, tol=None, verbosity=None)")]
             pub fn par_pure(
                 eos: &$py_eos,
                 min_temperature: PySINumber,
                 npoints: usize,
                 chunksize: usize,
+                nthreads: usize,
                 critical_temperature: Option<PySINumber>,
                 max_iter: Option<usize>,
                 tol: Option<f64>,
                 verbosity: Option<Verbosity>,
             ) -> PyResult<Self> {
+                let thread_pool = rayon_::ThreadPoolBuilder::new()
+                    .num_threads(nthreads)
+                    .build()
+                    .map_err(|e| pyo3::exceptions::PyValueError::new_err(e.to_string()))?;
                 let dia = PhaseDiagram::par_pure(
                     &eos.0,
                     min_temperature.into(),
                     npoints,
-                    critical_temperature.map(|t| t.into()),
                     chunksize,
+                    thread_pool,
+                    critical_temperature.map(|t| t.into()),
                     (max_iter, tol, verbosity).into(),
                 )?;
                 Ok(Self(dia))

diff --git a/src/estimator/diffusion.rs b/src/estimator/diffusion.rs
@@ -2,8 +2,6 @@ use super::{DataSet, EstimatorError};
 use feos_core::{DensityInitialization, EntropyScaling, EosUnit, EquationOfState, State};
 use ndarray::{arr1, Array1};
 use quantity::{QuantityArray1, QuantityScalar};
-#[cfg(feature = "rayon")]
-use rayon_::prelude::*;
 use std::collections::HashMap;
 use std::sync::Arc;
 

diff --git a/src/estimator/liquid_density.rs b/src/estimator/liquid_density.rs
@@ -5,8 +5,6 @@ use feos_core::{
 };
 use ndarray::arr1;
 use quantity::{QuantityArray1, QuantityScalar};
-#[cfg(feature = "rayon")]
-use rayon_::prelude::*;
 use std::collections::HashMap;
 use std::sync::Arc;
 

diff --git a/src/estimator/surface_tension.rs b/src/estimator/surface_tension.rs
@@ -0,0 +1,116 @@
+use super::{DataSet, EstimatorError};
+use feos_core::{Contributions, EosUnit, EquationOfState, PhaseEquilibrium, SolverOptions, State};
+use feos_dft::{HelmholtzEnergyFunctional, DFT};
+use feos_dft::{DFTProfile, DFTSolver, interface::PlanarInterface};
+use ndarray::Array1;
+use quantity::{QuantityArray1, QuantityScalar};
+use std::collections::HashMap;
+use std::rc::Rc;
+
+/// Store experimental surface tension data.
+#[derive(Clone)]
+pub struct SurfaceTension<U: EosUnit> {
+    pub target: QuantityArray1<U>,
+    temperature: QuantityArray1<U>,
+    max_temperature: QuantityScalar<U>,
+    datapoints: usize,
+    grid: usize,
+    width: QuantityScalar<U>,
+    solver_options: SolverOptions,
+}
+
+impl<U: EosUnit> SurfaceTension<U> {
+    /// Create a new data set for surface tension.
+    ///
+    /// If the equation of state fails to compute the surface tension
+    /// (e.g. when it underestimates the critical point) the vapor
+    /// pressure can be estimated.
+    /// If `extrapolate` is `true`, the surface tension is estimated by
+    /// calculating the slope of ln(p) over 1/T.
+    /// If `extrapolate` is `false`, it is set to `NAN`.
+    pub fn new(
+        target: QuantityArray1<U>,
+        temperature: QuantityArray1<U>,
+        grid: usize,
+        width: QuantityScalar<U>,
+        solver_options: Option<SolverOptions>,
+    ) -> Result<Self, EstimatorError> {
+        let datapoints = target.len();
+        let max_temperature = temperature
+            .to_reduced(U::reference_temperature())?
+            .into_iter()
+            .reduce(|a, b| a.max(b))
+            .unwrap()
+            * U::reference_temperature();
+        Ok(Self {
+            target,
+            temperature,
+            max_temperature,
+            datapoints,
+            grid,
+            width,
+            solver_options: solver_options.unwrap_or_default(),
+        })
+    }
+
+    /// Return temperature.
+    pub fn temperature(&self) -> QuantityArray1<U> {
+        self.temperature.clone()
+    }
+}
+
+impl<U: EosUnit, M: DFT<HelmholtzEnergyFunctional>> DataSet<U, M> for SurfaceTension<U> {
+    fn target(&self) -> &QuantityArray1<U> {
+        &self.target
+    }
+
+    fn target_str(&self) -> &str {
+        "surface tension"
+    }
+
+    fn input_str(&self) -> Vec<&str> {
+        vec!["temperature"]
+    }
+
+    fn predict(&self, eos: &Rc<M>) -> Result<QuantityArray1<U>, EstimatorError>
+    where
+        QuantityScalar<U>: std::fmt::Display + std::fmt::LowerExp,
+    {
+        let critical_point =
+            State::critical_point(eos, None, Some(self.max_temperature), self.solver_options)?;
+        let tc = critical_point.temperature;
+        let pc = critical_point.pressure(Contributions::Total);
+
+        let t0 = 0.9 * tc;
+        let p0 = PhaseEquilibrium::pure(eos, t0, None, self.solver_options)?
+            .vapor()
+            .pressure(Contributions::Total);
+
+        let b = pc.to_reduced(p0)?.ln() / (1.0 / tc - 1.0 / t0);
+        let a = pc.to_reduced(U::reference_pressure())?.ln() - b.to_reduced(tc)?;
+
+        let unit = self.target.get(0);
+        let mut prediction = Array1::zeros(self.datapoints) * unit;
+        for i in 0..self.datapoints {
+            let t = self.temperature.get(i);
+            let vle = PhaseEquilibrium::pure(&Rc::new(eos.into()), t, None, SolverOptions::default());
+            let gamma = vle
+                .and_then(|vle| PlanarInterface::from_tanh(&vle, self.grid, self.width, tc))
+                .and_then(|interface| interface.solve(None))
+                .and_then(|interface| interface.surface_tension());
+
+            if let Ok(gamma) = gamma {
+                prediction.try_set(i, gamma)?
+            } else {
+                prediction.try_set(i, f64::NAN * U::reference_surface_tension())?
+            }
+        }
+        Ok(prediction)
+    }
+
+    fn get_input(&self) -> HashMap<String, QuantityArray1<U>> {
+        let mut m = HashMap::with_capacity(1);
+        m.insert("temperature".to_owned(), self.temperature());
+        m
+    }
+}
diff --git a/src/estimator/thermal_conductivity.rs b/src/estimator/thermal_conductivity.rs
@@ -2,8 +2,6 @@ use super::{DataSet, EstimatorError};
 use feos_core::{DensityInitialization, EntropyScaling, EosUnit, EquationOfState, State};
 use ndarray::{arr1, Array1};
 use quantity::{QuantityArray1, QuantityScalar};
-#[cfg(feature = "rayon")]
-use rayon_::prelude::*;
 use std::collections::HashMap;
 use std::sync::Arc;
 

diff --git a/src/estimator/viscosity.rs b/src/estimator/viscosity.rs
@@ -2,8 +2,6 @@ use super::{DataSet, EstimatorError};
 use feos_core::{DensityInitialization, EntropyScaling, EosUnit, EquationOfState, State};
 use ndarray::arr1;
 use quantity::{QuantityArray1, QuantityScalar};
-#[cfg(feature = "rayon")]
-use rayon_::prelude::*;
 use std::collections::HashMap;
 use std::sync::Arc;
 

diff --git a/src/python/mod.rs b/src/python/mod.rs
@@ -25,9 +25,6 @@ pub fn feos(py: Python<'_>, m: &PyModule) -> PyResult<()> {
     m.add("__version__", env!("CARGO_PKG_VERSION"))?;
     m.add_wrapped(wrap_pymodule!(quantity))?;
 
-    #[cfg(feature = "rayon")]
-    m.add_wrapped(wrap_pyfunction!(feos_core::initialize_global_thread_pool_py))?;
-
     m.add_wrapped(wrap_pymodule!(eos))?;
     #[cfg(feature = "dft")]
     m.add_wrapped(wrap_pymodule!(dft))?;