feos_dft/interface/
surface_tension_diagram.rs

1use super::PlanarInterface;
2use crate::functional::{HelmholtzEnergyFunctional, DFT};
3use crate::solver::DFTSolver;
4use feos_core::{PhaseEquilibrium, ReferenceSystem, StateVec};
5use ndarray::{Array1, Array2};
6use quantity::{Length, Moles, SurfaceTension, Temperature};
7
8const DEFAULT_GRID_POINTS: usize = 2048;
9
10/// Container structure for the efficient calculation of surface tension diagrams.
11pub struct SurfaceTensionDiagram<F: HelmholtzEnergyFunctional> {
12    pub profiles: Vec<PlanarInterface<F>>,
13}
14
15// #[expect(clippy::ptr_arg)]
16impl<F: HelmholtzEnergyFunctional> SurfaceTensionDiagram<F> {
17    pub fn new(
18        dia: &[PhaseEquilibrium<DFT<F>, 2>],
19        init_densities: Option<bool>,
20        n_grid: Option<usize>,
21        l_grid: Option<Length>,
22        critical_temperature: Option<Temperature>,
23        fix_equimolar_surface: Option<bool>,
24        solver: Option<&DFTSolver>,
25    ) -> Self {
26        let n_grid = n_grid.unwrap_or(DEFAULT_GRID_POINTS);
27        let mut profiles: Vec<PlanarInterface<F>> = Vec::with_capacity(dia.len());
28        for vle in dia.iter() {
29            // check for a critical point
30            let profile = if PhaseEquilibrium::is_trivial_solution(vle.vapor(), vle.liquid()) {
31                Ok(PlanarInterface::from_tanh(
32                    vle,
33                    10,
34                    Length::from_reduced(100.0),
35                    Temperature::from_reduced(500.0),
36                    fix_equimolar_surface.unwrap_or(false),
37                ))
38            } else {
39                // initialize with pDGT for single segments and tanh for mixtures and segment DFT
40                if vle.vapor().eos.component_index().len() == 1 {
41                    PlanarInterface::from_pdgt(vle, n_grid, false)
42                } else {
43                    Ok(PlanarInterface::from_tanh(
44                        vle,
45                        n_grid,
46                        l_grid.unwrap_or(Length::from_reduced(100.0)),
47                        critical_temperature.unwrap_or(Temperature::from_reduced(500.0)),
48                        fix_equimolar_surface.unwrap_or(false),
49                    ))
50                }
51                .map(|mut profile| {
52                    if let Some(init) = profiles.last() {
53                        if init.profile.density.shape() == profile.profile.density.shape() {
54                            if let Some(scale) = init_densities {
55                                profile.set_density_inplace(&init.profile.density, scale)
56                            }
57                        }
58                    }
59                    profile
60                })
61            }
62            .and_then(|profile| profile.solve(solver));
63            if let Ok(profile) = profile {
64                profiles.push(profile);
65            }
66        }
67        Self { profiles }
68    }
69
70    pub fn vapor(&self) -> StateVec<'_, DFT<F>> {
71        self.profiles.iter().map(|p| p.vle.vapor()).collect()
72    }
73
74    pub fn liquid(&self) -> StateVec<'_, DFT<F>> {
75        self.profiles.iter().map(|p| p.vle.liquid()).collect()
76    }
77
78    pub fn surface_tension(&mut self) -> SurfaceTension<Array1<f64>> {
79        SurfaceTension::from_shape_fn(self.profiles.len(), |i| {
80            self.profiles[i].surface_tension.unwrap()
81        })
82    }
83
84    pub fn relative_adsorption(&self) -> Vec<Moles<Array2<f64>>> {
85        self.profiles
86            .iter()
87            .map(|planar_interf| planar_interf.relative_adsorption())
88            .collect()
89    }
90
91    pub fn interfacial_enrichment(&self) -> Vec<Array1<f64>> {
92        self.profiles
93            .iter()
94            .map(|planar_interf| planar_interf.interfacial_enrichment())
95            .collect()
96    }
97
98    pub fn interfacial_thickness(&self) -> Length<Array1<f64>> {
99        self.profiles
100            .iter()
101            .map(|planar_interf| planar_interf.interfacial_thickness().unwrap())
102            .collect()
103    }
104}
feos_dft/interface/surface_tension_diagram.rs

feos_dft/interface/
surface_tension_diagram.rs
