Skip to content

Commit 3d57b71

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

4 files changed

Lines changed: 191 additions & 50 deletions

File tree

‎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 {

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

Lines changed: 111 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -2,16 +2,20 @@ use super::PhaseEquilibrium;
22
use crate::equation_of_state::Residual;
33
use crate::errors::{FeosError, FeosResult};
44
use crate::state::{Contributions, State};
5-
use crate::{SolverOptions, Verbosity};
6-
use nalgebra::{DVector, Matrix3, Matrix4xX};
7-
use num_dual::{Dual, DualNum, first_derivative};
8-
use quantity::{Dimensionless, Moles, Pressure, Temperature};
5+
use crate::{ReferenceSystem, SolverOptions, Verbosity};
6+
use nalgebra::allocator::Allocator;
7+
use nalgebra::{DefaultAllocator, Dim, Matrix3, OVector, SVector, U1, U2, vector};
8+
use num_dual::{Dual, DualNum, DualStruct, Gradients, first_derivative, implicit_derivative_sp};
9+
use quantity::{Dimensionless, MOL, MolarVolume, Moles, Pressure, Quantity, Temperature};
910

1011
const MAX_ITER_TP: usize = 400;
1112
const TOL_TP: f64 = 1e-8;
1213

1314
/// # Flash calculations
14-
impl<E: Residual> PhaseEquilibrium<E, 2> {
15+
impl<E: Residual<N>, N: Gradients> PhaseEquilibrium<E, 2, N>
16+
where
17+
DefaultAllocator: Allocator<N> + Allocator<N, N>,
18+
{
1519
/// Perform a Tp-flash calculation. If no initial values are
1620
/// given, the solution is initialized using a stability analysis.
1721
///
@@ -21,8 +25,8 @@ impl<E: Residual> PhaseEquilibrium<E, 2> {
2125
eos: &E,
2226
temperature: Temperature,
2327
pressure: Pressure,
24-
feed: &Moles<DVector<f64>>,
25-
initial_state: Option<&PhaseEquilibrium<E, 2>>,
28+
feed: &Moles<OVector<f64, N>>,
29+
initial_state: Option<&PhaseEquilibrium<E, 2, N>>,
2630
options: SolverOptions,
2731
non_volatile_components: Option<Vec<usize>>,
2832
) -> FeosResult<Self> {
@@ -34,8 +38,69 @@ impl<E: Residual> PhaseEquilibrium<E, 2> {
3438
}
3539
}
3640

41+
impl<E: Residual<U2, D>, D: DualNum<f64> + Copy> PhaseEquilibrium<E, 2, U2, D> {
42+
/// Perform a Tp-flash calculation. If no initial values are
43+
/// given, the solution is initialized using a stability analysis.
44+
///
45+
/// The algorithm can be use to calculate phase equilibria of systems
46+
/// containing non-volatile components (e.g. ions).
47+
pub fn tp_flash_binary(
48+
eos: &E,
49+
temperature: Temperature<D>,
50+
pressure: Pressure<D>,
51+
feed: &Moles<SVector<D, 2>>,
52+
options: SolverOptions,
53+
) -> FeosResult<Self> {
54+
let z = feed.get(0).convert_into(feed.get(0) + feed.get(1));
55+
let total_moles = feed.sum();
56+
let moles = vector![z.re(), 1.0 - z.re()] * MOL;
57+
let vle_re = State::new_npt(&eos.re(), temperature.re(), pressure.re(), &moles, None)?
58+
.tp_flash(None, options, None)?;
59+
60+
// implicit differentiation
61+
let t = temperature.into_reduced();
62+
let p = pressure.into_reduced();
63+
let v_l = vle_re.liquid().density.into_reduced().recip();
64+
let v_v = vle_re.vapor().density.into_reduced().recip();
65+
let x = vle_re.liquid().molefracs[0];
66+
let y = vle_re.vapor().molefracs[0];
67+
let x = implicit_derivative_sp(
68+
|x, args: &SVector<_, _>| {
69+
let [[v_l, v_v, x, y]] = x.data.0;
70+
let [[t, p, z]] = args.data.0;
71+
let beta = (z - x) / (y - x);
72+
let molefracs_l = vector![x, -x + 1.0];
73+
let molefracs_v = vector![y, -y + 1.0];
74+
let eos = eos.lift();
75+
let a_l_res = eos.residual_molar_helmholtz_energy(t, v_l, &molefracs_l);
76+
let a_l_ig = (x * (x / v_l).ln() - (x - 1.0) * ((-x + 1.0) / v_l).ln() - 1.0) * t;
77+
let a_v_res = eos.residual_molar_helmholtz_energy(t, v_v, &molefracs_v);
78+
let a_v_ig = (y * (y / v_v).ln() - (y - 1.0) * ((-y + 1.0) / v_v).ln() - 1.0) * t;
79+
let a_res = (a_v_res + a_v_ig) * beta - (a_l_res + a_l_ig) * (beta - 1.0);
80+
let v = v_v * beta - v_l * (beta - 1.0);
81+
a_res + v * p
82+
},
83+
SVector::from([v_l, v_v, x, y]),
84+
&SVector::from([t, p, z]),
85+
);
86+
let [[v_l, v_v, x, y]] = x.data.0;
87+
let beta = (z - x) / (y - x);
88+
let volume = MolarVolume::from_reduced(v_l * (-beta + 1.0)) * total_moles;
89+
let moles =
90+
Quantity::new(vector![x, -x + 1.0] * (-beta + 1.0) * total_moles.convert_into(MOL));
91+
let liquid = State::new_nvt(eos, temperature, volume, &moles)?;
92+
let volume = MolarVolume::from_reduced(v_v * beta) * total_moles;
93+
let moles = Quantity::new(vector![y, -y + 1.0] * beta * total_moles.convert_into(MOL));
94+
let vapor = State::new_nvt(eos, temperature, volume, &moles)?;
95+
Ok(Self([vapor, liquid]))
96+
}
97+
}
98+
3799
/// # Flash calculations
38-
impl<E: Residual> State<E> {
100+
impl<E: Residual<N>, N: Gradients> State<E, N>
101+
where
102+
DefaultAllocator: Allocator<N> + Allocator<N, N>,
103+
{
39104
/// Perform a Tp-flash calculation using the [State] as feed.
40105
/// If no initial values are given, the solution is initialized
41106
/// using a stability analysis.
@@ -44,10 +109,10 @@ impl<E: Residual> State<E> {
44109
/// containing non-volatile components (e.g. ions).
45110
pub fn tp_flash(
46111
&self,
47-
initial_state: Option<&PhaseEquilibrium<E, 2>>,
112+
initial_state: Option<&PhaseEquilibrium<E, 2, N>>,
48113
options: SolverOptions,
49114
non_volatile_components: Option<Vec<usize>>,
50-
) -> FeosResult<PhaseEquilibrium<E, 2>> {
115+
) -> FeosResult<PhaseEquilibrium<E, 2, N>> {
51116
// initialization
52117
if let Some(init) = initial_state {
53118
let vle = self.tp_flash_(
@@ -76,10 +141,10 @@ impl<E: Residual> State<E> {
76141

77142
pub fn tp_flash_(
78143
&self,
79-
mut new_vle_state: PhaseEquilibrium<E, 2>,
144+
mut new_vle_state: PhaseEquilibrium<E, 2, N>,
80145
options: SolverOptions,
81146
non_volatile_components: Option<Vec<usize>>,
82-
) -> FeosResult<PhaseEquilibrium<E, 2>> {
147+
) -> FeosResult<PhaseEquilibrium<E, 2, N>> {
83148
// set options
84149
let (max_iter, tol, verbosity) = options.unwrap_or(MAX_ITER_TP, TOL_TP);
85150

@@ -92,8 +157,8 @@ impl<E: Residual> State<E> {
92157
verbosity,
93158
" {:4} | | {:10.8?} | {:10.8?}",
94159
0,
95-
new_vle_state.vapor().molefracs.data.as_vec(),
96-
new_vle_state.liquid().molefracs.data.as_vec(),
160+
new_vle_state.vapor().molefracs.as_slice(),
161+
new_vle_state.liquid().molefracs.as_slice(),
97162
);
98163

99164
let mut iter = 0;
@@ -169,7 +234,7 @@ impl<E: Residual> State<E> {
169234
Ok(new_vle_state)
170235
}
171236

172-
fn tangent_plane_distance(&self, trial_state: &State<E>) -> f64 {
237+
fn tangent_plane_distance(&self, trial_state: &State<E, N>) -> f64 {
173238
let ln_phi_z = self.ln_phi();
174239
let ln_phi_w = trial_state.ln_phi();
175240
let z = &self.molefracs;
@@ -178,20 +243,23 @@ impl<E: Residual> State<E> {
178243
}
179244
}
180245

181-
impl<E: Residual> PhaseEquilibrium<E, 2> {
246+
impl<E: Residual<N>, N: Gradients> PhaseEquilibrium<E, 2, N>
247+
where
248+
DefaultAllocator: Allocator<N> + Allocator<N, N>,
249+
{
182250
fn accelerated_successive_substitution(
183251
&mut self,
184-
feed_state: &State<E>,
252+
feed_state: &State<E, N>,
185253
iter: &mut usize,
186254
max_iter: usize,
187255
tol: f64,
188256
verbosity: Verbosity,
189257
non_volatile_components: &Option<Vec<usize>>,
190258
) -> FeosResult<()> {
259+
let (n, _) = feed_state.molefracs.shape_generic();
191260
for _ in 0..max_iter {
192261
// do 5 successive substitution steps and check for convergence
193-
let mut k_vec = Matrix4xX::zeros(self.vapor().eos.components());
194-
// let mut k_vec = Array::zeros((4, self.vapor().eos.components()));
262+
let mut k_vec = std::array::repeat(OVector::zeros_generic(n, U1));
195263
if self.successive_substitution(
196264
feed_state,
197265
5,
@@ -213,16 +281,19 @@ impl<E: Residual> PhaseEquilibrium<E, 2> {
213281
let gibbs = self.total_gibbs_energy();
214282

215283
// extrapolate K values
216-
let delta_vec = k_vec.rows_range(1..) - k_vec.rows_range(..3);
217-
let delta = Matrix3::from_fn(|i, j| delta_vec.row(i).dot(&delta_vec.row(j)));
284+
let delta_vec = [
285+
&k_vec[1] - &k_vec[0],
286+
&k_vec[2] - &k_vec[1],
287+
&k_vec[3] - &k_vec[2],
288+
];
289+
let delta = Matrix3::from_fn(|i, j| delta_vec[i].dot(&delta_vec[j]));
218290
let d = delta[(0, 1)] * delta[(0, 1)] - delta[(0, 0)] * delta[(1, 1)];
219291
let a = (delta[(0, 2)] * delta[(0, 1)] - delta[(1, 2)] * delta[(0, 0)]) / d;
220292
let b = (delta[(1, 2)] * delta[(0, 1)] - delta[(0, 2)] * delta[(1, 1)]) / d;
221293

222-
let mut k = (k_vec.row(3)
223-
+ ((b * delta_vec.row(1) + (a + b) * delta_vec.row(2)) / (1.0 - a - b)))
224-
.map(f64::exp)
225-
.transpose();
294+
let mut k = (&k_vec[3]
295+
+ ((b * &delta_vec[1] + (a + b) * &delta_vec[2]) / (1.0 - a - b)))
296+
.map(f64::exp);
226297

227298
// Set k = 0 for non-volatile components
228299
if let Some(nvc) = non_volatile_components.as_ref() {
@@ -245,10 +316,10 @@ impl<E: Residual> PhaseEquilibrium<E, 2> {
245316
#[expect(clippy::too_many_arguments)]
246317
fn successive_substitution(
247318
&mut self,
248-
feed_state: &State<E>,
319+
feed_state: &State<E, N>,
249320
iterations: usize,
250321
iter: &mut usize,
251-
k_vec: &mut Option<&mut Matrix4xX<f64>>,
322+
k_vec: &mut Option<&mut [OVector<f64, N>; 4]>,
252323
abs_tol: f64,
253324
verbosity: Verbosity,
254325
non_volatile_components: &Option<Vec<usize>>,
@@ -278,8 +349,8 @@ impl<E: Residual> PhaseEquilibrium<E, 2> {
278349
" {:4} | {:14.8e} | {:.8?} | {:.8?}",
279350
iter,
280351
res,
281-
self.vapor().molefracs.data.as_vec(),
282-
self.liquid().molefracs.data.as_vec(),
352+
self.vapor().molefracs.as_slice(),
353+
self.liquid().molefracs.as_slice(),
283354
);
284355
if res < abs_tol {
285356
return Ok(true);
@@ -289,16 +360,13 @@ impl<E: Residual> PhaseEquilibrium<E, 2> {
289360
if let Some(k_vec) = k_vec
290361
&& i >= iterations - 3
291362
{
292-
k_vec.set_row(
293-
i + 3 - iterations,
294-
&k.map(|ki| if ki > 0.0 { ki.ln() } else { 0.0 }).transpose(),
295-
);
363+
k_vec[i + 3 - iterations] = k.map(|ki| if ki > 0.0 { ki.ln() } else { 0.0 });
296364
}
297365
}
298366
Ok(false)
299367
}
300368

301-
fn update_states(&mut self, feed_state: &State<E>, k: &DVector<f64>) -> FeosResult<()> {
369+
fn update_states(&mut self, feed_state: &State<E, N>, k: &OVector<f64, N>) -> FeosResult<()> {
302370
// calculate vapor phase fraction using Rachford-Rice algorithm
303371
let mut beta = self.vapor_phase_fraction();
304372
beta = rachford_rice(&feed_state.molefracs, k, Some(beta))?;
@@ -314,7 +382,7 @@ impl<E: Residual> PhaseEquilibrium<E, 2> {
314382
Ok(())
315383
}
316384

317-
fn vle_init_stability(feed_state: &State<E>) -> FeosResult<(Self, Option<Self>)> {
385+
fn vle_init_stability(feed_state: &State<E, N>) -> FeosResult<(Self, Option<Self>)> {
318386
let mut stable_states = feed_state.stability_analysis(SolverOptions::default())?;
319387
let state1 = stable_states.pop();
320388
let state2 = stable_states.pop();
@@ -331,7 +399,14 @@ impl<E: Residual> PhaseEquilibrium<E, 2> {
331399
}
332400
}
333401

334-
fn rachford_rice(feed: &DVector<f64>, k: &DVector<f64>, beta_in: Option<f64>) -> FeosResult<f64> {
402+
fn rachford_rice<N: Dim>(
403+
feed: &OVector<f64, N>,
404+
k: &OVector<f64, N>,
405+
beta_in: Option<f64>,
406+
) -> FeosResult<f64>
407+
where
408+
DefaultAllocator: Allocator<N>,
409+
{
335410
const MAX_ITER: usize = 10;
336411
const ABS_TOL: f64 = 1e-6;
337412

0 commit comments

Comments
 (0)