From b09ea331bdf34d7c58c966ff2618a672a37f9f56 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Wed, 17 Dec 2025 13:19:28 +0100 Subject: [PATCH 1/5] Fix critical point evaluation if one component has 0 composition --- crates/feos-core/src/state/critical_point.rs | 9 +++++- crates/feos/tests/pcsaft/critical_point.rs | 31 +++++++++++++++++++- 2 files changed, 38 insertions(+), 2 deletions(-) diff --git a/crates/feos-core/src/state/critical_point.rs b/crates/feos-core/src/state/critical_point.rs index 1f18093e0..e77b95299 100644 --- a/crates/feos-core/src/state/critical_point.rs +++ b/crates/feos-core/src/state/critical_point.rs @@ -10,6 +10,7 @@ use num_dual::{ implicit_derivative_binary, implicit_derivative_vec, jacobian, partial, partial2, third_derivative, }; +use num_traits::Zero; use quantity::{Density, Pressure, Temperature}; const MAX_ITER_CRIT_POINT: usize = 50; @@ -537,7 +538,13 @@ where let n = n + sqrt_z.component_mul(&u).map(Dual3::from_re) * s; let t = Dual3::from_re(temperature); let v = Dual3::from_re(molar_volume); - let ig = n.dot(&n.map(|n| (n / v).ln() - 1.0)); + let ig = n.dot(&n.map(|n| { + if n.is_zero() { + Dual3::zero() + } else { + (n / v).ln() - 1.0 + } + })); eos.lift().residual_helmholtz_energy(t, v, &n) / t + ig }, D::from(0.0), diff --git a/crates/feos/tests/pcsaft/critical_point.rs b/crates/feos/tests/pcsaft/critical_point.rs index 6e4ee3156..9183e5e5d 100644 --- a/crates/feos/tests/pcsaft/critical_point.rs +++ b/crates/feos/tests/pcsaft/critical_point.rs @@ -1,10 +1,11 @@ use approx::assert_relative_eq; use feos::pcsaft::{PcSaft, PcSaftParameters}; -use feos_core::State; use feos_core::parameter::IdentifierOption; +use feos_core::{SolverOptions, State}; use nalgebra::dvector; use quantity::*; use std::error::Error; +use std::sync::Arc; #[test] fn test_critical_point_pure() -> Result<(), Box> { @@ -46,3 +47,31 @@ fn test_critical_point_mix() -> Result<(), Box> { ); Ok(()) } + +#[test] +fn test_critical_point_limits() -> Result<(), Box> { + let params = PcSaftParameters::from_json( + vec!["propane", "butane"], + "tests/pcsaft/test_parameters.json", + None, + IdentifierOption::Name, + )?; + let options = SolverOptions { + verbosity: feos_core::Verbosity::Iter, + ..Default::default() + }; + let saft = Arc::new(PcSaft::new(params)); + let cp_pure = State::critical_point_pure(&saft, None, None, options)?; + println!("{} {}", cp_pure[0], cp_pure[1]); + let molefracs = dvector![0.0, 1.0]; + let cp_2 = State::critical_point(&saft, Some(&molefracs), None, None, options)?; + println!("{}", cp_2); + let molefracs = dvector![1.0, 0.0]; + let cp_1 = State::critical_point(&saft, Some(&molefracs), None, None, options)?; + println!("{}", cp_1); + assert_eq!(cp_pure[0].temperature, cp_1.temperature); + assert_eq!(cp_pure[0].density, cp_1.density); + assert_eq!(cp_pure[1].temperature, cp_2.temperature); + assert_eq!(cp_pure[1].density, cp_2.density); + Ok(()) +} From 7eef180dcdf66c6185a50d2232a930ad8bf4c023 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Fri, 9 Jan 2026 15:48:00 +0100 Subject: [PATCH 2/5] fix tp-flash and add test --- .../src/phase_equilibria/tp_flash.rs | 11 +++- crates/feos/tests/pcsaft/tp_flash.rs | 51 ++++++++++++++++++- 2 files changed, 59 insertions(+), 3 deletions(-) diff --git a/crates/feos-core/src/phase_equilibria/tp_flash.rs b/crates/feos-core/src/phase_equilibria/tp_flash.rs index 1c9c3072a..0a1361dbf 100644 --- a/crates/feos-core/src/phase_equilibria/tp_flash.rs +++ b/crates/feos-core/src/phase_equilibria/tp_flash.rs @@ -5,6 +5,7 @@ use crate::state::{Contributions, State}; use crate::{SolverOptions, Verbosity}; use nalgebra::{DVector, Matrix3, Matrix4xX}; use num_dual::{Dual, DualNum, first_derivative}; +use num_traits::Zero; use quantity::{Dimensionless, Moles, Pressure, Temperature}; const MAX_ITER_TP: usize = 400; @@ -270,12 +271,20 @@ impl PhaseEquilibrium { .liquid() .molefracs .component_div(&self.vapor().molefracs) - .map(|i| if i > 0.0 { i.ln() } else { 0.0 }); + .map(|i| i.ln()); // Set residuum to 0 for non-volatile components if let Some(nvc) = non_volatile_components.as_ref() { nvc.iter().for_each(|&c| res_vec[c] = 0.0); } + + // Set residuum to 0 for non-present components + for (i, &x) in self.liquid().molefracs.iter().enumerate() { + if x.is_zero() { + res_vec[i] = 0.0 + } + } + let res = res_vec.norm(); log_iter!( verbosity, diff --git a/crates/feos/tests/pcsaft/tp_flash.rs b/crates/feos/tests/pcsaft/tp_flash.rs index ffe45e0b4..b31a14998 100644 --- a/crates/feos/tests/pcsaft/tp_flash.rs +++ b/crates/feos/tests/pcsaft/tp_flash.rs @@ -4,7 +4,6 @@ use feos_core::parameter::IdentifierOption; use feos_core::{Contributions, FeosResult, PhaseEquilibrium, SolverOptions}; use nalgebra::dvector; use quantity::*; -use std::error::Error; fn read_params(components: Vec<&str>) -> FeosResult { PcSaftParameters::from_json( @@ -16,7 +15,7 @@ fn read_params(components: Vec<&str>) -> FeosResult { } #[test] -fn test_tp_flash() -> Result<(), Box> { +fn test_tp_flash() -> FeosResult<()> { let propane = PcSaft::new(read_params(vec!["propane"])?); let butane = PcSaft::new(read_params(vec!["butane"])?); let t = 250.0 * KELVIN; @@ -63,3 +62,51 @@ fn test_tp_flash() -> Result<(), Box> { ); Ok(()) } + +#[test] +fn test_tp_flash_zero_component() -> FeosResult<()> { + let eos_full = PcSaft::new(PcSaftParameters::from_json( + vec!["propane", "butane", "hexane"], + "tests/pcsaft/test_parameters.json", + None, + IdentifierOption::Name, + )?); + let eos_binary = PcSaft::new(PcSaftParameters::from_json( + vec!["butane", "hexane"], + "tests/pcsaft/test_parameters.json", + None, + IdentifierOption::Name, + )?); + let options = SolverOptions { + verbosity: feos_core::Verbosity::Iter, + ..Default::default() + }; + let vle_full = PhaseEquilibrium::tp_flash( + &&eos_full, + 300.0 * KELVIN, + 1.2 * BAR, + &(dvector![0.0, 0.5, 0.5] * MOL), + None, + options, + None, + )?; + let vle_binary = PhaseEquilibrium::tp_flash( + &&eos_binary, + 300.0 * KELVIN, + 1.2 * BAR, + &(dvector![0.5, 0.5] * MOL), + None, + options, + None, + )?; + println!("{vle_full}\n{vle_binary}"); + assert_eq!( + vle_full.liquid().molefracs[1], + vle_binary.liquid().molefracs[0] + ); + assert_eq!( + vle_full.vapor().molefracs[1], + vle_binary.vapor().molefracs[0] + ); + Ok(()) +} From 5c2f6063f767ba411be97aded0c0486040c9945b Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Fri, 9 Jan 2026 15:51:51 +0100 Subject: [PATCH 3/5] makes slightly more sense to check for zeros in the feed composition than the liquid composition --- crates/feos-core/src/phase_equilibria/tp_flash.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/crates/feos-core/src/phase_equilibria/tp_flash.rs b/crates/feos-core/src/phase_equilibria/tp_flash.rs index 0a1361dbf..6f4afb1fa 100644 --- a/crates/feos-core/src/phase_equilibria/tp_flash.rs +++ b/crates/feos-core/src/phase_equilibria/tp_flash.rs @@ -279,7 +279,7 @@ impl PhaseEquilibrium { } // Set residuum to 0 for non-present components - for (i, &x) in self.liquid().molefracs.iter().enumerate() { + for (i, &x) in feed_state.molefracs.iter().enumerate() { if x.is_zero() { res_vec[i] = 0.0 } From 899a66d3d4335f78debeffa16b699b36e64690f3 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Fri, 9 Jan 2026 17:27:27 +0100 Subject: [PATCH 4/5] change residual in tp flash to avoid division by 0 --- crates/feos-core/src/phase_equilibria/tp_flash.rs | 15 +-------------- 1 file changed, 1 insertion(+), 14 deletions(-) diff --git a/crates/feos-core/src/phase_equilibria/tp_flash.rs b/crates/feos-core/src/phase_equilibria/tp_flash.rs index 6f4afb1fa..baefea80b 100644 --- a/crates/feos-core/src/phase_equilibria/tp_flash.rs +++ b/crates/feos-core/src/phase_equilibria/tp_flash.rs @@ -5,7 +5,6 @@ use crate::state::{Contributions, State}; use crate::{SolverOptions, Verbosity}; use nalgebra::{DVector, Matrix3, Matrix4xX}; use num_dual::{Dual, DualNum, first_derivative}; -use num_traits::Zero; use quantity::{Dimensionless, Moles, Pressure, Temperature}; const MAX_ITER_TP: usize = 400; @@ -266,25 +265,13 @@ impl PhaseEquilibrium { // check for convergence *iter += 1; - let mut res_vec = ln_phi_l - ln_phi_v - + self - .liquid() - .molefracs - .component_div(&self.vapor().molefracs) - .map(|i| i.ln()); + let mut res_vec = k.component_mul(&self.liquid().molefracs) - &self.vapor().molefracs; // Set residuum to 0 for non-volatile components if let Some(nvc) = non_volatile_components.as_ref() { nvc.iter().for_each(|&c| res_vec[c] = 0.0); } - // Set residuum to 0 for non-present components - for (i, &x) in feed_state.molefracs.iter().enumerate() { - if x.is_zero() { - res_vec[i] = 0.0 - } - } - let res = res_vec.norm(); log_iter!( verbosity, From bbc84941455db6687d14e9edf333fdbbcc6ef474 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 15 Jan 2026 13:57:35 +0100 Subject: [PATCH 5/5] add changelog entry --- CHANGELOG.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index c0502f7ee..93e16be91 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). ## [Unreleased] +## Fixed +- Fixed the calculation of critical points and tp-flashes when one or more of the mixture components have zero composition. [#331](https://github.com/feos-org/feos/pull/331) + ### Packaging - Updated `quantity` dependency to 0.13 and removed the `typenum` dependency. [#323](https://github.com/feos-org/feos/pull/323)