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
Prev Previous commit
Next Next commit
Fix other tests
  • Loading branch information
@prehner
prehner committed Oct 13, 2025
commit 4cba8823e16800c18f1ef784e643d6af5e5790ee
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 @@ -102,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
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/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