Skip to content

Pass initial density as optional argument to critical point algorithms #300

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
Oct 13, 2025
Merged

Pass initial density as optional argument to critical point algorithms #300

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
430300a
Pass initial density as optional argument to critical point algorithms
prehner Oct 13, 2025
4cba882
Fix other tests
prehner Oct 13, 2025
4f0e6c9
Fix doctest
prehner Oct 13, 2025
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
Next Next commit
Pass initial density as optional argument to critical point algorithms
  • Loading branch information
@prehner
prehner committed Oct 13, 2025
commit 430300a2a3753faa54d2cb43560aa7844379c42f
2 changes: 1 addition & 1 deletion crates/feos-core/src/cubic.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -208,7 +208,7 @@ mod tests {
let parameters = PengRobinsonParameters::new_pure(propane)?;
let pr = PengRobinson::new(parameters);
let options = SolverOptions::new().verbosity(Verbosity::Iter);
let cp = State::critical_point(&&pr, None, None, options)?;
let cp = State::critical_point(&&pr, None, None, None, options)?;
println!("{} {}", cp.temperature, cp.pressure(Contributions::Total));
assert_relative_eq!(cp.temperature, tc * KELVIN, max_relative = 1e-4);
assert_relative_eq!(
Expand Down
2 changes: 2 additions & 0 deletions crates/feos-core/src/phase_equilibria/phase_diagram_binary.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,7 @@ impl<E: Residual + Subset> PhaseDiagram<E, 2> {
temperature_or_pressure,
None,
None,
None,
SolverOptions::default(),
)?;
let cp_vle = PhaseEquilibrium::from_states(cp.clone(), cp.clone());
Expand All @@ -71,6 +72,7 @@ impl<E: Residual + Subset> PhaseDiagram<E, 2> {
temperature_or_pressure,
None,
None,
None,
SolverOptions::default(),
)?;
let cp_vle = PhaseEquilibrium::from_states(cp.clone(), cp.clone());
Expand Down
8 changes: 7 additions & 1 deletion crates/feos-core/src/phase_equilibria/phase_diagram_pure.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,13 @@ impl<E: Residual> PhaseDiagram<E, 2> {
) -> FeosResult<Self> {
let mut states = Vec::with_capacity(npoints);

let sc = State::critical_point(eos, None, critical_temperature, SolverOptions::default())?;
let sc = State::critical_point(
eos,
None,
critical_temperature,
None,
SolverOptions::default(),
)?;

let max_temperature = min_temperature
+ (sc.temperature - min_temperature) * ((npoints - 2) as f64 / (npoints - 1) as f64);
Expand Down
3 changes: 3 additions & 0 deletions crates/feos-core/src/phase_equilibria/phase_envelope.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@ impl<E: Residual> PhaseDiagram<E, 2> {
eos,
Some(molefracs),
critical_temperature,
None,
SolverOptions::default(),
)?;

Expand Down Expand Up @@ -70,6 +71,7 @@ impl<E: Residual> PhaseDiagram<E, 2> {
eos,
Some(molefracs),
critical_temperature,
None,
SolverOptions::default(),
)?;

Expand Down Expand Up @@ -134,6 +136,7 @@ impl<E: Residual> PhaseDiagram<E, 2> {
eos,
Some(molefracs),
critical_temperature,
None,
SolverOptions::default(),
)?;

Expand Down
2 changes: 1 addition & 1 deletion crates/feos-core/src/phase_equilibria/vle_pure.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -376,7 +376,7 @@ impl<E: Residual> PhaseEquilibrium<E, 2> {
vle = Some(_vle);
}

let cp = State::critical_point(eos, None, None, SolverOptions::default())?;
let cp = State::critical_point(eos, None, None, None, SolverOptions::default())?;
if pressure > cp.pressure(Contributions::Total) {
return Err(FeosError::SuperCritical);
};
Expand Down
41 changes: 32 additions & 9 deletions crates/feos-core/src/state/critical_point.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,12 +22,19 @@ impl<R: Residual + Subset> State<R> {
pub fn critical_point_pure(
eos: &R,
initial_temperature: Option<Temperature>,
initial_density: Option<Density>,
options: SolverOptions,
) -> FeosResult<Vec<Self>> {
(0..eos.components())
.map(|i| {
let pure_eos = eos.subset(&[i]);
let cp = State::critical_point(&pure_eos, None, initial_temperature, options)?;
let cp = State::critical_point(
&pure_eos,
None,
initial_temperature,
initial_density,
options,
)?;
let mut molefracs = DVector::zeros(eos.components());
molefracs[i] = 1.0;
State::new_intensive(eos, cp.temperature, cp.density, &molefracs)
Expand All @@ -45,15 +52,21 @@ where
temperature_or_pressure: TP,
initial_temperature: Option<Temperature>,
initial_molefracs: Option<[f64; 2]>,
initial_density: Option<Density>,
options: SolverOptions,
) -> FeosResult<Self> {
let eos_re = eos.re();
let n = N::from_usize(2);
let initial_molefracs = initial_molefracs.unwrap_or([0.5; 2]);
let initial_molefracs = OVector::from_fn_generic(n, U1, |i, _| initial_molefracs[i]);
if let Some(t) = temperature_or_pressure.temperature() {
let [rho0, rho1] =
critical_point_binary_t(&eos_re, t.re(), initial_molefracs, options)?;
let [rho0, rho1] = critical_point_binary_t(
&eos_re,
t.re(),
initial_molefracs,
initial_density,
options,
)?;
let rho = implicit_derivative_binary(
|rho0, rho1, &temperature| {
let rho = [rho0, rho1];
Expand All @@ -74,6 +87,7 @@ where
p.re(),
initial_temperature,
initial_molefracs,
initial_density,
options,
)?;
let trho = implicit_derivative_vec::<_, _, _, _, U3>(
Expand All @@ -99,21 +113,24 @@ where
eos: &E,
molefracs: Option<&OVector<D, N>>,
initial_temperature: Option<Temperature>,
initial_density: Option<Density>,
options: SolverOptions,
) -> FeosResult<Self> {
let eos_re = eos.re();
let molefracs = molefracs.map_or_else(E::pure_molefracs, |x| x.clone());
let x = &molefracs.map(|x| x.re());
let rho_init = initial_density.map(|r| r.into_reduced());
let trial_temperatures = [300.0, 700.0, 500.0];
let mut t_rho = None;
if let Some(t) = initial_temperature {
t_rho = Some(critical_point_hkm(&eos_re, x, t.into_reduced(), options)?);
let t = t.into_reduced();
t_rho = Some(critical_point_hkm(&eos_re, x, t, rho_init, options)?);
}
for &t in trial_temperatures.iter() {
if t_rho.is_some() {
break;
}
t_rho = critical_point_hkm(&eos_re, x, t, options).ok();
t_rho = critical_point_hkm(&eos_re, x, t, rho_init, options).ok();
}
let Some(t_rho) = t_rho else {
return Err(FeosError::NotConverged(String::from("Critical point")));
Expand All @@ -134,10 +151,12 @@ where
)
}
}

fn critical_point_hkm<E: Residual<N>, N: Gradients>(
eos: &E,
molefracs: &OVector<f64, N>,
initial_temperature: f64,
initial_density: Option<f64>,
options: SolverOptions,
) -> FeosResult<[f64; 2]>
where
Expand All @@ -147,7 +166,7 @@ where

let mut t = initial_temperature;
let max_density = eos.compute_max_density(molefracs);
let mut rho = 0.3 * max_density;
let mut rho = initial_density.unwrap_or(0.3 * max_density);

log_iter!(
verbosity,
Expand Down Expand Up @@ -220,6 +239,7 @@ fn critical_point_binary_t<E: Residual<N>, N: Gradients>(
eos: &E,
temperature: Temperature,
initial_molefracs: OVector<f64, N>,
initial_density: Option<Density>,
options: SolverOptions,
) -> FeosResult<[f64; 2]>
where
Expand All @@ -230,7 +250,8 @@ where
let t = temperature.to_reduced();
let n = N::from_usize(2);
let max_density = eos.compute_max_density(&initial_molefracs);
let mut rho = SVector::from([initial_molefracs[0], initial_molefracs[1]]) * 0.3 * max_density;
let rho_init = initial_density.map_or(0.3 * max_density, |r| r.into_reduced());
let mut rho = SVector::from([initial_molefracs[0], initial_molefracs[1]]) * rho_init;

log_iter!(
verbosity,
Expand Down Expand Up @@ -302,6 +323,7 @@ fn critical_point_binary_p<E: Residual<N>, N: Gradients>(
pressure: Pressure,
initial_temperature: Option<Temperature>,
initial_molefracs: OVector<f64, N>,
initial_density: Option<Density>,
options: SolverOptions,
) -> FeosResult<[f64; 3]>
where
Expand All @@ -312,7 +334,8 @@ where
let p = pressure.to_reduced();
let mut t = initial_temperature.map(|t| t.to_reduced()).unwrap_or(300.0);
let max_density = eos.compute_max_density(&initial_molefracs);
let mut rho = SVector::from([initial_molefracs[0], initial_molefracs[1]]) * 0.3 * max_density;
let rho_init = initial_density.map_or(0.3 * max_density, |r| r.into_reduced());
let mut rho = SVector::from([initial_molefracs[0], initial_molefracs[1]]) * rho_init;

log_iter!(
verbosity,
Expand Down Expand Up @@ -393,7 +416,7 @@ where
molefracs: Option<&OVector<f64, N>>,
options: SolverOptions,
) -> FeosResult<[Self; 2]> {
let critical_point = Self::critical_point(eos, molefracs, None, options)?;
let critical_point = Self::critical_point(eos, molefracs, None, None, options)?;
let molefracs = molefracs.map_or_else(E::pure_molefracs, |x| x.clone());
let spinodal_vapor = Self::calculate_spinodal(
eos,
Expand Down
20 changes: 14 additions & 6 deletions py-feos/src/state.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -190,19 +190,21 @@ impl PyState {
/// State : State at critical conditions
#[staticmethod]
#[pyo3(
text_signature = "(eos, initial_temperature=None, max_iter=None, tol=None, verbosity=None)"
text_signature = "(eos, initial_temperature=None, initial_density=None, max_iter=None, tol=None, verbosity=None)"
)]
#[pyo3(signature = (eos, initial_temperature=None, max_iter=None, tol=None, verbosity=None))]
#[pyo3(signature = (eos, initial_temperature=None, initial_density=None, max_iter=None, tol=None, verbosity=None))]
fn critical_point_pure(
eos: &PyEquationOfState,
initial_temperature: Option<Temperature>,
initial_density: Option<Density>,
max_iter: Option<usize>,
tol: Option<f64>,
verbosity: Option<PyVerbosity>,
) -> PyResult<Vec<Self>> {
let cp = State::critical_point_pure(
&eos.0,
initial_temperature,
initial_density,
(max_iter, tol, verbosity.map(|v| v.into())).into(),
)
.map_err(PyFeosError::from)?;
Expand Down Expand Up @@ -232,13 +234,14 @@ impl PyState {
/// State : State at critical conditions.
#[staticmethod]
#[pyo3(
text_signature = "(eos, molefracs=None, initial_temperature=None, max_iter=None, tol=None, verbosity=None)"
text_signature = "(eos, molefracs=None, initial_temperature=None, initial_density=None, max_iter=None, tol=None, verbosity=None)"
)]
#[pyo3(signature = (eos, molefracs=None, initial_temperature=None, max_iter=None, tol=None, verbosity=None))]
#[pyo3(signature = (eos, molefracs=None, initial_temperature=None, initial_density=None, max_iter=None, tol=None, verbosity=None))]
fn critical_point<'py>(
eos: &PyEquationOfState,
molefracs: Option<PyReadonlyArray1<'py, f64>>,
initial_temperature: Option<Temperature>,
initial_density: Option<Density>,
max_iter: Option<usize>,
tol: Option<f64>,
verbosity: Option<PyVerbosity>,
Expand All @@ -248,6 +251,7 @@ impl PyState {
&eos.0,
parse_molefracs(molefracs).as_ref(),
initial_temperature,
initial_density,
(max_iter, tol, verbosity.map(|v| v.into())).into(),
)
.map_err(PyFeosError::from)?,
Expand Down Expand Up @@ -278,14 +282,16 @@ impl PyState {
/// State : State at critical conditions.
#[staticmethod]
#[pyo3(
text_signature = "(eos, temperature_or_pressure, initial_temperature=None, initial_molefracs=None, max_iter=None, tol=None, verbosity=None)"
text_signature = "(eos, temperature_or_pressure, initial_temperature=None, initial_molefracs=None, initial_density=None, max_iter=None, tol=None, verbosity=None)"
)]
#[pyo3(signature = (eos, temperature_or_pressure, initial_temperature=None, initial_molefracs=None, max_iter=None, tol=None, verbosity=None))]
#[pyo3(signature = (eos, temperature_or_pressure, initial_temperature=None, initial_molefracs=None, initial_density=None, max_iter=None, tol=None, verbosity=None))]
#[expect(clippy::too_many_arguments)]
fn critical_point_binary(
eos: &PyEquationOfState,
temperature_or_pressure: Bound<'_, PyAny>,
initial_temperature: Option<Temperature>,
initial_molefracs: Option<[f64; 2]>,
initial_density: Option<Density>,
max_iter: Option<usize>,
tol: Option<f64>,
verbosity: Option<PyVerbosity>,
Expand All @@ -298,6 +304,7 @@ impl PyState {
t,
initial_temperature,
initial_molefracs,
initial_density,
(max_iter, tol, v).into(),
)
.map_err(PyFeosError::from)?,
Expand All @@ -309,6 +316,7 @@ impl PyState {
p,
initial_temperature,
initial_molefracs,
initial_density,
(max_iter, tol, v).into(),
)
.map_err(PyFeosError::from)?,
Expand Down
Loading