Skip to content

Remove unit generics #115

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.

Sign up for GitHub

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

Merged
prehner merged 5 commits into from
Jan 20, 2023
Merged

Remove unit generics #115

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
7b545d7
removed unit generics from feos-core
g-bauer Jan 18, 2023
58778ae
removed unit generics from feos-core python and fixed formats
g-bauer Jan 18, 2023
dd4872d
removed unit generics from feos-dft
g-bauer Jan 18, 2023
5e94f76
removed unit generics from feos, fixes in core and dft in python macros
g-bauer Jan 19, 2023
e4e7b8c
ran cargo fmt
g-bauer Jan 19, 2023
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
Next Next commit
removed unit generics from feos-core
  • Loading branch information
@g-bauer
g-bauer committed Jan 20, 2023
commit 7b545d73fee6364abe44d559e8a7742de7bfb41c
4 changes: 2 additions & 2 deletions feos-core/src/cubic.rs
View file
Open in desktop
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
40 changes: 20 additions & 20 deletions feos-core/src/density_iteration.rs
View file
Open in desktop
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::{SINumber, SIArray1, 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,8 @@ 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 +141,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 +174,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
View file
Open in desktop
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>;
}
12 changes: 6 additions & 6 deletions feos-core/src/joback.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ use crate::{
use conv::ValueInto;
use ndarray::Array1;
use num_dual::*;
use quantity::QuantityScalar;
use quantity::si::{SINumber, SIUnit};
use serde::{Deserialize, Serialize};
use std::fmt;

Expand Down Expand Up @@ -85,17 +85,17 @@ impl Joback {
}

/// Directly calculates the ideal gas heat capacity from the Joback model.
pub fn c_p<U: EosUnit>(
pub fn c_p(
&self,
temperature: QuantityScalar<U>,
temperature: SINumber,
molefracs: &Array1<f64>,
) -> EosResult<QuantityScalar<U>> {
let t = temperature.to_reduced(U::reference_temperature())?;
) -> EosResult<SINumber> {
let t = temperature.to_reduced(SIUnit::reference_temperature())?;
let mut c_p = 0.0;
for (j, &x) in self.records.iter().zip(molefracs.iter()) {
c_p += x * (j.a + j.b * t + j.c * t.powi(2) + j.d * t.powi(3) + j.e * t.powi(4));
}
Ok(c_p / RGAS * U::gas_constant())
Ok(c_p / RGAS * SIUnit::gas_constant())
}
}

Expand Down
Loading