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 all commits
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
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
16 changes: 14 additions & 2 deletions 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 Expand Up @@ -96,7 +102,13 @@ impl<E: Residual> PhaseDiagram<E, 2> {
where
E: Send + Sync,
{
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
2 changes: 1 addition & 1 deletion crates/feos/benches/dual_numbers.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ use typenum::P3;
fn state_pcsaft(n: usize, eos: &PcSaft) -> State<&PcSaft> {
let moles = DVector::from_element(n, 1.0 / n as f64) * 10.0 * MOL;
let molefracs = (&moles / moles.sum()).into_value();
let cp = State::critical_point(&eos, Some(&molefracs), None, Default::default()).unwrap();
let cp = State::critical_point(&eos, Some(&molefracs), None, None, Default::default()).unwrap();
let temperature = 0.8 * cp.temperature;
State::new_nvt(&eos, temperature, cp.volume, &moles).unwrap()
}
Expand Down
2 changes: 1 addition & 1 deletion crates/feos/benches/dual_numbers_saftvrmie.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ use quantity::*;
/// - molefracs (or moles) for equimolar mixture.
fn state_saftvrmie(n: usize, eos: &SaftVRMie) -> State<&SaftVRMie> {
let molefracs = DVector::from_element(n, 1.0 / n as f64);
let cp = State::critical_point(&eos, Some(&molefracs), None, Default::default()).unwrap();
let cp = State::critical_point(&eos, Some(&molefracs), None, None, Default::default()).unwrap();
let temperature = 0.8 * cp.temperature;
State::new_nvt(&eos, temperature, cp.volume, &(molefracs * 10. * MOL)).unwrap()
}
Expand Down
6 changes: 3 additions & 3 deletions crates/feos/benches/state_creation.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,12 +23,12 @@ fn npt<E: Residual>(

/// Evaluate critical point constructor
fn critical_point<E: Residual>((eos, n): (&E, Option<&DVector<f64>>)) {
State::critical_point(eos, n, None, Default::default()).unwrap();
State::critical_point(eos, n, None, None, Default::default()).unwrap();
}

/// Evaluate critical point constructor for binary systems at given T or p
fn critical_point_binary<E: Residual, TP: TemperatureOrPressure>((eos, tp): (&E, TP)) {
State::critical_point_binary(eos, tp, None, None, Default::default()).unwrap();
State::critical_point_binary(eos, tp, None, None, None, Default::default()).unwrap();
}

/// VLE for pure substance for given temperature or pressure
Expand Down Expand Up @@ -69,7 +69,7 @@ fn bench_states<E: Residual>(c: &mut Criterion, group_name: &str, eos: &E) {
let ncomponents = eos.components();
let x = DVector::from_element(ncomponents, 1.0 / ncomponents as f64);
let n = &x * 100.0 * MOL;
let crit = State::critical_point(eos, Some(&x), None, Default::default()).unwrap();
let crit = State::critical_point(eos, Some(&x), None, None, Default::default()).unwrap();
let vle = if ncomponents == 1 {
PhaseEquilibrium::pure(eos, crit.temperature * 0.95, None, Default::default()).unwrap()
} else {
Expand Down
2 changes: 1 addition & 1 deletion crates/feos/src/epcsaft/eos/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -273,7 +273,7 @@ mod tests {
fn critical_point() {
let e = ElectrolytePcSaft::new(propane_parameters()).unwrap();
let t = 300.0 * KELVIN;
let cp = State::critical_point(&&e, None, Some(t), Default::default());
let cp = State::critical_point(&&e, None, Some(t), None, Default::default());
if let Ok(v) = cp {
assert_relative_eq!(v.temperature, 375.1244078318015 * KELVIN, epsilon = 1e-8)
}
Expand Down
2 changes: 1 addition & 1 deletion crates/feos/src/lib.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@
//! let saft = PcSaft::new(parameters);
//!
//! // Define thermodynamic conditions.
//! let critical_point = State::critical_point(&&saft, Some(&dvector![1.0]), None, Default::default())?;
//! let critical_point = State::critical_point(&&saft, Some(&dvector![1.0]), None, None, Default::default())?;
//!
//! // Compute properties.
//! let p = critical_point.pressure(Contributions::Total);
Expand Down
2 changes: 1 addition & 1 deletion crates/feos/src/pcsaft/eos/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -506,7 +506,7 @@ mod tests {
fn critical_point() {
let e = &propane_parameters();
let t = 300.0 * KELVIN;
let cp = State::critical_point(&e, None, Some(t), Default::default());
let cp = State::critical_point(&e, None, Some(t), None, Default::default());
if let Ok(v) = cp {
assert_relative_eq!(v.temperature, 375.1244078318015 * KELVIN, epsilon = 1e-8)
}
Expand Down
4 changes: 2 additions & 2 deletions crates/feos/tests/gc_pcsaft/binary.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,9 +31,9 @@ fn test_binary() -> FeosResult<()> {
#[cfg(feature = "dft")]
let func = &GcPcSaftFunctional::new(parameters_func);
let molefracs = dvector![0.5, 0.5];
let cp = State::critical_point(&eos, Some(&molefracs), None, Default::default())?;
let cp = State::critical_point(&eos, Some(&molefracs), None, None, Default::default())?;
#[cfg(feature = "dft")]
let cp_func = State::critical_point(&func, Some(&molefracs), None, Default::default())?;
let cp_func = State::critical_point(&func, Some(&molefracs), None, None, Default::default())?;
println!("{}", cp.temperature);
#[cfg(feature = "dft")]
println!("{}", cp_func.temperature);
Expand Down
4 changes: 2 additions & 2 deletions crates/feos/tests/gc_pcsaft/dft.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -159,7 +159,7 @@ fn test_dft() -> Result<(), Box<dyn Error>> {
let t = 200.0 * KELVIN;
let w = 150.0 * ANGSTROM;
let points = 2048;
let tc = State::critical_point(&&func, None, None, Default::default())?.temperature;
let tc = State::critical_point(&&func, None, None, None, Default::default())?.temperature;
let vle = PhaseEquilibrium::pure(&&func, t, None, Default::default())?;
let profile = PlanarInterface::from_tanh(&vle, points, w, tc, false).solve(None)?;
println!(
Expand Down Expand Up @@ -253,7 +253,7 @@ fn test_dft_newton() -> Result<(), Box<dyn Error>> {
let t = 200.0 * KELVIN;
let w = 150.0 * ANGSTROM;
let points = 512;
let tc = State::critical_point(&&func, None, None, Default::default())?.temperature;
let tc = State::critical_point(&&func, None, None, None, Default::default())?.temperature;
let vle = PhaseEquilibrium::pure(&&func, t, None, Default::default())?;
let solver = DFTSolver::new(Some(Verbosity::Iter))
.picard_iteration(None, Some(10), None, None)
Expand Down
4 changes: 2 additions & 2 deletions crates/feos/tests/pcsaft/critical_point.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ fn test_critical_point_pure() -> Result<(), Box<dyn Error>> {
)?;
let saft = PcSaft::new(params);
let t = 300.0 * KELVIN;
let cp = State::critical_point(&&saft, None, Some(t), Default::default())?;
let cp = State::critical_point(&&saft, None, Some(t), None, Default::default())?;
assert_relative_eq!(cp.temperature, 375.12441 * KELVIN, max_relative = 1e-8);
assert_relative_eq!(
cp.density,
Expand All @@ -38,7 +38,7 @@ fn test_critical_point_mix() -> Result<(), Box<dyn Error>> {
let saft = PcSaft::new(params);
let t = 300.0 * KELVIN;
let molefracs = dvector![0.5, 0.5];
let cp = State::critical_point(&&saft, Some(&molefracs), Some(t), Default::default())?;
let cp = State::critical_point(&&saft, Some(&molefracs), Some(t), None, Default::default())?;
assert_relative_eq!(cp.temperature, 407.93481 * KELVIN, max_relative = 1e-8);
assert_relative_eq!(
cp.density,
Expand Down
6 changes: 3 additions & 3 deletions crates/feos/tests/pcsaft/dft.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,7 +104,7 @@ fn test_dft_propane() -> Result<(), Box<dyn Error>> {
let t = 200.0 * KELVIN;
let w = 150.0 * ANGSTROM;
let points = 2048;
let tc = State::critical_point(&&func_pure, None, None, Default::default())?.temperature;
let tc = State::critical_point(&&func_pure, None, None, None, Default::default())?.temperature;
let vle_pure = PhaseEquilibrium::pure(&&func_pure, t, None, Default::default())?;
let vle_full = PhaseEquilibrium::pure(&&func_full, t, None, Default::default())?;
let vle_full_vec = PhaseEquilibrium::pure(&&func_full_vec, t, None, Default::default())?;
Expand Down Expand Up @@ -214,7 +214,7 @@ fn test_dft_propane_newton() -> Result<(), Box<dyn Error>> {
let t = 200.0 * KELVIN;
let w = 150.0 * ANGSTROM;
let points = 512;
let tc = State::critical_point(&&func, None, None, Default::default())?.temperature;
let tc = State::critical_point(&&func, None, None, None, Default::default())?.temperature;
let vle = PhaseEquilibrium::pure(&&func, t, None, Default::default())?;
let solver = DFTSolver::new(Some(Verbosity::Iter)).newton(None, None, None, None);
PlanarInterface::from_tanh(&vle, points, w, tc, false).solve(Some(&solver))?;
Expand All @@ -235,7 +235,7 @@ fn test_dft_water() -> Result<(), Box<dyn Error>> {
let t = 400.0 * KELVIN;
let w = 120.0 * ANGSTROM;
let points = 2048;
let tc = State::critical_point(&&func_pure, None, None, Default::default())?.temperature;
let tc = State::critical_point(&&func_pure, None, None, None, Default::default())?.temperature;
let vle_pure = PhaseEquilibrium::pure(&&func_pure, t, None, Default::default())?;
let vle_full_vec = PhaseEquilibrium::pure(&&func_full_vec, t, None, Default::default())?;
let profile_pure = PlanarInterface::from_tanh(&vle_pure, points, w, tc, false).solve(None)?;
Expand Down
2 changes: 1 addition & 1 deletion crates/feos/tests/saftvrmie/critical_properties.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,7 @@ fn critical_properties_pure() {
let option = SolverOptions::default();
let p = parameters.remove(name).unwrap();
let eos = SaftVRMie::new(p);
let cp = State::critical_point(&&eos, None, t0, option).unwrap();
let cp = State::critical_point(&&eos, None, t0, None, option).unwrap();
assert_relative_eq!(cp.temperature, data.0, max_relative = 2e-3);
assert_relative_eq!(
cp.pressure(feos_core::Contributions::Total),
Expand Down
Loading