From 8213aade5f6151ac43793c811b47211f68c1010b Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Tue, 27 Jan 2026 09:50:15 +0100 Subject: [PATCH] Extend tp-flash to static arrays and AD --- .github/workflows/test.yml | 4 +- .github/workflows/wheels.yml | 4 +- CHANGELOG.md | 3 + crates/feos-core/src/phase_equilibria/mod.rs | 11 +- .../phase_equilibria/stability_analysis.rs | 25 ++- .../src/phase_equilibria/tp_flash.rs | 156 ++++++++++++++---- crates/feos/src/pcsaft/eos/mod.rs | 58 ++++++- 7 files changed, 207 insertions(+), 54 deletions(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 7c7b1ad4a..6ea6192e7 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -2,9 +2,9 @@ name: Test on: push: - branches: [main] + branches: [main, development] pull_request: - branches: [main] + branches: [main, development] env: CARGO_TERM_COLOR: always diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml index 44f16c092..92a98abb5 100644 --- a/.github/workflows/wheels.yml +++ b/.github/workflows/wheels.yml @@ -1,9 +1,9 @@ name: Build Wheels on: push: - branches: [main] + branches: [main, development] pull_request: - branches: [main] + branches: [main, development] jobs: linux: runs-on: ubuntu-latest diff --git a/CHANGELOG.md b/CHANGELOG.md index 31fa64b15..7e7427a82 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,9 @@ 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] +### Added +- Extended tp-flash algorithm to static numbers of components and enabled automatic differentiation for binary systems. [#336](https://github.com/feos-org/feos/pull/336) + ### Packaging - Updated `quantity` dependency to 0.13 and removed the `typenum` dependency. [#323](https://github.com/feos-org/feos/pull/323) diff --git a/crates/feos-core/src/phase_equilibria/mod.rs b/crates/feos-core/src/phase_equilibria/mod.rs index b27d97fbb..9683e49a1 100644 --- a/crates/feos-core/src/phase_equilibria/mod.rs +++ b/crates/feos-core/src/phase_equilibria/mod.rs @@ -3,8 +3,8 @@ use crate::errors::{FeosError, FeosResult}; use crate::state::{DensityInitialization, State}; use crate::{Contributions, ReferenceSystem}; use nalgebra::allocator::Allocator; -use nalgebra::{DVector, DefaultAllocator, Dim, Dyn, OVector}; -use num_dual::{DualNum, DualStruct}; +use nalgebra::{DefaultAllocator, Dim, Dyn, OVector}; +use num_dual::{DualNum, DualStruct, Gradients}; use quantity::{Energy, Moles, Pressure, RGAS, Temperature}; use std::fmt; use std::fmt::Write; @@ -168,7 +168,10 @@ where } } -impl PhaseEquilibrium { +impl, N: Gradients, const P: usize> PhaseEquilibrium +where + DefaultAllocator: Allocator, +{ pub(super) fn update_pressure( mut self, temperature: Temperature, @@ -189,7 +192,7 @@ impl PhaseEquilibrium { pub(super) fn update_moles( &mut self, pressure: Pressure, - moles: [&Moles>; P], + moles: [&Moles>; P], ) -> FeosResult<()> { for (i, s) in self.0.iter_mut().enumerate() { *s = State::new_npt( diff --git a/crates/feos-core/src/phase_equilibria/stability_analysis.rs b/crates/feos-core/src/phase_equilibria/stability_analysis.rs index e9af96e8e..69988f8c2 100644 --- a/crates/feos-core/src/phase_equilibria/stability_analysis.rs +++ b/crates/feos-core/src/phase_equilibria/stability_analysis.rs @@ -3,7 +3,9 @@ use crate::equation_of_state::Residual; use crate::errors::{FeosError, FeosResult}; use crate::state::{Contributions, DensityInitialization, State}; use crate::{ReferenceSystem, SolverOptions, Verbosity}; -use nalgebra::{DMatrix, DVector}; +use nalgebra::allocator::Allocator; +use nalgebra::{DefaultAllocator, OMatrix, OVector, U1}; +use num_dual::Gradients; use num_dual::linalg::LU; use num_dual::linalg::smallest_ev; use quantity::Moles; @@ -16,7 +18,10 @@ const MINIMIZE_KMAX: usize = 100; const ZERO_TPD: f64 = -1E-08; /// # Stability analysis -impl State { +impl, N: Gradients> State +where + DefaultAllocator: Allocator + Allocator, +{ /// Determine if the state is stable, i.e. if a phase split should /// occur or not. pub fn is_stable(&self, options: SolverOptions) -> FeosResult { @@ -26,7 +31,7 @@ impl State { /// Perform a stability analysis. The result is a list of [State]s with /// negative tangent plane distance (i.e. lower Gibbs energy) that can be /// used as initial estimates for a phase equilibrium calculation. - pub fn stability_analysis(&self, options: SolverOptions) -> FeosResult>> { + pub fn stability_analysis(&self, options: SolverOptions) -> FeosResult>> { let mut result = Vec::new(); for i_trial in 0..self.eos.components() + 1 { let phase = if i_trial == self.eos.components() { @@ -59,8 +64,9 @@ impl State { Ok(result) } - fn define_trial_state(&self, dominant_component: usize) -> FeosResult> { + fn define_trial_state(&self, dominant_component: usize) -> FeosResult> { let x_feed = &self.molefracs; + let (n, _) = x_feed.shape_generic(); let (x_trial, phase) = if dominant_component == self.eos.components() { // try an ideal vapor phase @@ -70,7 +76,7 @@ impl State { // try each component as nearly pure phase let factor = (1.0 - X_DOMINANT) / (x_feed.sum() - x_feed[dominant_component]); ( - DVector::from_fn(self.eos.components(), |i, _| { + OVector::from_fn_generic(n, U1, |i, _| { if i == dominant_component { X_DOMINANT } else { @@ -92,7 +98,7 @@ impl State { fn minimize_tpd( &self, - trial: &mut State, + trial: &mut State, options: SolverOptions, ) -> FeosResult<(Option, usize)> { let (max_iter, tol, verbosity) = options.unwrap_or(MINIMIZE_KMAX, MINIMIZE_TOL); @@ -154,9 +160,10 @@ impl State { Err(FeosError::NotConverged(String::from("stability analysis"))) } - fn stability_newton_step(&mut self, di: &DVector, tpd: &mut f64) -> FeosResult { + fn stability_newton_step(&mut self, di: &OVector, tpd: &mut f64) -> FeosResult { // save old values let tpd_old = *tpd; + let (n, _) = di.shape_generic(); // calculate residual and ideal hesse matrix let mut hesse = (self.dln_phi_dnj() * Moles::from_reduced(1.0)).into_value(); @@ -166,7 +173,7 @@ impl State { let sq_y = y.map(f64::sqrt); let gradient = (&ln_y + &lnphi - di).component_mul(&sq_y); - let hesse_ig = DMatrix::identity(self.eos.components(), self.eos.components()); + let hesse_ig = OMatrix::identity_generic(n, n); for i in 0..self.eos.components() { hesse.column_mut(i).component_mul_assign(&(sq_y[i] * &sq_y)); if y[i] > f64::EPSILON { @@ -181,7 +188,7 @@ impl State { // ! (3) objective function (tpd) does not descent // !----------------------------------------------------------------------------- let mut adjust_hessian = true; - let mut hessian: DMatrix; + let mut hessian: OMatrix; let mut eta_h = 1.0; while adjust_hessian { diff --git a/crates/feos-core/src/phase_equilibria/tp_flash.rs b/crates/feos-core/src/phase_equilibria/tp_flash.rs index baefea80b..cd9d9bea3 100644 --- a/crates/feos-core/src/phase_equilibria/tp_flash.rs +++ b/crates/feos-core/src/phase_equilibria/tp_flash.rs @@ -2,16 +2,22 @@ use super::PhaseEquilibrium; use crate::equation_of_state::Residual; use crate::errors::{FeosError, FeosResult}; use crate::state::{Contributions, State}; -use crate::{SolverOptions, Verbosity}; -use nalgebra::{DVector, Matrix3, Matrix4xX}; -use num_dual::{Dual, DualNum, first_derivative}; -use quantity::{Dimensionless, Moles, Pressure, Temperature}; +use crate::{ReferenceSystem, SolverOptions, Verbosity}; +use nalgebra::allocator::Allocator; +use nalgebra::{DefaultAllocator, Dim, Matrix3, OVector, SVector, U1, U2, vector}; +use num_dual::{ + Dual, Dual2Vec, DualNum, DualStruct, Gradients, first_derivative, implicit_derivative_sp, +}; +use quantity::{Dimensionless, MOL, MolarVolume, Moles, Pressure, Quantity, Temperature}; const MAX_ITER_TP: usize = 400; const TOL_TP: f64 = 1e-8; /// # Flash calculations -impl PhaseEquilibrium { +impl, N: Gradients> PhaseEquilibrium +where + DefaultAllocator: Allocator + Allocator, +{ /// Perform a Tp-flash calculation. If no initial values are /// given, the solution is initialized using a stability analysis. /// @@ -21,8 +27,8 @@ impl PhaseEquilibrium { eos: &E, temperature: Temperature, pressure: Pressure, - feed: &Moles>, - initial_state: Option<&PhaseEquilibrium>, + feed: &Moles>, + initial_state: Option<&PhaseEquilibrium>, options: SolverOptions, non_volatile_components: Option>, ) -> FeosResult { @@ -34,8 +40,76 @@ impl PhaseEquilibrium { } } +impl, D: DualNum + Copy> PhaseEquilibrium { + /// Perform a Tp-flash calculation for a binary mixture. + /// Compared to the version of the algorithm for a generic + /// number of components ([tp_flash](PhaseEquilibrium::tp_flash)), + /// this can be used in combination with automatic differentiation. + pub fn tp_flash_binary( + eos: &E, + temperature: Temperature, + pressure: Pressure, + feed: &Moles>, + options: SolverOptions, + ) -> FeosResult { + let z = feed.get(0).convert_into(feed.get(0) + feed.get(1)); + let total_moles = feed.sum(); + let moles = vector![z.re(), 1.0 - z.re()] * MOL; + let vle_re = State::new_npt(&eos.re(), temperature.re(), pressure.re(), &moles, None)? + .tp_flash(None, options, None)?; + + // implicit differentiation + + // specifications + let t = temperature.into_reduced(); + let p = pressure.into_reduced(); + + // molar volume and composition of the two phases + let variables = SVector::from([ + vle_re.liquid().density.into_reduced().recip(), + vle_re.vapor().density.into_reduced().recip(), + vle_re.liquid().molefracs[0], + vle_re.vapor().molefracs[0], + ]); + + // calculate derivatives for molar volumes and compositions (first component) + // with respect to t, p, or z or equation of state parameters + // using implicit differentiation of the minimum in the Gibbs energy + let [[v_l, v_v, x, y]] = implicit_derivative_sp( + |variables, &[t, p, z]: &[_; 3]| { + let [[v_l, v_v, x, y]] = variables.data.0; + let beta = (z - x) / (y - x); + let eos = eos.lift(); + let molar_gibbs_energy = |x: Dual2Vec<_, _, _>, v| { + let molefracs = vector![x, -x + 1.0]; + let a_res = eos.residual_molar_helmholtz_energy(t, v, &molefracs); + let a_ig = (x * (x / v).ln() - (x - 1.0) * ((-x + 1.0) / v).ln() - 1.0) * t; + a_res + a_ig + v * p + }; + // g = a + pv is the potential function for a tp flash using a Helmholtz energy model + // see https://www.sciencedirect.com/science/article/pii/S0378381299000928 + molar_gibbs_energy(y, v_v) * beta - molar_gibbs_energy(x, v_l) * (beta - 1.0) + }, + variables, + &[t, p, z], + ) + .data + .0; + let beta = (z - x) / (y - x); + let state = |x: D, v, phi| { + let volume = MolarVolume::from_reduced(v * phi) * total_moles; + let moles = Quantity::new(vector![x, -x + 1.0] * phi * total_moles.convert_into(MOL)); + State::new_nvt(eos, temperature, volume, &moles) + }; + Ok(Self([state(y, v_v, beta)?, state(x, v_l, -beta + 1.0)?])) + } +} + /// # Flash calculations -impl State { +impl, N: Gradients> State +where + DefaultAllocator: Allocator + Allocator, +{ /// Perform a Tp-flash calculation using the [State] as feed. /// If no initial values are given, the solution is initialized /// using a stability analysis. @@ -44,10 +118,10 @@ impl State { /// containing non-volatile components (e.g. ions). pub fn tp_flash( &self, - initial_state: Option<&PhaseEquilibrium>, + initial_state: Option<&PhaseEquilibrium>, options: SolverOptions, non_volatile_components: Option>, - ) -> FeosResult> { + ) -> FeosResult> { // initialization if let Some(init) = initial_state { let vle = self.tp_flash_( @@ -76,10 +150,10 @@ impl State { pub fn tp_flash_( &self, - mut new_vle_state: PhaseEquilibrium, + mut new_vle_state: PhaseEquilibrium, options: SolverOptions, non_volatile_components: Option>, - ) -> FeosResult> { + ) -> FeosResult> { // set options let (max_iter, tol, verbosity) = options.unwrap_or(MAX_ITER_TP, TOL_TP); @@ -92,8 +166,8 @@ impl State { verbosity, " {:4} | | {:10.8?} | {:10.8?}", 0, - new_vle_state.vapor().molefracs.data.as_vec(), - new_vle_state.liquid().molefracs.data.as_vec(), + new_vle_state.vapor().molefracs.as_slice(), + new_vle_state.liquid().molefracs.as_slice(), ); let mut iter = 0; @@ -169,7 +243,7 @@ impl State { Ok(new_vle_state) } - fn tangent_plane_distance(&self, trial_state: &State) -> f64 { + fn tangent_plane_distance(&self, trial_state: &State) -> f64 { let ln_phi_z = self.ln_phi(); let ln_phi_w = trial_state.ln_phi(); let z = &self.molefracs; @@ -178,20 +252,23 @@ impl State { } } -impl PhaseEquilibrium { +impl, N: Gradients> PhaseEquilibrium +where + DefaultAllocator: Allocator + Allocator, +{ fn accelerated_successive_substitution( &mut self, - feed_state: &State, + feed_state: &State, iter: &mut usize, max_iter: usize, tol: f64, verbosity: Verbosity, non_volatile_components: &Option>, ) -> FeosResult<()> { + let (n, _) = feed_state.molefracs.shape_generic(); for _ in 0..max_iter { // do 5 successive substitution steps and check for convergence - let mut k_vec = Matrix4xX::zeros(self.vapor().eos.components()); - // let mut k_vec = Array::zeros((4, self.vapor().eos.components())); + let mut k_vec = std::array::repeat(OVector::zeros_generic(n, U1)); if self.successive_substitution( feed_state, 5, @@ -213,16 +290,19 @@ impl PhaseEquilibrium { let gibbs = self.total_gibbs_energy(); // extrapolate K values - let delta_vec = k_vec.rows_range(1..) - k_vec.rows_range(..3); - let delta = Matrix3::from_fn(|i, j| delta_vec.row(i).dot(&delta_vec.row(j))); + let delta_vec = [ + &k_vec[1] - &k_vec[0], + &k_vec[2] - &k_vec[1], + &k_vec[3] - &k_vec[2], + ]; + let delta = Matrix3::from_fn(|i, j| delta_vec[i].dot(&delta_vec[j])); let d = delta[(0, 1)] * delta[(0, 1)] - delta[(0, 0)] * delta[(1, 1)]; let a = (delta[(0, 2)] * delta[(0, 1)] - delta[(1, 2)] * delta[(0, 0)]) / d; let b = (delta[(1, 2)] * delta[(0, 1)] - delta[(0, 2)] * delta[(1, 1)]) / d; - let mut k = (k_vec.row(3) - + ((b * delta_vec.row(1) + (a + b) * delta_vec.row(2)) / (1.0 - a - b))) - .map(f64::exp) - .transpose(); + let mut k = (&k_vec[3] + + ((b * &delta_vec[1] + (a + b) * &delta_vec[2]) / (1.0 - a - b))) + .map(f64::exp); // Set k = 0 for non-volatile components if let Some(nvc) = non_volatile_components.as_ref() { @@ -245,10 +325,10 @@ impl PhaseEquilibrium { #[expect(clippy::too_many_arguments)] fn successive_substitution( &mut self, - feed_state: &State, + feed_state: &State, iterations: usize, iter: &mut usize, - k_vec: &mut Option<&mut Matrix4xX>, + k_vec: &mut Option<&mut [OVector; 4]>, abs_tol: f64, verbosity: Verbosity, non_volatile_components: &Option>, @@ -278,8 +358,8 @@ impl PhaseEquilibrium { " {:4} | {:14.8e} | {:.8?} | {:.8?}", iter, res, - self.vapor().molefracs.data.as_vec(), - self.liquid().molefracs.data.as_vec(), + self.vapor().molefracs.as_slice(), + self.liquid().molefracs.as_slice(), ); if res < abs_tol { return Ok(true); @@ -289,16 +369,13 @@ impl PhaseEquilibrium { if let Some(k_vec) = k_vec && i >= iterations - 3 { - k_vec.set_row( - i + 3 - iterations, - &k.map(|ki| if ki > 0.0 { ki.ln() } else { 0.0 }).transpose(), - ); + k_vec[i + 3 - iterations] = k.map(|ki| if ki > 0.0 { ki.ln() } else { 0.0 }); } } Ok(false) } - fn update_states(&mut self, feed_state: &State, k: &DVector) -> FeosResult<()> { + fn update_states(&mut self, feed_state: &State, k: &OVector) -> FeosResult<()> { // calculate vapor phase fraction using Rachford-Rice algorithm let mut beta = self.vapor_phase_fraction(); beta = rachford_rice(&feed_state.molefracs, k, Some(beta))?; @@ -314,7 +391,7 @@ impl PhaseEquilibrium { Ok(()) } - fn vle_init_stability(feed_state: &State) -> FeosResult<(Self, Option)> { + fn vle_init_stability(feed_state: &State) -> FeosResult<(Self, Option)> { let mut stable_states = feed_state.stability_analysis(SolverOptions::default())?; let state1 = stable_states.pop(); let state2 = stable_states.pop(); @@ -331,7 +408,14 @@ impl PhaseEquilibrium { } } -fn rachford_rice(feed: &DVector, k: &DVector, beta_in: Option) -> FeosResult { +fn rachford_rice( + feed: &OVector, + k: &OVector, + beta_in: Option, +) -> FeosResult +where + DefaultAllocator: Allocator, +{ const MAX_ITER: usize = 10; const ABS_TOL: f64 = 1e-6; diff --git a/crates/feos/src/pcsaft/eos/mod.rs b/crates/feos/src/pcsaft/eos/mod.rs index b81d3ac7e..3a20c90e4 100644 --- a/crates/feos/src/pcsaft/eos/mod.rs +++ b/crates/feos/src/pcsaft/eos/mod.rs @@ -596,7 +596,7 @@ mod tests_parameter_fit { use super::*; use approx::assert_relative_eq; use feos_core::DensityInitialization::Liquid; - use feos_core::{Contributions, PropertiesAD, ReferenceSystem}; + use feos_core::{Contributions, PropertiesAD, ReferenceSystem, SolverOptions}; use feos_core::{FeosResult, ParametersAD, PhaseEquilibrium, State}; use nalgebra::{U1, U3, U8, vector}; use num_dual::{DualStruct, DualVec}; @@ -935,4 +935,60 @@ mod tests_parameter_fit { assert_relative_eq!(grad, dt_h, max_relative = 1e-6); Ok(()) } + + #[test] + fn test_tp_flash() -> FeosResult<()> { + let (pcsaft, _) = pcsaft_binary()?; + let pcsaft_ad = pcsaft.named_derivatives(["k_ij"]); + let temperature = 500.0 * KELVIN; + let pressure = 44.6 * BAR; + let x = vector![0.5, 0.5]; + let vle = PhaseEquilibrium::tp_flash_binary( + &pcsaft_ad, + Temperature::from_inner(&temperature), + Pressure::from_inner(&pressure), + &Moles::from_inner(&(x * MOL)), + SolverOptions { + verbosity: feos_core::Verbosity::Iter, + tol: Some(1e-10), + ..Default::default() + }, + )?; + let beta = vle + .vapor() + .total_moles + .convert_into(vle.vapor().total_moles + vle.liquid().total_moles); + let (beta, [[grad]]) = (beta.re, beta.eps.unwrap_generic(U1, U1).data.0); + + println!("{beta:.5}"); + println!("{grad:.5?}"); + + let (params, mut kij) = pcsaft.0; + let h = 1e-7; + kij += h; + let pcsaft_h = PcSaftBinary::new(params, kij); + let vle = PhaseEquilibrium::tp_flash_binary( + &pcsaft_h, + temperature, + pressure, + &(x * MOL), + SolverOptions { + tol: Some(1e-10), + ..Default::default() + }, + )?; + let beta_h = vle + .vapor() + .total_moles + .convert_into(vle.vapor().total_moles + vle.liquid().total_moles); + let dbeta_h = (beta_h - beta) / h; + println!( + "k_ij: {:11.5} {:11.5} {:.3e}", + dbeta_h, + grad, + ((dbeta_h - grad) / grad).abs() + ); + assert_relative_eq!(grad, dbeta_h, max_relative = 1e-4); + Ok(()) + } }