feos_core/phase_equilibria/
stability_analysis.rs

1use super::PhaseEquilibrium;
2use crate::equation_of_state::Residual;
3use crate::errors::{EosError, EosResult};
4use crate::state::{Contributions, DensityInitialization, State};
5use crate::{ReferenceSystem, SolverOptions, Verbosity};
6use ndarray::*;
7use num_dual::linalg::smallest_ev;
8use num_dual::linalg::LU;
9use quantity::Moles;
10use std::ops::MulAssign;
11
12const X_DOMINANT: f64 = 0.99;
13const MINIMIZE_TOL: f64 = 1E-06;
14const MIN_EIGENVAL: f64 = 1E-03;
15const ETA_STEP: f64 = 0.25;
16const MINIMIZE_KMAX: usize = 100;
17const ZERO_TPD: f64 = -1E-08;
18
19/// # Stability analysis
20impl<E: Residual> State<E> {
21    /// Determine if the state is stable, i.e. if a phase split should
22    /// occur or not.
23    pub fn is_stable(&self, options: SolverOptions) -> EosResult<bool> {
24        Ok(self.stability_analysis(options)?.is_empty())
25    }
26
27    /// Perform a stability analysis. The result is a list of [State]s with
28    /// negative tangent plane distance (i.e. lower Gibbs energy) that can be
29    /// used as initial estimates for a phase equilibrium calculation.
30    pub fn stability_analysis(&self, options: SolverOptions) -> EosResult<Vec<State<E>>> {
31        let mut result = Vec::new();
32        for i_trial in 0..self.eos.components() + 1 {
33            let phase = if i_trial == self.eos.components() {
34                "Vapor phase".to_string()
35            } else {
36                format!("Liquid phase {}", i_trial + 1)
37            };
38            if let Ok(mut trial_state) = self.define_trial_state(i_trial) {
39                let (tpd, i) = self.minimize_tpd(&mut trial_state, options)?;
40                let msg = if let Some(tpd) = tpd {
41                    if tpd < ZERO_TPD {
42                        if result
43                            .iter()
44                            .any(|s| PhaseEquilibrium::is_trivial_solution(s, &trial_state))
45                        {
46                            "Found already identified minimum"
47                        } else {
48                            result.push(trial_state);
49                            "Found candidate"
50                        }
51                    } else {
52                        "Found minimum > 0"
53                    }
54                } else {
55                    "Found trivial solution"
56                };
57                log_result!(options.verbosity, "{}: {} in {} step(s)\n", phase, msg, i);
58            }
59        }
60        Ok(result)
61    }
62
63    fn define_trial_state(&self, dominant_component: usize) -> EosResult<State<E>> {
64        let x_feed = &self.molefracs;
65
66        let (x_trial, phase) = if dominant_component == self.eos.components() {
67            // try an ideal vapor phase
68            let x_trial = self.ln_phi().mapv(f64::exp) * x_feed;
69            (&x_trial / x_trial.sum(), DensityInitialization::Vapor)
70        } else {
71            // try each component as nearly pure phase
72            let factor = (1.0 - X_DOMINANT) / (x_feed.sum() - x_feed[dominant_component]);
73            (
74                Array1::from_shape_fn(self.eos.components(), |i| {
75                    if i == dominant_component {
76                        X_DOMINANT
77                    } else {
78                        x_feed[i] * factor
79                    }
80                }),
81                DensityInitialization::Liquid,
82            )
83        };
84
85        State::new_npt(
86            &self.eos,
87            self.temperature,
88            self.pressure(Contributions::Total),
89            &Moles::from_reduced(x_trial),
90            phase,
91        )
92    }
93
94    fn minimize_tpd(
95        &self,
96        trial: &mut State<E>,
97        options: SolverOptions,
98    ) -> EosResult<(Option<f64>, usize)> {
99        let (max_iter, tol, verbosity) = options.unwrap_or(MINIMIZE_KMAX, MINIMIZE_TOL);
100        let mut newton = false;
101        let mut scaled_tol = tol;
102        let mut tpd = 1E10;
103        let di = self.molefracs.mapv(f64::ln) + self.ln_phi();
104
105        log_iter!(verbosity, " iter |    residual    |     tpd     | Newton");
106        log_iter!(verbosity, "{:-<46}", "");
107
108        for i in 1..=max_iter {
109            let error = if !newton {
110                // case: direct substitution
111                let y = (&di - &trial.ln_phi()).mapv(f64::exp);
112                let tpd_old = tpd;
113                tpd = 1.0 - y.sum();
114                let error = (&y / y.sum() - &trial.molefracs).mapv(f64::abs).sum();
115
116                *trial = State::new_npt(
117                    &trial.eos,
118                    trial.temperature,
119                    trial.pressure(Contributions::Total),
120                    &Moles::from_reduced(y),
121                    DensityInitialization::InitialDensity(trial.density),
122                )?;
123                if (i > 4 && error > scaled_tol) || (tpd > tpd_old + 1E-05 && i > 2) {
124                    newton = true; // switch to newton scheme
125                }
126                error
127            } else {
128                // case: newton step
129                trial.stability_newton_step(&di, &mut tpd)?
130            };
131            log_iter!(
132                verbosity,
133                " {:4} | {:14.8e} | {:11.8} | {}",
134                i,
135                error,
136                tpd,
137                newton
138            );
139            if PhaseEquilibrium::is_trivial_solution(self, &*trial) {
140                return Ok((None, i));
141            }
142            if tpd < -1E-02 {
143                scaled_tol = tol * 1E01
144            }
145            if tpd < -1E-01 {
146                scaled_tol = tol * 1E02
147            }
148            if tpd < -1E-01 && i > 5 {
149                scaled_tol = tol * 1E03
150            }
151            if error < scaled_tol {
152                return Ok((Some(tpd), i));
153            }
154        }
155        Err(EosError::NotConverged(String::from("stability analysis")))
156    }
157
158    fn stability_newton_step(&mut self, di: &Array1<f64>, tpd: &mut f64) -> EosResult<f64> {
159        // save old values
160        let tpd_old = *tpd;
161
162        // calculate residual and ideal hesse matrix
163        let mut hesse = (self.dln_phi_dnj() * Moles::from_reduced(1.0)).into_value();
164        let lnphi = self.ln_phi();
165        let y = self.moles.to_reduced();
166        let ln_y = Zip::from(&y).map_collect(|&y| if y > f64::EPSILON { y.ln() } else { 0.0 });
167        let sq_y = y.mapv(f64::sqrt);
168        let gradient = (&ln_y + &lnphi - di) * &sq_y;
169
170        let hesse_ig = Array2::eye(self.eos.components());
171        for i in 0..self.eos.components() {
172            hesse
173                .index_axis_mut(Axis(0), i)
174                .mul_assign(&(sq_y[i] * &sq_y));
175            if y[i] > f64::EPSILON {
176                hesse[[i, i]] += ln_y[i] + lnphi[i] - di[i];
177            }
178        }
179
180        // !-----------------------------------------------------------------------------
181        // ! use method of Murray, by adding a unity matrix to Hessian, if:
182        // ! (1) H is not positive definite
183        // ! (2) step size is too large
184        // ! (3) objective function (tpd) does not descent
185        // !-----------------------------------------------------------------------------
186        let mut adjust_hessian = true;
187        let mut hessian: Array2<f64>;
188        let mut eta_h = 1.0;
189
190        while adjust_hessian {
191            adjust_hessian = false;
192            hessian = &hesse + &(eta_h * &hesse_ig);
193
194            let (min_eigenval, _) = smallest_ev(hessian.clone());
195            if min_eigenval < MIN_EIGENVAL && eta_h < 20.0 {
196                eta_h += 2.0 * ETA_STEP;
197                adjust_hessian = true;
198                continue; // continue, because of Hessian-criterion (1): H not positive definite
199            }
200
201            // solve: hessian * delta_y = gradient
202            let delta_y = LU::new(hessian)?.solve(&gradient);
203            if delta_y
204                .iter()
205                .zip(y.iter())
206                .any(|(dy, y)| ((0.5 * dy).powi(2) / y).abs() > 5.0)
207            {
208                adjust_hessian = true;
209                eta_h += 2.0 * ETA_STEP;
210                continue; //  continue, because of Hessian-criterion (2): too large step-size
211            }
212
213            let y = (&sq_y - &(delta_y / 2.0)).mapv(|v| v.powi(2));
214            let ln_y = Zip::from(&y).map_collect(|&y| if y > f64::EPSILON { y.ln() } else { 0.0 });
215            *tpd = 1.0 + (&y * &(&ln_y + &lnphi - di - 1.0)).sum();
216            if *tpd > tpd_old + 0.0 * 1E-03 && eta_h < 30.0 {
217                eta_h += ETA_STEP;
218                adjust_hessian = true;
219                continue; // continue, because of Hessian-criterion (3): tpd does not descent
220            }
221
222            // accept step and update state
223            *self = State::new_npt(
224                &self.eos,
225                self.temperature,
226                self.pressure(Contributions::Total),
227                &Moles::from_reduced(y),
228                DensityInitialization::InitialDensity(self.density),
229            )?;
230        }
231        Ok(gradient.mapv(f64::abs).sum())
232    }
233}
feos_core/phase_equilibria/stability_analysis.rs

feos_core/phase_equilibria/
stability_analysis.rs
