diff --git a/Cargo.toml b/Cargo.toml index 1e93d4acd..4e2a63053 100644 --- 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 } [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 = ["dep: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 index 589c8db08..0268e309f 100644 --- 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 index c2e8d9d50..86e912327 100644 --- 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 } [dev-dependencies] approx = "0.4" [features] default = [] -python = ["pyo3", "numpy", "quantity/python", "num-dual/python"] +rayon = ["dep: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 index e741cecdb..d78774b4d 100644 --- 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, + parameters: Arc, } impl> HelmholtzEnergyDual 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, + parameters: Arc, /// 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) -> Self { + pub fn new(parameters: Arc) -> 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 { @@ -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> { 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 index 67c8cd4b7..c1ce504f9 100644 --- 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( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, pressure: QuantityScalar, moles: &QuantityArray1, @@ -142,7 +142,7 @@ pub fn density_iteration( } fn pressure_spinodal( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, rho_init: QuantityScalar, moles: &QuantityArray1, diff --git a/feos-core/src/equation_of_state.rs b/feos-core/src/equation_of_state.rs index 621874759..e2efca5c7 100644 --- a/feos-core/src/equation_of_state.rs +++ b/feos-core/src/equation_of_state.rs @@ -36,6 +36,8 @@ pub trait HelmholtzEnergy: + HelmholtzEnergyDual, f64>> + HelmholtzEnergyDual, f64>> + fmt::Display + + Send + + Sync { } @@ -52,6 +54,8 @@ impl HelmholtzEnergy for T where + HelmholtzEnergyDual, f64>> + HelmholtzEnergyDual, f64>> + fmt::Display + + Send + + Sync { } @@ -144,7 +148,7 @@ pub trait MolarWeight { } /// 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/errors.rs b/feos-core/src/errors.rs index 9561d5c06..0da42ea27 100644 --- a/feos-core/src/errors.rs +++ b/feos-core/src/errors.rs @@ -30,6 +30,9 @@ pub enum EosError { ParameterError(#[from] ParameterError), #[error(transparent)] LinAlgError(#[from] LinAlgError), + #[cfg(feature = "rayon")] + #[error(transparent)] + RayonError(#[from] rayon::ThreadPoolBuildError), } /// Convenience type for `Result`. diff --git a/feos-core/src/joback.rs b/feos-core/src/joback.rs index 1ec090249..fa1c2bd23 100644 --- 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 index 45e16b151..c8639edf7 100644 --- 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; fn reference_length() -> QuantityScalar; fn reference_density() -> QuantityScalar; diff --git a/feos-core/src/phase_equilibria/bubble_dew.rs b/feos-core/src/phase_equilibria/bubble_dew.rs index d61cb0810..c2217793d 100644 --- 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 PhaseEquilibrium { /// Calculate a phase equilibrium for a given temperature /// or pressure and composition of the liquid phase. pub fn bubble_point( - eos: &Rc, + eos: &Arc, temperature_or_pressure: QuantityScalar, liquid_molefracs: &Array1, tp_init: Option>, @@ -85,7 +85,7 @@ impl PhaseEquilibrium { /// Calculate a phase equilibrium for a given temperature /// or pressure and composition of the vapor phase. pub fn dew_point( - eos: &Rc, + eos: &Arc, temperature_or_pressure: QuantityScalar, vapor_molefracs: &Array1, tp_init: Option>, @@ -107,7 +107,7 @@ impl PhaseEquilibrium { } pub(super) fn bubble_dew_point( - eos: &Rc, + eos: &Arc, tp_spec: TPSpec, tp_init: Option>, molefracs_spec: &Array1, @@ -183,7 +183,7 @@ impl PhaseEquilibrium { } fn iterate_bubble_dew( - eos: &Rc, + eos: &Arc, tp_spec: TPSpec, tp_init: QuantityScalar, molefracs_spec: &Array1, @@ -204,7 +204,7 @@ impl PhaseEquilibrium { } fn starting_pressure_ideal_gas( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, molefracs_spec: &Array1, bubble: bool, @@ -220,7 +220,7 @@ impl PhaseEquilibrium { } pub(super) fn starting_pressure_ideal_gas_bubble( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, liquid_molefracs: &Array1, ) -> EosResult<(QuantityScalar, Array1)> { @@ -239,7 +239,7 @@ impl PhaseEquilibrium { } fn starting_pressure_ideal_gas_dew( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, vapor_molefracs: &Array1, ) -> EosResult<(QuantityScalar, Array1)> @@ -274,7 +274,7 @@ impl PhaseEquilibrium { } pub(super) fn starting_pressure_spinodal( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, molefracs: &Array1, ) -> EosResult> @@ -290,7 +290,7 @@ impl PhaseEquilibrium { } fn starting_x2_bubble( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, pressure: QuantityScalar, liquid_molefracs: &Array1, @@ -318,7 +318,7 @@ fn starting_x2_bubble( } fn starting_x2_dew( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, pressure: QuantityScalar, vapor_molefracs: &Array1, diff --git a/feos-core/src/phase_equilibria/mod.rs b/feos-core/src/phase_equilibria/mod.rs index b5cf20255..91c3d75d4 100644 --- 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 PhaseEquilibrium { } pub(super) fn new_npt( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, pressure: QuantityScalar, vapor_moles: &QuantityArray1, diff --git a/feos-core/src/phase_equilibria/phase_diagram_binary.rs b/feos-core/src/phase_equilibria/phase_diagram_binary.rs index eb865192e..b0a171fad 100644 --- 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 PhaseDiagram { /// phases are known, they can be passed as `x_lle` to avoid /// the calculation of unstable branches. pub fn binary_vle( - eos: &Rc, + eos: &Arc, temperature_or_pressure: QuantityScalar, npoints: Option, x_lle: Option<(f64, f64)>, @@ -99,7 +99,7 @@ impl PhaseDiagram { #[allow(clippy::type_complexity)] fn calculate_vlle( - eos: &Rc, + eos: &Arc, tp: TPSpec, npoints: usize, x_lle: (f64, f64), @@ -147,7 +147,7 @@ impl PhaseDiagram { /// liquid diagrams as well, as long as the feed composition is /// in a two phase region. pub fn lle( - eos: &Rc, + eos: &Arc, temperature_or_pressure: QuantityScalar, feed: &QuantityArray1, min_tp: QuantityScalar, @@ -184,7 +184,7 @@ impl PhaseDiagram { } fn iterate_vle( - eos: &Rc, + eos: &Arc, tp: TPSpec, x_lim: &[f64], vle_0: PhaseEquilibrium, @@ -266,7 +266,7 @@ impl PhaseDiagram { /// The `x_lle` parameter is used as initial values for the calculation /// of the heteroazeotrope. pub fn binary_vlle( - eos: &Rc, + eos: &Arc, temperature_or_pressure: QuantityScalar, x_lle: (f64, f64), tp_lim_lle: Option>, @@ -369,7 +369,7 @@ where /// Calculate a heteroazeotrope (three phase equilbrium) for a binary /// system and given pressure. pub fn heteroazeotrope( - eos: &Rc, + eos: &Arc, temperature_or_pressure: QuantityScalar, x_init: (f64, f64), tp_init: Option>, @@ -389,7 +389,7 @@ where /// Calculate a heteroazeotrope (three phase equilbrium) for a binary /// system and given temperature. fn heteroazeotrope_t( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, x_init: (f64, f64), p_init: Option>, @@ -541,7 +541,7 @@ where /// Calculate a heteroazeotrope (three phase equilbrium) for a binary /// system and given pressure. fn heteroazeotrope_p( - eos: &Rc, + eos: &Arc, pressure: QuantityScalar, x_init: (f64, f64), t_init: Option>, diff --git a/feos-core/src/phase_equilibria/phase_diagram_pure.rs b/feos-core/src/phase_equilibria/phase_diagram_pure.rs index 31ccdfb17..4dba91e47 100644 --- a/feos-core/src/phase_equilibria/phase_diagram_pure.rs +++ b/feos-core/src/phase_equilibria/phase_diagram_pure.rs @@ -3,8 +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}; -use std::rc::Rc; +#[cfg(feature = "rayon")] +use rayon::{prelude::*, ThreadPool}; +use std::sync::Arc; /// Pure component and binary mixture phase diagrams. pub struct PhaseDiagram { @@ -22,7 +26,7 @@ impl Clone for PhaseDiagram { impl PhaseDiagram { /// Calculate a phase diagram for a pure component. pub fn pure( - eos: &Rc, + eos: &Arc, min_temperature: QuantityScalar, npoints: usize, critical_temperature: Option>, @@ -61,3 +65,66 @@ impl PhaseDiagram { self.states.iter().map(|s| s.liquid()).collect() } } + +#[cfg(feature = "rayon")] +impl PhaseDiagram { + fn solve_temperatures( + eos: &Arc, + temperatures: ArrayView1, + options: SolverOptions, + ) -> EosResult>> + where + QuantityScalar: std::fmt::Display + std::fmt::LowerExp, + { + let mut states = Vec::with_capacity(temperatures.len()); + let mut vle = None; + for ti in temperatures { + 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 par_pure( + eos: &Arc, + min_temperature: QuantityScalar, + npoints: usize, + chunksize: usize, + thread_pool: ThreadPool, + critical_temperature: Option>, + options: SolverOptions, + ) -> EosResult + where + QuantityScalar: 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> = 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)); + Ok(PhaseDiagram { states }) + } +} diff --git a/feos-core/src/phase_equilibria/phase_envelope.rs b/feos-core/src/phase_equilibria/phase_envelope.rs index 161fb09a2..a51c96eef 100644 --- a/feos-core/src/phase_equilibria/phase_envelope.rs +++ b/feos-core/src/phase_equilibria/phase_envelope.rs @@ -4,12 +4,12 @@ use crate::errors::EosResult; use crate::state::State; use crate::{Contributions, EosUnit}; use quantity::{QuantityArray1, QuantityScalar}; -use std::rc::Rc; +use std::sync::Arc; impl PhaseDiagram { /// Calculate the bubble point line of a mixture with given composition. pub fn bubble_point_line( - eos: &Rc, + eos: &Arc, moles: &QuantityArray1, min_temperature: QuantityScalar, npoints: usize, @@ -61,7 +61,7 @@ impl PhaseDiagram { /// Calculate the dew point line of a mixture with given composition. pub fn dew_point_line( - eos: &Rc, + eos: &Arc, moles: &QuantityArray1, min_temperature: QuantityScalar, npoints: usize, @@ -129,7 +129,7 @@ impl PhaseDiagram { /// Calculate the spinodal lines for a mixture with fixed composition. pub fn spinodal( - eos: &Rc, + eos: &Arc, moles: &QuantityArray1, min_temperature: QuantityScalar, npoints: usize, diff --git a/feos-core/src/phase_equilibria/tp_flash.rs b/feos-core/src/phase_equilibria/tp_flash.rs index 05546e255..305416828 100644 --- a/feos-core/src/phase_equilibria/tp_flash.rs +++ b/feos-core/src/phase_equilibria/tp_flash.rs @@ -6,7 +6,7 @@ use crate::EosUnit; use ndarray::*; use num_dual::linalg::norm; use quantity::{QuantityArray1, QuantityScalar}; -use std::rc::Rc; +use std::sync::Arc; const MAX_ITER_TP: usize = 400; const TOL_TP: f64 = 1e-8; @@ -19,7 +19,7 @@ impl PhaseEquilibrium { /// The algorithm can be use to calculate phase equilibria of systems /// containing non-volatile components (e.g. ions). pub fn tp_flash( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, pressure: QuantityScalar, feed: &QuantityArray1, diff --git a/feos-core/src/phase_equilibria/vle_pure.rs b/feos-core/src/phase_equilibria/vle_pure.rs index 8b3130d3c..d5cdafd7d 100644 --- a/feos-core/src/phase_equilibria/vle_pure.rs +++ b/feos-core/src/phase_equilibria/vle_pure.rs @@ -6,7 +6,7 @@ use crate::EosUnit; use ndarray::{arr1, Array1}; use quantity::QuantityScalar; use std::convert::TryFrom; -use std::rc::Rc; +use std::sync::Arc; const SCALE_T_NEW: f64 = 0.7; @@ -17,7 +17,7 @@ const TOL_PURE: f64 = 1e-12; impl PhaseEquilibrium { /// Calculate a phase equilibrium for a pure component. pub fn pure( - eos: &Rc, + eos: &Arc, temperature_or_pressure: QuantityScalar, initial_state: Option<&PhaseEquilibrium>, options: SolverOptions, @@ -34,7 +34,7 @@ impl PhaseEquilibrium { /// Calculate a phase equilibrium for a pure component /// and given temperature. fn pure_t( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, initial_state: Option<&PhaseEquilibrium>, options: SolverOptions, @@ -168,7 +168,7 @@ impl PhaseEquilibrium { /// Calculate a phase equilibrium for a pure component /// and given pressure. fn pure_p( - eos: &Rc, + eos: &Arc, pressure: QuantityScalar, initial_state: Option<&Self>, options: SolverOptions, @@ -271,13 +271,13 @@ impl PhaseEquilibrium { Ok(Self([vapor, liquid])) } - fn init_pure_ideal_gas(eos: &Rc, temperature: QuantityScalar) -> EosResult { + fn init_pure_ideal_gas(eos: &Arc, temperature: QuantityScalar) -> EosResult { let m = arr1(&[1.0]) * U::reference_moles(); let p = Self::starting_pressure_ideal_gas_bubble(eos, temperature, &arr1(&[1.0]))?.0; PhaseEquilibrium::new_npt(eos, temperature, p, &m, &m)?.check_trivial_solution() } - fn init_pure_spinodal(eos: &Rc, temperature: QuantityScalar) -> EosResult + fn init_pure_spinodal(eos: &Arc, temperature: QuantityScalar) -> EosResult where QuantityScalar: std::fmt::Display, { @@ -287,7 +287,7 @@ impl PhaseEquilibrium { } /// Initialize a new VLE for a pure substance for a given pressure. - fn init_pure_p(eos: &Rc, pressure: QuantityScalar) -> EosResult + fn init_pure_p(eos: &Arc, pressure: QuantityScalar) -> EosResult where QuantityScalar: std::fmt::Display, { @@ -363,7 +363,7 @@ impl PhaseEquilibrium { /// Calculate the pure component vapor pressures of all /// components in the system for the given temperature. pub fn vapor_pressure( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, ) -> Vec>> where @@ -371,7 +371,7 @@ impl PhaseEquilibrium { { (0..eos.components()) .map(|i| { - let pure_eos = Rc::new(eos.subset(&[i])); + let pure_eos = Arc::new(eos.subset(&[i])); PhaseEquilibrium::pure_t(&pure_eos, temperature, None, SolverOptions::default()) .map(|vle| vle.vapor().pressure(Contributions::Total)) .ok() @@ -382,7 +382,7 @@ impl PhaseEquilibrium { /// Calculate the pure component boiling temperatures of all /// components in the system for the given pressure. pub fn boiling_temperature( - eos: &Rc, + eos: &Arc, pressure: QuantityScalar, ) -> Vec>> where @@ -390,7 +390,7 @@ impl PhaseEquilibrium { { (0..eos.components()) .map(|i| { - let pure_eos = Rc::new(eos.subset(&[i])); + let pure_eos = Arc::new(eos.subset(&[i])); PhaseEquilibrium::pure_p(&pure_eos, pressure, None, SolverOptions::default()) .map(|vle| vle.vapor().temperature) .ok() @@ -401,7 +401,7 @@ impl PhaseEquilibrium { /// Calculate the pure component phase equilibria of all /// components in the system. pub fn vle_pure_comps( - eos: &Rc, + eos: &Arc, temperature_or_pressure: QuantityScalar, ) -> Vec>> where @@ -409,7 +409,7 @@ impl PhaseEquilibrium { { (0..eos.components()) .map(|i| { - let pure_eos = Rc::new(eos.subset(&[i])); + let pure_eos = Arc::new(eos.subset(&[i])); PhaseEquilibrium::pure( &pure_eos, temperature_or_pressure, diff --git a/feos-core/src/python/cubic.rs b/feos-core/src/python/cubic.rs index a6f51d4ec..8096842d3 100644 --- a/feos-core/src/python/cubic.rs +++ b/feos-core/src/python/cubic.rs @@ -11,10 +11,10 @@ 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)] +#[pyclass(name = "PengRobinsonRecord")] #[derive(Clone)] pub struct PyPengRobinsonRecord(PengRobinsonRecord); @@ -59,9 +59,9 @@ impl_binary_record!(); /// Returns /// ------- /// PengRobinsonParameters -#[pyclass(name = "PengRobinsonParameters", unsendable)] +#[pyclass(name = "PengRobinsonParameters")] #[derive(Clone)] -pub struct PyPengRobinsonParameters(pub Rc); +pub struct PyPengRobinsonParameters(pub Arc); impl_parameter!(PengRobinsonParameters, PyPengRobinsonParameters); diff --git a/feos-core/src/python/parameter.rs b/feos-core/src/python/parameter.rs index e486db91a..d6562c42b 100644 --- a/feos-core/src/python/parameter.rs +++ b/feos-core/src/python/parameter.rs @@ -570,10 +570,13 @@ 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 +592,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 +624,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 +647,7 @@ macro_rules! impl_parameter { binary_path: Option, search_option: Option, ) -> Result { - Ok(Self(Rc::new(<$parameter>::from_json( + Ok(Self(Arc::new(<$parameter>::from_json( substances, pure_path, binary_path, @@ -670,7 +673,7 @@ macro_rules! impl_parameter { binary_path: Option<&str>, search_option: Option, ) -> Result { - 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 +716,7 @@ macro_rules! impl_parameter_from_segments { segment_records: Vec, binary_segment_records: Option>, ) -> Result { - 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 +748,7 @@ macro_rules! impl_parameter_from_segments { binary_path: Option, search_option: Option, ) -> Result { - 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 index d1f8a3293..0f42288fb 100644 --- a/feos-core/src/python/phase_equilibria.rs +++ b/feos-core/src/python/phase_equilibria.rs @@ -1,949 +1,1007 @@ -#[macro_export] -macro_rules! impl_phase_equilibrium { - ($eos:ty, $py_eos:ty) => { - /// A thermodynamic two phase equilibrium state. - #[pyclass(name = "PhaseEquilibrium", unsendable)] - #[derive(Clone)] - pub struct PyPhaseEquilibrium(PhaseEquilibrium); - - #[pymethods] - impl PyPhaseEquilibrium { - /// Create a liquid and vapor state in equilibrium - /// for a pure substance. - /// - /// Parameters - /// ---------- - /// eos : EquationOfState - /// The equation of state. - /// temperature_or_pressure : SINumber - /// The system temperature or pressure. - /// initial_state : PhaseEquilibrium, optional - /// A phase equilibrium used as initial guess. - /// Can speed up convergence. - /// max_iter : int, optional - /// The maximum number of iterations. - /// tol: float, optional - /// The solution tolerance. - /// verbosity : Verbosity, optional - /// The verbosity. - /// - /// Returns - /// ------- - /// PhaseEquilibrium - /// - /// Raises - /// ------ - /// RuntimeError - /// When pressure iteration fails or no phase equilibrium is found. - #[staticmethod] - #[pyo3(text_signature = "(eos, temperature_or_pressure, initial_state=None, max_iter=None, tol=None, verbosity=None)")] - pub fn pure( - eos: $py_eos, - temperature_or_pressure: PySINumber, - initial_state: Option<&PyPhaseEquilibrium>, - max_iter: Option, - tol: Option, - verbosity: Option, - ) -> PyResult { - Ok(Self(PhaseEquilibrium::pure( - &eos.0, - temperature_or_pressure.into(), - initial_state.and_then(|s| Some(&s.0)), - (max_iter, tol, verbosity).into(), - )?)) - } - - /// Create a liquid and vapor state in equilibrium - /// for given temperature, pressure and feed composition. - /// - /// Can also be used to calculate liquid liquid phase separation. - /// - /// Parameters - /// ---------- - /// eos : EquationOfState - /// The equation of state. - /// temperature : SINumber - /// The system temperature. - /// pressure : SINumber - /// The system pressure. - /// feed : SIArray1 - /// Feed composition (units of amount of substance). - /// initial_state : PhaseEquilibrium, optional - /// A phase equilibrium used as initial guess. - /// Can speed up convergence. - /// max_iter : int, optional - /// The maximum number of iterations. - /// tol: float, optional - /// The solution tolerance. - /// verbosity : Verbosity, optional - /// The verbosity. - /// - /// Returns - /// ------- - /// PhaseEquilibrium - /// - /// Raises - /// ------ - /// RuntimeError - /// When pressure iteration fails or no phase equilibrium is found. - #[staticmethod] - #[pyo3(text_signature = "(eos, temperature, pressure, feed, initial_state=None, max_iter=None, tol=None, verbosity=None, non_volatile_components=None)")] - pub fn tp_flash( - eos: $py_eos, - temperature: PySINumber, - pressure: PySINumber, - feed: &PySIArray1, - initial_state: Option<&PyPhaseEquilibrium>, - max_iter: Option, - tol: Option, - verbosity: Option, - non_volatile_components: Option>, - ) -> PyResult { - Ok(Self(PhaseEquilibrium::tp_flash( - &eos.0, - temperature.into(), - pressure.into(), - feed, - initial_state.and_then(|s| Some(&s.0)), - (max_iter, tol, verbosity).into(), non_volatile_components - )?)) - } - - /// Compute a phase equilibrium for given temperature - /// or pressure and liquid mole fractions. - /// - /// Parameters - /// ---------- - /// eos : EquationOfState - /// The equation of state. - /// temperature_or_pressure : SINumber - /// The system temperature_or_pressure. - /// liquid_molefracs : numpy.ndarray - /// The mole fraction of the liquid phase. - /// tp_init : SINumber, optional - /// The system pressure/temperature used as starting - /// condition for the iteration. - /// vapor_molefracs : numpy.ndarray, optional - /// The mole fraction of the vapor phase used as - /// starting condition for iteration. - /// max_iter_inner : int, optional - /// The maximum number of inner iterations. - /// max_iter_outer : int, optional - /// The maximum number of outer iterations. - /// tol_inner : float, optional - /// The solution tolerance in the inner loop. - /// tol_outer : float, optional - /// The solution tolerance in the outer loop. - /// verbosity : Verbosity, optional - /// The verbosity. - /// - /// Returns - /// ------- - /// PhaseEquilibrium - #[staticmethod] - #[pyo3(text_signature = "(eos, temperature_or_pressure, liquid_molefracs, tp_init=None, vapor_molefracs=None, max_iter_inner=None, max_iter_outer=None, tol_inner=None, tol_outer=None, verbosity=None)")] - pub fn bubble_point( - eos: $py_eos, - temperature_or_pressure: PySINumber, - liquid_molefracs: &PyArray1, - tp_init: Option, - vapor_molefracs: Option<&PyArray1>, - max_iter_inner: Option, - max_iter_outer: Option, - tol_inner: Option, - tol_outer: Option, - verbosity: Option, - ) -> PyResult { - let x = vapor_molefracs.and_then(|m| Some(m.to_owned_array())); - Ok(Self(PhaseEquilibrium::bubble_point( - &eos.0, - temperature_or_pressure.into(), - &liquid_molefracs.to_owned_array(), - tp_init.map(|p| p.into()), - x.as_ref(), - ( - (max_iter_inner, tol_inner, verbosity).into(), - (max_iter_outer, tol_outer, verbosity).into() - ) - )?)) - } - - /// Compute a phase equilibrium for given temperature - /// or pressure and vapor mole fractions. - /// - /// Parameters - /// ---------- - /// eos : EquationOfState - /// The equation of state. - /// temperature_or_pressure : SINumber - /// The system temperature or pressure. - /// vapor_molefracs : numpy.ndarray - /// The mole fraction of the vapor phase. - /// tp_init : SINumber, optional - /// The system pressure/temperature used as starting - /// condition for the iteration. - /// liquid_molefracs : numpy.ndarray, optional - /// The mole fraction of the liquid phase used as - /// starting condition for iteration. - /// max_iter_inner : int, optional - /// The maximum number of inner iterations. - /// max_iter_outer : int, optional - /// The maximum number of outer iterations. - /// tol_inner : float, optional - /// The solution tolerance in the inner loop. - /// tol_outer : float, optional - /// The solution tolerance in the outer loop. - /// verbosity : Verbosity, optional - /// The verbosity. - /// - /// Returns - /// ------- - /// PhaseEquilibrium - #[staticmethod] - #[pyo3(text_signature = "(eos, temperature_or_pressure, vapor_molefracs, tp_init=None, liquid_molefracs=None, max_iter_inner=None, max_iter_outer=None, tol_inner=None, tol_outer=None, verbosity=None)")] - pub fn dew_point( - eos: $py_eos, - temperature_or_pressure: PySINumber, - vapor_molefracs: &PyArray1, - tp_init: Option, - liquid_molefracs: Option<&PyArray1>, - max_iter_inner: Option, - max_iter_outer: Option, - tol_inner: Option, - tol_outer: Option, - verbosity: Option, - ) -> PyResult { - let x = liquid_molefracs.and_then(|m| Some(m.to_owned_array())); - Ok(Self(PhaseEquilibrium::dew_point( - &eos.0, - temperature_or_pressure.into(), - &vapor_molefracs.to_owned_array(), - tp_init.map(|p| p.into()), - x.as_ref(), - ( - (max_iter_inner, tol_inner, verbosity).into(), - (max_iter_outer, tol_outer, verbosity).into() - ) - )?)) - } - - #[getter] - fn get_vapor(&self) -> PyState { - PyState(self.0.vapor().clone()) - } - - #[getter] - fn get_liquid(&self) -> PyState { - PyState(self.0.liquid().clone()) - } - - /// Calculate a new PhaseEquilibrium with the given chemical potential. - /// The temperature remains constant, but the states are not in - /// a mechanical equilibrium anymore. - /// - /// Parameters - /// ---------- - /// chemical_potential: SIArray1 - /// The new chemical potential - /// - #[pyo3(text_signature = "(chemical_potential)")] - fn update_chemical_potential(slf: &PyCell, chemical_potential: &PySIArray1) -> PyResult<()> { - slf.borrow_mut().0.update_chemical_potential(chemical_potential)?; - Ok(()) - } - - /// Calculate the pure component vapor-liquid equilibria for all - /// components in the system. - /// - /// Parameters - /// ---------- - /// eos : EquationOfState - /// The equation of state. - /// temperature_or_pressure : SINumber - /// The system temperature or pressure. - /// - /// Returns - /// ------- - /// list[PhaseEquilibrium] - #[staticmethod] - #[pyo3(text_signature = "(eos, temperature_or_pressure)")] - fn vle_pure_comps(eos: $py_eos, temperature_or_pressure: PySINumber) -> Vec> { - PhaseEquilibrium::vle_pure_comps(&eos.0, temperature_or_pressure.into()) - .into_iter() - .map(|o| o.map(Self)) - .collect() - } - - /// Calculate the pure component vapor pressures for all the - /// components in the system. - /// - /// Parameters - /// ---------- - /// eos : EquationOfState - /// The equation of state. - /// temperature : SINumber - /// The system temperature. - /// - /// Returns - /// ------- - /// list[SINumber] - #[staticmethod] - #[pyo3(text_signature = "(eos, temperature)")] - fn vapor_pressure(eos: $py_eos, temperature: PySINumber) -> Vec> { - PhaseEquilibrium::vapor_pressure(&eos.0, temperature.into()) - .into_iter() - .map(|o| o.map(|n| n.into())) - .collect() - } - - /// Calculate the pure component boiling temperatures for all the - /// components in the system. - /// - /// Parameters - /// ---------- - /// eos : EquationOfState - /// The equation of state. - /// pressure : SINumber - /// The system pressure. - /// - /// Returns - /// ------- - /// list[SINumber] - #[staticmethod] - #[pyo3(text_signature = "(eos, pressure)")] - fn boiling_temperature(eos: $py_eos, pressure: PySINumber) -> Vec> { - PhaseEquilibrium::boiling_temperature(&eos.0, pressure.into()) - .into_iter() - .map(|o| o.map(|n| n.into())) - .collect() - } - - fn _repr_markdown_(&self) -> String { - self.0._repr_markdown_() - } - - fn __repr__(&self) -> PyResult { - Ok(self.0.to_string()) - } - } - - /// A thermodynamic three phase equilibrium state. - #[pyclass(name = "ThreePhaseEquilibrium", unsendable)] - #[derive(Clone)] - struct PyThreePhaseEquilibrium(PhaseEquilibrium); - - #[pymethods] - impl PyPhaseEquilibrium { - /// Calculate a heteroazeotrope in a binary mixture for a given temperature - /// or pressure. - /// - /// Parameters - /// ---------- - /// eos : EquationOfState - /// The equation of state. - /// temperature_or_pressure : SINumber - /// The system temperature or pressure. - /// x_init : list[float] - /// Initial guesses for the liquid molefracs of component 1 - /// at the heteroazeotropic point. - /// tp_init : SINumber, optional - /// Initial guess for the temperature/pressure at the - /// heteroszeotropic point. - /// max_iter : int, optional - /// The maximum number of iterations. - /// tol: float, optional - /// The solution tolerance. - /// verbosity : Verbosity, optional - /// The verbosity. - /// max_iter_bd_inner : int, optional - /// The maximum number of inner iterations in the bubble/dew point iteration. - /// max_iter_bd_outer : int, optional - /// The maximum number of outer iterations in the bubble/dew point iteration. - /// tol_bd_inner : float, optional - /// The solution tolerance in the inner loop of the bubble/dew point iteration. - /// tol_bd_outer : float, optional - /// The solution tolerance in the outer loop of the bubble/dew point iteration. - /// verbosity_bd : Verbosity, optional - /// The verbosity of the bubble/dew point iteration. - #[staticmethod] - #[pyo3(text_signature = "(eos, temperature_or_pressure, x_init, tp_init=None, max_iter=None, tol=None, verbosity=None, max_iter_bd_inner=None, max_iter_bd_outer=None, tol_bd_inner=None, tol_bd_outer=None, verbosity_bd=None)")] - fn heteroazeotrope( - eos: $py_eos, - temperature_or_pressure: PySINumber, - x_init: (f64, f64), - tp_init: Option, - max_iter: Option, - tol: Option, - verbosity: Option, - max_iter_bd_inner: Option, - max_iter_bd_outer: Option, - tol_bd_inner: Option, - tol_bd_outer: Option, - verbosity_bd: Option, - ) -> PyResult { - Ok(PyThreePhaseEquilibrium(PhaseEquilibrium::heteroazeotrope( - &eos.0, - temperature_or_pressure.into(), - x_init, - tp_init.map(|t| t.into()), - (max_iter, tol, verbosity).into(), - ( - (max_iter_bd_inner, tol_bd_inner, verbosity_bd).into(), - (max_iter_bd_outer, tol_bd_outer, verbosity_bd).into(), - ) - )?)) - } - } - - #[pymethods] - impl PyThreePhaseEquilibrium { - #[getter] - fn get_vapor(&self) -> PyState { - PyState(self.0.vapor().clone()) - } - - #[getter] - fn get_liquid1(&self) -> PyState { - PyState(self.0.liquid1().clone()) - } - - #[getter] - fn get_liquid2(&self) -> PyState { - PyState(self.0.liquid2().clone()) - } - - fn _repr_markdown_(&self) -> String { - self.0._repr_markdown_() - } - - fn __repr__(&self) -> PyResult { - Ok(self.0.to_string()) - } - } - - #[pymethods] - impl PyState { - /// Calculates a two phase Tp-flash with the state as feed. - /// - /// Parameters - /// ---------- - /// initial_state : PhaseEquilibrium, optional - /// A phase equilibrium used as initial guess. - /// Can speed up convergence. - /// max_iter : int, optional - /// The maximum number of iterations. - /// tol: float, optional - /// The solution tolerance. - /// verbosity : Verbosity, optional - /// The verbosity. - /// - /// Returns - /// ------- - /// PhaseEquilibrium - /// - /// Raises - /// ------ - /// RuntimeError - /// When pressure iteration fails or no phase equilibrium is found. - #[pyo3(text_signature = "($self, initial_state=None, max_iter=None, tol=None, verbosity=None, non_volatile_components=None)")] - pub fn tp_flash( - &self, - initial_state: Option<&PyPhaseEquilibrium>, - max_iter: Option, - tol: Option, - verbosity: Option, - non_volatile_components: Option>, - ) -> PyResult { - Ok(PyPhaseEquilibrium(self.0.tp_flash( - initial_state.and_then(|s| Some(&s.0)), - (max_iter, tol, verbosity).into(), - non_volatile_components - )?)) - } - } - - /// Phase diagram for a pure component or a binary mixture. - #[pyclass(name = "PhaseDiagram", unsendable)] - pub struct PyPhaseDiagram(PhaseDiagram); - - #[pymethods] - impl PyPhaseDiagram { - /// Calculate a pure component phase diagram. - /// - /// Parameters - /// ---------- - /// eos: Eos - /// The equation of state. - /// min_temperature: SINumber - /// The lower limit for the temperature. - /// npoints: int - /// The number of points. - /// 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 - #[staticmethod] - #[pyo3(text_signature = "(eos, min_temperature, npoints, critical_temperature=None, max_iter=None, tol=None, verbosity=None)")] - pub fn pure( - eos: &$py_eos, - min_temperature: PySINumber, - npoints: usize, - critical_temperature: Option, - max_iter: Option, - tol: Option, - verbosity: Option, - ) -> PyResult { - let dia = PhaseDiagram::pure( - &eos.0, - min_temperature.into(), - npoints, - critical_temperature.map(|t| t.into()), - (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 - /// bubble point line while the vapor states contain the corresponding - /// equilibrium states at different compositions. - /// - /// Parameters - /// ---------- - /// eos: Eos - /// The equation of state. - /// moles: SIArray1 - /// The moles of the individual components - /// min_temperature: SINumber - /// The lower limit for the temperature. - /// npoints: int - /// The number of points. - /// 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_inner : int, optional - /// The maximum number of inner iterations in the bubble/dew point iteration. - /// max_iter_outer : int, optional - /// The maximum number of outer iterations in the bubble/dew point iteration. - /// tol_inner : float, optional - /// The solution tolerance in the inner loop of the bubble/dew point iteration. - /// tol_outer : float, optional - /// The solution tolerance in the outer loop of the bubble/dew point iteration. - /// verbosity : Verbosity, optional - /// The verbosity of the bubble/dew point iteration. - /// - /// Returns - /// ------- - /// PhaseDiagram - #[staticmethod] - #[pyo3(text_signature = "(eos, moles, min_temperature, npoints, critical_temperature=None, max_iter_inner=None, max_iter_outer=None, tol_inner=None, tol_outer=None, verbosity=None)")] - pub fn bubble_point_line( - eos: &$py_eos, - moles: PySIArray1, - min_temperature: PySINumber, - npoints: usize, - critical_temperature: Option, - max_iter_inner: Option, - max_iter_outer: Option, - tol_inner: Option, - tol_outer: Option, - verbosity: Option, - ) -> PyResult { - let dia = PhaseDiagram::bubble_point_line( - &eos.0, - &moles, - min_temperature.into(), - npoints, - critical_temperature.map(|t| t.into()), - ( - (max_iter_inner, tol_inner, verbosity).into(), - (max_iter_outer, tol_outer, verbosity).into(), - ) - )?; - Ok(Self(dia)) - } - - /// Calculate the dew point line of a mixture with given composition. - /// - /// In the resulting phase diagram, the vapor states correspond to the - /// dew point line while the liquid states contain the corresponding - /// equilibrium states at different compositions. - /// - /// Parameters - /// ---------- - /// eos: Eos - /// The equation of state. - /// moles: SIArray1 - /// The moles of the individual components - /// min_temperature: SINumber - /// The lower limit for the temperature. - /// npoints: int - /// The number of points. - /// 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_inner : int, optional - /// The maximum number of inner iterations in the bubble/dew point iteration. - /// max_iter_outer : int, optional - /// The maximum number of outer iterations in the bubble/dew point iteration. - /// tol_inner : float, optional - /// The solution tolerance in the inner loop of the bubble/dew point iteration. - /// tol_outer : float, optional - /// The solution tolerance in the outer loop of the bubble/dew point iteration. - /// verbosity : Verbosity, optional - /// The verbosity of the bubble/dew point iteration. - /// - /// Returns - /// ------- - /// PhaseDiagram - #[staticmethod] - #[pyo3(text_signature = "(eos, moles, min_temperature, npoints, critical_temperature=None, max_iter_inner=None, max_iter_outer=None, tol_inner=None, tol_outer=None, verbosity=None)")] - pub fn dew_point_line( - eos: &$py_eos, - moles: PySIArray1, - min_temperature: PySINumber, - npoints: usize, - critical_temperature: Option, - max_iter_inner: Option, - max_iter_outer: Option, - tol_inner: Option, - tol_outer: Option, - verbosity: Option, - ) -> PyResult { - let dia = PhaseDiagram::dew_point_line( - &eos.0, - &moles, - min_temperature.into(), - npoints, - critical_temperature.map(|t| t.into()), - ( - (max_iter_inner, tol_inner, verbosity).into(), - (max_iter_outer, tol_outer, verbosity).into(), - ) - )?; - Ok(Self(dia)) - } - - /// Calculate the spinodal lines for a mixture with fixed composition. - /// - /// Parameters - /// ---------- - /// eos: Eos - /// The equation of state. - /// moles: SIArray1 - /// The moles of the individual components - /// min_temperature: SINumber - /// The lower limit for the temperature. - /// npoints: int - /// The number of points. - /// 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 - #[staticmethod] - #[pyo3(text_signature = "(eos, moles, min_temperature, npoints, critical_temperature=None, max_iter=None, tol=None, verbosity=None)")] - pub fn spinodal( - eos: &$py_eos, - moles: PySIArray1, - min_temperature: PySINumber, - npoints: usize, - critical_temperature: Option, - max_iter: Option, - tol: Option, - verbosity: Option, - ) -> PyResult { - let dia = PhaseDiagram::spinodal( - &eos.0, - &moles, - min_temperature.into(), - npoints, - critical_temperature.map(|t| t.into()), - (max_iter, tol, verbosity).into(), - )?; - Ok(Self(dia)) - } - - #[getter] - pub fn get_states(&self) -> Vec { - self.0 - .states - .iter() - .map(|vle| PyPhaseEquilibrium(vle.clone())) - .collect() - } - - #[getter] - pub fn get_vapor(&self) -> PyStateVec { - self.0.vapor().into() - } - - #[getter] - pub fn get_liquid(&self) -> PyStateVec { - self.0.liquid().into() - } - - /// Returns the phase diagram as dictionary. - /// - /// Units - /// ----- - /// temperature : K - /// pressure : Pa - /// densities : mol / m³ - /// molar enthalpies : kJ / mol - /// molar entropies : kJ / mol / K - /// - /// Returns - /// ------- - /// dict[str, list[float]] - /// Keys: property names. Values: property for each state. - /// - /// Notes - /// ----- - /// xi: liquid molefraction of component i - /// yi: vapor molefraction of component i - /// i: component index according to order in parameters. - pub fn to_dict(&self) -> PyResult>> { - let n = self.0.states[0].liquid().eos.components(); - let mut dict = HashMap::with_capacity(8 + 2 * n); - if n != 1 { - let xs = self.0.liquid().molefracs(); - let ys = self.0.vapor().molefracs(); - for i in 0..n { - dict.insert(String::from(format!("x{}", i)), xs.column(i).to_vec()); - dict.insert(String::from(format!("y{}", i)), ys.column(i).to_vec()); - } - } - dict.insert(String::from("temperature"), (self.0.vapor().temperature() / KELVIN).into_value()?.into_raw_vec()); - dict.insert(String::from("pressure"), (self.0.vapor().pressure() / PASCAL).into_value()?.into_raw_vec()); - dict.insert(String::from("density liquid"), (self.0.liquid().density() / (MOL / METER.powi(3))).into_value()?.into_raw_vec()); - dict.insert(String::from("density vapor"), (self.0.vapor().density() / (MOL / METER.powi(3))).into_value()?.into_raw_vec()); - dict.insert(String::from("molar enthalpy liquid"), (self.0.liquid().molar_enthalpy() / (KILO*JOULE / MOL)).into_value()?.into_raw_vec()); - dict.insert(String::from("molar enthalpy vapor"), (self.0.vapor().molar_enthalpy() / (KILO*JOULE / MOL)).into_value()?.into_raw_vec()); - dict.insert(String::from("molar entropy liquid"), (self.0.liquid().molar_entropy() / (KILO*JOULE / KELVIN / MOL)).into_value()?.into_raw_vec()); - dict.insert(String::from("molar entropy vapor"), (self.0.vapor().molar_entropy() / (KILO*JOULE / KELVIN / MOL)).into_value()?.into_raw_vec()); - Ok(dict) - } - - /// Binary phase diagram calculated using bubble/dew point iterations. - /// - /// Parameters - /// ---------- - /// eos : EquationOfState - /// The equation of state. - /// temperature_or_pressure: SINumber - /// The constant temperature or pressure. - /// npoints: int, optional - /// The number of points (default 51). - /// x_lle: (float, float), optional - /// An estimate for the molefractions of component 1 - /// at the heteroazeotrop - /// max_iter_inner : int, optional - /// The maximum number of inner iterations in the bubble/dew point iteration. - /// max_iter_outer : int, optional - /// The maximum number of outer iterations in the bubble/dew point iteration. - /// tol_inner : float, optional - /// The solution tolerance in the inner loop of the bubble/dew point iteration. - /// tol_outer : float, optional - /// The solution tolerance in the outer loop of the bubble/dew point iteration. - /// verbosity : Verbosity, optional - /// The verbosity of the bubble/dew point iteration. - /// - /// Returns - /// ------- - /// PhaseDiagram - #[staticmethod] - #[pyo3(text_signature = "(eos, temperature_or_pressure, npoints=None, x_lle=None, max_iter_inner=None, max_iter_outer=None, tol_inner=None, tol_outer=None, verbosity=None)")] - pub fn binary_vle( - eos: $py_eos, - temperature_or_pressure: PySINumber, - npoints: Option, - x_lle: Option<(f64, f64)>, - max_iter_inner: Option, - max_iter_outer: Option, - tol_inner: Option, - tol_outer: Option, - verbosity: Option, - ) -> PyResult { - let dia = PhaseDiagram::binary_vle( - &eos.0, - temperature_or_pressure.into(), - npoints, - x_lle, - ( - (max_iter_inner, tol_inner, verbosity).into(), - (max_iter_outer, tol_outer, verbosity).into(), - ) - )?; - Ok(Self(dia)) - } - - /// Create a new phase diagram using Tp flash calculations. - /// - /// The usual use case for this function is the calculation of - /// liquid-liquid phase diagrams, but it can be used for vapor- - /// liquid diagrams as well, as long as the feed composition is - /// in a two phase region. - /// - /// Parameters - /// ---------- - /// eos : EquationOfState - /// The equation of state. - /// temperature_or_pressure: SINumber - /// The consant temperature or pressure. - /// feed: SIArray1 - /// Mole numbers in the (unstable) feed state. - /// min_tp: - /// The lower limit of the temperature/pressure range. - /// max_tp: - /// The upper limit of the temperature/pressure range. - /// npoints: int, optional - /// The number of points (default 51). - /// - /// Returns - /// ------- - /// PhaseDiagram - #[staticmethod] - #[pyo3(text_signature = "(eos, temperature_or_pressure, feed, min_tp, max_tp, npoints=None)")] - pub fn lle( - eos: $py_eos, - temperature_or_pressure: PySINumber, - feed: PySIArray1, - min_tp: PySINumber, - max_tp: PySINumber, - npoints: Option, - ) -> PyResult { - let dia = PhaseDiagram::lle( - &eos.0, - temperature_or_pressure.into(), - &feed, - min_tp.into(), - max_tp.into(), - npoints, - )?; - Ok(Self(dia)) - } - } - - /// Phase diagram for a binary mixture exhibiting a heteroazeotrope. - #[pyclass(name = "PhaseDiagramHetero", unsendable)] - pub struct PyPhaseDiagramHetero(PhaseDiagramHetero); - - #[pymethods] - impl PyPhaseDiagram { - /// Phase diagram for a binary mixture exhibiting a heteroazeotrope. - /// - /// Parameters - /// ---------- - /// eos: SaftFunctional - /// The SAFT Helmholtz energy functional. - /// temperature_or_pressure: SINumber - /// The temperature_or_pressure. - /// x_lle: SINumber - /// Initial values for the molefractions of component 1 - /// at the heteroazeotrop. - /// tp_lim_lle: SINumber, optional - /// The minimum temperature up to which the LLE is calculated. - /// If it is not provided, no LLE is calcualted. - /// tp_init_vlle: SINumber, optional - /// Initial value for the calculation of the VLLE. - /// npoints_vle: int, optional - /// The number of points for the VLE (default 51). - /// npoints_lle: int, optional - /// The number of points for the LLE (default 51). - /// max_iter_inner : int, optional - /// The maximum number of inner iterations in the bubble/dew point iteration. - /// max_iter_outer : int, optional - /// The maximum number of outer iterations in the bubble/dew point iteration. - /// tol_inner : float, optional - /// The solution tolerance in the inner loop of the bubble/dew point iteration. - /// tol_outer : float, optional - /// The solution tolerance in the outer loop of the bubble/dew point iteration. - /// verbosity : Verbosity, optional - /// The verbosity of the bubble/dew point iteration. - /// - /// Returns - /// ------- - /// PhaseDiagramHetero - #[staticmethod] - #[pyo3(text_signature = "(eos, temperature_or_pressure, x_lle, tp_lim_lle=None, tp_init_vlle=None, npoints_vle=None, npoints_lle=None, max_iter_bd_inner=None, max_iter_bd_outer=None, tol_bd_inner=None, tol_bd_outer=None, verbosity_bd=None)")] - pub fn binary_vlle( - eos: $py_eos, - temperature_or_pressure: PySINumber, - x_lle: (f64, f64), - tp_lim_lle: Option, - tp_init_vlle: Option, - npoints_vle: Option, - npoints_lle: Option, - max_iter_inner: Option, - max_iter_outer: Option, - tol_inner: Option, - tol_outer: Option, - verbosity: Option, - ) -> PyResult { - let dia = PhaseDiagram::binary_vlle( - &eos.0, - temperature_or_pressure.into(), - x_lle, - tp_lim_lle.map(|t| t.into()), - tp_init_vlle.map(|t| t.into()), - npoints_vle, - npoints_lle, - ( - (max_iter_inner, tol_inner, verbosity).into(), - (max_iter_outer, tol_outer, verbosity).into(), - ) - )?; - Ok(PyPhaseDiagramHetero(dia)) - } - } - - #[pymethods] - impl PyPhaseDiagramHetero { - #[getter] - pub fn get_vle(&self) -> PyPhaseDiagram { - PyPhaseDiagram(self.0.vle().clone()) - } - - #[getter] - pub fn get_vle1(&self) -> PyPhaseDiagram { - PyPhaseDiagram(self.0.vle1.clone()) - } - - #[getter] - pub fn get_vle2(&self) -> PyPhaseDiagram { - PyPhaseDiagram(self.0.vle2.clone()) - } - - #[getter] - pub fn get_lle(&self) -> Option { - self.0 - .lle - .as_ref() - .map(|d| PyPhaseDiagram(d.clone())) - } - } - } -} +#[macro_export] +macro_rules! impl_phase_equilibrium { + ($eos:ty, $py_eos:ty) => { + /// A thermodynamic two phase equilibrium state. + #[pyclass(name = "PhaseEquilibrium")] + #[derive(Clone)] + pub struct PyPhaseEquilibrium(PhaseEquilibrium); + + #[pymethods] + impl PyPhaseEquilibrium { + /// Create a liquid and vapor state in equilibrium + /// for a pure substance. + /// + /// Parameters + /// ---------- + /// eos : EquationOfState + /// The equation of state. + /// temperature_or_pressure : SINumber + /// The system temperature or pressure. + /// initial_state : PhaseEquilibrium, optional + /// A phase equilibrium used as initial guess. + /// Can speed up convergence. + /// max_iter : int, optional + /// The maximum number of iterations. + /// tol: float, optional + /// The solution tolerance. + /// verbosity : Verbosity, optional + /// The verbosity. + /// + /// Returns + /// ------- + /// PhaseEquilibrium + /// + /// Raises + /// ------ + /// RuntimeError + /// When pressure iteration fails or no phase equilibrium is found. + #[staticmethod] + #[pyo3(text_signature = "(eos, temperature_or_pressure, initial_state=None, max_iter=None, tol=None, verbosity=None)")] + pub fn pure( + eos: $py_eos, + temperature_or_pressure: PySINumber, + initial_state: Option<&PyPhaseEquilibrium>, + max_iter: Option, + tol: Option, + verbosity: Option, + ) -> PyResult { + Ok(Self(PhaseEquilibrium::pure( + &eos.0, + temperature_or_pressure.into(), + initial_state.and_then(|s| Some(&s.0)), + (max_iter, tol, verbosity).into(), + )?)) + } + + /// Create a liquid and vapor state in equilibrium + /// for given temperature, pressure and feed composition. + /// + /// Can also be used to calculate liquid liquid phase separation. + /// + /// Parameters + /// ---------- + /// eos : EquationOfState + /// The equation of state. + /// temperature : SINumber + /// The system temperature. + /// pressure : SINumber + /// The system pressure. + /// feed : SIArray1 + /// Feed composition (units of amount of substance). + /// initial_state : PhaseEquilibrium, optional + /// A phase equilibrium used as initial guess. + /// Can speed up convergence. + /// max_iter : int, optional + /// The maximum number of iterations. + /// tol: float, optional + /// The solution tolerance. + /// verbosity : Verbosity, optional + /// The verbosity. + /// + /// Returns + /// ------- + /// PhaseEquilibrium + /// + /// Raises + /// ------ + /// RuntimeError + /// When pressure iteration fails or no phase equilibrium is found. + #[staticmethod] + #[pyo3(text_signature = "(eos, temperature, pressure, feed, initial_state=None, max_iter=None, tol=None, verbosity=None, non_volatile_components=None)")] + pub fn tp_flash( + eos: $py_eos, + temperature: PySINumber, + pressure: PySINumber, + feed: &PySIArray1, + initial_state: Option<&PyPhaseEquilibrium>, + max_iter: Option, + tol: Option, + verbosity: Option, + non_volatile_components: Option>, + ) -> PyResult { + Ok(Self(PhaseEquilibrium::tp_flash( + &eos.0, + temperature.into(), + pressure.into(), + feed, + initial_state.and_then(|s| Some(&s.0)), + (max_iter, tol, verbosity).into(), non_volatile_components + )?)) + } + + /// Compute a phase equilibrium for given temperature + /// or pressure and liquid mole fractions. + /// + /// Parameters + /// ---------- + /// eos : EquationOfState + /// The equation of state. + /// temperature_or_pressure : SINumber + /// The system temperature_or_pressure. + /// liquid_molefracs : numpy.ndarray + /// The mole fraction of the liquid phase. + /// tp_init : SINumber, optional + /// The system pressure/temperature used as starting + /// condition for the iteration. + /// vapor_molefracs : numpy.ndarray, optional + /// The mole fraction of the vapor phase used as + /// starting condition for iteration. + /// max_iter_inner : int, optional + /// The maximum number of inner iterations. + /// max_iter_outer : int, optional + /// The maximum number of outer iterations. + /// tol_inner : float, optional + /// The solution tolerance in the inner loop. + /// tol_outer : float, optional + /// The solution tolerance in the outer loop. + /// verbosity : Verbosity, optional + /// The verbosity. + /// + /// Returns + /// ------- + /// PhaseEquilibrium + #[staticmethod] + #[pyo3(text_signature = "(eos, temperature_or_pressure, liquid_molefracs, tp_init=None, vapor_molefracs=None, max_iter_inner=None, max_iter_outer=None, tol_inner=None, tol_outer=None, verbosity=None)")] + pub fn bubble_point( + eos: $py_eos, + temperature_or_pressure: PySINumber, + liquid_molefracs: &PyArray1, + tp_init: Option, + vapor_molefracs: Option<&PyArray1>, + max_iter_inner: Option, + max_iter_outer: Option, + tol_inner: Option, + tol_outer: Option, + verbosity: Option, + ) -> PyResult { + let x = vapor_molefracs.and_then(|m| Some(m.to_owned_array())); + Ok(Self(PhaseEquilibrium::bubble_point( + &eos.0, + temperature_or_pressure.into(), + &liquid_molefracs.to_owned_array(), + tp_init.map(|p| p.into()), + x.as_ref(), + ( + (max_iter_inner, tol_inner, verbosity).into(), + (max_iter_outer, tol_outer, verbosity).into() + ) + )?)) + } + + /// Compute a phase equilibrium for given temperature + /// or pressure and vapor mole fractions. + /// + /// Parameters + /// ---------- + /// eos : EquationOfState + /// The equation of state. + /// temperature_or_pressure : SINumber + /// The system temperature or pressure. + /// vapor_molefracs : numpy.ndarray + /// The mole fraction of the vapor phase. + /// tp_init : SINumber, optional + /// The system pressure/temperature used as starting + /// condition for the iteration. + /// liquid_molefracs : numpy.ndarray, optional + /// The mole fraction of the liquid phase used as + /// starting condition for iteration. + /// max_iter_inner : int, optional + /// The maximum number of inner iterations. + /// max_iter_outer : int, optional + /// The maximum number of outer iterations. + /// tol_inner : float, optional + /// The solution tolerance in the inner loop. + /// tol_outer : float, optional + /// The solution tolerance in the outer loop. + /// verbosity : Verbosity, optional + /// The verbosity. + /// + /// Returns + /// ------- + /// PhaseEquilibrium + #[staticmethod] + #[pyo3(text_signature = "(eos, temperature_or_pressure, vapor_molefracs, tp_init=None, liquid_molefracs=None, max_iter_inner=None, max_iter_outer=None, tol_inner=None, tol_outer=None, verbosity=None)")] + pub fn dew_point( + eos: $py_eos, + temperature_or_pressure: PySINumber, + vapor_molefracs: &PyArray1, + tp_init: Option, + liquid_molefracs: Option<&PyArray1>, + max_iter_inner: Option, + max_iter_outer: Option, + tol_inner: Option, + tol_outer: Option, + verbosity: Option, + ) -> PyResult { + let x = liquid_molefracs.and_then(|m| Some(m.to_owned_array())); + Ok(Self(PhaseEquilibrium::dew_point( + &eos.0, + temperature_or_pressure.into(), + &vapor_molefracs.to_owned_array(), + tp_init.map(|p| p.into()), + x.as_ref(), + ( + (max_iter_inner, tol_inner, verbosity).into(), + (max_iter_outer, tol_outer, verbosity).into() + ) + )?)) + } + + #[getter] + fn get_vapor(&self) -> PyState { + PyState(self.0.vapor().clone()) + } + + #[getter] + fn get_liquid(&self) -> PyState { + PyState(self.0.liquid().clone()) + } + + /// Calculate a new PhaseEquilibrium with the given chemical potential. + /// The temperature remains constant, but the states are not in + /// a mechanical equilibrium anymore. + /// + /// Parameters + /// ---------- + /// chemical_potential: SIArray1 + /// The new chemical potential + /// + #[pyo3(text_signature = "(chemical_potential)")] + fn update_chemical_potential(slf: &PyCell, chemical_potential: &PySIArray1) -> PyResult<()> { + slf.borrow_mut().0.update_chemical_potential(chemical_potential)?; + Ok(()) + } + + /// Calculate the pure component vapor-liquid equilibria for all + /// components in the system. + /// + /// Parameters + /// ---------- + /// eos : EquationOfState + /// The equation of state. + /// temperature_or_pressure : SINumber + /// The system temperature or pressure. + /// + /// Returns + /// ------- + /// list[PhaseEquilibrium] + #[staticmethod] + #[pyo3(text_signature = "(eos, temperature_or_pressure)")] + fn vle_pure_comps(eos: $py_eos, temperature_or_pressure: PySINumber) -> Vec> { + PhaseEquilibrium::vle_pure_comps(&eos.0, temperature_or_pressure.into()) + .into_iter() + .map(|o| o.map(Self)) + .collect() + } + + /// Calculate the pure component vapor pressures for all the + /// components in the system. + /// + /// Parameters + /// ---------- + /// eos : EquationOfState + /// The equation of state. + /// temperature : SINumber + /// The system temperature. + /// + /// Returns + /// ------- + /// list[SINumber] + #[staticmethod] + #[pyo3(text_signature = "(eos, temperature)")] + fn vapor_pressure(eos: $py_eos, temperature: PySINumber) -> Vec> { + PhaseEquilibrium::vapor_pressure(&eos.0, temperature.into()) + .into_iter() + .map(|o| o.map(|n| n.into())) + .collect() + } + + /// Calculate the pure component boiling temperatures for all the + /// components in the system. + /// + /// Parameters + /// ---------- + /// eos : EquationOfState + /// The equation of state. + /// pressure : SINumber + /// The system pressure. + /// + /// Returns + /// ------- + /// list[SINumber] + #[staticmethod] + #[pyo3(text_signature = "(eos, pressure)")] + fn boiling_temperature(eos: $py_eos, pressure: PySINumber) -> Vec> { + PhaseEquilibrium::boiling_temperature(&eos.0, pressure.into()) + .into_iter() + .map(|o| o.map(|n| n.into())) + .collect() + } + + fn _repr_markdown_(&self) -> String { + self.0._repr_markdown_() + } + + fn __repr__(&self) -> PyResult { + Ok(self.0.to_string()) + } + } + + /// A thermodynamic three phase equilibrium state. + #[pyclass(name = "ThreePhaseEquilibrium")] + #[derive(Clone)] + struct PyThreePhaseEquilibrium(PhaseEquilibrium); + + #[pymethods] + impl PyPhaseEquilibrium { + /// Calculate a heteroazeotrope in a binary mixture for a given temperature + /// or pressure. + /// + /// Parameters + /// ---------- + /// eos : EquationOfState + /// The equation of state. + /// temperature_or_pressure : SINumber + /// The system temperature or pressure. + /// x_init : list[float] + /// Initial guesses for the liquid molefracs of component 1 + /// at the heteroazeotropic point. + /// tp_init : SINumber, optional + /// Initial guess for the temperature/pressure at the + /// heteroszeotropic point. + /// max_iter : int, optional + /// The maximum number of iterations. + /// tol: float, optional + /// The solution tolerance. + /// verbosity : Verbosity, optional + /// The verbosity. + /// max_iter_bd_inner : int, optional + /// The maximum number of inner iterations in the bubble/dew point iteration. + /// max_iter_bd_outer : int, optional + /// The maximum number of outer iterations in the bubble/dew point iteration. + /// tol_bd_inner : float, optional + /// The solution tolerance in the inner loop of the bubble/dew point iteration. + /// tol_bd_outer : float, optional + /// The solution tolerance in the outer loop of the bubble/dew point iteration. + /// verbosity_bd : Verbosity, optional + /// The verbosity of the bubble/dew point iteration. + #[staticmethod] + #[pyo3(text_signature = "(eos, temperature_or_pressure, x_init, tp_init=None, max_iter=None, tol=None, verbosity=None, max_iter_bd_inner=None, max_iter_bd_outer=None, tol_bd_inner=None, tol_bd_outer=None, verbosity_bd=None)")] + fn heteroazeotrope( + eos: $py_eos, + temperature_or_pressure: PySINumber, + x_init: (f64, f64), + tp_init: Option, + max_iter: Option, + tol: Option, + verbosity: Option, + max_iter_bd_inner: Option, + max_iter_bd_outer: Option, + tol_bd_inner: Option, + tol_bd_outer: Option, + verbosity_bd: Option, + ) -> PyResult { + Ok(PyThreePhaseEquilibrium(PhaseEquilibrium::heteroazeotrope( + &eos.0, + temperature_or_pressure.into(), + x_init, + tp_init.map(|t| t.into()), + (max_iter, tol, verbosity).into(), + ( + (max_iter_bd_inner, tol_bd_inner, verbosity_bd).into(), + (max_iter_bd_outer, tol_bd_outer, verbosity_bd).into(), + ) + )?)) + } + } + + #[pymethods] + impl PyThreePhaseEquilibrium { + #[getter] + fn get_vapor(&self) -> PyState { + PyState(self.0.vapor().clone()) + } + + #[getter] + fn get_liquid1(&self) -> PyState { + PyState(self.0.liquid1().clone()) + } + + #[getter] + fn get_liquid2(&self) -> PyState { + PyState(self.0.liquid2().clone()) + } + + fn _repr_markdown_(&self) -> String { + self.0._repr_markdown_() + } + + fn __repr__(&self) -> PyResult { + Ok(self.0.to_string()) + } + } + + #[pymethods] + impl PyState { + /// Calculates a two phase Tp-flash with the state as feed. + /// + /// Parameters + /// ---------- + /// initial_state : PhaseEquilibrium, optional + /// A phase equilibrium used as initial guess. + /// Can speed up convergence. + /// max_iter : int, optional + /// The maximum number of iterations. + /// tol: float, optional + /// The solution tolerance. + /// verbosity : Verbosity, optional + /// The verbosity. + /// + /// Returns + /// ------- + /// PhaseEquilibrium + /// + /// Raises + /// ------ + /// RuntimeError + /// When pressure iteration fails or no phase equilibrium is found. + #[pyo3(text_signature = "($self, initial_state=None, max_iter=None, tol=None, verbosity=None, non_volatile_components=None)")] + pub fn tp_flash( + &self, + initial_state: Option<&PyPhaseEquilibrium>, + max_iter: Option, + tol: Option, + verbosity: Option, + non_volatile_components: Option>, + ) -> PyResult { + Ok(PyPhaseEquilibrium(self.0.tp_flash( + initial_state.and_then(|s| Some(&s.0)), + (max_iter, tol, verbosity).into(), + non_volatile_components + )?)) + } + } + + /// Phase diagram for a pure component or a binary mixture. + #[pyclass(name = "PhaseDiagram")] + pub struct PyPhaseDiagram(PhaseDiagram); + + #[pymethods] + impl PyPhaseDiagram { + /// Calculate a pure component phase diagram. + /// + /// Parameters + /// ---------- + /// eos: Eos + /// The equation of state. + /// min_temperature: SINumber + /// The lower limit for the temperature. + /// npoints: int + /// The number of points. + /// 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 + #[staticmethod] + #[pyo3(text_signature = "(eos, min_temperature, npoints, critical_temperature=None, max_iter=None, tol=None, verbosity=None)")] + pub fn pure( + eos: &$py_eos, + min_temperature: PySINumber, + npoints: usize, + critical_temperature: Option, + max_iter: Option, + tol: Option, + verbosity: Option, + ) -> PyResult { + let dia = PhaseDiagram::pure( + &eos.0, + min_temperature.into(), + npoints, + critical_temperature.map(|t| t.into()), + (max_iter, tol, verbosity).into(), + )?; + 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, 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, + max_iter: Option, + tol: Option, + verbosity: Option, + ) -> EosResult { + let thread_pool = rayon::ThreadPoolBuilder::new() + .num_threads(nthreads) + .build()?; + let dia = PhaseDiagram::par_pure( + &eos.0, + min_temperature.into(), + npoints, + chunksize, + thread_pool, + critical_temperature.map(|t| t.into()), + (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 + /// bubble point line while the vapor states contain the corresponding + /// equilibrium states at different compositions. + /// + /// Parameters + /// ---------- + /// eos: Eos + /// The equation of state. + /// moles: SIArray1 + /// The moles of the individual components + /// min_temperature: SINumber + /// The lower limit for the temperature. + /// npoints: int + /// The number of points. + /// 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_inner : int, optional + /// The maximum number of inner iterations in the bubble/dew point iteration. + /// max_iter_outer : int, optional + /// The maximum number of outer iterations in the bubble/dew point iteration. + /// tol_inner : float, optional + /// The solution tolerance in the inner loop of the bubble/dew point iteration. + /// tol_outer : float, optional + /// The solution tolerance in the outer loop of the bubble/dew point iteration. + /// verbosity : Verbosity, optional + /// The verbosity of the bubble/dew point iteration. + /// + /// Returns + /// ------- + /// PhaseDiagram + #[staticmethod] + #[pyo3(text_signature = "(eos, moles, min_temperature, npoints, critical_temperature=None, max_iter_inner=None, max_iter_outer=None, tol_inner=None, tol_outer=None, verbosity=None)")] + pub fn bubble_point_line( + eos: &$py_eos, + moles: PySIArray1, + min_temperature: PySINumber, + npoints: usize, + critical_temperature: Option, + max_iter_inner: Option, + max_iter_outer: Option, + tol_inner: Option, + tol_outer: Option, + verbosity: Option, + ) -> PyResult { + let dia = PhaseDiagram::bubble_point_line( + &eos.0, + &moles, + min_temperature.into(), + npoints, + critical_temperature.map(|t| t.into()), + ( + (max_iter_inner, tol_inner, verbosity).into(), + (max_iter_outer, tol_outer, verbosity).into(), + ) + )?; + Ok(Self(dia)) + } + + /// Calculate the dew point line of a mixture with given composition. + /// + /// In the resulting phase diagram, the vapor states correspond to the + /// dew point line while the liquid states contain the corresponding + /// equilibrium states at different compositions. + /// + /// Parameters + /// ---------- + /// eos: Eos + /// The equation of state. + /// moles: SIArray1 + /// The moles of the individual components + /// min_temperature: SINumber + /// The lower limit for the temperature. + /// npoints: int + /// The number of points. + /// 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_inner : int, optional + /// The maximum number of inner iterations in the bubble/dew point iteration. + /// max_iter_outer : int, optional + /// The maximum number of outer iterations in the bubble/dew point iteration. + /// tol_inner : float, optional + /// The solution tolerance in the inner loop of the bubble/dew point iteration. + /// tol_outer : float, optional + /// The solution tolerance in the outer loop of the bubble/dew point iteration. + /// verbosity : Verbosity, optional + /// The verbosity of the bubble/dew point iteration. + /// + /// Returns + /// ------- + /// PhaseDiagram + #[staticmethod] + #[pyo3(text_signature = "(eos, moles, min_temperature, npoints, critical_temperature=None, max_iter_inner=None, max_iter_outer=None, tol_inner=None, tol_outer=None, verbosity=None)")] + pub fn dew_point_line( + eos: &$py_eos, + moles: PySIArray1, + min_temperature: PySINumber, + npoints: usize, + critical_temperature: Option, + max_iter_inner: Option, + max_iter_outer: Option, + tol_inner: Option, + tol_outer: Option, + verbosity: Option, + ) -> PyResult { + let dia = PhaseDiagram::dew_point_line( + &eos.0, + &moles, + min_temperature.into(), + npoints, + critical_temperature.map(|t| t.into()), + ( + (max_iter_inner, tol_inner, verbosity).into(), + (max_iter_outer, tol_outer, verbosity).into(), + ) + )?; + Ok(Self(dia)) + } + + /// Calculate the spinodal lines for a mixture with fixed composition. + /// + /// Parameters + /// ---------- + /// eos: Eos + /// The equation of state. + /// moles: SIArray1 + /// The moles of the individual components + /// min_temperature: SINumber + /// The lower limit for the temperature. + /// npoints: int + /// The number of points. + /// 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 + #[staticmethod] + #[pyo3(text_signature = "(eos, moles, min_temperature, npoints, critical_temperature=None, max_iter=None, tol=None, verbosity=None)")] + pub fn spinodal( + eos: &$py_eos, + moles: PySIArray1, + min_temperature: PySINumber, + npoints: usize, + critical_temperature: Option, + max_iter: Option, + tol: Option, + verbosity: Option, + ) -> PyResult { + let dia = PhaseDiagram::spinodal( + &eos.0, + &moles, + min_temperature.into(), + npoints, + critical_temperature.map(|t| t.into()), + (max_iter, tol, verbosity).into(), + )?; + Ok(Self(dia)) + } + + #[getter] + pub fn get_states(&self) -> Vec { + self.0 + .states + .iter() + .map(|vle| PyPhaseEquilibrium(vle.clone())) + .collect() + } + + #[getter] + pub fn get_vapor(&self) -> PyStateVec { + self.0.vapor().into() + } + + #[getter] + pub fn get_liquid(&self) -> PyStateVec { + self.0.liquid().into() + } + + /// Returns the phase diagram as dictionary. + /// + /// Units + /// ----- + /// temperature : K + /// pressure : Pa + /// densities : mol / m³ + /// molar enthalpies : kJ / mol + /// molar entropies : kJ / mol / K + /// + /// Returns + /// ------- + /// dict[str, list[float]] + /// Keys: property names. Values: property for each state. + /// + /// Notes + /// ----- + /// xi: liquid molefraction of component i + /// yi: vapor molefraction of component i + /// i: component index according to order in parameters. + pub fn to_dict(&self) -> PyResult>> { + let n = self.0.states[0].liquid().eos.components(); + let mut dict = HashMap::with_capacity(8 + 2 * n); + if n != 1 { + let xs = self.0.liquid().molefracs(); + let ys = self.0.vapor().molefracs(); + for i in 0..n { + dict.insert(String::from(format!("x{}", i)), xs.column(i).to_vec()); + dict.insert(String::from(format!("y{}", i)), ys.column(i).to_vec()); + } + } + dict.insert(String::from("temperature"), (self.0.vapor().temperature() / KELVIN).into_value()?.into_raw_vec()); + dict.insert(String::from("pressure"), (self.0.vapor().pressure() / PASCAL).into_value()?.into_raw_vec()); + dict.insert(String::from("density liquid"), (self.0.liquid().density() / (MOL / METER.powi(3))).into_value()?.into_raw_vec()); + dict.insert(String::from("density vapor"), (self.0.vapor().density() / (MOL / METER.powi(3))).into_value()?.into_raw_vec()); + dict.insert(String::from("molar enthalpy liquid"), (self.0.liquid().molar_enthalpy() / (KILO*JOULE / MOL)).into_value()?.into_raw_vec()); + dict.insert(String::from("molar enthalpy vapor"), (self.0.vapor().molar_enthalpy() / (KILO*JOULE / MOL)).into_value()?.into_raw_vec()); + dict.insert(String::from("molar entropy liquid"), (self.0.liquid().molar_entropy() / (KILO*JOULE / KELVIN / MOL)).into_value()?.into_raw_vec()); + dict.insert(String::from("molar entropy vapor"), (self.0.vapor().molar_entropy() / (KILO*JOULE / KELVIN / MOL)).into_value()?.into_raw_vec()); + Ok(dict) + } + + /// Binary phase diagram calculated using bubble/dew point iterations. + /// + /// Parameters + /// ---------- + /// eos : EquationOfState + /// The equation of state. + /// temperature_or_pressure: SINumber + /// The constant temperature or pressure. + /// npoints: int, optional + /// The number of points (default 51). + /// x_lle: (float, float), optional + /// An estimate for the molefractions of component 1 + /// at the heteroazeotrop + /// max_iter_inner : int, optional + /// The maximum number of inner iterations in the bubble/dew point iteration. + /// max_iter_outer : int, optional + /// The maximum number of outer iterations in the bubble/dew point iteration. + /// tol_inner : float, optional + /// The solution tolerance in the inner loop of the bubble/dew point iteration. + /// tol_outer : float, optional + /// The solution tolerance in the outer loop of the bubble/dew point iteration. + /// verbosity : Verbosity, optional + /// The verbosity of the bubble/dew point iteration. + /// + /// Returns + /// ------- + /// PhaseDiagram + #[staticmethod] + #[pyo3(text_signature = "(eos, temperature_or_pressure, npoints=None, x_lle=None, max_iter_inner=None, max_iter_outer=None, tol_inner=None, tol_outer=None, verbosity=None)")] + pub fn binary_vle( + eos: $py_eos, + temperature_or_pressure: PySINumber, + npoints: Option, + x_lle: Option<(f64, f64)>, + max_iter_inner: Option, + max_iter_outer: Option, + tol_inner: Option, + tol_outer: Option, + verbosity: Option, + ) -> PyResult { + let dia = PhaseDiagram::binary_vle( + &eos.0, + temperature_or_pressure.into(), + npoints, + x_lle, + ( + (max_iter_inner, tol_inner, verbosity).into(), + (max_iter_outer, tol_outer, verbosity).into(), + ) + )?; + Ok(Self(dia)) + } + + /// Create a new phase diagram using Tp flash calculations. + /// + /// The usual use case for this function is the calculation of + /// liquid-liquid phase diagrams, but it can be used for vapor- + /// liquid diagrams as well, as long as the feed composition is + /// in a two phase region. + /// + /// Parameters + /// ---------- + /// eos : EquationOfState + /// The equation of state. + /// temperature_or_pressure: SINumber + /// The consant temperature or pressure. + /// feed: SIArray1 + /// Mole numbers in the (unstable) feed state. + /// min_tp: + /// The lower limit of the temperature/pressure range. + /// max_tp: + /// The upper limit of the temperature/pressure range. + /// npoints: int, optional + /// The number of points (default 51). + /// + /// Returns + /// ------- + /// PhaseDiagram + #[staticmethod] + #[pyo3(text_signature = "(eos, temperature_or_pressure, feed, min_tp, max_tp, npoints=None)")] + pub fn lle( + eos: $py_eos, + temperature_or_pressure: PySINumber, + feed: PySIArray1, + min_tp: PySINumber, + max_tp: PySINumber, + npoints: Option, + ) -> PyResult { + let dia = PhaseDiagram::lle( + &eos.0, + temperature_or_pressure.into(), + &feed, + min_tp.into(), + max_tp.into(), + npoints, + )?; + Ok(Self(dia)) + } + } + + /// Phase diagram for a binary mixture exhibiting a heteroazeotrope. + #[pyclass(name = "PhaseDiagramHetero")] + pub struct PyPhaseDiagramHetero(PhaseDiagramHetero); + + #[pymethods] + impl PyPhaseDiagram { + /// Phase diagram for a binary mixture exhibiting a heteroazeotrope. + /// + /// Parameters + /// ---------- + /// eos: SaftFunctional + /// The SAFT Helmholtz energy functional. + /// temperature_or_pressure: SINumber + /// The temperature_or_pressure. + /// x_lle: SINumber + /// Initial values for the molefractions of component 1 + /// at the heteroazeotrop. + /// tp_lim_lle: SINumber, optional + /// The minimum temperature up to which the LLE is calculated. + /// If it is not provided, no LLE is calcualted. + /// tp_init_vlle: SINumber, optional + /// Initial value for the calculation of the VLLE. + /// npoints_vle: int, optional + /// The number of points for the VLE (default 51). + /// npoints_lle: int, optional + /// The number of points for the LLE (default 51). + /// max_iter_inner : int, optional + /// The maximum number of inner iterations in the bubble/dew point iteration. + /// max_iter_outer : int, optional + /// The maximum number of outer iterations in the bubble/dew point iteration. + /// tol_inner : float, optional + /// The solution tolerance in the inner loop of the bubble/dew point iteration. + /// tol_outer : float, optional + /// The solution tolerance in the outer loop of the bubble/dew point iteration. + /// verbosity : Verbosity, optional + /// The verbosity of the bubble/dew point iteration. + /// + /// Returns + /// ------- + /// PhaseDiagramHetero + #[staticmethod] + #[pyo3(text_signature = "(eos, temperature_or_pressure, x_lle, tp_lim_lle=None, tp_init_vlle=None, npoints_vle=None, npoints_lle=None, max_iter_bd_inner=None, max_iter_bd_outer=None, tol_bd_inner=None, tol_bd_outer=None, verbosity_bd=None)")] + pub fn binary_vlle( + eos: $py_eos, + temperature_or_pressure: PySINumber, + x_lle: (f64, f64), + tp_lim_lle: Option, + tp_init_vlle: Option, + npoints_vle: Option, + npoints_lle: Option, + max_iter_inner: Option, + max_iter_outer: Option, + tol_inner: Option, + tol_outer: Option, + verbosity: Option, + ) -> PyResult { + let dia = PhaseDiagram::binary_vlle( + &eos.0, + temperature_or_pressure.into(), + x_lle, + tp_lim_lle.map(|t| t.into()), + tp_init_vlle.map(|t| t.into()), + npoints_vle, + npoints_lle, + ( + (max_iter_inner, tol_inner, verbosity).into(), + (max_iter_outer, tol_outer, verbosity).into(), + ) + )?; + Ok(PyPhaseDiagramHetero(dia)) + } + } + + #[pymethods] + impl PyPhaseDiagramHetero { + #[getter] + pub fn get_vle(&self) -> PyPhaseDiagram { + PyPhaseDiagram(self.0.vle().clone()) + } + + #[getter] + pub fn get_vle1(&self) -> PyPhaseDiagram { + PyPhaseDiagram(self.0.vle1.clone()) + } + + #[getter] + pub fn get_vle2(&self) -> PyPhaseDiagram { + PyPhaseDiagram(self.0.vle2.clone()) + } + + #[getter] + pub fn get_lle(&self) -> Option { + self.0 + .lle + .as_ref() + .map(|d| PyPhaseDiagram(d.clone())) + } + } + } +} diff --git a/feos-core/src/python/state.rs b/feos-core/src/python/state.rs index 6d36e512e..3c88e4001 100644 --- a/feos-core/src/python/state.rs +++ b/feos-core/src/python/state.rs @@ -46,7 +46,7 @@ macro_rules! impl_state { /// ------ /// Error /// When the state cannot be created using the combination of input. - #[pyclass(name = "State", unsendable)] + #[pyclass(name = "State")] #[derive(Clone)] #[pyo3(text_signature = "(eos, temperature=None, volume=None, density=None, partial_density=None, total_moles=None, moles=None, molefracs=None, pressure=None, molar_enthalpy=None, molar_entropy=None, molar_internal_energy=None, density_initialization=None, initial_temperature=None)")] pub struct PyState(pub State); @@ -1012,7 +1012,7 @@ macro_rules! impl_state { } - #[pyclass(name = "StateVec", unsendable)] + #[pyclass(name = "StateVec")] pub struct PyStateVec(Vec>); impl From> for PyStateVec { diff --git a/feos-core/src/state/builder.rs b/feos-core/src/state/builder.rs index 52af1e454..406cc3d04 100644 --- a/feos-core/src/state/builder.rs +++ b/feos-core/src/state/builder.rs @@ -4,7 +4,7 @@ use crate::errors::EosResult; use crate::EosUnit; use ndarray::Array1; use quantity::{QuantityArray1, QuantityScalar}; -use std::rc::Rc; +use std::sync::Arc; /// A simple tool to construct [State]s with arbitrary input parameters. /// @@ -13,12 +13,12 @@ use std::rc::Rc; /// # use feos_core::{EosResult, StateBuilder}; /// # use feos_core::cubic::{PengRobinson, PengRobinsonParameters}; /// # use quantity::si::*; -/// # use std::rc::Rc; +/// # use std::sync::Arc; /// # use ndarray::arr1; /// # use approx::assert_relative_eq; /// # fn main() -> EosResult<()> { /// // Create a state for given T,V,N -/// let eos = Rc::new(PengRobinson::new(Rc::new(PengRobinsonParameters::new_simple(&[369.8], &[41.9 * 1e5], &[0.15], &[15.0])?))); +/// let eos = Arc::new(PengRobinson::new(Arc::new(PengRobinsonParameters::new_simple(&[369.8], &[41.9 * 1e5], &[0.15], &[15.0])?))); /// let state = StateBuilder::new(&eos) /// .temperature(300.0 * KELVIN) /// .volume(12.5 * METER.powi(3)) @@ -27,7 +27,7 @@ use std::rc::Rc; /// assert_eq!(state.density, 0.2 * MOL / METER.powi(3)); /// /// // For a pure component, the composition does not need to be specified. -/// let eos = Rc::new(PengRobinson::new(Rc::new(PengRobinsonParameters::new_simple(&[369.8], &[41.9 * 1e5], &[0.15], &[15.0])?))); +/// let eos = Arc::new(PengRobinson::new(Arc::new(PengRobinsonParameters::new_simple(&[369.8], &[41.9 * 1e5], &[0.15], &[15.0])?))); /// let state = StateBuilder::new(&eos) /// .temperature(300.0 * KELVIN) /// .volume(12.5 * METER.powi(3)) @@ -36,8 +36,8 @@ use std::rc::Rc; /// assert_eq!(state.density, 0.2 * MOL / METER.powi(3)); /// /// // The state can be constructed without providing any extensive property. -/// let eos = Rc::new(PengRobinson::new( -/// Rc::new(PengRobinsonParameters::new_simple( +/// let eos = Arc::new(PengRobinson::new( +/// Arc::new(PengRobinsonParameters::new_simple( /// &[369.8, 305.4], /// &[41.9 * 1e5, 48.2 * 1e5], /// &[0.15, 0.10], @@ -54,7 +54,7 @@ use std::rc::Rc; /// # } /// ``` pub struct StateBuilder<'a, U: EosUnit, E: EquationOfState> { - eos: Rc, + eos: Arc, temperature: Option>, volume: Option>, density: Option>, @@ -72,7 +72,7 @@ pub struct StateBuilder<'a, U: EosUnit, E: EquationOfState> { impl<'a, U: EosUnit, E: EquationOfState> StateBuilder<'a, U, E> { /// Create a new `StateBuilder` for the given equation of state. - pub fn new(eos: &Rc) -> Self { + pub fn new(eos: &Arc) -> Self { StateBuilder { eos: eos.clone(), temperature: None, diff --git a/feos-core/src/state/critical_point.rs b/feos-core/src/state/critical_point.rs index 34f01d8a3..2976b6f59 100644 --- a/feos-core/src/state/critical_point.rs +++ b/feos-core/src/state/critical_point.rs @@ -9,7 +9,7 @@ use num_dual::{Dual, Dual3, Dual64, DualNum, DualVec64, HyperDual, StaticVec}; use num_traits::{One, Zero}; use quantity::{QuantityArray1, QuantityScalar}; use std::convert::TryFrom; -use std::rc::Rc; +use std::sync::Arc; const MAX_ITER_CRIT_POINT: usize = 50; const MAX_ITER_CRIT_POINT_BINARY: usize = 200; @@ -19,7 +19,7 @@ const TOL_CRIT_POINT: f64 = 1e-8; impl State { /// Calculate the pure component critical point of all components. pub fn critical_point_pure( - eos: &Rc, + eos: &Arc, initial_temperature: Option>, options: SolverOptions, ) -> EosResult> @@ -29,7 +29,7 @@ impl State { (0..eos.components()) .map(|i| { Self::critical_point( - &Rc::new(eos.subset(&[i])), + &Arc::new(eos.subset(&[i])), None, initial_temperature, options, @@ -39,7 +39,7 @@ impl State { } pub fn critical_point_binary( - eos: &Rc, + eos: &Arc, temperature_or_pressure: QuantityScalar, initial_temperature: Option>, initial_molefracs: Option<[f64; 2]>, @@ -64,7 +64,7 @@ impl State { /// Calculate the critical point of a system for given moles. pub fn critical_point( - eos: &Rc, + eos: &Arc, moles: Option<&QuantityArray1>, initial_temperature: Option>, options: SolverOptions, @@ -91,7 +91,7 @@ impl State { } fn critical_point_hkm( - eos: &Rc, + eos: &Arc, moles: &QuantityArray1, initial_temperature: QuantityScalar, options: SolverOptions, @@ -178,7 +178,7 @@ impl State { /// Calculate the critical point of a binary system for given temperature. fn critical_point_binary_t( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, initial_molefracs: Option<[f64; 2]>, options: SolverOptions, @@ -264,7 +264,7 @@ impl State { /// Calculate the critical point of a binary system for given pressure. fn critical_point_binary_p( - eos: &Rc, + eos: &Arc, pressure: QuantityScalar, initial_temperature: Option>, initial_molefracs: Option<[f64; 2]>, @@ -364,7 +364,7 @@ impl State { } pub fn spinodal( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, moles: Option<&QuantityArray1>, options: SolverOptions, @@ -393,7 +393,7 @@ impl State { } fn calculate_spinodal( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, moles: &QuantityArray1, density_initialization: DensityInitialization, @@ -469,7 +469,7 @@ impl State { } fn critical_point_objective( - eos: &Rc, + eos: &Arc, temperature: Dual64, density: Dual64, moles: &Array1, @@ -509,7 +509,7 @@ fn critical_point_objective( } fn critical_point_objective_t( - eos: &Rc, + eos: &Arc, temperature: f64, density: StaticVec, 2>, ) -> EosResult, 2>> { @@ -544,7 +544,7 @@ fn critical_point_objective_t( } fn critical_point_objective_p( - eos: &Rc, + eos: &Arc, pressure: f64, temperature: DualVec64<3>, density: StaticVec, 2>, @@ -591,7 +591,7 @@ fn critical_point_objective_p( } fn spinodal_objective( - eos: &Rc, + eos: &Arc, temperature: Dual64, density: Dual64, moles: &Array1, diff --git a/feos-core/src/state/mod.rs b/feos-core/src/state/mod.rs index 8baad3097..5b061686c 100644 --- 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::rc::Rc; +use std::sync::{Arc, Mutex}; mod builder; mod cache; @@ -123,7 +122,7 @@ impl> StateHD { #[derive(Debug)] pub struct State { /// Equation of state - pub eos: Rc, + pub eos: Arc, /// Temperature $T$ pub temperature: QuantityScalar, /// Volume $V$ @@ -145,7 +144,7 @@ pub struct State { /// Reduced moles reduced_moles: Array1, /// Cache - cache: RefCell, + cache: Mutex, } impl Clone for State { @@ -162,7 +161,7 @@ impl Clone for State { 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()), } } } @@ -224,7 +223,7 @@ impl State { /// and if values are finite. It will **not** validate physics, i.e. if the resulting /// densities are below the maximum packing fraction. pub fn new_nvt( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, volume: QuantityScalar, moles: &QuantityArray1, @@ -236,7 +235,7 @@ impl State { } pub(super) fn new_nvt_unchecked( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, volume: QuantityScalar, moles: &QuantityArray1, @@ -262,7 +261,7 @@ impl State { 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())), } } @@ -273,7 +272,7 @@ impl State { /// and if values are finite. It will **not** validate physics, i.e. if the resulting /// densities are below the maximum packing fraction. pub fn new_pure( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, density: QuantityScalar, ) -> EosResult { @@ -296,7 +295,7 @@ impl State { /// /// When the state cannot be created using the combination of inputs. pub fn new( - eos: &Rc, + eos: &Arc, temperature: Option>, volume: Option>, density: Option>, @@ -428,7 +427,7 @@ impl State { /// Return a new `State` using a density iteration. [DensityInitialization] is used to /// influence the calculation with respect to the possible solutions. pub fn new_npt( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, pressure: QuantityScalar, moles: &QuantityArray1, @@ -495,7 +494,7 @@ impl State { /// Return a new `State` for given pressure $p$, volume $V$, temperature $T$ and composition $x_i$. pub fn new_npvx( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, pressure: QuantityScalar, volume: QuantityScalar, @@ -510,7 +509,7 @@ impl State { /// Return a new `State` for given pressure $p$ and molar enthalpy $h$. pub fn new_nph( - eos: &Rc, + eos: &Arc, pressure: QuantityScalar, molar_enthalpy: QuantityScalar, moles: &QuantityArray1, @@ -531,7 +530,7 @@ impl State { /// Return a new `State` for given temperature $T$ and molar enthalpy $h$. pub fn new_nth( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, molar_enthalpy: QuantityScalar, moles: &QuantityArray1, @@ -558,7 +557,7 @@ impl State { /// Return a new `State` for given temperature $T$ and molar entropy $s$. pub fn new_nts( - eos: &Rc, + eos: &Arc, temperature: QuantityScalar, molar_entropy: QuantityScalar, moles: &QuantityArray1, @@ -582,7 +581,7 @@ impl State { /// Return a new `State` for given pressure $p$ and molar entropy $s$. pub fn new_nps( - eos: &Rc, + eos: &Arc, pressure: QuantityScalar, molar_entropy: QuantityScalar, moles: &QuantityArray1, @@ -603,7 +602,7 @@ impl State { /// Return a new `State` for given volume $V$ and molar internal energy $u$. pub fn new_nvu( - eos: &Rc, + eos: &Arc, volume: QuantityScalar, molar_internal_energy: QuantityScalar, moles: &QuantityArray1, diff --git a/feos-core/src/state/properties.rs b/feos-core/src/state/properties.rs index 13d0ea400..e84be6645 100644 --- a/feos-core/src/state/properties.rs +++ b/feos-core/src/state/properties.rs @@ -7,7 +7,7 @@ use num_dual::DualNum; use quantity::{QuantityArray, QuantityArray1, QuantityArray2, QuantityScalar}; use std::iter::FromIterator; use std::ops::{Add, Deref, Sub}; -use std::rc::Rc; +use std::sync::Arc; #[derive(Clone, Copy)] pub(crate) enum Evaluate { @@ -64,7 +64,7 @@ impl State { }; } - let mut cache = self.cache.borrow_mut(); + let mut cache = self.cache.lock().unwrap(); let residual = match evaluate { Evaluate::IdealGas => None, @@ -305,7 +305,7 @@ impl State { let pressure = self.pressure(Contributions::Total); (0..self.eos.components()) .map(|i| { - let eos = Rc::new(self.eos.subset(&[i])); + let eos = Arc::new(self.eos.subset(&[i])); let state = Self::new_npt( &eos, self.temperature, diff --git a/feos-dft/src/adsorption/mod.rs b/feos-dft/src/adsorption/mod.rs index e03d005bd..3a6ac50f2 100644 --- 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 ::Larger: Dimension, { fn new>( - functional: &Rc>, + functional: &Arc>, pore: &S, profiles: Vec>>, ) -> Self { @@ -62,7 +62,7 @@ where /// Calculate an adsorption isotherm (starting at low pressure) pub fn adsorption_isotherm>( - functional: &Rc>, + functional: &Arc>, temperature: QuantityScalar, pressure: &QuantityArray1, pore: &S, @@ -82,7 +82,7 @@ where /// Calculate an desorption isotherm (starting at high pressure) pub fn desorption_isotherm>( - functional: &Rc>, + functional: &Arc>, temperature: QuantityScalar, pressure: &QuantityArray1, pore: &S, @@ -108,7 +108,7 @@ where /// Calculate an equilibrium isotherm pub fn equilibrium_isotherm>( - functional: &Rc>, + functional: &Arc>, temperature: QuantityScalar, pressure: &QuantityArray1, pore: &S, @@ -195,7 +195,7 @@ where } fn isotherm>( - functional: &Rc>, + functional: &Arc>, temperature: QuantityScalar, pressure: &QuantityArray1, pore: &S, @@ -257,7 +257,7 @@ where /// Calculate the phase transition from an empty to a filled pore. pub fn phase_equilibrium>( - functional: &Rc>, + functional: &Arc>, temperature: QuantityScalar, p_min: QuantityScalar, p_max: QuantityScalar, diff --git a/feos-dft/src/adsorption/pore.rs b/feos-dft/src/adsorption/pore.rs index 4bbc20a48..2a384d8cc 100644 --- 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 { where D::Larger: Dimension, { - 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/convolver/mod.rs b/feos-dft/src/convolver/mod.rs index 9c49b0920..a558145b7 100644 --- a/feos-dft/src/convolver/mod.rs +++ b/feos-dft/src/convolver/mod.rs @@ -5,7 +5,7 @@ use ndarray::{Axis as Axis_nd, RemoveAxis, ScalarOperand, Slice}; use num_dual::*; use rustdct::DctNum; use std::ops::{AddAssign, MulAssign, SubAssign}; -use std::rc::Rc; +use std::sync::Arc; mod periodic_convolver; mod transform; @@ -19,7 +19,7 @@ use transform::*; /// Helmholtz energy functional. /// /// Parametrized over data types `T` and dimension of the problem `D`. -pub trait Convolver { +pub trait Convolver: Send + Sync { /// Convolve the profile with the given weight function. fn convolve(&self, profile: Array, weight_function: &WeightFunction) -> Array; @@ -79,9 +79,9 @@ pub struct ConvolverFFT { /// Lanczos sigma factor lanczos_sigma: Option>, /// Possibly curvilinear Fourier transform in the first dimension - transform: Rc>, + transform: Arc>, /// Vector of additional cartesian Fourier transforms in the other dimensions - cartesian_transforms: Vec>>, + cartesian_transforms: Vec>>, } impl ConvolverFFT @@ -96,7 +96,7 @@ where grid: &Grid, weight_functions: &[WeightFunctionInfo], lanczos: Option, - ) -> Rc> { + ) -> Arc> { match grid { Grid::Polar(r) => CurvilinearConvolver::new(r, &[], weight_functions, lanczos), Grid::Spherical(r) => CurvilinearConvolver::new(r, &[], weight_functions, lanczos), @@ -125,7 +125,7 @@ where cartesian_axes: &[&Axis], weight_functions: &[WeightFunctionInfo], lanczos: Option, - ) -> Rc> { + ) -> Arc> { // initialize the Fourier transform let mut cartesian_transforms = Vec::with_capacity(cartesian_axes.len()); let mut k_vec = Vec::with_capacity(cartesian_axes.len() + 1); @@ -222,7 +222,7 @@ where } // Return `FFTConvolver` - Rc::new(Self { + Arc::new(Self { k_abs, weight_functions: fft_weight_functions, lanczos_sigma, @@ -506,8 +506,8 @@ where /// The curvilinear convolver accounts for the shift that has to be performed /// for spherical and polar transforms. struct CurvilinearConvolver { - convolver: Rc>, - convolver_boundary: Rc>, + convolver: Arc>, + convolver_boundary: Arc>, } impl CurvilinearConvolver @@ -522,8 +522,8 @@ where z: &[&Axis], weight_functions: &[WeightFunctionInfo], lanczos: Option, - ) -> Rc> { - Rc::new(Self { + ) -> Arc> { + Arc::new(Self { convolver: ConvolverFFT::new(Some(r), z, weight_functions, lanczos), convolver_boundary: ConvolverFFT::new(None, z, weight_functions, lanczos), }) diff --git a/feos-dft/src/convolver/periodic_convolver.rs b/feos-dft/src/convolver/periodic_convolver.rs index 1d2ae1cfc..cc8a623b7 100644 --- a/feos-dft/src/convolver/periodic_convolver.rs +++ b/feos-dft/src/convolver/periodic_convolver.rs @@ -8,7 +8,6 @@ use rustfft::num_complex::Complex; use rustfft::{Fft, FftDirection, FftNum, FftPlanner}; use std::f64::consts::PI; use std::ops::AddAssign; -use std::rc::Rc; use std::sync::Arc; #[derive(Clone)] @@ -35,7 +34,7 @@ where axes: &[&Axis], weight_functions: &[WeightFunctionInfo], lanczos: Option, - ) -> Rc> { + ) -> Arc> { // initialize the Fourier transform let mut planner = FftPlanner::new(); let mut forward_transforms = Vec::with_capacity(axes.len()); @@ -130,7 +129,7 @@ where }); } - Rc::new(Self { + Arc::new(Self { k_abs, weight_functions: fft_weight_functions, lanczos_sigma, diff --git a/feos-dft/src/convolver/transform.rs b/feos-dft/src/convolver/transform.rs index 3ec46f8c6..de4a2379a 100644 --- a/feos-dft/src/convolver/transform.rs +++ b/feos-dft/src/convolver/transform.rs @@ -6,7 +6,6 @@ use rustdct::{DctNum, DctPlanner, TransformType2And3}; use rustfft::{num_complex::Complex, Fft, FftPlanner}; use std::f64::consts::PI; use std::ops::{DivAssign, SubAssign}; -use std::rc::Rc; use std::sync::Arc; #[derive(Clone, Copy)] @@ -26,7 +25,7 @@ impl SinCosTransform { } } -pub(super) trait FourierTransform> { +pub(super) trait FourierTransform>: Send + Sync { fn forward_transform(&self, f_r: ArrayView1, f_k: ArrayViewMut1, scalar: bool); fn back_transform(&self, f_k: ArrayViewMut1, f_r: ArrayViewMut1, scalar: bool); @@ -37,14 +36,14 @@ pub(super) struct CartesianTransform { } impl + DctNum + ScalarOperand> CartesianTransform { - pub(super) fn new(axis: &Axis) -> (Rc>, Array1) { + pub(super) fn new(axis: &Axis) -> (Arc>, Array1) { let (s, k) = Self::init(axis); - (Rc::new(s), k) + (Arc::new(s), k) } - pub(super) fn new_cartesian(axis: &Axis) -> (Rc, Array1) { + pub(super) fn new_cartesian(axis: &Axis) -> (Arc, Array1) { let (s, k) = Self::init(axis); - (Rc::new(s), k) + (Arc::new(s), k) } fn init(axis: &Axis) -> (Self, Array1) { @@ -130,12 +129,12 @@ pub(super) struct SphericalTransform { } impl + DctNum + ScalarOperand> SphericalTransform { - pub(super) fn new(axis: &Axis) -> (Rc>, Array1) { + pub(super) fn new(axis: &Axis) -> (Arc>, Array1) { let points = axis.grid.len(); let length = axis.length(); let k_grid: Array1<_> = (0..=points).map(|v| PI * v as f64 / length).collect(); ( - Rc::new(Self { + Arc::new(Self { r_grid: axis.grid.clone(), k_grid: k_grid.clone(), dct: DctPlanner::new().plan_dct2(points), @@ -224,7 +223,7 @@ pub(super) struct PolarTransform { } impl + DctNum + ScalarOperand> PolarTransform { - pub(super) fn new(axis: &Axis) -> (Rc>, Array1) { + pub(super) fn new(axis: &Axis) -> (Arc>, Array1) { let points = axis.grid.len(); let mut alpha = 0.002_f64; @@ -263,7 +262,7 @@ impl + DctNum + ScalarOperand> PolarTransform { ifft.process(jv.as_slice_mut().unwrap()); ( - Rc::new(Self { + Arc::new(Self { r_grid: axis.grid.clone(), k_grid: k_grid.clone(), fft, @@ -329,8 +328,8 @@ impl + DctNum + ScalarOperand> FourierTransform for PolarTran pub(super) struct NoTransform(); impl NoTransform { - pub(super) fn new>() -> (Rc>, Array1) { - (Rc::new(Self()), arr1(&[0.0])) + pub(super) fn new>() -> (Arc>, Array1) { + (Arc::new(Self()), arr1(&[0.0])) } } diff --git a/feos-dft/src/functional.rs b/feos-dft/src/functional.rs index 33b9fcc2b..ea07db2da 100644 --- a/feos-dft/src/functional.rs +++ b/feos-dft/src/functional.rs @@ -15,7 +15,7 @@ use quantity::{QuantityArray, QuantityArray1, QuantityScalar}; use std::borrow::Cow; use std::fmt; use std::ops::{AddAssign, Deref, MulAssign}; -use std::rc::Rc; +use std::sync::Arc; /// Wrapper struct for the [HelmholtzEnergyFunctional] trait. /// @@ -132,7 +132,7 @@ pub enum MoleculeShape<'a> { } /// A general Helmholtz energy functional. -pub trait HelmholtzEnergyFunctional: Sized { +pub trait HelmholtzEnergyFunctional: Sized + Send + Sync { /// Return a slice of [FunctionalContribution]s. fn contributions(&self) -> &[Box]; @@ -203,7 +203,7 @@ impl DFT { &self, temperature: QuantityScalar, density: &QuantityArray, - convolver: &Rc>, + convolver: &Arc>, ) -> EosResult> where U: EosUnit, @@ -275,7 +275,7 @@ impl DFT { &self, temperature: N, density: &Array, - convolver: &Rc>, + convolver: &Arc>, ) -> EosResult> where N: DualNum + ScalarOperand, @@ -311,7 +311,7 @@ impl DFT { &self, temperature: f64, density: &Array, - convolver: &Rc>, + convolver: &Arc>, contributions: Contributions, ) -> EosResult> where @@ -338,7 +338,7 @@ impl DFT { &self, temperature: f64, density: &Array, - convolver: &Rc>, + convolver: &Arc>, ) -> EosResult>> where D: Dimension, @@ -382,7 +382,7 @@ impl DFT { temperature: f64, density: &Array, external_potential: &Array, - convolver: &Rc>, + convolver: &Arc>, contributions: Contributions, ) -> EosResult> where @@ -411,7 +411,7 @@ impl DFT { &self, temperature: f64, density: &Array, - convolver: &Rc>, + convolver: &Arc>, ) -> EosResult<(Array, Array)> where D: Dimension, @@ -446,7 +446,7 @@ impl DFT { &self, temperature: f64, functional_derivative: &Array, - convolver: &Rc>, + convolver: &Arc>, ) -> Array where D: Dimension, diff --git a/feos-dft/src/functional_contribution.rs b/feos-dft/src/functional_contribution.rs index b07882015..7621c63e2 100644 --- a/feos-dft/src/functional_contribution.rs +++ b/feos-dft/src/functional_contribution.rs @@ -84,6 +84,8 @@ pub trait FunctionalContribution: + FunctionalContributionDual, f64>> + FunctionalContributionDual, f64>> + Display + + Sync + + Send { fn first_partial_derivatives( &self, @@ -167,5 +169,7 @@ impl FunctionalContribution for T where + FunctionalContributionDual, f64>> + FunctionalContributionDual, f64>> + Display + + Sync + + Send { } diff --git a/feos-dft/src/profile.rs b/feos-dft/src/profile.rs index 870b2b142..af24b20f4 100644 --- a/feos-dft/src/profile.rs +++ b/feos-dft/src/profile.rs @@ -13,7 +13,7 @@ use ndarray::{ use num_dual::Dual64; use quantity::{Quantity, QuantityArray, QuantityArray1, QuantityScalar}; use std::ops::MulAssign; -use std::rc::Rc; +use std::sync::Arc; pub(crate) const MAX_POTENTIAL: f64 = 50.0; #[cfg(feature = "3d_dft")] @@ -24,7 +24,7 @@ pub(crate) const CUTOFF_RADIUS: f64 = 14.0; /// In the most basic case, the chemical potential is specified in a DFT calculation, /// for more general systems, this trait provides the possibility to declare additional /// equations for the calculation of the chemical potential during the iteration. -pub trait DFTSpecification { +pub trait DFTSpecification: Send + Sync { fn calculate_chemical_potential( &self, profile: &DFTProfile, @@ -58,12 +58,12 @@ impl DFTSpecifications { /// particles constant in systems, where the number itself is difficult to obtain. pub fn moles_from_profile( profile: &DFTProfile, - ) -> EosResult> + ) -> EosResult> where ::Larger: Dimension, { let rho = profile.density.to_reduced(U::reference_density())?; - Ok(Rc::new(Self::Moles { + Ok(Arc::new(Self::Moles { moles: profile.integrate_reduced_comp(&rho), })) } @@ -74,7 +74,7 @@ impl DFTSpecifications { /// particles constant in systems, e.g. to fix the equimolar dividing surface. pub fn total_moles_from_profile( profile: &DFTProfile, - ) -> EosResult> + ) -> EosResult> where ::Larger: Dimension, { @@ -82,7 +82,7 @@ impl DFTSpecifications { .density .to_reduced(U::reference_density())? .sum_axis(Axis_nd(0)); - Ok(Rc::new(Self::TotalMoles { + Ok(Arc::new(Self::TotalMoles { total_moles: profile.integrate_reduced(rho), chemical_potential: profile.reduced_chemical_potential()?, })) @@ -117,11 +117,11 @@ impl DFTSpecification { pub grid: Grid, - pub convolver: Rc>, - pub dft: Rc>, + pub convolver: Arc>, + pub dft: Arc>, pub temperature: QuantityScalar, pub density: QuantityArray, - pub specification: Rc>, + pub specification: Arc>, pub external_potential: Array, pub bulk: State>, } @@ -187,7 +187,7 @@ where /// after this call if something else is required. pub fn new( grid: Grid, - convolver: Rc>, + convolver: Arc>, bulk: &State>, external_potential: Option>, density: Option<&QuantityArray>, @@ -228,7 +228,7 @@ where dft: bulk.eos.clone(), temperature: bulk.temperature, density, - specification: Rc::new(DFTSpecifications::ChemicalPotential), + specification: Arc::new(DFTSpecifications::ChemicalPotential), external_potential, bulk: bulk.clone(), }) diff --git a/feos-dft/src/python/adsorption/mod.rs b/feos-dft/src/python/adsorption/mod.rs index 89a1983e4..bb67130ab 100644 --- a/feos-dft/src/python/adsorption/mod.rs +++ b/feos-dft/src/python/adsorption/mod.rs @@ -7,11 +7,11 @@ pub use external_potential::PyExternalPotential; macro_rules! impl_adsorption { ($func:ty, $py_func:ty) => { /// Container structure for adsorption isotherms in 1D pores. - #[pyclass(name = "Adsorption1D", unsendable)] + #[pyclass(name = "Adsorption1D")] pub struct PyAdsorption1D(Adsorption1D); /// Container structure for adsorption isotherms in 3D pores. - #[pyclass(name = "Adsorption3D", unsendable)] + #[pyclass(name = "Adsorption3D")] pub struct PyAdsorption3D(Adsorption3D); impl_adsorption_isotherm!($func, $py_func, PyAdsorption1D, PyPore1D, PyPoreProfile1D); diff --git a/feos-dft/src/python/adsorption/pore.rs b/feos-dft/src/python/adsorption/pore.rs index 1da59fcaa..72aa616e4 100644 --- a/feos-dft/src/python/adsorption/pore.rs +++ b/feos-dft/src/python/adsorption/pore.rs @@ -25,7 +25,7 @@ macro_rules! impl_pore { #[pyo3(text_signature = "(geometry, pore_size, potential, n_grid=None, potential_cutoff=None)")] pub struct PyPore1D(Pore1D); - #[pyclass(name = "PoreProfile1D", unsendable)] + #[pyclass(name = "PoreProfile1D")] pub struct PyPoreProfile1D(PoreProfile1D); impl_1d_profile!(PyPoreProfile1D, [get_r, get_z]); @@ -151,7 +151,7 @@ macro_rules! impl_pore { #[pyo3(text_signature = "(system_size, n_grid, coordinates, sigma_ss, epsilon_k_ss, potential_cutoff=None, cutoff_radius=None)")] pub struct PyPore3D(Pore3D); - #[pyclass(name = "PoreProfile3D", unsendable)] + #[pyclass(name = "PoreProfile3D")] pub struct PyPoreProfile3D(PoreProfile3D); impl_3d_profile!(PyPoreProfile3D, get_x, get_y, get_z); diff --git a/feos-dft/src/python/interface/mod.rs b/feos-dft/src/python/interface/mod.rs index c8f23fd93..f0787b4ab 100644 --- a/feos-dft/src/python/interface/mod.rs +++ b/feos-dft/src/python/interface/mod.rs @@ -4,7 +4,7 @@ mod surface_tension_diagram; macro_rules! impl_planar_interface { ($func:ty) => { /// A one-dimensional density profile of a vapor-liquid or liquid-liquid interface. - #[pyclass(name = "PlanarInterface", unsendable)] + #[pyclass(name = "PlanarInterface")] pub struct PyPlanarInterface(PlanarInterface); impl_1d_profile!(PyPlanarInterface, [get_z]); diff --git a/feos-dft/src/python/interface/surface_tension_diagram.rs b/feos-dft/src/python/interface/surface_tension_diagram.rs index f45b9c52f..0f2687b8a 100644 --- a/feos-dft/src/python/interface/surface_tension_diagram.rs +++ b/feos-dft/src/python/interface/surface_tension_diagram.rs @@ -26,7 +26,7 @@ macro_rules! impl_surface_tension_diagram { /// ------- /// SurfaceTensionDiagram /// - #[pyclass(name = "SurfaceTensionDiagram", unsendable)] + #[pyclass(name = "SurfaceTensionDiagram")] #[pyo3(text_signature = "(dia, init_densities=None, n_grid=None, l_grid=None, critical_temperature=None, solver=None)")] pub struct PySurfaceTensionDiagram(SurfaceTensionDiagram); diff --git a/feos-dft/src/python/solvation.rs b/feos-dft/src/python/solvation.rs index 42ff1d633..6bd6d1959 100644 --- a/feos-dft/src/python/solvation.rs +++ b/feos-dft/src/python/solvation.rs @@ -26,7 +26,7 @@ macro_rules! impl_solvation_profile { /// ------- /// SolvationProfile /// - #[pyclass(name = "SolvationProfile", unsendable)] + #[pyclass(name = "SolvationProfile")] #[pyo3(text_signature = "(bulk, n_grid, coordinates, sigma, epsilon_k, system_size=None, cutoff_radius=None, potential_cutoff=None)")] pub struct PySolvationProfile(SolvationProfile); @@ -92,7 +92,7 @@ macro_rules! impl_pair_correlation { /// ------- /// PairCorrelation /// - #[pyclass(name = "PairCorrelation", unsendable)] + #[pyclass(name = "PairCorrelation")] #[pyo3(text_signature = "(bulk, test_particle, n_grid, width)")] pub struct PyPairCorrelation(PairCorrelation); diff --git a/src/association/mod.rs b/src/association/mod.rs index 3c6e2c39c..eb0e304f3 100644 --- a/src/association/mod.rs +++ b/src/association/mod.rs @@ -7,7 +7,7 @@ use num_dual::linalg::{norm, LU}; use num_dual::*; use serde::{Deserialize, Serialize}; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; #[cfg(feature = "dft")] mod dft; @@ -114,7 +114,7 @@ impl AssociationParameters { /// Implementation of the SAFT association Helmholtz energy /// contribution and functional. pub struct Association

{ - parameters: Rc

, + parameters: Arc

, association_parameters: AssociationParameters, max_iter: usize, tol: f64, @@ -123,7 +123,7 @@ pub struct Association

{ impl Association

{ pub fn new( - parameters: &Rc

, + parameters: &Arc

, association_parameters: &AssociationParameters, max_iter: usize, tol: f64, @@ -138,7 +138,7 @@ impl Association

{ } pub fn new_cross_association( - parameters: &Rc

, + parameters: &Arc

, association_parameters: &AssociationParameters, max_iter: usize, tol: f64, @@ -387,7 +387,7 @@ mod tests_pcsaft { #[test] fn helmholtz_energy() { - let params = Rc::new(water_parameters()); + let params = Arc::new(water_parameters()); let assoc = Association::new(¶ms, ¶ms.association, 50, 1e-10); let t = 350.0; let v = 41.248289328513216; @@ -399,7 +399,7 @@ mod tests_pcsaft { #[test] fn helmholtz_energy_cross() { - let params = Rc::new(water_parameters()); + let params = Arc::new(water_parameters()); let assoc = Association::new_cross_association(¶ms, ¶ms.association, 50, 1e-10); let t = 350.0; let v = 41.248289328513216; @@ -423,7 +423,7 @@ mod tests_gc_pcsaft { #[test] fn test_assoc_propanol() { - let params = Rc::new(propanol()); + let params = Arc::new(propanol()); let contrib = Association::new(¶ms, ¶ms.association, 50, 1e-10); let temperature = 300.0; let volume = METER @@ -443,7 +443,7 @@ mod tests_gc_pcsaft { #[test] fn test_cross_assoc_propanol() { - let params = Rc::new(propanol()); + let params = Arc::new(propanol()); let contrib = Association::new_cross_association(¶ms, ¶ms.association, 50, 1e-10); let temperature = 300.0; let volume = METER @@ -463,7 +463,7 @@ mod tests_gc_pcsaft { #[test] fn test_cross_assoc_ethanol_propanol() { - let params = Rc::new(ethanol_propanol(false)); + let params = Arc::new(ethanol_propanol(false)); let contrib = Association::new(¶ms, ¶ms.association, 50, 1e-10); let temperature = 300.0; let volume = METER diff --git a/src/estimator/binary_vle.rs b/src/estimator/binary_vle.rs index fe50e8e31..1053533df 100644 --- a/src/estimator/binary_vle.rs +++ b/src/estimator/binary_vle.rs @@ -7,7 +7,7 @@ use ndarray::{arr1, s, Array1, ArrayView1, Axis}; use quantity::{Quantity, QuantityArray1, QuantityScalar}; use std::collections::HashMap; use std::fmt::LowerExp; -use std::rc::Rc; +use std::sync::Arc; /// Different phases of experimental data points in the `BinaryVlePressure` data set. #[derive(Clone, Copy)] @@ -66,7 +66,7 @@ where ] } - fn predict(&self, eos: &Rc) -> Result, EstimatorError> + fn predict(&self, eos: &Arc) -> Result, EstimatorError> where QuantityScalar: std::fmt::Display + std::fmt::LowerExp, { @@ -155,7 +155,7 @@ where vec } - fn predict(&self, eos: &Rc) -> Result, EstimatorError> + fn predict(&self, eos: &Arc) -> Result, EstimatorError> where QuantityScalar: std::fmt::Display + std::fmt::LowerExp, { @@ -265,7 +265,7 @@ where vec } - fn predict(&self, eos: &Rc) -> Result, EstimatorError> + fn predict(&self, eos: &Arc) -> Result, EstimatorError> where QuantityScalar: std::fmt::Display + std::fmt::LowerExp, { diff --git a/src/estimator/dataset.rs b/src/estimator/dataset.rs index 49d7880ff..657b5eee5 100644 --- a/src/estimator/dataset.rs +++ b/src/estimator/dataset.rs @@ -9,13 +9,13 @@ use ndarray::Array1; use quantity::{QuantityArray1, QuantityScalar}; use std::collections::HashMap; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; /// Utilities for working with experimental data. /// /// Functionalities in the context of optimizations of /// parameters of equations of state. -pub trait DataSet { +pub trait DataSet: Send + Sync { /// Return target quantity. fn target(&self) -> &QuantityArray1; @@ -26,12 +26,12 @@ pub trait DataSet { fn input_str(&self) -> Vec<&str>; /// Evaluation of the equation of state for the target quantity. - fn predict(&self, eos: &Rc) -> Result, EstimatorError> + fn predict(&self, eos: &Arc) -> Result, EstimatorError> where QuantityScalar: std::fmt::Display + std::fmt::LowerExp; /// Evaluate the cost function. - fn cost(&self, eos: &Rc, loss: Loss) -> Result, EstimatorError> + fn cost(&self, eos: &Arc, loss: Loss) -> Result, EstimatorError> where QuantityScalar: std::fmt::Display + std::fmt::LowerExp, { @@ -50,7 +50,7 @@ pub trait DataSet { } /// Returns the relative difference between the equation of state and the experimental values. - fn relative_difference(&self, eos: &Rc) -> Result, EstimatorError> + fn relative_difference(&self, eos: &Arc) -> Result, EstimatorError> where QuantityScalar: std::fmt::Display + std::fmt::LowerExp, { @@ -60,7 +60,7 @@ pub trait DataSet { } /// Returns the mean of the absolute relative difference between the equation of state and the experimental values. - fn mean_absolute_relative_difference(&self, eos: &Rc) -> Result + fn mean_absolute_relative_difference(&self, eos: &Arc) -> Result where QuantityScalar: std::fmt::Display + std::fmt::LowerExp, { diff --git a/src/estimator/diffusion.rs b/src/estimator/diffusion.rs index e8d1bc3a1..a86dff3ce 100644 --- a/src/estimator/diffusion.rs +++ b/src/estimator/diffusion.rs @@ -3,7 +3,7 @@ use feos_core::{DensityInitialization, EntropyScaling, EosUnit, EquationOfState, use ndarray::{arr1, Array1}; use quantity::{QuantityArray1, QuantityScalar}; use std::collections::HashMap; -use std::rc::Rc; +use std::sync::Arc; /// Store experimental diffusion data. #[derive(Clone)] @@ -11,7 +11,6 @@ pub struct Diffusion { pub target: QuantityArray1, temperature: QuantityArray1, pressure: QuantityArray1, - datapoints: usize, } impl Diffusion { @@ -21,12 +20,10 @@ impl Diffusion { temperature: QuantityArray1, pressure: QuantityArray1, ) -> Result { - let datapoints = target.len(); Ok(Self { target, temperature, pressure, - datapoints, }) } @@ -54,20 +51,34 @@ impl> DataSet for Diffu vec!["temperature", "pressure"] } - fn predict(&self, eos: &Rc) -> Result, EstimatorError> + fn predict(&self, eos: &Arc) -> Result, EstimatorError> where QuantityScalar: std::fmt::Display + std::fmt::LowerExp, { - let unit = self.target.get(0); - let mut prediction = Array1::zeros(self.datapoints) * unit; let moles = arr1(&[1.0]) * U::reference_moles(); - for i in 0..self.datapoints { - let t = self.temperature.get(i); - let p = self.pressure.get(i); - let state = State::new_npt(eos, t, p, &moles, DensityInitialization::None)?; - prediction.try_set(i, state.diffusion()?)?; - } - Ok(prediction) + let ts = self + .temperature + .to_reduced(U::reference_temperature()) + .unwrap(); + let ps = self.pressure.to_reduced(U::reference_pressure()).unwrap(); + + let res = ts + .iter() + .zip(ps.iter()) + .map(|(&t, &p)| { + State::new_npt( + eos, + t * U::reference_temperature(), + p * U::reference_pressure(), + &moles, + DensityInitialization::None, + )? + .diffusion()? + .to_reduced(U::reference_diffusion()) + .map_err(EstimatorError::from) + }) + .collect::, EstimatorError>>(); + Ok(Array1::from_vec(res?) * U::reference_diffusion()) } fn get_input(&self) -> HashMap> { diff --git a/src/estimator/estimator.rs b/src/estimator/estimator.rs index fd7aaa972..6e7cb4212 100644 --- a/src/estimator/estimator.rs +++ b/src/estimator/estimator.rs @@ -9,12 +9,12 @@ use quantity::QuantityScalar; use std::fmt; use std::fmt::Display; use std::fmt::Write; -use std::rc::Rc; +use std::sync::Arc; /// A collection of [`DataSet`]s and weights that can be used to /// evaluate an equation of state versus experimental data. pub struct Estimator { - data: Vec>>, + data: Vec>>, weights: Vec, losses: Vec, } @@ -27,7 +27,7 @@ where /// /// The weights are normalized and used as multiplicator when the /// cost function across all `DataSet`s is evaluated. - pub fn new(data: Vec>>, weights: Vec, losses: Vec) -> Self { + pub fn new(data: Vec>>, weights: Vec, losses: Vec) -> Self { Self { data, weights, @@ -36,7 +36,7 @@ where } /// Add a `DataSet` and its weight. - pub fn add_data(&mut self, data: &Rc>, weight: f64, loss: Loss) { + pub fn add_data(&mut self, data: &Arc>, weight: f64, loss: Loss) { self.data.push(data.clone()); self.weights.push(weight); self.losses.push(loss); @@ -45,7 +45,7 @@ where /// Returns the cost of each `DataSet`. /// /// Each cost contains the inverse weight. - pub fn cost(&self, eos: &Rc) -> Result, EstimatorError> { + pub fn cost(&self, eos: &Arc) -> Result, EstimatorError> { let w = arr1(&self.weights) / self.weights.iter().sum::(); let predictions = self .data @@ -58,12 +58,12 @@ where } /// Returns the properties as computed by the equation of state for each `DataSet`. - pub fn predict(&self, eos: &Rc) -> Result>, EstimatorError> { + pub fn predict(&self, eos: &Arc) -> Result>, EstimatorError> { self.data.iter().map(|d| d.predict(eos)).collect() } /// Returns the relative difference for each `DataSet`. - pub fn relative_difference(&self, eos: &Rc) -> Result>, EstimatorError> { + pub fn relative_difference(&self, eos: &Arc) -> Result>, EstimatorError> { self.data .iter() .map(|d| d.relative_difference(eos)) @@ -73,7 +73,7 @@ where /// Returns the mean absolute relative difference for each `DataSet`. pub fn mean_absolute_relative_difference( &self, - eos: &Rc, + eos: &Arc, ) -> Result, EstimatorError> { self.data .iter() @@ -82,7 +82,7 @@ where } /// Returns the stored `DataSet`s. - pub fn datasets(&self) -> Vec>> { + pub fn datasets(&self) -> Vec>> { self.data.to_vec() } diff --git a/src/estimator/liquid_density.rs b/src/estimator/liquid_density.rs index 4a0446df0..c82066c1c 100644 --- a/src/estimator/liquid_density.rs +++ b/src/estimator/liquid_density.rs @@ -3,10 +3,10 @@ use feos_core::{ DensityInitialization, EosUnit, EquationOfState, MolarWeight, PhaseEquilibrium, SolverOptions, State, }; -use ndarray::{arr1, Array1}; +use ndarray::arr1; use quantity::{QuantityArray1, QuantityScalar}; use std::collections::HashMap; -use std::rc::Rc; +use std::sync::Arc; /// Liquid mass density data as function of pressure and temperature. #[derive(Clone)] @@ -17,8 +17,6 @@ pub struct LiquidDensity { temperature: QuantityArray1, /// pressure pressure: QuantityArray1, - /// number of data points - datapoints: usize, } impl LiquidDensity { @@ -28,12 +26,10 @@ impl LiquidDensity { temperature: QuantityArray1, pressure: QuantityArray1, ) -> Result { - let datapoints = target.len(); Ok(Self { target, temperature, pressure, - datapoints, }) } @@ -61,25 +57,21 @@ impl> DataSet for LiquidDe vec!["temperature", "pressure"] } - fn predict(&self, eos: &Rc) -> Result, EstimatorError> { + fn predict(&self, eos: &Arc) -> Result, EstimatorError> { let moles = arr1(&[1.0]) * U::reference_moles(); - let unit = self.target.get(0); - let mut prediction = Array1::zeros(self.datapoints) * unit; - for i in 0..self.datapoints { - let state = State::new_npt( - eos, - self.temperature.get(i), - self.pressure.get(i), - &moles, - DensityInitialization::Liquid, - ); - if let Ok(s) = state { - prediction.try_set(i, s.mass_density())?; - } else { - prediction.try_set(i, f64::NAN * unit)?; - } - } - Ok(prediction) + Ok(self + .temperature + .into_iter() + .zip(self.pressure.into_iter()) + .map(|(t, p)| { + let state = State::new_npt(eos, t, p, &moles, DensityInitialization::Liquid); + if let Ok(s) = state { + s.mass_density() + } else { + f64::NAN * U::reference_mass() / U::reference_volume() + } + }) + .collect()) } fn get_input(&self) -> HashMap> { @@ -95,7 +87,6 @@ impl> DataSet for LiquidDe pub struct EquilibriumLiquidDensity { pub target: QuantityArray1, temperature: QuantityArray1, - datapoints: usize, solver_options: SolverOptions, } @@ -106,11 +97,9 @@ impl EquilibriumLiquidDensity { temperature: QuantityArray1, vle_options: Option, ) -> Result { - let datapoints = target.len(); Ok(Self { target, temperature, - datapoints, solver_options: vle_options.unwrap_or_default(), }) } @@ -136,22 +125,21 @@ impl> DataSet vec!["temperature"] } - fn predict(&self, eos: &Rc) -> Result, EstimatorError> + fn predict(&self, eos: &Arc) -> Result, EstimatorError> where QuantityScalar: std::fmt::Display + std::fmt::LowerExp, { - 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); - if let Ok(state) = PhaseEquilibrium::pure(eos, t, None, self.solver_options) { - prediction.try_set(i, state.liquid().mass_density())?; - } else { - prediction.try_set(i, f64::NAN * U::reference_mass() / U::reference_volume())? - } - } - Ok(prediction) + Ok(self + .temperature + .into_iter() + .map(|t| { + if let Ok(state) = PhaseEquilibrium::pure(eos, t, None, self.solver_options) { + state.liquid().mass_density() + } else { + f64::NAN * U::reference_mass() / U::reference_volume() + } + }) + .collect()) } fn get_input(&self) -> HashMap> { diff --git a/src/estimator/mod.rs b/src/estimator/mod.rs index 0fe0d63b1..4f5894df9 100644 --- a/src/estimator/mod.rs +++ b/src/estimator/mod.rs @@ -6,16 +6,18 @@ use thiserror::Error; mod dataset; pub use dataset::DataSet; -mod binary_vle; -pub use binary_vle::{BinaryPhaseDiagram, BinaryVleChemicalPotential, BinaryVlePressure, Phase}; mod estimator; pub use estimator::Estimator; mod loss; pub use loss::Loss; + +// Properties mod vapor_pressure; pub use vapor_pressure::VaporPressure; mod liquid_density; pub use liquid_density::{EquilibriumLiquidDensity, LiquidDensity}; +mod binary_vle; +pub use binary_vle::{BinaryPhaseDiagram, BinaryVleChemicalPotential, BinaryVlePressure, Phase}; mod viscosity; pub use viscosity::Viscosity; mod thermal_conductivity; diff --git a/src/estimator/python.rs b/src/estimator/python.rs index 5d4314a38..4a88e69ef 100644 --- a/src/estimator/python.rs +++ b/src/estimator/python.rs @@ -11,7 +11,7 @@ impl From for PyErr { #[macro_export] macro_rules! impl_estimator { ($eos:ty, $py_eos:ty) => { - #[pyclass(name = "Loss", unsendable)] + #[pyclass(name = "Loss")] #[derive(Clone)] pub struct PyLoss(Loss); @@ -113,9 +113,9 @@ macro_rules! impl_estimator { /// A collection of experimental data that can be used to compute /// cost functions and make predictions using an equation of state. - #[pyclass(name = "DataSet", unsendable)] + #[pyclass(name = "DataSet")] #[derive(Clone)] - pub struct PyDataSet(Rc>); + pub struct PyDataSet(Arc>); #[pymethods] impl PyDataSet { @@ -244,7 +244,7 @@ macro_rules! impl_estimator { tol: Option, verbosity: Option, ) -> PyResult { - Ok(Self(Rc::new(VaporPressure::::new( + Ok(Self(Arc::new(VaporPressure::::new( target.clone().into(), temperature.clone().into(), extrapolate.unwrap_or(false), @@ -273,7 +273,7 @@ macro_rules! impl_estimator { temperature: &PySIArray1, pressure: &PySIArray1, ) -> PyResult { - Ok(Self(Rc::new(LiquidDensity::::new( + Ok(Self(Arc::new(LiquidDensity::::new( target.clone().into(), temperature.clone().into(), pressure.clone().into(), @@ -311,7 +311,7 @@ macro_rules! impl_estimator { tol: Option, verbosity: Option, ) -> PyResult { - Ok(Self(Rc::new(EquilibriumLiquidDensity::::new( + Ok(Self(Arc::new(EquilibriumLiquidDensity::::new( target.clone().into(), temperature.clone().into(), Some((max_iter, tol, verbosity).into()), @@ -343,7 +343,7 @@ macro_rules! impl_estimator { liquid_molefracs: &PyArray1, vapor_molefracs: &PyArray1, ) -> Self { - Self(Rc::new(BinaryVleChemicalPotential::new( + Self(Arc::new(BinaryVleChemicalPotential::new( temperature.clone().into(), pressure.clone().into(), liquid_molefracs.to_owned_array(), @@ -376,7 +376,7 @@ macro_rules! impl_estimator { molefracs: &PyArray1, phase: Phase, ) -> Self { - Self(Rc::new(BinaryVlePressure::new( + Self(Arc::new(BinaryVlePressure::new( temperature.clone().into(), pressure.clone().into(), molefracs.to_owned_array(), @@ -414,7 +414,7 @@ macro_rules! impl_estimator { vapor_molefracs: Option<&PyArray1>, npoints: Option, ) -> Self { - Self(Rc::new(BinaryPhaseDiagram::new( + Self(Arc::new(BinaryPhaseDiagram::new( specification.into(), temperature_or_pressure.clone().into(), liquid_molefracs.map(|x| x.to_owned_array()), @@ -466,7 +466,7 @@ macro_rules! impl_estimator { /// Returns /// ------- /// Estimator - #[pyclass(name = "Estimator", unsendable)] + #[pyclass(name = "Estimator")] #[pyo3(text_signature = "(data, weights, losses)")] pub struct PyEstimator(Estimator); @@ -633,7 +633,7 @@ macro_rules! impl_estimator_entropy_scaling { temperature: &PySIArray1, pressure: &PySIArray1, ) -> PyResult { - Ok(Self(Rc::new(Viscosity::::new( + Ok(Self(Arc::new(Viscosity::::new( target.clone().into(), temperature.clone().into(), pressure.clone().into(), @@ -661,7 +661,7 @@ macro_rules! impl_estimator_entropy_scaling { temperature: &PySIArray1, pressure: &PySIArray1, ) -> PyResult { - Ok(Self(Rc::new(ThermalConductivity::::new( + Ok(Self(Arc::new(ThermalConductivity::::new( target.clone().into(), temperature.clone().into(), pressure.clone().into(), @@ -689,7 +689,7 @@ macro_rules! impl_estimator_entropy_scaling { temperature: &PySIArray1, pressure: &PySIArray1, ) -> PyResult { - Ok(Self(Rc::new(Diffusion::::new( + Ok(Self(Arc::new(Diffusion::::new( target.clone().into(), temperature.clone().into(), pressure.clone().into(), diff --git a/src/estimator/thermal_conductivity.rs b/src/estimator/thermal_conductivity.rs index 281b35a45..3d64465c9 100644 --- a/src/estimator/thermal_conductivity.rs +++ b/src/estimator/thermal_conductivity.rs @@ -3,7 +3,7 @@ use feos_core::{DensityInitialization, EntropyScaling, EosUnit, EquationOfState, use ndarray::{arr1, Array1}; use quantity::{QuantityArray1, QuantityScalar}; use std::collections::HashMap; -use std::rc::Rc; +use std::sync::Arc; /// Store experimental thermal conductivity data. #[derive(Clone)] @@ -11,7 +11,6 @@ pub struct ThermalConductivity { pub target: QuantityArray1, temperature: QuantityArray1, pressure: QuantityArray1, - datapoints: usize, } impl ThermalConductivity { @@ -21,12 +20,10 @@ impl ThermalConductivity { temperature: QuantityArray1, pressure: QuantityArray1, ) -> Result { - let datapoints = target.len(); Ok(Self { target, temperature, pressure, - datapoints, }) } @@ -54,20 +51,38 @@ impl> DataSet for Therm vec!["temperature", "pressure"] } - fn predict(&self, eos: &Rc) -> Result, EstimatorError> + fn predict(&self, eos: &Arc) -> Result, EstimatorError> where QuantityScalar: std::fmt::Display + std::fmt::LowerExp, { - let unit = self.target.get(0); - let mut prediction = Array1::zeros(self.datapoints) * unit; let moles = arr1(&[1.0]) * U::reference_moles(); - for i in 0..self.datapoints { - let t = self.temperature.get(i); - let p = self.pressure.get(i); - let state = State::new_npt(eos, t, p, &moles, DensityInitialization::None)?; - prediction.try_set(i, state.thermal_conductivity()?)?; - } - Ok(prediction) + let ts = self + .temperature + .to_reduced(U::reference_temperature()) + .unwrap(); + let ps = self.pressure.to_reduced(U::reference_pressure()).unwrap(); + let unit = U::reference_energy() + / U::reference_time() + / U::reference_temperature() + / U::reference_length(); + + let res = ts + .iter() + .zip(ps.iter()) + .map(|(&t, &p)| { + State::new_npt( + eos, + t * U::reference_temperature(), + p * U::reference_pressure(), + &moles, + DensityInitialization::None, + )? + .thermal_conductivity()? + .to_reduced(unit) + .map_err(EstimatorError::from) + }) + .collect::, EstimatorError>>(); + Ok(Array1::from_vec(res?) * unit) } fn get_input(&self) -> HashMap> { diff --git a/src/estimator/vapor_pressure.rs b/src/estimator/vapor_pressure.rs index cd9021725..69707a82f 100644 --- a/src/estimator/vapor_pressure.rs +++ b/src/estimator/vapor_pressure.rs @@ -3,7 +3,7 @@ use feos_core::{Contributions, EosUnit, EquationOfState, PhaseEquilibrium, Solve use ndarray::Array1; use quantity::{QuantityArray1, QuantityScalar}; use std::collections::HashMap; -use std::rc::Rc; +use std::sync::Arc; /// Store experimental vapor pressure data. #[derive(Clone)] @@ -67,7 +67,7 @@ impl DataSet for VaporPressure { vec!["temperature"] } - fn predict(&self, eos: &Rc) -> Result, EstimatorError> + fn predict(&self, eos: &Arc) -> Result, EstimatorError> where QuantityScalar: std::fmt::Display + std::fmt::LowerExp, { diff --git a/src/estimator/viscosity.rs b/src/estimator/viscosity.rs index 48b74f1d7..2677c36ef 100644 --- a/src/estimator/viscosity.rs +++ b/src/estimator/viscosity.rs @@ -1,9 +1,9 @@ use super::{DataSet, EstimatorError}; use feos_core::{DensityInitialization, EntropyScaling, EosUnit, EquationOfState, State}; -use ndarray::{arr1, Array1}; +use ndarray::arr1; use quantity::{QuantityArray1, QuantityScalar}; use std::collections::HashMap; -use std::rc::Rc; +use std::sync::Arc; /// Store experimental viscosity data. #[derive(Clone)] @@ -11,7 +11,6 @@ pub struct Viscosity { pub target: QuantityArray1, temperature: QuantityArray1, pressure: QuantityArray1, - datapoints: usize, } impl Viscosity { @@ -21,12 +20,10 @@ impl Viscosity { temperature: QuantityArray1, pressure: QuantityArray1, ) -> Result { - let datapoints = target.len(); Ok(Self { target, temperature, pressure, - datapoints, }) } @@ -54,20 +51,20 @@ impl> DataSet for Visco vec!["temperature", "pressure"] } - fn predict(&self, eos: &Rc) -> Result, EstimatorError> + fn predict(&self, eos: &Arc) -> Result, EstimatorError> where QuantityScalar: std::fmt::Display + std::fmt::LowerExp, { - let unit = self.target.get(0); - let mut prediction = Array1::zeros(self.datapoints) * unit; let moles = arr1(&[1.0]) * U::reference_moles(); - for i in 0..self.datapoints { - let t = self.temperature.get(i); - let p = self.pressure.get(i); - let state = State::new_npt(eos, t, p, &moles, DensityInitialization::None)?; - prediction.try_set(i, state.viscosity()?)?; - } - Ok(prediction) + self.temperature + .into_iter() + .zip(self.pressure.into_iter()) + .map(|(t, p)| { + State::new_npt(eos, t, p, &moles, DensityInitialization::None)? + .viscosity() + .map_err(EstimatorError::from) + }) + .collect() } fn get_input(&self) -> HashMap> { diff --git a/src/gc_pcsaft/dft/dispersion.rs b/src/gc_pcsaft/dft/dispersion.rs index c1ed13872..4eade834f 100644 --- a/src/gc_pcsaft/dft/dispersion.rs +++ b/src/gc_pcsaft/dft/dispersion.rs @@ -9,15 +9,15 @@ use ndarray::*; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_6, PI}; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; #[derive(Clone)] pub struct AttractiveFunctional { - parameters: Rc, + parameters: Arc, } impl AttractiveFunctional { - pub fn new(parameters: &Rc) -> Self { + pub fn new(parameters: &Arc) -> Self { Self { parameters: parameters.clone(), } diff --git a/src/gc_pcsaft/dft/hard_chain.rs b/src/gc_pcsaft/dft/hard_chain.rs index c2fdc93b3..b7dc74793 100644 --- a/src/gc_pcsaft/dft/hard_chain.rs +++ b/src/gc_pcsaft/dft/hard_chain.rs @@ -8,15 +8,15 @@ use ndarray::*; use num_dual::DualNum; use petgraph::visit::EdgeRef; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; #[derive(Clone)] pub struct ChainFunctional { - parameters: Rc, + parameters: Arc, } impl ChainFunctional { - pub fn new(parameters: &Rc) -> Self { + pub fn new(parameters: &Arc) -> Self { Self { parameters: parameters.clone(), } diff --git a/src/gc_pcsaft/dft/mod.rs b/src/gc_pcsaft/dft/mod.rs index 5fff0d02e..c201c90c7 100644 --- a/src/gc_pcsaft/dft/mod.rs +++ b/src/gc_pcsaft/dft/mod.rs @@ -10,7 +10,7 @@ use num_dual::DualNum; use petgraph::graph::UnGraph; use quantity::si::{SIArray1, SIUnit, GRAM, MOL}; use std::f64::consts::FRAC_PI_6; -use std::rc::Rc; +use std::sync::Arc; mod dispersion; mod hard_chain; @@ -21,14 +21,14 @@ pub use parameter::GcPcSaftFunctionalParameters; /// gc-PC-SAFT Helmholtz energy functional. pub struct GcPcSaftFunctional { - pub parameters: Rc, + pub parameters: Arc, fmt_version: FMTVersion, options: GcPcSaftOptions, contributions: Vec>, } impl GcPcSaftFunctional { - pub fn new(parameters: Rc) -> DFT { + pub fn new(parameters: Arc) -> DFT { Self::with_options( parameters, FMTVersion::WhiteBear, @@ -37,7 +37,7 @@ impl GcPcSaftFunctional { } pub fn with_options( - parameters: Rc, + parameters: Arc, fmt_version: FMTVersion, saft_options: GcPcSaftOptions, ) -> DFT { @@ -83,7 +83,7 @@ impl HelmholtzEnergyFunctional for GcPcSaftFunctional { fn subset(&self, component_list: &[usize]) -> DFT { Self::with_options( - Rc::new(self.parameters.subset(component_list)), + Arc::new(self.parameters.subset(component_list)), self.fmt_version, self.options, ) diff --git a/src/gc_pcsaft/eos/dispersion.rs b/src/gc_pcsaft/eos/dispersion.rs index a0f9e7d05..a104b35c7 100644 --- a/src/gc_pcsaft/eos/dispersion.rs +++ b/src/gc_pcsaft/eos/dispersion.rs @@ -4,7 +4,7 @@ use feos_core::{HelmholtzEnergyDual, StateHD}; use num_dual::DualNum; use std::f64::consts::PI; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; pub const A0: [f64; 7] = [ 0.91056314451539, @@ -63,7 +63,7 @@ pub const B2: [f64; 7] = [ #[derive(Clone)] pub struct Dispersion { - pub parameters: Rc, + pub parameters: Arc, } impl> HelmholtzEnergyDual for Dispersion { @@ -141,7 +141,7 @@ mod test { fn test_dispersion_propane() { let parameters = propane(); let contrib = Dispersion { - parameters: Rc::new(parameters), + parameters: Arc::new(parameters), }; let temperature = 300.0; let volume = METER @@ -163,7 +163,7 @@ mod test { fn test_dispersion_propanol() { let parameters = propanol(); let contrib = Dispersion { - parameters: Rc::new(parameters), + parameters: Arc::new(parameters), }; let temperature = 300.0; let volume = METER diff --git a/src/gc_pcsaft/eos/hard_chain.rs b/src/gc_pcsaft/eos/hard_chain.rs index 3ce6a03b9..419ac92bc 100644 --- a/src/gc_pcsaft/eos/hard_chain.rs +++ b/src/gc_pcsaft/eos/hard_chain.rs @@ -3,11 +3,11 @@ use crate::hard_sphere::HardSphereProperties; use feos_core::{HelmholtzEnergyDual, StateHD}; use num_dual::*; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; #[derive(Clone)] pub struct HardChain { - pub parameters: Rc, + pub parameters: Arc, } impl> HelmholtzEnergyDual for HardChain { @@ -56,7 +56,7 @@ mod test { fn test_hc_propane() { let parameters = propane(); let contrib = HardChain { - parameters: Rc::new(parameters), + parameters: Arc::new(parameters), }; let temperature = 300.0; let volume = METER @@ -82,7 +82,7 @@ mod test { fn test_hc_propanol() { let parameters = propanol(); let contrib = HardChain { - parameters: Rc::new(parameters), + parameters: Arc::new(parameters), }; let temperature = 300.0; let volume = METER diff --git a/src/gc_pcsaft/eos/mod.rs b/src/gc_pcsaft/eos/mod.rs index aeb407ddc..58d32cd38 100644 --- a/src/gc_pcsaft/eos/mod.rs +++ b/src/gc_pcsaft/eos/mod.rs @@ -6,7 +6,7 @@ use feos_core::{EquationOfState, HelmholtzEnergy, IdealGasContribution, MolarWei use ndarray::Array1; use quantity::si::*; use std::f64::consts::FRAC_PI_6; -use std::rc::Rc; +use std::sync::Arc; pub(crate) mod dispersion; mod hard_chain; @@ -40,18 +40,18 @@ impl Default for GcPcSaftOptions { /// gc-PC-SAFT equation of state pub struct GcPcSaft { - pub parameters: Rc, + pub parameters: Arc, options: GcPcSaftOptions, contributions: Vec>, joback: Joback, } impl GcPcSaft { - pub fn new(parameters: Rc) -> Self { + pub fn new(parameters: Arc) -> Self { Self::with_options(parameters, GcPcSaftOptions::default()) } - pub fn with_options(parameters: Rc, options: GcPcSaftOptions) -> Self { + pub fn with_options(parameters: Arc, options: GcPcSaftOptions) -> Self { let mut contributions: Vec> = Vec::with_capacity(7); contributions.push(Box::new(HardSphere::new(¶meters))); contributions.push(Box::new(HardChain { @@ -90,7 +90,7 @@ impl EquationOfState for GcPcSaft { fn subset(&self, component_list: &[usize]) -> Self { Self::with_options( - Rc::new(self.parameters.subset(component_list)), + Arc::new(self.parameters.subset(component_list)), self.options, ) } @@ -130,7 +130,7 @@ mod test { #[test] fn hs_propane() { let parameters = propane(); - let contrib = HardSphere::new(&Rc::new(parameters)); + let contrib = HardSphere::new(&Arc::new(parameters)); let temperature = 300.0; let volume = METER .powi(3) @@ -150,7 +150,7 @@ mod test { #[test] fn hs_propanol() { let parameters = propanol(); - let contrib = HardSphere::new(&Rc::new(parameters)); + let contrib = HardSphere::new(&Arc::new(parameters)); let temperature = 300.0; let volume = METER .powi(3) @@ -169,7 +169,7 @@ mod test { #[test] fn assoc_propanol() { - let parameters = Rc::new(propanol()); + let parameters = Arc::new(propanol()); let contrib = Association::new(¶meters, ¶meters.association, 50, 1e-10); let temperature = 300.0; let volume = METER @@ -189,7 +189,7 @@ mod test { #[test] fn cross_assoc_propanol() { - let parameters = Rc::new(propanol()); + let parameters = Arc::new(propanol()); let contrib = Association::new_cross_association(¶meters, ¶meters.association, 50, 1e-10); let temperature = 300.0; @@ -210,7 +210,7 @@ mod test { #[test] fn cross_assoc_ethanol_propanol() { - let parameters = Rc::new(ethanol_propanol(false)); + let parameters = Arc::new(ethanol_propanol(false)); let contrib = Association::new(¶meters, ¶meters.association, 50, 1e-10); let temperature = 300.0; let volume = METER diff --git a/src/gc_pcsaft/eos/polar.rs b/src/gc_pcsaft/eos/polar.rs index 301a12a8d..9256524d1 100644 --- a/src/gc_pcsaft/eos/polar.rs +++ b/src/gc_pcsaft/eos/polar.rs @@ -5,7 +5,7 @@ use ndarray::prelude::*; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_3, PI}; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; // Dipole parameters pub const AD: [[f64; 3]; 5] = [ @@ -54,7 +54,7 @@ fn triplet_integral_ijk>(mijk1: f64, mijk2: f64, eta: D) -> D { } pub struct Dipole { - parameters: Rc, + parameters: Arc, mij1: Array2, mij2: Array2, mijk1: Array3, @@ -64,7 +64,7 @@ pub struct Dipole { } impl Dipole { - pub fn new(parameters: &Rc) -> Self { + pub fn new(parameters: &Arc) -> Self { let ndipole = parameters.dipole_comp.len(); let f2_term = Array2::from_shape_fn([ndipole; 2], |(i, j)| { diff --git a/src/gc_pcsaft/micelles.rs b/src/gc_pcsaft/micelles.rs index 2a25fc920..f68af8201 100644 --- a/src/gc_pcsaft/micelles.rs +++ b/src/gc_pcsaft/micelles.rs @@ -5,7 +5,7 @@ use feos_dft::{ }; use ndarray::prelude::*; use quantity::{QuantityArray1, QuantityArray2, QuantityScalar}; -use std::rc::Rc; +use std::sync::Arc; pub enum MicelleInitialization { ExternalPotential(f64, f64), @@ -170,7 +170,7 @@ impl MicelleProfile { )?; // specify specification - profile.specification = Rc::new(specification); + profile.specification = Arc::new(specification); Ok(Self { profile, @@ -211,7 +211,7 @@ impl MicelleProfile { pub fn update_specification(&self, specification: MicelleSpecification) -> Self { let mut profile = self.clone(); - profile.profile.specification = Rc::new(specification); + profile.profile.specification = Arc::new(specification); profile.delta_omega = None; profile.delta_n = None; profile @@ -233,7 +233,7 @@ impl MicelleProfile { let pressure = self.profile.bulk.pressure(Contributions::Total); let eos = self.profile.bulk.eos.clone(); let indices = self.profile.bulk.eos.component_index().into_owned(); - self.profile.specification = Rc::new(MicelleSpecification::ChemicalPotential); + self.profile.specification = Arc::new(MicelleSpecification::ChemicalPotential); for _ in 0..options.max_iter.unwrap_or(MAX_ITER_MICELLE) { // check for convergence diff --git a/src/gc_pcsaft/mod.rs b/src/gc_pcsaft/mod.rs index 863737dac..7ddcd87a3 100644 --- a/src/gc_pcsaft/mod.rs +++ b/src/gc_pcsaft/mod.rs @@ -1,5 +1,5 @@ //! Heterosegmented Group Contribution PC-SAFT -//! +//! //! - [Gross et al. (2003)](https://doi.org/10.1021/ie020509y) //! - [Sauer et al. (2014)](https://doi.org/10.1021/ie502203w) #![warn(clippy::all)] diff --git a/src/gc_pcsaft/python/micelles.rs b/src/gc_pcsaft/python/micelles.rs index 873d45725..4bc658005 100644 --- a/src/gc_pcsaft/python/micelles.rs +++ b/src/gc_pcsaft/python/micelles.rs @@ -2,7 +2,7 @@ macro_rules! impl_micelle_profile { ($func:ty) => { /// A one-dimensional profile of a spherical or cylindrical micelle. - #[pyclass(name = "MicelleProfile", unsendable)] + #[pyclass(name = "MicelleProfile")] pub struct PyMicelleProfile(MicelleProfile); impl_1d_profile!(PyMicelleProfile, [get_r]); diff --git a/src/gc_pcsaft/python/mod.rs b/src/gc_pcsaft/python/mod.rs index 0fa938bb7..a92edd5b0 100644 --- a/src/gc_pcsaft/python/mod.rs +++ b/src/gc_pcsaft/python/mod.rs @@ -13,12 +13,12 @@ use feos_core::{impl_json_handling, impl_parameter_from_segments, impl_segment_r #[cfg(feature = "dft")] use numpy::{PyArray2, ToPyArray}; use pyo3::prelude::*; -use std::rc::Rc; +use std::sync::Arc; #[cfg(feature = "micelles")] mod micelles; -#[pyclass(name = "GcPcSaftRecord", unsendable)] +#[pyclass(name = "GcPcSaftRecord")] #[pyo3(text_signature = "(m, sigma, epsilon_k, mu=None, association_record=None, psi_dft=None)")] #[derive(Clone)] pub struct PyGcPcSaftRecord(GcPcSaftRecord); @@ -83,12 +83,12 @@ impl_segment_record!( PyJobackRecord ); -#[pyclass(name = "GcPcSaftEosParameters", unsendable)] +#[pyclass(name = "GcPcSaftEosParameters")] #[pyo3( text_signature = "(pure_records, segmentbinary_records=None, substances=None, search_option='Name')" )] #[derive(Clone)] -pub struct PyGcPcSaftEosParameters(pub Rc); +pub struct PyGcPcSaftEosParameters(pub Arc); impl_parameter_from_segments!(GcPcSaftEosParameters, PyGcPcSaftEosParameters); @@ -104,12 +104,12 @@ impl PyGcPcSaftEosParameters { } #[cfg(feature = "dft")] -#[pyclass(name = "GcPcSaftFunctionalParameters", unsendable)] +#[pyclass(name = "GcPcSaftFunctionalParameters")] #[pyo3( text_signature = "(pure_records, segmentbinary_records=None, substances=None, search_option='Name')" )] #[derive(Clone)] -pub struct PyGcPcSaftFunctionalParameters(pub Rc); +pub struct PyGcPcSaftFunctionalParameters(pub Arc); #[cfg(feature = "dft")] impl_parameter_from_segments!(GcPcSaftFunctionalParameters, PyGcPcSaftFunctionalParameters); diff --git a/src/hard_sphere/dft.rs b/src/hard_sphere/dft.rs index 50b275f8f..b261b41dd 100644 --- a/src/hard_sphere/dft.rs +++ b/src/hard_sphere/dft.rs @@ -9,7 +9,7 @@ use ndarray::*; use num_dual::DualNum; use std::f64::consts::PI; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; use super::{HardSphereProperties, MonomerShape}; @@ -62,7 +62,7 @@ pub enum FMTVersion { /// /// The geometry coefficients $C_{k,\alpha}$ and the segment diameters $d_\alpha$ are specified via the [HardSphereProperties] trait. pub struct FMTContribution

{ - pub properties: Rc

, + pub properties: Arc

, version: FMTVersion, } @@ -76,7 +76,7 @@ impl

Clone for FMTContribution

{ } impl

FMTContribution

{ - pub fn new(properties: &Rc

, version: FMTVersion) -> Self { + pub fn new(properties: &Arc

, version: FMTVersion) -> Self { Self { properties: properties.clone(), version, @@ -310,14 +310,14 @@ impl HardSphereProperties for HardSphereParameters { /// [HelmholtzEnergyFunctional] for hard sphere systems. pub struct FMTFunctional { - properties: Rc, + properties: Arc, contributions: Vec>, version: FMTVersion, } impl FMTFunctional { pub fn new(sigma: &Array1, version: FMTVersion) -> DFT { - let properties = Rc::new(HardSphereParameters { + let properties = Arc::new(HardSphereParameters { sigma: sigma.clone(), }); let contributions: Vec> = diff --git a/src/hard_sphere/mod.rs b/src/hard_sphere/mod.rs index 4926c24aa..92d96e7a2 100644 --- a/src/hard_sphere/mod.rs +++ b/src/hard_sphere/mod.rs @@ -3,10 +3,9 @@ use feos_core::{HelmholtzEnergyDual, StateHD}; use ndarray::*; use num_dual::DualNum; -use std::borrow::Cow; use std::f64::consts::FRAC_PI_6; use std::fmt; -use std::rc::Rc; +use std::{borrow::Cow, sync::Arc}; #[cfg(feature = "dft")] mod dft; @@ -106,11 +105,11 @@ pub trait HardSphereProperties { /// /// The geometry coefficients $C_{k,\alpha}$ and the segment diameters $d_\alpha$ are specified via the [HardSphereProperties] trait. pub struct HardSphere

{ - parameters: Rc

, + parameters: Arc

, } impl

HardSphere

{ - pub fn new(parameters: &Rc

) -> Self { + pub fn new(parameters: &Arc

) -> Self { Self { parameters: parameters.clone(), } diff --git a/src/lib.rs b/src/lib.rs index 0cb72b6ec..069356ff3 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -8,7 +8,7 @@ //! use feos_core::parameter::{IdentifierOption, Parameter}; //! use feos_core::{Contributions, State}; //! use quantity::si::KELVIN; -//! use std::rc::Rc; +//! use std::sync::Arc; //! //! // Read parameters from json file. //! let parameters = PcSaftParameters::from_json( @@ -19,7 +19,7 @@ //! )?; //! //! // Define equation of state. -//! let saft = Rc::new(PcSaft::new(Rc::new(parameters))); +//! let saft = Arc::new(PcSaft::new(Arc::new(parameters))); //! //! // Define thermodynamic conditions. //! let critical_point = State::critical_point(&saft, None, None, Default::default())?; @@ -33,9 +33,9 @@ #![warn(clippy::all)] #![allow(clippy::too_many_arguments)] -pub mod eos; #[cfg(feature = "dft")] pub mod dft; +pub mod eos; #[cfg(feature = "estimator")] pub mod estimator; diff --git a/src/pcsaft/dft/dispersion.rs b/src/pcsaft/dft/dispersion.rs index f368845e1..05908daac 100644 --- a/src/pcsaft/dft/dispersion.rs +++ b/src/pcsaft/dft/dispersion.rs @@ -10,7 +10,7 @@ use ndarray::*; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_3, PI}; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; /// psi Parameter for DFT (Sauer2017) const PSI_DFT: f64 = 1.3862; @@ -19,11 +19,11 @@ const PSI_PDGT: f64 = 1.3286; #[derive(Clone)] pub struct AttractiveFunctional { - parameters: Rc, + parameters: Arc, } impl AttractiveFunctional { - pub fn new(parameters: Rc) -> Self { + pub fn new(parameters: Arc) -> Self { Self { parameters } } } diff --git a/src/pcsaft/dft/hard_chain.rs b/src/pcsaft/dft/hard_chain.rs index 1a1ea8893..df8d2013a 100644 --- a/src/pcsaft/dft/hard_chain.rs +++ b/src/pcsaft/dft/hard_chain.rs @@ -7,15 +7,15 @@ use feos_dft::{ use ndarray::*; use num_dual::DualNum; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; #[derive(Clone)] pub struct ChainFunctional { - parameters: Rc, + parameters: Arc, } impl ChainFunctional { - pub fn new(parameters: Rc) -> Self { + pub fn new(parameters: Arc) -> Self { Self { parameters } } } diff --git a/src/pcsaft/dft/mod.rs b/src/pcsaft/dft/mod.rs index df6d18d34..593e1dc78 100644 --- a/src/pcsaft/dft/mod.rs +++ b/src/pcsaft/dft/mod.rs @@ -12,7 +12,7 @@ use ndarray::{Array1, Array2}; use num_traits::One; use quantity::si::*; use std::f64::consts::FRAC_PI_6; -use std::rc::Rc; +use std::sync::Arc; mod dispersion; mod hard_chain; @@ -24,7 +24,7 @@ use pure_saft_functional::*; /// PC-SAFT Helmholtz energy functional. pub struct PcSaftFunctional { - pub parameters: Rc, + pub parameters: Arc, fmt_version: FMTVersion, options: PcSaftOptions, contributions: Vec>, @@ -32,16 +32,16 @@ pub struct PcSaftFunctional { } impl PcSaftFunctional { - pub fn new(parameters: Rc) -> DFT { + pub fn new(parameters: Arc) -> DFT { Self::with_options(parameters, FMTVersion::WhiteBear, PcSaftOptions::default()) } - pub fn new_full(parameters: Rc, fmt_version: FMTVersion) -> DFT { + pub fn new_full(parameters: Arc, fmt_version: FMTVersion) -> DFT { Self::with_options(parameters, fmt_version, PcSaftOptions::default()) } pub fn with_options( - parameters: Rc, + parameters: Arc, fmt_version: FMTVersion, saft_options: PcSaftOptions, ) -> DFT { @@ -106,7 +106,7 @@ impl PcSaftFunctional { impl HelmholtzEnergyFunctional for PcSaftFunctional { fn subset(&self, component_list: &[usize]) -> DFT { Self::with_options( - Rc::new(self.parameters.subset(component_list)), + Arc::new(self.parameters.subset(component_list)), self.fmt_version, self.options, ) diff --git a/src/pcsaft/dft/pure_saft_functional.rs b/src/pcsaft/dft/pure_saft_functional.rs index 464eb672a..66279ef3d 100644 --- a/src/pcsaft/dft/pure_saft_functional.rs +++ b/src/pcsaft/dft/pure_saft_functional.rs @@ -12,7 +12,7 @@ use ndarray::*; use num_dual::*; use std::f64::consts::{FRAC_PI_6, PI}; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; const PI36M1: f64 = 1.0 / (36.0 * PI); const N3_CUTOFF: f64 = 1e-5; @@ -20,12 +20,12 @@ const N0_CUTOFF: f64 = 1e-9; #[derive(Clone)] pub struct PureFMTAssocFunctional { - parameters: Rc, + parameters: Arc, version: FMTVersion, } impl PureFMTAssocFunctional { - pub fn new(parameters: Rc, version: FMTVersion) -> Self { + pub fn new(parameters: Arc, version: FMTVersion) -> Self { Self { parameters, version, @@ -154,11 +154,11 @@ impl fmt::Display for PureFMTAssocFunctional { #[derive(Clone)] pub struct PureChainFunctional { - parameters: Rc, + parameters: Arc, } impl PureChainFunctional { - pub fn new(parameters: Rc) -> Self { + pub fn new(parameters: Arc) -> Self { Self { parameters } } } @@ -206,11 +206,11 @@ impl fmt::Display for PureChainFunctional { #[derive(Clone)] pub struct PureAttFunctional { - parameters: Rc, + parameters: Arc, } impl PureAttFunctional { - pub fn new(parameters: Rc) -> Self { + pub fn new(parameters: Arc) -> Self { Self { parameters } } } diff --git a/src/pcsaft/eos/dispersion.rs b/src/pcsaft/eos/dispersion.rs index 4a3c67c87..a8a4f0c9d 100644 --- a/src/pcsaft/eos/dispersion.rs +++ b/src/pcsaft/eos/dispersion.rs @@ -4,7 +4,7 @@ use feos_core::{HelmholtzEnergyDual, StateHD}; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_3, PI}; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; pub const A0: [f64; 7] = [ 0.91056314451539, @@ -62,7 +62,7 @@ pub const B2: [f64; 7] = [ ]; pub struct Dispersion { - pub parameters: Rc, + pub parameters: Arc, } impl> HelmholtzEnergyDual for Dispersion { diff --git a/src/pcsaft/eos/hard_chain.rs b/src/pcsaft/eos/hard_chain.rs index f4da72970..ea604e26d 100644 --- a/src/pcsaft/eos/hard_chain.rs +++ b/src/pcsaft/eos/hard_chain.rs @@ -4,10 +4,10 @@ use feos_core::{HelmholtzEnergyDual, StateHD}; use ndarray::Array; use num_dual::*; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; pub struct HardChain { - pub parameters: Rc, + pub parameters: Arc, } impl> HelmholtzEnergyDual for HardChain { diff --git a/src/pcsaft/eos/mod.rs b/src/pcsaft/eos/mod.rs index a10476d3d..3e6d2b384 100644 --- a/src/pcsaft/eos/mod.rs +++ b/src/pcsaft/eos/mod.rs @@ -10,7 +10,7 @@ use feos_core::{ use ndarray::Array1; use quantity::si::*; use std::f64::consts::{FRAC_PI_6, PI}; -use std::rc::Rc; +use std::sync::Arc; pub(crate) mod dispersion; pub(crate) mod hard_chain; @@ -50,18 +50,18 @@ impl Default for PcSaftOptions { /// PC-SAFT equation of state. pub struct PcSaft { - parameters: Rc, + parameters: Arc, options: PcSaftOptions, contributions: Vec>, ideal_gas: IdealGasContributions, } impl PcSaft { - pub fn new(parameters: Rc) -> Self { + pub fn new(parameters: Arc) -> Self { Self::with_options(parameters, PcSaftOptions::default()) } - pub fn with_options(parameters: Rc, options: PcSaftOptions) -> Self { + pub fn with_options(parameters: Arc, options: PcSaftOptions) -> Self { let mut contributions: Vec> = Vec::with_capacity(7); contributions.push(Box::new(HardSphere::new(¶meters))); contributions.push(Box::new(HardChain { @@ -116,7 +116,7 @@ impl EquationOfState for PcSaft { fn subset(&self, component_list: &[usize]) -> Self { Self::with_options( - Rc::new(self.parameters.subset(component_list)), + Arc::new(self.parameters.subset(component_list)), self.options, ) } @@ -286,7 +286,7 @@ impl EntropyScaling for PcSaft { } let p = &self.parameters; let mws = self.molar_weight(); - let state = State::new_nvt(&Rc::new(Self::new(p.clone())), temperature, volume, moles)?; + let state = State::new_nvt(&Arc::new(Self::new(p.clone())), temperature, volume, moles)?; let res: Array1 = (0..self.components()) .map(|i| { let tr = (temperature / p.epsilon_k[i] / KELVIN) @@ -350,7 +350,7 @@ mod tests { #[test] fn ideal_gas_pressure() { - let e = Rc::new(PcSaft::new(propane_parameters())); + let e = Arc::new(PcSaft::new(propane_parameters())); let t = 200.0 * KELVIN; let v = 1e-3 * METER.powi(3); let n = arr1(&[1.0]) * MOL; @@ -366,7 +366,7 @@ mod tests { #[test] fn ideal_gas_heat_capacity_joback() { - let e = Rc::new(PcSaft::new(propane_parameters())); + let e = Arc::new(PcSaft::new(propane_parameters())); let t = 200.0 * KELVIN; let v = 1e-3 * METER.powi(3); let n = arr1(&[1.0]) * MOL; @@ -412,7 +412,7 @@ mod tests { #[test] fn association() { - let parameters = Rc::new(water_parameters()); + let parameters = Arc::new(water_parameters()); let assoc = Association::new(¶meters, ¶meters.association, 50, 1e-10); let t = 350.0; let v = 41.248289328513216; @@ -424,7 +424,7 @@ mod tests { #[test] fn cross_association() { - let parameters = Rc::new(water_parameters()); + let parameters = Arc::new(water_parameters()); let assoc = Association::new_cross_association(¶meters, ¶meters.association, 50, 1e-10); let t = 350.0; @@ -437,7 +437,7 @@ mod tests { #[test] fn new_tpn() { - let e = Rc::new(PcSaft::new(propane_parameters())); + let e = Arc::new(PcSaft::new(propane_parameters())); let t = 300.0 * KELVIN; let p = BAR; let m = arr1(&[1.0]) * MOL; @@ -452,7 +452,7 @@ mod tests { #[test] fn vle_pure() { - let e = Rc::new(PcSaft::new(propane_parameters())); + let e = Arc::new(PcSaft::new(propane_parameters())); let t = 300.0 * KELVIN; let vle = PhaseEquilibrium::pure(&e, t, None, Default::default()); if let Ok(v) = vle { @@ -466,7 +466,7 @@ mod tests { #[test] fn critical_point() { - let e = Rc::new(PcSaft::new(propane_parameters())); + let e = Arc::new(PcSaft::new(propane_parameters())); let t = 300.0 * KELVIN; let cp = State::critical_point(&e, None, Some(t), Default::default()); if let Ok(v) = cp { @@ -476,7 +476,7 @@ mod tests { #[test] fn speed_of_sound() { - let e = Rc::new(PcSaft::new(propane_parameters())); + let e = Arc::new(PcSaft::new(propane_parameters())); let t = 300.0 * KELVIN; let p = BAR; let m = arr1(&[1.0]) * MOL; @@ -490,9 +490,9 @@ mod tests { #[test] fn mix_single() { - let e1 = Rc::new(PcSaft::new(propane_parameters())); - let e2 = Rc::new(PcSaft::new(butane_parameters())); - let e12 = Rc::new(PcSaft::new(propane_butane_parameters())); + let e1 = Arc::new(PcSaft::new(propane_parameters())); + let e2 = Arc::new(PcSaft::new(butane_parameters())); + let e12 = Arc::new(PcSaft::new(propane_butane_parameters())); let t = 300.0 * KELVIN; let v = 0.02456883872966545 * METER.powi(3); let m1 = arr1(&[2.0]) * MOL; @@ -516,7 +516,7 @@ mod tests { #[test] fn viscosity() -> EosResult<()> { - let e = Rc::new(PcSaft::new(propane_parameters())); + let e = Arc::new(PcSaft::new(propane_parameters())); let t = 300.0 * KELVIN; let p = BAR; let n = arr1(&[1.0]) * MOL; @@ -539,7 +539,7 @@ mod tests { #[test] fn diffusion() -> EosResult<()> { - let e = Rc::new(PcSaft::new(propane_parameters())); + let e = Arc::new(PcSaft::new(propane_parameters())); let t = 300.0 * KELVIN; let p = BAR; let n = arr1(&[1.0]) * MOL; diff --git a/src/pcsaft/eos/polar.rs b/src/pcsaft/eos/polar.rs index 501fd189a..704da0059 100644 --- a/src/pcsaft/eos/polar.rs +++ b/src/pcsaft/eos/polar.rs @@ -5,7 +5,7 @@ use ndarray::prelude::*; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_3, PI}; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; pub const ALPHA: f64 = 1.1937350; @@ -183,7 +183,7 @@ fn triplet_integral_ijk_dq>(mijk: f64, eta: D, c: &[[f64; 2]]) - } pub struct Dipole { - pub parameters: Rc, + pub parameters: Arc, } impl> HelmholtzEnergyDual for Dipole { @@ -279,7 +279,7 @@ impl fmt::Display for Dipole { } pub struct Quadrupole { - pub parameters: Rc, + pub parameters: Arc, } impl> HelmholtzEnergyDual for Quadrupole { @@ -384,7 +384,7 @@ pub enum DQVariants { } pub struct DipoleQuadrupole { - pub parameters: Rc, + pub parameters: Arc, pub variant: DQVariants, } @@ -505,7 +505,7 @@ mod tests { #[test] fn test_dipolar_contribution() { let dp = Dipole { - parameters: Rc::new(dme_parameters()), + parameters: Arc::new(dme_parameters()), }; let t = 350.0; let v = 1000.0; @@ -518,7 +518,7 @@ mod tests { #[test] fn test_quadrupolar_contribution() { let qp = Quadrupole { - parameters: Rc::new(carbon_dioxide_parameters()), + parameters: Arc::new(carbon_dioxide_parameters()), }; let t = 350.0; let v = 1000.0; @@ -531,17 +531,17 @@ mod tests { #[test] fn test_dipolar_quadrupolar_contribution() { let dp = Dipole { - parameters: Rc::new(dme_co2_parameters()), + parameters: Arc::new(dme_co2_parameters()), }; let qp = Quadrupole { - parameters: Rc::new(dme_co2_parameters()), + parameters: Arc::new(dme_co2_parameters()), }; // let dpqp = DipoleQuadrupole { - // parameters: Rc::new(dme_co2_parameters()), + // parameters: Arc::new(dme_co2_parameters()), // variant: DQVariants::DQ35, // }; let disp = Dispersion { - parameters: Rc::new(dme_co2_parameters()), + parameters: Arc::new(dme_co2_parameters()), }; let t = 350.0; let v = 1000.0; diff --git a/src/pcsaft/eos/qspr.rs b/src/pcsaft/eos/qspr.rs index 98cf41975..baa2e6436 100644 --- a/src/pcsaft/eos/qspr.rs +++ b/src/pcsaft/eos/qspr.rs @@ -3,7 +3,7 @@ use feos_core::IdealGasContributionDual; use ndarray::Array1; use num_dual::*; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; const RGAS: f64 = 6.022140857 * 1.38064852; const KB: f64 = 1.38064852e-23; @@ -65,7 +65,7 @@ const AP_400: [f64; 6] = [ #[allow(clippy::upper_case_acronyms)] pub struct QSPR { - pub parameters: Rc, + pub parameters: Arc, } impl> IdealGasContributionDual for QSPR { diff --git a/src/pcsaft/parameters.rs b/src/pcsaft/parameters.rs index 20b7bf3d0..f2619c1f2 100644 --- a/src/pcsaft/parameters.rs +++ b/src/pcsaft/parameters.rs @@ -532,7 +532,7 @@ impl std::fmt::Display for PcSaftParameters { pub mod utils { use super::*; use feos_core::joback::JobackRecord; - use std::rc::Rc; + use std::sync::Arc; // use feos_core::parameter::SegmentRecord; // pub fn pure_record_vec() -> Vec> { @@ -696,7 +696,7 @@ pub mod utils { // PcSaftParameters::from_records(vec![methane_record], None).unwrap() // } - pub fn propane_parameters() -> Rc { + pub fn propane_parameters() -> Arc { let propane_json = r#" { "identifier": { @@ -719,7 +719,7 @@ pub mod utils { }"#; let propane_record: PureRecord = serde_json::from_str(propane_json).expect("Unable to parse json."); - Rc::new(PcSaftParameters::new_pure(propane_record)) + Arc::new(PcSaftParameters::new_pure(propane_record)) } // pub fn propane_homogc_parameters() -> PcSaftParameters { @@ -793,7 +793,7 @@ pub mod utils { PcSaftParameters::new_pure(co2_record) } - pub fn butane_parameters() -> Rc { + pub fn butane_parameters() -> Arc { let butane_json = r#" { "identifier": { @@ -813,7 +813,7 @@ pub mod utils { }"#; let butane_record: PureRecord = serde_json::from_str(butane_json).expect("Unable to parse json."); - Rc::new(PcSaftParameters::new_pure(butane_record)) + Arc::new(PcSaftParameters::new_pure(butane_record)) } pub fn dme_parameters() -> PcSaftParameters { @@ -907,7 +907,7 @@ pub mod utils { PcSaftParameters::new_binary(binary_record, None) } - pub fn propane_butane_parameters() -> Rc { + pub fn propane_butane_parameters() -> Arc { let binary_json = r#"[ { "identifier": { @@ -949,7 +949,7 @@ pub mod utils { ]"#; let binary_record: Vec> = serde_json::from_str(binary_json).expect("Unable to parse json."); - Rc::new(PcSaftParameters::new_binary(binary_record, None)) + Arc::new(PcSaftParameters::new_binary(binary_record, None)) } // pub fn water_hexane_parameters() -> PcSaftParameters { diff --git a/src/pcsaft/python.rs b/src/pcsaft/python.rs index 49757767d..50694f3c8 100644 --- a/src/pcsaft/python.rs +++ b/src/pcsaft/python.rs @@ -12,7 +12,7 @@ use numpy::{PyArray2, PyReadonlyArray2, ToPyArray}; use pyo3::exceptions::PyTypeError; use pyo3::prelude::*; use std::convert::{TryFrom, TryInto}; -use std::rc::Rc; +use std::sync::Arc; /// Create a set of PC-Saft parameters from records. #[pyclass(name = "PcSaftRecord")] @@ -125,7 +125,7 @@ impl_json_handling!(PyPcSaftRecord); impl_pure_record!(PcSaftRecord, PyPcSaftRecord, JobackRecord, PyJobackRecord); impl_segment_record!(PcSaftRecord, PyPcSaftRecord, JobackRecord, PyJobackRecord); -#[pyclass(name = "PcSaftBinaryRecord", unsendable)] +#[pyclass(name = "PcSaftBinaryRecord")] #[pyo3( text_signature = "(pure_records, binary_records=None, substances=None, search_option='Name')" )] @@ -151,12 +151,12 @@ impl_binary_record!(PcSaftBinaryRecord, PyPcSaftBinaryRecord); /// Returns /// ------- /// PcSaftParameters -#[pyclass(name = "PcSaftParameters", unsendable)] +#[pyclass(name = "PcSaftParameters")] #[pyo3( text_signature = "(pure_records, binary_records=None, substances=None, search_option='Name')" )] #[derive(Clone)] -pub struct PyPcSaftParameters(pub Rc); +pub struct PyPcSaftParameters(pub Arc); impl_parameter!(PcSaftParameters, PyPcSaftParameters); impl_parameter_from_segments!(PcSaftParameters, PyPcSaftParameters); diff --git a/src/pets/dft/dispersion.rs b/src/pets/dft/dispersion.rs index f53cd254b..79c648497 100644 --- a/src/pets/dft/dispersion.rs +++ b/src/pets/dft/dispersion.rs @@ -9,7 +9,7 @@ use ndarray::*; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_3, PI}; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; /// psi Parameter for DFT (Heier2018) const PSI_DFT: f64 = 1.21; @@ -18,11 +18,11 @@ const PSI_PDGT: f64 = 1.21; #[derive(Clone)] pub struct AttractiveFunctional { - parameters: Rc, + parameters: Arc, } impl AttractiveFunctional { - pub fn new(parameters: Rc) -> Self { + pub fn new(parameters: Arc) -> Self { Self { parameters } } } diff --git a/src/pets/dft/mod.rs b/src/pets/dft/mod.rs index 252852630..33e350ead 100644 --- a/src/pets/dft/mod.rs +++ b/src/pets/dft/mod.rs @@ -13,7 +13,7 @@ use num_dual::DualNum; use pure_pets_functional::*; use quantity::si::*; use std::f64::consts::FRAC_PI_6; -use std::rc::Rc; +use std::sync::Arc; mod dispersion; mod pure_pets_functional; @@ -21,7 +21,7 @@ mod pure_pets_functional; /// PeTS Helmholtz energy functional. pub struct PetsFunctional { /// PeTS parameters of all substances in the system - pub parameters: Rc, + pub parameters: Arc, fmt_version: FMTVersion, options: PetsOptions, contributions: Vec>, @@ -33,19 +33,19 @@ impl PetsFunctional { /// /// # Defaults /// `FMTVersion`: `FMTVersion::WhiteBear` - pub fn new(parameters: Rc) -> DFT { + pub fn new(parameters: Arc) -> DFT { Self::with_options(parameters, FMTVersion::WhiteBear, PetsOptions::default()) } /// PeTS functional with default options for and provided FMT version. #[allow(non_snake_case)] - pub fn new_full(parameters: Rc, fmt_Version: FMTVersion) -> DFT { + pub fn new_full(parameters: Arc, fmt_Version: FMTVersion) -> DFT { Self::with_options(parameters, fmt_Version, PetsOptions::default()) } /// PeTS functional with provided options for FMT and equation of state options. pub fn with_options( - parameters: Rc, + parameters: Arc, fmt_version: FMTVersion, pets_options: PetsOptions, ) -> DFT { @@ -93,7 +93,7 @@ impl PetsFunctional { impl HelmholtzEnergyFunctional for PetsFunctional { fn subset(&self, component_list: &[usize]) -> DFT { Self::with_options( - Rc::new(self.parameters.subset(component_list)), + Arc::new(self.parameters.subset(component_list)), self.fmt_version, self.options, ) diff --git a/src/pets/dft/pure_pets_functional.rs b/src/pets/dft/pure_pets_functional.rs index d17034984..b7aaaa798 100644 --- a/src/pets/dft/pure_pets_functional.rs +++ b/src/pets/dft/pure_pets_functional.rs @@ -9,19 +9,19 @@ use ndarray::*; use num_dual::*; use std::f64::consts::{FRAC_PI_6, PI}; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; const PI36M1: f64 = 1.0 / (36.0 * PI); const N3_CUTOFF: f64 = 1e-5; #[derive(Clone)] pub struct PureFMTFunctional { - parameters: Rc, + parameters: Arc, version: FMTVersion, } impl PureFMTFunctional { - pub fn new(parameters: Rc, version: FMTVersion) -> Self { + pub fn new(parameters: Arc, version: FMTVersion) -> Self { Self { parameters, version, @@ -116,11 +116,11 @@ impl fmt::Display for PureFMTFunctional { #[derive(Clone)] pub struct PureAttFunctional { - parameters: Rc, + parameters: Arc, } impl PureAttFunctional { - pub fn new(parameters: Rc) -> Self { + pub fn new(parameters: Arc) -> Self { Self { parameters } } } diff --git a/src/pets/eos/dispersion.rs b/src/pets/eos/dispersion.rs index 00e555209..09db1c17c 100644 --- a/src/pets/eos/dispersion.rs +++ b/src/pets/eos/dispersion.rs @@ -4,7 +4,7 @@ use feos_core::{HelmholtzEnergyDual, StateHD}; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_3, PI}; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; pub const A: [f64; 7] = [ 0.690603404, @@ -27,7 +27,7 @@ pub const B: [f64; 7] = [ #[derive(Debug, Clone)] pub struct Dispersion { - pub parameters: Rc, + pub parameters: Arc, } impl> HelmholtzEnergyDual for Dispersion { diff --git a/src/pets/eos/mod.rs b/src/pets/eos/mod.rs index c878da495..e1f8c9d4b 100644 --- a/src/pets/eos/mod.rs +++ b/src/pets/eos/mod.rs @@ -9,7 +9,7 @@ use feos_core::{ use ndarray::Array1; use quantity::si::*; use std::f64::consts::{FRAC_PI_6, PI}; -use std::rc::Rc; +use std::sync::Arc; pub(crate) mod dispersion; mod qspr; @@ -40,7 +40,7 @@ impl Default for PetsOptions { /// PeTS equation of state. pub struct Pets { - parameters: Rc, + parameters: Arc, options: PetsOptions, contributions: Vec>, ideal_gas: IdealGasContributions, @@ -48,12 +48,12 @@ pub struct Pets { impl Pets { /// PeTS equation of state with default options. - pub fn new(parameters: Rc) -> Self { + pub fn new(parameters: Arc) -> Self { Self::with_options(parameters, PetsOptions::default()) } /// PeTS equation of state with provided options. - pub fn with_options(parameters: Rc, options: PetsOptions) -> Self { + pub fn with_options(parameters: Arc, options: PetsOptions) -> Self { let contributions: Vec> = vec![ Box::new(HardSphere::new(¶meters)), Box::new(Dispersion { @@ -82,7 +82,7 @@ impl EquationOfState for Pets { fn subset(&self, component_list: &[usize]) -> Self { Self::with_options( - Rc::new(self.parameters.subset(component_list)), + Arc::new(self.parameters.subset(component_list)), self.options, ) } @@ -270,7 +270,7 @@ impl EntropyScaling for Pets { } let p = &self.parameters; let state = State::new_nvt( - &Rc::new(Self::new(self.parameters.clone())), + &Arc::new(Self::new(self.parameters.clone())), temperature, volume, moles, @@ -337,7 +337,7 @@ mod tests { #[test] fn ideal_gas_pressure() { - let e = Rc::new(Pets::new(argon_parameters())); + let e = Arc::new(Pets::new(argon_parameters())); let t = 200.0 * KELVIN; let v = 1e-3 * METER.powi(3); let n = arr1(&[1.0]) * MOL; @@ -353,7 +353,7 @@ mod tests { #[test] fn ideal_gas_heat_capacity_joback() { - let e = Rc::new(Pets::new(argon_parameters())); + let e = Arc::new(Pets::new(argon_parameters())); let t = 200.0 * KELVIN; let v = 1e-3 * METER.powi(3); let n = arr1(&[1.0]) * MOL; @@ -388,7 +388,7 @@ mod tests { #[test] fn new_tpn() { - let e = Rc::new(Pets::new(argon_parameters())); + let e = Arc::new(Pets::new(argon_parameters())); let t = 300.0 * KELVIN; let p = BAR; let m = arr1(&[1.0]) * MOL; @@ -403,7 +403,7 @@ mod tests { #[test] fn vle_pure_t() { - let e = Rc::new(Pets::new(argon_parameters())); + let e = Arc::new(Pets::new(argon_parameters())); let t = 300.0 * KELVIN; let vle = PhaseEquilibrium::pure(&e, t, None, Default::default()); if let Ok(v) = vle { @@ -417,7 +417,7 @@ mod tests { // #[test] // fn critical_point() { - // let e = Rc::new(Pets::new(argon_parameters())); + // let e = Arc::new(Pets::new(argon_parameters())); // let t = 300.0 * KELVIN; // let cp = State::critical_point(&e, None, Some(t), Default::default()); // if let Ok(v) = cp { @@ -427,7 +427,7 @@ mod tests { // #[test] // fn speed_of_sound() { - // let e = Rc::new(Pets::new(argon_parameters())); + // let e = Arc::new(Pets::new(argon_parameters())); // let t = 300.0 * KELVIN; // let p = BAR; // let m = arr1(&[1.0]) * MOL; @@ -441,9 +441,9 @@ mod tests { #[test] fn mix_single() { - let e1 = Rc::new(Pets::new(argon_parameters())); - let e2 = Rc::new(Pets::new(krypton_parameters())); - let e12 = Rc::new(Pets::new(argon_krypton_parameters())); + let e1 = Arc::new(Pets::new(argon_parameters())); + let e2 = Arc::new(Pets::new(krypton_parameters())); + let e12 = Arc::new(Pets::new(argon_krypton_parameters())); let t = 300.0 * KELVIN; let v = 0.02456883872966545 * METER.powi(3); let m1 = arr1(&[2.0]) * MOL; @@ -467,7 +467,7 @@ mod tests { // #[test] // fn viscosity() -> EosResult<()> { - // let e = Rc::new(Pets::new(argon_parameters())); + // let e = Arc::new(Pets::new(argon_parameters())); // let t = 300.0 * KELVIN; // let p = BAR; // let n = arr1(&[1.0]) * MOL; @@ -490,7 +490,7 @@ mod tests { // #[test] // fn diffusion() -> EosResult<()> { - // let e = Rc::new(Pets::new(argon_parameters())); + // let e = Arc::new(Pets::new(argon_parameters())); // let t = 300.0 * KELVIN; // let p = BAR; // let n = arr1(&[1.0]) * MOL; diff --git a/src/pets/eos/qspr.rs b/src/pets/eos/qspr.rs index e4345213b..78f9e8572 100644 --- a/src/pets/eos/qspr.rs +++ b/src/pets/eos/qspr.rs @@ -3,7 +3,7 @@ use feos_core::IdealGasContributionDual; use ndarray::Array1; use num_dual::*; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; const RGAS: f64 = 6.022140857 * 1.38064852; const KB: f64 = 1.38064852e-23; @@ -65,7 +65,7 @@ const NA_NP_400: [f64; 6] = [ #[allow(clippy::upper_case_acronyms)] pub struct QSPR { - pub parameters: Rc, + pub parameters: Arc, } impl> IdealGasContributionDual for QSPR { diff --git a/src/pets/mod.rs b/src/pets/mod.rs index 373ef8726..c79c12580 100644 --- a/src/pets/mod.rs +++ b/src/pets/mod.rs @@ -1,9 +1,9 @@ //! Perturbed Truncated and Shifted (PeTS) equation of state -//! +//! //! [Heier et al. (2018)](https://doi.org/10.1080/00268976.2018.1447153) -//! +//! //! PeTS is an equation of state for the truncated and shifted Lennar-Jones fluid with cut-off -//! distance 2.5 $\sigma$. +//! distance 2.5 $\sigma$. //! It utilizes a hard-sphere fluid as reference with an attactive perturbation. #![warn(clippy::all)] #![allow(clippy::too_many_arguments)] @@ -18,7 +18,7 @@ mod parameters; #[cfg(feature = "dft")] pub use dft::PetsFunctional; pub use eos::{Pets, PetsOptions}; -pub use parameters::{PetsParameters, PetsRecord, PetsBinaryRecord}; +pub use parameters::{PetsBinaryRecord, PetsParameters, PetsRecord}; #[cfg(feature = "python")] pub mod python; diff --git a/src/pets/parameters.rs b/src/pets/parameters.rs index 39874e665..096f6e724 100644 --- a/src/pets/parameters.rs +++ b/src/pets/parameters.rs @@ -286,9 +286,9 @@ impl std::fmt::Display for PetsParameters { pub mod utils { use super::*; use feos_core::joback::JobackRecord; - use std::rc::Rc; + use std::sync::Arc; - pub fn argon_parameters() -> Rc { + pub fn argon_parameters() -> Arc { let argon_json = r#" { "identifier": { @@ -310,10 +310,10 @@ pub mod utils { }"#; let argon_record: PureRecord = serde_json::from_str(argon_json).expect("Unable to parse json."); - Rc::new(PetsParameters::new_pure(argon_record)) + Arc::new(PetsParameters::new_pure(argon_record)) } - pub fn krypton_parameters() -> Rc { + pub fn krypton_parameters() -> Arc { let krypton_json = r#" { "identifier": { @@ -332,10 +332,10 @@ pub mod utils { }"#; let krypton_record: PureRecord = serde_json::from_str(krypton_json).expect("Unable to parse json."); - Rc::new(PetsParameters::new_pure(krypton_record)) + Arc::new(PetsParameters::new_pure(krypton_record)) } - pub fn argon_krypton_parameters() -> Rc { + pub fn argon_krypton_parameters() -> Arc { let binary_json = r#"[ { "identifier": { @@ -376,6 +376,6 @@ pub mod utils { ]"#; let binary_record: Vec> = serde_json::from_str(binary_json).expect("Unable to parse json."); - Rc::new(PetsParameters::new_binary(binary_record, None)) + Arc::new(PetsParameters::new_binary(binary_record, None)) } } diff --git a/src/pets/python.rs b/src/pets/python.rs index 3f0ebac4f..f0b365fa1 100644 --- a/src/pets/python.rs +++ b/src/pets/python.rs @@ -9,10 +9,10 @@ use numpy::{PyArray2, PyReadonlyArray2, ToPyArray}; use pyo3::exceptions::{PyTypeError, PyValueError}; use pyo3::prelude::*; use std::convert::{TryFrom, TryInto}; -use std::rc::Rc; +use std::sync::Arc; /// Create a set of PeTS parameters from records. -#[pyclass(name = "PetsRecord", unsendable)] +#[pyclass(name = "PetsRecord")] #[pyo3( text_signature = "(sigma, epsilon_k, viscosity=None, diffusion=None, thermal_conductivity=None)" )] @@ -71,7 +71,7 @@ impl PyPetsRecord { impl_json_handling!(PyPetsRecord); impl_pure_record!(PetsRecord, PyPetsRecord, JobackRecord, PyJobackRecord); -#[pyclass(name = "PetsBinaryRecord", unsendable)] +#[pyclass(name = "PetsBinaryRecord")] #[pyo3( text_signature = "(pure_records, binary_records=None, substances=None, search_option='Name')" )] @@ -93,12 +93,12 @@ impl_binary_record!(PetsBinaryRecord, PyPetsBinaryRecord); /// When not provided, all entries of `pure_records` are used. /// search_option : {'Name', 'Cas', 'Inchi', 'IupacName', 'Formula', 'Smiles'}, optional, defaults to 'Name'. /// Identifier that is used to search substance. -#[pyclass(name = "PetsParameters", unsendable)] +#[pyclass(name = "PetsParameters")] #[pyo3( text_signature = "(pure_records, binary_records=None, substances=None, search_option='Name')" )] #[derive(Clone)] -pub struct PyPetsParameters(pub Rc); +pub struct PyPetsParameters(pub Arc); #[pymethods] impl PyPetsParameters { @@ -191,7 +191,7 @@ impl PyPetsParameters { None => Array2::from_shape_fn((n, n), |(_, _)| PetsBinaryRecord::from(0.0)), }; - Ok(Self(Rc::new(PetsParameters::from_records( + Ok(Self(Arc::new(PetsParameters::from_records( pure_records, binary, )))) @@ -241,7 +241,7 @@ impl PyPetsParameters { PetsRecord::new(sigma, epsilon_k, viscosity, diffusion, thermal_conductivity), None, ); - Self(Rc::new(PetsParameters::new_pure(pure_record))) + Self(Arc::new(PetsParameters::new_pure(pure_record))) } #[getter] diff --git a/src/python/dft.rs b/src/python/dft.rs index 46426db24..1965a6cc3 100644 --- a/src/python/dft.rs +++ b/src/python/dft.rs @@ -31,11 +31,11 @@ use pyo3::wrap_pymodule; use quantity::python::*; use quantity::si::*; use std::collections::HashMap; -use std::rc::Rc; +use std::sync::Arc; -#[pyclass(name = "HelmholtzEnergyFunctional", unsendable)] +#[pyclass(name = "HelmholtzEnergyFunctional")] #[derive(Clone)] -pub struct PyFunctionalVariant(pub Rc>); +pub struct PyFunctionalVariant(pub Arc>); #[pymethods] impl PyFunctionalVariant { @@ -85,7 +85,7 @@ impl PyFunctionalVariant { tol_cross_assoc, dq_variant, }; - Self(Rc::new( + Self(Arc::new( PcSaftFunctional::with_options(parameters.0, fmt_version, options).into(), )) } @@ -131,7 +131,7 @@ impl PyFunctionalVariant { max_iter_cross_assoc, tol_cross_assoc, }; - Self(Rc::new( + Self(Arc::new( GcPcSaftFunctional::with_options(parameters.0, fmt_version, options).into(), )) } @@ -157,7 +157,7 @@ impl PyFunctionalVariant { #[pyo3(text_signature = "(parameters, fmt_version, max_eta)")] fn pets(parameters: PyPetsParameters, fmt_version: FMTVersion, max_eta: f64) -> Self { let options = PetsOptions { max_eta }; - Self(Rc::new( + Self(Arc::new( PetsFunctional::with_options(parameters.0, fmt_version, options).into(), )) } @@ -177,7 +177,7 @@ impl PyFunctionalVariant { #[staticmethod] #[pyo3(text_signature = "(sigma, version)")] fn fmt(sigma: &PyArray1, fmt_version: FMTVersion) -> Self { - Self(Rc::new( + Self(Arc::new( FMTFunctional::new(&sigma.to_owned_array(), fmt_version).into(), )) } diff --git a/src/python/eos.rs b/src/python/eos.rs index c2605210b..1c44ebd83 100644 --- a/src/python/eos.rs +++ b/src/python/eos.rs @@ -34,11 +34,11 @@ use pyo3::wrap_pymodule; use quantity::python::*; use quantity::si::*; use std::collections::HashMap; -use std::rc::Rc; +use std::sync::Arc; -#[pyclass(name = "EquationOfState", unsendable)] +#[pyclass(name = "EquationOfState")] #[derive(Clone)] -pub struct PyEosVariant(pub Rc); +pub struct PyEosVariant(pub Arc); #[pymethods] impl PyEosVariant { @@ -86,7 +86,7 @@ impl PyEosVariant { tol_cross_assoc, dq_variant, }; - Self(Rc::new(EosVariant::PcSaft(PcSaft::with_options( + Self(Arc::new(EosVariant::PcSaft(PcSaft::with_options( parameters.0, options, )))) @@ -129,7 +129,7 @@ impl PyEosVariant { max_iter_cross_assoc, tol_cross_assoc, }; - Self(Rc::new(EosVariant::GcPcSaft(GcPcSaft::with_options( + Self(Arc::new(EosVariant::GcPcSaft(GcPcSaft::with_options( parameters.0, options, )))) @@ -150,7 +150,7 @@ impl PyEosVariant { #[staticmethod] #[pyo3(text_signature = "(parameters)")] pub fn peng_robinson(parameters: PyPengRobinsonParameters) -> Self { - Self(Rc::new(EosVariant::PengRobinson(PengRobinson::new( + Self(Arc::new(EosVariant::PengRobinson(PengRobinson::new( parameters.0, )))) } @@ -169,7 +169,7 @@ impl PyEosVariant { #[staticmethod] #[pyo3(text_signature = "(obj)")] fn python(obj: Py) -> PyResult { - Ok(Self(Rc::new(EosVariant::Python(PyEoSObj::new(obj)?)))) + Ok(Self(Arc::new(EosVariant::Python(PyEoSObj::new(obj)?)))) } /// PeTS equation of state. @@ -192,7 +192,7 @@ impl PyEosVariant { #[pyo3(text_signature = "(parameters, max_eta)")] fn pets(parameters: PyPetsParameters, max_eta: f64) -> Self { let options = PetsOptions { max_eta }; - Self(Rc::new(EosVariant::Pets(Pets::with_options( + Self(Arc::new(EosVariant::Pets(Pets::with_options( parameters.0, options, )))) @@ -223,7 +223,7 @@ impl PyEosVariant { max_eta, perturbation, }; - Self(Rc::new(EosVariant::UVTheory(UVTheory::with_options( + Self(Arc::new(EosVariant::UVTheory(UVTheory::with_options( parameters.0, options, )))) diff --git a/src/python/mod.rs b/src/python/mod.rs index 9c8c94f04..bc0142571 100644 --- a/src/python/mod.rs +++ b/src/python/mod.rs @@ -24,6 +24,7 @@ use dft::__PYO3_PYMODULE_DEF_DFT; pub fn feos(py: Python<'_>, m: &PyModule) -> PyResult<()> { m.add("__version__", env!("CARGO_PKG_VERSION"))?; m.add_wrapped(wrap_pymodule!(quantity))?; + m.add_wrapped(wrap_pymodule!(eos))?; #[cfg(feature = "dft")] m.add_wrapped(wrap_pymodule!(dft))?; diff --git a/src/uvtheory/eos/attractive_perturbation_bh.rs b/src/uvtheory/eos/attractive_perturbation_bh.rs index 861b59445..a157724fc 100644 --- a/src/uvtheory/eos/attractive_perturbation_bh.rs +++ b/src/uvtheory/eos/attractive_perturbation_bh.rs @@ -6,7 +6,7 @@ use num_dual::DualNum; use std::{ f64::consts::{FRAC_PI_3, PI}, fmt, - rc::Rc, + sync::Arc, }; const C_BH: [[f64; 4]; 2] = [ @@ -41,7 +41,7 @@ const C2: [[f64; 2]; 3] = [ #[derive(Debug, Clone)] pub struct AttractivePerturbationBH { - pub parameters: Rc, + pub parameters: Arc, } impl fmt::Display for AttractivePerturbationBH { @@ -260,7 +260,7 @@ mod test { let p = methane_parameters(24.0, 6.0); let pt = AttractivePerturbationBH { - parameters: Rc::new(p.clone()), + parameters: Arc::new(p.clone()), }; let state = StateHD::new( reduced_temperature * p.epsilon_k[0], diff --git a/src/uvtheory/eos/attractive_perturbation_wca.rs b/src/uvtheory/eos/attractive_perturbation_wca.rs index b591a3645..99d304a46 100644 --- a/src/uvtheory/eos/attractive_perturbation_wca.rs +++ b/src/uvtheory/eos/attractive_perturbation_wca.rs @@ -3,7 +3,7 @@ use crate::uvtheory::parameters::*; use feos_core::{HelmholtzEnergyDual, StateHD}; use ndarray::Array1; use num_dual::DualNum; -use std::{f64::consts::PI, fmt, rc::Rc}; +use std::{f64::consts::PI, fmt, sync::Arc}; const C_WCA: [[f64; 6]; 6] = [ [ @@ -68,7 +68,7 @@ const C2: [[f64; 2]; 3] = [ #[derive(Debug, Clone)] pub struct AttractivePerturbationWCA { - pub parameters: Rc, + pub parameters: Arc, } impl fmt::Display for AttractivePerturbationWCA { @@ -306,7 +306,7 @@ mod test { let p = methane_parameters(24.0, 6.0); let pt = AttractivePerturbationWCA { - parameters: Rc::new(p.clone()), + parameters: Arc::new(p.clone()), }; let state = StateHD::new( reduced_temperature * p.epsilon_k[0], @@ -423,7 +423,7 @@ mod test { assert_relative_eq!(delta_b2, -4.7846399638747954, epsilon = 1e-6); // Full attractive contribution let pt = AttractivePerturbationWCA { - parameters: Rc::new(p), + parameters: Arc::new(p), }; let a = pt.helmholtz_energy(&state) / (moles[0] + moles[1]); @@ -489,7 +489,7 @@ mod test { // Full attractive contribution let pt = AttractivePerturbationWCA { - parameters: Rc::new(p), + parameters: Arc::new(p), }; let a = pt.helmholtz_energy(&state) / (moles[0] + moles[1]); assert_relative_eq!(a, -1.3318659166866607, epsilon = 1e-5); diff --git a/src/uvtheory/eos/hard_sphere_bh.rs b/src/uvtheory/eos/hard_sphere_bh.rs index 50e9a7ad9..4453310b9 100644 --- a/src/uvtheory/eos/hard_sphere_bh.rs +++ b/src/uvtheory/eos/hard_sphere_bh.rs @@ -4,7 +4,7 @@ use lazy_static::lazy_static; use ndarray::prelude::*; use num_dual::DualNum; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; lazy_static! { static ref BH_CONSTANTS_ETA_B: Array2 = arr2(&[ @@ -22,7 +22,7 @@ lazy_static! { #[derive(Debug, Clone)] pub struct HardSphere { - pub parameters: Rc, + pub parameters: Arc, } impl> HelmholtzEnergyDual for HardSphere { @@ -176,7 +176,7 @@ mod test { // let reduced_temperature = 4.0; // let hs = HardSphere { - // parameters: Rc::new(p.clone()), + // parameters: Arc::new(p.clone()), // }; // let particles = arr1(&[1000.0]); // let n = &particles / 6.02214076e23; diff --git a/src/uvtheory/eos/hard_sphere_wca.rs b/src/uvtheory/eos/hard_sphere_wca.rs index c20500179..5d256fa2b 100644 --- a/src/uvtheory/eos/hard_sphere_wca.rs +++ b/src/uvtheory/eos/hard_sphere_wca.rs @@ -5,7 +5,7 @@ use ndarray::prelude::*; use num_dual::DualNum; use std::f64::consts::PI; use std::fmt; -use std::rc::Rc; +use std::sync::Arc; lazy_static! { static ref WCA_CONSTANTS_ETA_B: Array2 = arr2(&[ @@ -38,7 +38,7 @@ lazy_static! { #[derive(Debug, Clone)] pub struct HardSphereWCA { - pub parameters: Rc, + pub parameters: Arc, } impl> HelmholtzEnergyDual for HardSphereWCA { @@ -295,7 +295,7 @@ mod test { ); let pt = HardSphereWCA { - parameters: Rc::new(p), + parameters: Arc::new(p), }; let state = StateHD::new(reduced_temperature, reduced_volume, moles.clone()); let a = pt.helmholtz_energy(&state) / (moles[0] + moles[1]); diff --git a/src/uvtheory/eos/mod.rs b/src/uvtheory/eos/mod.rs index d19e1a3de..8b59fd532 100644 --- a/src/uvtheory/eos/mod.rs +++ b/src/uvtheory/eos/mod.rs @@ -5,7 +5,7 @@ use super::parameters::UVParameters; use feos_core::{parameter::Parameter, EquationOfState, HelmholtzEnergy}; use ndarray::Array1; use std::f64::consts::FRAC_PI_6; -use std::rc::Rc; +use std::sync::Arc; pub(crate) mod attractive_perturbation_bh; pub(crate) mod attractive_perturbation_wca; @@ -47,19 +47,19 @@ impl Default for UVTheoryOptions { /// uv-theory equation of state pub struct UVTheory { - parameters: Rc, + parameters: Arc, options: UVTheoryOptions, contributions: Vec>, } impl UVTheory { /// uv-theory with default options (WCA). - pub fn new(parameters: Rc) -> Self { + pub fn new(parameters: Arc) -> Self { Self::with_options(parameters, UVTheoryOptions::default()) } /// uv-theory with provided options. - pub fn with_options(parameters: Rc, options: UVTheoryOptions) -> Self { + pub fn with_options(parameters: Arc, options: UVTheoryOptions) -> Self { let mut contributions: Vec> = Vec::with_capacity(3); match options.perturbation { @@ -102,7 +102,7 @@ impl EquationOfState for UVTheory { fn subset(&self, component_list: &[usize]) -> Self { Self::with_options( - Rc::new(self.parameters.subset(component_list)), + Arc::new(self.parameters.subset(component_list)), self.options.clone(), ) } @@ -138,7 +138,7 @@ mod test { let i = Identifier::new(None, None, None, None, None, None); let pr = PureRecord::new(i, 1.0, r, None); let parameters = UVParameters::new_pure(pr); - let eos = Rc::new(UVTheory::new(Rc::new(parameters))); + let eos = Arc::new(UVTheory::new(Arc::new(parameters))); let reduced_temperature = 4.0; //let reduced_temperature = 1.0; @@ -168,7 +168,7 @@ mod test { max_eta: 0.5, perturbation: Perturbation::BarkerHenderson, }; - let eos = Rc::new(UVTheory::with_options(Rc::new(parameters), options)); + let eos = Arc::new(UVTheory::with_options(Arc::new(parameters), options)); let reduced_temperature = 4.0; let reduced_density = 1.0; @@ -221,7 +221,7 @@ mod test { perturbation: Perturbation::BarkerHenderson, }; - let eos_bh = Rc::new(UVTheory::with_options(Rc::new(uv_parameters), options)); + let eos_bh = Arc::new(UVTheory::with_options(Arc::new(uv_parameters), options)); let state_bh = State::new_nvt(&eos_bh, t_x, volume, &moles).unwrap(); let a_bh = state_bh @@ -249,7 +249,7 @@ mod test { let volume = (p.sigma[0] * ANGSTROM).powi(3) / reduced_density * NAV * total_moles; // EoS - let eos_wca = Rc::new(UVTheory::new(Rc::new(p))); + let eos_wca = Arc::new(UVTheory::new(Arc::new(p))); let state_wca = State::new_nvt(&eos_wca, t_x, volume, &moles).unwrap(); let a_wca = state_wca .molar_helmholtz_energy(Contributions::ResidualNvt) @@ -278,7 +278,7 @@ mod test { let volume = NAV * total_moles / density; // EoS - let eos_wca = Rc::new(UVTheory::new(Rc::new(p))); + let eos_wca = Arc::new(UVTheory::new(Arc::new(p))); let state_wca = State::new_nvt(&eos_wca, t_x, volume, &moles).unwrap(); let a_wca = state_wca .molar_helmholtz_energy(Contributions::ResidualNvt) diff --git a/src/uvtheory/eos/reference_perturbation_bh.rs b/src/uvtheory/eos/reference_perturbation_bh.rs index f09f0cda0..e54ca0c2e 100644 --- a/src/uvtheory/eos/reference_perturbation_bh.rs +++ b/src/uvtheory/eos/reference_perturbation_bh.rs @@ -5,11 +5,11 @@ use crate::uvtheory::parameters::*; use feos_core::{HelmholtzEnergyDual, StateHD}; use num_dual::DualNum; use std::fmt; -use std::{f64::consts::PI, rc::Rc}; +use std::{f64::consts::PI, sync::Arc}; #[derive(Debug, Clone)] pub struct ReferencePerturbationBH { - pub parameters: Rc, + pub parameters: Arc, } impl fmt::Display for ReferencePerturbationBH { @@ -62,7 +62,7 @@ mod test { let p = test_parameters(24.0, 6.0, 1.0, 1.0); let pt = ReferencePerturbationBH { - parameters: Rc::new(p), + parameters: Arc::new(p), }; let state = StateHD::new(reduced_temperature, reduced_volume, moles.clone()); let a = pt.helmholtz_energy(&state) / moles[0]; diff --git a/src/uvtheory/eos/reference_perturbation_wca.rs b/src/uvtheory/eos/reference_perturbation_wca.rs index 1a34bb85e..89fb58801 100644 --- a/src/uvtheory/eos/reference_perturbation_wca.rs +++ b/src/uvtheory/eos/reference_perturbation_wca.rs @@ -1,18 +1,16 @@ use super::hard_sphere_wca::{ - diameter_wca, dimensionless_diameter_q_wca, packing_fraction, packing_fraction_a, packing_fraction_b, - }; use crate::uvtheory::parameters::*; use feos_core::{HelmholtzEnergyDual, StateHD}; use num_dual::DualNum; use std::fmt; -use std::{f64::consts::PI, rc::Rc}; +use std::{f64::consts::PI, sync::Arc}; #[derive(Debug, Clone)] pub struct ReferencePerturbationWCA { - pub parameters: Rc, + pub parameters: Arc, } impl fmt::Display for ReferencePerturbationWCA { @@ -47,7 +45,6 @@ impl> HelmholtzEnergyDual for ReferencePerturbationWCA { let q_ij = dimensionless_diameter_q_wca(t_ij, D::from(rep_ij), D::from(att_ij)) * p.sigma_ij[[i, j]]; - a += x[i] * x[j] * ((-eta_a[[i, j]] * 0.5 + 1.0) / (-eta_a[[i, j]] + 1.0).powi(3) @@ -80,7 +77,7 @@ mod test { let p = test_parameters(24.0, 6.0, 1.0, 1.0); let pt = ReferencePerturbationWCA { - parameters: Rc::new(p), + parameters: Arc::new(p), }; let state = StateHD::new(reduced_temperature, reduced_volume, moles.clone()); let a = pt.helmholtz_energy(&state) / moles[0]; @@ -88,7 +85,6 @@ mod test { } #[test] fn test_delta_a0_wca_mixture() { - let moles = arr1(&[0.40000000000000002, 0.59999999999999998]); let reduced_temperature = 1.0; let reduced_density = 0.90000000000000002; @@ -102,12 +98,11 @@ mod test { ); let pt = ReferencePerturbationWCA { - parameters: Rc::new(p), + parameters: Arc::new(p), }; let state = StateHD::new(reduced_temperature, reduced_volume, moles.clone()); let a = pt.helmholtz_energy(&state) / (moles[0] + moles[1]); assert_relative_eq!(a, 0.308268896386771, epsilon = 1e-6); - } } diff --git a/src/uvtheory/eos/ufraction.rs b/src/uvtheory/eos/ufraction.rs index 62ec965dc..21e789809 100644 --- a/src/uvtheory/eos/ufraction.rs +++ b/src/uvtheory/eos/ufraction.rs @@ -3,4 +3,4 @@ use num_dual::DualNum; trait UFraction> { fn u_fraction(&self, state: &StateHD) -> D; -} \ No newline at end of file +} diff --git a/src/uvtheory/mod.rs b/src/uvtheory/mod.rs index ca0740641..481932928 100644 --- a/src/uvtheory/mod.rs +++ b/src/uvtheory/mod.rs @@ -1,5 +1,5 @@ //! uv-theory equation of state -//! +//! //! [van Westen et al. (2021)](https://doi.org/10.1063/5.0073572) mod eos; mod parameters; diff --git a/src/uvtheory/python.rs b/src/uvtheory/python.rs index b66c0def0..b442971d7 100644 --- a/src/uvtheory/python.rs +++ b/src/uvtheory/python.rs @@ -10,15 +10,15 @@ use numpy::PyReadonlyArray2; use pyo3::exceptions::PyTypeError; use pyo3::prelude::*; use std::convert::{TryFrom, TryInto}; -use std::rc::Rc; +use std::sync::Arc; /// Create a set of UV Theory parameters from records. -#[pyclass(name = "NoRecord", unsendable)] +#[pyclass(name = "NoRecord")] #[derive(Clone)] struct PyNoRecord(NoRecord); /// Create a set of UV Theory parameters from records. -#[pyclass(name = "UVRecord", unsendable)] +#[pyclass(name = "UVRecord")] #[pyo3(text_signature = "(rep, att, sigma, epsilon_k)")] #[derive(Clone)] pub struct PyUVRecord(UVRecord); @@ -37,7 +37,7 @@ impl PyUVRecord { impl_json_handling!(PyUVRecord); -#[pyclass(name = "UVBinaryRecord", unsendable)] +#[pyclass(name = "UVBinaryRecord")] #[derive(Clone)] pub struct PyUVBinaryRecord(UVBinaryRecord); impl_binary_record!(UVBinaryRecord, PyUVBinaryRecord); @@ -56,10 +56,10 @@ impl_binary_record!(UVBinaryRecord, PyUVBinaryRecord); /// When not provided, all entries of `pure_records` are used. /// search_option : IdentifierOption, optional, defaults to IdentifierOption.Name /// Identifier that is used to search binary records. -#[pyclass(name = "UVParameters", unsendable)] +#[pyclass(name = "UVParameters")] #[pyo3(text_signature = "(pure_records, binary_records, substances, search_option)")] #[derive(Clone)] -pub struct PyUVParameters(pub Rc); +pub struct PyUVParameters(pub Arc); #[pymethods] impl PyUVParameters { @@ -98,7 +98,7 @@ impl PyUVParameters { }) .collect(); let binary = Array2::from_shape_fn((n, n), |(_, _)| UVBinaryRecord { k_ij: 0.0 }); - Self(Rc::new(UVParameters::from_records(pure_records, binary))) + Self(Arc::new(UVParameters::from_records(pure_records, binary))) } } diff --git a/tests/gc_pcsaft/binary.rs b/tests/gc_pcsaft/binary.rs index 07976b6f7..a9430ca7c 100644 --- a/tests/gc_pcsaft/binary.rs +++ b/tests/gc_pcsaft/binary.rs @@ -6,7 +6,7 @@ use feos_core::parameter::{IdentifierOption, ParameterHetero}; use feos_core::{EosResult, State}; use ndarray::arr1; use quantity::si::{KELVIN, MOL}; -use std::rc::Rc; +use std::sync::Arc; #[test] fn test_binary() -> EosResult<()> { @@ -27,9 +27,9 @@ fn test_binary() -> EosResult<()> { IdentifierOption::Name, ) .unwrap(); - let eos = Rc::new(GcPcSaft::new(Rc::new(parameters))); + let eos = Arc::new(GcPcSaft::new(Arc::new(parameters))); #[cfg(feature = "dft")] - let func = Rc::new(GcPcSaftFunctional::new(Rc::new(parameters_func))); + let func = Arc::new(GcPcSaftFunctional::new(Arc::new(parameters_func))); let moles = arr1(&[0.5, 0.5]) * MOL; let cp = State::critical_point(&eos, Some(&moles), None, Default::default())?; #[cfg(feature = "dft")] diff --git a/tests/gc_pcsaft/dft.rs b/tests/gc_pcsaft/dft.rs index cf0a77a98..6cee26de1 100644 --- a/tests/gc_pcsaft/dft.rs +++ b/tests/gc_pcsaft/dft.rs @@ -14,7 +14,7 @@ use feos_dft::{DFTSolver, Geometry}; use ndarray::arr1; use quantity::si::*; use std::error::Error; -use std::rc::Rc; +use std::sync::Arc; #[test] #[allow(non_snake_case)] @@ -41,8 +41,8 @@ fn test_bulk_implementation() -> Result<(), Box> { ) .unwrap(); - let eos = Rc::new(GcPcSaft::new(Rc::new(parameters))); - let func = Rc::new(GcPcSaftFunctional::new(Rc::new(parameters_func))); + let eos = Arc::new(GcPcSaft::new(Arc::new(parameters))); + let func = Arc::new(GcPcSaftFunctional::new(Arc::new(parameters_func))); let t = 200.0 * KELVIN; let v = 0.002 * METER.powi(3) * NAV / NAV_old; let n = arr1(&[1.5]) * MOL; @@ -115,18 +115,18 @@ fn test_bulk_association() -> Result<(), Box> { vec!["OH".into(), "CH2".into(), "CH2".into(), "OH".into()], None, ); - let eos_parameters = Rc::new(GcPcSaftEosParameters::from_segments( + let eos_parameters = Arc::new(GcPcSaftEosParameters::from_segments( vec![ethylene_glycol.clone()], segment_records.clone(), None, )?); - let eos = Rc::new(GcPcSaft::new(eos_parameters.clone())); - let func_parameters = Rc::new(GcPcSaftFunctionalParameters::from_segments( + let eos = Arc::new(GcPcSaft::new(eos_parameters.clone())); + let func_parameters = Arc::new(GcPcSaftFunctionalParameters::from_segments( vec![ethylene_glycol], segment_records, None, )?); - let func = Rc::new(GcPcSaftFunctional::new(func_parameters.clone())); + let func = Arc::new(GcPcSaftFunctional::new(func_parameters.clone())); let t = 200.0 * KELVIN; let v = 0.002 * METER.powi(3); @@ -187,7 +187,7 @@ fn test_dft() -> Result<(), Box> { ) .unwrap(); - let func = Rc::new(GcPcSaftFunctional::new(Rc::new(parameters))); + let func = Arc::new(GcPcSaftFunctional::new(Arc::new(parameters))); let t = 200.0 * KELVIN; let w = 150.0 * ANGSTROM; let points = 2048; @@ -233,7 +233,7 @@ fn test_dft_assoc() -> Result<(), Box> { ) .unwrap(); - let func = Rc::new(GcPcSaftFunctional::new(Rc::new(parameters))); + let func = Arc::new(GcPcSaftFunctional::new(Arc::new(parameters))); let t = 300.0 * KELVIN; let w = 100.0 * ANGSTROM; let points = 4096; diff --git a/tests/pcsaft/critical_point.rs b/tests/pcsaft/critical_point.rs index ad75809e4..d9425ecbb 100644 --- a/tests/pcsaft/critical_point.rs +++ b/tests/pcsaft/critical_point.rs @@ -5,7 +5,7 @@ use feos_core::State; use ndarray::arr1; use quantity::si::*; use std::error::Error; -use std::rc::Rc; +use std::sync::Arc; #[test] fn test_critical_point_pure() -> Result<(), Box> { @@ -15,7 +15,7 @@ fn test_critical_point_pure() -> Result<(), Box> { None, IdentifierOption::Name, )?; - let saft = Rc::new(PcSaft::new(Rc::new(params))); + let saft = Arc::new(PcSaft::new(Arc::new(params))); let t = 300.0 * KELVIN; let cp = State::critical_point(&saft, None, Some(t), Default::default())?; assert_relative_eq!(cp.temperature, 375.12441 * KELVIN, max_relative = 1e-8); @@ -35,7 +35,7 @@ fn test_critical_point_mix() -> Result<(), Box> { None, IdentifierOption::Name, )?; - let saft = Rc::new(PcSaft::new(Rc::new(params))); + let saft = Arc::new(PcSaft::new(Arc::new(params))); let t = 300.0 * KELVIN; let moles = arr1(&[1.5, 1.5]) * MOL; let cp = State::critical_point(&saft, Some(&moles), Some(t), Default::default())?; diff --git a/tests/pcsaft/dft.rs b/tests/pcsaft/dft.rs index 7ffafed65..45ce81c55 100644 --- a/tests/pcsaft/dft.rs +++ b/tests/pcsaft/dft.rs @@ -9,7 +9,7 @@ use feos_dft::interface::PlanarInterface; use ndarray::{arr1, Axis}; use quantity::si::*; use std::error::Error; -use std::rc::Rc; +use std::sync::Arc; #[test] #[allow(non_snake_case)] @@ -18,15 +18,15 @@ fn test_bulk_implementations() -> Result<(), Box> { let KB_old = 1.38064852e-23 * JOULE / KELVIN; let NAV_old = 6.022140857e23 / MOL; - let params = Rc::new(PcSaftParameters::from_json( + let params = Arc::new(PcSaftParameters::from_json( vec!["water_np"], "tests/pcsaft/test_parameters.json", None, IdentifierOption::Name, )?); - let eos = Rc::new(PcSaft::new(params.clone())); - let func_pure = Rc::new(PcSaftFunctional::new(params.clone())); - let func_full = Rc::new(PcSaftFunctional::new_full( + let eos = Arc::new(PcSaft::new(params.clone())); + let func_pure = Arc::new(PcSaftFunctional::new(params.clone())); + let func_full = Arc::new(PcSaftFunctional::new_full( params, FMTVersion::KierlikRosinberg, )); @@ -93,18 +93,18 @@ fn test_dft_propane() -> Result<(), Box> { let KB_old = 1.38064852e-23 * JOULE / KELVIN; let NAV_old = 6.022140857e23 / MOL; - let params = Rc::new(PcSaftParameters::from_json( + let params = Arc::new(PcSaftParameters::from_json( vec!["propane"], "tests/pcsaft/test_parameters.json", None, IdentifierOption::Name, )?); - let func_pure = Rc::new(PcSaftFunctional::new(params.clone())); - let func_full = Rc::new(PcSaftFunctional::new_full( + let func_pure = Arc::new(PcSaftFunctional::new(params.clone())); + let func_full = Arc::new(PcSaftFunctional::new_full( params.clone(), FMTVersion::KierlikRosinberg, )); - let func_full_vec = Rc::new(PcSaftFunctional::new_full(params, FMTVersion::WhiteBear)); + let func_full_vec = Arc::new(PcSaftFunctional::new_full(params, FMTVersion::WhiteBear)); let t = 200.0 * KELVIN; let w = 150.0 * ANGSTROM; let points = 2048; @@ -217,14 +217,14 @@ fn test_dft_water() -> Result<(), Box> { let KB_old = 1.38064852e-23 * JOULE / KELVIN; let NAV_old = 6.022140857e23 / MOL; - let params = Rc::new(PcSaftParameters::from_json( + let params = Arc::new(PcSaftParameters::from_json( vec!["water_np"], "tests/pcsaft/test_parameters.json", None, IdentifierOption::Name, )?); - let func_pure = Rc::new(PcSaftFunctional::new(params.clone())); - let func_full_vec = Rc::new(PcSaftFunctional::new_full(params, FMTVersion::WhiteBear)); + let func_pure = Arc::new(PcSaftFunctional::new(params.clone())); + let func_full_vec = Arc::new(PcSaftFunctional::new_full(params, FMTVersion::WhiteBear)); let t = 400.0 * KELVIN; let w = 120.0 * ANGSTROM; let points = 2048; @@ -304,7 +304,7 @@ fn test_entropy_bulk_values() -> Result<(), Box> { None, IdentifierOption::Name, )?; - let func = Rc::new(PcSaftFunctional::new(Rc::new(params))); + let func = Arc::new(PcSaftFunctional::new(Arc::new(params))); let vle = PhaseEquilibrium::pure(&func, 350.0 * KELVIN, None, Default::default())?; let profile = PlanarInterface::from_pdgt(&vle, 2048)?.solve(None)?; let s_res = profile diff --git a/tests/pcsaft/properties.rs b/tests/pcsaft/properties.rs index a39605c72..f27c33712 100644 --- a/tests/pcsaft/properties.rs +++ b/tests/pcsaft/properties.rs @@ -5,7 +5,7 @@ use feos_core::{EquationOfState, StateBuilder}; use ndarray::*; use quantity::si::*; use std::error::Error; -use std::rc::Rc; +use std::sync::Arc; #[test] fn test_dln_phi_dp() -> Result<(), Box> { @@ -15,7 +15,7 @@ fn test_dln_phi_dp() -> Result<(), Box> { None, IdentifierOption::Name, )?; - let saft = Rc::new(PcSaft::new(Rc::new(params))); + let saft = Arc::new(PcSaft::new(Arc::new(params))); let t = 300.0 * KELVIN; let p = BAR; let h = 1e-1 * PASCAL; @@ -48,7 +48,7 @@ fn test_virial_is_not_nan() -> Result<(), Box> { None, IdentifierOption::Name, )?; - let saft = Rc::new(PcSaft::new(Rc::new(params))); + let saft = Arc::new(PcSaft::new(Arc::new(params))); let virial_b = saft.second_virial_coefficient(300.0 * KELVIN, None)?; assert!(!virial_b.is_nan()); Ok(()) diff --git a/tests/pcsaft/stability_analysis.rs b/tests/pcsaft/stability_analysis.rs index b6083b4dd..a722f527c 100644 --- a/tests/pcsaft/stability_analysis.rs +++ b/tests/pcsaft/stability_analysis.rs @@ -4,7 +4,7 @@ use feos_core::{DensityInitialization, PhaseEquilibrium, State}; use ndarray::arr1; use quantity::si::*; use std::error::Error; -use std::rc::Rc; +use std::sync::Arc; #[test] fn test_stability_analysis() -> Result<(), Box> { @@ -14,7 +14,7 @@ fn test_stability_analysis() -> Result<(), Box> { None, IdentifierOption::Name, )?; - let mix = Rc::new(PcSaft::new(Rc::new(params))); + let mix = Arc::new(PcSaft::new(Arc::new(params))); let unstable = State::new_npt( &mix, 300.0 * KELVIN, @@ -31,7 +31,7 @@ fn test_stability_analysis() -> Result<(), Box> { None, IdentifierOption::Name, )?; - let mix = Rc::new(PcSaft::new(Rc::new(params))); + let mix = Arc::new(PcSaft::new(Arc::new(params))); let vle = PhaseEquilibrium::bubble_point( &mix, 300.0 * KELVIN, diff --git a/tests/pcsaft/state_creation_mixture.rs b/tests/pcsaft/state_creation_mixture.rs index ee4078196..8b40d8cb6 100644 --- a/tests/pcsaft/state_creation_mixture.rs +++ b/tests/pcsaft/state_creation_mixture.rs @@ -6,10 +6,10 @@ use ndarray::prelude::*; use ndarray::Zip; use quantity::si::*; use std::error::Error; -use std::rc::Rc; +use std::sync::Arc; -fn propane_butane_parameters() -> Result, ParameterError> { - Ok(Rc::new(PcSaftParameters::from_json( +fn propane_butane_parameters() -> Result, ParameterError> { + Ok(Arc::new(PcSaftParameters::from_json( vec!["propane", "butane"], "tests/pcsaft/test_parameters.json", None, @@ -19,7 +19,7 @@ fn propane_butane_parameters() -> Result, ParameterError> { #[test] fn pressure_entropy_molefracs() -> Result<(), Box> { - let saft = Rc::new(PcSaft::new(propane_butane_parameters()?)); + let saft = Arc::new(PcSaft::new(propane_butane_parameters()?)); let pressure = BAR; let temperature = 300.0 * KELVIN; let x = arr1(&[0.3, 0.7]); @@ -50,7 +50,7 @@ fn pressure_entropy_molefracs() -> Result<(), Box> { #[test] fn volume_temperature_molefracs() -> Result<(), Box> { - let saft = Rc::new(PcSaft::new(propane_butane_parameters()?)); + let saft = Arc::new(PcSaft::new(propane_butane_parameters()?)); let temperature = 300.0 * KELVIN; let volume = 1.5e-3 * METER.powi(3); let moles = MOL; @@ -67,7 +67,7 @@ fn volume_temperature_molefracs() -> Result<(), Box> { #[test] fn temperature_partial_density() -> Result<(), Box> { - let saft = Rc::new(PcSaft::new(propane_butane_parameters()?)); + let saft = Arc::new(PcSaft::new(propane_butane_parameters()?)); let temperature = 300.0 * KELVIN; let x = arr1(&[0.3, 0.7]); let partial_density = x.clone() * MOL / METER.powi(3); @@ -86,7 +86,7 @@ fn temperature_partial_density() -> Result<(), Box> { #[test] fn temperature_density_molefracs() -> Result<(), Box> { - let saft = Rc::new(PcSaft::new(propane_butane_parameters()?)); + let saft = Arc::new(PcSaft::new(propane_butane_parameters()?)); let temperature = 300.0 * KELVIN; let x = arr1(&[0.3, 0.7]); let density = MOL / METER.powi(3); @@ -104,7 +104,7 @@ fn temperature_density_molefracs() -> Result<(), Box> { #[test] fn temperature_pressure_molefracs() -> Result<(), Box> { - let saft = Rc::new(PcSaft::new(propane_butane_parameters()?)); + let saft = Arc::new(PcSaft::new(propane_butane_parameters()?)); let temperature = 300.0 * KELVIN; let pressure = BAR; let x = arr1(&[0.3, 0.7]); diff --git a/tests/pcsaft/state_creation_pure.rs b/tests/pcsaft/state_creation_pure.rs index 750c3dd58..a4a8015e5 100644 --- a/tests/pcsaft/state_creation_pure.rs +++ b/tests/pcsaft/state_creation_pure.rs @@ -8,10 +8,10 @@ use feos_core::{ use quantity::si::*; use quantity::QuantityScalar; use std::error::Error; -use std::rc::Rc; +use std::sync::Arc; -fn propane_parameters() -> Result, ParameterError> { - Ok(Rc::new(PcSaftParameters::from_json( +fn propane_parameters() -> Result, ParameterError> { + Ok(Arc::new(PcSaftParameters::from_json( vec!["propane"], "tests/pcsaft/test_parameters.json", None, @@ -21,7 +21,7 @@ fn propane_parameters() -> Result, ParameterError> { #[test] fn temperature_volume() -> Result<(), Box> { - let saft = Rc::new(PcSaft::new(propane_parameters()?)); + let saft = Arc::new(PcSaft::new(propane_parameters()?)); let temperature = 300.0 * KELVIN; let volume = 1.5e-3 * METER.powi(3); let moles = MOL; @@ -36,7 +36,7 @@ fn temperature_volume() -> Result<(), Box> { #[test] fn temperature_density() -> Result<(), Box> { - let saft = Rc::new(PcSaft::new(propane_parameters()?)); + let saft = Arc::new(PcSaft::new(propane_parameters()?)); let temperature = 300.0 * KELVIN; let density = MOL / METER.powi(3); let state = StateBuilder::new(&saft) @@ -49,7 +49,7 @@ fn temperature_density() -> Result<(), Box> { #[test] fn temperature_total_moles_volume() -> Result<(), Box> { - let saft = Rc::new(PcSaft::new(propane_parameters()?)); + let saft = Arc::new(PcSaft::new(propane_parameters()?)); let temperature = 300.0 * KELVIN; let total_moles = MOL; let volume = METER.powi(3); @@ -65,7 +65,7 @@ fn temperature_total_moles_volume() -> Result<(), Box> { #[test] fn temperature_total_moles_density() -> Result<(), Box> { - let saft = Rc::new(PcSaft::new(propane_parameters()?)); + let saft = Arc::new(PcSaft::new(propane_parameters()?)); let temperature = 300.0 * KELVIN; let total_moles = MOL; let density = MOL / METER.powi(3); @@ -84,7 +84,7 @@ fn temperature_total_moles_density() -> Result<(), Box> { #[test] fn pressure_temperature() -> Result<(), Box> { - let saft = Rc::new(PcSaft::new(propane_parameters()?)); + let saft = Arc::new(PcSaft::new(propane_parameters()?)); let pressure = BAR; let temperature = 300.0 * KELVIN; let state = StateBuilder::new(&saft) @@ -101,7 +101,7 @@ fn pressure_temperature() -> Result<(), Box> { #[test] fn pressure_temperature_phase() -> Result<(), Box> { - let saft = Rc::new(PcSaft::new(propane_parameters()?)); + let saft = Arc::new(PcSaft::new(propane_parameters()?)); let pressure = BAR; let temperature = 300.0 * KELVIN; let state = StateBuilder::new(&saft) @@ -119,7 +119,7 @@ fn pressure_temperature_phase() -> Result<(), Box> { #[test] fn pressure_temperature_initial_density() -> Result<(), Box> { - let saft = Rc::new(PcSaft::new(propane_parameters()?)); + let saft = Arc::new(PcSaft::new(propane_parameters()?)); let pressure = BAR; let temperature = 300.0 * KELVIN; let state = StateBuilder::new(&saft) @@ -137,7 +137,7 @@ fn pressure_temperature_initial_density() -> Result<(), Box> { #[test] fn pressure_enthalpy_vapor() -> Result<(), Box> { - let saft = Rc::new(PcSaft::new(propane_parameters()?)); + let saft = Arc::new(PcSaft::new(propane_parameters()?)); let pressure = 0.3 * BAR; let molar_enthalpy = 2000.0 * JOULE / MOL; let state = StateBuilder::new(&saft) @@ -176,7 +176,7 @@ fn pressure_enthalpy_vapor() -> Result<(), Box> { #[test] fn density_internal_energy() -> Result<(), Box> { - let saft = Rc::new(PcSaft::new(propane_parameters()?)); + let saft = Arc::new(PcSaft::new(propane_parameters()?)); let pressure = 5.0 * BAR; let temperature = 315.0 * KELVIN; let total_moles = 2.5 * MOL; @@ -203,7 +203,7 @@ fn density_internal_energy() -> Result<(), Box> { #[test] fn pressure_enthalpy_total_moles_vapor() -> Result<(), Box> { - let saft = Rc::new(PcSaft::new(propane_parameters()?)); + let saft = Arc::new(PcSaft::new(propane_parameters()?)); let pressure = 0.3 * BAR; let molar_enthalpy = 2000.0 * JOULE / MOL; let total_moles = 2.5 * MOL; @@ -244,7 +244,7 @@ fn pressure_enthalpy_total_moles_vapor() -> Result<(), Box> { #[test] fn pressure_entropy_vapor() -> Result<(), Box> { - let saft = Rc::new(PcSaft::new(propane_parameters()?)); + let saft = Arc::new(PcSaft::new(propane_parameters()?)); let pressure = 0.3 * BAR; let molar_entropy = -2.0 * JOULE / MOL / KELVIN; let state = StateBuilder::new(&saft) @@ -283,7 +283,7 @@ fn pressure_entropy_vapor() -> Result<(), Box> { #[test] fn temperature_entropy_vapor() -> Result<(), Box> { - let saft = Rc::new(PcSaft::new(propane_parameters()?)); + let saft = Arc::new(PcSaft::new(propane_parameters()?)); let pressure = 3.0 * BAR; let temperature = 315.15 * KELVIN; let total_moles = 3.0 * MOL; @@ -341,7 +341,7 @@ fn assert_multiple_states( #[test] fn test_consistency() -> Result<(), Box> { let p = propane_parameters()?; - let saft = Rc::new(PcSaft::new(p)); + let saft = Arc::new(PcSaft::new(p)); let temperatures = [350.0 * KELVIN, 400.0 * KELVIN, 450.0 * KELVIN]; let pressures = [1.0 * BAR, 2.0 * BAR, 3.0 * BAR]; diff --git a/tests/pcsaft/tp_flash.rs b/tests/pcsaft/tp_flash.rs index 41b27e155..ff57226ea 100644 --- a/tests/pcsaft/tp_flash.rs +++ b/tests/pcsaft/tp_flash.rs @@ -5,10 +5,10 @@ use feos_core::{Contributions, PhaseEquilibrium, SolverOptions}; use ndarray::*; use quantity::si::*; use std::error::Error; -use std::rc::Rc; +use std::sync::Arc; -fn read_params(components: Vec<&str>) -> Result, ParameterError> { - Ok(Rc::new(PcSaftParameters::from_json( +fn read_params(components: Vec<&str>) -> Result, ParameterError> { + Ok(Arc::new(PcSaftParameters::from_json( components, "tests/pcsaft/test_parameters.json", None, @@ -18,8 +18,8 @@ fn read_params(components: Vec<&str>) -> Result, ParameterE #[test] fn test_tp_flash() -> Result<(), Box> { - let propane = Rc::new(PcSaft::new(read_params(vec!["propane"])?)); - let butane = Rc::new(PcSaft::new(read_params(vec!["butane"])?)); + let propane = Arc::new(PcSaft::new(read_params(vec!["propane"])?)); + let butane = Arc::new(PcSaft::new(read_params(vec!["butane"])?)); let t = 250.0 * KELVIN; let p_propane = PhaseEquilibrium::pure(&propane, t, None, Default::default())? .vapor() @@ -32,7 +32,7 @@ fn test_tp_flash() -> Result<(), Box> { let y1 = (x1 * p_propane / p).into_value()?; let z1 = 0.5 * (x1 + y1); println!("{} {} {} {} {}", p_propane, p_butane, x1, y1, z1); - let mix = Rc::new(PcSaft::new(read_params(vec!["propane", "butane"])?)); + let mix = Arc::new(PcSaft::new(read_params(vec!["propane", "butane"])?)); let options = SolverOptions::new().max_iter(100).tol(1e-12); let vle = PhaseEquilibrium::tp_flash( &mix, diff --git a/tests/pcsaft/vle_pure.rs b/tests/pcsaft/vle_pure.rs index 19381bb49..eb31d5078 100644 --- a/tests/pcsaft/vle_pure.rs +++ b/tests/pcsaft/vle_pure.rs @@ -4,7 +4,7 @@ use feos_core::parameter::{IdentifierOption, Parameter}; use feos_core::{Contributions, PhaseEquilibrium}; use quantity::si::*; use std::error::Error; -use std::rc::Rc; +use std::sync::Arc; #[test] fn vle_pure_temperature() -> Result<(), Box> { @@ -14,7 +14,7 @@ fn vle_pure_temperature() -> Result<(), Box> { None, IdentifierOption::Name, )?; - let saft = Rc::new(PcSaft::new(Rc::new(params))); + let saft = Arc::new(PcSaft::new(Arc::new(params))); let temperatures = [ 170.0 * KELVIN, 200.0 * KELVIN, @@ -42,7 +42,7 @@ fn vle_pure_pressure() -> Result<(), Box> { None, IdentifierOption::Name, )?; - let saft = Rc::new(PcSaft::new(Rc::new(params))); + let saft = Arc::new(PcSaft::new(Arc::new(params))); let pressures = [0.1 * BAR, 1.0 * BAR, 10.0 * BAR, 30.0 * BAR, 44.0 * BAR]; for &p in pressures.iter() { let state = PhaseEquilibrium::pure(&saft, p, None, Default::default())?;