Extend tp-flash to static arrays and AD

.github/workflows/test.yml

@@ -2,9 +2,9 @@ name: Test
 
 on:
   push:
-    branches: [main]
+    branches: [main, development]
   pull_request:
-    branches: [main]
+    branches: [main, development]
 
 env:
   CARGO_TERM_COLOR: always

.github/workflows/wheels.yml

@@ -1,9 +1,9 @@
 name: Build Wheels
 on:
   push:
-    branches: [main]
+    branches: [main, development]
   pull_request:
-    branches: [main]
+    branches: [main, development]
 jobs:
   linux:
     runs-on: ubuntu-latest

CHANGELOG.md

@@ -5,6 +5,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Breaking]
+### Added
+- Extended tp-flash algorithm to static numbers of components and enabled automatic differentiation for binary systems. [#336](https://github.com/feos-org/feos/pull/336)
+
 ### Packaging
 - Updated `quantity` dependency to 0.13 and removed the `typenum` dependency. [#323](https://github.com/feos-org/feos/pull/323)
 

crates/feos-core/src/phase_equilibria/mod.rs

@@ -3,8 +3,8 @@ use crate::errors::{FeosError, FeosResult};
 use crate::state::{DensityInitialization, State};
 use crate::{Contributions, ReferenceSystem};
 use nalgebra::allocator::Allocator;
-use nalgebra::{DVector, DefaultAllocator, Dim, Dyn, OVector};
-use num_dual::{DualNum, DualStruct};
+use nalgebra::{DefaultAllocator, Dim, Dyn, OVector};
+use num_dual::{DualNum, DualStruct, Gradients};
 use quantity::{Energy, Moles, Pressure, RGAS, Temperature};
 use std::fmt;
 use std::fmt::Write;
@@ -168,7 +168,10 @@ where
     }
 }
 
-impl<E: Residual, const P: usize> PhaseEquilibrium<E, P> {
+impl<E: Residual<N>, N: Gradients, const P: usize> PhaseEquilibrium<E, P, N>
+where
+    DefaultAllocator: Allocator<N>,
+{
     pub(super) fn update_pressure(
         mut self,
         temperature: Temperature,
@@ -189,7 +192,7 @@ impl<E: Residual, const P: usize> PhaseEquilibrium<E, P> {
     pub(super) fn update_moles(
         &mut self,
         pressure: Pressure,
-        moles: [&Moles<DVector<f64>>; P],
+        moles: [&Moles<OVector<f64, N>>; P],
     ) -> FeosResult<()> {
         for (i, s) in self.0.iter_mut().enumerate() {
             *s = State::new_npt(

crates/feos-core/src/phase_equilibria/stability_analysis.rs

@@ -3,7 +3,9 @@ use crate::equation_of_state::Residual;
 use crate::errors::{FeosError, FeosResult};
 use crate::state::{Contributions, DensityInitialization, State};
 use crate::{ReferenceSystem, SolverOptions, Verbosity};
-use nalgebra::{DMatrix, DVector};
+use nalgebra::allocator::Allocator;
+use nalgebra::{DefaultAllocator, OMatrix, OVector, U1};
+use num_dual::Gradients;
 use num_dual::linalg::LU;
 use num_dual::linalg::smallest_ev;
 use quantity::Moles;
@@ -16,7 +18,10 @@ const MINIMIZE_KMAX: usize = 100;
 const ZERO_TPD: f64 = -1E-08;
 
 /// # Stability analysis
-impl<E: Residual> State<E> {
+impl<E: Residual<N>, N: Gradients> State<E, N>
+where
+    DefaultAllocator: Allocator<N> + Allocator<N, N>,
+{
     /// Determine if the state is stable, i.e. if a phase split should
     /// occur or not.
     pub fn is_stable(&self, options: SolverOptions) -> FeosResult<bool> {
@@ -26,7 +31,7 @@ impl<E: Residual> State<E> {
     /// Perform a stability analysis. The result is a list of [State]s with
     /// negative tangent plane distance (i.e. lower Gibbs energy) that can be
     /// used as initial estimates for a phase equilibrium calculation.
-    pub fn stability_analysis(&self, options: SolverOptions) -> FeosResult<Vec<State<E>>> {
+    pub fn stability_analysis(&self, options: SolverOptions) -> FeosResult<Vec<State<E, N>>> {
         let mut result = Vec::new();
         for i_trial in 0..self.eos.components() + 1 {
             let phase = if i_trial == self.eos.components() {
@@ -59,8 +64,9 @@ impl<E: Residual> State<E> {
         Ok(result)
     }
 
-    fn define_trial_state(&self, dominant_component: usize) -> FeosResult<State<E>> {
+    fn define_trial_state(&self, dominant_component: usize) -> FeosResult<State<E, N>> {
         let x_feed = &self.molefracs;
+        let (n, _) = x_feed.shape_generic();
 
         let (x_trial, phase) = if dominant_component == self.eos.components() {
             // try an ideal vapor phase
@@ -70,7 +76,7 @@ impl<E: Residual> State<E> {
             // try each component as nearly pure phase
             let factor = (1.0 - X_DOMINANT) / (x_feed.sum() - x_feed[dominant_component]);
             (
-                DVector::from_fn(self.eos.components(), |i, _| {
+                OVector::from_fn_generic(n, U1, |i, _| {
                     if i == dominant_component {
                         X_DOMINANT
                     } else {
@@ -92,7 +98,7 @@ impl<E: Residual> State<E> {
 
     fn minimize_tpd(
         &self,
-        trial: &mut State<E>,
+        trial: &mut State<E, N>,
         options: SolverOptions,
     ) -> FeosResult<(Option<f64>, usize)> {
         let (max_iter, tol, verbosity) = options.unwrap_or(MINIMIZE_KMAX, MINIMIZE_TOL);
@@ -154,9 +160,10 @@ impl<E: Residual> State<E> {
         Err(FeosError::NotConverged(String::from("stability analysis")))
     }
 
-    fn stability_newton_step(&mut self, di: &DVector<f64>, tpd: &mut f64) -> FeosResult<f64> {
+    fn stability_newton_step(&mut self, di: &OVector<f64, N>, tpd: &mut f64) -> FeosResult<f64> {
         // save old values
         let tpd_old = *tpd;
+        let (n, _) = di.shape_generic();
 
         // calculate residual and ideal hesse matrix
         let mut hesse = (self.dln_phi_dnj() * Moles::from_reduced(1.0)).into_value();
@@ -166,7 +173,7 @@ impl<E: Residual> State<E> {
         let sq_y = y.map(f64::sqrt);
         let gradient = (&ln_y + &lnphi - di).component_mul(&sq_y);
 
-        let hesse_ig = DMatrix::identity(self.eos.components(), self.eos.components());
+        let hesse_ig = OMatrix::identity_generic(n, n);
         for i in 0..self.eos.components() {
             hesse.column_mut(i).component_mul_assign(&(sq_y[i] * &sq_y));
             if y[i] > f64::EPSILON {
@@ -181,7 +188,7 @@ impl<E: Residual> State<E> {
         // ! (3) objective function (tpd) does not descent
         // !-----------------------------------------------------------------------------
         let mut adjust_hessian = true;
-        let mut hessian: DMatrix<f64>;
+        let mut hessian: OMatrix<f64, N, N>;
         let mut eta_h = 1.0;
 
         while adjust_hessian {