Skip to content
This repository was archived by the owner on Jul 28, 2022. It is now read-only.

Added DFTSpecification interface to PoreProfile.initialize #4

Open
JohannesEller wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Added DFTSpecification interface to PoreProfile.initialize #4

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
9ebee6e
Added DFTSpecification interface to PoreProfile.initialize; started i…
Nov 13, 2021
7af3607
Added isostere calculation
Nov 13, 2021
cd9559e
Update isostere calculation; TODO: python interface
Nov 28, 2021
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
Prev Previous commit
Next Next commit
Added isostere calculation
  • Loading branch information
Johannes Eller committed Nov 13, 2021
commit 7af3607bbc5157eae856c0ac8131afc9450bfe51
182 changes: 84 additions & 98 deletions src/adsorption/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,12 +1,11 @@
//! Adsorption profiles and isotherms.
use super::functional::{HelmholtzEnergyFunctional, DFT};
// use super::profile::DFTSpecifications;
use super::profile::DFTSpecifications;
use super::solver::DFTSolver;
use feos_core::{
Contributions, EosError, EosResult, EosUnit, EquationOfState, StateBuilder, VLEOptions,
};
use ndarray::{arr1, Array1, Dimension, Ix1, Ix3};
// use quantity::si::{KB, NAV};
use quantity::{QuantityArray1, QuantityArray2, QuantityScalar};
use std::rc::Rc;

Expand All @@ -19,69 +18,41 @@ 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 temperature grid of adsorption isostere.
// pub enum TemperatureSpecification<U> {
// /// Specify the minimal and maximal pressure, and the number of points
// Tlim {
// t_min: QuantityScalar<U>,
// t_max: QuantityScalar<U>,
// points: usize,
// },
// /// Provide a custom array of pressure points.
// Tvec(QuantityArray1<U>),
// }

// impl<U: EosUnit> TemperatureSpecification<U>
// where
// QuantityScalar<U>: std::fmt::Display,
// {
// fn to_vec(&self) -> EosResult<QuantityArray1<U>> {
// Ok(match self {
// Self::Tlim {
// t_min,
// t_max,
// points,
// } => QuantityArray1::linspace(*t_min, *t_max, *points)?,
// Self::Tvec(temperature) => temperature.clone(),
// })
// }
// }

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

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

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

Expand All @@ -94,15 +65,15 @@ where
{
let p_eq = equilibrium.pressure().get(0);
match self {
Self::Plim {
p_min,
p_max,
Self::Qlim {
q_min,
q_max,
points,
} => Ok((
QuantityArray1::linspace(*p_min, p_eq, points / 2)?,
QuantityArray1::linspace(*p_max, p_eq, points - points / 2)?,
QuantityArray1::linspace(*q_min, p_eq, points / 2)?,
QuantityArray1::linspace(*q_max, p_eq, points - points / 2)?,
)),
Self::Pvec(pressure) => {
Self::Qvec(pressure) => {
let index = (0..pressure.len()).find(|&i| pressure.get(i) > p_eq);
Ok(if let Some(index) = index {
(
Expand Down Expand Up @@ -146,7 +117,7 @@ where
pub fn adsorption_isotherm<S: PoreSpecification<U, D, F>>(
functional: &Rc<DFT<F>>,
temperature: QuantityScalar<U>,
pressure: &PressureSpecification<U>,
pressure: &QuantitySpecification<U>,
pore: &S,
molefracs: Option<&Array1<f64>>,
solver: Option<&DFTSolver>,
Expand All @@ -159,7 +130,7 @@ where
pub fn desorption_isotherm<S: PoreSpecification<U, D, F>>(
functional: &Rc<DFT<F>>,
temperature: QuantityScalar<U>,
pressure: &PressureSpecification<U>,
pressure: &QuantitySpecification<U>,
pore: &S,
molefracs: Option<&Array1<f64>>,
solver: Option<&DFTSolver>,
Expand All @@ -178,12 +149,12 @@ where
pub fn equilibrium_isotherm<S: PoreSpecification<U, D, F>>(
functional: &Rc<DFT<F>>,
temperature: QuantityScalar<U>,
pressure: &PressureSpecification<U>,
pressure: &QuantitySpecification<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.q_min_max();
let equilibrium = Self::phase_equilibrium(
functional,
temperature,
Expand Down Expand Up @@ -405,46 +376,61 @@ where
))
}

// pub fn isostere<S: PoreSpecification<U, D, F>>(
// &self,
// functional: &Rc<DFT<F>>,
// initial_pressure: QuantityScalar<U>,
// temperature: &TemperatureSpecification<U>,
// pore: &S,
// molefracs: Option<&Array1<f64>>,
// solver: Option<&DFTSolver>,
// ) -> EosResult<Adsorption<U, D, F>> {
// let temperature = temperature.to_vec()?;
// 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(temperature.len());

// // Create the initial bulk state at the start of the isostere
// let mut initial_bulk = StateBuilder::new(functional)
// .temperature(temperature.get(0))
// .pressure(initial_pressure)
// .moles(&moles)
// .build()?;
// if functional.components() > 1 && !initial_bulk.is_stable(VLEOptions::default())? {
// initial_bulk = initial_bulk
// .tp_flash(None, VLEOptions::default(), None)?
// .vapor()
// .clone();
// }

// // Initial PoreProfile for given mu,V,T
// let initial_profile = pore.initialize(&initial_bulk, None, None)?.solve(solver);

// // Specify the number of particles for the isostere calculation
// let n = initial_profile?.profile.moles().to_reduced(NAV)?;
// let specification = Rc::new(DFTSpecifications::Moles { moles: n });

// for i in 1..temperature.len() {}
// let dummy_state =

// unimplemented!()
// }
pub fn isostere<S: PoreSpecification<U, D, F>>(
&self,
functional: &Rc<DFT<F>>,
initial_pressure: QuantityScalar<U>,
temperature: &QuantitySpecification<U>,
pore: &S,
molefracs: Option<&Array1<f64>>,
solver: Option<&DFTSolver>,
) -> EosResult<Adsorption<U, D, F>> {
let temperature = temperature.to_vec()?;
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(temperature.len());

// Create the initial bulk state at the start of the isostere
let mut initial_bulk = StateBuilder::new(functional)
.temperature(temperature.get(0))
.pressure(initial_pressure)
.moles(&moles)
.build()?;
if functional.components() > 1 && !initial_bulk.is_stable(VLEOptions::default())? {
initial_bulk = initial_bulk
.tp_flash(None, VLEOptions::default(), None)?
.vapor()
.clone();
}

// Initial PoreProfile for given mu,V,T
let initial_profile = pore.initialize(&initial_bulk, None, None)?.solve(solver)?;

let specification = DFTSpecifications::moles_from_profile(&initial_profile.profile)?;

profiles.push(Ok(initial_profile.clone()));

let partial_density = initial_profile.profile.moles() / initial_profile.profile.volume();

for i in 1..temperature.len() {
let dummy_state = StateBuilder::new(functional)
.temperature(temperature.get(i))
.partial_density(&partial_density)
.build()?;

let mut p = pore
.initialize(&dummy_state, None, Some(specification.clone()))?
.solve(solver)?;
p.profile
.bulk
.update_chemical_potential(&p.profile.chemical_potential)?;

profiles.push(Ok(p));
}

Ok(Adsorption(profiles, functional.components()))
}

pub fn pressure(&self) -> QuantityArray1<U> {
QuantityArray1::from_shape_fn(self.0.len(), |i| match &self.0[i] {
Expand Down
6 changes: 3 additions & 3 deletions src/profile.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,14 +63,14 @@ impl DFTSpecifications {
/// particles constant in systems, where the number itself is difficult to obtain.
pub fn moles_from_profile<U: EosUnit, D: Dimension, F: HelmholtzEnergyFunctional>(
profile: &DFTProfile<U, D, F>,
) -> EosResult<Rc<Self>>
) -> EosResult<Self>
where
<D as Dimension>::Larger: Dimension<Smaller = D>,
{
let rho = profile.density.to_reduced(U::reference_density())?;
Ok(Rc::new(Self::Moles {
Ok(Self::Moles {
moles: profile.integrate_reduced_comp(&rho),
}))
})
}

/// Calculate the number of particles from the profile.
Expand Down