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

feos_core/phase_equilibria/
stability_analysis.rs
