Fix critical point evaluation if one component has 0 composition

prehner · prehner · commit 8e1722369ac7 · 2025-12-17T13:19:28.000+01:00
diff --git a/crates/feos-core/src/state/critical_point.rs b/crates/feos-core/src/state/critical_point.rs
@@ -10,6 +10,7 @@ use num_dual::{
     implicit_derivative_binary, implicit_derivative_vec, jacobian, partial, partial2,
     third_derivative,
 };
+use num_traits::Zero;
 use quantity::{Density, Pressure, Temperature};
 
 const MAX_ITER_CRIT_POINT: usize = 50;
@@ -537,7 +538,13 @@ where
             let n = n + sqrt_z.component_mul(&u).map(Dual3::from_re) * s;
             let t = Dual3::from_re(temperature);
             let v = Dual3::from_re(molar_volume);
-            let ig = n.dot(&n.map(|n| (n / v).ln() - 1.0));
+            let ig = n.dot(&n.map(|n| {
+                if n.is_zero() {
+                    Dual3::zero()
+                } else {
+                    (n / v).ln() - 1.0
+                }
+            }));
             eos.lift().residual_helmholtz_energy(t, v, &n) / t + ig
         },
         D::from(0.0),
diff --git a/crates/feos/tests/pcsaft/critical_point.rs b/crates/feos/tests/pcsaft/critical_point.rs
@@ -1,10 +1,11 @@
 use approx::assert_relative_eq;
 use feos::pcsaft::{PcSaft, PcSaftParameters};
-use feos_core::State;
 use feos_core::parameter::IdentifierOption;
+use feos_core::{SolverOptions, State};
 use nalgebra::dvector;
 use quantity::*;
 use std::error::Error;
+use std::sync::Arc;
 use typenum::P3;
 
 #[test]
@@ -47,3 +48,31 @@ fn test_critical_point_mix() -> Result<(), Box<dyn Error>> {
     );
     Ok(())
 }
+
+#[test]
+fn test_critical_point_limits() -> Result<(), Box<dyn Error>> {
+    let params = PcSaftParameters::from_json(
+        vec!["propane", "butane"],
+        "tests/pcsaft/test_parameters.json",
+        None,
+        IdentifierOption::Name,
+    )?;
+    let options = SolverOptions {
+        verbosity: feos_core::Verbosity::Iter,
+        ..Default::default()
+    };
+    let saft = Arc::new(PcSaft::new(params));
+    let cp_pure = State::critical_point_pure(&saft, None, None, options)?;
+    println!("{} {}", cp_pure[0], cp_pure[1]);
+    let molefracs = dvector![0.0, 1.0];
+    let cp_2 = State::critical_point(&saft, Some(&molefracs), None, None, options)?;
+    println!("{}", cp_2);
+    let molefracs = dvector![1.0, 0.0];
+    let cp_1 = State::critical_point(&saft, Some(&molefracs), None, None, options)?;
+    println!("{}", cp_1);
+    assert_eq!(cp_pure[0].temperature, cp_1.temperature);
+    assert_eq!(cp_pure[0].density, cp_1.density);
+    assert_eq!(cp_pure[1].temperature, cp_2.temperature);
+    assert_eq!(cp_pure[1].density, cp_2.density);
+    Ok(())
+}
