Skip to content

Fixed transposed jacobian layout in binary association #386

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
g-bauer merged 3 commits into from
Sep 8, 2026
Merged

Fixed transposed jacobian layout in binary association #386

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d76fc8a
Fixed transposed jacobian layout in binary association, fixed inconsi…
g-bauer Sep 3, 2026
a94257c
Changed spinodal to return second pressure derivative and got rid of …
g-bauer Sep 5, 2026
42e8cc4
Updated changelog
g-bauer Sep 8, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 off-by-one error for the convergence-check in `_density_iteration`. [#386](https://github.com/feos-org/feos/pull/386)
- Fixed erroneously transposed jacobian layout in binary association of `PcSaftBinary`. [#386](https://github.com/feos-org/feos/pull/386)

## [0.10.1] - 2026-07-24
### Fixed
Expand Down
80 changes: 62 additions & 18 deletions crates/feos-core/src/density_iteration.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,9 +119,8 @@ where
}

let maxiter = 50;
let mut iterations = 0;
let mut converged = false;
'iteration: for k in 0..maxiter {
iterations += 1;
let (_, mut p, mut dp_drho) = eos.p_dpdrho(temperature, rho, molefracs);

// attempt to correct for poor initial density rho_init
Expand Down Expand Up @@ -149,7 +148,7 @@ where
let (_, _, d2pdrho2) = eos.p_dpdrho_d2pdrho2(temperature, rho, molefracs);

if rho > 0.85 * maxdensity {
let (sp_p, sp_rho) =
let (sp_p, sp_rho, _) =
_pressure_spinodal(eos, temperature, initial_density, molefracs)?;
rho = sp_rho;
error = sp_p - pressure;
Expand All @@ -167,7 +166,7 @@ where
rho = (rho * 1.1).min(maxdensity)
}
} else if error.is_sign_positive() && d2pdrho2.is_sign_positive() {
let (sp_p, sp_rho) =
let (sp_p, sp_rho, _) =
_pressure_spinodal(eos, temperature, initial_density, molefracs)?;
rho = sp_rho;
error = sp_p - pressure;
Expand All @@ -177,7 +176,7 @@ where
rho = (rho * 1.1).min(maxdensity)
}
} else if error.is_sign_negative() && d2pdrho2.is_sign_negative() {
let (sp_p, sp_rho) =
let (sp_p, sp_rho, _) =
_pressure_spinodal(eos, temperature, initial_density, molefracs)?;
rho = sp_rho;
error = sp_p - pressure;
Expand All @@ -187,9 +186,9 @@ where
rho *= 0.8
}
} else if error.is_sign_negative() && d2pdrho2.is_sign_positive() {
let (_, rho_l) = _pressure_spinodal(eos, temperature, 0.8 * maxdensity, molefracs)?;
let (_, rho_l) = _pressure_spinodal_branch(eos, temperature, molefracs, Liquid)?;
let (sp_v_p, rho_v) =
_pressure_spinodal(eos, temperature, 0.001 * maxdensity, molefracs)?;
_pressure_spinodal_branch(eos, temperature, molefracs, Vapor)?;
error = sp_v_p - pressure;
if error.is_sign_positive()
&& (initial_density - rho_v).abs() < (initial_density - rho_l).abs()
Expand All @@ -199,9 +198,9 @@ where
rho = (rho_l * 1.1).min(maxdensity)
}
} else if error.is_sign_positive() && d2pdrho2.is_sign_negative() {
let (_, rho_l) = _pressure_spinodal(eos, temperature, 0.8 * maxdensity, molefracs)?;
let (_, rho_l) = _pressure_spinodal_branch(eos, temperature, molefracs, Liquid)?;
let (sp_v_p, rho_v) =
_pressure_spinodal(eos, temperature, 0.001 * maxdensity, molefracs)?;
_pressure_spinodal_branch(eos, temperature, molefracs, Vapor)?;
error = sp_v_p - pressure;
if error.is_sign_negative()
&& (initial_density - rho_v).abs() > (initial_density - rho_l).abs()
Expand All @@ -221,27 +220,33 @@ where
// Newton step
rho += delta_rho;
if error.abs() < f64::max(abstol, rho * reltol) {
converged = true;
break 'iteration;
}
}
if iterations == maxiter + 1 {
Err(FeosError::NotConverged("density_iteration".to_owned()))
} else {
if converged {
Ok(rho)
} else {
Err(FeosError::NotConverged("density_iteration".to_owned()))
}
}

/// Spinodal (dp/drho = 0) closest to `rho_init`. Returns `(p, rho)`.
///
/// Which spinodal is found depends on the initial density.
/// If vapor/liquid branch is required, [`_pressure_spinodal_branch`]
/// can be used, which verifies the result.
pub(crate) fn _pressure_spinodal<E: Residual<N>, N: Dim>(
eos: &E,
temperature: f64,
rho_init: f64,
molefracs: &OVector<f64, N>,
) -> FeosResult<(f64, f64)>
) -> FeosResult<(f64, f64, f64)>
where
DefaultAllocator: Allocator<N>,
{
let maxiter = 30;
let abstol = 1e-8;
let tol = 1e-8;

let maxdensity = eos.compute_max_density(molefracs);
let mut rho = rho_init;
Expand All @@ -256,18 +261,57 @@ where

for _ in 0..maxiter {
let (p, dpdrho, d2pdrho2) = eos.p_dpdrho_d2pdrho2(temperature, rho, molefracs);

if dpdrho.abs() < tol {
return Ok((p, rho, d2pdrho2));
}
let mut delta_rho = -dpdrho / d2pdrho2;
// Check failure mode: d2pdrho2 is zero which makes delta_rho infinite.
if !delta_rho.is_finite() {
return Err(FeosError::NotConverged("pressure_spinodal".to_owned()));
}
if delta_rho.abs() > 0.05 * maxdensity {
delta_rho = 0.05 * maxdensity * delta_rho.signum()
}
delta_rho = delta_rho.max(-rho * 0.95); // prevent stepping to rho < 0.0
delta_rho = delta_rho.min(maxdensity - rho); // prevent stepping to rho > maxdensity
rho += delta_rho;
}
Err(FeosError::NotConverged("pressure_spinodal".to_owned()))
}

if dpdrho.abs() < abstol {
return Ok((p, rho));
/// Spinodal of the requested branch (liquid or vapor). Returns `(p, rho)`.
///
/// Errors if the iteration converged to the other branch.
pub(crate) fn _pressure_spinodal_branch<E: Residual<N>, N: Dim>(
eos: &E,
temperature: f64,
molefracs: &OVector<f64, N>,
branch: DensityInitialization<f64>,
) -> FeosResult<(f64, f64)>
where
DefaultAllocator: Allocator<N>,
{
let maxdensity = eos.compute_max_density(molefracs);
let rho_init = match branch {
Liquid => 0.8 * maxdensity,
Vapor => 0.001 * maxdensity,
InitialDensity(_) => {
return Err(FeosError::Error(String::from(
"`_pressure_spinodal_branch`: branch must be Liquid or Vapor",
)));
}
};
let (p, rho, d2pdrho2) = _pressure_spinodal(eos, temperature, rho_init, molefracs)?;
let (on_branch, name) = match branch {
Liquid => (d2pdrho2 > 0.0, "liquid"),
Vapor => (d2pdrho2 < 0.0, "vapor"),
InitialDensity(_) => unreachable!(),
};
if on_branch {
Ok((p, rho))
} else {
Err(FeosError::UndeterminedState(format!(
"pressure_spinodal: iteration converged to the wrong branch (requested {name} spinodal)"
)))
}
Err(FeosError::NotConverged("pressure_spinodal".to_owned()))
}
7 changes: 3 additions & 4 deletions crates/feos-core/src/phase_equilibria/vle_pure.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
use super::{PhaseEquilibrium, TRIVIAL_REL_DEVIATION};
use crate::density_iteration::{_density_iteration, _pressure_spinodal};
use crate::density_iteration::{_density_iteration, _pressure_spinodal_branch};
use crate::equation_of_state::{Residual, Subset};
use crate::errors::{FeosError, FeosResult};
use crate::state::{Contributions, DensityInitialization, State};
Expand Down Expand Up @@ -228,10 +228,9 @@ where
DefaultAllocator: Allocator<N>,
{
let x = E::pure_molefracs();
let maxdensity = eos.compute_max_density(&x);
let t = temperature.into_reduced();
let (p_l, _) = _pressure_spinodal(eos, t, 0.8 * maxdensity, &x)?;
let (p_v, _) = _pressure_spinodal(eos, t, 0.001 * maxdensity, &x)?;
let (p_l, _) = _pressure_spinodal_branch(eos, t, &x, DensityInitialization::Liquid)?;
let (p_v, _) = _pressure_spinodal_branch(eos, t, &x, DensityInitialization::Vapor)?;
let p = 0.5 * (0.0_f64.max(p_l) + p_v);
let rho_l = _density_iteration(eos, t, p, &x, DensityInitialization::Liquid)?;
let rho_v = _density_iteration(eos, t, p, &x, DensityInitialization::Vapor)?;
Expand Down
Loading