Skip to content

Commit

Permalink
QoL improvements for working with ideal gas models
Browse files Browse the repository at this point in the history
  • Loading branch information
prehner committed Nov 28, 2023
1 parent 01d8238 commit d0a8c82
Show file tree
Hide file tree
Showing 17 changed files with 149 additions and 102 deletions.
13 changes: 9 additions & 4 deletions feos-core/src/cubic.rs
Original file line number Diff line number Diff line change
Expand Up @@ -26,15 +26,18 @@ pub struct PengRobinsonRecord {
pc: f64,
/// acentric factor
acentric_factor: f64,
/// molar weight
molarweight: f64,
}

impl PengRobinsonRecord {
/// Create a new pure substance record for the Peng-Robinson equation of state.
pub fn new(tc: f64, pc: f64, acentric_factor: f64) -> Self {
pub fn new(tc: f64, pc: f64, acentric_factor: f64, molarweight: f64) -> Self {
Self {
tc,
pc,
acentric_factor,
molarweight,
}
}
}
Expand All @@ -43,7 +46,8 @@ impl std::fmt::Display for PengRobinsonRecord {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
write!(f, "PengRobinsonRecord(tc={} K", self.tc)?;
write!(f, ", pc={} Pa", self.pc)?;
write!(f, ", acentric factor={}", self.acentric_factor)
write!(f, ", acentric factor={}", self.acentric_factor)?;
write!(f, ", molar weight={})", self.molarweight)
}
}

Expand Down Expand Up @@ -93,9 +97,10 @@ impl PengRobinsonParameters {
tc: tc[i],
pc: pc[i],
acentric_factor: acentric_factor[i],
molarweight: molarweight[i],
};
let id = Identifier::default();
PureRecord::new(id, molarweight[i], record)
PureRecord::new(id, record)
})
.collect();
PengRobinsonParameters::from_records(records, None)
Expand All @@ -120,8 +125,8 @@ impl Parameter for PengRobinsonParameters {
let mut kappa = Array1::zeros(n);

for (i, record) in pure_records.iter().enumerate() {
molarweight[i] = record.molarweight;
let r = &record.model_record;
molarweight[i] = r.molarweight;
tc[i] = r.tc;
a[i] = 0.45724 * r.tc.powi(2) * KB_A3 / r.pc;
b[i] = 0.07780 * r.tc * KB_A3 / r.pc;
Expand Down
17 changes: 14 additions & 3 deletions feos-core/src/equation_of_state/mod.rs
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ mod residual;

pub use helmholtz_energy::{HelmholtzEnergy, HelmholtzEnergyDual};
pub use ideal_gas::{DeBroglieWavelength, DeBroglieWavelengthDual, IdealGas};
use residual::NoResidual;
pub use residual::{EntropyScaling, Residual};

/// The number of components that the model is initialized for.
Expand All @@ -27,14 +28,24 @@ pub trait Components {
/// and a residual Helmholtz energy model.
#[derive(Clone)]
pub struct EquationOfState<I, R> {
pub ideal_gas: Arc<I>,
pub ideal_gas: I,
pub residual: Arc<R>,
}

impl<I, R> EquationOfState<I, R> {
/// Return a new [EquationOfState] with the given ideal gas
/// and residual models.
pub fn new(ideal_gas: Arc<I>, residual: Arc<R>) -> Self {
pub fn new(ideal_gas: I, residual: Arc<R>) -> Self {
Self {
ideal_gas,
residual,
}
}
}

impl<I: IdealGas> EquationOfState<I, NoResidual> {
pub fn ideal_gas(ideal_gas: I) -> Self {
let residual = Arc::new(NoResidual(ideal_gas.components()));
Self {
ideal_gas,
residual,
Expand All @@ -54,7 +65,7 @@ impl<I: Components, R: Components> Components for EquationOfState<I, R> {

fn subset(&self, component_list: &[usize]) -> Self {
Self::new(
Arc::new(self.ideal_gas.subset(component_list)),
self.ideal_gas.subset(component_list),
Arc::new(self.residual.subset(component_list)),
)
}
Expand Down
26 changes: 26 additions & 0 deletions feos-core/src/equation_of_state/residual.rs
Original file line number Diff line number Diff line change
Expand Up @@ -178,3 +178,29 @@ pub trait EntropyScaling {
) -> EosResult<ThermalConductivity>;
fn thermal_conductivity_correlation(&self, s_res: f64, x: &Array1<f64>) -> EosResult<f64>;
}

pub struct NoResidual(pub usize);

impl Components for NoResidual {
fn components(&self) -> usize {
self.0
}

fn subset(&self, component_list: &[usize]) -> Self {
Self(component_list.len())
}
}

impl Residual for NoResidual {
fn compute_max_density(&self, _: &Array1<f64>) -> f64 {
1.0
}

fn contributions(&self) -> &[Box<dyn HelmholtzEnergy>] {
&[]
}

fn molar_weight(&self) -> MolarWeight<Array1<f64>> {
unimplemented!()
}
}
34 changes: 9 additions & 25 deletions feos-core/src/joback.rs
Original file line number Diff line number Diff line change
Expand Up @@ -240,32 +240,15 @@ const KB: f64 = 1.38064852e-23;

#[cfg(test)]
mod tests {
use crate::si::*;
use crate::{Contributions, Residual, State, StateBuilder};
use crate::{si::*, EquationOfState};
use crate::{Contributions, State, StateBuilder};
use approx::assert_relative_eq;
use ndarray::arr1;
use std::sync::Arc;
use typenum::P3;

use super::*;

// implement Residual to test Joback as equation of state
impl Residual for Joback {
fn compute_max_density(&self, _moles: &Array1<f64>) -> f64 {
1.0
}

fn contributions(&self) -> &[Box<dyn crate::HelmholtzEnergy>] {
&[]
}

fn molar_weight(&self) -> MolarWeight<Array1<f64>> {
MolarWeight::from_shape_fn(self.components(), |i| {
self.parameters.pure_records[i].molarweight * GRAM / MOL
})
}
}

#[test]
fn paper_example() -> EosResult<()> {
let segments_json = r#"[
Expand Down Expand Up @@ -350,8 +333,9 @@ mod tests {
);
assert_relative_eq!(jr.e, 0.0);

let pr = PureRecord::new(Identifier::default(), 1.0, jr);
let eos = Arc::new(Joback::new(Arc::new(JobackParameters::new_pure(pr)?)));
let pr = PureRecord::new(Identifier::default(), jr);
let joback = Joback::new(Arc::new(JobackParameters::new_pure(pr)?));
let eos = Arc::new(EquationOfState::ideal_gas(joback));
let state = State::new_nvt(
&eos,
1000.0 * KELVIN,
Expand All @@ -373,20 +357,20 @@ mod tests {
fn c_p_comparison() -> EosResult<()> {
let record1 = PureRecord::new(
Identifier::default(),
1.0,
JobackRecord::new(1.0, 0.2, 0.03, 0.004, 0.005),
);
let record2 = PureRecord::new(
Identifier::default(),
1.0,
JobackRecord::new(-5.0, 0.4, 0.03, 0.002, 0.001),
);
let parameters = Arc::new(JobackParameters::new_binary(vec![record1, record2], None)?);
let joback = Arc::new(Joback::new(parameters));
let joback = Joback::new(parameters.clone());
let ideal_gas = Joback::new(parameters);
let eos = Arc::new(EquationOfState::ideal_gas(ideal_gas));
let temperature = 300.0 * KELVIN;
let volume = METER.powi::<P3>();
let moles = &arr1(&[1.0, 3.0]) * MOL;
let state = StateBuilder::new(&joback)
let state = StateBuilder::new(&eos)
.temperature(temperature)
.volume(volume)
.moles(&moles)
Expand Down
7 changes: 3 additions & 4 deletions feos-core/src/lib.rs
Original file line number Diff line number Diff line change
Expand Up @@ -173,14 +173,13 @@ mod tests {
fn validate_residual_properties() -> EosResult<()> {
let mixture = pure_record_vec();
let propane = mixture[0].clone();
let parameters = PengRobinsonParameters::new_pure(propane)?;
let residual = Arc::new(PengRobinson::new(Arc::new(parameters)));
let parameters = Arc::new(PengRobinsonParameters::new_pure(propane)?);
let residual = Arc::new(PengRobinson::new(parameters));
let joback_parameters = Arc::new(JobackParameters::new_pure(PureRecord::new(
Identifier::default(),
1.0,
JobackRecord::new(0.0, 0.0, 0.0, 0.0, 0.0),
))?);
let ideal_gas = Arc::new(Joback::new(joback_parameters));
let ideal_gas = Joback::new(joback_parameters);
let eos = Arc::new(EquationOfState::new(ideal_gas, residual.clone()));

let sr = StateBuilder::new(&residual)
Expand Down
2 changes: 1 addition & 1 deletion feos-core/src/parameter/mod.rs
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ where
fn from_model_records(model_records: Vec<Self::Pure>) -> Result<Self, ParameterError> {
let pure_records = model_records
.into_iter()
.map(|r| PureRecord::new(Default::default(), Default::default(), r))
.map(|r| PureRecord::new(Default::default(), r))
.collect();
Self::from_records(pure_records, None)
}
Expand Down
19 changes: 6 additions & 13 deletions feos-core/src/parameter/model_record.rs
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
use super::identifier::Identifier;
use super::segment::SegmentRecord;
use super::{IdentifierOption, ParameterError};
use conv::ValueInto;
use serde::de::DeserializeOwned;
use serde::{Deserialize, Serialize};
use std::collections::{HashMap, HashSet};
Expand All @@ -13,16 +12,14 @@ use std::path::Path;
#[derive(Serialize, Deserialize, Debug, Clone)]
pub struct PureRecord<M> {
pub identifier: Identifier,
pub molarweight: f64,
pub model_record: M,
}

impl<M> PureRecord<M> {
/// Create a new `PureRecord`.
pub fn new(identifier: Identifier, molarweight: f64, model_record: M) -> Self {
pub fn new(identifier: Identifier, model_record: M) -> Self {
Self {
identifier,
molarweight,
model_record,
}
}
Expand All @@ -33,19 +30,16 @@ impl<M> PureRecord<M> {
/// and the ideal gas record.
pub fn from_segments<S, T>(identifier: Identifier, segments: S) -> Result<Self, ParameterError>
where
T: Copy + ValueInto<f64>,
M: FromSegments<T>,
S: IntoIterator<Item = (SegmentRecord<M>, T)>,
{
let mut molarweight = 0.0;
let mut model_segments = Vec::new();
for (s, n) in segments {
molarweight += s.molarweight * n.value_into().unwrap();
model_segments.push((s.model_record, n));
}
let model_segments: Vec<_> = segments
.into_iter()
.map(|(s, n)| (s.model_record, n))
.collect();
let model_record = M::from_segments(&model_segments)?;

Ok(Self::new(identifier, molarweight, model_record))
Ok(Self::new(identifier, model_record))
}

/// Create pure substance parameters from a json file.
Expand Down Expand Up @@ -104,7 +98,6 @@ where
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
write!(f, "PureRecord(")?;
write!(f, "\n\tidentifier={},", self.identifier)?;
write!(f, "\n\tmolarweight={},", self.molarweight)?;
write!(f, "\n\tmodel_record={},", self.model_record)?;
write!(f, "\n)")
}
Expand Down
9 changes: 7 additions & 2 deletions feos-core/src/python/cubic.rs
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,13 @@ pub struct PyPengRobinsonRecord(PengRobinsonRecord);
#[pymethods]
impl PyPengRobinsonRecord {
#[new]
fn new(tc: f64, pc: f64, acentric_factor: f64) -> Self {
Self(PengRobinsonRecord::new(tc, pc, acentric_factor))
fn new(tc: f64, pc: f64, acentric_factor: f64, molarweight: f64) -> Self {
Self(PengRobinsonRecord::new(
tc,
pc,
acentric_factor,
molarweight,
))
}

fn __repr__(&self) -> PyResult<String> {
Expand Down
22 changes: 2 additions & 20 deletions feos-core/src/python/parameter/mod.rs
Original file line number Diff line number Diff line change
Expand Up @@ -399,16 +399,8 @@ macro_rules! impl_pure_record {
impl PyPureRecord {
#[new]
#[pyo3(text_signature = "(identifier, molarweight, model_record)")]
fn new(
identifier: PyIdentifier,
molarweight: f64,
model_record: $py_model_record,
) -> PyResult<Self> {
Ok(Self(PureRecord::new(
identifier.0,
molarweight,
model_record.0,
)))
fn new(identifier: PyIdentifier, model_record: $py_model_record) -> PyResult<Self> {
Ok(Self(PureRecord::new(identifier.0, model_record.0)))
}

#[getter]
Expand All @@ -421,16 +413,6 @@ macro_rules! impl_pure_record {
self.0.identifier = identifier.0;
}

#[getter]
fn get_molarweight(&self) -> f64 {
self.0.molarweight
}

#[setter]
fn set_molarweight(&mut self, molarweight: f64) {
self.0.molarweight = molarweight;
}

#[getter]
fn get_model_record(&self) -> $py_model_record {
$py_model_record(self.0.model_record.clone())
Expand Down
30 changes: 24 additions & 6 deletions feos-derive/src/residual.rs
Original file line number Diff line number Diff line change
Expand Up @@ -24,20 +24,38 @@ fn impl_residual(
) -> proc_macro2::TokenStream {
let compute_max_density = variants.iter().map(|v| {
let name = &v.ident;
quote! {
Self::#name(residual) => residual.compute_max_density(moles)
if name == "NoModel" {
quote! {
Self::#name(_) => 0.0
}
} else {
quote! {
Self::#name(residual) => residual.compute_max_density(moles)
}
}
});
let contributions = variants.iter().map(|v| {
let name = &v.ident;
quote! {
Self::#name(residual) => residual.contributions()
if name == "NoModel" {
quote! {
Self::#name(_) => &[]
}
} else {
quote! {
Self::#name(residual) => residual.contributions()
}
}
});
let molar_weight = variants.iter().map(|v| {
let name = &v.ident;
quote! {
Self::#name(residual) => residual.molar_weight()
if name == "NoModel" {
quote! {
Self::#name(_) => panic!("OH NO")
}
} else {
quote! {
Self::#name(residual) => residual.molar_weight()
}
}
});

Expand Down
2 changes: 1 addition & 1 deletion feos-dft/src/functional.rs
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,7 @@ impl<F> Deref for DFT<F> {

impl<F> DFT<F> {
pub fn ideal_gas<I>(self, ideal_gas: I) -> DFT<EquationOfState<I, F>> {
DFT(EquationOfState::new(Arc::new(ideal_gas), Arc::new(self.0)))
DFT(EquationOfState::new(ideal_gas, Arc::new(self.0)))
}
}

Expand Down
1 change: 1 addition & 0 deletions src/eos.rs
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@ use ndarray::Array1;
/// are undesirable (e.g. FFI).
#[derive(Components, Residual)]
pub enum ResidualModel {
NoModel(usize),
#[cfg(feature = "pcsaft")]
#[implement(entropy_scaling)]
PcSaft(PcSaft),
Expand Down
Loading

0 comments on commit d0a8c82

Please sign in to comment.