From 118789e0b5328204d516ff17eacd64664d987f5c Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Tue, 21 Jul 2026 17:58:24 +0200 Subject: [PATCH] Clean the implementation of temperature_or_pressure arguments somewhat --- CHANGELOG.md | 4 + .../src/phase_equilibria/bubble_dew.rs | 206 ++++++++---------- .../phase_equilibria/phase_diagram_binary.rs | 65 +++--- crates/feos/src/pcsaft/eos/mod.rs | 4 +- 4 files changed, 132 insertions(+), 147 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7a05e9188..adc0b1945 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,10 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Breaking] +### Changed +- Changed data type of initial temperatures or pressure for phase equilibrium calculations (`TemperatureOrPressure::Other`) from `D` to `f64`. [#369](https://github.com/feos-org/feos/pull/369) + ## [Unreleased] ## [0.10.0] - 2026-07-15 diff --git a/crates/feos-core/src/phase_equilibria/bubble_dew.rs b/crates/feos-core/src/phase_equilibria/bubble_dew.rs index f956d4d7d..6d2c9333d 100644 --- a/crates/feos-core/src/phase_equilibria/bubble_dew.rs +++ b/crates/feos-core/src/phase_equilibria/bubble_dew.rs @@ -31,10 +31,7 @@ pub trait TemperatureOrPressure + Copy = f64>: Copy { fn temperature(&self) -> Option>; fn pressure(&self) -> Option>; - fn temperature_pressure( - &self, - tp_init: Option, - ) -> (Option>, Option>, bool); + fn specification(&self, tp_init: Option) -> TemperatureOrPressureSpecification; fn from_state, N: Gradients>(state: &State) -> Self::Other where @@ -50,7 +47,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> { @@ -61,30 +58,27 @@ impl + Copy> TemperatureOrPressure for Temperature { None } - fn temperature_pressure( - &self, - tp_init: Option, - ) -> (Option>, Option>, bool) { - (Some(*self), tp_init, true) + fn specification(&self, tp_init: Option) -> TemperatureOrPressureSpecification { + TemperatureOrPressureSpecification::Temperature(*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 +89,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> { @@ -106,34 +100,43 @@ impl + Copy> TemperatureOrPressure Some(*self) } - fn temperature_pressure( - &self, - tp_init: Option, - ) -> (Option>, Option>, bool) { - (tp_init, Some(*self), false) + fn specification(&self, tp_init: Option) -> TemperatureOrPressureSpecification { + TemperatureOrPressureSpecification::Pressure(*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), ) } } +/// Specification and initial value for a phase equilibrium calculation. +pub enum TemperatureOrPressureSpecification { + Temperature(Temperature, Option), + Pressure(Pressure, Option), +} + +impl TemperatureOrPressureSpecification { + fn is_temperature(&self) -> bool { + matches!(self, Self::Temperature(_, _)) + } +} + /// # Bubble and dew point calculations impl, N: Gradients, D: DualNum + Copy> PhaseEquilibrium where @@ -197,134 +200,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 + .specification(tp_init) + .is_temperature(); 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.specification(tp_init) { + TemperatureOrPressureSpecification::Temperature(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() + TemperatureOrPressureSpecification::Pressure(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..f33b6b75d 100644 --- a/crates/feos-core/src/phase_equilibria/phase_diagram_binary.rs +++ b/crates/feos-core/src/phase_equilibria/phase_diagram_binary.rs @@ -1,6 +1,7 @@ use super::bubble_dew::TemperatureOrPressure; use super::{PhaseDiagram, PhaseEquilibrium}; use crate::errors::{FeosError, FeosResult}; +use crate::phase_equilibria::bubble_dew::TemperatureOrPressureSpecification; use crate::state::{Contributions, DensityInitialization::Vapor, State}; use crate::{ReferenceSystem, Residual, SolverOptions, Subset}; use nalgebra::{DVector, dvector, matrix, stack, vector}; @@ -503,31 +504,27 @@ 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( - eos, - temperature.ok_or(FeosError::Error( - "Temperature information is expected for heteroazeotrope calculation." - .to_string(), - ))?, - x_init, - pressure, - options, - bubble_dew_options, - ) - } else { - PhaseEquilibrium::heteroazeotrope_p( - eos, - pressure.ok_or(FeosError::Error( - "Pressure information is expected for heteroazeotrope calculation.".to_string(), - ))?, - x_init, - temperature, - options, - bubble_dew_options, - ) + match temperature_or_pressure.specification(tp_init) { + TemperatureOrPressureSpecification::Temperature(temperature, pressure) => { + PhaseEquilibrium::heteroazeotrope_t( + eos, + temperature, + x_init, + pressure, + options, + bubble_dew_options, + ) + } + TemperatureOrPressureSpecification::Pressure(pressure, temperature) => { + PhaseEquilibrium::heteroazeotrope_p( + eos, + pressure, + x_init, + temperature, + options, + bubble_dew_options, + ) + } } } @@ -826,12 +823,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.specification(None) { + TemperatureOrPressureSpecification::Temperature(temperature, _) => { + Self::iterate_azeotrope_t(eos, temperature, x0, 10, 1e-10) + } + TemperatureOrPressureSpecification::Pressure(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 +848,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 +883,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,