Skip to content

Commit 8213aad

Browse files
prehnerprehner
committed
Extend tp-flash to static arrays and AD
1 parent 711b625 commit 8213aad

File tree

7 files changed

+207
-54
lines changed

7 files changed

+207
-54
lines changed

‎.github/workflows/test.yml‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,9 +2,9 @@ name: Test
22

33
on:
44
push:
5-
branches: [main]
5+
branches: [main, development]
66
pull_request:
7-
branches: [main]
7+
branches: [main, development]
88

99
env:
1010
CARGO_TERM_COLOR: always

‎.github/workflows/wheels.yml‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,9 @@
11
name: Build Wheels
22
on:
33
push:
4-
branches: [main]
4+
branches: [main, development]
55
pull_request:
6-
branches: [main]
6+
branches: [main, development]
77
jobs:
88
linux:
99
runs-on: ubuntu-latest

‎CHANGELOG.md‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
55
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
66

77
## [Breaking]
8+
### Added
9+
- 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)
10+
811
### Packaging
912
- Updated `quantity` dependency to 0.13 and removed the `typenum` dependency. [#323](https://github.com/feos-org/feos/pull/323)
1013

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

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -3,8 +3,8 @@ use crate::errors::{FeosError, FeosResult};
33
use crate::state::{DensityInitialization, State};
44
use crate::{Contributions, ReferenceSystem};
55
use nalgebra::allocator::Allocator;
6-
use nalgebra::{DVector, DefaultAllocator, Dim, Dyn, OVector};
7-
use num_dual::{DualNum, DualStruct};
6+
use nalgebra::{DefaultAllocator, Dim, Dyn, OVector};
7+
use num_dual::{DualNum, DualStruct, Gradients};
88
use quantity::{Energy, Moles, Pressure, RGAS, Temperature};
99
use std::fmt;
1010
use std::fmt::Write;
@@ -168,7 +168,10 @@ where
168168
}
169169
}
170170

171-
impl<E: Residual, const P: usize> PhaseEquilibrium<E, P> {
171+
impl<E: Residual<N>, N: Gradients, const P: usize> PhaseEquilibrium<E, P, N>
172+
where
173+
DefaultAllocator: Allocator<N>,
174+
{
172175
pub(super) fn update_pressure(
173176
mut self,
174177
temperature: Temperature,
@@ -189,7 +192,7 @@ impl<E: Residual, const P: usize> PhaseEquilibrium<E, P> {
189192
pub(super) fn update_moles(
190193
&mut self,
191194
pressure: Pressure,
192-
moles: [&Moles<DVector<f64>>; P],
195+
moles: [&Moles<OVector<f64, N>>; P],
193196
) -> FeosResult<()> {
194197
for (i, s) in self.0.iter_mut().enumerate() {
195198
*s = State::new_npt(

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

Lines changed: 16 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,9 @@ use crate::equation_of_state::Residual;
33
use crate::errors::{FeosError, FeosResult};
44
use crate::state::{Contributions, DensityInitialization, State};
55
use crate::{ReferenceSystem, SolverOptions, Verbosity};
6-
use nalgebra::{DMatrix, DVector};
6+
use nalgebra::allocator::Allocator;
7+
use nalgebra::{DefaultAllocator, OMatrix, OVector, U1};
8+
use num_dual::Gradients;
79
use num_dual::linalg::LU;
810
use num_dual::linalg::smallest_ev;
911
use quantity::Moles;
@@ -16,7 +18,10 @@ const MINIMIZE_KMAX: usize = 100;
1618
const ZERO_TPD: f64 = -1E-08;
1719

1820
/// # Stability analysis
19-
impl<E: Residual> State<E> {
21+
impl<E: Residual<N>, N: Gradients> State<E, N>
22+
where
23+
DefaultAllocator: Allocator<N> + Allocator<N, N>,
24+
{
2025
/// Determine if the state is stable, i.e. if a phase split should
2126
/// occur or not.
2227
pub fn is_stable(&self, options: SolverOptions) -> FeosResult<bool> {
@@ -26,7 +31,7 @@ impl<E: Residual> State<E> {
2631
/// Perform a stability analysis. The result is a list of [State]s with
2732
/// negative tangent plane distance (i.e. lower Gibbs energy) that can be
2833
/// used as initial estimates for a phase equilibrium calculation.
29-
pub fn stability_analysis(&self, options: SolverOptions) -> FeosResult<Vec<State<E>>> {
34+
pub fn stability_analysis(&self, options: SolverOptions) -> FeosResult<Vec<State<E, N>>> {
3035
let mut result = Vec::new();
3136
for i_trial in 0..self.eos.components() + 1 {
3237
let phase = if i_trial == self.eos.components() {
@@ -59,8 +64,9 @@ impl<E: Residual> State<E> {
5964
Ok(result)
6065
}
6166

62-
fn define_trial_state(&self, dominant_component: usize) -> FeosResult<State<E>> {
67+
fn define_trial_state(&self, dominant_component: usize) -> FeosResult<State<E, N>> {
6368
let x_feed = &self.molefracs;
69+
let (n, _) = x_feed.shape_generic();
6470

6571
let (x_trial, phase) = if dominant_component == self.eos.components() {
6672
// try an ideal vapor phase
@@ -70,7 +76,7 @@ impl<E: Residual> State<E> {
7076
// try each component as nearly pure phase
7177
let factor = (1.0 - X_DOMINANT) / (x_feed.sum() - x_feed[dominant_component]);
7278
(
73-
DVector::from_fn(self.eos.components(), |i, _| {
79+
OVector::from_fn_generic(n, U1, |i, _| {
7480
if i == dominant_component {
7581
X_DOMINANT
7682
} else {
@@ -92,7 +98,7 @@ impl<E: Residual> State<E> {
9298

9399
fn minimize_tpd(
94100
&self,
95-
trial: &mut State<E>,
101+
trial: &mut State<E, N>,
96102
options: SolverOptions,
97103
) -> FeosResult<(Option<f64>, usize)> {
98104
let (max_iter, tol, verbosity) = options.unwrap_or(MINIMIZE_KMAX, MINIMIZE_TOL);
@@ -154,9 +160,10 @@ impl<E: Residual> State<E> {
154160
Err(FeosError::NotConverged(String::from("stability analysis")))
155161
}
156162

157-
fn stability_newton_step(&mut self, di: &DVector<f64>, tpd: &mut f64) -> FeosResult<f64> {
163+
fn stability_newton_step(&mut self, di: &OVector<f64, N>, tpd: &mut f64) -> FeosResult<f64> {
158164
// save old values
159165
let tpd_old = *tpd;
166+
let (n, _) = di.shape_generic();
160167

161168
// calculate residual and ideal hesse matrix
162169
let mut hesse = (self.dln_phi_dnj() * Moles::from_reduced(1.0)).into_value();
@@ -166,7 +173,7 @@ impl<E: Residual> State<E> {
166173
let sq_y = y.map(f64::sqrt);
167174
let gradient = (&ln_y + &lnphi - di).component_mul(&sq_y);
168175

169-
let hesse_ig = DMatrix::identity(self.eos.components(), self.eos.components());
176+
let hesse_ig = OMatrix::identity_generic(n, n);
170177
for i in 0..self.eos.components() {
171178
hesse.column_mut(i).component_mul_assign(&(sq_y[i] * &sq_y));
172179
if y[i] > f64::EPSILON {
@@ -181,7 +188,7 @@ impl<E: Residual> State<E> {
181188
// ! (3) objective function (tpd) does not descent
182189
// !-----------------------------------------------------------------------------
183190
let mut adjust_hessian = true;
184-
let mut hessian: DMatrix<f64>;
191+
let mut hessian: OMatrix<f64, N, N>;
185192
let mut eta_h = 1.0;
186193

187194
while adjust_hessian {

0 commit comments

Comments
 (0)