Made State sync + send, added parallel execution to pure PhaseDiagram…

… constructor
feos-org · g-bauer · Oct 18, 2022 · Sep 27, 2022 · Sep 28, 2022 · Sep 30, 2022
commit 47d1a7fbe7f2931be93e2731828299ca56011afe
diff --git a/Cargo.toml b/Cargo.toml
@@ -28,7 +28,7 @@ feos-core = { version = "0.3", path = "feos-core" }
 feos-dft = { version = "0.3", path = "feos-dft", optional = true }
 feos-derive = { version = "0.1", path = "feos-derive" }
 numpy = { version = "0.16", optional = true }
-ndarray = { version = "0.15", features = ["approx"] }
+ndarray = { version = "0.15", features = ["approx", "rayon"] }
 petgraph = { version = "0.6", optional = true }
 thiserror = "1.0"
 conv = "0.3"
@@ -37,6 +37,7 @@ serde = "1.0"
 serde_json = "1.0"
 lazy_static = { version = "1.4", optional = true }
 indexmap = "1.8"
+rayon = "1.5"
 
 [dependencies.pyo3]
 version = "0.16"

diff --git a/feos-core/Cargo.toml b/feos-core/Cargo.toml
@@ -19,7 +19,7 @@ rustdoc-args = [ "--html-in-header", "./docs-header.html" ]
 [dependencies]
 quantity = "0.5"
 num-dual = { version = "0.5", features = ["linalg"] }
-ndarray = { version = "0.15", features = ["serde"] }
+ndarray = { version = "0.15", features = ["serde", "rayon"] }
 num-traits = "0.2"
 thiserror = "1.0"
 serde = { version = "1.0", features = ["derive"] }
@@ -28,6 +28,7 @@ indexmap = "1.8"
 conv = "0.3"
 numpy = { version = "0.16", optional = true }
 pyo3 = { version = "0.16", optional = true }
+rayon = "1.5"
 
 [dev-dependencies]
 approx = "0.4"

diff --git a/feos-core/src/lib.rs b/feos-core/src/lib.rs
@@ -48,7 +48,7 @@ pub use state::{Contributions, DensityInitialization, State, StateBuilder, State
 pub mod python;
 
 /// Consistent conversions between quantities and reduced properties.
-pub trait EosUnit: Unit {
+pub trait EosUnit: Unit + Sync + Send {
     fn reference_temperature() -> QuantityScalar<Self>;
     fn reference_length() -> QuantityScalar<Self>;
     fn reference_density() -> 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,7 +3,10 @@ use crate::equation_of_state::EquationOfState;
 use crate::errors::EosResult;
 use crate::state::{State, StateVec};
 use crate::EosUnit;
+use ndarray::{Array1, ArrayView1};
+use numpy::ndarray::Axis;
 use quantity::{QuantityArray1, QuantityScalar};
+use rayon::prelude::*;
 use std::sync::Arc;
 
 /// Pure component and binary mixture phase diagrams.
@@ -51,6 +54,68 @@ impl<U: EosUnit, E: EquationOfState> PhaseDiagram<U, E> {
         Ok(PhaseDiagram { states })
     }
 
+    fn solve_range(
+        eos: &Arc<E>,
+        range: 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 vle = None;
+        for ti in range {
+            vle = PhaseEquilibrium::pure(eos, *ti * U::reference_temperature(), vle.as_ref(), options).ok();
+            if let Some(vle) = vle.as_ref() {
+                states.push(vle.clone());
+            }
+        }
+        Ok(states)
+    }
+
+    pub fn pure_par(
+        eos: &Arc<E>,
+        min_temperature: QuantityScalar<U>,
+        npoints: usize,
+        critical_temperature: Option<QuantityScalar<U>>,
+        chunksize: usize,
+        options: SolverOptions,
+    ) -> EosResult<Self>
+    where
+        QuantityScalar<U>: std::fmt::Display + std::fmt::LowerExp,
+    {
+        let sc = State::critical_point(eos, None, critical_temperature, SolverOptions::default())?;
+
+        let max_temperature = min_temperature
+            + (sc.temperature - min_temperature) * ((npoints - 2) as f64 / (npoints - 1) as f64);
+        let temperatures = Array1::linspace(
+            min_temperature.to_reduced(U::reference_temperature())?,
+            max_temperature.to_reduced(U::reference_temperature())?,
+            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();
+
+        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 })
+    }
+
     /// Return the vapor states of the diagram.
     pub fn vapor(&self) -> StateVec<'_, U, E> {
         self.states.iter().map(|s| s.vapor()).collect()

diff --git a/feos-core/src/python/cubic.rs b/feos-core/src/python/cubic.rs
@@ -11,7 +11,7 @@ use numpy::PyReadonlyArray2;
 use pyo3::exceptions::PyTypeError;
 use pyo3::prelude::*;
 use std::convert::{TryFrom, TryInto};
-use std::rc::Rc;
+use std::sync::Arc;
 
 /// A pure substance parameter for the Peng-Robinson equation of state.
 #[pyclass(name = "PengRobinsonRecord", unsendable)]
@@ -61,7 +61,7 @@ impl_binary_record!();
 /// PengRobinsonParameters
 #[pyclass(name = "PengRobinsonParameters", unsendable)]
 #[derive(Clone)]
-pub struct PyPengRobinsonParameters(pub Rc<PengRobinsonParameters>);
+pub struct PyPengRobinsonParameters(pub Arc<PengRobinsonParameters>);
 
 impl_parameter!(PengRobinsonParameters, PyPengRobinsonParameters);
 

diff --git a/feos-core/src/python/parameter.rs b/feos-core/src/python/parameter.rs
@@ -570,10 +570,10 @@ macro_rules! impl_parameter {
                             "Could not parse binary input!"
                         )))
                     };
-                    Ok(Self(Rc::new(<$parameter>::from_records(prs, brs.unwrap()))))
+                    Ok(Self(Arc::new(<$parameter>::from_records(prs, brs.unwrap()))))
                 } else {
                     let n = prs.len();
-                    Ok(Self(Rc::new(<$parameter>::from_records(
+                    Ok(Self(Arc::new(<$parameter>::from_records(
                         prs,
                         Array2::from_elem([n, n], <$parameter as Parameter>::Binary::default()),
                     ))))
@@ -589,7 +589,7 @@ macro_rules! impl_parameter {
             #[staticmethod]
             #[pyo3(text_signature = "(pure_record)")]
             fn new_pure(pure_record: PyPureRecord) -> Self {
-                Self(Rc::new(<$parameter>::new_pure(pure_record.0)))
+                Self(Arc::new(<$parameter>::new_pure(pure_record.0)))
             }
 
             /// Creates parameters for a binary system from pure records and an optional
@@ -621,7 +621,7 @@ macro_rules! impl_parameter {
                         }
                     })
                     .transpose()?;
-                Ok(Self(Rc::new(<$parameter>::new_binary(prs, br))))
+                Ok(Self(Arc::new(<$parameter>::new_binary(prs, br))))
             }
 
             /// Creates parameters from json files.
@@ -644,7 +644,7 @@ macro_rules! impl_parameter {
                 binary_path: Option<String>,
                 search_option: Option<IdentifierOption>,
             ) -> Result<Self, ParameterError> {
-                Ok(Self(Rc::new(<$parameter>::from_json(
+                Ok(Self(Arc::new(<$parameter>::from_json(
                     substances,
                     pure_path,
                     binary_path,
@@ -670,7 +670,7 @@ macro_rules! impl_parameter {
                 binary_path: Option<&str>,
                 search_option: Option<IdentifierOption>,
             ) -> Result<Self, ParameterError> {
-                Ok(Self(Rc::new(<$parameter>::from_multiple_json(
+                Ok(Self(Arc::new(<$parameter>::from_multiple_json(
                     &input,
                     binary_path,
                     search_option.unwrap_or(IdentifierOption::Name),
@@ -713,7 +713,7 @@ macro_rules! impl_parameter_from_segments {
                 segment_records: Vec<PySegmentRecord>,
                 binary_segment_records: Option<Vec<PyBinarySegmentRecord>>,
             ) -> Result<Self, ParameterError> {
-                Ok(Self(Rc::new(<$parameter>::from_segments(
+                Ok(Self(Arc::new(<$parameter>::from_segments(
                     chemical_records.into_iter().map(|cr| cr.0).collect(),
                     segment_records.into_iter().map(|sr| sr.0).collect(),
                     binary_segment_records.map(|r| r.into_iter().map(|r| BinaryRecord{id1:r.0.id1,id2:r.0.id2,model_record:r.0.model_record.into()}).collect()),
@@ -745,7 +745,7 @@ macro_rules! impl_parameter_from_segments {
                 binary_path: Option<String>,
                 search_option: Option<IdentifierOption>,
             ) -> Result<Self, ParameterError> {
-                Ok(Self(Rc::new(<$parameter>::from_json_segments(
+                Ok(Self(Arc::new(<$parameter>::from_json_segments(
                     &substances,
                     pure_path,
                     segments_path,

diff --git a/feos-core/src/python/phase_equilibria.rs b/feos-core/src/python/phase_equilibria.rs
@@ -513,6 +513,29 @@ macro_rules! impl_phase_equilibrium {
                 Ok(Self(dia))
             }
 
+            #[staticmethod]
+            #[pyo3(text_signature = "(eos, min_temperature, npoints, critical_temperature=None, max_iter=None, tol=None, verbosity=None)")]
+            pub fn pure_par(
+                eos: &$py_eos,
+                min_temperature: PySINumber,
+                npoints: usize,
+                chunksize: usize,
+                critical_temperature: Option<PySINumber>,
+                max_iter: Option<usize>,
+                tol: Option<f64>,
+                verbosity: Option<Verbosity>,
+            ) -> PyResult<Self> {
+                let dia = PhaseDiagram::pure_par(
+                    &eos.0,
+                    min_temperature.into(),
+                    npoints,
+                    critical_temperature.map(|t| t.into()),
+                    chunksize,
+                    (max_iter, tol, verbosity).into(),
+                )?;
+                Ok(Self(dia))
+            }
+
             /// Calculate the bubble point line of a mixture with given composition.
             ///
             /// In the resulting phase diagram, the liquid states correspond to the

diff --git a/feos-core/src/state/mod.rs b/feos-core/src/state/mod.rs
@@ -15,10 +15,9 @@ use ndarray::prelude::*;
 use num_dual::linalg::{norm, LU};
 use num_dual::*;
 use quantity::{QuantityArray1, QuantityScalar};
-use std::cell::RefCell;
 use std::convert::TryFrom;
 use std::fmt;
-use std::sync::Arc;
+use std::sync::{Arc, Mutex};
 
 mod builder;
 mod cache;
@@ -145,7 +144,7 @@ pub struct State<U, E> {
     /// Reduced moles
     reduced_moles: Array1<f64>,
     /// Cache
-    cache: RefCell<Cache>,
+    cache: Mutex<Cache>,
 }
 
 impl<U: Clone, E> Clone for State<U, E> {
@@ -162,7 +161,7 @@ impl<U: Clone, E> Clone for State<U, E> {
             reduced_temperature: self.reduced_temperature,
             reduced_volume: self.reduced_volume,
             reduced_moles: self.reduced_moles.clone(),
-            cache: self.cache.clone(),
+            cache: Mutex::new(self.cache.lock().unwrap().clone()),
         }
     }
 }
@@ -262,7 +261,7 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
             reduced_temperature: t,
             reduced_volume: v,
             reduced_moles: m,
-            cache: RefCell::new(Cache::with_capacity(eos.components())),
+            cache: Mutex::new(Cache::with_capacity(eos.components())),
         }
     }
 

diff --git a/feos-core/src/state/properties.rs b/feos-core/src/state/properties.rs
@@ -64,7 +64,7 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
             };
         }
 
-        let mut cache = self.cache.borrow_mut();
+        let mut cache = self.cache.lock().unwrap();
 
         let residual = match evaluate {
             Evaluate::IdealGas => None,

diff --git a/feos-dft/src/adsorption/mod.rs b/feos-dft/src/adsorption/mod.rs
@@ -8,7 +8,7 @@ use feos_core::{
 use ndarray::{Array1, Dimension, Ix1, Ix3, RemoveAxis};
 use quantity::{QuantityArray1, QuantityArray2, QuantityScalar};
 use std::iter;
-use std::rc::Rc;
+use std::sync::Arc;
 
 mod external_potential;
 #[cfg(feature = "3d_dft")]
@@ -49,7 +49,7 @@ where
     <D::Larger as Dimension>::Larger: Dimension<Smaller = D::Larger>,
 {
     fn new<S: PoreSpecification<U, D>>(
-        functional: &Rc<DFT<F>>,
+        functional: &Arc<DFT<F>>,
         pore: &S,
         profiles: Vec<EosResult<PoreProfile<U, D, F>>>,
     ) -> Self {
@@ -62,7 +62,7 @@ where
 
     /// Calculate an adsorption isotherm (starting at low pressure)
     pub fn adsorption_isotherm<S: PoreSpecification<U, D>>(
-        functional: &Rc<DFT<F>>,
+        functional: &Arc<DFT<F>>,
         temperature: QuantityScalar<U>,
         pressure: &QuantityArray1<U>,
         pore: &S,
@@ -82,7 +82,7 @@ where
 
     /// Calculate an desorption isotherm (starting at high pressure)
     pub fn desorption_isotherm<S: PoreSpecification<U, D>>(
-        functional: &Rc<DFT<F>>,
+        functional: &Arc<DFT<F>>,
         temperature: QuantityScalar<U>,
         pressure: &QuantityArray1<U>,
         pore: &S,
@@ -108,7 +108,7 @@ where
 
     /// Calculate an equilibrium isotherm
     pub fn equilibrium_isotherm<S: PoreSpecification<U, D>>(
-        functional: &Rc<DFT<F>>,
+        functional: &Arc<DFT<F>>,
         temperature: QuantityScalar<U>,
         pressure: &QuantityArray1<U>,
         pore: &S,
@@ -195,7 +195,7 @@ where
     }
 
     fn isotherm<S: PoreSpecification<U, D>>(
-        functional: &Rc<DFT<F>>,
+        functional: &Arc<DFT<F>>,
         temperature: QuantityScalar<U>,
         pressure: &QuantityArray1<U>,
         pore: &S,
@@ -257,7 +257,7 @@ where
 
     /// Calculate the phase transition from an empty to a filled pore.
     pub fn phase_equilibrium<S: PoreSpecification<U, D>>(
-        functional: &Rc<DFT<F>>,
+        functional: &Arc<DFT<F>>,
         temperature: QuantityScalar<U>,
         p_min: QuantityScalar<U>,
         p_max: QuantityScalar<U>,

diff --git a/feos-dft/src/adsorption/pore.rs b/feos-dft/src/adsorption/pore.rs
@@ -10,7 +10,7 @@ use ndarray::prelude::*;
 use ndarray::Axis as Axis_nd;
 use ndarray::RemoveAxis;
 use quantity::{QuantityArray, QuantityArray2, QuantityScalar};
-use std::rc::Rc;
+use std::sync::Arc;
 
 const POTENTIAL_OFFSET: f64 = 2.0;
 const DEFAULT_GRID_POINTS: usize = 2048;
@@ -60,7 +60,7 @@ pub trait PoreSpecification<U: EosUnit, D: Dimension> {
     where
         D::Larger: Dimension<Smaller = D>,
     {
-        let bulk = StateBuilder::new(&Rc::new(Helium::new()))
+        let bulk = StateBuilder::new(&Arc::new(Helium::new()))
             .temperature(298.0 * U::reference_temperature())
             .density(U::reference_density())
             .build()?;

diff --git a/feos-dft/src/functional.rs b/feos-dft/src/functional.rs
@@ -132,7 +132,7 @@ pub enum MoleculeShape<'a> {
 }
 
 /// A general Helmholtz energy functional.
-pub trait HelmholtzEnergyFunctional: Sized {
+pub trait HelmholtzEnergyFunctional: Sized + Sync + Send {
     /// Return a slice of [FunctionalContribution]s.
     fn contributions(&self) -> &[Box<dyn FunctionalContribution>];
 

diff --git a/feos-dft/src/functional_contribution.rs b/feos-dft/src/functional_contribution.rs
@@ -84,6 +84,8 @@ pub trait FunctionalContribution:
     + FunctionalContributionDual<Dual3<DualVec64<2>, f64>>
     + FunctionalContributionDual<Dual3<DualVec64<3>, f64>>
     + Display
+    + Sync
+    + Send
 {
     fn first_partial_derivatives(
         &self,
@@ -167,5 +169,7 @@ impl<T> FunctionalContribution for T where
         + FunctionalContributionDual<Dual3<DualVec64<2>, f64>>
         + FunctionalContributionDual<Dual3<DualVec64<3>, f64>>
         + Display
+        + Sync
+        + Send
 {
 }