fix tp-flash and add test

prehner · prehner · commit 7eef180dcdf6 · 2026-01-09T15:48:00.000+01:00
diff --git a/crates/feos-core/src/phase_equilibria/tp_flash.rs b/crates/feos-core/src/phase_equilibria/tp_flash.rs
@@ -5,6 +5,7 @@ use crate::state::{Contributions, State};
 use crate::{SolverOptions, Verbosity};
 use nalgebra::{DVector, Matrix3, Matrix4xX};
 use num_dual::{Dual, DualNum, first_derivative};
+use num_traits::Zero;
 use quantity::{Dimensionless, Moles, Pressure, Temperature};
 
 const MAX_ITER_TP: usize = 400;
@@ -270,12 +271,20 @@ impl<E: Residual> PhaseEquilibrium<E, 2> {
                     .liquid()
                     .molefracs
                     .component_div(&self.vapor().molefracs)
-                    .map(|i| if i > 0.0 { i.ln() } else { 0.0 });
+                    .map(|i| i.ln());
 
             // Set residuum to 0 for non-volatile components
             if let Some(nvc) = non_volatile_components.as_ref() {
                 nvc.iter().for_each(|&c| res_vec[c] = 0.0);
             }
+
+            // Set residuum to 0 for non-present components
+            for (i, &x) in self.liquid().molefracs.iter().enumerate() {
+                if x.is_zero() {
+                    res_vec[i] = 0.0
+                }
+            }
+
             let res = res_vec.norm();
             log_iter!(
                 verbosity,
diff --git a/crates/feos/tests/pcsaft/tp_flash.rs b/crates/feos/tests/pcsaft/tp_flash.rs
@@ -4,7 +4,6 @@ use feos_core::parameter::IdentifierOption;
 use feos_core::{Contributions, FeosResult, PhaseEquilibrium, SolverOptions};
 use nalgebra::dvector;
 use quantity::*;
-use std::error::Error;
 
 fn read_params(components: Vec<&str>) -> FeosResult<PcSaftParameters> {
     PcSaftParameters::from_json(
@@ -16,7 +15,7 @@ fn read_params(components: Vec<&str>) -> FeosResult<PcSaftParameters> {
 }
 
 #[test]
-fn test_tp_flash() -> Result<(), Box<dyn Error>> {
+fn test_tp_flash() -> FeosResult<()> {
     let propane = PcSaft::new(read_params(vec!["propane"])?);
     let butane = PcSaft::new(read_params(vec!["butane"])?);
     let t = 250.0 * KELVIN;
@@ -63,3 +62,51 @@ fn test_tp_flash() -> Result<(), Box<dyn Error>> {
     );
     Ok(())
 }
+
+#[test]
+fn test_tp_flash_zero_component() -> FeosResult<()> {
+    let eos_full = PcSaft::new(PcSaftParameters::from_json(
+        vec!["propane", "butane", "hexane"],
+        "tests/pcsaft/test_parameters.json",
+        None,
+        IdentifierOption::Name,
+    )?);
+    let eos_binary = PcSaft::new(PcSaftParameters::from_json(
+        vec!["butane", "hexane"],
+        "tests/pcsaft/test_parameters.json",
+        None,
+        IdentifierOption::Name,
+    )?);
+    let options = SolverOptions {
+        verbosity: feos_core::Verbosity::Iter,
+        ..Default::default()
+    };
+    let vle_full = PhaseEquilibrium::tp_flash(
+        &&eos_full,
+        300.0 * KELVIN,
+        1.2 * BAR,
+        &(dvector![0.0, 0.5, 0.5] * MOL),
+        None,
+        options,
+        None,
+    )?;
+    let vle_binary = PhaseEquilibrium::tp_flash(
+        &&eos_binary,
+        300.0 * KELVIN,
+        1.2 * BAR,
+        &(dvector![0.5, 0.5] * MOL),
+        None,
+        options,
+        None,
+    )?;
+    println!("{vle_full}\n{vle_binary}");
+    assert_eq!(
+        vle_full.liquid().molefracs[1],
+        vle_binary.liquid().molefracs[0]
+    );
+    assert_eq!(
+        vle_full.vapor().molefracs[1],
+        vle_binary.vapor().molefracs[0]
+    );
+    Ok(())
+}
