Skip to content

Commit bd4280b

Browse files
prehnerprehner
authored
Improve the robustness of adsorption isotherm calculations (feos-org#50)
* Improve the robustness of adsorption isotherm calculations * Update notebooks
1 parent a318e44 commit bd4280b

File tree

8 files changed

+136
-209
lines changed

8 files changed

+136
-209
lines changed

‎examples/pcsaft/adsorption_3D.ipynb‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -192,7 +192,7 @@
192192
{
193193
"data": {
194194
"text/plain": [
195-
"<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7f69a8a18be0>"
195+
"<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7f8a755c5220>"
196196
]
197197
},
198198
"execution_count": 6,
@@ -294,14 +294,15 @@
294294
"name": "stdout",
295295
"output_type": "stream",
296296
"text": [
297-
"CPU times: user 17min 23s, sys: 2min 31s, total: 19min 55s\n",
298-
"Wall time: 19min 47s\n"
297+
"CPU times: user 17min 11s, sys: 2min 13s, total: 19min 25s\n",
298+
"Wall time: 19min 18s\n"
299299
]
300300
}
301301
],
302302
"source": [
303303
"%%time\n",
304-
"isotherm = Adsorption3D.adsorption_isotherm(func, temperature = T, pressure = (100*PASCAL, 50*BAR, 20), pore=pore, solver=solver)"
304+
"pressure = SIArray1.linspace(100*PASCAL, 50*BAR, 20)\n",
305+
"isotherm = Adsorption3D.adsorption_isotherm(func, temperature = T, pressure = pressure, pore=pore, solver=solver)"
305306
]
306307
},
307308
{

‎examples/pcsaft/adsorption_isotherms.ipynb‎

Lines changed: 31 additions & 30 deletions
Large diffs are not rendered by default.

‎examples/pcsaft/fea_adsorption.ipynb‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@
88
"source": [
99
"from feos.dft import ExternalPotential, HelmholtzEnergyFunctional, Adsorption1D, Geometry, Pore1D, State\n",
1010
"from feos.pcsaft import PcSaftParameters\n",
11-
"from feos.si import ANGSTROM, KELVIN, BAR, NAV, KILO, METER, MOL\n",
11+
"from feos.si import ANGSTROM, KELVIN, BAR, NAV, KILO, METER, MOL, SIArray1\n",
1212
"\n",
1313
"import pandas as pd\n",
1414
"import numpy as np\n",
@@ -94,8 +94,8 @@
9494
"text": [
9595
"6.84534814234455 3220\n",
9696
"[5.9595, 5.9594999999999985, 5.9595]\n",
97-
"CPU times: user 321 ms, sys: 232 µs, total: 321 ms\n",
98-
"Wall time: 320 ms\n"
97+
"CPU times: user 328 ms, sys: 3.58 ms, total: 331 ms\n",
98+
"Wall time: 331 ms\n"
9999
]
100100
}
101101
],
@@ -140,7 +140,7 @@
140140
"metadata": {},
141141
"outputs": [],
142142
"source": [
143-
"isotherm = Adsorption1D.adsorption_isotherm(func, 298.0 * KELVIN, (1.0e-3 * BAR, 5.0 * BAR, 51), pore)"
143+
"isotherm = Adsorption1D.adsorption_isotherm(func, 298.0 * KELVIN, SIArray1.linspace(1.0e-3 * BAR, 5.0 * BAR, 51), pore)"
144144
]
145145
},
146146
{
@@ -151,7 +151,7 @@
151151
{
152152
"data": {
153153
"text/plain": [
154-
"<matplotlib.legend.Legend at 0x7f68dce9abb0>"
154+
"<matplotlib.legend.Legend at 0x7f539de574c0>"
155155
]
156156
},
157157
"execution_count": 8,

‎examples/pcsaft/pore_geometry.ipynb‎

Lines changed: 8 additions & 8 deletions
Large diffs are not rendered by default.

‎feos-dft/CHANGELOG.md‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
66

77
## [Unreleased]
88

9+
## [0.3.1] - 2022-08-26
10+
### Changed
11+
- `Adsorption::equilibrium_isotherm`, `Adsorption::adsorption_isotherm`, and `Adsorption::desorption_isotherm` only accept `SIArray1`s as input for the pressure discretization. [#50](https://github.com/feos-org/feos/pull/50)
12+
- For the calculation of desorption isotherms and the filled branches of equilibrium isotherms, the density profiles are initialized using a liquid bulk density. [#50](https://github.com/feos-org/feos/pull/50)
13+
914
## [0.3.0] - 2022-07-13
1015
### Added
1116
- Added getters for the fields of `Pore1D` in Python. [#30](https://github.com/feos-org/feos-dft/pull/30)
@@ -16,7 +21,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
1621
- Changed interface of `PairCorrelationFunction` to facilitate the calculation of pair correlation functions in mixtures. [#29](https://github.com/feos-org/feos-dft/pull/29)
1722

1823
### Removed
19-
- Moved the implementation of fundamental measure theory to the new `feos-saft` crate. [#27](https://github.com/feos-org/feos/pull/27)
24+
- Moved the implementation of fundamental measure theory to the `feos` crate. [#27](https://github.com/feos-org/feos/pull/27)
2025

2126
## [0.2.0] - 2022-04-12
2227
### Added

‎feos-dft/Cargo.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
[package]
22
name = "feos-dft"
3-
version = "0.3.0"
3+
version = "0.3.1"
44
authors = ["Philipp Rehner <prehner@ethz.ch>", "Gernot Bauer <bauer@itt.uni-stuttgart.de>"]
55
edition = "2021"
66
license = "MIT OR Apache-2.0"

‎feos-dft/src/adsorption/mod.rs‎

Lines changed: 68 additions & 126 deletions
Original file line numberDiff line numberDiff line change
@@ -2,10 +2,12 @@
22
use super::functional::{HelmholtzEnergyFunctional, DFT};
33
use super::solver::DFTSolver;
44
use feos_core::{
5-
Contributions, EosError, EosResult, EosUnit, EquationOfState, SolverOptions, StateBuilder,
5+
Contributions, DensityInitialization, EosError, EosResult, EosUnit, EquationOfState,
6+
SolverOptions, State, StateBuilder,
67
};
7-
use ndarray::{arr1, Array1, Dimension, Ix1, Ix3, RemoveAxis};
8+
use ndarray::{Array1, Dimension, Ix1, Ix3, RemoveAxis};
89
use quantity::{QuantityArray1, QuantityArray2, QuantityScalar};
10+
use std::iter;
911
use std::rc::Rc;
1012

1113
mod external_potential;
@@ -17,93 +19,6 @@ pub use pore::{Pore1D, Pore3D, PoreProfile, PoreProfile1D, PoreProfile3D, PoreSp
1719
const MAX_ITER_ADSORPTION_EQUILIBRIUM: usize = 50;
1820
const TOL_ADSORPTION_EQUILIBRIUM: f64 = 1e-8;
1921

20-
/// Possible inputs for the pressure grid of adsorption isotherms.
21-
pub enum PressureSpecification<U> {
22-
/// Specify the minimal and maximal pressure, and the number of points
23-
Plim {
24-
p_min: QuantityScalar<U>,
25-
p_max: QuantityScalar<U>,
26-
points: usize,
27-
},
28-
/// Provide a custom array of pressure points.
29-
Pvec(QuantityArray1<U>),
30-
}
31-
32-
impl<U: EosUnit> PressureSpecification<U>
33-
where
34-
QuantityScalar<U>: std::fmt::Display,
35-
{
36-
fn p_min_max(&self) -> (QuantityScalar<U>, QuantityScalar<U>) {
37-
match self {
38-
Self::Plim {
39-
p_min,
40-
p_max,
41-
points: _,
42-
} => (*p_min, *p_max),
43-
Self::Pvec(pressure) => (pressure.get(0), pressure.get(pressure.len() - 1)),
44-
}
45-
}
46-
47-
fn to_vec(&self) -> EosResult<QuantityArray1<U>> {
48-
Ok(match self {
49-
Self::Plim {
50-
p_min,
51-
p_max,
52-
points,
53-
} => QuantityArray1::linspace(*p_min, *p_max, *points)?,
54-
Self::Pvec(pressure) => pressure.clone(),
55-
})
56-
}
57-
58-
fn equilibrium<
59-
D: Dimension + RemoveAxis + 'static,
60-
F: HelmholtzEnergyFunctional + FluidParameters,
61-
>(
62-
&self,
63-
equilibrium: &Adsorption<U, D, F>,
64-
) -> EosResult<(QuantityArray1<U>, QuantityArray1<U>)>
65-
where
66-
D::Larger: Dimension<Smaller = D>,
67-
D::Smaller: Dimension<Larger = D>,
68-
<D::Larger as Dimension>::Larger: Dimension<Smaller = D::Larger>,
69-
{
70-
let p_eq = equilibrium.pressure().get(0);
71-
match self {
72-
Self::Plim {
73-
p_min,
74-
p_max,
75-
points,
76-
} => Ok((
77-
QuantityArray1::linspace(*p_min, p_eq, points / 2)?,
78-
QuantityArray1::linspace(*p_max, p_eq, points - points / 2)?,
79-
)),
80-
Self::Pvec(pressure) => {
81-
let index = (0..pressure.len()).find(|&i| pressure.get(i) > p_eq);
82-
Ok(if let Some(index) = index {
83-
(
84-
QuantityArray1::from_shape_fn(index + 1, |i| {
85-
if i == index {
86-
p_eq
87-
} else {
88-
pressure.get(i)
89-
}
90-
}),
91-
QuantityArray1::from_shape_fn(pressure.len() - index + 1, |i| {
92-
if i == pressure.len() - index {
93-
p_eq
94-
} else {
95-
pressure.get(pressure.len() - i - 1)
96-
}
97-
}),
98-
)
99-
} else {
100-
(pressure.clone(), arr1(&[]) * U::reference_pressure())
101-
})
102-
}
103-
}
104-
}
105-
}
106-
10722
/// Container structure for the calculation of adsorption isotherms.
10823
pub struct Adsorption<U, D: Dimension, F> {
10924
components: usize,
@@ -143,28 +58,41 @@ where
14358
pub fn adsorption_isotherm<S: PoreSpecification<U, D>>(
14459
functional: &Rc<DFT<F>>,
14560
temperature: QuantityScalar<U>,
146-
pressure: &PressureSpecification<U>,
61+
pressure: &QuantityArray1<U>,
14762
pore: &S,
14863
molefracs: Option<&Array1<f64>>,
14964
solver: Option<&DFTSolver>,
15065
) -> EosResult<Adsorption<U, D, F>> {
151-
let pressure = pressure.to_vec()?;
152-
Self::isotherm(functional, temperature, &pressure, pore, molefracs, solver)
66+
Self::isotherm(
67+
functional,
68+
temperature,
69+
pressure,
70+
pore,
71+
molefracs,
72+
DensityInitialization::Vapor,
73+
solver,
74+
)
15375
}
15476

15577
/// Calculate an desorption isotherm (starting at high pressure)
15678
pub fn desorption_isotherm<S: PoreSpecification<U, D>>(
15779
functional: &Rc<DFT<F>>,
15880
temperature: QuantityScalar<U>,
159-
pressure: &PressureSpecification<U>,
81+
pressure: &QuantityArray1<U>,
16082
pore: &S,
16183
molefracs: Option<&Array1<f64>>,
16284
solver: Option<&DFTSolver>,
16385
) -> EosResult<Adsorption<U, D, F>> {
164-
let pressure = pressure.to_vec()?;
165-
let pressure =
166-
QuantityArray1::from_shape_fn(pressure.len(), |i| pressure.get(pressure.len() - i - 1));
167-
let isotherm = Self::isotherm(functional, temperature, &pressure, pore, molefracs, solver)?;
86+
let pressure = pressure.into_iter().rev().collect();
87+
let isotherm = Self::isotherm(
88+
functional,
89+
temperature,
90+
&pressure,
91+
pore,
92+
molefracs,
93+
DensityInitialization::Liquid,
94+
solver,
95+
)?;
16896
Ok(Adsorption::new(
16997
functional,
17098
pore,
@@ -176,12 +104,12 @@ where
176104
pub fn equilibrium_isotherm<S: PoreSpecification<U, D>>(
177105
functional: &Rc<DFT<F>>,
178106
temperature: QuantityScalar<U>,
179-
pressure: &PressureSpecification<U>,
107+
pressure: &QuantityArray1<U>,
180108
pore: &S,
181109
molefracs: Option<&Array1<f64>>,
182110
solver: Option<&DFTSolver>,
183111
) -> EosResult<Adsorption<U, D, F>> {
184-
let (p_min, p_max) = pressure.p_min_max();
112+
let (p_min, p_max) = (pressure.get(0), pressure.get(pressure.len() - 1));
185113
let equilibrium = Self::phase_equilibrium(
186114
functional,
187115
temperature,
@@ -193,20 +121,28 @@ where
193121
SolverOptions::default(),
194122
);
195123
if let Ok(equilibrium) = equilibrium {
196-
let pressure = pressure.equilibrium(&equilibrium)?;
197-
let adsorption = Self::isotherm(
124+
let p_eq = equilibrium.pressure().get(0);
125+
let p_ads = pressure
126+
.into_iter()
127+
.filter(|&p| p <= p_eq)
128+
.chain(iter::once(p_eq))
129+
.collect();
130+
let p_des = iter::once(p_eq)
131+
.chain(pressure.into_iter().filter(|&p| p > p_eq))
132+
.collect();
133+
let adsorption = Self::adsorption_isotherm(
198134
functional,
199135
temperature,
200-
&pressure.0,
136+
&p_ads,
201137
pore,
202138
molefracs,
203139
solver,
204140
)?
205141
.profiles;
206-
let desorption = Self::isotherm(
142+
let desorption = Self::desorption_isotherm(
207143
functional,
208144
temperature,
209-
&pressure.1,
145+
&p_des,
210146
pore,
211147
molefracs,
212148
solver,
@@ -215,7 +151,7 @@ where
215151
Ok(Adsorption {
216152
profiles: adsorption
217153
.into_iter()
218-
.chain(desorption.into_iter().rev())
154+
.chain(desorption.into_iter())
219155
.collect(),
220156
components: functional.components(),
221157
dimension: pore.dimension(),
@@ -258,28 +194,33 @@ where
258194
pressure: &QuantityArray1<U>,
259195
pore: &S,
260196
molefracs: Option<&Array1<f64>>,
197+
density_initialization: DensityInitialization<U>,
261198
solver: Option<&DFTSolver>,
262199
) -> EosResult<Adsorption<U, D, F>> {
263200
let moles =
264201
functional.validate_moles(molefracs.map(|x| x * U::reference_moles()).as_ref())?;
265202
let mut profiles: Vec<EosResult<PoreProfile<U, D, F>>> = Vec::with_capacity(pressure.len());
266203

267-
// Calculate the external potential once
268-
let mut bulk = StateBuilder::new(functional)
269-
.temperature(temperature)
270-
.pressure(pressure.get(0))
271-
.moles(&moles)
272-
.build()?;
204+
// On the first iteration, initialize the density profile according to the direction
205+
// and calculate the external potential once.
206+
let mut bulk = State::new_npt(
207+
functional,
208+
temperature,
209+
pressure.get(0),
210+
&moles,
211+
density_initialization,
212+
)?;
273213
if functional.components() > 1 && !bulk.is_stable(SolverOptions::default())? {
274-
bulk = bulk
275-
.tp_flash(None, SolverOptions::default(), None)?
276-
.vapor()
277-
.clone();
214+
let vle = bulk.tp_flash(None, SolverOptions::default(), None)?;
215+
bulk = match density_initialization {
216+
DensityInitialization::Liquid => vle.liquid().clone(),
217+
DensityInitialization::Vapor => vle.vapor().clone(),
218+
_ => unreachable!(),
219+
};
278220
}
279-
let external_potential = pore
280-
.initialize(&bulk, None, None)?
281-
.profile
282-
.external_potential;
221+
let profile = pore.initialize(&bulk, None, None)?.solve(solver)?.profile;
222+
let external_potential = Some(&profile.external_potential);
223+
let mut old_density = Some(&profile.density);
283224

284225
for i in 0..pressure.len() {
285226
let mut bulk = StateBuilder::new(functional)
@@ -293,15 +234,16 @@ where
293234
.vapor()
294235
.clone();
295236
}
296-
let old_density = if let Some(Ok(l)) = profiles.last() {
237+
238+
let p = pore.initialize(&bulk, old_density, external_potential)?;
239+
let p2 = pore.initialize(&bulk, None, external_potential)?;
240+
profiles.push(p.solve(solver).or_else(|_| p2.solve(solver)));
241+
242+
old_density = if let Some(Ok(l)) = profiles.last() {
297243
Some(&l.profile.density)
298244
} else {
299245
None
300246
};
301-
302-
let p = pore.initialize(&bulk, old_density, Some(&external_potential))?;
303-
let p2 = pore.initialize(&bulk, None, Some(&external_potential))?;
304-
profiles.push(p.solve(solver).or_else(|_| p2.solve(solver)));
305247
}
306248

307249
Ok(Adsorption::new(functional, pore, profiles))
@@ -344,11 +286,11 @@ where
344286
let nl = liquid.profile.bulk.density
345287
* (liquid.profile.moles() * liquid.profile.bulk.molar_volume(Contributions::Total))
346288
.sum();
347-
let f = |s: &mut PoreProfile<U, D, F>, n: QuantityScalar<U>| -> EosResult<_> {
289+
let f = |s: &PoreProfile<U, D, F>, n: QuantityScalar<U>| -> EosResult<_> {
348290
Ok(s.grand_potential.unwrap()
349291
+ s.profile.bulk.molar_gibbs_energy(Contributions::Total) * n)
350292
};
351-
let mut g = (f(&mut liquid, nl)? - f(&mut vapor, nv)?) / (nl - nv);
293+
let mut g = (f(&liquid, nl)? - f(&vapor, nv)?) / (nl - nv);
352294

353295
// update filled pore with limited step size
354296
let mut bulk = StateBuilder::new(functional)

0 commit comments

Comments
 (0)