From 55a40ef2ca0db2dbb37e03d542319e7669ea789c Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Tue, 27 Jan 2026 15:57:59 +0100 Subject: [PATCH] Enable AD for pure-component VLEs with given pressure --- crates/feos-core/src/ad/mod.rs | 66 ++- .../src/equation_of_state/residual.rs | 36 +- crates/feos-core/src/phase_equilibria/mod.rs | 16 +- .../src/phase_equilibria/vle_pure.rs | 432 ++++++++++-------- .../src/state/residual_properties.rs | 9 - crates/feos/src/pcsaft/eos/mod.rs | 44 ++ py-feos/src/ad/mod.rs | 28 +- py-feos/src/lib.rs | 1 + 8 files changed, 411 insertions(+), 221 deletions(-) diff --git a/crates/feos-core/src/ad/mod.rs b/crates/feos-core/src/ad/mod.rs index ce92d05d3..308e722d1 100644 --- a/crates/feos-core/src/ad/mod.rs +++ b/crates/feos-core/src/ad/mod.rs @@ -4,7 +4,7 @@ use crate::{FeosResult, PhaseEquilibrium, ReferenceSystem, Residual}; use nalgebra::{Const, SVector, U1, U2}; #[cfg(feature = "rayon")] use ndarray::{Array1, Array2, ArrayView2, Zip}; -use num_dual::{Derivative, DualSVec, DualStruct}; +use num_dual::{Derivative, DualNum, DualSVec, DualStruct, first_derivative, partial2}; use quantity::{Density, Pressure, Temperature}; #[cfg(feature = "rayon")] use quantity::{KELVIN, KILO, METER, MOL, PASCAL}; @@ -66,6 +66,50 @@ pub trait PropertiesAD { Ok(Pressure::from_reduced(p)) } + fn boiling_temperature( + &self, + pressure: Pressure, + ) -> FeosResult>> + where + Self: Residual>, + { + let eos_f64 = self.re(); + let (temperature, [vapor_density, liquid_density]) = + PhaseEquilibrium::pure_p(&eos_f64, pressure, None, Default::default())?; + + // implicit differentiation is implemented here instead of just calling pure_t with dual + // numbers, because for the first derivative, we can avoid calculating density derivatives. + let t = temperature.into_reduced(); + let v1 = 1.0 / liquid_density.to_reduced(); + let v2 = 1.0 / vapor_density.to_reduced(); + let p = pressure.into_reduced(); + let t = Gradient::from(t); + let t = t + { + let v1 = Gradient::from(v1); + let v2 = Gradient::from(v2); + let p = Gradient::from(p); + let x = Self::pure_molefracs(); + + let residual_entropy = |v| { + let (a, s) = first_derivative( + partial2( + |t, &v, x| self.lift().residual_molar_helmholtz_energy(t, v, x), + &v, + &x, + ), + t, + ); + (a, -s) + }; + let (a1, s1) = residual_entropy(v1); + let (a2, s2) = residual_entropy(v2); + + let ln_rho = (v1 / v2).ln(); + (p * (v2 - v1) + (a2 - a1 + t * ln_rho)) / (s2 - s1 - ln_rho) + }; + Ok(Temperature::from_reduced(t)) + } + fn equilibrium_liquid_density( &self, temperature: Temperature, @@ -111,6 +155,26 @@ pub trait PropertiesAD { ) } + #[cfg(feature = "rayon")] + fn boiling_temperature_parallel( + parameter_names: [String; P], + parameters: ArrayView2, + input: ArrayView2, + ) -> (Array1, Array2, Array1) + where + Self: ParametersAD<1>, + { + parallelize::<_, Self, _, _>( + parameter_names, + parameters, + input, + |eos: &Self::Lifted>, inp| { + eos.boiling_temperature(inp[0] * PASCAL) + .map(|p| p.convert_into(KELVIN)) + }, + ) + } + #[cfg(feature = "rayon")] fn liquid_density_parallel( parameter_names: [String; P], diff --git a/crates/feos-core/src/equation_of_state/residual.rs b/crates/feos-core/src/equation_of_state/residual.rs index 326aff883..9de69d3db 100644 --- a/crates/feos-core/src/equation_of_state/residual.rs +++ b/crates/feos-core/src/equation_of_state/residual.rs @@ -1,6 +1,9 @@ use crate::{FeosError, FeosResult, ReferenceSystem, state::StateHD}; +use nalgebra::SVector; use nalgebra::{DVector, DefaultAllocator, Dim, Dyn, OMatrix, OVector, U1, allocator::Allocator}; -use num_dual::{DualNum, Gradients, partial, partial2, second_derivative, third_derivative}; +use num_dual::{ + DualNum, Gradients, hessian, partial, partial2, second_derivative, third_derivative, +}; use quantity::ad::first_derivative; use quantity::*; use std::ops::{Deref, Div}; @@ -313,12 +316,41 @@ where molar_volume, ); ( - a * density, + a, -da + temperature * density, molar_volume * molar_volume * d2a + temperature, ) } + /// calculates a_res, p, s_res, dp_drho, dp_dt + fn p_dpdrho_dpdt( + &self, + temperature: D, + density: D, + molefracs: &OVector, + ) -> (D, D, D, D, D) { + let molar_volume = density.recip(); + let (a, da, d2a) = hessian::<_, _, _, nalgebra::U2, _>( + partial( + |vt: SVector<_, 2>, x: &OVector<_, N>| { + let [[v, t]] = vt.data.0; + self.lift().residual_molar_helmholtz_energy(t, v, x) + }, + molefracs, + ), + &SVector::from([molar_volume, temperature]), + ); + let [[da_dv, da_dt]] = da.data.0; + let [[d2a_dv2, d2a_dvdt], _] = d2a.data.0; + ( + a, + -da_dv + temperature * density, + -da_dt, + molar_volume * molar_volume * d2a_dv2 + temperature, + -d2a_dvdt + density, + ) + } + /// calculates p, dp_drho, d2p_drho2 fn p_dpdrho_d2pdrho2( &self, diff --git a/crates/feos-core/src/phase_equilibria/mod.rs b/crates/feos-core/src/phase_equilibria/mod.rs index 9683e49a1..4a10aa0fe 100644 --- a/crates/feos-core/src/phase_equilibria/mod.rs +++ b/crates/feos-core/src/phase_equilibria/mod.rs @@ -1,5 +1,5 @@ use crate::equation_of_state::Residual; -use crate::errors::{FeosError, FeosResult}; +use crate::errors::FeosResult; use crate::state::{DensityInitialization, State}; use crate::{Contributions, ReferenceSystem}; use nalgebra::allocator::Allocator; @@ -44,12 +44,6 @@ pub struct PhaseEquilibrium + C where DefaultAllocator: Allocator; -// impl Clone for PhaseEquilibrium { -// fn clone(&self) -> Self { -// Self(self.0.clone()) -// } -// } - impl fmt::Display for PhaseEquilibrium { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { for (i, s) in self.0.iter().enumerate() { @@ -224,14 +218,6 @@ impl, N: Dim> PhaseEquilibrium where DefaultAllocator: Allocator, { - pub(super) fn check_trivial_solution(self) -> FeosResult { - if Self::is_trivial_solution(self.vapor(), self.liquid()) { - Err(FeosError::TrivialSolution) - } else { - Ok(self) - } - } - /// Check if the two states form a trivial solution pub fn is_trivial_solution(state1: &State, state2: &State) -> bool { let rho1 = state1.molefracs.clone() * state1.density.into_reduced(); diff --git a/crates/feos-core/src/phase_equilibria/vle_pure.rs b/crates/feos-core/src/phase_equilibria/vle_pure.rs index 5e58601af..e3575c094 100644 --- a/crates/feos-core/src/phase_equilibria/vle_pure.rs +++ b/crates/feos-core/src/phase_equilibria/vle_pure.rs @@ -5,40 +5,37 @@ use crate::errors::{FeosError, FeosResult}; use crate::state::{Contributions, DensityInitialization, State}; use crate::{ReferenceSystem, SolverOptions, TemperatureOrPressure, Verbosity}; use nalgebra::allocator::Allocator; -use nalgebra::{DVector, DefaultAllocator, Dim, dvector}; -use num_dual::{DualNum, DualStruct, Gradients}; -use quantity::{Density, Pressure, RGAS, Temperature}; +use nalgebra::{DVector, DefaultAllocator, Dim, SVector, U1, U2}; +use num_dual::{DualNum, DualStruct, Gradients, gradient, partial}; +use quantity::{Density, Pressure, Temperature}; const SCALE_T_NEW: f64 = 0.7; const MAX_ITER_PURE: usize = 50; const TOL_PURE: f64 = 1e-12; /// # Pure component phase equilibria -impl PhaseEquilibrium { +impl, N: Gradients, D: DualNum + Copy> PhaseEquilibrium +where + DefaultAllocator: Allocator + Allocator + Allocator, +{ /// Calculate a phase equilibrium for a pure component. - pub fn pure( + pub fn pure>( eos: &E, temperature_or_pressure: TP, initial_state: Option<&Self>, options: SolverOptions, ) -> FeosResult { - if let Some(t) = temperature_or_pressure.temperature() { + let (t, rho) = if let Some(t) = temperature_or_pressure.temperature() { let (_, rho) = Self::pure_t(eos, t, initial_state, options)?; - Ok(Self(rho.map(|r| { - State::new_intensive(eos, t, r, &dvector![1.0]).unwrap() - }))) + (t, rho) } else if let Some(p) = temperature_or_pressure.pressure() { - Self::pure_p(eos, p, initial_state, options) + Self::pure_p(eos, p, initial_state, options)? } else { unreachable!() - } + }; + Ok(Self(rho.map(|r| State::new_pure(eos, t, r).unwrap()))) } -} -impl, N: Gradients, D: DualNum + Copy> PhaseEquilibrium -where - DefaultAllocator: Allocator, -{ /// Calculate a phase equilibrium for a pure component /// and given temperature. pub fn pure_t( @@ -89,9 +86,9 @@ where for _ in 0..D::NDERIV { let v_l = liquid_density.recip(); let v_v = vapor_density.recip(); - let (f_l, p_l, dp_l) = eos.p_dpdrho(t, liquid_density, &x); - let (f_v, p_v, dp_v) = eos.p_dpdrho(t, vapor_density, &x); - pressure = -(f_l * v_l - f_v * v_v + t * (v_v / v_l).ln()) / (v_l - v_v); + let (a_l, p_l, dp_l) = eos.p_dpdrho(t, liquid_density, &x); + let (a_v, p_v, dp_v) = eos.p_dpdrho(t, vapor_density, &x); + pressure = -(a_l - a_v + t * (v_v / v_l).ln()) / (v_l - v_v); liquid_density += (pressure - p_l) / dp_l; vapor_density += (pressure - p_v) / dp_v; } @@ -133,15 +130,14 @@ where for i in 1..=max_iter { // calculate properties - let (f_l_res, p_l, p_rho_l) = eos.p_dpdrho(temperature, liquid_density, &x); - let (f_v_res, p_v, p_rho_v) = eos.p_dpdrho(temperature, vapor_density, &x); + let (a_l_res, p_l, p_rho_l) = eos.p_dpdrho(temperature, liquid_density, &x); + let (a_v_res, p_v, p_rho_v) = eos.p_dpdrho(temperature, vapor_density, &x); // Estimate the new pressure let v_v = vapor_density.recip(); let v_l = liquid_density.recip(); let delta_v = v_v - v_l; - let delta_a = - f_v_res * v_v - f_l_res * v_l + temperature * (vapor_density / liquid_density).ln(); + let delta_a = a_v_res - a_l_res + temperature * (vapor_density / liquid_density).ln(); let mut p_new = -delta_a / delta_v; // If the pressure becomes negative, assume the gas phase is ideal. The @@ -238,198 +234,248 @@ where Ok((p, [rho_v, rho_l])) } -impl PhaseEquilibrium { - fn new_pt(eos: &E, temperature: Temperature, pressure: Pressure) -> FeosResult { - let liquid = State::new_xpt( - eos, - temperature, - pressure, - &dvector![1.0], - Some(DensityInitialization::Liquid), - )?; - let vapor = State::new_xpt( - eos, - temperature, - pressure, - &dvector![1.0], - Some(DensityInitialization::Vapor), - )?; - Ok(Self([vapor, liquid])) - } - +impl, N: Gradients, D: DualNum + Copy> PhaseEquilibrium +where + DefaultAllocator: Allocator + Allocator + Allocator, +{ /// Calculate a phase equilibrium for a pure component /// and given pressure. - fn pure_p( + pub fn pure_p( eos: &E, - pressure: Pressure, + pressure: Pressure, initial_state: Option<&Self>, options: SolverOptions, - ) -> FeosResult { - let (max_iter, tol, verbosity) = options.unwrap_or(MAX_ITER_PURE, TOL_PURE); + ) -> FeosResult<(Temperature, [Density; 2])> { + let eos_f64 = eos.re(); + let p = pressure.into_reduced(); // Initialize the phase equilibrium - let mut vle = match initial_state { - Some(init) => init - .clone() - .update_pressure(init.vapor().temperature, pressure)?, - None => PhaseEquilibrium::init_pure_p(eos, pressure)?, + let vle = match initial_state { + Some(init) => ( + init.vapor().temperature.into_reduced().re(), + [ + init.vapor().density.into_reduced().re(), + init.liquid().density.into_reduced().re(), + ], + ), + None => init_pure_p(&eos_f64, pressure.re())?, }; + let (t, [rho_v, rho_l]) = iterate_pure_p(&eos_f64, p.re(), vle, options)?; + // Implicit differentiation + let mut temperature = D::from(t); + let mut vapor_density = D::from(rho_v); + let mut liquid_density = D::from(rho_l); + let x = E::pure_molefracs(); + for _ in 0..D::NDERIV { + let v_l = liquid_density.recip(); + let v_v = vapor_density.recip(); + let (a_l, p_l, s_l, p_rho_l, p_t_l) = + eos.p_dpdrho_dpdt(temperature, liquid_density, &x); + let (a_v, p_v, s_v, p_rho_v, p_t_v) = eos.p_dpdrho_dpdt(temperature, vapor_density, &x); + let ln_rho = (v_l / v_v).ln(); + let delta_t = + (p * (v_v - v_l) + (a_v - a_l + temperature * ln_rho)) / (s_v - s_l - ln_rho); + temperature += delta_t; + liquid_density += (p - p_l - p_t_l * delta_t) / p_rho_l; + vapor_density += (p - p_v - p_t_v * delta_t) / p_rho_v; + } + Ok(( + Temperature::from_reduced(temperature), + [ + Density::from_reduced(vapor_density), + Density::from_reduced(liquid_density), + ], + )) + } +} + +/// Calculate a phase equilibrium for a pure component +/// and given pressure. +fn iterate_pure_p, N: Dim>( + eos: &E, + pressure: f64, + (mut temperature, [mut vapor_density, mut liquid_density]): (f64, [f64; 2]), + options: SolverOptions, +) -> FeosResult<(f64, [f64; 2])> +where + DefaultAllocator: Allocator, +{ + let (max_iter, tol, verbosity) = options.unwrap_or(MAX_ITER_PURE, TOL_PURE); + let x = E::pure_molefracs(); + + log_iter!( + verbosity, + " iter | residual | temperature | liquid density | vapor density " + ); + log_iter!(verbosity, "{:-<89}", ""); + log_iter!( + verbosity, + " {:4} | | {:13.8} | {:12.8} | {:12.8}", + 0, + Temperature::from_reduced(temperature), + Density::from_reduced(liquid_density), + Density::from_reduced(vapor_density) + ); + for i in 1..=max_iter { + // calculate properties + let (a_l_res, p_l, s_l_res, p_rho_l, p_t_l) = + eos.p_dpdrho_dpdt(temperature, liquid_density, &x); + let (a_v_res, p_v, s_v_res, p_rho_v, p_t_v) = + eos.p_dpdrho_dpdt(temperature, vapor_density, &x); + + // calculate the molar volumes + let v_l = liquid_density.recip(); + let v_v = vapor_density.recip(); + + // estimate the temperature steps + let ln_rho = (v_l / v_v).ln(); + let delta_t = (pressure * (v_v - v_l) + (a_v_res - a_l_res + temperature * ln_rho)) + / (s_v_res - s_l_res - ln_rho); + temperature += delta_t; + + // calculate Newton steps for the densities and update state. + let rho_l = liquid_density + (pressure - p_l - p_t_l * delta_t) / p_rho_l; + let rho_v = vapor_density + (pressure - p_v - p_t_v * delta_t) / p_rho_v; + + if rho_l.is_sign_negative() || rho_v.is_sign_negative() || delta_t.abs() > 1.0 { + // if densities are negative or the temperature step is large use density iteration instead + liquid_density = _density_iteration( + eos, + temperature, + pressure, + &x, + DensityInitialization::InitialDensity(liquid_density), + )?; + vapor_density = _density_iteration( + eos, + temperature, + pressure, + &x, + DensityInitialization::InitialDensity(vapor_density), + )?; + } else { + liquid_density = rho_l; + vapor_density = rho_v; + } + + // check for trivial solution + if (vapor_density / liquid_density - 1.0).abs() < TRIVIAL_REL_DEVIATION { + return Err(FeosError::TrivialSolution); + } + + // check for convergence + let res = delta_t.abs(); log_iter!( verbosity, - " iter | residual | temperature | liquid density | vapor density " - ); - log_iter!(verbosity, "{:-<89}", ""); - log_iter!( - verbosity, - " {:4} | | {:13.8} | {:12.8} | {:12.8}", - 0, - vle.vapor().temperature, - vle.liquid().density, - vle.vapor().density + " {:4} | {:14.8e} | {:13.8} | {:12.8} | {:12.8}", + i, + res, + Temperature::from_reduced(temperature), + Density::from_reduced(liquid_density), + Density::from_reduced(vapor_density) ); - for i in 1..=max_iter { - // calculate the pressures and derivatives - let (p_l, p_rho_l) = vle.liquid().p_dpdrho(); - let (p_v, p_rho_v) = vle.vapor().p_dpdrho(); - let p_t_l = vle.liquid().dp_dt(Contributions::Total); - let p_t_v = vle.vapor().dp_dt(Contributions::Total); - - // calculate the residual molar entropies (already cached) - let s_l_res = vle.liquid().residual_molar_entropy(); - let s_v_res = vle.vapor().residual_molar_entropy(); - - // calculate the residual molar Helmholtz energies (already cached) - let a_l_res = vle.liquid().residual_molar_helmholtz_energy(); - let a_v_res = vle.vapor().residual_molar_helmholtz_energy(); - - // calculate the molar volumes - let v_l = 1.0 / vle.liquid().density; - let v_v = 1.0 / vle.vapor().density; - - // estimate the temperature steps - let kt = RGAS * vle.vapor().temperature; - let ln_rho = (v_l / v_v).into_value().ln(); - let delta_t = (pressure * (v_v - v_l) + (a_v_res - a_l_res + kt * ln_rho)) - / (s_v_res - s_l_res - RGAS * ln_rho); - let t_new = vle.vapor().temperature + delta_t; - - // calculate Newton steps for the densities and update state. - let rho_l = vle.liquid().density + (pressure - p_l - p_t_l * delta_t) / p_rho_l; - let rho_v = vle.vapor().density + (pressure - p_v - p_t_v * delta_t) / p_rho_v; - - if rho_l.is_sign_negative() - || rho_v.is_sign_negative() - || delta_t.abs() > Temperature::from_reduced(1.0) - { - // if densities are negative or the temperature step is large use density iteration instead - vle = vle - .update_pressure(t_new, pressure)? - .check_trivial_solution()?; - } else { - // update state - vle = Self([ - State::new_pure(eos, t_new, rho_v)?, - State::new_pure(eos, t_new, rho_l)?, - ]); - } - - // check for convergence - let res = delta_t.abs(); - log_iter!( + if res < temperature * tol { + log_result!( verbosity, - " {:4} | {:14.8e} | {:13.8} | {:12.8} | {:12.8}", - i, - res, - vle.vapor().temperature, - vle.liquid().density, - vle.vapor().density + "PhaseEquilibrium::pure_p: calculation converged in {} step(s)\n", + i ); - if res < vle.vapor().temperature * tol { - log_result!( - verbosity, - "PhaseEquilibrium::pure_p: calculation converged in {} step(s)\n", - i - ); - return Ok(vle); - } + return Ok((temperature, [vapor_density, liquid_density])); } - Err(FeosError::NotConverged("pure_p".to_owned())) } + Err(FeosError::NotConverged("pure_p".to_owned())) +} - /// Initialize a new VLE for a pure substance for a given pressure. - fn init_pure_p(eos: &E, pressure: Pressure) -> FeosResult { - let trial_temperatures = [ - Temperature::from_reduced(300.0), - Temperature::from_reduced(500.0), - Temperature::from_reduced(200.0), - ]; - let x = dvector![1.0]; - let mut vle = None; - let mut t0 = Temperature::from_reduced(1.0); - for t in trial_temperatures.iter() { - t0 = *t; - let _vle = PhaseEquilibrium::new_pt(eos, *t, pressure)?; - if !Self::is_trivial_solution(_vle.vapor(), _vle.liquid()) { - return Ok(_vle); +/// Initialize a new VLE for a pure substance for a given pressure. +fn init_pure_p, N: Gradients>( + eos: &E, + pressure: Pressure, +) -> FeosResult<(f64, [f64; 2])> +where + DefaultAllocator: Allocator + Allocator + Allocator, +{ + let trial_temperatures = [300.0, 500.0, 200.0]; + let p = pressure.into_reduced(); + let x = E::pure_molefracs(); + let mut vle = None; + for t in trial_temperatures { + let liquid_density = _density_iteration(eos, t, p, &x, DensityInitialization::Liquid)?; + let vapor_density = _density_iteration(eos, t, p, &x, DensityInitialization::Vapor)?; + let _vle = (t, [vapor_density, liquid_density]); + if (vapor_density / liquid_density - 1.0).abs() >= TRIVIAL_REL_DEVIATION { + return Ok(_vle); + } + vle = Some(_vle); + } + let Some((t0, [mut rho_v, mut rho_l])) = vle else { + unreachable!() + }; + let [mut t_v, mut t_l] = [t0, t0]; + + let cp = State::critical_point(eos, None, None, None, SolverOptions::default())?; + let cp_density = cp.density.into_reduced(); + if pressure > cp.pressure(Contributions::Total) { + return Err(FeosError::SuperCritical); + }; + + if rho_v < cp_density { + // reduce temperature of liquid phase... + for _ in 0..8 { + t_l *= SCALE_T_NEW; + rho_l = _density_iteration(eos, t_l, p, &x, DensityInitialization::Liquid)?; + if rho_l > cp_density { + break; } - vle = Some(_vle); } - - let cp = State::critical_point(eos, None, None, None, SolverOptions::default())?; - if pressure > cp.pressure(Contributions::Total) { - return Err(FeosError::SuperCritical); - }; - if let Some(mut e) = vle { - if e.vapor().density < cp.density { - for _ in 0..8 { - t0 *= SCALE_T_NEW; - e.0[1] = - State::new_xpt(eos, t0, pressure, &x, Some(DensityInitialization::Liquid))?; - if e.liquid().density > cp.density { - break; - } - } - } else { - for _ in 0..8 { - t0 /= SCALE_T_NEW; - e.0[0] = - State::new_xpt(eos, t0, pressure, &x, Some(DensityInitialization::Vapor))?; - if e.vapor().density < cp.density { - break; - } - } + } else { + // ...or increase temperature of vapor phase + for _ in 0..8 { + t_v /= SCALE_T_NEW; + rho_v = _density_iteration(eos, t_v, p, &x, DensityInitialization::Vapor)?; + if rho_v < cp_density { + break; } + } + } - for _ in 0..20 { - let h = |s: &State<_>| s.residual_enthalpy() + s.total_moles * RGAS * s.temperature; - t0 = (h(e.vapor()) - h(e.liquid())) - / (e.vapor().residual_entropy() - - e.liquid().residual_entropy() - - RGAS - * e.vapor().total_moles - * ((e.vapor().density / e.liquid().density).into_value().ln())); - let trial_state = - State::new_xpt(eos, t0, pressure, &x, Some(DensityInitialization::Vapor))?; - if trial_state.density < cp.density { - e.0[0] = trial_state; - } - let trial_state = - State::new_xpt(eos, t0, pressure, &x, Some(DensityInitialization::Liquid))?; - if trial_state.density > cp.density { - e.0[1] = trial_state; - } - if e.liquid().temperature == e.vapor().temperature { - return Ok(e); - } - } - Err(FeosError::IterationFailed( - "new_init_p: could not find proper initial state".to_owned(), - )) - } else { - unreachable!() + // determine new temperatures and assign them to either the liquid or the vapor phase until + // both phases have the same temperature + for _ in 0..20 { + let h_s = |t, v| { + let (a_res, da_res) = gradient::<_, _, _, U2, _>( + partial( + |t_v: SVector<_, _>, x| { + let [[t, v]] = t_v.data.0; + eos.lift().residual_molar_helmholtz_energy(t, v, x) + }, + &x, + ), + &SVector::from([t, v]), + ); + let [[da_res_dt, da_res_dv]] = da_res.data.0; + (a_res - t * da_res_dt - v * da_res_dv + t, -da_res_dt) + }; + let (h_l, s_l_res) = h_s(t_l, rho_l.recip()); + let (h_v, s_v_res) = h_s(t_v, rho_v.recip()); + let t = (h_v - h_l) / (s_v_res - s_l_res - (rho_v / rho_l).ln()); + let trial_density = _density_iteration(eos, t, p, &x, DensityInitialization::Vapor)?; + if trial_density < cp_density { + rho_v = trial_density; + t_v = t; + } + let trial_density = _density_iteration(eos, t, p, &x, DensityInitialization::Liquid)?; + if trial_density > cp_density { + rho_l = trial_density; + t_l = t; + } + if t_l == t_v { + return Ok((t_l, [rho_v, rho_l])); } } + Err(FeosError::IterationFailed( + "new_init_p: could not find proper initial state".to_owned(), + )) } impl PhaseEquilibrium { @@ -453,7 +499,7 @@ impl PhaseEquilibrium { .map(|i| { let pure_eos = eos.subset(&[i]); PhaseEquilibrium::pure_p(&pure_eos, pressure, None, SolverOptions::default()) - .map(|vle| vle.vapor().temperature) + .map(|(t, _)| t) .ok() }) .collect() diff --git a/crates/feos-core/src/state/residual_properties.rs b/crates/feos-core/src/state/residual_properties.rs index 179ff6d19..174b3697a 100644 --- a/crates/feos-core/src/state/residual_properties.rs +++ b/crates/feos-core/src/state/residual_properties.rs @@ -233,15 +233,6 @@ where .into_value() } - // This function is designed specifically for use in density iterations - pub(crate) fn p_dpdrho(&self) -> (Pressure, as Div>>::Output) { - let dp_dv = self.dp_dv(Contributions::Total); - ( - self.pressure(Contributions::Total), - (-self.volume * dp_dv / self.density), - ) - } - /// Partial molar volume: $v_i=\left(\frac{\partial V}{\partial N_i}\right)_{T,p,N_j}$ pub fn partial_molar_volume(&self) -> MolarVolume> { -self.dp_dni(Contributions::Total) / self.dp_dv(Contributions::Total) diff --git a/crates/feos/src/pcsaft/eos/mod.rs b/crates/feos/src/pcsaft/eos/mod.rs index 3a20c90e4..5edd0f0eb 100644 --- a/crates/feos/src/pcsaft/eos/mod.rs +++ b/crates/feos/src/pcsaft/eos/mod.rs @@ -685,6 +685,50 @@ mod tests_parameter_fit { Ok(()) } + #[test] + fn test_boiling_temperature_derivatives_fit() -> FeosResult<()> { + let pcsaft = pcsaft_non_assoc(); + let pcsaft_ad = pcsaft.named_derivatives(["m", "sigma", "epsilon_k"]); + let pressure = BAR; + let t = pcsaft_ad.boiling_temperature(pressure)?; + let t = t.convert_into(KELVIN); + let (t, grad) = (t.re, t.eps.unwrap_generic(U3, U1)); + + println!("{t:.5}"); + println!("{grad:.5?}"); + + let (t_check, _) = PhaseEquilibrium::pure_p( + &pcsaft_ad, + Pressure::from_inner(&pressure), + None, + Default::default(), + )?; + let t_check = t_check.convert_into(KELVIN); + let (t_check, grad_check) = (t_check.re, t_check.eps.unwrap_generic(U3, U1)); + println!("{t_check:.5}"); + println!("{grad_check:.5?}"); + assert_relative_eq!(t, t_check, max_relative = 1e-15); + assert_relative_eq!(grad, grad_check, max_relative = 1e-15); + + for (i, par) in ["m", "sigma", "epsilon_k"].into_iter().enumerate() { + let mut params = pcsaft.0; + let h = params[i] * 1e-8; + params[i] += h; + let pcsaft_h = PcSaftPure(params); + let (t_h, _) = PhaseEquilibrium::pure_p(&pcsaft_h, pressure, None, Default::default())?; + let dt_h = (t_h.convert_into(KELVIN) - t) / h; + let dt = grad[i]; + println!( + "{par:12}: {:11.5} {:11.5} {:.3e}", + dt_h, + dt, + ((dt_h - dt) / dt).abs() + ); + assert_relative_eq!(dt, dt_h, max_relative = 1e-6); + } + Ok(()) + } + #[test] fn test_equilibrium_liquid_density_derivatives_fit() -> FeosResult<()> { let pcsaft = pcsaft_non_assoc(); diff --git a/py-feos/src/ad/mod.rs b/py-feos/src/ad/mod.rs index 1286f6c14..c18d1eacb 100644 --- a/py-feos/src/ad/mod.rs +++ b/py-feos/src/ad/mod.rs @@ -57,6 +57,32 @@ pub fn vapor_pressure_derivatives<'py>( _vapor_pressure_derivatives(model, parameter_names, parameters, input) } +/// Calculate boiling temperatures and derivatives w.r.t. model parameters. +/// +/// Parameters +/// ---------- +/// model: EquationOfStateAD +/// The equation of state to use. +/// parameter_names: List[string] +/// The name of the parameters for which derivatives are calculated. +/// parameters: np.ndarray[float] +/// The parameters for every data point. +/// input: np.ndarray[float] +/// The pressure (in Pa) for every data point. +/// +/// Returns +/// ------- +/// (np.ndarray[float], np.ndarray[float], np.ndarray[bool]): The boiling temperature (in K), gradients, and convergence status. +#[pyfunction] +pub fn boiling_temperature_derivatives<'py>( + model: PyEquationOfStateAD, + parameter_names: Bound<'py, PyAny>, + parameters: PyReadonlyArray2, + input: PyReadonlyArray2, +) -> GradResult<'py> { + _boiling_temperature_derivatives(model, parameter_names, parameters, input) +} + /// Calculate liquid densities and derivatives w.r.t. model parameters. /// /// Parameters @@ -222,7 +248,7 @@ macro_rules! impl_evaluate_gradients { impl_evaluate_gradients!( pure, - [vapor_pressure, liquid_density, equilibrium_liquid_density], + [vapor_pressure, boiling_temperature, liquid_density, equilibrium_liquid_density], {PcSaftNonAssoc: PcSaftPure, PcSaftFull: PcSaftPure} ); diff --git a/py-feos/src/lib.rs b/py-feos/src/lib.rs index b4f94441a..7ae4d269c 100644 --- a/py-feos/src/lib.rs +++ b/py-feos/src/lib.rs @@ -91,6 +91,7 @@ fn feos(m: &Bound<'_, PyModule>) -> PyResult<()> { #[cfg(feature = "ad")] { m.add_function(wrap_pyfunction!(ad::vapor_pressure_derivatives, m)?)?; + m.add_function(wrap_pyfunction!(ad::boiling_temperature_derivatives, m)?)?; m.add_function(wrap_pyfunction!(ad::liquid_density_derivatives, m)?)?; m.add_function(wrap_pyfunction!( ad::equilibrium_liquid_density_derivatives,