Skip to content
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.

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

Remove unit generics #115

Merged
merged 5 commits into from
Jan 20, 2023
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
1 change: 1 addition & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
- Export `EosVariant` and `FunctionalVariant` directly in the crate root instead of their own modules. [#62](https://github.com/feos-org/feos/pull/62)
- Changed constructors `VaporPressure::new` and `DataSet.vapor_pressure` (Python) to take a new optional argument `critical_temperature`. [#86](https://github.com/feos-org/feos/pull/86)
- The limitations of the homo gc method for PC-SAFT are enforced more strictly. [#88](https://github.com/feos-org/feos/pull/88)
- Removed generics for units in all structs and traits in favor of static SI units. [#115](https://github.com/feos-org/feos/pull/115)

## [0.3.0] - 2022-09-14
- Major restructuring of the entire `feos` project. All individual models are reunited in the `feos` crate. `feos-core` and `feos-dft` still live as individual crates within the `feos` workspace.
Expand Down
8 changes: 2 additions & 6 deletions benches/dual_numbers.rs
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ use std::sync::Arc;
/// - temperature is 80% of critical temperature,
/// - volume is critical volume,
/// - molefracs (or moles) for equimolar mixture.
fn state_pcsaft(parameters: PcSaftParameters) -> State<SIUnit, PcSaft> {
fn state_pcsaft(parameters: PcSaftParameters) -> State<PcSaft> {
let n = parameters.pure_records.len();
let eos = Arc::new(PcSaft::new(Arc::new(parameters)));
let moles = Array::from_elem(n, 1.0 / n as f64) * 10.0 * MOL;
Expand All @@ -36,11 +36,7 @@ where
}

/// Benchmark for evaluation of the Helmholtz energy for different dual number types.
fn bench_dual_numbers<E: EquationOfState>(
c: &mut Criterion,
group_name: &str,
state: State<SIUnit, E>,
) {
fn bench_dual_numbers<E: EquationOfState>(c: &mut Criterion, group_name: &str, state: State<E>) {
let mut group = c.benchmark_group(group_name);
group.bench_function("a_f64", |b| {
b.iter(|| a_res((&state.eos, &state.derive0())))
Expand Down
2 changes: 1 addition & 1 deletion benches/state_creation.rs
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ fn npt<E: EquationOfState>(
SINumber,
SINumber,
&SIArray1,
DensityInitialization<SIUnit>,
DensityInitialization,
),
) {
State::new_npt(eos, t, p, n, rho0).unwrap();
Expand Down
6 changes: 3 additions & 3 deletions benches/state_properties.rs
Original file line number Diff line number Diff line change
Expand Up @@ -8,11 +8,11 @@ use ndarray::arr1;
use quantity::si::*;
use std::sync::Arc;

type S = State<SIUnit, PcSaft>;
type S = State<PcSaft>;

/// Evaluate a property of a state given the EoS, the property to compute,
/// temperature, volume, moles, and the contributions to consider.
fn property<E: EquationOfState, T, F: Fn(&State<SIUnit, E>, Contributions) -> T>(
fn property<E: EquationOfState, T, F: Fn(&State<E>, Contributions) -> T>(
(eos, property, t, v, n, contributions): (
&Arc<E>,
F,
Expand All @@ -28,7 +28,7 @@ fn property<E: EquationOfState, T, F: Fn(&State<SIUnit, E>, Contributions) -> T>

/// Evaluate a property with of a state given the EoS, the property to compute,
/// temperature, volume, moles.
fn property_no_contributions<E: EquationOfState, T, F: Fn(&State<SIUnit, E>) -> T>(
fn property_no_contributions<E: EquationOfState, T, F: Fn(&State<E>) -> T>(
(eos, property, t, v, n): (&Arc<E>, F, SINumber, SINumber, &SIArray1),
) -> T {
let state = State::new_nvt(eos, t, v, n).unwrap();
Expand Down
1 change: 1 addition & 0 deletions feos-core/CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
- Added `StateVec::moles` getter. [#113](https://github.com/feos-org/feos/pull/113)
- Added public constructors `PhaseDiagram::new` and `StateVec::new` that allow the creation of the respective structs from a list of `PhaseEquilibrium`s or `State`s in Rust and Python. [#113](https://github.com/feos-org/feos/pull/113)
- Added new variant `EosError::Error` which allows dispatching generic errors that are not covered by the existing variants. [#98](https://github.com/feos-org/feos/pull/98)
- Removed generics for units in all structs and traits in favor of static SI units. [#115](https://github.com/feos-org/feos/pull/115)

### Changed
- Added `Sync` and `Send` as supertraits to `EquationOfState`. [#57](https://github.com/feos-org/feos/pull/57)
Expand Down
4 changes: 2 additions & 2 deletions feos-core/src/cubic.rs
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ use crate::state::StateHD;
use crate::MolarWeight;
use ndarray::{Array1, Array2};
use num_dual::DualNum;
use quantity::si::{SIArray1, SIUnit};
use quantity::si::SIArray1;
use serde::{Deserialize, Serialize};
use std::f64::consts::SQRT_2;
use std::fmt;
Expand Down Expand Up @@ -255,7 +255,7 @@ impl EquationOfState for PengRobinson {
}
}

impl MolarWeight<SIUnit> for PengRobinson {
impl MolarWeight for PengRobinson {
fn molar_weight(&self) -> SIArray1 {
self.parameters.molarweight.clone() * GRAM / MOL
}
Expand Down
43 changes: 23 additions & 20 deletions feos-core/src/density_iteration.rs
Original file line number Diff line number Diff line change
Expand Up @@ -2,26 +2,26 @@ use crate::equation_of_state::EquationOfState;
use crate::errors::{EosError, EosResult};
use crate::state::State;
use crate::EosUnit;
use quantity::{QuantityArray1, QuantityScalar};
use quantity::si::{SIArray1, SINumber, SIUnit};
use std::sync::Arc;

pub fn density_iteration<U: EosUnit, E: EquationOfState>(
pub fn density_iteration<E: EquationOfState>(
eos: &Arc<E>,
temperature: QuantityScalar<U>,
pressure: QuantityScalar<U>,
moles: &QuantityArray1<U>,
initial_density: QuantityScalar<U>,
) -> EosResult<State<U, E>> {
temperature: SINumber,
pressure: SINumber,
moles: &SIArray1,
initial_density: SINumber,
) -> EosResult<State<E>> {
let maxdensity = eos.max_density(Some(moles))?;
let (abstol, reltol) = (1e-12, 1e-14);
let n = moles.sum();

let mut rho = initial_density;
if rho <= 0.0 * U::reference_density() {
if rho <= 0.0 * SIUnit::reference_density() {
return Err(EosError::InvalidState(
String::from("density iteration"),
String::from("density"),
rho.to_reduced(U::reference_density())?,
rho.to_reduced(SIUnit::reference_density())?,
));
}

Expand Down Expand Up @@ -117,7 +117,7 @@ pub fn density_iteration<U: EosUnit, E: EquationOfState>(
} else {
rho = (rho + initial_density) * 0.5;
if (rho - initial_density)
.to_reduced(U::reference_density())?
.to_reduced(SIUnit::reference_density())?
.abs()
< 1e-8
{
Expand All @@ -128,8 +128,11 @@ pub fn density_iteration<U: EosUnit, E: EquationOfState>(
}
// Newton step
rho += delta_rho;
if error.to_reduced(U::reference_pressure())?.abs()
< f64::max(abstol, (rho * reltol).to_reduced(U::reference_density())?)
if error.to_reduced(SIUnit::reference_pressure())?.abs()
< f64::max(
abstol,
(rho * reltol).to_reduced(SIUnit::reference_density())?,
)
{
break 'iteration;
}
Expand All @@ -141,24 +144,24 @@ pub fn density_iteration<U: EosUnit, E: EquationOfState>(
}
}

fn pressure_spinodal<U: EosUnit, E: EquationOfState>(
fn pressure_spinodal<E: EquationOfState>(
eos: &Arc<E>,
temperature: QuantityScalar<U>,
rho_init: QuantityScalar<U>,
moles: &QuantityArray1<U>,
) -> EosResult<[QuantityScalar<U>; 2]> {
temperature: SINumber,
rho_init: SINumber,
moles: &SIArray1,
) -> EosResult<[SINumber; 2]> {
let maxiter = 30;
let abstol = 1e-8;

let maxdensity = eos.max_density(Some(moles))?;
let n = moles.sum();
let mut rho = rho_init;

if rho <= 0.0 * U::reference_density() {
if rho <= 0.0 * SIUnit::reference_density() {
return Err(EosError::InvalidState(
String::from("pressure spinodal"),
String::from("density"),
rho.to_reduced(U::reference_density())?,
rho.to_reduced(SIUnit::reference_density())?,
));
}

Expand All @@ -174,7 +177,7 @@ fn pressure_spinodal<U: EosUnit, E: EquationOfState>(
rho += delta_rho;

if dpdrho
.to_reduced(U::reference_pressure() / U::reference_density())?
.to_reduced(SIUnit::reference_pressure() / SIUnit::reference_density())?
.abs()
< abstol
{
Expand Down
96 changes: 45 additions & 51 deletions feos-core/src/equation_of_state.rs
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ use num_dual::{
Dual, Dual2_64, Dual3, Dual3_64, Dual64, DualNum, DualVec64, HyperDual, HyperDual64,
};
use num_traits::{One, Zero};
use quantity::{QuantityArray1, QuantityScalar};
use quantity::si::{SIArray1, SINumber, SIUnit};
use std::fmt;

/// Individual Helmholtz energy contribution that can
Expand Down Expand Up @@ -149,8 +149,8 @@ impl fmt::Display for DefaultIdealGasContribution {
///
/// The trait is required to be able to calculate (mass)
/// specific properties.
pub trait MolarWeight<U: EosUnit> {
fn molar_weight(&self) -> QuantityArray1<U>;
pub trait MolarWeight {
fn molar_weight(&self) -> SIArray1;
}

/// A general equation of state.
Expand Down Expand Up @@ -219,15 +219,12 @@ pub trait EquationOfState: Send + Sync {
/// of components of the equation of state. For a pure component, however,
/// no moles need to be provided. In that case, it is set to the constant
/// reference value.
fn validate_moles<U: EosUnit>(
&self,
moles: Option<&QuantityArray1<U>>,
) -> EosResult<QuantityArray1<U>> {
fn validate_moles(&self, moles: Option<&SIArray1>) -> EosResult<SIArray1> {
let l = moles.map_or(1, |m| m.len());
if self.components() == l {
match moles {
Some(m) => Ok(m.to_owned()),
None => Ok(Array::ones(1) * U::reference_moles()),
None => Ok(Array::ones(1) * SIUnit::reference_moles()),
}
} else {
Err(EosError::IncompatibleComponents(self.components(), l))
Expand All @@ -240,105 +237,102 @@ pub trait EquationOfState: Send + Sync {
/// equilibria and other iterations. It is not explicitly meant to
/// be a mathematical limit for the density (if those exist in the
/// equation of state anyways).
fn max_density<U: EosUnit>(
&self,
moles: Option<&QuantityArray1<U>>,
) -> EosResult<QuantityScalar<U>> {
fn max_density(&self, moles: Option<&SIArray1>) -> EosResult<SINumber> {
let mr = self
.validate_moles(moles)?
.to_reduced(U::reference_moles())?;
Ok(self.compute_max_density(&mr) * U::reference_density())
.to_reduced(SIUnit::reference_moles())?;
Ok(self.compute_max_density(&mr) * SIUnit::reference_density())
}

/// Calculate the second virial coefficient $B(T)$
fn second_virial_coefficient<U: EosUnit>(
fn second_virial_coefficient(
&self,
temperature: QuantityScalar<U>,
moles: Option<&QuantityArray1<U>>,
) -> EosResult<QuantityScalar<U>> {
temperature: SINumber,
moles: Option<&SIArray1>,
) -> EosResult<SINumber> {
let mr = self.validate_moles(moles)?;
let x = mr.to_reduced(mr.sum())?;
let mut rho = HyperDual64::zero();
rho.eps1[0] = 1.0;
rho.eps2[0] = 1.0;
let t = HyperDual64::from(temperature.to_reduced(U::reference_temperature())?);
let t = HyperDual64::from(temperature.to_reduced(SIUnit::reference_temperature())?);
let s = StateHD::new_virial(t, rho, x);
Ok(self.evaluate_residual(&s).eps1eps2[(0, 0)] * 0.5 / U::reference_density())
Ok(self.evaluate_residual(&s).eps1eps2[(0, 0)] * 0.5 / SIUnit::reference_density())
}

/// Calculate the third virial coefficient $C(T)$
fn third_virial_coefficient<U: EosUnit>(
fn third_virial_coefficient(
&self,
temperature: QuantityScalar<U>,
moles: Option<&QuantityArray1<U>>,
) -> EosResult<QuantityScalar<U>> {
temperature: SINumber,
moles: Option<&SIArray1>,
) -> EosResult<SINumber> {
let mr = self.validate_moles(moles)?;
let x = mr.to_reduced(mr.sum())?;
let rho = Dual3_64::zero().derive();
let t = Dual3_64::from(temperature.to_reduced(U::reference_temperature())?);
let t = Dual3_64::from(temperature.to_reduced(SIUnit::reference_temperature())?);
let s = StateHD::new_virial(t, rho, x);
Ok(self.evaluate_residual(&s).v3 / 3.0 / U::reference_density().powi(2))
Ok(self.evaluate_residual(&s).v3 / 3.0 / SIUnit::reference_density().powi(2))
}

/// Calculate the temperature derivative of the second virial coefficient $B'(T)$
fn second_virial_coefficient_temperature_derivative<U: EosUnit>(
fn second_virial_coefficient_temperature_derivative(
&self,
temperature: QuantityScalar<U>,
moles: Option<&QuantityArray1<U>>,
) -> EosResult<QuantityScalar<U>> {
temperature: SINumber,
moles: Option<&SIArray1>,
) -> EosResult<SINumber> {
let mr = self.validate_moles(moles)?;
let x = mr.to_reduced(mr.sum())?;
let mut rho = HyperDual::zero();
rho.eps1[0] = Dual64::one();
rho.eps2[0] = Dual64::one();
let t = HyperDual::from_re(
Dual64::from(temperature.to_reduced(U::reference_temperature())?).derive(),
Dual64::from(temperature.to_reduced(SIUnit::reference_temperature())?).derive(),
);
let s = StateHD::new_virial(t, rho, x);
Ok(self.evaluate_residual(&s).eps1eps2[(0, 0)].eps[0] * 0.5
/ (U::reference_density() * U::reference_temperature()))
/ (SIUnit::reference_density() * SIUnit::reference_temperature()))
}

/// Calculate the temperature derivative of the third virial coefficient $C'(T)$
fn third_virial_coefficient_temperature_derivative<U: EosUnit>(
fn third_virial_coefficient_temperature_derivative(
&self,
temperature: QuantityScalar<U>,
moles: Option<&QuantityArray1<U>>,
) -> EosResult<QuantityScalar<U>> {
temperature: SINumber,
moles: Option<&SIArray1>,
) -> EosResult<SINumber> {
let mr = self.validate_moles(moles)?;
let x = mr.to_reduced(mr.sum())?;
let rho = Dual3::zero().derive();
let t = Dual3::from_re(
Dual64::from(temperature.to_reduced(U::reference_temperature())?).derive(),
Dual64::from(temperature.to_reduced(SIUnit::reference_temperature())?).derive(),
);
let s = StateHD::new_virial(t, rho, x);
Ok(self.evaluate_residual(&s).v3.eps[0]
/ 3.0
/ (U::reference_density().powi(2) * U::reference_temperature()))
/ (SIUnit::reference_density().powi(2) * SIUnit::reference_temperature()))
}
}

/// Reference values and residual entropy correlations for entropy scaling.
pub trait EntropyScaling<U: EosUnit> {
pub trait EntropyScaling {
fn viscosity_reference(
&self,
temperature: QuantityScalar<U>,
volume: QuantityScalar<U>,
moles: &QuantityArray1<U>,
) -> EosResult<QuantityScalar<U>>;
temperature: SINumber,
volume: SINumber,
moles: &SIArray1,
) -> EosResult<SINumber>;
fn viscosity_correlation(&self, s_res: f64, x: &Array1<f64>) -> EosResult<f64>;
fn diffusion_reference(
&self,
temperature: QuantityScalar<U>,
volume: QuantityScalar<U>,
moles: &QuantityArray1<U>,
) -> EosResult<QuantityScalar<U>>;
temperature: SINumber,
volume: SINumber,
moles: &SIArray1,
) -> EosResult<SINumber>;
fn diffusion_correlation(&self, s_res: f64, x: &Array1<f64>) -> EosResult<f64>;
fn thermal_conductivity_reference(
&self,
temperature: QuantityScalar<U>,
volume: QuantityScalar<U>,
moles: &QuantityArray1<U>,
) -> EosResult<QuantityScalar<U>>;
temperature: SINumber,
volume: SINumber,
moles: &SIArray1,
) -> EosResult<SINumber>;
fn thermal_conductivity_correlation(&self, s_res: f64, x: &Array1<f64>) -> EosResult<f64>;
}
Loading