Skip to content

Improve the robustness of adsorption isotherm calculations #50

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
prehner merged 2 commits into from
Aug 26, 2022
Merged

Improve the robustness of adsorption isotherm calculations #50

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
855351b
Improve the robustness of adsorption isotherm calculations
prehner Aug 26, 2022
869ab1c
Update notebooks
prehner Aug 26, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Next Next commit
Improve the robustness of adsorption isotherm calculations
  • Loading branch information
@prehner
prehner committed Aug 26, 2022
commit 855351bda54cf8e59c0c7708f99327138992cb65
7 changes: 6 additions & 1 deletion feos-dft/CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0

## [Unreleased]

## [0.3.1] - 2022-08-26
### Changed
- `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)
- 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)

## [0.3.0] - 2022-07-13
### Added
- Added getters for the fields of `Pore1D` in Python. [#30](https://github.com/feos-org/feos-dft/pull/30)
Expand All @@ -16,7 +21,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
- Changed interface of `PairCorrelationFunction` to facilitate the calculation of pair correlation functions in mixtures. [#29](https://github.com/feos-org/feos-dft/pull/29)

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

## [0.2.0] - 2022-04-12
### Added
Expand Down
194 changes: 68 additions & 126 deletions feos-dft/src/adsorption/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,10 +2,12 @@
use super::functional::{HelmholtzEnergyFunctional, DFT};
use super::solver::DFTSolver;
use feos_core::{
Contributions, EosError, EosResult, EosUnit, EquationOfState, SolverOptions, StateBuilder,
Contributions, DensityInitialization, EosError, EosResult, EosUnit, EquationOfState,
SolverOptions, State, StateBuilder,
};
use ndarray::{arr1, Array1, Dimension, Ix1, Ix3, RemoveAxis};
use ndarray::{Array1, Dimension, Ix1, Ix3, RemoveAxis};
use quantity::{QuantityArray1, QuantityArray2, QuantityScalar};
use std::iter;
use std::rc::Rc;

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

/// Possible inputs for the pressure grid of adsorption isotherms.
pub enum PressureSpecification<U> {
/// Specify the minimal and maximal pressure, and the number of points
Plim {
p_min: QuantityScalar<U>,
p_max: QuantityScalar<U>,
points: usize,
},
/// Provide a custom array of pressure points.
Pvec(QuantityArray1<U>),
}

impl<U: EosUnit> PressureSpecification<U>
where
QuantityScalar<U>: std::fmt::Display,
{
fn p_min_max(&self) -> (QuantityScalar<U>, QuantityScalar<U>) {
match self {
Self::Plim {
p_min,
p_max,
points: _,
} => (*p_min, *p_max),
Self::Pvec(pressure) => (pressure.get(0), pressure.get(pressure.len() - 1)),
}
}

fn to_vec(&self) -> EosResult<QuantityArray1<U>> {
Ok(match self {
Self::Plim {
p_min,
p_max,
points,
} => QuantityArray1::linspace(*p_min, *p_max, *points)?,
Self::Pvec(pressure) => pressure.clone(),
})
}

fn equilibrium<
D: Dimension + RemoveAxis + 'static,
F: HelmholtzEnergyFunctional + FluidParameters,
>(
&self,
equilibrium: &Adsorption<U, D, F>,
) -> EosResult<(QuantityArray1<U>, QuantityArray1<U>)>
where
D::Larger: Dimension<Smaller = D>,
D::Smaller: Dimension<Larger = D>,
<D::Larger as Dimension>::Larger: Dimension<Smaller = D::Larger>,
{
let p_eq = equilibrium.pressure().get(0);
match self {
Self::Plim {
p_min,
p_max,
points,
} => Ok((
QuantityArray1::linspace(*p_min, p_eq, points / 2)?,
QuantityArray1::linspace(*p_max, p_eq, points - points / 2)?,
)),
Self::Pvec(pressure) => {
let index = (0..pressure.len()).find(|&i| pressure.get(i) > p_eq);
Ok(if let Some(index) = index {
(
QuantityArray1::from_shape_fn(index + 1, |i| {
if i == index {
p_eq
} else {
pressure.get(i)
}
}),
QuantityArray1::from_shape_fn(pressure.len() - index + 1, |i| {
if i == pressure.len() - index {
p_eq
} else {
pressure.get(pressure.len() - i - 1)
}
}),
)
} else {
(pressure.clone(), arr1(&[]) * U::reference_pressure())
})
}
}
}
}

/// Container structure for the calculation of adsorption isotherms.
pub struct Adsorption<U, D: Dimension, F> {
components: usize,
Expand Down Expand Up @@ -143,28 +58,41 @@ where
pub fn adsorption_isotherm<S: PoreSpecification<U, D>>(
functional: &Rc<DFT<F>>,
temperature: QuantityScalar<U>,
pressure: &PressureSpecification<U>,
pressure: &QuantityArray1<U>,
pore: &S,
molefracs: Option<&Array1<f64>>,
solver: Option<&DFTSolver>,
) -> EosResult<Adsorption<U, D, F>> {
let pressure = pressure.to_vec()?;
Self::isotherm(functional, temperature, &pressure, pore, molefracs, solver)
Self::isotherm(
functional,
temperature,
pressure,
pore,
molefracs,
DensityInitialization::Vapor,
solver,
)
}

/// Calculate an desorption isotherm (starting at high pressure)
pub fn desorption_isotherm<S: PoreSpecification<U, D>>(
functional: &Rc<DFT<F>>,
temperature: QuantityScalar<U>,
pressure: &PressureSpecification<U>,
pressure: &QuantityArray1<U>,
pore: &S,
molefracs: Option<&Array1<f64>>,
solver: Option<&DFTSolver>,
) -> EosResult<Adsorption<U, D, F>> {
let pressure = pressure.to_vec()?;
let pressure =
QuantityArray1::from_shape_fn(pressure.len(), |i| pressure.get(pressure.len() - i - 1));
let isotherm = Self::isotherm(functional, temperature, &pressure, pore, molefracs, solver)?;
let pressure = pressure.into_iter().rev().collect();
let isotherm = Self::isotherm(
functional,
temperature,
&pressure,
pore,
molefracs,
DensityInitialization::Liquid,
solver,
)?;
Ok(Adsorption::new(
functional,
pore,
Expand All @@ -176,12 +104,12 @@ where
pub fn equilibrium_isotherm<S: PoreSpecification<U, D>>(
functional: &Rc<DFT<F>>,
temperature: QuantityScalar<U>,
pressure: &PressureSpecification<U>,
pressure: &QuantityArray1<U>,
pore: &S,
molefracs: Option<&Array1<f64>>,
solver: Option<&DFTSolver>,
) -> EosResult<Adsorption<U, D, F>> {
let (p_min, p_max) = pressure.p_min_max();
let (p_min, p_max) = (pressure.get(0), pressure.get(pressure.len() - 1));
let equilibrium = Self::phase_equilibrium(
functional,
temperature,
Expand All @@ -193,20 +121,28 @@ where
SolverOptions::default(),
);
if let Ok(equilibrium) = equilibrium {
let pressure = pressure.equilibrium(&equilibrium)?;
let adsorption = Self::isotherm(
let p_eq = equilibrium.pressure().get(0);
let p_ads = pressure
.into_iter()
.filter(|&p| p <= p_eq)
.chain(iter::once(p_eq))
.collect();
let p_des = iter::once(p_eq)
.chain(pressure.into_iter().filter(|&p| p > p_eq))
.collect();
let adsorption = Self::adsorption_isotherm(
functional,
temperature,
&pressure.0,
&p_ads,
pore,
molefracs,
solver,
)?
.profiles;
let desorption = Self::isotherm(
let desorption = Self::desorption_isotherm(
functional,
temperature,
&pressure.1,
&p_des,
pore,
molefracs,
solver,
Expand All @@ -215,7 +151,7 @@ where
Ok(Adsorption {
profiles: adsorption
.into_iter()
.chain(desorption.into_iter().rev())
.chain(desorption.into_iter())
.collect(),
components: functional.components(),
dimension: pore.dimension(),
Expand Down Expand Up @@ -258,28 +194,33 @@ where
pressure: &QuantityArray1<U>,
pore: &S,
molefracs: Option<&Array1<f64>>,
density_initialization: DensityInitialization<U>,
solver: Option<&DFTSolver>,
) -> EosResult<Adsorption<U, D, F>> {
let moles =
functional.validate_moles(molefracs.map(|x| x * U::reference_moles()).as_ref())?;
let mut profiles: Vec<EosResult<PoreProfile<U, D, F>>> = Vec::with_capacity(pressure.len());

// Calculate the external potential once
let mut bulk = StateBuilder::new(functional)
.temperature(temperature)
.pressure(pressure.get(0))
.moles(&moles)
.build()?;
// On the first iteration, initialize the density profile according to the direction
// and calculate the external potential once.
let mut bulk = State::new_npt(
functional,
temperature,
pressure.get(0),
&moles,
density_initialization,
)?;
if functional.components() > 1 && !bulk.is_stable(SolverOptions::default())? {
bulk = bulk
.tp_flash(None, SolverOptions::default(), None)?
.vapor()
.clone();
let vle = bulk.tp_flash(None, SolverOptions::default(), None)?;
bulk = match density_initialization {
DensityInitialization::Liquid => vle.liquid().clone(),
DensityInitialization::Vapor => vle.vapor().clone(),
_ => unreachable!(),
};
}
let external_potential = pore
.initialize(&bulk, None, None)?
.profile
.external_potential;
let profile = pore.initialize(&bulk, None, None)?.solve(solver)?.profile;
let external_potential = Some(&profile.external_potential);
let mut old_density = Some(&profile.density);

for i in 0..pressure.len() {
let mut bulk = StateBuilder::new(functional)
Expand All @@ -293,15 +234,16 @@ where
.vapor()
.clone();
}
let old_density = if let Some(Ok(l)) = profiles.last() {

let p = pore.initialize(&bulk, old_density, external_potential)?;
let p2 = pore.initialize(&bulk, None, external_potential)?;
profiles.push(p.solve(solver).or_else(|_| p2.solve(solver)));

old_density = if let Some(Ok(l)) = profiles.last() {
Some(&l.profile.density)
} else {
None
};

let p = pore.initialize(&bulk, old_density, Some(&external_potential))?;
let p2 = pore.initialize(&bulk, None, Some(&external_potential))?;
profiles.push(p.solve(solver).or_else(|_| p2.solve(solver)));
}

Ok(Adsorption::new(functional, pore, profiles))
Expand Down Expand Up @@ -344,11 +286,11 @@ where
let nl = liquid.profile.bulk.density
* (liquid.profile.moles() * liquid.profile.bulk.molar_volume(Contributions::Total))
.sum();
let f = |s: &mut PoreProfile<U, D, F>, n: QuantityScalar<U>| -> EosResult<_> {
let f = |s: &PoreProfile<U, D, F>, n: QuantityScalar<U>| -> EosResult<_> {
Ok(s.grand_potential.unwrap()
+ s.profile.bulk.molar_gibbs_energy(Contributions::Total) * n)
};
let mut g = (f(&mut liquid, nl)? - f(&mut vapor, nv)?) / (nl - nv);
let mut g = (f(&liquid, nl)? - f(&vapor, nv)?) / (nl - nv);

// update filled pore with limited step size
let mut bulk = StateBuilder::new(functional)
Expand Down
Loading