Migrate FeOs to new quantity crate

feos-org · prehner · Oct 28, 2024 · Sep 23, 2024 · Sep 23, 2024 · Sep 24, 2024
commit a1c975d06184cf69e85ca267793bf5aaec504847
diff --git a/feos-core/Cargo.toml b/feos-core/Cargo.toml
@@ -19,7 +19,7 @@ rustdoc-args = [ "--html-in-header", "./docs-header.html" ]
 features = [ "rayon" ]
 
 [dependencies]
-quantity = { version = "0.8", optional = true }
+quantity = { path = "../../../quantity" }
 num-dual = { version = "0.9", features = ["linalg"] }
 ndarray = { version = "0.15", features = ["serde", "approx-0_5"] }
 nalgebra = "0.32"

diff --git a/feos-core/src/cubic.rs b/feos-core/src/cubic.rs
@@ -6,10 +6,10 @@
 //! [this wikipedia article](https://en.wikipedia.org/wiki/Cubic_equations_of_state#Peng%E2%80%93Robinson_equation_of_state).
 use crate::equation_of_state::{Components, Residual};
 use crate::parameter::{Identifier, Parameter, ParameterError, PureRecord};
-use crate::si::{MolarWeight, GRAM, MOL};
 use crate::state::StateHD;
 use ndarray::{Array1, Array2, ScalarOperand};
 use num_dual::DualNum;
+use quantity::{MolarWeight, GRAM, MOL};
 use serde::{Deserialize, Serialize};
 use std::f64::consts::SQRT_2;
 use std::fmt;
@@ -223,10 +223,10 @@ impl Residual for PengRobinson {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::si::{KELVIN, PASCAL};
     use crate::state::{Contributions, State};
     use crate::{EosResult, SolverOptions, Verbosity};
     use approx::*;
+    use quantity::{KELVIN, PASCAL};
     use std::sync::Arc;
 
     fn pure_record_vec() -> Vec<PureRecord<PengRobinsonRecord>> {

diff --git a/feos-core/src/density_iteration.rs b/feos-core/src/density_iteration.rs
@@ -1,8 +1,9 @@
 use crate::equation_of_state::Residual;
 use crate::errors::{EosError, EosResult};
-use crate::si::{Density, Moles, Pressure, Temperature};
 use crate::state::State;
+use crate::ReferenceSystem;
 use ndarray::Array1;
+use quantity::{Density, Moles, Pressure, Temperature};
 use std::sync::Arc;
 
 pub fn density_iteration<E: Residual>(

diff --git a/feos-core/src/equation_of_state/mod.rs b/feos-core/src/equation_of_state/mod.rs
@@ -1,8 +1,8 @@
-use crate::{
-    si::{Diffusivity, MolarWeight, Moles, Temperature, ThermalConductivity, Viscosity, Volume},
-    EosResult,
-};
+use crate::EosResult;
 use ndarray::{Array1, ScalarOperand};
+use quantity::{
+    Diffusivity, MolarWeight, Moles, Temperature, ThermalConductivity, Viscosity, Volume,
+};
 use std::sync::Arc;
 
 mod ideal_gas;

diff --git a/feos-core/src/equation_of_state/residual.rs b/feos-core/src/equation_of_state/residual.rs
@@ -1,11 +1,10 @@
 use super::Components;
-use crate::si::*;
-use crate::StateHD;
-use crate::{EosError, EosResult};
+use crate::{EosError, EosResult, ReferenceSystem, StateHD};
 use ndarray::prelude::*;
 use ndarray::ScalarOperand;
 use num_dual::*;
 use num_traits::{One, Zero};
+use quantity::*;
 use std::ops::Div;
 
 /// A reisdual Helmholtz energy model.

diff --git a/feos-core/src/lib.rs b/feos-core/src/lib.rs
@@ -1,6 +1,9 @@
 #![warn(clippy::all)]
 #![allow(clippy::reversed_empty_ranges)]
 #![warn(clippy::allow_attributes)]
+use quantity::{Quantity, SIUnit};
+use std::ops::{Div, Mul};
+use typenum::Integer;
 
 /// Print messages with level `Verbosity::Iter` or higher.
 #[macro_export]
@@ -28,7 +31,7 @@ mod equation_of_state;
 mod errors;
 pub mod parameter;
 mod phase_equilibria;
-pub mod si;
+// pub mod si;
 mod state;
 pub use equation_of_state::{
     Components, EntropyScaling, EquationOfState, IdealGas, NoResidual, Residual,
@@ -115,18 +118,82 @@ impl SolverOptions {
     }
 }
 
+/// Reference values used for reduced properties in feos
+const REFERENCE_VALUES: [f64; 7] = [
+    1e-12,               // 1 ps
+    1e-10,               // 1 Å
+    1.380649e-27,        // Fixed through k_B
+    1.0,                 // 1 A
+    1.0,                 // 1 K
+    1.0 / 6.02214076e23, // 1/N_AV
+    1.0,                 // 1 Cd
+];
+
+pub trait ReferenceSystem<Inner> {
+    fn from_reduced(value: Inner) -> Self
+    where
+        Inner: Mul<f64, Output = Inner>;
+    fn to_reduced(&self) -> Inner
+    where
+        for<'a> &'a Inner: Div<f64, Output = Inner>;
+}
+
+/// Conversion to and from reduced units
+impl<
+        Inner,
+        T: Integer,
+        L: Integer,
+        M: Integer,
+        I: Integer,
+        THETA: Integer,
+        N: Integer,
+        J: Integer,
+    > ReferenceSystem<Inner> for Quantity<Inner, SIUnit<T, L, M, I, THETA, N, J>>
+{
+    fn from_reduced(value: Inner) -> Self
+    where
+        Inner: Mul<f64, Output = Inner>,
+    {
+        Self::new(
+            value
+                * REFERENCE_VALUES[0].powi(T::I32)
+                * REFERENCE_VALUES[1].powi(L::I32)
+                * REFERENCE_VALUES[2].powi(M::I32)
+                * REFERENCE_VALUES[3].powi(I::I32)
+                * REFERENCE_VALUES[4].powi(THETA::I32)
+                * REFERENCE_VALUES[5].powi(N::I32)
+                * REFERENCE_VALUES[6].powi(J::I32),
+        )
+    }
+
+    fn to_reduced(&self) -> Inner
+    where
+        for<'a> &'a Inner: Div<f64, Output = Inner>,
+    {
+        self.convert_to(Quantity::new(
+            REFERENCE_VALUES[0].powi(T::I32)
+                * REFERENCE_VALUES[1].powi(L::I32)
+                * REFERENCE_VALUES[2].powi(M::I32)
+                * REFERENCE_VALUES[3].powi(I::I32)
+                * REFERENCE_VALUES[4].powi(THETA::I32)
+                * REFERENCE_VALUES[5].powi(N::I32)
+                * REFERENCE_VALUES[6].powi(J::I32),
+        ))
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use crate::cubic::*;
     use crate::equation_of_state::{Components, EquationOfState, IdealGas};
     use crate::parameter::*;
-    use crate::si::{BAR, KELVIN, MOL, RGAS};
     use crate::Contributions;
     use crate::EosResult;
     use crate::StateBuilder;
     use approx::*;
     use ndarray::Array1;
     use num_dual::DualNum;
+    use quantity::{BAR, KELVIN, MOL, RGAS};
     use std::sync::Arc;
 
     // Only to be able to instantiate an `EquationOfState`

diff --git a/feos-core/src/phase_equilibria/bubble_dew.rs b/feos-core/src/phase_equilibria/bubble_dew.rs
@@ -1,16 +1,15 @@
 use super::PhaseEquilibrium;
 use crate::equation_of_state::Residual;
 use crate::errors::{EosError, EosResult};
-use crate::si::{Density, Dimensionless, Moles, Pressure, Quantity, SIUnit, Temperature, RGAS};
-use crate::state::TPSpec;
 use crate::state::{
     Contributions,
     DensityInitialization::{InitialDensity, Liquid, Vapor},
-    State, StateBuilder,
+    State, StateBuilder, TPSpec,
 };
-use crate::{SolverOptions, Verbosity};
+use crate::{ReferenceSystem, SolverOptions, Verbosity};
 use ndarray::*;
 use num_dual::linalg::{norm, LU};
+use quantity::{Density, Dimensionless, Moles, Pressure, Quantity, SIUnit, Temperature, RGAS};
 use std::fmt;
 use std::sync::Arc;
 use typenum::{N1, N2, P1, Z0};
@@ -323,7 +322,7 @@ impl TemperatureOrPressure for Quantity<f64, SIUnit<N2, N1, P1, Z0, Z0, Z0, Z0>>
         // Derivative w.r.t. temperature
         let ln_phi_1_dt = state1.dln_phi_dt();
         let ln_phi_2_dt = state2.dln_phi_dt();
-        let df = ((ln_phi_1_dt - ln_phi_2_dt) * Dimensionless::from(&state1.molefracs * &k)).sum();
+        let df = ((ln_phi_1_dt - ln_phi_2_dt) * Dimensionless::new(&state1.molefracs * &k)).sum();
         let mut tstep = -f / df;
 
         // catch too big t-steps
@@ -550,7 +549,7 @@ impl<E: Residual> PhaseEquilibrium<E, 2> {
         let p_l = liquid.pressure(Contributions::Total);
         let mu_l = liquid.residual_chemical_potential();
         let p_i = (liquid_molefracs * temperature * density * RGAS)
-            * Dimensionless::from(
+            * Dimensionless::new(
                 ((mu_l - p_l * v_l) / (RGAS * temperature))
                     .into_value()
                     .mapv(f64::exp),

diff --git a/feos-core/src/phase_equilibria/mod.rs b/feos-core/src/phase_equilibria/mod.rs
@@ -1,9 +1,9 @@
 use crate::equation_of_state::Residual;
 use crate::errors::{EosError, EosResult};
-use crate::si::{Dimensionless, Energy, Moles, Pressure, Temperature, RGAS};
 use crate::state::{DensityInitialization, State};
-use crate::Contributions;
+use crate::{Contributions, ReferenceSystem};
 use ndarray::Array1;
+use quantity::{Dimensionless, Energy, Moles, Pressure, Temperature, RGAS};
 use std::fmt;
 use std::fmt::Write;
 use std::sync::Arc;
@@ -193,7 +193,7 @@ impl<E: Residual, const N: usize> PhaseEquilibrium<E, N> {
             let ln_rho = s.partial_density.to_reduced().mapv(f64::ln);
             acc + s.residual_helmholtz_energy()
                 + s.pressure(Contributions::Total) * s.volume
-                + RGAS * s.temperature * (s.moles.clone() * Dimensionless::from(ln_rho - 1.0)).sum()
+                + RGAS * s.temperature * (s.moles.clone() * Dimensionless::new(ln_rho - 1.0)).sum()
         })
     }
 }

diff --git a/feos-core/src/phase_equilibria/phase_diagram_binary.rs b/feos-core/src/phase_equilibria/phase_diagram_binary.rs
@@ -2,11 +2,11 @@ use super::bubble_dew::TemperatureOrPressure;
 use super::{PhaseDiagram, PhaseEquilibrium};
 use crate::equation_of_state::Residual;
 use crate::errors::{EosError, EosResult};
-use crate::si::{Density, Moles, Pressure, Temperature, RGAS};
 use crate::state::{Contributions, DensityInitialization, State, StateBuilder, TPSpec};
-use crate::SolverOptions;
+use crate::{ReferenceSystem, SolverOptions};
 use ndarray::{arr1, arr2, concatenate, s, Array1, Array2, Axis};
 use num_dual::linalg::{norm, LU};
+use quantity::{Density, Moles, Pressure, Temperature, RGAS};
 use std::sync::Arc;
 
 const DEFAULT_POINTS: usize = 51;

diff --git a/feos-core/src/phase_equilibria/phase_diagram_pure.rs b/feos-core/src/phase_equilibria/phase_diagram_pure.rs
@@ -1,7 +1,7 @@
 use super::PhaseEquilibrium;
 use crate::equation_of_state::Residual;
 use crate::errors::EosResult;
-use crate::si::Temperature;
+use quantity::Temperature;
 use crate::state::{State, StateVec};
 use crate::SolverOptions;
 #[cfg(feature = "rayon")]

diff --git a/feos-core/src/phase_equilibria/phase_envelope.rs b/feos-core/src/phase_equilibria/phase_envelope.rs
@@ -1,7 +1,7 @@
 use super::{PhaseDiagram, PhaseEquilibrium};
 use crate::equation_of_state::Residual;
 use crate::errors::EosResult;
-use crate::si::{Moles, Pressure, Temperature};
+use quantity::{Moles, Pressure, Temperature};
 use crate::state::{Contributions, State};
 use crate::SolverOptions;
 use ndarray::Array1;

diff --git a/feos-core/src/phase_equilibria/stability_analysis.rs b/feos-core/src/phase_equilibria/stability_analysis.rs
@@ -1,12 +1,12 @@
 use super::PhaseEquilibrium;
 use crate::equation_of_state::Residual;
 use crate::errors::{EosError, EosResult};
-use crate::si::Moles;
 use crate::state::{Contributions, DensityInitialization, State};
-use crate::{SolverOptions, Verbosity};
+use crate::{ReferenceSystem, SolverOptions, Verbosity};
 use ndarray::*;
 use num_dual::linalg::smallest_ev;
 use num_dual::linalg::LU;
+use quantity::Moles;
 use std::ops::MulAssign;
 
 const X_DOMINANT: f64 = 0.99;

diff --git a/feos-core/src/phase_equilibria/tp_flash.rs b/feos-core/src/phase_equilibria/tp_flash.rs
@@ -1,11 +1,11 @@
 use super::PhaseEquilibrium;
 use crate::equation_of_state::Residual;
 use crate::errors::{EosError, EosResult};
-use crate::si::{Dimensionless, Moles, Pressure, Temperature};
 use crate::state::{Contributions, DensityInitialization, State};
 use crate::{SolverOptions, Verbosity};
 use ndarray::*;
 use num_dual::linalg::norm;
+use quantity::{Dimensionless, Moles, Pressure, Temperature};
 use std::sync::Arc;
 
 const MAX_ITER_TP: usize = 400;
@@ -315,9 +315,9 @@ impl<E: Residual> PhaseEquilibrium<E, 2> {
         beta = rachford_rice(&feed_state.molefracs, k, Some(beta))?;
 
         // update VLE
-        let v = feed_state.moles.clone() * Dimensionless::from(beta * k / (1.0 - beta + beta * k));
+        let v = feed_state.moles.clone() * Dimensionless::new(beta * k / (1.0 - beta + beta * k));
         let l =
-            feed_state.moles.clone() * Dimensionless::from((1.0 - beta) / (1.0 - beta + beta * k));
+            feed_state.moles.clone() * Dimensionless::new((1.0 - beta) / (1.0 - beta + beta * k));
         self.update_moles(feed_state.pressure(Contributions::Total), [&v, &l])?;
         Ok(())
     }

diff --git a/feos-core/src/phase_equilibria/vle_pure.rs b/feos-core/src/phase_equilibria/vle_pure.rs
@@ -1,10 +1,10 @@
 use super::PhaseEquilibrium;
 use crate::equation_of_state::Residual;
 use crate::errors::{EosError, EosResult};
-use crate::si::{Moles, Pressure, Temperature, RGAS};
 use crate::state::{Contributions, DensityInitialization, State, TPSpec};
-use crate::{SolverOptions, TemperatureOrPressure, Verbosity};
+use crate::{ReferenceSystem, SolverOptions, TemperatureOrPressure, Verbosity};
 use ndarray::{arr1, Array1};
+use quantity::{Moles, Pressure, Temperature, RGAS};
 use std::sync::Arc;
 
 const SCALE_T_NEW: f64 = 0.7;

diff --git a/feos-core/src/python/user_defined.rs b/feos-core/src/python/user_defined.rs
@@ -1,5 +1,5 @@
 #![allow(non_snake_case)]
-use crate::si::MolarWeight;
+use quantity::MolarWeight;
 use crate::{Components, IdealGas, Residual, StateHD};
 use ndarray::{Array1, ScalarOperand};
 use num_dual::*;