Convergence improvement for tp_flash (#219)

Author: prehner

CHANGELOG.md

@@ -6,6 +6,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+## [0.6.1] - 2024-01-11
+- Python only: Release the changes introduced in `feos-core` 0.6.1.
+
 ## [0.6.0] - 2023-12-19
 ### Added
 - Added `EquationOfState.ideal_gas()` to initialize an equation of state that only consists of an ideal gas contribution. [#204](https://github.com/feos-org/feos/pull/204)

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "feos"
-version = "0.6.0"
+version = "0.6.1"
 authors = ["Gernot Bauer <[email protected]>", "Philipp Rehner <[email protected]>"]
 edition = "2021"
 readme = "README.md"

feos-core/CHANGELOG.md

@@ -6,6 +6,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## Unreleased
 
+## [0.6.1] 2024-01-11
+### Fixed
+- Improved convergence of `tp_flash` for certain edge cases. [#219](https://github.com/feos-org/feos/pull/219)
+
 ## [0.6.0] 2023-12-19
 ### Added
 - Added `EquationOfState::ideal_gas` to initialize an equation of state that only consists of an ideal gas contribution. [#204](https://github.com/feos-org/feos/pull/204)

feos-core/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "feos-core"
-version = "0.6.0"
+version = "0.6.1"
 authors = ["Gernot Bauer <[email protected]>",
            "Philipp Rehner <[email protected]"]
 edition = "2021"

feos-core/src/phase_equilibria/tp_flash.rs

@@ -51,18 +51,42 @@ impl<E: Residual> State<E> {
         initial_state: Option<&PhaseEquilibrium<E, 2>>,
         options: SolverOptions,
         non_volatile_components: Option<Vec<usize>>,
+    ) -> EosResult<PhaseEquilibrium<E, 2>> {
+        // initialization
+        if let Some(init) = initial_state {
+            let vle = self.tp_flash_(
+                init.clone()
+                    .update_pressure(self.temperature, self.pressure(Contributions::Total))?,
+                options,
+                non_volatile_components.clone(),
+            );
+            if vle.is_ok() {
+                return vle;
+            }
+        }
+
+        let (init1, init2) = PhaseEquilibrium::vle_init_stability(self)?;
+        let vle = self.tp_flash_(init1, options, non_volatile_components.clone());
+        if vle.is_ok() {
+            return vle;
+        }
+
+        if let Some(init2) = init2 {
+            self.tp_flash_(init2, options, non_volatile_components)
+        } else {
+            vle
+        }
+    }
+
+    pub fn tp_flash_(
+        &self,
+        mut new_vle_state: PhaseEquilibrium<E, 2>,
+        options: SolverOptions,
+        non_volatile_components: Option<Vec<usize>>,
     ) -> EosResult<PhaseEquilibrium<E, 2>> {
         // set options
         let (max_iter, tol, verbosity) = options.unwrap_or(MAX_ITER_TP, TOL_TP);
 
-        // initialization
-        let mut new_vle_state = match initial_state {
-            Some(init) => init
-                .clone()
-                .update_pressure(self.temperature, self.pressure(Contributions::Total))?,
-            None => PhaseEquilibrium::vle_init_stability(self)?,
-        };
-
         log_iter!(
             verbosity,
             " iter |    residual    |  phase I mole fractions  |  phase II mole fractions  "
@@ -297,14 +321,19 @@ impl<E: Residual> PhaseEquilibrium<E, 2> {
         Ok(())
     }
 
-    fn vle_init_stability(feed_state: &State<E>) -> EosResult<Self> {
+    fn vle_init_stability(feed_state: &State<E>) -> EosResult<(Self, Option<Self>)> {
         let mut stable_states = feed_state.stability_analysis(SolverOptions::default())?;
         let state1 = stable_states.pop();
         let state2 = stable_states.pop();
-        match (state1, state2) {
-            (Some(s1), Some(s2)) => Ok(Self::from_states(s1, s2)),
-            (Some(s1), None) => Ok(Self::from_states(s1, feed_state.clone())),
-            _ => Err(EosError::NoPhaseSplit),
+        if let Some(s1) = state1 {
+            let init1 = Self::from_states(s1.clone(), feed_state.clone());
+            if let Some(s2) = state2 {
+                Ok((Self::from_states(s1, s2), Some(init1)))
+            } else {
+                Ok((init1, None))
+            }
+        } else {
+            Err(EosError::NoPhaseSplit)
         }
     }
 }