From 767c2fb53002f60d40637402a47114e3254473ed Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Tue, 21 Jul 2026 15:08:31 +0200 Subject: [PATCH 1/2] Changed some unsafe unwraps to proper error handling --- .../src/phase_equilibria/bubble_dew.rs | 44 +++++++++++-------- .../phase_equilibria/phase_diagram_binary.rs | 31 +++++++------ .../phase_equilibria/stability_analysis.rs | 4 +- crates/feos-core/src/state/mod.rs | 2 +- 4 files changed, 47 insertions(+), 34 deletions(-) diff --git a/crates/feos-core/src/phase_equilibria/bubble_dew.rs b/crates/feos-core/src/phase_equilibria/bubble_dew.rs index 05dd42849..f956d4d7d 100644 --- a/crates/feos-core/src/phase_equilibria/bubble_dew.rs +++ b/crates/feos-core/src/phase_equilibria/bubble_dew.rs @@ -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() { @@ -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, @@ -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, @@ -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, @@ -306,6 +313,7 @@ where }; // implicit differentiation + // unwraps here are safe let (mut t, mut p) = if iterate_p { ( temperature.unwrap().into_reduced(), @@ -329,7 +337,7 @@ where &mut molar_volume, &mut rho2, Verbosity::None, - ) + )? } else { Self::newton_step_p( eos, @@ -339,7 +347,7 @@ where &mut molar_volume, &mut rho2, Verbosity::None, - ) + )? }; } let state1 = State::new( @@ -372,7 +380,7 @@ where molar_volume: &mut D, partial_density_other_phase: &mut OVector, verbosity: Verbosity, - ) -> f64 { + ) -> FeosResult { // 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); @@ -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 { @@ -430,7 +438,7 @@ where x.as_slice(), true, ); - error + Ok(error) } fn newton_step_p( @@ -441,7 +449,7 @@ where molar_volume: &mut D, partial_density_other_phase: &mut OVector, verbosity: Verbosity, - ) -> f64 { + ) -> FeosResult { // 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); @@ -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 { @@ -506,7 +514,7 @@ where x.as_slice(), true, ); - error + Ok(error) } } @@ -598,7 +606,7 @@ where &mut molar_volume, &mut rho2, options_outer.verbosity, - ) + )? } else { Self::newton_step_p( &state1.eos, @@ -608,7 +616,7 @@ where &mut molar_volume, &mut rho2, options_outer.verbosity, - ) + )? }; *temperature = Temperature::from_reduced(t); *pressure = Pressure::from_reduced(p); diff --git a/crates/feos-core/src/phase_equilibria/phase_diagram_binary.rs b/crates/feos-core/src/phase_equilibria/phase_diagram_binary.rs index 03ac6e1cb..a3b23cbcc 100644 --- a/crates/feos-core/src/phase_equilibria/phase_diagram_binary.rs +++ b/crates/feos-core/src/phase_equilibria/phase_diagram_binary.rs @@ -97,7 +97,7 @@ impl PhaseDiagram { 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 }) } @@ -236,7 +236,7 @@ fn check_for_vlle( states: Vec>, npoints: usize, bubble_dew_options: (SolverOptions, SolverOptions), -) -> Vec> { +) -> FeosResult>> { let n = states.len(); let p: Vec<_> = states .iter() @@ -290,7 +290,7 @@ fn check_for_vlle( 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]), @@ -304,17 +304,17 @@ fn check_for_vlle( 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() @@ -368,7 +368,7 @@ fn check_for_vlle( 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]), @@ -382,21 +382,21 @@ fn check_for_vlle( 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. @@ -508,7 +508,10 @@ impl PhaseEquilibrium { 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, @@ -517,7 +520,9 @@ impl PhaseEquilibrium { } 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, diff --git a/crates/feos-core/src/phase_equilibria/stability_analysis.rs b/crates/feos-core/src/phase_equilibria/stability_analysis.rs index ec2329da6..31f248dfa 100644 --- a/crates/feos-core/src/phase_equilibria/stability_analysis.rs +++ b/crates/feos-core/src/phase_equilibria/stability_analysis.rs @@ -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); diff --git a/crates/feos-core/src/state/mod.rs b/crates/feos-core/src/state/mod.rs index daa7a7ff4..18e7190dc 100644 --- a/crates/feos-core/src/state/mod.rs +++ b/crates/feos-core/src/state/mod.rs @@ -6,7 +6,7 @@ //! * the volume //! //! Internally, all properties are computed using such states as input. -use crate::density_iteration::density_iteration; +use crate::density_iteration::{_density_iteration, density_iteration}; use crate::equation_of_state::Residual; use crate::errors::{FeosError, FeosResult}; use crate::{ReferenceSystem, Total}; From 665f16df459bfa72dc2c94fd3a0cd3bbf71222f6 Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Tue, 21 Jul 2026 15:18:55 +0200 Subject: [PATCH 2/2] Removed unused import --- crates/feos-core/src/state/mod.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/crates/feos-core/src/state/mod.rs b/crates/feos-core/src/state/mod.rs index 18e7190dc..daa7a7ff4 100644 --- a/crates/feos-core/src/state/mod.rs +++ b/crates/feos-core/src/state/mod.rs @@ -6,7 +6,7 @@ //! * the volume //! //! Internally, all properties are computed using such states as input. -use crate::density_iteration::{_density_iteration, density_iteration}; +use crate::density_iteration::density_iteration; use crate::equation_of_state::Residual; use crate::errors::{FeosError, FeosResult}; use crate::{ReferenceSystem, Total};