Skip to content

Fixed some unwraps in VLE routines #367

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 2 commits into from
Jul 21, 2026
Merged

Fixed some unwraps in VLE routines #367

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
767c2fb
Changed some unsafe unwraps to proper error handling
g-bauer Jul 21, 2026
665f16d
Removed unused import
g-bauer Jul 21, 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
44 changes: 26 additions & 18 deletions crates/feos-core/src/phase_equilibria/bubble_dew.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -230,7 +230,9 @@ where
let molefracs_spec_re = molefracs_spec.map(|x| x.re());
let (v1, rho2) = if iterate_p {
// temperature is specified
let temperature_re = temperature_re.as_mut().unwrap();
let temperature_re = temperature_re.as_mut().ok_or(FeosError::Error(
"Temperature information is expected for bubble/dew calculation.".to_string(),
))?;

// First use given initial pressure if applicable
if let Some(p) = pressure_re.as_mut() {
Expand All @@ -253,11 +255,11 @@ where
bubble,
)
.and_then(|(p, x)| {
pressure_re = Some(p);
let p = pressure_re.insert(p);
PhaseEquilibrium::iterate_bubble_dew(
&eos_re,
temperature_re,
pressure_re.as_mut().unwrap(),
p,
&molefracs_spec_re,
molefracs_init.or(Some(&x)),
bubble,
Expand All @@ -274,11 +276,11 @@ where
&molefracs_spec_re,
)
.and_then(|p| {
pressure_re = Some(p);
let p = pressure_re.insert(p);
PhaseEquilibrium::iterate_bubble_dew(
&eos_re,
temperature_re,
pressure_re.as_mut().unwrap(),
p,
&molefracs_spec_re,
molefracs_init,
bubble,
Expand All @@ -290,9 +292,14 @@ where
}
} else {
// pressure is specified
let pressure_re = pressure_re.as_mut().unwrap();

let temperature_re = temperature_re.as_mut().expect("An initial temperature is required for the calculation of bubble/dew points at given pressure!");
let pressure_re = pressure_re.as_mut().ok_or(FeosError::Error(
"Pressure information is expected for bubble/dew calculation.".to_string(),
))?;

let temperature_re = temperature_re
.as_mut()
.ok_or(FeosError::Error(
"An initial temperature is required for the calculation of bubble/dew points at given pressure.".to_string()))?;
PhaseEquilibrium::iterate_bubble_dew(
&eos.re(),
temperature_re,
Expand All @@ -306,6 +313,7 @@ where
};

// implicit differentiation
// unwraps here are safe
let (mut t, mut p) = if iterate_p {
(
temperature.unwrap().into_reduced(),
Expand All @@ -329,7 +337,7 @@ where
&mut molar_volume,
&mut rho2,
Verbosity::None,
)
)?
} else {
Self::newton_step_p(
eos,
Expand All @@ -339,7 +347,7 @@ where
&mut molar_volume,
&mut rho2,
Verbosity::None,
)
)?
};
}
let state1 = State::new(
Expand Down Expand Up @@ -372,7 +380,7 @@ where
molar_volume: &mut D,
partial_density_other_phase: &mut OVector<D, N>,
verbosity: Verbosity,
) -> f64 {
) -> FeosResult<f64> {
// calculate properties
let (p_1, mu_res_1, dp_1, dmu_1) = eos.dmu_drho(temperature, partial_density_other_phase);
let (p_2, mu_res_2, dp_2, dmu_2) = eos.dmu_dv(temperature, *molar_volume, molefracs);
Expand Down Expand Up @@ -409,7 +417,7 @@ where
});

// calculate Newton step
let dx = LU::<_, _, Dyn>::new(jac).unwrap().solve(&f);
let dx = LU::<_, _, Dyn>::new(jac)?.solve(&f);

// apply Newton step
for i in 0..n {
Expand All @@ -430,7 +438,7 @@ where
x.as_slice(),
true,
);
error
Ok(error)
}

fn newton_step_p(
Expand All @@ -441,7 +449,7 @@ where
molar_volume: &mut D,
partial_density_other_phase: &mut OVector<D, N>,
verbosity: Verbosity,
) -> f64 {
) -> FeosResult<f64> {
// calculate properties
let (p_1, mu_res_1, dp_1, dmu_1) = eos.dmu_drho(*temperature, partial_density_other_phase);
let (p_2, mu_res_2, dp_2, dmu_2) = eos.dmu_dv(*temperature, *molar_volume, molefracs);
Expand Down Expand Up @@ -485,7 +493,7 @@ where
});

// calculate Newton step
let dx = LU::<_, _, Dyn>::new(jac).unwrap().solve(&f);
let dx = LU::<_, _, Dyn>::new(jac)?.solve(&f);

// apply Newton step
for i in 0..n {
Expand All @@ -506,7 +514,7 @@ where
x.as_slice(),
true,
);
error
Ok(error)
}
}

Expand Down Expand Up @@ -598,7 +606,7 @@ where
&mut molar_volume,
&mut rho2,
options_outer.verbosity,
)
)?
} else {
Self::newton_step_p(
&state1.eos,
Expand All @@ -608,7 +616,7 @@ where
&mut molar_volume,
&mut rho2,
options_outer.verbosity,
)
)?
};
*temperature = Temperature::from_reduced(t);
*pressure = Pressure::from_reduced(p);
Expand Down
31 changes: 18 additions & 13 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 @@ -97,7 +97,7 @@ impl<E: Residual + Subset> PhaseDiagram<E, 2> {
if !bubble {
states = states.into_iter().rev().collect();
}
let states = check_for_vlle(temperature_or_pressure, states, npoints, bubble_dew_options);
let states = check_for_vlle(temperature_or_pressure, states, npoints, bubble_dew_options)?;
Ok(Self { states })
}

Expand Down Expand Up @@ -236,7 +236,7 @@ fn check_for_vlle<E: Residual + Subset, TP: TemperatureOrPressure>(
states: Vec<PhaseEquilibrium<E, 2>>,
npoints: usize,
bubble_dew_options: (SolverOptions, SolverOptions),
) -> Vec<PhaseEquilibrium<E, 2>> {
) -> FeosResult<Vec<PhaseEquilibrium<E, 2>>> {
let n = states.len();
let p: Vec<_> = states
.iter()
Expand Down Expand Up @@ -290,7 +290,7 @@ fn check_for_vlle<E: Residual + Subset, TP: TemperatureOrPressure>(
let a = matrix![yi[1] - yi[0], yj[0] - yj[1];
(pi[1] - pi[0]).into_reduced(), (pj[0] - pj[1]).into_reduced()];
let b = vector![yj[0] - yi[0], (pj[0] - pi[0]).into_reduced()];
let [[r, s]] = LU::new(a).unwrap().solve(&b).data.0;
let [[r, s]] = LU::new(a)?.solve(&b).data.0;
let (xi, xj, p) = (
xi[0] + r * (xi[1] - xi[0]),
xj[0] + s * (xj[1] - xj[0]),
Expand All @@ -304,17 +304,17 @@ fn check_for_vlle<E: Residual + Subset, TP: TemperatureOrPressure>(
Default::default(),
bubble_dew_options,
) else {
return states;
return Ok(states);
};
let x_hetero = (vlle.liquid1().molefracs[0], vlle.liquid2().molefracs[0]);
return PhaseDiagram::binary_vle(
return Ok(PhaseDiagram::binary_vle(
&states[0].liquid().eos,
tp,
Some(npoints),
Some(x_hetero),
bubble_dew_options,
)
.map_or(states, |dia| dia.states);
.map_or(states, |dia| dia.states));
}
}
} else if let Some(p) = tp.pressure()
Expand Down Expand Up @@ -368,7 +368,7 @@ fn check_for_vlle<E: Residual + Subset, TP: TemperatureOrPressure>(
let a = matrix![yi[1] - yi[0], yj[0] - yj[1];
(ti[1] - ti[0]).into_reduced(), (tj[0] - tj[1]).into_reduced()];
let b = vector![yj[0] - yi[0], (tj[0] - ti[0]).into_reduced()];
let [[r, s]] = LU::new(a).unwrap().solve(&b).data.0;
let [[r, s]] = LU::new(a)?.solve(&b).data.0;
let (xi, xj, t) = (
xi[0] + r * (xi[1] - xi[0]),
xj[0] + s * (xj[1] - xj[0]),
Expand All @@ -382,21 +382,21 @@ fn check_for_vlle<E: Residual + Subset, TP: TemperatureOrPressure>(
Default::default(),
bubble_dew_options,
) else {
return states;
return Ok(states);
};
let x_hetero = (vlle.liquid1().molefracs[0], vlle.liquid2().molefracs[0]);
return PhaseDiagram::binary_vle(
return Ok(PhaseDiagram::binary_vle(
&states[0].liquid().eos,
tp,
Some(npoints),
Some(x_hetero),
bubble_dew_options,
)
.map_or(states, |dia| dia.states);
.map_or(states, |dia| dia.states));
}
}
}
states
Ok(states)
}

/// Phase diagram (Txy or pxy) for a system with heteroazeotropic phase behavior.
Expand Down Expand Up @@ -508,7 +508,10 @@ impl<E: Residual> PhaseEquilibrium<E, 3> {
if iterate_p {
PhaseEquilibrium::heteroazeotrope_t(
eos,
temperature.unwrap(),
temperature.ok_or(FeosError::Error(
"Temperature information is expected for heteroazeotrope calculation."
.to_string(),
))?,
x_init,
pressure,
options,
Expand All @@ -517,7 +520,9 @@ impl<E: Residual> PhaseEquilibrium<E, 3> {
} else {
PhaseEquilibrium::heteroazeotrope_p(
eos,
pressure.unwrap(),
pressure.ok_or(FeosError::Error(
"Pressure information is expected for heteroazeotrope calculation.".to_string(),
))?,
x_init,
temperature,
options,
Expand Down
4 changes: 2 additions & 2 deletions crates/feos-core/src/phase_equilibria/stability_analysis.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -168,9 +168,9 @@ where
// calculate residual and ideal hesse matrix
// TODO: this should not require extensive properties, but I couldn't rewrite it
// quickly without breaking it.
let mut hesse = self.n_dln_phi_dnj() / self.total_moles().unwrap().into_reduced();
let mut hesse = self.n_dln_phi_dnj() / self.total_moles()?.into_reduced();
let lnphi = self.ln_phi();
let y = self.moles().unwrap().into_reduced();
let y = self.moles()?.into_reduced();
let ln_y = y.map(|y| if y > f64::EPSILON { y.ln() } else { 0.0 });
let sq_y = y.map(f64::sqrt);
let gradient = (&ln_y + &lnphi - di).component_mul(&sq_y);
Expand Down
Loading