diff --git a/Cargo.toml b/Cargo.toml index 32771d66a..d45aa8b22 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -44,6 +44,7 @@ criterion = "0.8" paste = "1.0" rusqlite = "0.40" csv = "1.0" +either = "1.0" feos-core = { version = "0.10", path = "crates/feos-core" } feos-dft = { version = "0.10", path = "crates/feos-dft" } diff --git a/crates/feos-core/Cargo.toml b/crates/feos-core/Cargo.toml index c1831a1b7..8ee3138a6 100644 --- a/crates/feos-core/Cargo.toml +++ b/crates/feos-core/Cargo.toml @@ -28,6 +28,7 @@ rayon = { workspace = true, optional = true } csv = { workspace = true } itertools = { workspace = true } rusqlite = { workspace = true, features = ["bundled"], optional = true } +either = { workspace = true } [dev-dependencies] approx = { workspace = true } diff --git a/crates/feos-core/src/phase_equilibria/bubble_dew.rs b/crates/feos-core/src/phase_equilibria/bubble_dew.rs index f956d4d7d..b014c4e03 100644 --- a/crates/feos-core/src/phase_equilibria/bubble_dew.rs +++ b/crates/feos-core/src/phase_equilibria/bubble_dew.rs @@ -5,6 +5,7 @@ use crate::state::{ DensityInitialization::{InitialDensity, Liquid, Vapor}, }; use crate::{Composition, ReferenceSystem, Residual, SolverOptions, State, Verbosity}; +use either::Either; use nalgebra::allocator::Allocator; use nalgebra::{DMatrix, DVector, DefaultAllocator, Dim, Dyn, OVector, U1}; #[cfg(feature = "ndarray")] @@ -31,10 +32,11 @@ pub trait TemperatureOrPressure + Copy = f64>: Copy { fn temperature(&self) -> Option>; fn pressure(&self) -> Option>; + #[expect(clippy::type_complexity)] fn temperature_pressure( &self, tp_init: Option, - ) -> (Option>, Option>, bool); + ) -> Either<(Temperature, Option), (Pressure, Option)>; fn from_state, N: Gradients>(state: &State) -> Self::Other where @@ -50,7 +52,7 @@ pub trait TemperatureOrPressure + Copy = f64>: Copy { } impl + Copy> TemperatureOrPressure for Temperature { - type Other = Pressure; + type Other = Pressure; const IDENTIFIER: &'static str = "temperature"; fn temperature(&self) -> Option> { @@ -64,27 +66,27 @@ impl + Copy> TemperatureOrPressure for Temperature { fn temperature_pressure( &self, tp_init: Option, - ) -> (Option>, Option>, bool) { - (Some(*self), tp_init, true) + ) -> Either<(Temperature, Option), (Pressure, Option)> { + Either::Left((*self, tp_init)) } fn from_state, N: Gradients>(state: &State) -> Self::Other where DefaultAllocator: Allocator, { - state.pressure(Contributions::Total) + state.pressure(Contributions::Total).re() } #[cfg(feature = "ndarray")] fn linspace( &self, - start: Pressure, - end: Pressure, + start: Pressure, + end: Pressure, n: usize, ) -> (Temperature>, Pressure>) { ( Temperature::linspace(self.re(), self.re(), n), - Pressure::linspace(start.re(), end.re(), n), + Pressure::linspace(start, end, n), ) } } @@ -95,7 +97,7 @@ impl + Copy> TemperatureOrPressure for Temperature { impl + Copy> TemperatureOrPressure for Quantity> { - type Other = Temperature; + type Other = Temperature; const IDENTIFIER: &'static str = "pressure"; fn temperature(&self) -> Option> { @@ -109,26 +111,26 @@ impl + Copy> TemperatureOrPressure fn temperature_pressure( &self, tp_init: Option, - ) -> (Option>, Option>, bool) { - (tp_init, Some(*self), false) + ) -> Either<(Temperature, Option), (Pressure, Option)> { + Either::Right((*self, tp_init)) } fn from_state, N: Dim>(state: &State) -> Self::Other where DefaultAllocator: Allocator, { - state.temperature + state.temperature.re() } #[cfg(feature = "ndarray")] fn linspace( &self, - start: Temperature, - end: Temperature, + start: Temperature, + end: Temperature, n: usize, ) -> (Temperature>, Pressure>) { ( - Temperature::linspace(start.re(), end.re(), n), + Temperature::linspace(start, end, n), Pressure::linspace(self.re(), self.re(), n), ) } @@ -197,134 +199,113 @@ where } Ok(vle) } else { - let (temperature, pressure, iterate_p) = - temperature_or_pressure.temperature_pressure(tp_init); Self::bubble_dew_point_tp( eos, - temperature, - pressure, + temperature_or_pressure, + tp_init, vapor_molefracs, liquid_molefracs, bubble, - iterate_p, options, ) } } - #[expect(clippy::too_many_arguments)] - fn bubble_dew_point_tp>( + fn bubble_dew_point_tp, X: Composition>( eos: &E, - temperature: Option>, - pressure: Option>, + temperature_or_pressure: TP, + tp_init: Option, composition: X, molefracs_init: Option<&OVector>, bubble: bool, - iterate_p: bool, options: (SolverOptions, SolverOptions), ) -> FeosResult { let eos_re = eos.re(); - let mut temperature_re = temperature.map(|t| t.re()); - let mut pressure_re = pressure.map(|p| p.re()); + let iterate_p = temperature_or_pressure + .temperature_pressure(tp_init) + .is_left(); let (molefracs_spec, total_moles) = composition.into_molefracs(eos)?; let molefracs_spec_re = molefracs_spec.map(|x| x.re()); - let (v1, rho2) = if iterate_p { - // temperature is specified - let temperature_re = temperature_re.as_mut().ok_or(FeosError::Error( - "Temperature information is expected for bubble/dew calculation.".to_string(), - ))?; - - // First use given initial pressure if applicable - if let Some(p) = pressure_re.as_mut() { - PhaseEquilibrium::iterate_bubble_dew( - &eos_re, - temperature_re, - p, - &molefracs_spec_re, - molefracs_init, - bubble, - iterate_p, - options, - )? - } else { - // Next try to initialize with an ideal gas assumption - let x2 = PhaseEquilibrium::starting_pressure_ideal_gas( - &eos_re, - *temperature_re, - &molefracs_spec_re, - bubble, - ) - .and_then(|(p, x)| { - let p = pressure_re.insert(p); - PhaseEquilibrium::iterate_bubble_dew( + let (mut t, mut p, v1, rho2) = match temperature_or_pressure.temperature_pressure(tp_init) { + Either::Left((t, mut p)) => { + // First use given initial pressure if applicable + let (p, v1, rho2) = if let Some(p) = p.as_mut() { + let (v1, rho2) = PhaseEquilibrium::iterate_bubble_dew( &eos_re, - temperature_re, + &mut t.re(), p, &molefracs_spec_re, - molefracs_init.or(Some(&x)), + molefracs_init, bubble, iterate_p, options, - ) - }); - - // Finally use the spinodal to initialize the calculation - x2.or_else(|_| { - PhaseEquilibrium::starting_pressure_spinodal( + )?; + (*p, v1, rho2) + } else { + let x2 = PhaseEquilibrium::starting_pressure_ideal_gas( &eos_re, - *temperature_re, + t.re(), &molefracs_spec_re, + bubble, ) - .and_then(|p| { - let p = pressure_re.insert(p); - PhaseEquilibrium::iterate_bubble_dew( + .and_then(|(mut p, x)| { + let (v1, rho2) = PhaseEquilibrium::iterate_bubble_dew( &eos_re, - temperature_re, - p, + &mut t.re(), + &mut p, &molefracs_spec_re, - molefracs_init, + molefracs_init.or(Some(&x)), bubble, iterate_p, options, + )?; + Ok((p, v1, rho2)) + }); + + // Finally use the spinodal to initialize the calculation + x2.or_else(|_| { + PhaseEquilibrium::starting_pressure_spinodal( + &eos_re, + t.re(), + &molefracs_spec_re, ) - }) - })? + .and_then(|mut p| { + let (v1, rho2) = PhaseEquilibrium::iterate_bubble_dew( + &eos_re, + &mut t.re(), + &mut p, + &molefracs_spec_re, + molefracs_init, + bubble, + iterate_p, + options, + )?; + Ok((p, v1, rho2)) + }) + })? + }; + (t.into_reduced(), D::from(p.into_reduced()), v1, rho2) } - } else { - // pressure is specified - let pressure_re = pressure_re.as_mut().ok_or(FeosError::Error( - "Pressure information is expected for bubble/dew calculation.".to_string(), - ))?; - - let temperature_re = temperature_re - .as_mut() + Either::Right((p, mut t)) => { + let mut pressure_re = p.re(); + let t = t.as_mut() .ok_or(FeosError::Error( "An initial temperature is required for the calculation of bubble/dew points at given pressure.".to_string()))?; - PhaseEquilibrium::iterate_bubble_dew( - &eos.re(), - temperature_re, - pressure_re, - &molefracs_spec_re, - molefracs_init, - bubble, - iterate_p, - options, - )? + let (v1, rho2) = PhaseEquilibrium::iterate_bubble_dew( + &eos.re(), + t, + &mut pressure_re, + &molefracs_spec_re, + molefracs_init, + bubble, + iterate_p, + options, + )?; + (D::from(t.into_reduced()), p.into_reduced(), v1, rho2) + } }; // implicit differentiation - // unwraps here are safe - let (mut t, mut p) = if iterate_p { - ( - temperature.unwrap().into_reduced(), - D::from(pressure_re.unwrap().into_reduced()), - ) - } else { - ( - D::from(temperature_re.unwrap().into_reduced()), - pressure.unwrap().into_reduced(), - ) - }; let mut molar_volume = D::from(v1); let mut rho2 = rho2.map(D::from); for _ in 0..D::NDERIV { diff --git a/crates/feos-core/src/phase_equilibria/phase_diagram_binary.rs b/crates/feos-core/src/phase_equilibria/phase_diagram_binary.rs index a3b23cbcc..65b0cbf46 100644 --- a/crates/feos-core/src/phase_equilibria/phase_diagram_binary.rs +++ b/crates/feos-core/src/phase_equilibria/phase_diagram_binary.rs @@ -3,6 +3,7 @@ use super::{PhaseDiagram, PhaseEquilibrium}; use crate::errors::{FeosError, FeosResult}; use crate::state::{Contributions, DensityInitialization::Vapor, State}; use crate::{ReferenceSystem, Residual, SolverOptions, Subset}; +use either::Either; use nalgebra::{DVector, dvector, matrix, stack, vector}; use ndarray::{Array1, s}; use num_dual::linalg::LU; @@ -503,31 +504,23 @@ impl PhaseEquilibrium { options: SolverOptions, bubble_dew_options: (SolverOptions, SolverOptions), ) -> FeosResult { - let (temperature, pressure, iterate_p) = - temperature_or_pressure.temperature_pressure(tp_init); - if iterate_p { - PhaseEquilibrium::heteroazeotrope_t( + match temperature_or_pressure.temperature_pressure(tp_init) { + Either::Left((temperature, pressure)) => PhaseEquilibrium::heteroazeotrope_t( eos, - temperature.ok_or(FeosError::Error( - "Temperature information is expected for heteroazeotrope calculation." - .to_string(), - ))?, + temperature, x_init, pressure, options, bubble_dew_options, - ) - } else { - PhaseEquilibrium::heteroazeotrope_p( + ), + Either::Right((pressure, temperature)) => PhaseEquilibrium::heteroazeotrope_p( eos, - pressure.ok_or(FeosError::Error( - "Pressure information is expected for heteroazeotrope calculation.".to_string(), - ))?, + pressure, x_init, temperature, options, bubble_dew_options, - ) + ), } } @@ -826,12 +819,14 @@ impl PhaseEquilibrium { let x0 = -ln_alpha1 / (ln_alpha2 - ln_alpha1); // solve for the azeotropic composition and return the corresponding VLE state - let (temperature, pressure, iterate_t) = temperature_or_pressure.temperature_pressure(None); - (if iterate_t { - Self::iterate_azeotrope_t(eos, temperature.unwrap(), x0, 10, 1e-10) - } else { - let t_init = vle1.liquid().temperature.min(vle2.liquid().temperature); - Self::iterate_azeotrope_p(eos, pressure.unwrap(), x0, t_init, 10, 1e-10) + (match temperature_or_pressure.temperature_pressure(None) { + Either::Left((temperature, _)) => { + Self::iterate_azeotrope_t(eos, temperature, x0, 10, 1e-10) + } + Either::Right((pressure, _)) => { + let t_init = vle1.liquid().temperature.min(vle2.liquid().temperature); + Self::iterate_azeotrope_p(eos, pressure, x0, t_init, 10, 1e-10) + } }) .map(Some) } @@ -849,7 +844,7 @@ impl PhaseEquilibrium { PhaseEquilibrium::bubble_point( &eos.lift(), t, - &dvector![x, -x + 1.0], + dvector![x, -x + 1.0], None, None, Default::default(), @@ -884,7 +879,7 @@ impl PhaseEquilibrium { PhaseEquilibrium::bubble_point( &eos.lift(), p, - &dvector![x, -x + 1.0], + dvector![x, -x + 1.0], Some(t_init), None, Default::default(), diff --git a/crates/feos/src/pcsaft/eos/mod.rs b/crates/feos/src/pcsaft/eos/mod.rs index 6a68f67f7..aa8e05d0c 100644 --- a/crates/feos/src/pcsaft/eos/mod.rs +++ b/crates/feos/src/pcsaft/eos/mod.rs @@ -890,7 +890,7 @@ mod tests_parameter_fit { let pcsaft_ad = PcSaftBinary::::seed_derivatives(&flat_binary_params(&pcsaft), ["k_ij"]); let pressure = Pressure::from_reduced(DualVec::from(45. * BAR.into_reduced())); - let t_init = Temperature::from_reduced(DualVec::from(500.0)); + let t_init = Temperature::from_reduced(500.0); let x = DualVec::from(0.5); let t = PhaseEquilibrium::bubble_point( &pcsaft_ad, @@ -939,7 +939,7 @@ mod tests_parameter_fit { let pcsaft_ad = PcSaftBinary::::seed_derivatives(&flat_binary_params(&pcsaft), ["k_ij"]); let pressure = Pressure::from_reduced(DualVec::from(45. * BAR.into_reduced())); - let t_init = Temperature::from_reduced(DualVec::from(500.0)); + let t_init = Temperature::from_reduced(500.0); let x = DualVec::from(0.5); let t = PhaseEquilibrium::dew_point( &pcsaft_ad,