feos_core/phase_equilibria/
stability_analysis.rs

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

feos_core/phase_equilibria/
stability_analysis.rs
