feos-org · JohannesEller · Nov 13, 2021 · Nov 13, 2021 · Nov 28, 2021
diff --git a/src/adsorption/mod.rs b/src/adsorption/mod.rs
@@ -1,5 +1,8 @@
 //! Adsorption profiles and isotherms.
+use crate::DFTSpecification;
+
 use super::functional::{HelmholtzEnergyFunctional, DFT};
+use super::profile::DFTSpecifications;
 use super::solver::DFTSolver;
 use feos_core::{
     Contributions, EosError, EosResult, EosUnit, EquationOfState, StateBuilder, VLEOptions,
@@ -17,41 +20,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 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(),
         })
     }
 
@@ -64,15 +67,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 {
                     (
@@ -116,7 +119,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>,
@@ -129,7 +132,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>,
@@ -148,12 +151,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,
@@ -247,7 +250,10 @@ where
                 .vapor()
                 .clone();
         }
-        let external_potential = pore.initialize(&bulk, None)?.profile.external_potential;
+        let external_potential = pore
+            .initialize(&bulk, None, None)?
+            .profile
+            .external_potential;
 
         for i in 0..pressure.len() {
             let mut bulk = StateBuilder::new(functional)
@@ -261,7 +267,7 @@ where
                     .vapor()
                     .clone();
             }
-            let mut p = pore.initialize(&bulk, Some(&external_potential))?;
+            let mut p = pore.initialize(&bulk, Some(&external_potential), None)?;
             let p2 = p.clone();
             if let Some(Ok(l)) = profiles.last() {
                 p.profile.density = l.profile.density.clone();
@@ -300,8 +306,8 @@ where
             .liquid()
             .build()?;
 
-        let mut vapor = pore.initialize(&vapor_bulk, None)?.solve(None)?;
-        let mut liquid = pore.initialize(&liquid_bulk, None)?.solve(solver)?;
+        let mut vapor = pore.initialize(&vapor_bulk, None, None)?.solve(None)?;
+        let mut liquid = pore.initialize(&liquid_bulk, None, None)?.solve(solver)?;
 
         // calculate initial value for the molar gibbs energy
         let nv = vapor.profile.bulk.density
@@ -372,6 +378,62 @@ where
         ))
     }
 
+    /// Calculate the isostere starting from the lowest temperature.
+    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 = if let Ok(initial_profile) = &initial_profile {
+            DFTSpecifications::moles_from_profile(&initial_profile.profile)?
+        } else {
+            unimplemented!()
+        };
+
+        profiles.push(initial_profile);
+
+        for i in 1..temperature.len() {
+            let bulk = StateBuilder::new(functional)
+                .temperature(temperature.get(i))
+                .partial_density(&initial_bulk.partial_density)
+                .build()?;
+            let mut p = pore.initialize(&bulk, None, Some(specification.clone()))?;
+            if let Some(Ok(l)) = profiles.last() {
+                p.profile.density = l.profile.density.clone();
+            }
+            profiles.push(p.solve(solver));
+        }
+
+        Ok(Adsorption(profiles, functional.components()))
+    }
+
     pub fn pressure(&self) -> QuantityArray1<U> {
         QuantityArray1::from_shape_fn(self.0.len(), |i| match &self.0[i] {
             Ok(p) => {

diff --git a/src/adsorption/pore.rs b/src/adsorption/pore.rs
@@ -2,7 +2,7 @@ use crate::adsorption::{ExternalPotential, FluidParameters};
 use crate::convolver::ConvolverFFT;
 use crate::functional::{HelmholtzEnergyFunctional, DFT};
 use crate::geometry::{Axis, AxisGeometry, Grid};
-use crate::profile::{DFTProfile, CUTOFF_RADIUS, MAX_POTENTIAL};
+use crate::profile::{DFTProfile, DFTSpecifications, CUTOFF_RADIUS, MAX_POTENTIAL};
 use crate::solver::DFTSolver;
 use feos_core::{Contributions, EosResult, EosUnit, State};
 use ndarray::prelude::*;
@@ -88,6 +88,7 @@ pub trait PoreSpecification<U, D: Dimension, F> {
         &self,
         bulk: &State<U, DFT<F>>,
         external_potential: Option<&Array<f64, D::Larger>>,
+        specification: Option<DFTSpecifications>,
     ) -> EosResult<PoreProfile<U, D, F>>;
 }
 
@@ -150,6 +151,14 @@ where
         self.interfacial_tension = None;
         self
     }
+
+    pub fn update_specification(&self, specification: DFTSpecifications) -> Self {
+        let mut profile = self.clone();
+        profile.profile.specification = Rc::new(specification);
+        profile.grand_potential = None;
+        profile.interfacial_tension = None;
+        profile
+    }
 }
 
 impl<U: EosUnit, F: HelmholtzEnergyFunctional + FluidParameters> PoreSpecification<U, Ix1, F>
@@ -159,6 +168,7 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional + FluidParameters> PoreSpecificati
         &self,
         bulk: &State<U, DFT<F>>,
         external_potential: Option<&Array2<f64>>,
+        specification: Option<DFTSpecifications>,
     ) -> EosResult<PoreProfile1D<U, F>> {
         let dft = &bulk.eos;
         let n_grid = self.n_grid.unwrap_or(DEFAULT_GRID_POINTS);
@@ -195,7 +205,13 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional + FluidParameters> PoreSpecificati
         let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1));
 
         Ok(PoreProfile {
-            profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential))?,
+            profile: DFTProfile::new(
+                grid,
+                convolver,
+                bulk,
+                Some(external_potential),
+                specification,
+            )?,
             grand_potential: None,
             interfacial_tension: None,
         })
@@ -209,6 +225,7 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional, P: FluidParameters> PoreSpecifica
         &self,
         bulk: &State<U, DFT<F>>,
         external_potential: Option<&Array4<f64>>,
+        specification: Option<DFTSpecifications>,
     ) -> EosResult<PoreProfile3D<U, F>> {
         let dft = &bulk.eos;
 
@@ -251,7 +268,13 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional, P: FluidParameters> PoreSpecifica
         let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1));
 
         Ok(PoreProfile {
-            profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential))?,
+            profile: DFTProfile::new(
+                grid,
+                convolver,
+                bulk,
+                Some(external_potential),
+                specification,
+            )?,
             grand_potential: None,
             interfacial_tension: None,
         })

diff --git a/src/interface/mod.rs b/src/interface/mod.rs
@@ -82,7 +82,7 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional> PlanarInterface<U, F> {
         let convolver = ConvolverFFT::plan(&grid, &weight_functions, None);
 
         Ok(Self {
-            profile: DFTProfile::new(grid, convolver, vle.vapor(), None)?,
+            profile: DFTProfile::new(grid, convolver, vle.vapor(), None, None)?,
             vle: vle.clone(),
             surface_tension: None,
             equimolar_radius: None,

diff --git a/src/profile.rs b/src/profile.rs
@@ -33,6 +33,7 @@ pub trait DFTSpecification<U, D: Dimension, F> {
 }
 
 /// Common specifications for the grand potentials in a DFT calculation.
+#[derive(Clone)]
 pub enum DFTSpecifications {
     /// DFT with specified chemical potential.
     ChemicalPotential,
@@ -49,21 +50,27 @@ pub enum DFTSpecifications {
     },
 }
 
+impl Default for DFTSpecifications {
+    fn default() -> Self {
+        Self::ChemicalPotential
+    }
+}
+
 impl DFTSpecifications {
     /// Calculate the number of particles from the profile.
     ///
     /// Call this after initializing the density profile to keep the number of
     /// 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.
@@ -198,6 +205,7 @@ where
         convolver: Rc<dyn Convolver<f64, D>>,
         bulk: &State<U, DFT<F>>,
         external_potential: Option<Array<f64, D::Larger>>,
+        specification: Option<DFTSpecifications>,
     ) -> EosResult<Self> {
         let dft = bulk.eos.clone();
 
@@ -224,14 +232,17 @@ where
                 .assign(&(isaft.index_axis(Axis_nd(0), s).map(|is| is.min(1.0)) * bulk_density[c]));
         }
 
+        // initialize DFT specification
+        let specification = specification.unwrap_or_default();
+
         Ok(Self {
             grid,
             convolver,
             dft: bulk.eos.clone(),
             temperature: bulk.temperature,
             density: density * U::reference_density(),
             chemical_potential: bulk.chemical_potential(Contributions::Total),
-            specification: Rc::new(DFTSpecifications::ChemicalPotential),
+            specification: Rc::new(specification),
             external_potential,
             bulk: bulk.clone(),
         })