update num-dual dependency

feos-org · prehner · May 31, 2023 · May 30, 2023 · May 30, 2023 · May 31, 2023
commit 8c733eea79e3f9bc1573cec893728a3d4fb691a4
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
@@ -2,9 +2,9 @@ name: Test
 
 on:
   push:
-    branches: [main]
+    branches: [main, development]
   pull_request:
-    branches: [main]
+    branches: [main, development]
 
 env:
   CARGO_TERM_COLOR: always

diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml
@@ -1,9 +1,9 @@
 name: Build Wheels
 on:
   push:
-    branches: [main]
+    branches: [main, development]
   pull_request:
-    branches: [main]
+    branches: [main, development]
 jobs:
   linux:
     runs-on: ubuntu-latest

diff --git a/Cargo.toml b/Cargo.toml
@@ -23,7 +23,7 @@ crate-type = ["rlib", "cdylib"]
 
 [dependencies]
 quantity = "0.6"
-num-dual = "0.6"
+num-dual = "0.7"
 feos-core = { version = "0.4", path = "feos-core" }
 feos-dft = { version = "0.4", path = "feos-dft", optional = true }
 feos-derive = { version = "0.2", path = "feos-derive" }

diff --git a/benches/dual_numbers.rs b/benches/dual_numbers.rs
@@ -28,7 +28,7 @@ fn state_pcsaft(parameters: PcSaftParameters) -> State<PcSaft> {
 }
 
 /// Residual Helmholtz energy given an equation of state and a StateHD.
-fn a_res<D: DualNum<f64>, E: EquationOfState>(inp: (&Arc<E>, &StateHD<D>)) -> D
+fn a_res<D: DualNum<f64> + Copy, E: EquationOfState>(inp: (&Arc<E>, &StateHD<D>)) -> D
 where
     (dyn HelmholtzEnergy + 'static): HelmholtzEnergyDual<D>,
 {

diff --git a/feos-core/Cargo.toml b/feos-core/Cargo.toml
@@ -19,8 +19,9 @@ features = [ "rayon" ]
 
 [dependencies]
 quantity = "0.6"
-num-dual = { version = "0.6", features = ["linalg"] }
+num-dual = { version = "0.7", features = ["linalg"] }
 ndarray = { version = "0.15", features = ["serde"] }
+nalgebra = "0.32"
 num-traits = "0.2"
 thiserror = "1.0"
 serde = { version = "1.0", features = ["derive"] }
@@ -37,4 +38,4 @@ approx = "0.4"
 [features]
 default = []
 rayon = ["dep:rayon", "ndarray/rayon"]
-python = ["pyo3", "numpy", "quantity/python", "num-dual/python", "rayon"]
+python = ["pyo3", "numpy", "quantity/python", "num-dual/python_macro", "rayon"]
diff --git a/feos-core/src/cubic.rs b/feos-core/src/cubic.rs
@@ -170,7 +170,7 @@ struct PengRobinsonContribution {
     parameters: Arc<PengRobinsonParameters>,
 }
 
-impl<D: DualNum<f64>> HelmholtzEnergyDual<D> for PengRobinsonContribution {
+impl<D: DualNum<f64> + Copy> HelmholtzEnergyDual<D> for PengRobinsonContribution {
     fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
         // temperature dependent a parameter
         let p = &self.parameters;

diff --git a/feos-core/src/equation_of_state.rs b/feos-core/src/equation_of_state.rs
@@ -3,9 +3,10 @@ use crate::state::StateHD;
 use crate::EosUnit;
 use ndarray::prelude::*;
 use num_dual::{
-    Dual, Dual2_64, Dual3, Dual3_64, Dual64, DualNum, DualVec64, HyperDual, HyperDual64,
+    first_derivative, second_derivative, third_derivative, Dual, Dual2, Dual2_64, Dual3, Dual3_64,
+    Dual64, DualNum, DualSVec64, HyperDual, HyperDual64,
 };
-use num_traits::{One, Zero};
+use num_traits::Zero;
 use quantity::si::{SIArray1, SINumber, SIUnit};
 use std::fmt;
 
@@ -28,16 +29,17 @@ pub trait HelmholtzEnergyDual<D: DualNum<f64>> {
 pub trait HelmholtzEnergy:
     HelmholtzEnergyDual<f64>
     + HelmholtzEnergyDual<Dual64>
-    + HelmholtzEnergyDual<Dual<DualVec64<3>, f64>>
+    + HelmholtzEnergyDual<Dual<DualSVec64<3>, f64>>
     + HelmholtzEnergyDual<HyperDual64>
     + HelmholtzEnergyDual<Dual2_64>
     + HelmholtzEnergyDual<Dual3_64>
     + HelmholtzEnergyDual<HyperDual<Dual64, f64>>
-    + HelmholtzEnergyDual<HyperDual<DualVec64<2>, f64>>
-    + HelmholtzEnergyDual<HyperDual<DualVec64<3>, f64>>
+    + HelmholtzEnergyDual<HyperDual<DualSVec64<2>, f64>>
+    + HelmholtzEnergyDual<HyperDual<DualSVec64<3>, f64>>
+    + HelmholtzEnergyDual<Dual2<Dual64, f64>>
     + HelmholtzEnergyDual<Dual3<Dual64, f64>>
-    + HelmholtzEnergyDual<Dual3<DualVec64<2>, f64>>
-    + HelmholtzEnergyDual<Dual3<DualVec64<3>, f64>>
+    + HelmholtzEnergyDual<Dual3<DualSVec64<2>, f64>>
+    + HelmholtzEnergyDual<Dual3<DualSVec64<3>, f64>>
     + fmt::Display
     + Send
     + Sync
@@ -47,16 +49,17 @@ pub trait HelmholtzEnergy:
 impl<T> HelmholtzEnergy for T where
     T: HelmholtzEnergyDual<f64>
         + HelmholtzEnergyDual<Dual64>
-        + HelmholtzEnergyDual<Dual<DualVec64<3>, f64>>
+        + HelmholtzEnergyDual<Dual<DualSVec64<3>, f64>>
         + HelmholtzEnergyDual<HyperDual64>
         + HelmholtzEnergyDual<Dual2_64>
         + HelmholtzEnergyDual<Dual3_64>
         + HelmholtzEnergyDual<HyperDual<Dual64, f64>>
-        + HelmholtzEnergyDual<HyperDual<DualVec64<2>, f64>>
-        + HelmholtzEnergyDual<HyperDual<DualVec64<3>, f64>>
+        + HelmholtzEnergyDual<HyperDual<DualSVec64<2>, f64>>
+        + HelmholtzEnergyDual<HyperDual<DualSVec64<3>, f64>>
+        + HelmholtzEnergyDual<Dual2<Dual64, f64>>
         + HelmholtzEnergyDual<Dual3<Dual64, f64>>
-        + HelmholtzEnergyDual<Dual3<DualVec64<2>, f64>>
-        + HelmholtzEnergyDual<Dual3<DualVec64<3>, f64>>
+        + HelmholtzEnergyDual<Dual3<DualSVec64<2>, f64>>
+        + HelmholtzEnergyDual<Dual3<DualSVec64<3>, f64>>
         + fmt::Display
         + Send
         + Sync
@@ -70,7 +73,7 @@ impl<T> HelmholtzEnergy for T where
 /// the specific types in the supertraits of [IdealGasContribution]
 /// so that the implementor can be used as an ideal gas
 /// contribution in the equation of state.
-pub trait IdealGasContributionDual<D: DualNum<f64>> {
+pub trait IdealGasContributionDual<D: DualNum<f64> + Copy> {
     /// The thermal de Broglie wavelength of each component in the form $\ln\left(\frac{\Lambda^3}{\AA^3}\right)$
     fn de_broglie_wavelength(&self, temperature: D, components: usize) -> Array1<D>;
 
@@ -101,39 +104,41 @@ pub trait IdealGasContributionDual<D: DualNum<f64>> {
 pub trait IdealGasContribution:
     IdealGasContributionDual<f64>
     + IdealGasContributionDual<Dual64>
-    + IdealGasContributionDual<Dual<DualVec64<3>, f64>>
+    + IdealGasContributionDual<Dual<DualSVec64<3>, f64>>
     + IdealGasContributionDual<HyperDual64>
     + IdealGasContributionDual<Dual2_64>
     + IdealGasContributionDual<Dual3_64>
     + IdealGasContributionDual<HyperDual<Dual64, f64>>
-    + IdealGasContributionDual<HyperDual<DualVec64<2>, f64>>
-    + IdealGasContributionDual<HyperDual<DualVec64<3>, f64>>
+    + IdealGasContributionDual<HyperDual<DualSVec64<2>, f64>>
+    + IdealGasContributionDual<HyperDual<DualSVec64<3>, f64>>
+    + IdealGasContributionDual<Dual2<Dual64, f64>>
     + IdealGasContributionDual<Dual3<Dual64, f64>>
-    + IdealGasContributionDual<Dual3<DualVec64<2>, f64>>
-    + IdealGasContributionDual<Dual3<DualVec64<3>, f64>>
+    + IdealGasContributionDual<Dual3<DualSVec64<2>, f64>>
+    + IdealGasContributionDual<Dual3<DualSVec64<3>, f64>>
     + fmt::Display
 {
 }
 
 impl<T> IdealGasContribution for T where
     T: IdealGasContributionDual<f64>
         + IdealGasContributionDual<Dual64>
-        + IdealGasContributionDual<Dual<DualVec64<3>, f64>>
+        + IdealGasContributionDual<Dual<DualSVec64<3>, f64>>
         + IdealGasContributionDual<HyperDual64>
         + IdealGasContributionDual<Dual2_64>
         + IdealGasContributionDual<Dual3_64>
         + IdealGasContributionDual<HyperDual<Dual64, f64>>
-        + IdealGasContributionDual<HyperDual<DualVec64<2>, f64>>
-        + IdealGasContributionDual<HyperDual<DualVec64<3>, f64>>
+        + IdealGasContributionDual<HyperDual<DualSVec64<2>, f64>>
+        + IdealGasContributionDual<HyperDual<DualSVec64<3>, f64>>
+        + IdealGasContributionDual<Dual2<Dual64, f64>>
         + IdealGasContributionDual<Dual3<Dual64, f64>>
-        + IdealGasContributionDual<Dual3<DualVec64<2>, f64>>
-        + IdealGasContributionDual<Dual3<DualVec64<3>, f64>>
+        + IdealGasContributionDual<Dual3<DualSVec64<2>, f64>>
+        + IdealGasContributionDual<Dual3<DualSVec64<3>, f64>>
         + fmt::Display
 {
 }
 
 struct DefaultIdealGasContribution;
-impl<D: DualNum<f64>> IdealGasContributionDual<D> for DefaultIdealGasContribution {
+impl<D: DualNum<f64> + Copy> IdealGasContributionDual<D> for DefaultIdealGasContribution {
     fn de_broglie_wavelength(&self, _: D, components: usize) -> Array1<D> {
         Array1::zeros(components)
     }
@@ -175,7 +180,7 @@ pub trait EquationOfState: Send + Sync {
     fn residual(&self) -> &[Box<dyn HelmholtzEnergy>];
 
     /// Evaluate the residual reduced Helmholtz energy $\beta A^\mathrm{res}$.
-    fn evaluate_residual<D: DualNum<f64>>(&self, state: &StateHD<D>) -> D
+    fn evaluate_residual<D: DualNum<f64> + Copy>(&self, state: &StateHD<D>) -> D
     where
         dyn HelmholtzEnergy: HelmholtzEnergyDual<D>,
     {
@@ -187,7 +192,7 @@ pub trait EquationOfState: Send + Sync {
 
     /// Evaluate the reduced Helmholtz energy of each individual contribution
     /// and return them together with a string representation of the contribution.
-    fn evaluate_residual_contributions<D: DualNum<f64>>(
+    fn evaluate_residual_contributions<D: DualNum<f64> + Copy>(
         &self,
         state: &StateHD<D>,
     ) -> Vec<(String, D)>
@@ -252,12 +257,10 @@ pub trait EquationOfState: Send + Sync {
     ) -> 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(SIUnit::reference_temperature())?);
-        let s = StateHD::new_virial(t, rho, x);
-        Ok(self.evaluate_residual(&s).eps1eps2[(0, 0)] * 0.5 / SIUnit::reference_density())
+        let t = temperature.to_reduced(SIUnit::reference_temperature())?;
+        let a_res = |rho| self.evaluate_residual(&StateHD::new_virial(t.into(), rho, x));
+        let (_, _, b) = second_derivative(a_res, 0.0);
+        Ok(b * 0.5 / SIUnit::reference_density())
     }
 
     /// Calculate the third virial coefficient $C(T)$
@@ -268,10 +271,10 @@ pub trait EquationOfState: Send + Sync {
     ) -> 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(SIUnit::reference_temperature())?);
-        let s = StateHD::new_virial(t, rho, x);
-        Ok(self.evaluate_residual(&s).v3 / 3.0 / SIUnit::reference_density().powi(2))
+        let t = temperature.to_reduced(SIUnit::reference_temperature())?;
+        let a_res = |rho| self.evaluate_residual(&StateHD::new_virial(t.into(), rho, x));
+        let (_, _, _, c) = third_derivative(a_res, 0.0);
+        Ok(c / 3.0 / SIUnit::reference_density().powi(2))
     }
 
     /// Calculate the temperature derivative of the second virial coefficient $B'(T)$
@@ -282,15 +285,16 @@ pub trait EquationOfState: Send + Sync {
     ) -> 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(SIUnit::reference_temperature())?).derive(),
-        );
-        let s = StateHD::new_virial(t, rho, x);
-        Ok(self.evaluate_residual(&s).eps1eps2[(0, 0)].eps[0] * 0.5
-            / (SIUnit::reference_density() * SIUnit::reference_temperature()))
+        let t = temperature.to_reduced(SIUnit::reference_temperature())?;
+        let b = |t| {
+            let a_res = |rho: Dual2<Dual64, f64>| {
+                self.evaluate_residual(&StateHD::new_virial(Dual2::from_re(t), rho, x))
+            };
+            let (_, _, b) = second_derivative(a_res, Dual64::zero());
+            b
+        };
+        let (_, b_t) = first_derivative(b, t);
+        Ok(b_t * 0.5 / (SIUnit::reference_density() * SIUnit::reference_temperature()))
     }
 
     /// Calculate the temperature derivative of the third virial coefficient $C'(T)$
@@ -301,14 +305,15 @@ pub trait EquationOfState: Send + Sync {
     ) -> 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(SIUnit::reference_temperature())?).derive(),
-        );
-        let s = StateHD::new_virial(t, rho, x);
-        Ok(self.evaluate_residual(&s).v3.eps[0]
-            / 3.0
-            / (SIUnit::reference_density().powi(2) * SIUnit::reference_temperature()))
+        let t = temperature.to_reduced(SIUnit::reference_temperature())?;
+        let c = |t| {
+            let a_res =
+                |rho| self.evaluate_residual(&StateHD::new_virial(Dual3::from_re(t), rho, x));
+            let (_, _, _, c) = third_derivative(a_res, Dual64::zero());
+            c
+        };
+        let (_, c_t) = first_derivative(c, t);
+        Ok(c_t / 3.0 / (SIUnit::reference_density().powi(2) * SIUnit::reference_temperature()))
     }
 }
 

diff --git a/feos-core/src/joback.rs b/feos-core/src/joback.rs
@@ -107,7 +107,7 @@ const P0: f64 = 1.0e5;
 const A3: f64 = 1e-30;
 const KB: f64 = 1.38064852e-23;
 
-impl<D: DualNum<f64>> IdealGasContributionDual<D> for Joback {
+impl<D: DualNum<f64> + Copy> IdealGasContributionDual<D> for Joback {
     fn de_broglie_wavelength(&self, temperature: D, components: usize) -> Array1<D> {
         let t = temperature;
         let t2 = t * t;

diff --git a/feos-core/src/python/user_defined.rs b/feos-core/src/python/user_defined.rs
@@ -2,11 +2,11 @@ use crate::{EquationOfState, HelmholtzEnergy, HelmholtzEnergyDual, MolarWeight,
 use ndarray::Array1;
 use num_dual::*;
 use numpy::convert::IntoPyArray;
-use numpy::{PyReadonlyArrayDyn, PyArray};
+use numpy::{PyArray, PyReadonlyArrayDyn};
 use pyo3::exceptions::PyTypeError;
 use pyo3::prelude::*;
 use quantity::python::PySIArray1;
-use quantity::si::{SIArray1};
+use quantity::si::SIArray1;
 use std::fmt;
 
 struct PyHelmholtzEnergy(Py<PyAny>);
@@ -205,35 +205,36 @@ state!(PyStateF, f64, f64);
 helmholtz_energy!(PyStateF, f64, f64);
 
 impl_dual_state_helmholtz_energy!(PyStateD, PyDual64, Dual64, f64);
-dual_number!(PyDualVec3, DualVec64<3>, f64);
+dual_number!(PyDualVec3, DualSVec64<3>, f64);
 impl_dual_state_helmholtz_energy!(
     PyStateDualDualVec3,
     PyDualDualVec3,
-    Dual<DualVec64<3>, f64>,
+    Dual<DualSVec64<3>, f64>,
     PyDualVec3
 );
-impl_dual_state_helmholtz_energy!(
-    PyStateHD,
-    PyHyperDual64,
-    HyperDual64,
-    f64
-);
+impl_dual_state_helmholtz_energy!(PyStateHD, PyHyperDual64, HyperDual64, f64);
 impl_dual_state_helmholtz_energy!(PyStateD2, PyDual2_64, Dual2_64, f64);
 impl_dual_state_helmholtz_energy!(PyStateD3, PyDual3_64, Dual3_64, f64);
 impl_dual_state_helmholtz_energy!(PyStateHDD, PyHyperDualDual64, HyperDual<Dual64, f64>, PyDual64);
-dual_number!(PyDualVec2, DualVec64<2>, f64);
+dual_number!(PyDualVec2, DualSVec64<2>, f64);
 impl_dual_state_helmholtz_energy!(
     PyStateHDDVec2,
     PyHyperDualVec2,
-    HyperDual<DualVec64<2>, f64>,
+    HyperDual<DualSVec64<2>, f64>,
     PyDualVec2
 );
 impl_dual_state_helmholtz_energy!(
     PyStateHDDVec3,
     PyHyperDualVec3,
-    HyperDual<DualVec64<3>, f64>,
+    HyperDual<DualSVec64<3>, f64>,
     PyDualVec3
 );
+impl_dual_state_helmholtz_energy!(
+    PyStateD2D,
+    PyDual2Dual64,
+    Dual2<Dual64, f64>,
+    PyDual64
+);
 impl_dual_state_helmholtz_energy!(
     PyStateD3D,
     PyDual3Dual64,
@@ -243,12 +244,12 @@ impl_dual_state_helmholtz_energy!(
 impl_dual_state_helmholtz_energy!(
     PyStateD3DVec2,
     PyDual3DualVec2,
-    Dual3<DualVec64<2>, f64>,
+    Dual3<DualSVec64<2>, f64>,
     PyDualVec2
 );
 impl_dual_state_helmholtz_energy!(
     PyStateD3DVec3,
     PyDual3DualVec3,
-    Dual3<DualVec64<3>, f64>,
+    Dual3<DualSVec64<3>, f64>,
     PyDualVec3
 );