Compile-time checked units (#181)

Cargo.toml

@@ -22,7 +22,7 @@ members = ["feos-core", "feos-dft", "feos-derive"]
 crate-type = ["rlib", "cdylib"]
 
 [dependencies]
-quantity = "0.6"
+quantity = { version = "0.6", optional = true }
 num-dual = "0.7"
 feos-core = { version = "0.4", path = "feos-core" }
 feos-dft = { version = "0.4", path = "feos-dft", optional = true }
@@ -39,6 +39,7 @@ lazy_static = { version = "1.4", optional = true }
 indexmap = "1.8"
 rayon = { version = "1.5", optional = true }
 itertools = "0.10"
+typenum = "1.16"
 
 [dependencies.pyo3]
 version = "0.18"
@@ -64,7 +65,7 @@ uvtheory = ["lazy_static"]
 pets = []
 saftvrqmie = []
 rayon = ["dep:rayon", "ndarray/rayon", "feos-core/rayon", "feos-dft?/rayon"]
-python = ["pyo3", "numpy", "feos-core/python", "feos-dft?/python", "rayon"]
+python = ["pyo3", "numpy", "quantity/python", "feos-core/python", "feos-dft?/python", "rayon"]
 all_models = ["dft", "estimator", "pcsaft", "gc_pcsaft", "uvtheory", "pets", "saftvrqmie"]
 
 [[bench]]

benches/contributions.rs

@@ -8,9 +8,9 @@
 use criterion::{criterion_group, criterion_main, Criterion};
 use feos::pcsaft::{PcSaft, PcSaftParameters};
 use feos_core::parameter::{IdentifierOption, Parameter};
+use feos_core::si::*;
 use feos_core::{DensityInitialization, Derivative, Residual, State};
 use ndarray::arr1;
-use quantity::si::*;
 use std::sync::Arc;
 
 /// Benchmark for the PC-SAFT equation of state

benches/dft_pore.rs

@@ -5,12 +5,13 @@ use feos::gc_pcsaft::{GcPcSaftFunctional, GcPcSaftFunctionalParameters};
 use feos::hard_sphere::{FMTFunctional, FMTVersion};
 use feos::pcsaft::{PcSaftFunctional, PcSaftParameters};
 use feos_core::parameter::{IdentifierOption, Parameter, ParameterHetero};
+use feos_core::si::{ANGSTROM, KELVIN, NAV};
 use feos_core::{PhaseEquilibrium, State, StateBuilder};
 use feos_dft::adsorption::{ExternalPotential, Pore1D, PoreSpecification};
 use feos_dft::{DFTSolver, Geometry};
 use ndarray::arr1;
-use quantity::si::{ANGSTROM, KELVIN, NAV};
 use std::sync::Arc;
+use typenum::P3;
 
 fn fmt(c: &mut Criterion) {
     let mut group = c.benchmark_group("DFT_pore_fmt");
@@ -23,9 +24,9 @@ fn fmt(c: &mut Criterion) {
         None,
         None,
     );
-    let bulk = State::new_pure(&func, KELVIN, 0.75 / NAV / ANGSTROM.powi(3)).unwrap();
+    let bulk = State::new_pure(&func, KELVIN, 0.75 / NAV / ANGSTROM.powi::<P3>()).unwrap();
     group.bench_function("liquid", |b| {
-        b.iter(|| pore.initialize(&bulk, None, None).unwrap().solve(None))
+        b.iter(|| pore.initialize(&bulk, None, None).solve(None))
     });
 }
 
@@ -53,11 +54,11 @@ fn pcsaft(c: &mut Criterion) {
     let vle = PhaseEquilibrium::pure(&func, 300.0 * KELVIN, None, Default::default()).unwrap();
     let bulk = vle.liquid();
     group.bench_function("butane_liquid", |b| {
-        b.iter(|| pore.initialize(bulk, None, None).unwrap().solve(None))
+        b.iter(|| pore.initialize(bulk, None, None).solve(None))
     });
     let bulk = State::new_pure(&func, 300.0 * KELVIN, vle.vapor().density * 0.2).unwrap();
     group.bench_function("butane_vapor", |b| {
-        b.iter(|| pore.initialize(&bulk, None, None).unwrap().solve(None))
+        b.iter(|| pore.initialize(&bulk, None, None).solve(None))
     });
 
     let parameters = PcSaftParameters::from_json(
@@ -79,15 +80,15 @@ fn pcsaft(c: &mut Criterion) {
     .unwrap();
     let bulk = vle.liquid();
     group.bench_function("butane_pentane_liquid", |b| {
-        b.iter(|| pore.initialize(bulk, None, None).unwrap().solve(None))
+        b.iter(|| pore.initialize(bulk, None, None).solve(None))
     });
     let bulk = StateBuilder::new(&func)
         .temperature(300.0 * KELVIN)
         .partial_density(&(&vle.vapor().partial_density * 0.2))
         .build()
         .unwrap();
     group.bench_function("butane_pentane_vapor", |b| {
-        b.iter(|| pore.initialize(&bulk, None, None).unwrap().solve(None))
+        b.iter(|| pore.initialize(&bulk, None, None).solve(None))
     });
 }
 
@@ -121,11 +122,7 @@ fn gc_pcsaft(c: &mut Criterion) {
         .picard_iteration(None, None, Some(1e-5), None)
         .anderson_mixing(None, None, None, None, None);
     group.bench_function("butane_liquid", |b| {
-        b.iter(|| {
-            pore.initialize(bulk, None, None)
-                .unwrap()
-                .solve(Some(&solver))
-        })
+        b.iter(|| pore.initialize(bulk, None, None).solve(Some(&solver)))
     });
 }
 

benches/dual_numbers.rs

@@ -5,14 +5,15 @@
 //! creation.
 use criterion::{criterion_group, criterion_main, Criterion};
 use feos::pcsaft::{PcSaft, PcSaftParameters};
+use feos_core::si::*;
 use feos_core::{
     parameter::{IdentifierOption, Parameter},
     Derivative, HelmholtzEnergy, HelmholtzEnergyDual, Residual, State, StateHD,
 };
 use ndarray::{arr1, Array};
 use num_dual::DualNum;
-use quantity::si::*;
 use std::sync::Arc;
+use typenum::P3;
 
 /// Helper function to create a state for given parameters.
 /// - temperature is 80% of critical temperature,
@@ -123,7 +124,7 @@ fn methane_co2_pcsaft(c: &mut Criterion) {
 
     // 230 K, 50 bar, x0 = 0.15
     let temperature = 230.0 * KELVIN;
-    let density = 24.16896 * KILO * MOL / METER.powi(3);
+    let density = 24.16896 * KILO * MOL / METER.powi::<P3>();
     let volume = 10.0 * MOL / density;
     let x = arr1(&[0.15, 0.85]);
     let moles = &x * 10.0 * MOL;

benches/state_creation.rs

@@ -1,47 +1,61 @@
+#![allow(clippy::type_complexity)]
 use criterion::{criterion_group, criterion_main, Criterion};
 use feos::pcsaft::{PcSaft, PcSaftParameters};
+use feos_core::si::*;
 use feos_core::{
     parameter::{IdentifierOption, Parameter},
-    Contributions, DensityInitialization, PhaseEquilibrium, Residual, State,
+    Contributions, DensityInitialization, PhaseEquilibrium, Residual, State, TPSpec,
 };
 use ndarray::{Array, Array1};
-use quantity::si::*;
 use std::sync::Arc;
 
 /// Evaluate NPT constructor
 fn npt<E: Residual>(
     (eos, t, p, n, rho0): (
         &Arc<E>,
-        SINumber,
-        SINumber,
-        &SIArray1,
+        Temperature<f64>,
+        Pressure<f64>,
+        &Moles<Array1<f64>>,
         DensityInitialization,
     ),
 ) {
     State::new_npt(eos, t, p, n, rho0).unwrap();
 }
 
 /// Evaluate critical point constructor
-fn critical_point<E: Residual>((eos, n): (&Arc<E>, Option<&SIArray1>)) {
+fn critical_point<E: Residual>((eos, n): (&Arc<E>, Option<&Moles<Array1<f64>>>)) {
     State::critical_point(eos, n, None, Default::default()).unwrap();
 }
 
 /// Evaluate critical point constructor for binary systems at given T or p
-fn critical_point_binary<E: Residual>((eos, tp): (&Arc<E>, SINumber)) {
+fn critical_point_binary<E: Residual, TP>((eos, tp): (&Arc<E>, TP))
+where
+    TPSpec: From<TP>,
+{
     State::critical_point_binary(eos, tp, None, None, Default::default()).unwrap();
 }
 
 /// VLE for pure substance for given temperature or pressure
-fn pure<E: Residual>((eos, t_or_p): (&Arc<E>, SINumber)) {
+fn pure<E: Residual, TP>((eos, t_or_p): (&Arc<E>, TP))
+where
+    TPSpec: From<TP>,
+{
     PhaseEquilibrium::pure(eos, t_or_p, None, Default::default()).unwrap();
 }
 
 /// Evaluate temperature, pressure flash.
-fn tp_flash<E: Residual>((eos, t, p, feed): (&Arc<E>, SINumber, SINumber, &SIArray1)) {
+fn tp_flash<E: Residual>(
+    (eos, t, p, feed): (
+        &Arc<E>,
+        Temperature<f64>,
+        Pressure<f64>,
+        &Moles<Array1<f64>>,
+    ),
+) {
     PhaseEquilibrium::tp_flash(eos, t, p, feed, None, Default::default(), None).unwrap();
 }
 
-fn bubble_point<E: Residual>((eos, t, x): (&Arc<E>, SINumber, &Array1<f64>)) {
+fn bubble_point<E: Residual>((eos, t, x): (&Arc<E>, Temperature<f64>, &Array1<f64>)) {
     PhaseEquilibrium::bubble_point(
         eos,
         t,
@@ -53,7 +67,7 @@ fn bubble_point<E: Residual>((eos, t, x): (&Arc<E>, SINumber, &Array1<f64>)) {
     .unwrap();
 }
 
-fn dew_point<E: Residual>((eos, t, y): (&Arc<E>, SINumber, &Array1<f64>)) {
+fn dew_point<E: Residual>((eos, t, y): (&Arc<E>, Temperature<f64>, &Array1<f64>)) {
     PhaseEquilibrium::dew_point(
         eos,
         t,

benches/state_properties.rs

@@ -1,12 +1,14 @@
+#![allow(clippy::type_complexity)]
 use criterion::{criterion_group, criterion_main, Criterion};
 use feos::pcsaft::{PcSaft, PcSaftParameters};
+use feos_core::si::*;
 use feos_core::{
     parameter::{IdentifierOption, Parameter},
     Contributions, Residual, State,
 };
-use ndarray::arr1;
-use quantity::si::*;
+use ndarray::{arr1, Array1};
 use std::sync::Arc;
+use typenum::P3;
 
 type S = State<PcSaft>;
 
@@ -16,9 +18,9 @@ fn property<E: Residual, T, F: Fn(&State<E>, Contributions) -> T>(
     (eos, property, t, v, n, contributions): (
         &Arc<E>,
         F,
-        SINumber,
-        SINumber,
-        &SIArray1,
+        Temperature<f64>,
+        Volume<f64>,
+        &Moles<Array1<f64>>,
         Contributions,
     ),
 ) -> T {
@@ -29,7 +31,13 @@ fn property<E: Residual, T, F: Fn(&State<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: Residual, T, F: Fn(&State<E>) -> T>(
-    (eos, property, t, v, n): (&Arc<E>, F, SINumber, SINumber, &SIArray1),
+    (eos, property, t, v, n): (
+        &Arc<E>,
+        F,
+        Temperature<f64>,
+        Volume<f64>,
+        &Moles<Array1<f64>>,
+    ),
 ) -> T {
     let state = State::new_nvt(eos, t, v, n).unwrap();
     property(&state)
@@ -45,7 +53,7 @@ fn properties_pcsaft(c: &mut Criterion) {
     .unwrap();
     let eos = Arc::new(PcSaft::new(Arc::new(parameters)));
     let t = 300.0 * KELVIN;
-    let density = 71.18 * KILO * MOL / METER.powi(3);
+    let density = 71.18 * KILO * MOL / METER.powi::<P3>();
     let v = 100.0 * MOL / density;
     let x = arr1(&[1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0]);
     let m = &x * 100.0 * MOL;
@@ -80,7 +88,7 @@ fn properties_pcsaft_polar(c: &mut Criterion) {
     .unwrap();
     let eos = Arc::new(PcSaft::new(Arc::new(parameters)));
     let t = 300.0 * KELVIN;
-    let density = 71.18 * KILO * MOL / METER.powi(3);
+    let density = 71.18 * KILO * MOL / METER.powi::<P3>();
     let v = 100.0 * MOL / density;
     let x = arr1(&[1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0]);
     let m = &x * 100.0 * MOL;

feos-core/CHANGELOG.md

@@ -12,8 +12,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added `JobackParameters` struct that implements `Parameters` including Python bindings. [#158](https://github.com/feos-org/feos/pull/158)
 - Added `Parameter::from_model_records` as a simpler interface to generate parameters.  [#169](https://github.com/feos-org/feos/pull/169)
 - Added optional `Phase` argument for constructors of dynamic properties of `DataSet`s.  [#174](https://github.com/feos-org/feos/pull/174)
-- Added `molar_weight` method to `Residual` trait. [#177](https://github.com/feos-org/feos/pull/158)
-- Added molar versions for entropy, enthalpy, etc. for residual properties. [#177](https://github.com/feos-org/feos/pull/158)
+- Added `molar_weight` method to `Residual` trait. [#177](https://github.com/feos-org/feos/pull/177)
+- Added molar versions for entropy, enthalpy, etc. for residual properties. [#177](https://github.com/feos-org/feos/pull/177)
 
 ### Changed
 - Changed constructors of `Parameter` trait to return `Result`s. [#161](https://github.com/feos-org/feos/pull/161)
@@ -28,16 +28,17 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Bubble and dew point iterations will not attempt a second iteration if no solution is found for the given initial pressure. [#166](https://github.com/feos-org/feos/pull/166)
 - Made the binary records in the constructions and getters of the `Parameter` trait optional. [#169](https://github.com/feos-org/feos/pull/169)
 - Changed the second argument of `new_binary` in Python from a `BinaryRecord` to the corresponding binary model record (analogous to the Rust implementation). [#169](https://github.com/feos-org/feos/pull/169)
-- Renamed `c_v` and `c_p` to `isochoric_heat_capacity` and `isobaric_heat_capacity`, respectively, and added prefixes for molar and specific properties. [#177](https://github.com/feos-org/feos/pull/158)
-- `State.helmholtz_energy_contributions` in Python now accepts optional `Contributions`. [#177](https://github.com/feos-org/feos/pull/158)
-- Changed `StateVec` Python getters for entropy and enthalpy to functions that accept optional `Contributions`. [#177](https://github.com/feos-org/feos/pull/158)
-- `PhaseDiagram.to_dict` in Python now accepts optional `Contributions` and includes mass specific properties. [#177](https://github.com/feos-org/feos/pull/158)
+- Renamed `c_v` and `c_p` to `isochoric_heat_capacity` and `isobaric_heat_capacity`, respectively, and added prefixes for molar and specific properties. [#177](https://github.com/feos-org/feos/pull/177)
+- `State.helmholtz_energy_contributions` in Python now accepts optional `Contributions`. [#177](https://github.com/feos-org/feos/pull/177)
+- Changed `StateVec` Python getters for entropy and enthalpy to functions that accept optional `Contributions`. [#177](https://github.com/feos-org/feos/pull/177)
+- `PhaseDiagram.to_dict` in Python now accepts optional `Contributions` and includes mass specific properties. [#177](https://github.com/feos-org/feos/pull/177)
+- Replaced the run-time unit checks from the `quantity` crate with compile-time unit checks with custom implementations in the new `feos_core.si` module. [#181](https://github.com/feos-org/feos/pull/181)
 
 ### Removed
 - Removed `EquationOfState` trait. [#158](https://github.com/feos-org/feos/pull/158)
 - Removed ideal gas dependencies from `PureRecord` and `SegmentRecord`. [#158](https://github.com/feos-org/feos/pull/158)
 - Removed Python getter and setter functions and optional arguments for ideal gas records in macros. [#158](https://github.com/feos-org/feos/pull/158)
-- Removed `MolarWeight` trait. [#177](https://github.com/feos-org/feos/pull/158)
+- Removed `MolarWeight` trait. [#177](https://github.com/feos-org/feos/pull/177)
 
 ### Fixed
 - The vapor and liquid states in a bubble or dew point iteration are assigned correctly according to the inputs, rather than based on the mole density which can be incorrect for mixtures with large differences in molar weights. [#166](https://github.com/feos-org/feos/pull/166)

feos-core/Cargo.toml

@@ -18,7 +18,7 @@ rustdoc-args = [ "--html-in-header", "./docs-header.html" ]
 features = [ "rayon" ]
 
 [dependencies]
-quantity = "0.6"
+quantity = { version = "0.6", optional = true }
 num-dual = { version = "0.7", features = ["linalg"] }
 ndarray = { version = "0.15", features = ["serde"] }
 nalgebra = "0.32"
@@ -31,9 +31,10 @@ conv = "0.3"
 numpy = { version = "0.18", optional = true }
 pyo3 = { version = "0.18", optional = true }
 rayon = { version = "1.5", optional = true }
-
-[dev-dependencies]
+typenum = "1.16"
 approx = "0.4"
+regex = "1.9"
+once_cell = "1.18"
 
 [features]
 default = []

feos-core/src/cubic.rs

@@ -6,11 +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, HelmholtzEnergy, HelmholtzEnergyDual, Residual};
 use crate::parameter::{Identifier, Parameter, ParameterError, PureRecord};
-use crate::si::{GRAM, MOL};
+use crate::si::{MolarWeight, GRAM, MOL};
 use crate::state::StateHD;
 use ndarray::{Array1, Array2};
 use num_dual::DualNum;
-use quantity::si::SIArray1;
 use serde::{Deserialize, Serialize};
 use std::f64::consts::SQRT_2;
 use std::fmt;
@@ -232,18 +231,18 @@ impl Residual for PengRobinson {
         &self.contributions
     }
 
-    fn molar_weight(&self) -> SIArray1 {
-        self.parameters.molarweight.clone() * GRAM / MOL
+    fn molar_weight(&self) -> MolarWeight<Array1<f64>> {
+        &self.parameters.molarweight * (GRAM / MOL)
     }
 }
 
 #[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::si::*;
     use std::sync::Arc;
 
     fn pure_record_vec() -> Vec<PureRecord<PengRobinsonRecord>> {