Skip to content

Commit 43a2d3f

Browse files
prehnerprehner
authored
Pass initial density as optional argument to critical point algorithms (#300)
* Pass initial density as optional argument to critical point algorithms * Fix other tests * Fix doctest
1 parent b4c2976 commit 43a2d3f

File tree

18 files changed

+85
-37
lines changed

18 files changed

+85
-37
lines changed

‎crates/feos-core/src/cubic.rs‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -208,7 +208,7 @@ mod tests {
208208
let parameters = PengRobinsonParameters::new_pure(propane)?;
209209
let pr = PengRobinson::new(parameters);
210210
let options = SolverOptions::new().verbosity(Verbosity::Iter);
211-
let cp = State::critical_point(&&pr, None, None, options)?;
211+
let cp = State::critical_point(&&pr, None, None, None, options)?;
212212
println!("{} {}", cp.temperature, cp.pressure(Contributions::Total));
213213
assert_relative_eq!(cp.temperature, tc * KELVIN, max_relative = 1e-4);
214214
assert_relative_eq!(

‎crates/feos-core/src/phase_equilibria/phase_diagram_binary.rs‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -60,6 +60,7 @@ impl<E: Residual + Subset> PhaseDiagram<E, 2> {
6060
temperature_or_pressure,
6161
None,
6262
None,
63+
None,
6364
SolverOptions::default(),
6465
)?;
6566
let cp_vle = PhaseEquilibrium::from_states(cp.clone(), cp.clone());
@@ -71,6 +72,7 @@ impl<E: Residual + Subset> PhaseDiagram<E, 2> {
7172
temperature_or_pressure,
7273
None,
7374
None,
75+
None,
7476
SolverOptions::default(),
7577
)?;
7678
let cp_vle = PhaseEquilibrium::from_states(cp.clone(), cp.clone());

‎crates/feos-core/src/phase_equilibria/phase_diagram_pure.rs‎

Lines changed: 14 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,13 @@ impl<E: Residual> PhaseDiagram<E, 2> {
3535
) -> FeosResult<Self> {
3636
let mut states = Vec::with_capacity(npoints);
3737

38-
let sc = State::critical_point(eos, None, critical_temperature, SolverOptions::default())?;
38+
let sc = State::critical_point(
39+
eos,
40+
None,
41+
critical_temperature,
42+
None,
43+
SolverOptions::default(),
44+
)?;
3945

4046
let max_temperature = min_temperature
4147
+ (sc.temperature - min_temperature) * ((npoints - 2) as f64 / (npoints - 1) as f64);
@@ -96,7 +102,13 @@ impl<E: Residual> PhaseDiagram<E, 2> {
96102
where
97103
E: Send + Sync,
98104
{
99-
let sc = State::critical_point(eos, None, critical_temperature, SolverOptions::default())?;
105+
let sc = State::critical_point(
106+
eos,
107+
None,
108+
critical_temperature,
109+
None,
110+
SolverOptions::default(),
111+
)?;
100112

101113
let max_temperature = min_temperature
102114
+ (sc.temperature - min_temperature) * ((npoints - 2) as f64 / (npoints - 1) as f64);

‎crates/feos-core/src/phase_equilibria/phase_envelope.rs‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@ impl<E: Residual> PhaseDiagram<E, 2> {
2222
eos,
2323
Some(molefracs),
2424
critical_temperature,
25+
None,
2526
SolverOptions::default(),
2627
)?;
2728

@@ -70,6 +71,7 @@ impl<E: Residual> PhaseDiagram<E, 2> {
7071
eos,
7172
Some(molefracs),
7273
critical_temperature,
74+
None,
7375
SolverOptions::default(),
7476
)?;
7577

@@ -134,6 +136,7 @@ impl<E: Residual> PhaseDiagram<E, 2> {
134136
eos,
135137
Some(molefracs),
136138
critical_temperature,
139+
None,
137140
SolverOptions::default(),
138141
)?;
139142

‎crates/feos-core/src/phase_equilibria/vle_pure.rs‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -376,7 +376,7 @@ impl<E: Residual> PhaseEquilibrium<E, 2> {
376376
vle = Some(_vle);
377377
}
378378

379-
let cp = State::critical_point(eos, None, None, SolverOptions::default())?;
379+
let cp = State::critical_point(eos, None, None, None, SolverOptions::default())?;
380380
if pressure > cp.pressure(Contributions::Total) {
381381
return Err(FeosError::SuperCritical);
382382
};

‎crates/feos-core/src/state/critical_point.rs‎

Lines changed: 32 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -22,12 +22,19 @@ impl<R: Residual + Subset> State<R> {
2222
pub fn critical_point_pure(
2323
eos: &R,
2424
initial_temperature: Option<Temperature>,
25+
initial_density: Option<Density>,
2526
options: SolverOptions,
2627
) -> FeosResult<Vec<Self>> {
2728
(0..eos.components())
2829
.map(|i| {
2930
let pure_eos = eos.subset(&[i]);
30-
let cp = State::critical_point(&pure_eos, None, initial_temperature, options)?;
31+
let cp = State::critical_point(
32+
&pure_eos,
33+
None,
34+
initial_temperature,
35+
initial_density,
36+
options,
37+
)?;
3138
let mut molefracs = DVector::zeros(eos.components());
3239
molefracs[i] = 1.0;
3340
State::new_intensive(eos, cp.temperature, cp.density, &molefracs)
@@ -45,15 +52,21 @@ where
4552
temperature_or_pressure: TP,
4653
initial_temperature: Option<Temperature>,
4754
initial_molefracs: Option<[f64; 2]>,
55+
initial_density: Option<Density>,
4856
options: SolverOptions,
4957
) -> FeosResult<Self> {
5058
let eos_re = eos.re();
5159
let n = N::from_usize(2);
5260
let initial_molefracs = initial_molefracs.unwrap_or([0.5; 2]);
5361
let initial_molefracs = OVector::from_fn_generic(n, U1, |i, _| initial_molefracs[i]);
5462
if let Some(t) = temperature_or_pressure.temperature() {
55-
let [rho0, rho1] =
56-
critical_point_binary_t(&eos_re, t.re(), initial_molefracs, options)?;
63+
let [rho0, rho1] = critical_point_binary_t(
64+
&eos_re,
65+
t.re(),
66+
initial_molefracs,
67+
initial_density,
68+
options,
69+
)?;
5770
let rho = implicit_derivative_binary(
5871
|rho0, rho1, &temperature| {
5972
let rho = [rho0, rho1];
@@ -74,6 +87,7 @@ where
7487
p.re(),
7588
initial_temperature,
7689
initial_molefracs,
90+
initial_density,
7791
options,
7892
)?;
7993
let trho = implicit_derivative_vec::<_, _, _, _, U3>(
@@ -99,21 +113,24 @@ where
99113
eos: &E,
100114
molefracs: Option<&OVector<D, N>>,
101115
initial_temperature: Option<Temperature>,
116+
initial_density: Option<Density>,
102117
options: SolverOptions,
103118
) -> FeosResult<Self> {
104119
let eos_re = eos.re();
105120
let molefracs = molefracs.map_or_else(E::pure_molefracs, |x| x.clone());
106121
let x = &molefracs.map(|x| x.re());
122+
let rho_init = initial_density.map(|r| r.into_reduced());
107123
let trial_temperatures = [300.0, 700.0, 500.0];
108124
let mut t_rho = None;
109125
if let Some(t) = initial_temperature {
110-
t_rho = Some(critical_point_hkm(&eos_re, x, t.into_reduced(), options)?);
126+
let t = t.into_reduced();
127+
t_rho = Some(critical_point_hkm(&eos_re, x, t, rho_init, options)?);
111128
}
112129
for &t in trial_temperatures.iter() {
113130
if t_rho.is_some() {
114131
break;
115132
}
116-
t_rho = critical_point_hkm(&eos_re, x, t, options).ok();
133+
t_rho = critical_point_hkm(&eos_re, x, t, rho_init, options).ok();
117134
}
118135
let Some(t_rho) = t_rho else {
119136
return Err(FeosError::NotConverged(String::from("Critical point")));
@@ -134,10 +151,12 @@ where
134151
)
135152
}
136153
}
154+
137155
fn critical_point_hkm<E: Residual<N>, N: Gradients>(
138156
eos: &E,
139157
molefracs: &OVector<f64, N>,
140158
initial_temperature: f64,
159+
initial_density: Option<f64>,
141160
options: SolverOptions,
142161
) -> FeosResult<[f64; 2]>
143162
where
@@ -147,7 +166,7 @@ where
147166

148167
let mut t = initial_temperature;
149168
let max_density = eos.compute_max_density(molefracs);
150-
let mut rho = 0.3 * max_density;
169+
let mut rho = initial_density.unwrap_or(0.3 * max_density);
151170

152171
log_iter!(
153172
verbosity,
@@ -220,6 +239,7 @@ fn critical_point_binary_t<E: Residual<N>, N: Gradients>(
220239
eos: &E,
221240
temperature: Temperature,
222241
initial_molefracs: OVector<f64, N>,
242+
initial_density: Option<Density>,
223243
options: SolverOptions,
224244
) -> FeosResult<[f64; 2]>
225245
where
@@ -230,7 +250,8 @@ where
230250
let t = temperature.to_reduced();
231251
let n = N::from_usize(2);
232252
let max_density = eos.compute_max_density(&initial_molefracs);
233-
let mut rho = SVector::from([initial_molefracs[0], initial_molefracs[1]]) * 0.3 * max_density;
253+
let rho_init = initial_density.map_or(0.3 * max_density, |r| r.into_reduced());
254+
let mut rho = SVector::from([initial_molefracs[0], initial_molefracs[1]]) * rho_init;
234255

235256
log_iter!(
236257
verbosity,
@@ -302,6 +323,7 @@ fn critical_point_binary_p<E: Residual<N>, N: Gradients>(
302323
pressure: Pressure,
303324
initial_temperature: Option<Temperature>,
304325
initial_molefracs: OVector<f64, N>,
326+
initial_density: Option<Density>,
305327
options: SolverOptions,
306328
) -> FeosResult<[f64; 3]>
307329
where
@@ -312,7 +334,8 @@ where
312334
let p = pressure.to_reduced();
313335
let mut t = initial_temperature.map(|t| t.to_reduced()).unwrap_or(300.0);
314336
let max_density = eos.compute_max_density(&initial_molefracs);
315-
let mut rho = SVector::from([initial_molefracs[0], initial_molefracs[1]]) * 0.3 * max_density;
337+
let rho_init = initial_density.map_or(0.3 * max_density, |r| r.into_reduced());
338+
let mut rho = SVector::from([initial_molefracs[0], initial_molefracs[1]]) * rho_init;
316339

317340
log_iter!(
318341
verbosity,
@@ -393,7 +416,7 @@ where
393416
molefracs: Option<&OVector<f64, N>>,
394417
options: SolverOptions,
395418
) -> FeosResult<[Self; 2]> {
396-
let critical_point = Self::critical_point(eos, molefracs, None, options)?;
419+
let critical_point = Self::critical_point(eos, molefracs, None, None, options)?;
397420
let molefracs = molefracs.map_or_else(E::pure_molefracs, |x| x.clone());
398421
let spinodal_vapor = Self::calculate_spinodal(
399422
eos,

‎crates/feos/benches/dual_numbers.rs‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ use typenum::P3;
2222
fn state_pcsaft(n: usize, eos: &PcSaft) -> State<&PcSaft> {
2323
let moles = DVector::from_element(n, 1.0 / n as f64) * 10.0 * MOL;
2424
let molefracs = (&moles / moles.sum()).into_value();
25-
let cp = State::critical_point(&eos, Some(&molefracs), None, Default::default()).unwrap();
25+
let cp = State::critical_point(&eos, Some(&molefracs), None, None, Default::default()).unwrap();
2626
let temperature = 0.8 * cp.temperature;
2727
State::new_nvt(&eos, temperature, cp.volume, &moles).unwrap()
2828
}

‎crates/feos/benches/dual_numbers_saftvrmie.rs‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@ use quantity::*;
1818
/// - molefracs (or moles) for equimolar mixture.
1919
fn state_saftvrmie(n: usize, eos: &SaftVRMie) -> State<&SaftVRMie> {
2020
let molefracs = DVector::from_element(n, 1.0 / n as f64);
21-
let cp = State::critical_point(&eos, Some(&molefracs), None, Default::default()).unwrap();
21+
let cp = State::critical_point(&eos, Some(&molefracs), None, None, Default::default()).unwrap();
2222
let temperature = 0.8 * cp.temperature;
2323
State::new_nvt(&eos, temperature, cp.volume, &(molefracs * 10. * MOL)).unwrap()
2424
}

‎crates/feos/benches/state_creation.rs‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -23,12 +23,12 @@ fn npt<E: Residual>(
2323

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

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

3434
/// VLE for pure substance for given temperature or pressure
@@ -69,7 +69,7 @@ fn bench_states<E: Residual>(c: &mut Criterion, group_name: &str, eos: &E) {
6969
let ncomponents = eos.components();
7070
let x = DVector::from_element(ncomponents, 1.0 / ncomponents as f64);
7171
let n = &x * 100.0 * MOL;
72-
let crit = State::critical_point(eos, Some(&x), None, Default::default()).unwrap();
72+
let crit = State::critical_point(eos, Some(&x), None, None, Default::default()).unwrap();
7373
let vle = if ncomponents == 1 {
7474
PhaseEquilibrium::pure(eos, crit.temperature * 0.95, None, Default::default()).unwrap()
7575
} else {

‎crates/feos/src/epcsaft/eos/mod.rs‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -273,7 +273,7 @@ mod tests {
273273
fn critical_point() {
274274
let e = ElectrolytePcSaft::new(propane_parameters()).unwrap();
275275
let t = 300.0 * KELVIN;
276-
let cp = State::critical_point(&&e, None, Some(t), Default::default());
276+
let cp = State::critical_point(&&e, None, Some(t), None, Default::default());
277277
if let Ok(v) = cp {
278278
assert_relative_eq!(v.temperature, 375.1244078318015 * KELVIN, epsilon = 1e-8)
279279
}

0 commit comments

Comments
 (0)