feos_core/phase_equilibria/
phase_diagram_pure.rs

1use super::PhaseEquilibrium;
2#[cfg(feature = "rayon")]
3use crate::ReferenceSystem;
4use crate::SolverOptions;
5use crate::equation_of_state::Residual;
6use crate::errors::FeosResult;
7use crate::state::{State, StateVec};
8#[cfg(feature = "rayon")]
9use ndarray::{Array1, ArrayView1, Axis};
10use quantity::Temperature;
11#[cfg(feature = "rayon")]
12use rayon::{ThreadPool, prelude::*};
13
14/// Pure component and binary mixture phase diagrams.
15#[derive(Clone)]
16pub struct PhaseDiagram<E, const N: usize> {
17    pub states: Vec<PhaseEquilibrium<E, N>>,
18}
19
20impl<E, const N: usize> PhaseDiagram<E, N> {
21    /// Create a phase diagram from a list of phase equilibria.
22    pub fn new(states: Vec<PhaseEquilibrium<E, N>>) -> Self {
23        Self { states }
24    }
25}
26
27impl<E: Residual> PhaseDiagram<E, 2> {
28    /// Calculate a phase diagram for a pure component.
29    pub fn pure(
30        eos: &E,
31        min_temperature: Temperature,
32        npoints: usize,
33        critical_temperature: Option<Temperature>,
34        options: SolverOptions,
35    ) -> FeosResult<Self> {
36        let mut states = Vec::with_capacity(npoints);
37
38        let sc = State::critical_point(
39            eos,
40            None,
41            critical_temperature,
42            None,
43            SolverOptions::default(),
44        )?;
45
46        let max_temperature = min_temperature
47            + (sc.temperature - min_temperature) * ((npoints - 2) as f64 / (npoints - 1) as f64);
48        let temperatures = Temperature::linspace(min_temperature, max_temperature, npoints - 1);
49
50        let mut vle = None;
51        for ti in &temperatures {
52            vle = PhaseEquilibrium::pure(eos, ti, vle.as_ref(), options).ok();
53            if let Some(vle) = vle.as_ref() {
54                states.push(vle.clone());
55            }
56        }
57        states.push(PhaseEquilibrium::from_states(sc.clone(), sc));
58
59        Ok(PhaseDiagram::new(states))
60    }
61
62    /// Return the vapor states of the diagram.
63    pub fn vapor(&self) -> StateVec<'_, E> {
64        self.states.iter().map(|s| s.vapor()).collect()
65    }
66
67    /// Return the liquid states of the diagram.
68    pub fn liquid(&self) -> StateVec<'_, E> {
69        self.states.iter().map(|s| s.liquid()).collect()
70    }
71}
72
73#[cfg(feature = "rayon")]
74impl<E: Residual> PhaseDiagram<E, 2> {
75    fn solve_temperatures(
76        eos: &E,
77        temperatures: ArrayView1<f64>,
78        options: SolverOptions,
79    ) -> FeosResult<Vec<PhaseEquilibrium<E, 2>>> {
80        let mut states = Vec::with_capacity(temperatures.len());
81        let mut vle = None;
82        for ti in temperatures {
83            vle =
84                PhaseEquilibrium::pure(eos, Temperature::from_reduced(*ti), vle.as_ref(), options)
85                    .ok();
86            if let Some(vle) = vle.as_ref() {
87                states.push(vle.clone());
88            }
89        }
90        Ok(states)
91    }
92
93    pub fn par_pure(
94        eos: &E,
95        min_temperature: Temperature,
96        npoints: usize,
97        chunksize: usize,
98        thread_pool: ThreadPool,
99        critical_temperature: Option<Temperature>,
100        options: SolverOptions,
101    ) -> FeosResult<Self>
102    where
103        E: Send + Sync,
104    {
105        let sc = State::critical_point(
106            eos,
107            None,
108            critical_temperature,
109            None,
110            SolverOptions::default(),
111        )?;
112
113        let max_temperature = min_temperature
114            + (sc.temperature - min_temperature) * ((npoints - 2) as f64 / (npoints - 1) as f64);
115        let temperatures = Array1::linspace(
116            min_temperature.to_reduced(),
117            max_temperature.to_reduced(),
118            npoints - 1,
119        );
120
121        let mut states: Vec<PhaseEquilibrium<E, 2>> = thread_pool.install(|| {
122            temperatures
123                .axis_chunks_iter(Axis(0), chunksize)
124                .into_par_iter()
125                .filter_map(|t| Self::solve_temperatures(eos, t, options).ok())
126                .flatten()
127                .collect()
128        });
129
130        states.push(PhaseEquilibrium::from_states(sc.clone(), sc));
131        Ok(PhaseDiagram::new(states))
132    }
133}
feos_core/phase_equilibria/phase_diagram_pure.rs

feos_core/phase_equilibria/
phase_diagram_pure.rs
