feos-org · g-bauer · Oct 18, 2022 · Sep 27, 2022 · Sep 28, 2022 · Sep 30, 2022
diff --git a/Cargo.toml b/Cargo.toml
@@ -37,6 +37,7 @@ serde = "1.0"
 serde_json = "1.0"
 lazy_static = { version = "1.4", optional = true }
 indexmap = "1.8"
+rayon_ = { version = "1.5", optional = true, package = "rayon" }
 
 [dependencies.pyo3]
 version = "0.16"
@@ -55,5 +56,6 @@ pcsaft = ["association"]
 gc_pcsaft = ["association"]
 uvtheory = ["lazy_static"]
 pets = []
-python = ["pyo3", "numpy", "feos-core/python", "feos-dft?/python"]
+rayon = ["rayon_", "ndarray/rayon", "feos-core/rayon"]
+python = ["pyo3", "numpy", "feos-core/python", "feos-dft?/python", "rayon"]
 all_models = ["dft", "estimator", "pcsaft", "gc_pcsaft", "uvtheory", "pets"]
diff --git a/examples/pcsaft_phase_diagram.ipynb b/examples/pcsaft_phase_diagram.ipynb
@@ -490,7 +490,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.7"
+   "version": "3.9.12"
   }
  },
  "nbformat": 4,

diff --git a/feos-core/Cargo.toml b/feos-core/Cargo.toml
@@ -28,10 +28,12 @@ indexmap = "1.8"
 conv = "0.3"
 numpy = { version = "0.16", optional = true }
 pyo3 = { version = "0.16", optional = true }
+rayon_ = { version = "1.5", optional = true, package = "rayon" }
 
 [dev-dependencies]
 approx = "0.4"
 
 [features]
 default = []
-python = ["pyo3", "numpy", "quantity/python", "num-dual/python"]
+rayon = ["rayon_", "ndarray/rayon"]
+python = ["pyo3", "numpy", "quantity/python", "num-dual/python", "rayon"]
diff --git a/feos-core/src/cubic.rs b/feos-core/src/cubic.rs
@@ -18,7 +18,7 @@ use quantity::si::{SIArray1, SIUnit};
 use serde::{Deserialize, Serialize};
 use std::f64::consts::SQRT_2;
 use std::fmt;
-use std::rc::Rc;
+use std::sync::Arc;
 
 const KB_A3: f64 = 13806490.0;
 
@@ -167,7 +167,7 @@ impl Parameter for PengRobinsonParameters {
 }
 
 struct PengRobinsonContribution {
-    parameters: Rc<PengRobinsonParameters>,
+    parameters: Arc<PengRobinsonParameters>,
 }
 
 impl<D: DualNum<f64>> HelmholtzEnergyDual<D> for PengRobinsonContribution {
@@ -206,7 +206,7 @@ impl fmt::Display for PengRobinsonContribution {
 /// A simple version of the Peng-Robinson equation of state.
 pub struct PengRobinson {
     /// Parameters
-    parameters: Rc<PengRobinsonParameters>,
+    parameters: Arc<PengRobinsonParameters>,
     /// Ideal gas contributions to the Helmholtz energy
     ideal_gas: Joback,
     /// Non-ideal contributions to the Helmholtz energy
@@ -215,7 +215,7 @@ pub struct PengRobinson {
 
 impl PengRobinson {
     /// Create a new equation of state from a set of parameters.
-    pub fn new(parameters: Rc<PengRobinsonParameters>) -> Self {
+    pub fn new(parameters: Arc<PengRobinsonParameters>) -> Self {
         let ideal_gas = parameters.joback_records.as_ref().map_or_else(
             || Joback::default(parameters.tc.len()),
             |j| Joback::new(j.clone()),
@@ -238,7 +238,7 @@ impl EquationOfState for PengRobinson {
     }
 
     fn subset(&self, component_list: &[usize]) -> Self {
-        Self::new(Rc::new(self.parameters.subset(component_list)))
+        Self::new(Arc::new(self.parameters.subset(component_list)))
     }
 
     fn compute_max_density(&self, moles: &Array1<f64>) -> f64 {
@@ -270,7 +270,7 @@ mod tests {
     use crate::{EosResult, Verbosity};
     use approx::*;
     use quantity::si::*;
-    use std::rc::Rc;
+    use std::sync::Arc;
 
     fn pure_record_vec() -> Vec<PureRecord<PengRobinsonRecord, JobackRecord>> {
         let records = r#"[
@@ -317,7 +317,7 @@ mod tests {
         let tc = propane.model_record.tc;
         let pc = propane.model_record.pc;
         let parameters = PengRobinsonParameters::from_records(vec![propane], Array2::zeros((1, 1)));
-        let pr = Rc::new(PengRobinson::new(Rc::new(parameters)));
+        let pr = Arc::new(PengRobinson::new(Arc::new(parameters)));
         let options = SolverOptions::new().verbosity(Verbosity::Iter);
         let cp = State::critical_point(&pr, None, None, options)?;
         println!("{} {}", cp.temperature, cp.pressure(Contributions::Total));

diff --git a/feos-core/src/density_iteration.rs b/feos-core/src/density_iteration.rs
@@ -3,10 +3,10 @@ use crate::errors::{EosError, EosResult};
 use crate::state::State;
 use crate::EosUnit;
 use quantity::{QuantityArray1, QuantityScalar};
-use std::rc::Rc;
+use std::sync::Arc;
 
 pub fn density_iteration<U: EosUnit, E: EquationOfState>(
-    eos: &Rc<E>,
+    eos: &Arc<E>,
     temperature: QuantityScalar<U>,
     pressure: QuantityScalar<U>,
     moles: &QuantityArray1<U>,
@@ -142,7 +142,7 @@ pub fn density_iteration<U: EosUnit, E: EquationOfState>(
 }
 
 fn pressure_spinodal<U: EosUnit, E: EquationOfState>(
-    eos: &Rc<E>,
+    eos: &Arc<E>,
     temperature: QuantityScalar<U>,
     rho_init: QuantityScalar<U>,
     moles: &QuantityArray1<U>,

diff --git a/feos-core/src/equation_of_state.rs b/feos-core/src/equation_of_state.rs
@@ -36,6 +36,8 @@ pub trait HelmholtzEnergy:
     + HelmholtzEnergyDual<Dual3<DualVec64<2>, f64>>
     + HelmholtzEnergyDual<Dual3<DualVec64<3>, f64>>
     + fmt::Display
+    + Send
+    + Sync
 {
 }
 
@@ -52,6 +54,8 @@ impl<T> HelmholtzEnergy for T where
         + HelmholtzEnergyDual<Dual3<DualVec64<2>, f64>>
         + HelmholtzEnergyDual<Dual3<DualVec64<3>, f64>>
         + fmt::Display
+        + Send
+        + Sync
 {
 }
 
@@ -144,7 +148,7 @@ pub trait MolarWeight<U: EosUnit> {
 }
 
 /// A general equation of state.
-pub trait EquationOfState {
+pub trait EquationOfState: Send + Sync {
     /// Return the number of components of the equation of state.
     fn components(&self) -> usize;
 

diff --git a/feos-core/src/joback.rs b/feos-core/src/joback.rs
@@ -165,7 +165,7 @@ mod tests {
     use approx::assert_relative_eq;
     use ndarray::arr1;
     use quantity::si::*;
-    use std::rc::Rc;
+    use std::sync::Arc;
 
     use super::*;
 
@@ -259,7 +259,7 @@ mod tests {
         );
         assert_relative_eq!(jr.e, 0.0);
 
-        let eos = Rc::new(Joback::new(vec![jr]));
+        let eos = Arc::new(Joback::new(vec![jr]));
         let state = State::new_nvt(
             &eos,
             1000.0 * KELVIN,
@@ -281,7 +281,7 @@ mod tests {
     fn c_p_comparison() -> EosResult<()> {
         let record1 = JobackRecord::new(1.0, 0.2, 0.03, 0.004, 0.005);
         let record2 = JobackRecord::new(-5.0, 0.4, 0.03, 0.002, 0.001);
-        let joback = Rc::new(Joback::new(vec![record1, record2]));
+        let joback = Arc::new(Joback::new(vec![record1, record2]));
         let temperature = 300.0 * KELVIN;
         let volume = METER.powi(3);
         let moles = arr1(&[1.0, 3.0]) * MOL;

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 + Send + Sync {
     fn reference_temperature() -> QuantityScalar<Self>;
     fn reference_length() -> QuantityScalar<Self>;
     fn reference_density() -> QuantityScalar<Self>;

diff --git a/feos-core/src/phase_equilibria/bubble_dew.rs b/feos-core/src/phase_equilibria/bubble_dew.rs
@@ -10,7 +10,7 @@ use ndarray::*;
 use num_dual::linalg::{norm, LU};
 use quantity::{QuantityArray1, QuantityScalar};
 use std::convert::TryFrom;
-use std::rc::Rc;
+use std::sync::Arc;
 
 const MAX_ITER_INNER: usize = 5;
 const TOL_INNER: f64 = 1e-9;
@@ -61,7 +61,7 @@ impl<U: EosUnit, E: EquationOfState> PhaseEquilibrium<U, E, 2> {
     /// Calculate a phase equilibrium for a given temperature
     /// or pressure and composition of the liquid phase.
     pub fn bubble_point(
-        eos: &Rc<E>,
+        eos: &Arc<E>,
         temperature_or_pressure: QuantityScalar<U>,
         liquid_molefracs: &Array1<f64>,
         tp_init: Option<QuantityScalar<U>>,
@@ -85,7 +85,7 @@ impl<U: EosUnit, E: EquationOfState> PhaseEquilibrium<U, E, 2> {
     /// Calculate a phase equilibrium for a given temperature
     /// or pressure and composition of the vapor phase.
     pub fn dew_point(
-        eos: &Rc<E>,
+        eos: &Arc<E>,
         temperature_or_pressure: QuantityScalar<U>,
         vapor_molefracs: &Array1<f64>,
         tp_init: Option<QuantityScalar<U>>,
@@ -107,7 +107,7 @@ impl<U: EosUnit, E: EquationOfState> PhaseEquilibrium<U, E, 2> {
     }
 
     pub(super) fn bubble_dew_point(
-        eos: &Rc<E>,
+        eos: &Arc<E>,
         tp_spec: TPSpec<U>,
         tp_init: Option<QuantityScalar<U>>,
         molefracs_spec: &Array1<f64>,
@@ -183,7 +183,7 @@ impl<U: EosUnit, E: EquationOfState> PhaseEquilibrium<U, E, 2> {
     }
 
     fn iterate_bubble_dew(
-        eos: &Rc<E>,
+        eos: &Arc<E>,
         tp_spec: TPSpec<U>,
         tp_init: QuantityScalar<U>,
         molefracs_spec: &Array1<f64>,
@@ -204,7 +204,7 @@ impl<U: EosUnit, E: EquationOfState> PhaseEquilibrium<U, E, 2> {
     }
 
     fn starting_pressure_ideal_gas(
-        eos: &Rc<E>,
+        eos: &Arc<E>,
         temperature: QuantityScalar<U>,
         molefracs_spec: &Array1<f64>,
         bubble: bool,
@@ -220,7 +220,7 @@ impl<U: EosUnit, E: EquationOfState> PhaseEquilibrium<U, E, 2> {
     }
 
     pub(super) fn starting_pressure_ideal_gas_bubble(
-        eos: &Rc<E>,
+        eos: &Arc<E>,
         temperature: QuantityScalar<U>,
         liquid_molefracs: &Array1<f64>,
     ) -> EosResult<(QuantityScalar<U>, Array1<f64>)> {
@@ -239,7 +239,7 @@ impl<U: EosUnit, E: EquationOfState> PhaseEquilibrium<U, E, 2> {
     }
 
     fn starting_pressure_ideal_gas_dew(
-        eos: &Rc<E>,
+        eos: &Arc<E>,
         temperature: QuantityScalar<U>,
         vapor_molefracs: &Array1<f64>,
     ) -> EosResult<(QuantityScalar<U>, Array1<f64>)>
@@ -274,7 +274,7 @@ impl<U: EosUnit, E: EquationOfState> PhaseEquilibrium<U, E, 2> {
     }
 
     pub(super) fn starting_pressure_spinodal(
-        eos: &Rc<E>,
+        eos: &Arc<E>,
         temperature: QuantityScalar<U>,
         molefracs: &Array1<f64>,
     ) -> EosResult<QuantityScalar<U>>
@@ -290,7 +290,7 @@ impl<U: EosUnit, E: EquationOfState> PhaseEquilibrium<U, E, 2> {
 }
 
 fn starting_x2_bubble<U: EosUnit, E: EquationOfState>(
-    eos: &Rc<E>,
+    eos: &Arc<E>,
     temperature: QuantityScalar<U>,
     pressure: QuantityScalar<U>,
     liquid_molefracs: &Array1<f64>,
@@ -318,7 +318,7 @@ fn starting_x2_bubble<U: EosUnit, E: EquationOfState>(
 }
 
 fn starting_x2_dew<U: EosUnit, E: EquationOfState>(
-    eos: &Rc<E>,
+    eos: &Arc<E>,
     temperature: QuantityScalar<U>,
     pressure: QuantityScalar<U>,
     vapor_molefracs: &Array1<f64>,

diff --git a/feos-core/src/phase_equilibria/mod.rs b/feos-core/src/phase_equilibria/mod.rs
@@ -5,7 +5,7 @@ use crate::EosUnit;
 use quantity::{QuantityArray1, QuantityScalar};
 use std::fmt;
 use std::fmt::Write;
-use std::rc::Rc;
+use std::sync::Arc;
 
 mod bubble_dew;
 mod phase_diagram_binary;
@@ -196,7 +196,7 @@ impl<U: EosUnit, E: EquationOfState> PhaseEquilibrium<U, E, 2> {
     }
 
     pub(super) fn new_npt(
-        eos: &Rc<E>,
+        eos: &Arc<E>,
         temperature: QuantityScalar<U>,
         pressure: QuantityScalar<U>,
         vapor_moles: &QuantityArray1<U>,

diff --git a/feos-core/src/phase_equilibria/phase_diagram_binary.rs b/feos-core/src/phase_equilibria/phase_diagram_binary.rs
@@ -7,7 +7,7 @@ use ndarray::{arr1, arr2, concatenate, s, Array1, Array2, Axis};
 use num_dual::linalg::{norm, LU};
 use quantity::{QuantityArray1, QuantityScalar};
 use std::convert::{TryFrom, TryInto};
-use std::rc::Rc;
+use std::sync::Arc;
 
 const DEFAULT_POINTS: usize = 51;
 
@@ -19,7 +19,7 @@ impl<U: EosUnit, E: EquationOfState> PhaseDiagram<U, E> {
     /// phases are known, they can be passed as `x_lle` to avoid
     /// the calculation of unstable branches.
     pub fn binary_vle(
-        eos: &Rc<E>,
+        eos: &Arc<E>,
         temperature_or_pressure: QuantityScalar<U>,
         npoints: Option<usize>,
         x_lle: Option<(f64, f64)>,
@@ -99,7 +99,7 @@ impl<U: EosUnit, E: EquationOfState> PhaseDiagram<U, E> {
 
     #[allow(clippy::type_complexity)]
     fn calculate_vlle(
-        eos: &Rc<E>,
+        eos: &Arc<E>,
         tp: TPSpec<U>,
         npoints: usize,
         x_lle: (f64, f64),
@@ -147,7 +147,7 @@ impl<U: EosUnit, E: EquationOfState> PhaseDiagram<U, E> {
     /// liquid diagrams as well, as long as the feed composition is
     /// in a two phase region.
     pub fn lle(
-        eos: &Rc<E>,
+        eos: &Arc<E>,
         temperature_or_pressure: QuantityScalar<U>,
         feed: &QuantityArray1<U>,
         min_tp: QuantityScalar<U>,
@@ -184,7 +184,7 @@ impl<U: EosUnit, E: EquationOfState> PhaseDiagram<U, E> {
 }
 
 fn iterate_vle<U: EosUnit, E: EquationOfState>(
-    eos: &Rc<E>,
+    eos: &Arc<E>,
     tp: TPSpec<U>,
     x_lim: &[f64],
     vle_0: PhaseEquilibrium<U, E, 2>,
@@ -266,7 +266,7 @@ impl<U: EosUnit, E: EquationOfState> PhaseDiagram<U, E> {
     /// The `x_lle` parameter is used as initial values for the calculation
     /// of the heteroazeotrope.
     pub fn binary_vlle(
-        eos: &Rc<E>,
+        eos: &Arc<E>,
         temperature_or_pressure: QuantityScalar<U>,
         x_lle: (f64, f64),
         tp_lim_lle: Option<QuantityScalar<U>>,
@@ -369,7 +369,7 @@ where
     /// Calculate a heteroazeotrope (three phase equilbrium) for a binary
     /// system and given pressure.
     pub fn heteroazeotrope(
-        eos: &Rc<E>,
+        eos: &Arc<E>,
         temperature_or_pressure: QuantityScalar<U>,
         x_init: (f64, f64),
         tp_init: Option<QuantityScalar<U>>,
@@ -389,7 +389,7 @@ where
     /// Calculate a heteroazeotrope (three phase equilbrium) for a binary
     /// system and given temperature.
     fn heteroazeotrope_t(
-        eos: &Rc<E>,
+        eos: &Arc<E>,
         temperature: QuantityScalar<U>,
         x_init: (f64, f64),
         p_init: Option<QuantityScalar<U>>,
@@ -541,7 +541,7 @@ where
     /// Calculate a heteroazeotrope (three phase equilbrium) for a binary
     /// system and given pressure.
     fn heteroazeotrope_p(
-        eos: &Rc<E>,
+        eos: &Arc<E>,
         pressure: QuantityScalar<U>,
         x_init: (f64, f64),
         t_init: Option<QuantityScalar<U>>,