use derivatives with units for state properties

feos-org · g-bauer · Oct 13, 2025 · Aug 18, 2025 · Aug 18, 2025 · Aug 18, 2025
commit ba7d6ec80cd624e0484f703ab50a415e1563ea79
diff --git a/crates/feos-core/src/equation_of_state/mod.rs b/crates/feos-core/src/equation_of_state/mod.rs
@@ -1,11 +1,11 @@
 use std::ops::Deref;
 
-use crate::state::StateHD;
+use crate::{ReferenceSystem, state::StateHD};
 use nalgebra::{
     Const, DVector, DefaultAllocator, Dim, Dyn, OVector, SVector, U1, allocator::Allocator,
 };
 use num_dual::DualNum;
-use quantity::MolarWeight;
+use quantity::{Energy, MolarEnergy, MolarWeight, Moles, Temperature, Volume};
 
 mod ideal_gas;
 mod residual;
@@ -129,13 +129,13 @@ where
         )
     }
 
-    fn ideal_gas_helmholtz_energy<D2: DualNum<f64, Inner = D> + Copy>(
+    fn ideal_gas_molar_helmholtz_energy<D2: DualNum<f64, Inner = D> + Copy>(
         &self,
         temperature: D2,
-        volume: D2,
-        moles: &OVector<D2, N>,
+        molar_volume: D2,
+        molefracs: &OVector<D2, N>,
     ) -> D2 {
-        let partial_density = moles / volume;
+        let partial_density = molefracs / molar_volume;
         let mut res = D2::from(0.0);
         for (i, &r) in self.ideal_gas().zip(partial_density.iter()) {
             let ln_rho_m1 = if r.re() == 0.0 {
@@ -145,11 +145,32 @@ where
             };
             res += r * (i.lift().ln_lambda3(temperature) + ln_rho_m1)
         }
-        res * volume * temperature
+        res * molar_volume * temperature
     }
 
-    fn ideal_gas_helmholtz_energy_0(&self, temperature: D, volume: D, moles: &OVector<D, N>) -> D {
-        let partial_density = moles / volume;
+    fn ideal_gas_helmholtz_energy_unit<D2: DualNum<f64, Inner = D> + Copy>(
+        &self,
+        temperature: Temperature<D2>,
+        volume: Volume<D2>,
+        moles: &Moles<OVector<D2, N>>,
+    ) -> Energy<D2> {
+        let total_moles = moles.sum();
+        let molefracs = moles / total_moles;
+        let molar_volume = volume / total_moles;
+        MolarEnergy::from_reduced(self.ideal_gas_molar_helmholtz_energy(
+            temperature.into_reduced(),
+            molar_volume.into_reduced(),
+            &molefracs,
+        )) * total_moles
+    }
+
+    fn ideal_gas_molar_helmholtz_energy_0(
+        &self,
+        temperature: D,
+        molar_volume: D,
+        molefracs: &OVector<D, N>,
+    ) -> D {
+        let partial_density = molefracs / molar_volume;
         let mut res = D::from(0.0);
         for (i, &r) in self.ideal_gas().zip(partial_density.iter()) {
             let ln_rho_m1 = if r.re() == 0.0 {
@@ -159,7 +180,23 @@ where
             };
             res += r * (i.ln_lambda3(temperature) + ln_rho_m1)
         }
-        res * volume * temperature
+        res * molar_volume * temperature
+    }
+
+    fn ideal_gas_helmholtz_energy_0(
+        &self,
+        temperature: Temperature<D>,
+        volume: Volume<D>,
+        moles: &Moles<OVector<D, N>>,
+    ) -> Energy<D> {
+        let total_moles = moles.sum();
+        let molefracs = moles / total_moles;
+        let molar_volume = volume / total_moles;
+        MolarEnergy::from_reduced(self.ideal_gas_molar_helmholtz_energy_0(
+            temperature.into_reduced(),
+            molar_volume.into_reduced(),
+            &molefracs,
+        )) * total_moles
     }
 }
 

diff --git a/crates/feos-core/src/equation_of_state/residual.rs b/crates/feos-core/src/equation_of_state/residual.rs
@@ -198,6 +198,22 @@ where
         self.reduced_residual_helmholtz_energy_density(&state) * temperature * volume
     }
 
+    /// Evaluate the residual Helmholtz energy $A^\mathrm{res}$.
+    fn residual_helmholtz_energy_unit(
+        &self,
+        temperature: Temperature<D>,
+        volume: Volume<D>,
+        moles: &Moles<OVector<D, N>>,
+    ) -> Energy<D> {
+        let temperature = temperature.into_reduced();
+        let total_moles = moles.sum();
+        let molar_volume = (volume / total_moles).into_reduced();
+        let molefracs = moles / total_moles;
+        let state = StateHD::new(temperature, molar_volume, &molefracs);
+        Pressure::from_reduced(self.reduced_residual_helmholtz_energy_density(&state) * temperature)
+            * volume
+    }
+
     /// Check if the provided optional mole number is consistent with the
     /// equation of state.
     ///

diff --git a/crates/feos-core/src/state/cache.rs b/crates/feos-core/src/state/cache.rs
@@ -1,19 +1,26 @@
-use crate::state::VectorPartialDerivative;
-
-use super::{Derivative, PartialDerivative};
 use nalgebra::allocator::Allocator;
-use nalgebra::{DefaultAllocator, Dim, OMatrix, OVector, Scalar};
-use std::collections::HashMap;
+use nalgebra::{DefaultAllocator, Dim, OVector, Scalar};
+use quantity::*;
+use std::sync::OnceLock;
+use typenum::Diff;
 
 #[derive(Clone, Debug)]
+#[expect(clippy::type_complexity)]
 pub struct Cache<D: Scalar, N: Dim>
 where
     DefaultAllocator: Allocator<N>,
 {
-    pub map: HashMap<PartialDerivative, D>,
-    pub vec_map: HashMap<VectorPartialDerivative, OVector<D, N>>,
-    pub hit: u64,
-    pub miss: u64,
+    pub a: OnceLock<Energy<D>>,
+    pub da_dt: OnceLock<Entropy<D>>,
+    pub da_dv: OnceLock<Pressure<D>>,
+    pub da_dn: OnceLock<MolarEnergy<OVector<D, N>>>,
+    pub d2a_dt2: OnceLock<Quantity<D, Diff<_Entropy, _Temperature>>>,
+    pub d2a_dv2: OnceLock<Quantity<D, Diff<_Pressure, _Volume>>>,
+    pub d2a_dtdv: OnceLock<Quantity<D, Diff<_Pressure, _Temperature>>>,
+    pub d2a_dndt: OnceLock<Quantity<OVector<D, N>, Diff<_MolarEnergy, _Temperature>>>,
+    pub d2a_dndv: OnceLock<Quantity<OVector<D, N>, Diff<_MolarEnergy, _Volume>>>,
+    pub d3a_dt3: OnceLock<Quantity<D, Diff<Diff<_Entropy, _Temperature>, _Temperature>>>,
+    pub d3a_dv3: OnceLock<Quantity<D, Diff<Diff<_Pressure, _Volume>, _Volume>>>,
 }
 
 impl<D: Scalar + Copy, N: Dim> Cache<D, N>
@@ -22,154 +29,17 @@ where
 {
     pub fn new() -> Self {
         Self {
-            map: HashMap::with_capacity(8),
-            vec_map: HashMap::with_capacity(3),
-            hit: 0,
-            miss: 0,
-        }
-    }
-
-    pub fn get_or_insert_with_zeroth<F: FnOnce() -> D>(&mut self, f: F) -> D {
-        if let Some(&value) = self.map.get(&PartialDerivative::Zeroth) {
-            self.hit += 1;
-            value
-        } else {
-            self.miss += 1;
-            let value = f();
-            self.map.insert(PartialDerivative::Zeroth, value);
-            value
-        }
-    }
-
-    pub fn get_or_insert_with_first_scalar<F: FnOnce() -> (D, D)>(
-        &mut self,
-        derivative: Derivative,
-        f: F,
-    ) -> D {
-        if let Some(&value) = self.map.get(&PartialDerivative::First(derivative)) {
-            self.hit += 1;
-            value
-        } else {
-            self.miss += 1;
-            let (f, df) = f();
-            self.map.insert(PartialDerivative::Zeroth, f);
-            self.map.insert(PartialDerivative::First(derivative), df);
-            df
-        }
-    }
-
-    pub fn get_or_insert_with_second<F: FnOnce() -> (D, D, D)>(
-        &mut self,
-        derivative: Derivative,
-        f: F,
-    ) -> D {
-        if let Some(&value) = self.map.get(&PartialDerivative::Second(derivative)) {
-            self.hit += 1;
-            value
-        } else {
-            self.miss += 1;
-            let (f, df, d2f) = f();
-            self.map.insert(PartialDerivative::Zeroth, f);
-            self.map.insert(PartialDerivative::First(derivative), df);
-            self.map.insert(PartialDerivative::Second(derivative), d2f);
-            d2f
-        }
-    }
-
-    pub fn get_or_insert_with_second_mixed_scalar<F: FnOnce() -> (D, D, D, D)>(
-        &mut self,
-        derivative1: Derivative,
-        derivative2: Derivative,
-        f: F,
-    ) -> D {
-        if let Some(&value) = self.map.get(&PartialDerivative::SecondMixed) {
-            self.hit += 1;
-            value
-        } else {
-            self.miss += 1;
-            let (f, df1, df2, d2f) = f();
-            self.map.insert(PartialDerivative::Zeroth, f);
-            self.map.insert(PartialDerivative::First(derivative1), df1);
-            self.map.insert(PartialDerivative::First(derivative2), df2);
-            self.map.insert(PartialDerivative::SecondMixed, d2f);
-            d2f
-        }
-    }
-
-    pub fn get_or_insert_with_third<F: FnOnce() -> (D, D, D, D)>(
-        &mut self,
-        derivative: Derivative,
-        f: F,
-    ) -> D {
-        if let Some(&value) = self.map.get(&PartialDerivative::Third(derivative)) {
-            self.hit += 1;
-            value
-        } else {
-            self.miss += 1;
-            let (f, df, d2f, d3f) = f();
-            self.map.insert(PartialDerivative::Zeroth, f);
-            self.map.insert(PartialDerivative::First(derivative), df);
-            self.map.insert(PartialDerivative::Second(derivative), d2f);
-            self.map.insert(PartialDerivative::Third(derivative), d3f);
-            d3f
-        }
-    }
-
-    pub fn get_or_insert_with_first_vector<F: FnOnce() -> (D, OVector<D, N>)>(
-        &mut self,
-        f: F,
-    ) -> OVector<D, N> {
-        if let Some(value) = self.vec_map.get(&VectorPartialDerivative::First) {
-            self.hit += 1;
-            value.clone()
-        } else {
-            self.miss += 1;
-            let (f, df) = f();
-            self.map.insert(PartialDerivative::Zeroth, f);
-            self.vec_map
-                .insert(VectorPartialDerivative::First, df.clone());
-            df
-        }
-    }
-
-    pub fn get_or_insert_with_second_vector<F: FnOnce() -> (D, OVector<D, N>, OMatrix<D, N, N>)>(
-        &mut self,
-        f: F,
-    ) -> OMatrix<D, N, N>
-    where
-        DefaultAllocator: Allocator<N, N>,
-    {
-        self.miss += 1;
-        let (f, df, d2f) = f();
-        self.map.insert(PartialDerivative::Zeroth, f);
-        self.vec_map.insert(VectorPartialDerivative::First, df);
-        d2f
-    }
-
-    pub fn get_or_insert_with_second_mixed_vector<
-        F: FnOnce() -> (D, OVector<D, N>, D, OVector<D, N>),
-    >(
-        &mut self,
-        derivative: Derivative,
-        f: F,
-    ) -> OVector<D, N> {
-        if let Some(value) = self
-            .vec_map
-            .get(&VectorPartialDerivative::SecondMixed(derivative))
-        {
-            self.hit += 1;
-            value.clone()
-        } else {
-            self.miss += 1;
-            let (f, df1, df2, d2f) = f();
-            self.map.insert(PartialDerivative::Zeroth, f);
-            self.vec_map.insert(VectorPartialDerivative::First, df1);
-            self.map.insert(PartialDerivative::First(derivative), df2);
-            self.vec_map.insert(
-                VectorPartialDerivative::SecondMixed(derivative),
-                d2f.clone(),
-            );
-            d2f
+            a: OnceLock::new(),
+            da_dt: OnceLock::new(),
+            da_dv: OnceLock::new(),
+            da_dn: OnceLock::new(),
+            d2a_dt2: OnceLock::new(),
+            d2a_dv2: OnceLock::new(),
+            d2a_dtdv: OnceLock::new(),
+            d2a_dndt: OnceLock::new(),
+            d2a_dndv: OnceLock::new(),
+            d3a_dt3: OnceLock::new(),
+            d3a_dv3: OnceLock::new(),
         }
     }
 }