feos_dft/interface/
surface_tension_diagram.rs

1use super::PlanarInterface;
2use crate::functional::HelmholtzEnergyFunctional;
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<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                        && init.profile.density.shape() == profile.profile.density.shape()
54                        && let Some(scale) = init_densities
55                    {
56                        profile.set_density_inplace(&init.profile.density, scale)
57                    }
58                    profile
59                })
60            }
61            .and_then(|profile| profile.solve(solver));
62            if let Ok(profile) = profile {
63                profiles.push(profile);
64            }
65        }
66        Self { profiles }
67    }
68
69    pub fn vapor(&self) -> StateVec<'_, F> {
70        self.profiles.iter().map(|p| p.vle.vapor()).collect()
71    }
72
73    pub fn liquid(&self) -> StateVec<'_, F> {
74        self.profiles.iter().map(|p| p.vle.liquid()).collect()
75    }
76
77    pub fn surface_tension(&mut self) -> SurfaceTension<Array1<f64>> {
78        SurfaceTension::from_shape_fn(self.profiles.len(), |i| {
79            self.profiles[i].surface_tension.unwrap()
80        })
81    }
82
83    pub fn relative_adsorption(&self) -> Vec<Moles<Array2<f64>>> {
84        self.profiles
85            .iter()
86            .map(|planar_interf| planar_interf.relative_adsorption())
87            .collect()
88    }
89
90    pub fn interfacial_enrichment(&self) -> Vec<Array1<f64>> {
91        self.profiles
92            .iter()
93            .map(|planar_interf| planar_interf.interfacial_enrichment())
94            .collect()
95    }
96
97    pub fn interfacial_thickness(&self) -> Length<Array1<f64>> {
98        self.profiles
99            .iter()
100            .map(|planar_interf| planar_interf.interfacial_thickness().unwrap())
101            .collect()
102    }
103}
feos_dft/interface/surface_tension_diagram.rs

feos_dft/interface/
surface_tension_diagram.rs
