Skip to content
This repository was archived by the owner on Jun 14, 2022. It is now read-only.
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
3 changes: 3 additions & 0 deletions src/errors.rs
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
use crate::parameter::ParameterError;
use argmin::core::Error as ArgminError;
use num_dual::linalg::LinAlgError;
use quantity::QuantityError;
Expand Down Expand Up @@ -25,6 +26,8 @@ pub enum EosError {
#[error(transparent)]
QuantityError(#[from] QuantityError),
#[error(transparent)]
ParameterError(#[from] ParameterError),
#[error(transparent)]
ArgminError(#[from] ArgminError),
#[error(transparent)]
LinAlgError(#[from] LinAlgError),
Expand Down
14 changes: 6 additions & 8 deletions src/joback.rs
Original file line number Diff line number Diff line change
Expand Up @@ -171,7 +171,7 @@ mod tests {
struct ModelRecord;

#[test]
fn paper_example() {
fn paper_example() -> EosResult<()> {
let segments_json = r#"[
{
"identifier": "-Cl",
Expand Down Expand Up @@ -226,8 +226,7 @@ mod tests {
],
None,
)
.segment_count(&segment_records)
.unwrap();
.segment_count(&segment_records)?;
assert_eq!(segments.get(&segment_records[0]), Some(&2.0));
assert_eq!(segments.get(&segment_records[1]), Some(&4.0));
assert_eq!(segments.get(&segment_records[2]), Some(&2.0));
Expand Down Expand Up @@ -264,17 +263,16 @@ mod tests {
1000.0 * KELVIN,
1.0 * ANGSTROM.powi(3),
&(arr1(&[1.0]) * MOL),
)
.unwrap();
)?;
assert!(
(state
.c_p(Contributions::IdealGas)
.to_reduced(JOULE / MOL / KELVIN)
.unwrap()
.to_reduced(JOULE / MOL / KELVIN)?
- 224.6)
.abs()
< 1.0
)
);
Ok(())
}

#[test]
Expand Down
46 changes: 20 additions & 26 deletions src/python/utils.rs
Original file line number Diff line number Diff line change
Expand Up @@ -138,15 +138,12 @@ macro_rules! impl_estimator {
target: &PySIArray1,
temperature: &PySIArray1,
std_parameters: Option<Vec<f64>>,
) -> Self {
Self(Rc::new(
VaporPressure::<SIUnit>::new(
target.clone().into(),
temperature.clone().into(),
std_parameters.unwrap_or(vec![0.0, 0.0, 0.0]),
)
.unwrap(),
))
) -> PyResult<Self> {
Ok(Self(Rc::new(VaporPressure::<SIUnit>::new(
target.clone().into(),
temperature.clone().into(),
std_parameters.unwrap_or(vec![0.0, 0.0, 0.0]),
)?)))
}

/// Create a DataSet with experimental data for liquid density.
Expand Down Expand Up @@ -177,15 +174,12 @@ macro_rules! impl_estimator {
target: &PySIArray1,
temperature: &PySIArray1,
pressure: &PySIArray1,
) -> Self {
Self(Rc::new(
LiquidDensity::<SIUnit>::new(
target.clone().into(),
temperature.clone().into(),
pressure.clone().into(),
)
.unwrap(),
))
) -> PyResult<Self> {
Ok(Self(Rc::new(LiquidDensity::<SIUnit>::new(
target.clone().into(),
temperature.clone().into(),
pressure.clone().into(),
)?)))
}

/// Create a DataSet with experimental data for liquid density
Expand All @@ -211,14 +205,14 @@ macro_rules! impl_estimator {
/// eos_python.saft.estimator.DataSet.relative_difference
#[staticmethod]
#[pyo3(text_signature = "(target, temperature)")]
fn equilibrium_liquid_density(target: &PySIArray1, temperature: &PySIArray1) -> Self {
Self(Rc::new(
EquilibriumLiquidDensity::<SIUnit>::new(
target.clone().into(),
temperature.clone().into(),
)
.unwrap(),
))
fn equilibrium_liquid_density(
target: &PySIArray1,
temperature: &PySIArray1,
) -> PyResult<Self> {
Ok(Self(Rc::new(EquilibriumLiquidDensity::<SIUnit>::new(
target.clone().into(),
temperature.clone().into(),
)?)))
}

/// Return `input` as ``Dict[str, SIArray1]``.
Expand Down
6 changes: 3 additions & 3 deletions src/state/builder.rs
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ use std::rc::Rc;
/// # use approx::assert_relative_eq;
/// # fn main() -> EosResult<()> {
/// // Create a state for given T,V,N
/// let eos = Rc::new(PengRobinson::new(Rc::new(PengRobinsonParameters::new_simple(&[369.8], &[41.9 * 1e5], &[0.15], &[15.0]).unwrap())));
/// let eos = Rc::new(PengRobinson::new(Rc::new(PengRobinsonParameters::new_simple(&[369.8], &[41.9 * 1e5], &[0.15], &[15.0])?)));
/// let state = StateBuilder::new(&eos)
/// .temperature(300.0 * KELVIN)
/// .volume(12.5 * METER.powi(3))
Expand All @@ -27,7 +27,7 @@ use std::rc::Rc;
/// assert_eq!(state.density, 0.2 * MOL / METER.powi(3));
///
/// // For a pure component, the composition does not need to be specified.
/// let eos = Rc::new(PengRobinson::new(Rc::new(PengRobinsonParameters::new_simple(&[369.8], &[41.9 * 1e5], &[0.15], &[15.0]).unwrap())));
/// let eos = Rc::new(PengRobinson::new(Rc::new(PengRobinsonParameters::new_simple(&[369.8], &[41.9 * 1e5], &[0.15], &[15.0])?)));
/// let state = StateBuilder::new(&eos)
/// .temperature(300.0 * KELVIN)
/// .volume(12.5 * METER.powi(3))
Expand All @@ -42,7 +42,7 @@ use std::rc::Rc;
/// &[41.9 * 1e5, 48.2 * 1e5],
/// &[0.15, 0.10],
/// &[15.0, 30.0]
/// ).unwrap())
/// )?)
/// ));
/// let state = StateBuilder::new(&eos)
/// .temperature(300.0 * KELVIN)
Expand Down
82 changes: 44 additions & 38 deletions src/utils/dataset.rs
Original file line number Diff line number Diff line change
Expand Up @@ -100,8 +100,7 @@ impl<U: EosUnit> VaporPressure<U> {
) -> Result<Self, FitError> {
let datapoints = target.len();
let max_temperature = temperature
.to_reduced(U::reference_temperature())
.unwrap()
.to_reduced(U::reference_temperature())?
.into_iter()
.reduce(|a, b| a.max(b))
.unwrap()
Expand Down Expand Up @@ -147,10 +146,13 @@ impl<U: EosUnit, E: EquationOfState> DataSet<U, E> for VaporPressure<U> {
where
QuantityScalar<U>: std::fmt::Display + std::fmt::LowerExp,
{
let tc =
State::critical_point(eos, None, Some(self.max_temperature), SolverOptions::default())
.unwrap()
.temperature;
let tc = State::critical_point(
eos,
None,
Some(self.max_temperature),
SolverOptions::default(),
)?
.temperature;

let unit = self.target.get(0);
let mut prediction = Array1::zeros(self.datapoints) * unit;
Expand All @@ -160,12 +162,12 @@ impl<U: EosUnit, E: EquationOfState> DataSet<U, E> for VaporPressure<U> {
if let Some(pvap) =
PhaseEquilibrium::vapor_pressure(eos, self.temperature.get(i))[0]
{
prediction.try_set(i, pvap).unwrap();
prediction.try_set(i, pvap)?;
} else {
prediction.try_set(i, f64::NAN * unit).unwrap();
prediction.try_set(i, f64::NAN * unit)?;
}
} else {
prediction.try_set(i, f64::NAN * unit).unwrap();
prediction.try_set(i, f64::NAN * unit)?;
}
}
Ok(prediction)
Expand All @@ -176,14 +178,18 @@ impl<U: EosUnit, E: EquationOfState> DataSet<U, E> for VaporPressure<U> {
QuantityScalar<U>: std::fmt::Display + std::fmt::LowerExp,
{
let tc_inv = 1.0
/ State::critical_point(eos, None, Some(self.max_temperature), SolverOptions::default())
.unwrap()
.temperature;

let reduced_temperatures = (0..self.datapoints)
.map(|i| (self.temperature.get(i) * tc_inv).into_value().unwrap())
/ State::critical_point(
eos,
None,
Some(self.max_temperature),
SolverOptions::default(),
)?
.temperature;

let reduced_temperatures: Result<_, _> = (0..self.datapoints)
.map(|i| (self.temperature.get(i) * tc_inv).into_value())
.collect();
let mut weights = self.weight_from_std(&reduced_temperatures);
let mut weights = self.weight_from_std(&reduced_temperatures?);
weights /= weights.sum();

let prediction = &self.predict(eos)?;
Expand All @@ -193,8 +199,7 @@ impl<U: EosUnit, E: EquationOfState> DataSet<U, E> for VaporPressure<U> {
cost[i] = weights[i]
* 5.0
* (self.temperature.get(i) - 1.0 / tc_inv)
.to_reduced(U::reference_temperature())
.unwrap();
.to_reduced(U::reference_temperature())?;
} else {
cost[i] = weights[i]
* ((self.target.get(i) - prediction.get(i)) / self.target.get(i))
Expand Down Expand Up @@ -275,9 +280,9 @@ impl<U: EosUnit, E: EquationOfState + MolarWeight<U>> DataSet<U, E> for LiquidDe
DensityInitialization::Liquid,
);
if let Ok(s) = state {
prediction.try_set(i, s.mass_density()).unwrap();
prediction.try_set(i, s.mass_density())?;
} else {
prediction.try_set(i, 1.0e10 * unit).unwrap();
prediction.try_set(i, 1.0e10 * unit)?;
}
}
Ok(prediction)
Expand Down Expand Up @@ -322,8 +327,7 @@ impl<U: EosUnit> EquilibriumLiquidDensity<U> {
) -> Result<Self, FitError> {
let datapoints = target.len();
let max_temperature = temperature
.to_reduced(U::reference_temperature())
.unwrap()
.to_reduced(U::reference_temperature())?
.into_iter()
.reduce(|a, b| a.max(b))
.unwrap()
Expand Down Expand Up @@ -362,23 +366,24 @@ impl<U: EosUnit, E: EquationOfState + MolarWeight<U>> DataSet<U, E>
where
QuantityScalar<U>: std::fmt::Display + std::fmt::LowerExp,
{
let tc =
State::critical_point(eos, None, Some(self.max_temperature), SolverOptions::default())
.unwrap()
.temperature;
let tc = State::critical_point(
eos,
None,
Some(self.max_temperature),
SolverOptions::default(),
)?
.temperature;

let unit = self.target.get(0);
let mut prediction = Array1::zeros(self.datapoints) * unit;
for i in 0..self.datapoints {
let t: QuantityScalar<U> = self.temperature.get(i);
if t < tc {
let state: PhaseEquilibrium<U, E, 2> =
PhaseEquilibrium::pure_t(eos, t, None, SolverOptions::default()).unwrap();
prediction
.try_set(i, state.liquid().mass_density())
.unwrap();
PhaseEquilibrium::pure_t(eos, t, None, SolverOptions::default())?;
prediction.try_set(i, state.liquid().mass_density())?;
} else {
prediction.try_set(i, f64::NAN * unit).unwrap();
prediction.try_set(i, f64::NAN * unit)?;
}
}
Ok(prediction)
Expand All @@ -388,20 +393,21 @@ impl<U: EosUnit, E: EquationOfState + MolarWeight<U>> DataSet<U, E>
where
QuantityScalar<U>: std::fmt::Display + std::fmt::LowerExp,
{
let tc =
State::critical_point(eos, None, Some(self.max_temperature), SolverOptions::default())
.unwrap()
.temperature;
let tc = State::critical_point(
eos,
None,
Some(self.max_temperature),
SolverOptions::default(),
)?
.temperature;
let n_inv = 1.0 / self.datapoints as f64;
let prediction = &self.predict(eos)?;
let mut cost = Array1::zeros(self.datapoints);
for i in 0..self.datapoints {
if prediction.get(i).is_nan() {
cost[i] = n_inv
* 5.0
* (self.temperature.get(i) - tc)
.to_reduced(U::reference_temperature())
.unwrap();
* (self.temperature.get(i) - tc).to_reduced(U::reference_temperature())?;
} else {
cost[i] = n_inv
* ((self.target.get(i) - prediction.get(i)) / self.target.get(i))
Expand Down
4 changes: 3 additions & 1 deletion src/utils/estimator.rs
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
//! optimization.
use super::dataset::*;
use crate::equation_of_state::EquationOfState;
use crate::EosUnit;
use crate::{EosError, EosUnit};
use ndarray::{arr1, concatenate, Array1, ArrayView1, Axis};
use quantity::{QuantityArray1, QuantityError, QuantityScalar};
use std::fmt;
Expand All @@ -25,6 +25,8 @@ pub enum FitError {
ParseError(#[from] ParseFloatError),
#[error(transparent)]
QuantityError(#[from] QuantityError),
#[error(transparent)]
EosError(#[from] EosError),
}

/// A collection of [`DataSet`]s and weights that can be used to
Expand Down