From 8c733eea79e3f9bc1573cec893728a3d4fb691a4 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Tue, 30 May 2023 12:51:10 +0200 Subject: [PATCH] update num-dual dependency --- .github/workflows/test.yml | 4 +- .github/workflows/wheels.yml | 4 +- Cargo.toml | 2 +- benches/dual_numbers.rs | 2 +- feos-core/Cargo.toml | 5 +- feos-core/src/cubic.rs | 2 +- feos-core/src/equation_of_state.rs | 111 ++++++----- feos-core/src/joback.rs | 2 +- feos-core/src/python/user_defined.rs | 31 +-- feos-core/src/state/cache.rs | 18 +- feos-core/src/state/critical_point.rs | 186 +++++++++--------- feos-core/src/state/mod.rs | 32 +-- feos-core/src/state/properties.rs | 36 ++-- feos-dft/Cargo.toml | 3 +- feos-dft/src/convolver/mod.rs | 4 +- feos-dft/src/functional.rs | 22 +-- feos-dft/src/functional_contribution.rs | 59 +++--- feos-dft/src/ideal_chain_contribution.rs | 2 +- feos-dft/src/pdgt.rs | 13 +- feos-dft/src/profile.rs | 10 +- feos-dft/src/solver.rs | 9 +- feos-dft/src/weight_functions.rs | 4 +- src/association/dft.rs | 14 +- src/association/mod.rs | 41 ++-- src/gc_pcsaft/dft/dispersion.rs | 4 +- src/gc_pcsaft/dft/hard_chain.rs | 2 +- src/gc_pcsaft/dft/mod.rs | 2 +- src/gc_pcsaft/eos/dispersion.rs | 16 +- src/gc_pcsaft/eos/hard_chain.rs | 16 +- src/gc_pcsaft/eos/mod.rs | 35 ++-- src/gc_pcsaft/eos/parameter.rs | 2 +- src/gc_pcsaft/eos/polar.rs | 6 +- src/hard_sphere/dft.rs | 2 +- src/hard_sphere/mod.rs | 8 +- src/pcsaft/dft/dispersion.rs | 6 +- src/pcsaft/dft/hard_chain.rs | 2 +- src/pcsaft/dft/polar.rs | 10 +- src/pcsaft/dft/pure_saft_functional.rs | 8 +- src/pcsaft/eos/dispersion.rs | 2 +- src/pcsaft/eos/hard_chain.rs | 2 +- src/pcsaft/eos/polar.rs | 17 +- src/pcsaft/eos/qspr.rs | 2 +- src/pcsaft/parameters.rs | 2 +- src/pets/dft/dispersion.rs | 6 +- src/pets/dft/pure_pets_functional.rs | 4 +- src/pets/eos/dispersion.rs | 2 +- src/pets/eos/qspr.rs | 2 +- src/pets/parameters.rs | 2 +- src/saftvrqmie/dft/dispersion.rs | 6 +- src/saftvrqmie/dft/mod.rs | 2 +- src/saftvrqmie/dft/non_additive_hs.rs | 4 +- src/saftvrqmie/eos/dispersion.rs | 146 +++++++------- src/saftvrqmie/eos/hard_sphere.rs | 72 ++++--- src/saftvrqmie/eos/non_additive_hs.rs | 4 +- .../eos/attractive_perturbation_bh.rs | 18 +- .../eos/attractive_perturbation_uvb3.rs | 29 +-- .../eos/attractive_perturbation_wca.rs | 18 +- src/uvtheory/eos/hard_sphere_bh.rs | 23 ++- src/uvtheory/eos/hard_sphere_wca.rs | 29 +-- src/uvtheory/eos/reference_perturbation_bh.rs | 2 +- .../eos/reference_perturbation_uvb3.rs | 2 +- .../eos/reference_perturbation_wca.rs | 2 +- src/uvtheory/parameters.rs | 4 +- 63 files changed, 615 insertions(+), 522 deletions(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 73c300ab2..1f06c7322 100644 --- 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 index 8e167f82a..4ccc9a432 100644 --- 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 index 357ee68d0..ff44725de 100644 --- 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 index 0eff389fa..f78ca03f7 100644 --- a/benches/dual_numbers.rs +++ b/benches/dual_numbers.rs @@ -28,7 +28,7 @@ fn state_pcsaft(parameters: PcSaftParameters) -> State { } /// Residual Helmholtz energy given an equation of state and a StateHD. -fn a_res, E: EquationOfState>(inp: (&Arc, &StateHD)) -> D +fn a_res + Copy, E: EquationOfState>(inp: (&Arc, &StateHD)) -> D where (dyn HelmholtzEnergy + 'static): HelmholtzEnergyDual, { diff --git a/feos-core/Cargo.toml b/feos-core/Cargo.toml index 4b3dc2940..4788e6690 100644 --- 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 index 2732ea0a2..f0dc05862 100644 --- a/feos-core/src/cubic.rs +++ b/feos-core/src/cubic.rs @@ -170,7 +170,7 @@ struct PengRobinsonContribution { parameters: Arc, } -impl> HelmholtzEnergyDual for PengRobinsonContribution { +impl + Copy> HelmholtzEnergyDual for PengRobinsonContribution { fn helmholtz_energy(&self, state: &StateHD) -> 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 index abb500012..6d7b523fc 100644 --- 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> { pub trait HelmholtzEnergy: HelmholtzEnergyDual + HelmholtzEnergyDual - + HelmholtzEnergyDual, f64>> + + HelmholtzEnergyDual, f64>> + HelmholtzEnergyDual + HelmholtzEnergyDual + HelmholtzEnergyDual + HelmholtzEnergyDual> - + HelmholtzEnergyDual, f64>> - + HelmholtzEnergyDual, f64>> + + HelmholtzEnergyDual, f64>> + + HelmholtzEnergyDual, f64>> + + HelmholtzEnergyDual> + HelmholtzEnergyDual> - + HelmholtzEnergyDual, f64>> - + HelmholtzEnergyDual, f64>> + + HelmholtzEnergyDual, f64>> + + HelmholtzEnergyDual, f64>> + fmt::Display + Send + Sync @@ -47,16 +49,17 @@ pub trait HelmholtzEnergy: impl HelmholtzEnergy for T where T: HelmholtzEnergyDual + HelmholtzEnergyDual - + HelmholtzEnergyDual, f64>> + + HelmholtzEnergyDual, f64>> + HelmholtzEnergyDual + HelmholtzEnergyDual + HelmholtzEnergyDual + HelmholtzEnergyDual> - + HelmholtzEnergyDual, f64>> - + HelmholtzEnergyDual, f64>> + + HelmholtzEnergyDual, f64>> + + HelmholtzEnergyDual, f64>> + + HelmholtzEnergyDual> + HelmholtzEnergyDual> - + HelmholtzEnergyDual, f64>> - + HelmholtzEnergyDual, f64>> + + HelmholtzEnergyDual, f64>> + + HelmholtzEnergyDual, f64>> + fmt::Display + Send + Sync @@ -70,7 +73,7 @@ impl 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> { +pub trait IdealGasContributionDual + 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; @@ -101,16 +104,17 @@ pub trait IdealGasContributionDual> { pub trait IdealGasContribution: IdealGasContributionDual + IdealGasContributionDual - + IdealGasContributionDual, f64>> + + IdealGasContributionDual, f64>> + IdealGasContributionDual + IdealGasContributionDual + IdealGasContributionDual + IdealGasContributionDual> - + IdealGasContributionDual, f64>> - + IdealGasContributionDual, f64>> + + IdealGasContributionDual, f64>> + + IdealGasContributionDual, f64>> + + IdealGasContributionDual> + IdealGasContributionDual> - + IdealGasContributionDual, f64>> - + IdealGasContributionDual, f64>> + + IdealGasContributionDual, f64>> + + IdealGasContributionDual, f64>> + fmt::Display { } @@ -118,22 +122,23 @@ pub trait IdealGasContribution: impl IdealGasContribution for T where T: IdealGasContributionDual + IdealGasContributionDual - + IdealGasContributionDual, f64>> + + IdealGasContributionDual, f64>> + IdealGasContributionDual + IdealGasContributionDual + IdealGasContributionDual + IdealGasContributionDual> - + IdealGasContributionDual, f64>> - + IdealGasContributionDual, f64>> + + IdealGasContributionDual, f64>> + + IdealGasContributionDual, f64>> + + IdealGasContributionDual> + IdealGasContributionDual> - + IdealGasContributionDual, f64>> - + IdealGasContributionDual, f64>> + + IdealGasContributionDual, f64>> + + IdealGasContributionDual, f64>> + fmt::Display { } struct DefaultIdealGasContribution; -impl> IdealGasContributionDual for DefaultIdealGasContribution { +impl + Copy> IdealGasContributionDual for DefaultIdealGasContribution { fn de_broglie_wavelength(&self, _: D, components: usize) -> Array1 { Array1::zeros(components) } @@ -175,7 +180,7 @@ pub trait EquationOfState: Send + Sync { fn residual(&self) -> &[Box]; /// Evaluate the residual reduced Helmholtz energy $\beta A^\mathrm{res}$. - fn evaluate_residual>(&self, state: &StateHD) -> D + fn evaluate_residual + Copy>(&self, state: &StateHD) -> D where dyn HelmholtzEnergy: HelmholtzEnergyDual, { @@ -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>( + fn evaluate_residual_contributions + Copy>( &self, state: &StateHD, ) -> Vec<(String, D)> @@ -252,12 +257,10 @@ pub trait EquationOfState: Send + Sync { ) -> EosResult { 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 { 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 { 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| { + 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 { 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 index 643f7bc66..b8d73e90e 100644 --- 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> IdealGasContributionDual for Joback { +impl + Copy> IdealGasContributionDual for Joback { fn de_broglie_wavelength(&self, temperature: D, components: usize) -> Array1 { 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 index 52c34e10c..ace9713ac 100644 --- 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); @@ -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, f64>, + Dual, 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, PyDual64); -dual_number!(PyDualVec2, DualVec64<2>, f64); +dual_number!(PyDualVec2, DualSVec64<2>, f64); impl_dual_state_helmholtz_energy!( PyStateHDDVec2, PyHyperDualVec2, - HyperDual, f64>, + HyperDual, f64>, PyDualVec2 ); impl_dual_state_helmholtz_energy!( PyStateHDDVec3, PyHyperDualVec3, - HyperDual, f64>, + HyperDual, f64>, PyDualVec3 ); +impl_dual_state_helmholtz_energy!( + PyStateD2D, + PyDual2Dual64, + Dual2, + 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, f64>, + Dual3, f64>, PyDualVec2 ); impl_dual_state_helmholtz_energy!( PyStateD3DVec3, PyDual3DualVec3, - Dual3, f64>, + Dual3, f64>, PyDualVec3 ); diff --git a/feos-core/src/state/cache.rs b/feos-core/src/state/cache.rs index f0c41eabd..605e50459 100644 --- a/feos-core/src/state/cache.rs +++ b/feos-core/src/state/cache.rs @@ -45,8 +45,8 @@ impl Cache { let value = f(); self.map.insert(PartialDerivative::Zeroth, value.re); self.map - .insert(PartialDerivative::First(derivative), value.eps[0]); - value.eps[0] + .insert(PartialDerivative::First(derivative), value.eps.unwrap()); + value.eps.unwrap() } } @@ -66,12 +66,12 @@ impl Cache { let value = f(); self.map.insert(PartialDerivative::Zeroth, value.re); self.map - .insert(PartialDerivative::First(derivative), value.v1[0]); + .insert(PartialDerivative::First(derivative), value.v1.unwrap()); self.map.insert( PartialDerivative::SecondMixed(derivative, derivative), - value.v2[0], + value.v2.unwrap(), ); - value.v2[0] + value.v2.unwrap() } } @@ -91,14 +91,14 @@ impl Cache { let value = f(); self.map.insert(PartialDerivative::Zeroth, value.re); self.map - .insert(PartialDerivative::First(derivative1), value.eps1[0]); + .insert(PartialDerivative::First(derivative1), value.eps1.unwrap()); self.map - .insert(PartialDerivative::First(derivative2), value.eps2[0]); + .insert(PartialDerivative::First(derivative2), value.eps2.unwrap()); self.map.insert( PartialDerivative::SecondMixed(d1, d2), - value.eps1eps2[(0, 0)], + value.eps1eps2.unwrap(), ); - value.eps1eps2[(0, 0)] + value.eps1eps2.unwrap() } } diff --git a/feos-core/src/state/critical_point.rs b/feos-core/src/state/critical_point.rs index b2cb1b3c3..e1279c052 100644 --- a/feos-core/src/state/critical_point.rs +++ b/feos-core/src/state/critical_point.rs @@ -3,9 +3,12 @@ use crate::equation_of_state::EquationOfState; use crate::errors::{EosError, EosResult}; use crate::phase_equilibria::{SolverOptions, Verbosity}; use crate::{DensityInitialization, EosUnit}; -use ndarray::{arr1, arr2, Array1, Array2}; -use num_dual::linalg::{norm, smallest_ev, LU}; -use num_dual::{Dual, Dual3, Dual64, DualNum, DualVec64, HyperDual, StaticVec}; +use nalgebra::{DMatrix, DVector, SVector, SymmetricEigen}; +use ndarray::{arr1, Array1}; +use num_dual::{ + first_derivative, try_first_derivative, try_jacobian, Derivative, Dual, Dual3, Dual64, DualNum, + DualSVec64, DualVec, HyperDual, +}; use num_traits::{One, Zero}; use quantity::si::{SIArray1, SINumber, SIUnit}; use std::convert::TryFrom; @@ -123,16 +126,12 @@ impl State { for i in 1..=max_iter { // calculate residuals and derivatives w.r.t. temperature and density - let [t_dual, rho_dual] = *StaticVec::new_vec([t, rho]) - .map(DualVec64::<2>::from_re) - .derive() - .raw_array(); - let res = critical_point_objective(eos, t_dual, rho_dual, &n)?; - let h = arr2(res.jacobian().raw_data()); - let res = arr1(res.map(|r| r.re()).raw_array()); + let res = |x: SVector, 2>| critical_point_objective(eos, x[0], x[1], &n); + let (res, jac) = try_jacobian(res, SVector::from([t, rho]))?; // calculate Newton step - let mut delta = LU::new(h)?.solve(&res); + let delta = jac.lu().solve(&res); + let mut delta = delta.ok_or(EosError::IterationFailed("Critical point".into()))?; // reduce step if necessary if delta[0].abs() > 0.25 * t { @@ -151,13 +150,13 @@ impl State { verbosity, " {:4} | {:14.8e} | {:13.8} | {:12.8}", i, - norm(&res), + res.norm(), t * SIUnit::reference_temperature(), rho * SIUnit::reference_density(), ); // check convergence - if norm(&res) < tol { + if res.norm() < tol { log_result!( verbosity, "Critical point calculation converged in {} step(s)\n", @@ -188,9 +187,9 @@ impl State { options.unwrap_or(MAX_ITER_CRIT_POINT_BINARY, TOL_CRIT_POINT); let t = temperature.to_reduced(SIUnit::reference_temperature())?; - let x = StaticVec::new_vec(initial_molefracs.unwrap_or([0.5, 0.5])); + let x = SVector::from(initial_molefracs.unwrap_or([0.5, 0.5])); let max_density = eos - .max_density(Some(&(arr1(x.raw_array()) * SIUnit::reference_moles())))? + .max_density(Some(&(arr1(&x.data.0[0]) * SIUnit::reference_moles())))? .to_reduced(SIUnit::reference_density())?; let mut rho = x * 0.3 * max_density; @@ -209,17 +208,12 @@ impl State { for i in 1..=max_iter { // calculate residuals and derivatives w.r.t. partial densities - let r = StaticVec::new_vec([DualVec64::from_re(rho[0]), DualVec64::from_re(rho[1])]) - .derive(); - let res = critical_point_objective_t(eos, t, r)?; + let res = |rho| critical_point_objective_t(eos, t, rho); + let (res, jac) = try_jacobian(res, rho)?; // calculate Newton step - let h = res.jacobian(); - let res = res.map(|r| r.re); - let mut delta = StaticVec::new_vec([ - h[(1, 1)] * res[0] - h[(0, 1)] * res[1], - h[(0, 0)] * res[1] - h[(1, 0)] * res[0], - ]) / (h[(0, 0)] * h[(1, 1)] - h[(0, 1)] * h[(1, 0)]); + let delta = jac.lu().solve(&res); + let mut delta = delta.ok_or(EosError::IterationFailed("Critical point".into()))?; // reduce step if necessary for i in 0..2 { @@ -253,7 +247,7 @@ impl State { eos, t * SIUnit::reference_temperature(), SIUnit::reference_volume(), - &(arr1(rho.raw_array()) * SIUnit::reference_moles()), + &(arr1(&rho.data.0[0]) * SIUnit::reference_moles()), ); } } @@ -279,9 +273,9 @@ impl State { .map(|t| t.to_reduced(SIUnit::reference_temperature())) .transpose()? .unwrap_or(300.0); - let x = StaticVec::new_vec(initial_molefracs.unwrap_or([0.5, 0.5])); + let x = SVector::from(initial_molefracs.unwrap_or([0.5, 0.5])); let max_density = eos - .max_density(Some(&(arr1(x.raw_array()) * SIUnit::reference_moles())))? + .max_density(Some(&(arr1(&x.data.0[0]) * SIUnit::reference_moles())))? .to_reduced(SIUnit::reference_density())?; let mut rho = x * 0.3 * max_density; @@ -301,19 +295,15 @@ impl State { for i in 1..=max_iter { // calculate residuals and derivatives w.r.t. temperature and partial densities - let x = StaticVec::new_vec([ - DualVec64::from_re(t), - DualVec64::from_re(rho[0]), - DualVec64::from_re(rho[1]), - ]) - .derive(); - let r = StaticVec::new_vec([x[1], x[2]]); - let res = critical_point_objective_p(eos, p, x[0], r)?; + let res = |x: SVector, 3>| { + let r = SVector::from([x[1], x[2]]); + critical_point_objective_p(eos, p, x[0], r) + }; + let (res, jac) = try_jacobian(res, SVector::from([t, rho[0], rho[1]]))?; // calculate Newton step - let h = arr2(res.jacobian().raw_data()); - let res = arr1(res.map(|r| r.re).raw_array()); - let mut delta = LU::new(h)?.solve(&res); + let delta = jac.lu().solve(&res); + let mut delta = delta.ok_or(EosError::IterationFailed("Critical point".into()))?; // reduce step if necessary if delta[0].abs() > 0.25 * t { @@ -337,14 +327,14 @@ impl State { verbosity, " {:4} | {:14.8e} | {:13.8} | {:12.8} | {:12.8}", i, - norm(&res), + res.norm(), t * SIUnit::reference_temperature(), rho[0] * SIUnit::reference_density(), rho[1] * SIUnit::reference_density(), ); // check convergence - if norm(&res) < tol { + if res.norm() < tol { log_result!( verbosity, "Critical point calculation converged in {} step(s)\n", @@ -354,7 +344,7 @@ impl State { eos, t * SIUnit::reference_temperature(), SIUnit::reference_volume(), - &(arr1(rho.raw_array()) * SIUnit::reference_moles()), + &(arr1(&rho.data.0[0]) * SIUnit::reference_moles()), ); } } @@ -427,10 +417,11 @@ impl State { for i in 1..=max_iter { // calculate residuals and derivative w.r.t. density - let res = spinodal_objective(eos, Dual64::from(t), Dual64::from(rho).derive(), &n)?; + let (f, df) = + try_first_derivative(|rho| spinodal_objective(eos, t.into(), rho, &n), rho)?; // calculate Newton step - let mut delta = res.re / res.eps[0]; + let mut delta = f / df; // reduce step if necessary if delta.abs() > 0.03 * max_density { @@ -445,12 +436,12 @@ impl State { verbosity, " {:4} | {:14.8e} | {:12.8}", i, - res.re.abs(), + f.abs(), rho * SIUnit::reference_density(), ); // check convergence - if res.re.abs() < tol { + if f.abs() < tol { log_result!( verbosity, "Spinodal calculation converged in {} step(s)\n", @@ -470,20 +461,20 @@ impl State { fn critical_point_objective( eos: &Arc, - temperature: DualVec64<2>, - density: DualVec64<2>, + temperature: DualSVec64<2>, + density: DualSVec64<2>, moles: &Array1, -) -> EosResult, 2>> { +) -> EosResult, 2>> { // calculate second partial derivatives w.r.t. moles let t = HyperDual::from_re(temperature); let v = HyperDual::from_re(density.recip() * moles.sum()); - let qij = Array2::from_shape_fn((eos.components(), eos.components()), |(i, j)| { + let qij = DMatrix::from_fn(eos.components(), eos.components(), |i, j| { let mut m = moles.mapv(HyperDual::from); - m[i].eps1[0] = DualVec64::one(); - m[j].eps2[0] = DualVec64::one(); + m[i].eps1 = Derivative::some(SVector::from_element(DualSVec64::one())); + m[j].eps2 = Derivative::some(SVector::from_element(DualSVec64::one())); let state = StateHD::new(t, v, m); - (eos.evaluate_residual(&state).eps1eps2[(0, 0)] - + eos.ideal_gas().evaluate(&state).eps1eps2[(0, 0)]) + (eos.evaluate_residual(&state).eps1eps2.unwrap() + + eos.ideal_gas().evaluate(&state).eps1eps2.unwrap()) * (moles[i] * moles[j]).sqrt() }); @@ -493,10 +484,10 @@ fn critical_point_objective( // evaluate third partial derivative w.r.t. s let moles_hd = Array1::from_shape_fn(eos.components(), |i| { Dual3::new( - DualVec64::from(moles[i]), + DualSVec64::from(moles[i]), evec[i] * moles[i].sqrt(), - DualVec64::zero(), - DualVec64::zero(), + DualSVec64::zero(), + DualSVec64::zero(), ) }); let state_s = StateHD::new( @@ -505,24 +496,24 @@ fn critical_point_objective( moles_hd, ); let res = eos.evaluate_residual(&state_s) + eos.ideal_gas().evaluate(&state_s); - Ok(StaticVec::new_vec([eval, res.v3])) + Ok(SVector::from([eval, res.v3])) } fn critical_point_objective_t( eos: &Arc, temperature: f64, - density: StaticVec, 2>, -) -> EosResult, 2>> { + density: SVector, 2>, +) -> EosResult, 2>> { // calculate second partial derivatives w.r.t. moles let t = HyperDual::from(temperature); let v = HyperDual::from(1.0); - let qij = Array2::from_shape_fn((eos.components(), eos.components()), |(i, j)| { + let qij = DMatrix::from_fn(eos.components(), eos.components(), |i, j| { let mut m = density.map(HyperDual::from_re); - m[i].eps1[0] = DualVec64::one(); - m[j].eps2[0] = DualVec64::one(); + m[i].eps1 = Derivative::some(SVector::from_element(DualSVec64::one())); + m[j].eps2 = Derivative::some(SVector::from_element(DualSVec64::one())); let state = StateHD::new(t, v, arr1(&[m[0], m[1]])); - (eos.evaluate_residual(&state).eps1eps2[(0, 0)] - + eos.ideal_gas().evaluate(&state).eps1eps2[(0, 0)]) + (eos.evaluate_residual(&state).eps1eps2.unwrap() + + eos.ideal_gas().evaluate(&state).eps1eps2.unwrap()) * (density[i] * density[j]).sqrt() }); @@ -534,31 +525,31 @@ fn critical_point_objective_t( Dual3::new( density[i], evec[i] * density[i].sqrt(), - DualVec64::zero(), - DualVec64::zero(), + DualSVec64::zero(), + DualSVec64::zero(), ) }); let state_s = StateHD::new(Dual3::from(temperature), Dual3::from(1.0), moles_hd); let res = eos.evaluate_residual(&state_s) + eos.ideal_gas().evaluate(&state_s); - Ok(StaticVec::new_vec([eval, res.v3])) + Ok(SVector::from([eval, res.v3])) } fn critical_point_objective_p( eos: &Arc, pressure: f64, - temperature: DualVec64<3>, - density: StaticVec, 2>, -) -> EosResult, 3>> { + temperature: DualSVec64<3>, + density: SVector, 2>, +) -> EosResult, 3>> { // calculate second partial derivatives w.r.t. moles let t = HyperDual::from_re(temperature); let v = HyperDual::from(1.0); - let qij = Array2::from_shape_fn((eos.components(), eos.components()), |(i, j)| { + let qij = DMatrix::from_fn(eos.components(), eos.components(), |i, j| { let mut m = density.map(HyperDual::from_re); - m[i].eps1[0] = DualVec64::one(); - m[j].eps2[0] = DualVec64::one(); + m[i].eps1 = Derivative::some(SVector::from_element(DualSVec64::one())); + m[j].eps2 = Derivative::some(SVector::from_element(DualSVec64::one())); let state = StateHD::new(t, v, arr1(&[m[0], m[1]])); - (eos.evaluate_residual(&state).eps1eps2[(0, 0)] - + eos.ideal_gas().evaluate(&state).eps1eps2[(0, 0)]) + (eos.evaluate_residual(&state).eps1eps2.unwrap() + + eos.ideal_gas().evaluate(&state).eps1eps2.unwrap()) * (density[i] * density[j]).sqrt() }); @@ -570,24 +561,22 @@ fn critical_point_objective_p( Dual3::new( density[i], evec[i] * density[i].sqrt(), - DualVec64::zero(), - DualVec64::zero(), + DualSVec64::zero(), + DualSVec64::zero(), ) }); let state_s = StateHD::new(Dual3::from_re(temperature), Dual3::from(1.0), moles_hd); let res = eos.evaluate_residual(&state_s) + eos.ideal_gas().evaluate(&state_s); // calculate pressure - let v = Dual::from(1.0).derive(); - let m = arr1(&[Dual::from_re(density[0]), Dual::from_re(density[1])]); - let state_p = StateHD::new(Dual::from_re(temperature), v, m); - let p = eos.evaluate_residual(&state_p) + eos.ideal_gas().evaluate(&state_p); - - Ok(StaticVec::new_vec([ - eval, - res.v3, - p.eps[0] * temperature + pressure, - ])) + let a = |v| { + let m = arr1(&[Dual::from_re(density[0]), Dual::from_re(density[1])]); + let state_p = StateHD::new(Dual::from_re(temperature), v, m); + eos.evaluate_residual(&state_p) + eos.ideal_gas().evaluate(&state_p) + }; + let (_, p) = first_derivative(a, DualVec::one()); + + Ok(SVector::from([eval, res.v3, p * temperature + pressure])) } fn spinodal_objective( @@ -599,13 +588,13 @@ fn spinodal_objective( // calculate second partial derivatives w.r.t. moles let t = HyperDual::from_re(temperature); let v = HyperDual::from_re(density.recip() * moles.sum()); - let qij = Array2::from_shape_fn((eos.components(), eos.components()), |(i, j)| { + let qij = DMatrix::from_fn(eos.components(), eos.components(), |i, j| { let mut m = moles.mapv(HyperDual::from); - m[i].eps1[0] = Dual64::one(); - m[j].eps2[0] = Dual64::one(); + m[i].eps1 = Derivative::some(SVector::from_element(DualSVec64::one())); + m[j].eps2 = Derivative::some(SVector::from_element(DualSVec64::one())); let state = StateHD::new(t, v, m); - (eos.evaluate_residual(&state).eps1eps2[(0, 0)] - + eos.ideal_gas().evaluate(&state).eps1eps2[(0, 0)]) + (eos.evaluate_residual(&state).eps1eps2.unwrap() + + eos.ideal_gas().evaluate(&state).eps1eps2.unwrap()) * (moles[i] * moles[j]).sqrt() }); @@ -614,3 +603,16 @@ fn spinodal_objective( Ok(eval) } + +fn smallest_ev( + m: DMatrix>, +) -> (DualSVec64, DVector>) { + let eig = SymmetricEigen::new(m); + let (e, ev) = eig + .eigenvalues + .iter() + .zip(eig.eigenvectors.column_iter()) + .reduce(|e1, e2| if e1.0 < e2.0 { e1 } else { e2 }) + .unwrap(); + (*e, ev.into()) +} diff --git a/feos-core/src/state/mod.rs b/feos-core/src/state/mod.rs index 2133a0165..ca566add3 100644 --- a/feos-core/src/state/mod.rs +++ b/feos-core/src/state/mod.rs @@ -58,7 +58,7 @@ pub struct StateHD> { pub partial_density: Array1, } -impl> StateHD { +impl + Copy> StateHD { /// Create a new `StateHD` for given temperature volume and moles. pub fn new(temperature: D, volume: D, moles: Array1) -> Self { let total_moles = moles.sum(); @@ -686,9 +686,9 @@ impl State { let mut v = Dual64::from(self.reduced_volume); let mut n = self.reduced_moles.mapv(Dual64::from); match derivative { - Derivative::DT => t = t.derive(), - Derivative::DV => v = v.derive(), - Derivative::DN(i) => n[i] = n[i].derive(), + Derivative::DT => t = t.derivative(), + Derivative::DV => v = v.derivative(), + Derivative::DN(i) => n[i] = n[i].derivative(), } StateHD::new(t, v, n) } @@ -699,9 +699,9 @@ impl State { let mut v = Dual2_64::from(self.reduced_volume); let mut n = self.reduced_moles.mapv(Dual2_64::from); match derivative { - Derivative::DT => t = t.derive(), - Derivative::DV => v = v.derive(), - Derivative::DN(i) => n[i] = n[i].derive(), + Derivative::DT => t = t.derivative(), + Derivative::DV => v = v.derivative(), + Derivative::DN(i) => n[i] = n[i].derivative(), } StateHD::new(t, v, n) } @@ -716,14 +716,14 @@ impl State { let mut v = HyperDual64::from(self.reduced_volume); let mut n = self.reduced_moles.mapv(HyperDual64::from); match derivative1 { - Derivative::DT => t = t.derive1(), - Derivative::DV => v = v.derive1(), - Derivative::DN(i) => n[i] = n[i].derive1(), + Derivative::DT => t = t.derivative1(), + Derivative::DV => v = v.derivative1(), + Derivative::DN(i) => n[i] = n[i].derivative1(), } match derivative2 { - Derivative::DT => t = t.derive2(), - Derivative::DV => v = v.derive2(), - Derivative::DN(i) => n[i] = n[i].derive2(), + Derivative::DT => t = t.derivative2(), + Derivative::DV => v = v.derivative2(), + Derivative::DN(i) => n[i] = n[i].derivative2(), } StateHD::new(t, v, n) } @@ -734,9 +734,9 @@ impl State { let mut v = Dual3_64::from(self.reduced_volume); let mut n = self.reduced_moles.mapv(Dual3_64::from); match derivative { - Derivative::DT => t = t.derive(), - Derivative::DV => v = v.derive(), - Derivative::DN(i) => n[i] = n[i].derive(), + Derivative::DT => t = t.derivative(), + Derivative::DV => v = v.derivative(), + Derivative::DN(i) => n[i] = n[i].derivative(), }; StateHD::new(t, v, n) } diff --git a/feos-core/src/state/properties.rs b/feos-core/src/state/properties.rs index baf749718..e607be859 100644 --- a/feos-core/src/state/properties.rs +++ b/feos-core/src/state/properties.rs @@ -47,18 +47,23 @@ impl State { } PartialDerivative::First(v) => { let new_state = self.derive1(v); - -(new_state.moles.sum() * new_state.temperature * new_state.volume.ln()).eps[0] + -(new_state.moles.sum() * new_state.temperature * new_state.volume.ln()) + .eps + .unwrap() * (SIUnit::reference_energy() / v.reference()) } PartialDerivative::Second(v) => { let new_state = self.derive2(v); - -(new_state.moles.sum() * new_state.temperature * new_state.volume.ln()).v2[0] + -(new_state.moles.sum() * new_state.temperature * new_state.volume.ln()) + .v2 + .unwrap() * (SIUnit::reference_energy() / (v.reference() * v.reference())) } PartialDerivative::SecondMixed(v1, v2) => { let new_state = self.derive2_mixed(v1, v2); -(new_state.moles.sum() * new_state.temperature * new_state.volume.ln()) - .eps1eps2[(0, 0)] + .eps1eps2 + .unwrap() * (SIUnit::reference_energy() / (v1.reference() * v2.reference())) } PartialDerivative::Third(v) => { @@ -123,20 +128,25 @@ impl State { } PartialDerivative::First(v) => { let new_state = self.derive1(v); - (self.eos.ideal_gas().evaluate(&new_state) * new_state.temperature).eps[0] + (self.eos.ideal_gas().evaluate(&new_state) * new_state.temperature) + .eps + .unwrap() * SIUnit::reference_energy() / v.reference() } PartialDerivative::Second(v) => { let new_state = self.derive2(v); - (self.eos.ideal_gas().evaluate(&new_state) * new_state.temperature).v2[0] + (self.eos.ideal_gas().evaluate(&new_state) * new_state.temperature) + .v2 + .unwrap() * SIUnit::reference_energy() / (v.reference() * v.reference()) } PartialDerivative::SecondMixed(v1, v2) => { let new_state = self.derive2_mixed(v1, v2); - (self.eos.ideal_gas().evaluate(&new_state) * new_state.temperature).eps1eps2 - [(0, 0)] + (self.eos.ideal_gas().evaluate(&new_state) * new_state.temperature) + .eps1eps2 + .unwrap() * SIUnit::reference_energy() / (v1.reference() * v2.reference()) } @@ -541,13 +551,15 @@ impl State { let ig = self.eos.ideal_gas(); res.push(( ig.to_string(), - -(ig.evaluate(&new_state) * new_state.temperature).eps[0] + -(ig.evaluate(&new_state) * new_state.temperature) + .eps + .unwrap() * SIUnit::reference_pressure(), )); for (s, v) in contributions { res.push(( s, - -(v * new_state.temperature).eps[0] * SIUnit::reference_pressure(), + -(v * new_state.temperature).eps.unwrap() * SIUnit::reference_pressure(), )); } res @@ -561,13 +573,15 @@ impl State { let ig = self.eos.ideal_gas(); res.push(( ig.to_string(), - (ig.evaluate(&new_state) * new_state.temperature).eps[0] + (ig.evaluate(&new_state) * new_state.temperature) + .eps + .unwrap() * SIUnit::reference_molar_energy(), )); for (s, v) in contributions { res.push(( s, - (v * new_state.temperature).eps[0] * SIUnit::reference_molar_energy(), + (v * new_state.temperature).eps.unwrap() * SIUnit::reference_molar_energy(), )); } res diff --git a/feos-dft/Cargo.toml b/feos-dft/Cargo.toml index 02e7145d4..f9f1b7af2 100644 --- a/feos-dft/Cargo.toml +++ b/feos-dft/Cargo.toml @@ -18,9 +18,10 @@ features = [ "rayon" ] [dependencies] quantity = { version = "0.6", features = ["linalg"] } -num-dual = "0.6" +num-dual = "0.7" feos-core = { version = "0.4", path = "../feos-core" } ndarray = "0.15" +nalgebra = "0.32" rustdct = "0.7" rustfft = "6.0" ang = "0.6" diff --git a/feos-dft/src/convolver/mod.rs b/feos-dft/src/convolver/mod.rs index 5cdaa56d2..f609916cc 100644 --- a/feos-dft/src/convolver/mod.rs +++ b/feos-dft/src/convolver/mod.rs @@ -40,7 +40,7 @@ pub(crate) struct BulkConvolver { weight_constants: Vec>, } -impl> BulkConvolver { +impl + Copy + Send + Sync> BulkConvolver { pub(crate) fn new(weight_functions: Vec>) -> Arc> { let weight_constants = weight_functions .into_iter() @@ -50,7 +50,7 @@ impl> BulkConvolver { } } -impl> Convolver for BulkConvolver +impl + Copy + Send + Sync> Convolver for BulkConvolver where Array2: Dot, Output = Array1>, { diff --git a/feos-dft/src/functional.rs b/feos-dft/src/functional.rs index 2c037060a..88f06702f 100644 --- a/feos-dft/src/functional.rs +++ b/feos-dft/src/functional.rs @@ -51,7 +51,7 @@ impl MolarWeight for DFT { } struct DefaultIdealGasContribution(); -impl> IdealGasContributionDual for DefaultIdealGasContribution { +impl + Copy> IdealGasContributionDual for DefaultIdealGasContribution { fn de_broglie_wavelength(&self, _: D, components: usize) -> Array1 { Array1::zeros(components) } @@ -80,7 +80,7 @@ impl EquationOfState for DFT { unreachable!() } - fn evaluate_residual>(&self, state: &StateHD) -> D + fn evaluate_residual + Copy>(&self, state: &StateHD) -> D where dyn HelmholtzEnergy: HelmholtzEnergyDual, { @@ -91,7 +91,7 @@ impl EquationOfState for DFT { + self.ideal_chain_contribution().helmholtz_energy(state) } - fn evaluate_residual_contributions>( + fn evaluate_residual_contributions + Copy>( &self, state: &StateHD, ) -> Vec<(String, D)> @@ -278,7 +278,7 @@ impl DFT { convolver: &Arc>, ) -> EosResult> where - N: DualNum + ScalarOperand, + N: DualNum + Copy + ScalarOperand, dyn FunctionalContribution: FunctionalContributionDual, D: Dimension, D::Larger: Dimension, @@ -318,7 +318,7 @@ impl DFT { D: Dimension, D::Larger: Dimension, { - let temperature_dual = Dual64::from(temperature).derive(); + let temperature_dual = Dual64::from(temperature).derivative(); let mut helmholtz_energy_density = self.intrinsic_helmholtz_energy_density(temperature_dual, density, convolver)?; match contributions { @@ -328,7 +328,7 @@ impl DFT { Contributions::ResidualNpt|Contributions::IdealGas => panic!("Entropy density can only be calculated for Contributions::Residual or Contributions::Total"), Contributions::ResidualNvt => (), } - Ok(helmholtz_energy_density.mapv(|f| -f.eps[0])) + Ok(helmholtz_energy_density.mapv(|f| -f.eps.unwrap())) } /// Calculate the individual contributions to the entropy density. @@ -346,7 +346,7 @@ impl DFT { ::Larger: Dimension, { let density_dual = density.mapv(Dual64::from); - let temperature_dual = Dual64::from(temperature).derive(); + let temperature_dual = Dual64::from(temperature).derivative(); let weighted_densities = convolver.weighted_densities(&density_dual); let functional_contributions = self.contributions(); let mut helmholtz_energy_density: Vec> = @@ -370,7 +370,7 @@ impl DFT { } Ok(helmholtz_energy_density .iter() - .map(|v| v.mapv(|f| -(f * temperature_dual).eps[0])) + .map(|v| v.mapv(|f| -(f * temperature_dual).eps.unwrap())) .collect()) } @@ -389,7 +389,7 @@ impl DFT { D: Dimension, D::Larger: Dimension, { - let temperature_dual = Dual64::from(temperature).derive(); + let temperature_dual = Dual64::from(temperature).derivative(); let mut helmholtz_energy_density_dual = self.intrinsic_helmholtz_energy_density(temperature_dual, density, convolver)?; match contributions { @@ -400,7 +400,7 @@ impl DFT { Contributions::ResidualNvt => (), } let helmholtz_energy_density = helmholtz_energy_density_dual - .mapv(|f| f.re - f.eps[0] * temperature) + .mapv(|f| f.re - f.eps.unwrap() * temperature) + (external_potential * density).sum_axis(Axis(0)) * temperature; Ok(helmholtz_energy_density) } @@ -452,7 +452,7 @@ impl DFT { D: Dimension, D::Larger: Dimension, { - let temperature_dual = Dual64::from(temperature).derive(); + let temperature_dual = Dual64::from(temperature).derivative(); let density_dual = density.mapv(Dual64::from); let weighted_densities = convolver.weighted_densities(&density_dual); let contributions = self.contributions(); diff --git a/feos-dft/src/functional_contribution.rs b/feos-dft/src/functional_contribution.rs index 67541246d..33b75e7a5 100644 --- a/feos-dft/src/functional_contribution.rs +++ b/feos-dft/src/functional_contribution.rs @@ -3,7 +3,7 @@ use feos_core::{EosResult, HelmholtzEnergyDual, StateHD}; use ndarray::prelude::*; use ndarray::RemoveAxis; use num_dual::*; -use num_traits::{One, Zero}; +use num_traits::Zero; use std::fmt::Display; macro_rules! impl_helmholtz_energy { @@ -34,16 +34,17 @@ macro_rules! impl_helmholtz_energy { impl_helmholtz_energy!(f64); impl_helmholtz_energy!(Dual64); -impl_helmholtz_energy!(Dual, f64>); +impl_helmholtz_energy!(Dual, f64>); impl_helmholtz_energy!(HyperDual64); impl_helmholtz_energy!(Dual2_64); impl_helmholtz_energy!(Dual3_64); impl_helmholtz_energy!(HyperDual); -impl_helmholtz_energy!(HyperDual, f64>); -impl_helmholtz_energy!(HyperDual, f64>); +impl_helmholtz_energy!(HyperDual, f64>); +impl_helmholtz_energy!(HyperDual, f64>); +impl_helmholtz_energy!(Dual2); impl_helmholtz_energy!(Dual3); -impl_helmholtz_energy!(Dual3, f64>); -impl_helmholtz_energy!(Dual3, f64>); +impl_helmholtz_energy!(Dual3, f64>); +impl_helmholtz_energy!(Dual3, f64>); /// Individual functional contribution that can /// be evaluated using generalized (hyper) dual numbers. @@ -76,16 +77,17 @@ pub trait FunctionalContribution: FunctionalContributionDual + FunctionalContributionDual + FunctionalContributionDual> - + FunctionalContributionDual, f64>> + + FunctionalContributionDual, f64>> + FunctionalContributionDual + FunctionalContributionDual + FunctionalContributionDual + FunctionalContributionDual> - + FunctionalContributionDual, f64>> - + FunctionalContributionDual, f64>> + + FunctionalContributionDual, f64>> + + FunctionalContributionDual, f64>> + + FunctionalContributionDual> + FunctionalContributionDual> - + FunctionalContributionDual, f64>> - + FunctionalContributionDual, f64>> + + FunctionalContributionDual, f64>> + + FunctionalContributionDual, f64>> + Display + Sync + Send @@ -103,13 +105,13 @@ pub trait FunctionalContribution: for i in 0..wd.shape()[0] { wd.index_axis_mut(Axis(0), i) - .map_inplace(|x| x.eps[0] = 1.0); + .map_inplace(|x| x.eps = Derivative::derivative()); phi = self.calculate_helmholtz_energy_density(t, wd.view())?; first_partial_derivative .index_axis_mut(Axis(0), i) - .assign(&phi.mapv(|p| p.eps[0])); + .assign(&phi.mapv(|p| p.eps.unwrap())); wd.index_axis_mut(Axis(0), i) - .map_inplace(|x| x.eps[0] = 0.0); + .map_inplace(|x| x.eps = Derivative::none()); } helmholtz_energy_density.assign(&phi.mapv(|p| p.re)); Ok(()) @@ -128,13 +130,13 @@ pub trait FunctionalContribution: for i in 0..wd.shape()[0] { wd.index_axis_mut(Axis(0), i) - .map_inplace(|x| x.eps[0] = Dual64::one()); + .map_inplace(|x| x.eps = Derivative::derivative()); phi = self.calculate_helmholtz_energy_density(t, wd.view())?; first_partial_derivative .index_axis_mut(Axis(0), i) - .assign(&phi.mapv(|p| p.eps[0])); + .assign(&phi.mapv(|p| p.eps.unwrap())); wd.index_axis_mut(Axis(0), i) - .map_inplace(|x| x.eps[0] = Dual64::zero()); + .map_inplace(|x| x.eps = Derivative::none()); } helmholtz_energy_density.assign(&phi.mapv(|p| p.re)); Ok(()) @@ -154,12 +156,12 @@ pub trait FunctionalContribution: for i in 0..wd.shape()[0] { wd.index_axis_mut(Axis(0), i) - .map_inplace(|x| x.eps1[0] = 1.0); + .map_inplace(|x| x.eps1 = Derivative::derivative()); for j in 0..=i { wd.index_axis_mut(Axis(0), j) - .map_inplace(|x| x.eps2[0] = 1.0); + .map_inplace(|x| x.eps2 = Derivative::derivative()); phi = self.calculate_helmholtz_energy_density(t, wd.view())?; - let p = phi.mapv(|p| p.eps1eps2[(0, 0)]); + let p = phi.mapv(|p| p.eps1eps2.unwrap()); second_partial_derivative .index_axis_mut(Axis(0), i) .index_axis_mut(Axis(0), j) @@ -171,13 +173,13 @@ pub trait FunctionalContribution: .assign(&p); } wd.index_axis_mut(Axis(0), j) - .map_inplace(|x| x.eps2[0] = 0.0); + .map_inplace(|x| x.eps2 = Derivative::none()); } first_partial_derivative .index_axis_mut(Axis(0), i) - .assign(&phi.mapv(|p| p.eps1[0])); + .assign(&phi.mapv(|p| p.eps1.unwrap())); wd.index_axis_mut(Axis(0), i) - .map_inplace(|x| x.eps1[0] = 0.0); + .map_inplace(|x| x.eps1 = Derivative::none()); } helmholtz_energy_density.assign(&phi.mapv(|p| p.re)); Ok(()) @@ -188,16 +190,17 @@ impl FunctionalContribution for T where T: FunctionalContributionDual + FunctionalContributionDual + FunctionalContributionDual> - + FunctionalContributionDual, f64>> + + FunctionalContributionDual, f64>> + FunctionalContributionDual + FunctionalContributionDual + FunctionalContributionDual + FunctionalContributionDual> - + FunctionalContributionDual, f64>> - + FunctionalContributionDual, f64>> + + FunctionalContributionDual, f64>> + + FunctionalContributionDual, f64>> + + FunctionalContributionDual> + FunctionalContributionDual> - + FunctionalContributionDual, f64>> - + FunctionalContributionDual, f64>> + + FunctionalContributionDual, f64>> + + FunctionalContributionDual, f64>> + Display + Sync + Send diff --git a/feos-dft/src/ideal_chain_contribution.rs b/feos-dft/src/ideal_chain_contribution.rs index 9eece53aa..39847458e 100644 --- a/feos-dft/src/ideal_chain_contribution.rs +++ b/feos-dft/src/ideal_chain_contribution.rs @@ -19,7 +19,7 @@ impl IdealChainContribution { } } -impl> HelmholtzEnergyDual for IdealChainContribution { +impl + Copy> HelmholtzEnergyDual for IdealChainContribution { fn helmholtz_energy(&self, state: &StateHD) -> D { let segments = self.component_index.len(); if self.component_index[segments - 1] + 1 != segments { diff --git a/feos-dft/src/pdgt.rs b/feos-dft/src/pdgt.rs index bdea179f2..9c6bbf4b9 100644 --- a/feos-dft/src/pdgt.rs +++ b/feos-dft/src/pdgt.rs @@ -3,19 +3,18 @@ use super::functional_contribution::FunctionalContribution; use super::weight_functions::WeightFunctionInfo; use feos_core::{Contributions, EosResult, EosUnit, EquationOfState, PhaseEquilibrium}; use ndarray::*; -use num_dual::HyperDual64; -// use quantity::{SIArray2, SINumber}; +use num_dual::Dual2_64; use quantity::si::{SIArray1, SIArray2, SINumber, SIUnit}; use std::ops::AddAssign; -impl WeightFunctionInfo { +impl WeightFunctionInfo { fn pdgt_weight_constants(&self) -> (Array2, Array2, Array2) { - let k = HyperDual64::from(0.0).derive1().derive2(); + let k = Dual2_64::from(0.0).derivative(); let w = self.weight_constants(k, 1); ( w.mapv(|w| w.re), - w.mapv(|w| -w.eps1[0]), - w.mapv(|w| -0.5 * w.eps1eps2[(0, 0)]), + w.mapv(|w| -w.v1.unwrap()), + w.mapv(|w| -0.5 * w.v2.unwrap()), ) } } @@ -32,7 +31,7 @@ impl dyn FunctionalContribution { influence_matrix: Option<&mut Array3>, ) -> EosResult<()> { // calculate weighted densities - let weight_functions = self.weight_functions_pdgt(HyperDual64::from(temperature)); + let weight_functions = self.weight_functions_pdgt(Dual2_64::from(temperature)); let (w0, w1, w2) = weight_functions.pdgt_weight_constants(); let weighted_densities = w0.dot(density); diff --git a/feos-dft/src/profile.rs b/feos-dft/src/profile.rs index 86ab4479e..42902c0d4 100644 --- a/feos-dft/src/profile.rs +++ b/feos-dft/src/profile.rs @@ -485,7 +485,7 @@ where let functional_contributions = self.dft.contributions(); let weight_functions: Vec> = functional_contributions .iter() - .map(|c| c.weight_functions(Dual64::from(t).derive())) + .map(|c| c.weight_functions(Dual64::from(t).derivative())) .collect(); let convolver = ConvolverFFT::plan(&self.grid, &weight_functions, None); @@ -518,7 +518,7 @@ where let functional_contributions = self.dft.contributions(); let weight_functions: Vec> = functional_contributions .iter() - .map(|c| c.weight_functions(Dual64::from(t).derive())) + .map(|c| c.weight_functions(Dual64::from(t).derivative())) .collect(); let convolver = ConvolverFFT::plan(&self.grid, &weight_functions, None); @@ -617,7 +617,7 @@ where let functional_contributions = self.dft.contributions(); let weight_functions: Vec> = functional_contributions .iter() - .map(|c| c.weight_functions(Dual64::from(t).derive())) + .map(|c| c.weight_functions(Dual64::from(t).derivative())) .collect(); let convolver: Arc> = ConvolverFFT::plan(&self.grid, &weight_functions, None); @@ -638,11 +638,11 @@ where let x = (self.bulk.partial_molar_volume(Contributions::Total) * self.bulk.dp_dt(Contributions::Total)) .to_reduced(SIUnit::reference_molar_entropy())?; - let mut lhs = dfdrhodt.mapv(|d| d.eps[0]); + let mut lhs = dfdrhodt.mapv(|d| d.eps.unwrap()); lhs.outer_iter_mut() .zip(dfdrhodt_bulk.into_iter()) .zip(x.into_iter()) - .for_each(|((mut lhs, d), x)| lhs -= d.eps[0] - x); + .for_each(|((mut lhs, d), x)| lhs -= d.eps.unwrap() - x); lhs.outer_iter_mut() .zip(rho.outer_iter()) .zip(rho_bulk.into_iter()) diff --git a/feos-dft/src/solver.rs b/feos-dft/src/solver.rs index 3da93b1f8..47bd4851d 100644 --- a/feos-dft/src/solver.rs +++ b/feos-dft/src/solver.rs @@ -1,8 +1,8 @@ use crate::{DFTProfile, HelmholtzEnergyFunctional, WeightFunction, WeightFunctionShape}; use feos_core::{log_iter, log_result, EosError, EosResult, EosUnit, Verbosity}; +use nalgebra::{DMatrix, DVector}; use ndarray::prelude::*; use ndarray::RemoveAxis; -use num_dual::linalg::LU; use petgraph::graph::Graph; use petgraph::visit::EdgeRef; use petgraph::Directed; @@ -398,7 +398,7 @@ where } // calculate alpha - r = Array::from_shape_fn((m + 1, m + 1), |(i, j)| match (i == m, j == m) { + r = DMatrix::from_fn(m + 1, m + 1, |i, j| match (i == m, j == m) { (false, false) => { let (resi, resi_bulk, _) = &resm[i]; let (resj, resj_bulk, _) = &resm[j]; @@ -407,9 +407,10 @@ where (true, true) => 0.0, _ => 1.0, }); - alpha = Array::zeros(m + 1); + alpha = DVector::zeros(m + 1); alpha[m] = 1.0; - alpha = LU::new(r)?.solve(&alpha); + let alpha = r.lu().solve(&alpha); + let alpha = alpha.ok_or(EosError::Error("alpha matrix is not invertible".into()))?; // update solution rho.fill(0.0); diff --git a/feos-dft/src/weight_functions.rs b/feos-dft/src/weight_functions.rs index 871a0b59d..4f4c0665e 100644 --- a/feos-dft/src/weight_functions.rs +++ b/feos-dft/src/weight_functions.rs @@ -15,7 +15,7 @@ pub struct WeightFunction { pub shape: WeightFunctionShape, } -impl> WeightFunction { +impl + Copy> WeightFunction { /// Create a new weight function without prefactor pub fn new_unscaled(kernel_radius: Array1, shape: WeightFunctionShape) -> Self { Self { @@ -276,7 +276,7 @@ impl WeightFunctionInfo { } } -impl> WeightFunctionInfo { +impl + Copy> WeightFunctionInfo { /// calculates the matrix of weight constants for this set of weighted densities pub fn weight_constants(&self, k: T, dimensions: usize) -> Array2 { let segments = self.component_index.len(); diff --git a/src/association/dft.rs b/src/association/dft.rs index 05c92b378..5a8838a2c 100644 --- a/src/association/dft.rs +++ b/src/association/dft.rs @@ -12,7 +12,7 @@ pub const N0_CUTOFF: f64 = 1e-9; impl FunctionalContributionDual for Association

where - N: DualNum + ScalarOperand, + N: DualNum + Copy + ScalarOperand, P: HardSphereProperties, { fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { @@ -111,7 +111,7 @@ where impl Association

{ pub fn calculate_helmholtz_energy_density< - N: DualNum + ScalarOperand, + N: DualNum + Copy + ScalarOperand, S: Data, >( &self, @@ -183,7 +183,10 @@ impl Association

{ } } - fn helmholtz_energy_density_ab_analytic + ScalarOperand, S: Data>( + fn helmholtz_energy_density_ab_analytic< + N: DualNum + Copy + ScalarOperand, + S: Data, + >( &self, temperature: N, rho0: &Array2, @@ -217,7 +220,10 @@ impl Association

{ rhoa * xa.mapv(f) + rhob * xb.mapv(f) } - fn helmholtz_energy_density_cc_analytic + ScalarOperand, S: Data>( + fn helmholtz_energy_density_cc_analytic< + N: DualNum + Copy + ScalarOperand, + S: Data, + >( &self, temperature: N, rho0: &Array2, diff --git a/src/association/mod.rs b/src/association/mod.rs index 8dfcbc9de..3a9d7385c 100644 --- a/src/association/mod.rs +++ b/src/association/mod.rs @@ -212,7 +212,7 @@ impl Association

{ res } - fn association_strength>( + fn association_strength + Copy>( &self, temperature: D, diameter: &Array1, @@ -241,7 +241,7 @@ impl Association

{ } } -impl + ScalarOperand, P: HardSphereProperties> HelmholtzEnergyDual +impl + Copy + ScalarOperand, P: HardSphereProperties> HelmholtzEnergyDual for Association

{ fn helmholtz_energy(&self, state: &StateHD) -> D { @@ -305,7 +305,11 @@ impl

fmt::Display for Association

{ } impl Association

{ - fn helmholtz_energy_ab_analytic>(&self, state: &StateHD, delta: D) -> D { + fn helmholtz_energy_ab_analytic + Copy>( + &self, + state: &StateHD, + delta: D, + ) -> D { let a = &self.association_parameters; // site densities @@ -322,7 +326,11 @@ impl Association

{ (rhoa * (xa.ln() - xa * 0.5 + 0.5) + rhob * (xb.ln() - xb * 0.5 + 0.5)) * state.volume } - fn helmholtz_energy_cc_analytic>(&self, state: &StateHD, delta: D) -> D { + fn helmholtz_energy_cc_analytic + Copy>( + &self, + state: &StateHD, + delta: D, + ) -> D { let a = &self.association_parameters; // site density @@ -337,7 +345,7 @@ impl Association

{ #[allow(clippy::too_many_arguments)] fn helmholtz_energy_density_cross_association< - D: DualNum + ScalarOperand, + D: DualNum + Copy + ScalarOperand, S: Data, >( rho: &ArrayBase, @@ -396,7 +404,7 @@ impl Association

{ Ok((rho * x_dual.mapv(f)).sum()) } - fn newton_step_cross_association + ScalarOperand, S: Data>( + fn newton_step_cross_association + Copy + ScalarOperand, S: Data>( x: &mut Array1, delta_ab: &Array2, delta_cc: &Array2, @@ -541,11 +549,12 @@ mod tests_gc_pcsaft { let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); let state = StateHD::new( Dual64::from_re(temperature), - Dual64::from_re(volume).derive(), + Dual64::from_re(volume).derivative(), arr1(&[Dual64::from_re(moles)]), ); - let pressure = - -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + let pressure = -contrib.helmholtz_energy(&state).eps.unwrap() + * temperature + * EosUnit::reference_pressure(); assert_relative_eq!(pressure, -3.6819598891967344 * PASCAL, max_relative = 1e-10); } @@ -561,11 +570,12 @@ mod tests_gc_pcsaft { let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); let state = StateHD::new( Dual64::from_re(temperature), - Dual64::from_re(volume).derive(), + Dual64::from_re(volume).derivative(), arr1(&[Dual64::from_re(moles)]), ); - let pressure = - -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + let pressure = -contrib.helmholtz_energy(&state).eps.unwrap() + * temperature + * EosUnit::reference_pressure(); assert_relative_eq!(pressure, -3.6819598891967344 * PASCAL, max_relative = 1e-10); } @@ -583,11 +593,12 @@ mod tests_gc_pcsaft { .unwrap(); let state = StateHD::new( Dual64::from_re(temperature), - Dual64::from_re(volume).derive(), + Dual64::from_re(volume).derivative(), moles.mapv(Dual64::from_re), ); - let pressure = - -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + let pressure = -contrib.helmholtz_energy(&state).eps.unwrap() + * temperature + * EosUnit::reference_pressure(); assert_relative_eq!(pressure, -26.105606376765632 * PASCAL, max_relative = 1e-10); } } diff --git a/src/gc_pcsaft/dft/dispersion.rs b/src/gc_pcsaft/dft/dispersion.rs index 4eade834f..f9baea129 100644 --- a/src/gc_pcsaft/dft/dispersion.rs +++ b/src/gc_pcsaft/dft/dispersion.rs @@ -24,7 +24,9 @@ impl AttractiveFunctional { } } -impl + ScalarOperand> FunctionalContributionDual for AttractiveFunctional { +impl + Copy + ScalarOperand> FunctionalContributionDual + for AttractiveFunctional +{ fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { let p = &self.parameters; diff --git a/src/gc_pcsaft/dft/hard_chain.rs b/src/gc_pcsaft/dft/hard_chain.rs index b7dc74793..c3552821f 100644 --- a/src/gc_pcsaft/dft/hard_chain.rs +++ b/src/gc_pcsaft/dft/hard_chain.rs @@ -23,7 +23,7 @@ impl ChainFunctional { } } -impl + ScalarOperand> FunctionalContributionDual for ChainFunctional { +impl + Copy + ScalarOperand> FunctionalContributionDual for ChainFunctional { fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { let p = &self.parameters; let d = p.hs_diameter(temperature); diff --git a/src/gc_pcsaft/dft/mod.rs b/src/gc_pcsaft/dft/mod.rs index e8e502d65..2b0f43b38 100644 --- a/src/gc_pcsaft/dft/mod.rs +++ b/src/gc_pcsaft/dft/mod.rs @@ -128,7 +128,7 @@ impl HardSphereProperties for GcPcSaftFunctionalParameters { MonomerShape::Heterosegmented([m.clone(), m.clone(), m.clone(), m], &self.component_index) } - fn hs_diameter>(&self, temperature: D) -> Array1 { + fn hs_diameter + Copy>(&self, temperature: D) -> Array1 { let ti = temperature.recip() * -3.0; Array1::from_shape_fn(self.sigma.len(), |i| { -((ti * self.epsilon_k[i]).exp() * 0.12 - 1.0) * self.sigma[i] diff --git a/src/gc_pcsaft/eos/dispersion.rs b/src/gc_pcsaft/eos/dispersion.rs index a104b35c7..b79896d19 100644 --- a/src/gc_pcsaft/eos/dispersion.rs +++ b/src/gc_pcsaft/eos/dispersion.rs @@ -66,7 +66,7 @@ pub struct Dispersion { pub parameters: Arc, } -impl> HelmholtzEnergyDual for Dispersion { +impl + Copy> HelmholtzEnergyDual for Dispersion { fn helmholtz_energy(&self, state: &StateHD) -> D { // auxiliary variables let p = &self.parameters; @@ -151,11 +151,12 @@ mod test { let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); let state = StateHD::new( Dual64::from_re(temperature), - Dual64::from_re(volume).derive(), + Dual64::from_re(volume).derivative(), arr1(&[Dual64::from_re(moles)]), ); - let pressure = - -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + let pressure = -contrib.helmholtz_energy(&state).eps.unwrap() + * temperature + * EosUnit::reference_pressure(); assert_relative_eq!(pressure, -2.846724434944439 * PASCAL, max_relative = 1e-10); } @@ -173,11 +174,12 @@ mod test { let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); let state = StateHD::new( Dual64::from_re(temperature), - Dual64::from_re(volume).derive(), + Dual64::from_re(volume).derivative(), arr1(&[Dual64::from_re(moles)]), ); - let pressure = - -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + let pressure = -contrib.helmholtz_energy(&state).eps.unwrap() + * temperature + * EosUnit::reference_pressure(); assert_relative_eq!(pressure, -5.432173507270732 * PASCAL, max_relative = 1e-10); } } diff --git a/src/gc_pcsaft/eos/hard_chain.rs b/src/gc_pcsaft/eos/hard_chain.rs index 419ac92bc..bb996a37f 100644 --- a/src/gc_pcsaft/eos/hard_chain.rs +++ b/src/gc_pcsaft/eos/hard_chain.rs @@ -10,7 +10,7 @@ pub struct HardChain { pub parameters: Arc, } -impl> HelmholtzEnergyDual for HardChain { +impl + Copy> HelmholtzEnergyDual for HardChain { fn helmholtz_energy(&self, state: &StateHD) -> D { // temperature dependent segment diameter let diameter = self.parameters.hs_diameter(state.temperature); @@ -66,11 +66,12 @@ mod test { let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); let state = StateHD::new( Dual64::from_re(temperature), - Dual64::from_re(volume).derive(), + Dual64::new_scalar(volume, 1.0), arr1(&[Dual64::from_re(moles)]), ); - let pressure = - -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + let pressure = -contrib.helmholtz_energy(&state).eps.unwrap() + * temperature + * EosUnit::reference_pressure(); assert_relative_eq!( pressure, -7.991735636207462e-1 * PASCAL, @@ -92,11 +93,12 @@ mod test { let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); let state = StateHD::new( Dual64::from_re(temperature), - Dual64::from_re(volume).derive(), + Dual64::from_re(volume).derivative(), arr1(&[Dual64::from_re(moles)]), ); - let pressure = - -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + let pressure = -contrib.helmholtz_energy(&state).eps.unwrap() + * temperature + * EosUnit::reference_pressure(); assert_relative_eq!(pressure, -1.2831486124723626 * PASCAL, max_relative = 1e-10); } } diff --git a/src/gc_pcsaft/eos/mod.rs b/src/gc_pcsaft/eos/mod.rs index eec50fea4..3c821cb33 100644 --- a/src/gc_pcsaft/eos/mod.rs +++ b/src/gc_pcsaft/eos/mod.rs @@ -139,11 +139,12 @@ mod test { let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); let state = StateHD::new( Dual64::from_re(temperature), - Dual64::from_re(volume).derive(), + Dual64::from_re(volume).derivative(), arr1(&[Dual64::from_re(moles)]), ); - let pressure = - -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + let pressure = -contrib.helmholtz_energy(&state).eps.unwrap() + * temperature + * EosUnit::reference_pressure(); assert_relative_eq!(pressure, 1.5285037907989527 * PASCAL, max_relative = 1e-10); } @@ -159,11 +160,12 @@ mod test { let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); let state = StateHD::new( Dual64::from_re(temperature), - Dual64::from_re(volume).derive(), + Dual64::from_re(volume).derivative(), arr1(&[Dual64::from_re(moles)]), ); - let pressure = - -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + let pressure = -contrib.helmholtz_energy(&state).eps.unwrap() + * temperature + * EosUnit::reference_pressure(); assert_relative_eq!(pressure, 2.3168212018200243 * PASCAL, max_relative = 1e-10); } @@ -179,11 +181,12 @@ mod test { let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); let state = StateHD::new( Dual64::from_re(temperature), - Dual64::from_re(volume).derive(), + Dual64::from_re(volume).derivative(), arr1(&[Dual64::from_re(moles)]), ); - let pressure = - -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + let pressure = -contrib.helmholtz_energy(&state).eps.unwrap() + * temperature + * EosUnit::reference_pressure(); assert_relative_eq!(pressure, -3.6819598891967344 * PASCAL, max_relative = 1e-10); } @@ -200,11 +203,12 @@ mod test { let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); let state = StateHD::new( Dual64::from_re(temperature), - Dual64::from_re(volume).derive(), + Dual64::from_re(volume).derivative(), arr1(&[Dual64::from_re(moles)]), ); - let pressure = - -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + let pressure = -contrib.helmholtz_energy(&state).eps.unwrap() + * temperature + * EosUnit::reference_pressure(); assert_relative_eq!(pressure, -3.6819598891967344 * PASCAL, max_relative = 1e-10); } @@ -222,11 +226,12 @@ mod test { .unwrap(); let state = StateHD::new( Dual64::from_re(temperature), - Dual64::from_re(volume).derive(), + Dual64::from_re(volume).derivative(), moles.mapv(Dual64::from_re), ); - let pressure = - -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + let pressure = -contrib.helmholtz_energy(&state).eps.unwrap() + * temperature + * EosUnit::reference_pressure(); assert_relative_eq!(pressure, -26.105606376765632 * PASCAL, max_relative = 1e-10); } } diff --git a/src/gc_pcsaft/eos/parameter.rs b/src/gc_pcsaft/eos/parameter.rs index 795ea2792..c7567ed81 100644 --- a/src/gc_pcsaft/eos/parameter.rs +++ b/src/gc_pcsaft/eos/parameter.rs @@ -277,7 +277,7 @@ impl HardSphereProperties for GcPcSaftEosParameters { MonomerShape::Heterosegmented([m.clone(), m.clone(), m.clone(), m], &self.component_index) } - fn hs_diameter>(&self, temperature: D) -> Array1 { + fn hs_diameter + Copy>(&self, temperature: D) -> Array1 { let ti = temperature.recip() * -3.0; Array1::from_shape_fn(self.sigma.len(), |i| { -((ti * self.epsilon_k[i]).exp() * 0.12 - 1.0) * self.sigma[i] diff --git a/src/gc_pcsaft/eos/polar.rs b/src/gc_pcsaft/eos/polar.rs index e72dc2f37..b36f47b59 100644 --- a/src/gc_pcsaft/eos/polar.rs +++ b/src/gc_pcsaft/eos/polar.rs @@ -33,7 +33,7 @@ pub const CD: [[f64; 3]; 4] = [ pub const PI_SQ_43: f64 = 4.0 * PI * FRAC_PI_3; -fn pair_integral_ij>(mij1: f64, mij2: f64, eta: D, eps_ij_t: D) -> D { +fn pair_integral_ij + Copy>(mij1: f64, mij2: f64, eta: D, eps_ij_t: D) -> D { let eta2 = eta * eta; let etas = [D::one(), eta, eta2, eta2 * eta, eta2 * eta2]; (0..AD.len()) @@ -45,7 +45,7 @@ fn pair_integral_ij>(mij1: f64, mij2: f64, eta: D, eps_ij_t: D) .sum() } -fn triplet_integral_ijk>(mijk1: f64, mijk2: f64, eta: D) -> D { +fn triplet_integral_ijk + Copy>(mijk1: f64, mijk2: f64, eta: D) -> D { let eta2 = eta * eta; let etas = [D::one(), eta, eta2, eta2 * eta]; (0..CD.len()) @@ -118,7 +118,7 @@ impl Dipole { } } -impl> HelmholtzEnergyDual for Dipole { +impl + Copy> HelmholtzEnergyDual for Dipole { fn helmholtz_energy(&self, state: &StateHD) -> D { let p = &self.parameters; let ndipole = p.dipole_comp.len(); diff --git a/src/hard_sphere/dft.rs b/src/hard_sphere/dft.rs index b261b41dd..c886e2c95 100644 --- a/src/hard_sphere/dft.rs +++ b/src/hard_sphere/dft.rs @@ -84,7 +84,7 @@ impl

FMTContribution

{ } } -impl> FunctionalContributionDual +impl + Copy> FunctionalContributionDual for FMTContribution

{ fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { diff --git a/src/hard_sphere/mod.rs b/src/hard_sphere/mod.rs index 92d96e7a2..51a9da45d 100644 --- a/src/hard_sphere/mod.rs +++ b/src/hard_sphere/mod.rs @@ -31,7 +31,7 @@ pub trait HardSphereProperties { fn monomer_shape>(&self, temperature: D) -> MonomerShape; /// The temperature dependent hard-sphere diameters of every segment. - fn hs_diameter>(&self, temperature: D) -> Array1; + fn hs_diameter + Copy>(&self, temperature: D) -> Array1; /// For every segment, the index of the component that it is on. fn component_index(&self) -> Cow> { @@ -55,7 +55,7 @@ pub trait HardSphereProperties { } /// The packing fractions $\zeta_k$. - fn zeta, const N: usize>( + fn zeta + Copy, const N: usize>( &self, temperature: D, partial_density: &Array1, @@ -77,7 +77,7 @@ pub trait HardSphereProperties { } /// The fraction $\frac{\zeta_2}{\zeta_3}$ evaluated in a way to avoid a division by 0 when the density is 0. - fn zeta_23>(&self, temperature: D, molefracs: &Array1) -> D { + fn zeta_23 + Copy>(&self, temperature: D, molefracs: &Array1) -> D { let component_index = self.component_index(); let geometry_coefficients = self.geometry_coefficients(temperature); let diameter = self.hs_diameter(temperature); @@ -116,7 +116,7 @@ impl

HardSphere

{ } } -impl, P: HardSphereProperties> HelmholtzEnergyDual for HardSphere

{ +impl + Copy, P: HardSphereProperties> HelmholtzEnergyDual for HardSphere

{ fn helmholtz_energy(&self, state: &StateHD) -> D { let p = &self.parameters; let zeta = p.zeta(state.temperature, &state.partial_density, [0, 1, 2, 3]); diff --git a/src/pcsaft/dft/dispersion.rs b/src/pcsaft/dft/dispersion.rs index 05908daac..f3d7c7181 100644 --- a/src/pcsaft/dft/dispersion.rs +++ b/src/pcsaft/dft/dispersion.rs @@ -28,7 +28,7 @@ impl AttractiveFunctional { } } -fn att_weight_functions + ScalarOperand>( +fn att_weight_functions + Copy + ScalarOperand>( p: &PcSaftParameters, psi: f64, temperature: N, @@ -40,7 +40,9 @@ fn att_weight_functions + ScalarOperand>( ) } -impl + ScalarOperand> FunctionalContributionDual for AttractiveFunctional { +impl + Copy + ScalarOperand> FunctionalContributionDual + for AttractiveFunctional +{ fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { att_weight_functions(&self.parameters, PSI_DFT, temperature) } diff --git a/src/pcsaft/dft/hard_chain.rs b/src/pcsaft/dft/hard_chain.rs index df8d2013a..7119e2e40 100644 --- a/src/pcsaft/dft/hard_chain.rs +++ b/src/pcsaft/dft/hard_chain.rs @@ -20,7 +20,7 @@ impl ChainFunctional { } } -impl + ScalarOperand> FunctionalContributionDual for ChainFunctional { +impl + Copy + ScalarOperand> FunctionalContributionDual for ChainFunctional { fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { let p = &self.parameters; let d = p.hs_diameter(temperature); diff --git a/src/pcsaft/dft/polar.rs b/src/pcsaft/dft/polar.rs index 606e77c4d..efb20adf7 100644 --- a/src/pcsaft/dft/polar.rs +++ b/src/pcsaft/dft/polar.rs @@ -8,7 +8,7 @@ use ndarray::*; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_3, PI}; -pub(super) fn calculate_helmholtz_energy_density_polar + ScalarOperand>( +pub(super) fn calculate_helmholtz_energy_density_polar + Copy + ScalarOperand>( parameters: &PcSaftParameters, temperature: N, density: ArrayView2, @@ -38,7 +38,7 @@ pub(super) fn calculate_helmholtz_energy_density_polar + ScalarO Ok(phi) } -pub fn pair_integral_ij + ScalarOperand>( +pub fn pair_integral_ij + Copy + ScalarOperand>( mij1: f64, mij2: f64, eta: &Array1, @@ -93,7 +93,7 @@ fn triplet_integral_ijk_dq + ScalarOperand>( integral } -fn phi_polar_dipole + ScalarOperand>( +fn phi_polar_dipole + Copy + ScalarOperand>( p: &PcSaftParameters, temperature: N, density: ArrayView2, @@ -181,7 +181,7 @@ fn phi_polar_dipole + ScalarOperand>( Ok(result) } -fn phi_polar_quadrupole + ScalarOperand>( +fn phi_polar_quadrupole + Copy + ScalarOperand>( p: &PcSaftParameters, temperature: N, density: ArrayView2, @@ -269,7 +269,7 @@ fn phi_polar_quadrupole + ScalarOperand>( Ok(result) } -fn phi_polar_dipole_quadrupole + ScalarOperand>( +fn phi_polar_dipole_quadrupole + Copy + ScalarOperand>( p: &PcSaftParameters, temperature: N, density: ArrayView2, diff --git a/src/pcsaft/dft/pure_saft_functional.rs b/src/pcsaft/dft/pure_saft_functional.rs index 801a079a8..d8b0291a1 100644 --- a/src/pcsaft/dft/pure_saft_functional.rs +++ b/src/pcsaft/dft/pure_saft_functional.rs @@ -35,7 +35,9 @@ impl PureFMTAssocFunctional { } } -impl + ScalarOperand> FunctionalContributionDual for PureFMTAssocFunctional { +impl + Copy + ScalarOperand> FunctionalContributionDual + for PureFMTAssocFunctional +{ fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { let r = self.parameters.hs_diameter(temperature) * 0.5; WeightFunctionInfo::new(arr1(&[0]), false).extend( @@ -154,7 +156,7 @@ impl PureChainFunctional { } } -impl + ScalarOperand> FunctionalContributionDual for PureChainFunctional { +impl + Copy + ScalarOperand> FunctionalContributionDual for PureChainFunctional { fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { let d = self.parameters.hs_diameter(temperature); WeightFunctionInfo::new(arr1(&[0]), true) @@ -206,7 +208,7 @@ impl PureAttFunctional { } } -impl + ScalarOperand> FunctionalContributionDual for PureAttFunctional { +impl + Copy + ScalarOperand> FunctionalContributionDual for PureAttFunctional { fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { let d = self.parameters.hs_diameter(temperature); const PSI: f64 = 1.3862; // Homosegmented DFT (Sauer2017) diff --git a/src/pcsaft/eos/dispersion.rs b/src/pcsaft/eos/dispersion.rs index a8a4f0c9d..416778685 100644 --- a/src/pcsaft/eos/dispersion.rs +++ b/src/pcsaft/eos/dispersion.rs @@ -65,7 +65,7 @@ pub struct Dispersion { pub parameters: Arc, } -impl> HelmholtzEnergyDual for Dispersion { +impl + Copy> HelmholtzEnergyDual for Dispersion { fn helmholtz_energy(&self, state: &StateHD) -> D { // auxiliary variables let n = self.parameters.m.len(); diff --git a/src/pcsaft/eos/hard_chain.rs b/src/pcsaft/eos/hard_chain.rs index ea604e26d..a046c893c 100644 --- a/src/pcsaft/eos/hard_chain.rs +++ b/src/pcsaft/eos/hard_chain.rs @@ -10,7 +10,7 @@ pub struct HardChain { pub parameters: Arc, } -impl> HelmholtzEnergyDual for HardChain { +impl + Copy> HelmholtzEnergyDual for HardChain { fn helmholtz_energy(&self, state: &StateHD) -> D { let p = &self.parameters; let d = self.parameters.hs_diameter(state.temperature); diff --git a/src/pcsaft/eos/polar.rs b/src/pcsaft/eos/polar.rs index b94df9087..9efe1c428 100644 --- a/src/pcsaft/eos/polar.rs +++ b/src/pcsaft/eos/polar.rs @@ -127,7 +127,7 @@ impl MeanSegmentNumbers { } } -fn pair_integral_ij>( +fn pair_integral_ij + Copy>( mij1: f64, mij2: f64, etas: &[D], @@ -144,13 +144,18 @@ fn pair_integral_ij>( .sum() } -fn triplet_integral_ijk>(mijk1: f64, mijk2: f64, etas: &[D], c: &[[f64; 3]]) -> D { +fn triplet_integral_ijk + Copy>( + mijk1: f64, + mijk2: f64, + etas: &[D], + c: &[[f64; 3]], +) -> D { (0..c.len()) .map(|i| etas[i] * (c[i][0] + mijk1 * c[i][1] + mijk2 * c[i][2])) .sum() } -fn triplet_integral_ijk_dq>(mijk: f64, etas: &[D], c: &[[f64; 2]]) -> D { +fn triplet_integral_ijk_dq + Copy>(mijk: f64, etas: &[D], c: &[[f64; 2]]) -> D { (0..c.len()) .map(|i| etas[i] * (c[i][0] + mijk * c[i][1])) .sum() @@ -160,7 +165,7 @@ pub struct Dipole { pub parameters: Arc, } -impl> HelmholtzEnergyDual for Dipole { +impl + Copy> HelmholtzEnergyDual for Dipole { fn helmholtz_energy(&self, state: &StateHD) -> D { let m = MeanSegmentNumbers::new(&self.parameters, Multipole::Dipole); let p = &self.parameters; @@ -237,7 +242,7 @@ pub struct Quadrupole { pub parameters: Arc, } -impl> HelmholtzEnergyDual for Quadrupole { +impl + Copy> HelmholtzEnergyDual for Quadrupole { fn helmholtz_energy(&self, state: &StateHD) -> D { let m = MeanSegmentNumbers::new(&self.parameters, Multipole::Quadrupole); let p = &self.parameters; @@ -324,7 +329,7 @@ pub struct DipoleQuadrupole { pub variant: DQVariants, } -impl> HelmholtzEnergyDual for DipoleQuadrupole { +impl + Copy> HelmholtzEnergyDual for DipoleQuadrupole { fn helmholtz_energy(&self, state: &StateHD) -> D { let p = &self.parameters; diff --git a/src/pcsaft/eos/qspr.rs b/src/pcsaft/eos/qspr.rs index 0f2ce9b4b..844c21ce5 100644 --- a/src/pcsaft/eos/qspr.rs +++ b/src/pcsaft/eos/qspr.rs @@ -68,7 +68,7 @@ pub struct QSPR { pub parameters: Arc, } -impl> IdealGasContributionDual for QSPR { +impl + Copy> IdealGasContributionDual for QSPR { fn de_broglie_wavelength(&self, temperature: D, components: usize) -> Array1 { let (c_300, c_400) = if self.parameters.association.is_empty() { match self.parameters.ndipole + self.parameters.nquadpole { diff --git a/src/pcsaft/parameters.rs b/src/pcsaft/parameters.rs index 5cfef1564..892cbee6b 100644 --- a/src/pcsaft/parameters.rs +++ b/src/pcsaft/parameters.rs @@ -469,7 +469,7 @@ impl HardSphereProperties for PcSaftParameters { MonomerShape::NonSpherical(self.m.mapv(N::from)) } - fn hs_diameter>(&self, temperature: D) -> Array1 { + fn hs_diameter + Copy>(&self, temperature: D) -> Array1 { let ti = temperature.recip() * -3.0; Array::from_shape_fn(self.sigma.len(), |i| { -((ti * self.epsilon_k[i]).exp() * 0.12 - 1.0) * self.sigma[i] diff --git a/src/pets/dft/dispersion.rs b/src/pets/dft/dispersion.rs index 79c648497..89e556aee 100644 --- a/src/pets/dft/dispersion.rs +++ b/src/pets/dft/dispersion.rs @@ -27,7 +27,7 @@ impl AttractiveFunctional { } } -fn att_weight_functions + ScalarOperand>( +fn att_weight_functions + Copy + ScalarOperand>( p: &PetsParameters, psi: f64, temperature: N, @@ -39,7 +39,9 @@ fn att_weight_functions + ScalarOperand>( ) } -impl + ScalarOperand> FunctionalContributionDual for AttractiveFunctional { +impl + Copy + ScalarOperand> FunctionalContributionDual + for AttractiveFunctional +{ fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { att_weight_functions(&self.parameters, PSI_DFT, temperature) } diff --git a/src/pets/dft/pure_pets_functional.rs b/src/pets/dft/pure_pets_functional.rs index b7aaaa798..85ec7c910 100644 --- a/src/pets/dft/pure_pets_functional.rs +++ b/src/pets/dft/pure_pets_functional.rs @@ -29,7 +29,7 @@ impl PureFMTFunctional { } } -impl + ScalarOperand> FunctionalContributionDual for PureFMTFunctional { +impl + Copy + ScalarOperand> FunctionalContributionDual for PureFMTFunctional { fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { let r = self.parameters.hs_diameter(temperature) * 0.5; WeightFunctionInfo::new(arr1(&[0]), false).extend( @@ -125,7 +125,7 @@ impl PureAttFunctional { } } -impl + ScalarOperand> FunctionalContributionDual for PureAttFunctional { +impl + Copy + ScalarOperand> FunctionalContributionDual for PureAttFunctional { fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { let d = self.parameters.hs_diameter(temperature); const PSI: f64 = 1.21; // Homosegmented DFT (Heier2018) diff --git a/src/pets/eos/dispersion.rs b/src/pets/eos/dispersion.rs index 09db1c17c..4b656d412 100644 --- a/src/pets/eos/dispersion.rs +++ b/src/pets/eos/dispersion.rs @@ -30,7 +30,7 @@ pub struct Dispersion { pub parameters: Arc, } -impl> HelmholtzEnergyDual for Dispersion { +impl + Copy> HelmholtzEnergyDual for Dispersion { fn helmholtz_energy(&self, state: &StateHD) -> D { // auxiliary variables let n = self.parameters.sigma.len(); diff --git a/src/pets/eos/qspr.rs b/src/pets/eos/qspr.rs index 78f9e8572..e4d5aa6e4 100644 --- a/src/pets/eos/qspr.rs +++ b/src/pets/eos/qspr.rs @@ -68,7 +68,7 @@ pub struct QSPR { pub parameters: Arc, } -impl> IdealGasContributionDual for QSPR { +impl + Copy> IdealGasContributionDual for QSPR { fn de_broglie_wavelength(&self, temperature: D, components: usize) -> Array1 { let (c_300, c_400) = (NA_NP_300, NA_NP_400); diff --git a/src/pets/parameters.rs b/src/pets/parameters.rs index 096f6e724..436c9dc3a 100644 --- a/src/pets/parameters.rs +++ b/src/pets/parameters.rs @@ -237,7 +237,7 @@ impl HardSphereProperties for PetsParameters { MonomerShape::Spherical(self.sigma.len()) } - fn hs_diameter>(&self, temperature: D) -> Array1 { + fn hs_diameter + Copy>(&self, temperature: D) -> Array1 { let ti = temperature.recip() * -3.052785558; Array::from_shape_fn(self.sigma.len(), |i| { -((ti * self.epsilon_k[i]).exp() * 0.127112544 - 1.0) * self.sigma[i] diff --git a/src/saftvrqmie/dft/dispersion.rs b/src/saftvrqmie/dft/dispersion.rs index c824385bf..4b21d0538 100644 --- a/src/saftvrqmie/dft/dispersion.rs +++ b/src/saftvrqmie/dft/dispersion.rs @@ -26,7 +26,7 @@ impl AttractiveFunctional { } } -fn att_weight_functions + ScalarOperand>( +fn att_weight_functions + Copy + ScalarOperand>( p: &SaftVRQMieParameters, psi: f64, temperature: N, @@ -38,7 +38,9 @@ fn att_weight_functions + ScalarOperand>( ) } -impl + ScalarOperand> FunctionalContributionDual for AttractiveFunctional { +impl + Copy + ScalarOperand> FunctionalContributionDual + for AttractiveFunctional +{ fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { att_weight_functions(&self.parameters, PSI_DFT, temperature) } diff --git a/src/saftvrqmie/dft/mod.rs b/src/saftvrqmie/dft/mod.rs index da9d32855..d728fe9db 100644 --- a/src/saftvrqmie/dft/mod.rs +++ b/src/saftvrqmie/dft/mod.rs @@ -116,7 +116,7 @@ impl HardSphereProperties for SaftVRQMieParameters { MonomerShape::Spherical(self.m.len()) } - fn hs_diameter>(&self, temperature: D) -> Array1 { + fn hs_diameter + Copy>(&self, temperature: D) -> Array1 { self.hs_diameter(temperature) } } diff --git a/src/saftvrqmie/dft/non_additive_hs.rs b/src/saftvrqmie/dft/non_additive_hs.rs index fe94e362e..f7517c700 100644 --- a/src/saftvrqmie/dft/non_additive_hs.rs +++ b/src/saftvrqmie/dft/non_additive_hs.rs @@ -24,7 +24,7 @@ impl NonAddHardSphereFunctional { impl FunctionalContributionDual for NonAddHardSphereFunctional where - N: DualNum + ScalarOperand, + N: DualNum + Copy + ScalarOperand, { fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { let p = &self.parameters; @@ -146,7 +146,7 @@ where } } -pub fn non_additive_hs_energy_density + ScalarOperand>( +pub fn non_additive_hs_energy_density + Copy + ScalarOperand>( parameters: &SaftVRQMieParameters, d_hs_ij: &Array2, d_hs_add_ij: &Array2, diff --git a/src/saftvrqmie/eos/dispersion.rs b/src/saftvrqmie/eos/dispersion.rs index 1db59440e..8bff9d590 100644 --- a/src/saftvrqmie/eos/dispersion.rs +++ b/src/saftvrqmie/eos/dispersion.rs @@ -34,7 +34,7 @@ pub struct Alpha> { alpha_ij: Array2, } -impl> Alpha { +impl + Copy> Alpha { pub fn new( parameters: &SaftVRQMieParameters, sigma_eff_ij: &Array2, @@ -76,7 +76,7 @@ pub struct Dispersion { pub parameters: Arc, } -impl> HelmholtzEnergyDual for Dispersion { +impl + Copy> HelmholtzEnergyDual for Dispersion { fn helmholtz_energy(&self, state: &StateHD) -> D { // auxiliary variables let n = self.parameters.m.len(); @@ -136,7 +136,7 @@ impl> HelmholtzEnergyDual for Dispersion { } #[cfg(feature = "dft")] -pub fn dispersion_energy_density>( +pub fn dispersion_energy_density + Copy>( parameters: &SaftVRQMieParameters, d_hs_ij: &Array2, s_eff_ij: &Array2, @@ -172,7 +172,7 @@ pub fn dispersion_energy_density>( rho_s * (a1 * inv_t + a2 * inv_t.powi(2) + a3 * inv_t.powi(3)) } -fn zeta_saft_vrq_mie>( +fn zeta_saft_vrq_mie + Copy>( m: &Array1, x_s: &Array1, diameter: &Array2, @@ -187,7 +187,7 @@ fn zeta_saft_vrq_mie>( zeta * FRAC_PI_6 * rho_s } -fn first_order_perturbation>( +fn first_order_perturbation + Copy>( parameters: &SaftVRQMieParameters, x_s: &Array1, zeta: D, @@ -223,7 +223,7 @@ fn first_order_perturbation>( a1 } -fn first_order_perturbation_ij>( +fn first_order_perturbation_ij + Copy>( lambda_a: f64, lambda_r: f64, epsilon_k: f64, @@ -244,7 +244,7 @@ fn first_order_perturbation_ij>( (int_qa * qa1 - int_qr * qr1 + int_a - int_r) * c } -fn eta_eff>(lambda: f64, zeta: D) -> D { +fn eta_eff + Copy>(lambda: f64, zeta: D) -> D { let inv_lambda = Array1::from(vec![ 1.0, 1.0 / lambda, @@ -260,7 +260,7 @@ fn eta_eff>(lambda: f64, zeta: D) -> D { zeta * (zeta * (zeta * (zeta * c[3] + c[2]) + c[1]) + c[0]) } -fn sutherland>(lambda: f64, epsilon_k: f64, zeta: D, x0: D) -> D { +fn sutherland + Copy>(lambda: f64, epsilon_k: f64, zeta: D, x0: D) -> D { let ef = eta_eff(lambda, zeta); (-ef * 0.5 + 1.0) * -12.0 * x0.powf(lambda) * epsilon_k / (lambda - 3.0) / (-ef + 1.0).powi(3) } @@ -278,7 +278,7 @@ fn jlambda>(lambda: f64, x0: D) -> D { /// B is divided by the packing fraction /// /// \author Morten Hammer, February 2018 -fn b>(lambda: f64, epsilon_k: f64, zeta: D, x0: D, x0_eff: D) -> D { +fn b + Copy>(lambda: f64, epsilon_k: f64, zeta: D, x0: D, x0_eff: D) -> D { let ilambda = ilambda(lambda, x0_eff); let jlambda = jlambda(lambda, x0_eff); let denum = (-zeta + 1.0).powi(3); @@ -289,7 +289,7 @@ fn b>(lambda: f64, epsilon_k: f64, zeta: D, x0: D, x0_eff: D) -> } #[inline] -fn combine_sutherland_and_b>( +fn combine_sutherland_and_b + Copy>( lambda: f64, epsilon_k: f64, zeta: D, @@ -301,7 +301,7 @@ fn combine_sutherland_and_b>( int_as + int_b } -fn second_order_perturbation>( +fn second_order_perturbation + Copy>( parameters: &SaftVRQMieParameters, alpha: &Alpha, x_s: &Array1, @@ -353,7 +353,7 @@ fn quantum_prefactor(lambda: f64) -> f64 { lambda * (lambda - 1.0) } -fn second_order_perturbation_ij>( +fn second_order_perturbation_ij + Copy>( lambda_a: f64, lambda_r: f64, epsilon_k: f64, @@ -398,7 +398,7 @@ fn second_order_perturbation_ij>( a2_ij * 0.5 * epsilon_k * c.powi(2) } -fn third_order_perturbation>( +fn third_order_perturbation + Copy>( parameters: &SaftVRQMieParameters, alpha: &Alpha, x_s: &Array1, @@ -417,7 +417,7 @@ fn third_order_perturbation>( a3 } -fn third_order_perturbation_ij>( +fn third_order_perturbation_ij + Copy>( i: usize, j: usize, epsilon_k_eff: D, @@ -444,9 +444,6 @@ mod tests { use crate::saftvrqmie::parameters::utils::hydrogen_fh1; use approx::assert_relative_eq; use ndarray::arr1; - use num_dual::Dual2; - use num_dual::DualNum; - use num_traits::Zero; #[test] fn test_eta_eff() { @@ -466,7 +463,7 @@ mod tests { #[test] fn test_alpha() { - let temperature = Dual2::from_re(26.7060).derive(); + let temperature = 26.7060; let parameters = hydrogen_fh1(); let n = 1; let s_eff_ij = Array2::from_shape_fn((n, n), |(i, j)| { @@ -476,46 +473,42 @@ mod tests { parameters.calc_epsilon_k_eff_ij(i, j, temperature) }); let alpha = Alpha::new(¶meters, &s_eff_ij, &epsilon_k_eff_ij, temperature); - assert_relative_eq!( - alpha.alpha_ij[[0, 0]].re(), - 1.0239374984636636, - epsilon = 5e-8 - ); + assert_relative_eq!(alpha.alpha_ij[[0, 0]], 1.0239374984636636, epsilon = 5e-8); } #[test] fn test_sutherland() { - let x0 = Dual2::from_re(1.1).derive(); - let zeta = Dual2::from_re(0.333).derive(); + let x0 = 1.1; + let zeta = 0.333; let lambda = 13.77; let eps_div_k = 13.88; let asa = sutherland(lambda, eps_div_k, zeta, x0); - assert_relative_eq!(asa.re(), -122.12017536923423, epsilon = 1e-12); + assert_relative_eq!(asa, -122.12017536923423, epsilon = 1e-12); } #[test] fn test_b() { - let x0 = Dual2::from_re(1.1).derive(); - let zeta = Dual2::from_re(0.333).derive(); + let x0 = 1.1; + let zeta = 0.333; let lambda = 13.77; let eps_div_k = 13.88; let ba = b(lambda, eps_div_k, zeta, x0, x0); - assert_relative_eq!(ba.re(), 93.436438943866293, epsilon = 1e-12); + assert_relative_eq!(ba, 93.436438943866293, epsilon = 1e-12); } #[test] fn test_quantum_d_ij() { let p = hydrogen_fh1(); - let temperature = Dual2::from_re(26.7060).derive(); + let temperature = 26.7060; let dq_ij = p.quantum_d_ij(0, 0, temperature); - assert_relative_eq!(dq_ij.re(), 7.5092605940987542e-2, epsilon = 5e-8); + assert_relative_eq!(dq_ij, 7.5092605940987542e-2, epsilon = 5e-8); } #[test] fn test_first_order_perturbation_ij() { let p = hydrogen_fh1(); - let temperature = Dual2::from_re(26.7060).derive(); - let zeta = Dual2::from_re(0.333).derive(); + let temperature = 26.7060; + let zeta = 0.333; let dq_div_s2 = p.quantum_d_ij(0, 0, temperature) / p.sigma_ij[[0, 0]].powi(2); let s_eff = p.calc_sigma_eff_ij(0, 0, temperature); let d_hs = p.hs_diameter_ij(0, 0, temperature, s_eff); @@ -532,15 +525,15 @@ mod tests { p.c_ij[[0, 0]], dq_div_s2, ); - let rel_err = (a1_ij.re() + 332.00915966785539) / 332.00915966785539; + let rel_err = (a1_ij + 332.00915966785539) / 332.00915966785539; assert_relative_eq!(rel_err, 0.0, epsilon = 1e-7); } #[test] fn test_second_order_perturbation_ij() { let p = hydrogen_fh1(); - let temperature = Dual2::from_re(26.7060).derive(); - let zeta = Dual2::from_re(0.333).derive(); + let temperature = 26.7060; + let zeta = 0.333; let dq_div_s2 = p.quantum_d_ij(0, 0, temperature) / p.sigma_ij[[0, 0]].powi(2); let s_eff = p.calc_sigma_eff_ij(0, 0, temperature); let d_hs = p.hs_diameter_ij(0, 0, temperature, s_eff); @@ -557,15 +550,15 @@ mod tests { p.c_ij[[0, 0]], dq_div_s2, ); - let rel_err = (a2_ij.re() + 1907.5055256805874) / 1907.5055256805874; + let rel_err = (a2_ij + 1907.5055256805874) / 1907.5055256805874; assert_relative_eq!(rel_err, 0.0, epsilon = 1e-7); } #[test] fn test_third_order_perturbation_ij() { let p = hydrogen_fh1(); - let temperature = Dual2::from_re(26.7060).derive(); - let zeta_bar = Dual2::from_re(0.333).derive(); + let temperature = 26.7060; + let zeta_bar = 0.333; let n = 1; let s_eff_ij = Array2::from_shape_fn((n, n), |(i, j)| p.calc_sigma_eff_ij(i, j, temperature)); @@ -575,16 +568,16 @@ mod tests { let a3_ij = third_order_perturbation_ij(0, 0, epsilon_k_eff_ij[[0, 0]], &alpha, zeta_bar); - let rel_err = (a3_ij.re() + 25.807966819127916) / 25.807966819127916; + let rel_err = (a3_ij + 25.807966819127916) / 25.807966819127916; assert_relative_eq!(rel_err, 0.0, epsilon = 5e-7); } #[test] fn test_zeta_saft_vrq_mie() { let p = hydrogen_fh1(); - let t = Dual2::from_re(26.7060).derive(); - let v = Dual2::from_re(1.0e26).derive(); - let n = Dual2::from_re(6.02214076e23).derive(); + let t = 26.7060; + let v = 1.0e26; + let n = 6.02214076e23; let state = StateHD::new(t, v, arr1(&[n])); let nc = 1; // temperature dependent sigma @@ -603,23 +596,23 @@ mod tests { x_s[i] *= inv_x_s_sum; } // Segment density - let mut rho_s = Dual2::zero(); + let mut rho_s = 0.0; for i in 0..nc { rho_s += state.partial_density[i] * p.m[i]; } // packing fractions let zeta = zeta_saft_vrq_mie(&p.m, &x_s, &d_hs_ij, rho_s); let zeta_bar = zeta_saft_vrq_mie(&p.m, &x_s, &s_eff_ij, rho_s); - assert_relative_eq!(zeta.re(), 9.7717457994590765E-002, epsilon = 5e-9); - assert_relative_eq!(zeta_bar.re(), 0.10864364645845238, epsilon = 5e-9); + assert_relative_eq!(zeta, 9.7717457994590765E-002, epsilon = 5e-9); + assert_relative_eq!(zeta_bar, 0.10864364645845238, epsilon = 5e-9); } #[test] fn test_perturbation_terms() { let p = hydrogen_fh1(); - let t = Dual2::from_re(26.7060).derive(); - let v = Dual2::from_re(1.0e26).derive(); - let n = Dual2::from_re(6.02214076e23).derive(); + let t = 26.7060; + let v = 1.0e26; + let n = 6.02214076e23; let state = StateHD::new(t, v, arr1(&[n])); let nc = 1; // temperature dependent sigma @@ -639,7 +632,7 @@ mod tests { } // Segment density - let mut rho_s = Dual2::zero(); + let mut rho_s = 0.0; for i in 0..nc { rho_s += state.partial_density[i] * p.m[i]; } @@ -666,12 +659,12 @@ mod tests { ); let a3 = third_order_perturbation(&p, &alpha, &x_s, zeta_bar, &epsilon_k_eff_ij); - let rel_err_a1 = (a1.re() + 30.702499892515764) / 30.702499892515764; - let rel_err_a2 = (a2.re() + 67.046957636607587) / 67.046957636607587; - let rel_err_a3 = (a3.re() + 470.96241656623727) / 470.96241656623727; - assert_relative_eq!(rel_err_a1.re(), 0.0, epsilon = 5e-7); - assert_relative_eq!(rel_err_a2.re(), 0.0, epsilon = 5e-7); - assert_relative_eq!(rel_err_a3.re(), 0.0, epsilon = 5e-7); + let rel_err_a1 = (a1 + 30.702499892515764) / 30.702499892515764; + let rel_err_a2 = (a2 + 67.046957636607587) / 67.046957636607587; + let rel_err_a3 = (a3 + 470.96241656623727) / 470.96241656623727; + assert_relative_eq!(rel_err_a1, 0.0, epsilon = 5e-7); + assert_relative_eq!(rel_err_a2, 0.0, epsilon = 5e-7); + assert_relative_eq!(rel_err_a3, 0.0, epsilon = 5e-7); } #[test] @@ -688,19 +681,18 @@ mod tests { ]; let na = 6.02214076e23; for (it, &a) in a_ref.iter().enumerate() { - let t = Dual2::from_re(26.7060 * (it + 1) as f64).derive(); - let v = Dual2::from_re(1.0e26).derive(); - let n = Dual2::from_re(na).derive(); - let state = StateHD::new(t, v, arr1(&[n])); + let t = 26.7060 * (it + 1) as f64; + let v = 1.0e26; + let state = StateHD::new(t, v, arr1(&[na])); let a_disp = disp.helmholtz_energy(&state) / na; - assert_relative_eq!(a_disp.re(), a, epsilon = 1e-7); + assert_relative_eq!(a_disp, a, epsilon = 1e-7); } - let t = Dual2::from_re(26.7060).derive(); - let v = Dual2::from_re(1.0e26 * 2.0).derive(); - let n = Dual2::from_re(na * 2.0).derive(); + let t = 26.7060; + let v = 1.0e26 * 2.0; + let n = na * 2.0; let state = StateHD::new(t, v, arr1(&[n])); let a_disp = disp.helmholtz_energy(&state) / na; - assert_relative_eq!(a_disp.re(), a_ref[0] * 2.0, epsilon = 1e-7); + assert_relative_eq!(a_disp, a_ref[0] * 2.0, epsilon = 1e-7); } #[test] @@ -728,21 +720,21 @@ mod tests { -0.84210863940206726, ]; let na = 6.02214076e23; - let n = [Dual2::from_re(1.1 * na), Dual2::from_re(1.0 * na)]; - let v = Dual2::from_re(1.0e26).derive(); + let n = [1.1 * na, 1.0 * na]; + let v = 1.0e26; for (it, &a) in a_ref.iter().enumerate() { - let t = Dual2::from_re(30.0 * (it + 1) as f64).derive(); + let t = 30.0 * (it + 1) as f64; let state = StateHD::new(t, v, arr1(&n)); let a_disp = disp.helmholtz_energy(&state) / na; dbg!(it); - assert_relative_eq!(a_disp.re(), a, epsilon = 1e-7); + assert_relative_eq!(a_disp, a, epsilon = 1e-7); } - let t = Dual2::from_re(30.0).derive(); - let v = Dual2::from_re(1.0e26 * 2.0).derive(); - let n = [Dual2::from_re(2.2 * na), Dual2::from_re(2.0 * na)]; + let t = 30.0; + let v = 1.0e26 * 2.0; + let n = [2.2 * na, 2.0 * na]; let state = StateHD::new(t, v, arr1(&n)); let a_disp = disp.helmholtz_energy(&state) / na; - assert_relative_eq!(a_disp.re(), a_ref[0] * 2.0, epsilon = 1e-7); + assert_relative_eq!(a_disp, a_ref[0] * 2.0, epsilon = 1e-7); } #[cfg(feature = "dft")] @@ -753,8 +745,8 @@ mod tests { }; let p = &disp.parameters; let n = p.m.len(); - let rho = Array1::from_shape_fn(n, |_i| Dual2::from_re(0.01)); - let t = Dual2::from_re(25.0).derive(); + let rho = Array1::from_shape_fn(n, |_i| 0.01); + let t = 25.0; // temperature dependent segment radius // calc & store this in struct let s_eff_ij = Array2::from_shape_fn((n, n), |(i, j)| p.calc_sigma_eff_ij(i, j, t)); @@ -783,7 +775,7 @@ mod tests { ); dbg!(rho); - dbg!(a_disp.re()); - assert_relative_eq!(a_disp.re(), -0.022349175545184223, epsilon = 1e-7); + dbg!(a_disp); + assert_relative_eq!(a_disp, -0.022349175545184223, epsilon = 1e-7); } } diff --git a/src/saftvrqmie/eos/hard_sphere.rs b/src/saftvrqmie/eos/hard_sphere.rs index 452339a35..67e72b713 100644 --- a/src/saftvrqmie/eos/hard_sphere.rs +++ b/src/saftvrqmie/eos/hard_sphere.rs @@ -62,7 +62,7 @@ const W_K21: [f64; 21] = [ impl SaftVRQMieParameters { #[inline] - pub fn hs_diameter>(&self, temperature: D) -> Array1 { + pub fn hs_diameter + Copy>(&self, temperature: D) -> Array1 { Array1::from_shape_fn(self.m.len(), |i| -> D { let sigma_eff = self.calc_sigma_eff_ij(i, i, temperature); self.hs_diameter_ij(i, i, temperature, sigma_eff) @@ -70,7 +70,7 @@ impl SaftVRQMieParameters { } #[inline] - pub fn hs_diameter_ij>( + pub fn hs_diameter_ij + Copy>( &self, i: usize, j: usize, @@ -89,7 +89,7 @@ impl SaftVRQMieParameters { d_hs } - pub fn zero_integrand>( + pub fn zero_integrand + Copy>( &self, i: usize, j: usize, @@ -118,13 +118,18 @@ impl SaftVRQMieParameters { } #[inline] - pub fn epsilon_k_eff>(&self, temperature: D) -> Array1 { + pub fn epsilon_k_eff + Copy>(&self, temperature: D) -> Array1 { Array1::from_shape_fn(self.m.len(), |i| -> D { self.calc_epsilon_k_eff_ij(i, i, temperature) }) } - pub fn calc_epsilon_k_eff_ij>(&self, i: usize, j: usize, temperature: D) -> D { + pub fn calc_epsilon_k_eff_ij + Copy>( + &self, + i: usize, + j: usize, + temperature: D, + ) -> D { let mut r = D::one() * self.sigma_ij[[i, j]]; let mut u_vec = [D::zero(), D::zero(), D::zero()]; for _k in 1..20 { @@ -141,13 +146,18 @@ impl SaftVRQMieParameters { } #[inline] - pub fn sigma_eff>(&self, temperature: D) -> Array1 { + pub fn sigma_eff + Copy>(&self, temperature: D) -> Array1 { Array1::from_shape_fn(self.m.len(), |i| -> D { self.calc_sigma_eff_ij(i, i, temperature) }) } - pub fn calc_sigma_eff_ij>(&self, i: usize, j: usize, temperature: D) -> D { + pub fn calc_sigma_eff_ij + Copy>( + &self, + i: usize, + j: usize, + temperature: D, + ) -> D { let mut r = D::one() * self.sigma_ij[[i, j]]; let mut u_vec = [D::zero(), D::zero(), D::zero()]; for _k in 1..20 { @@ -169,7 +179,7 @@ impl SaftVRQMieParameters { } /// Feynman-Hibbs corrected potential - pub fn qmie_potential_ij>( + pub fn qmie_potential_ij + Copy>( &self, i: usize, j: usize, @@ -219,7 +229,7 @@ pub struct HardSphere { pub parameters: Arc, } -impl> HelmholtzEnergyDual for HardSphere { +impl + Copy> HelmholtzEnergyDual for HardSphere { fn helmholtz_energy(&self, state: &StateHD) -> D { let d = self.parameters.hs_diameter(state.temperature); let zeta = zeta(&self.parameters.m, &state.partial_density, &d); @@ -238,7 +248,7 @@ impl fmt::Display for HardSphere { } } -pub fn zeta>( +pub fn zeta + Copy>( m: &Array1, partial_density: &Array1, diameter: &Array1, @@ -254,7 +264,11 @@ pub fn zeta>( zeta } -pub fn zeta_23>(m: &Array1, molefracs: &Array1, diameter: &Array1) -> D { +pub fn zeta_23 + Copy>( + m: &Array1, + molefracs: &Array1, + diameter: &Array1, +) -> D { let mut zeta: [D; 2] = [D::zero(), D::zero()]; for i in 0..m.len() { for (k, z) in zeta.iter_mut().enumerate() { @@ -271,59 +285,59 @@ mod tests { use crate::saftvrqmie::parameters::utils::hydrogen_fh1; use approx::assert_relative_eq; use ndarray::arr1; - use num_dual::Dual2; + use num_dual::Dual64; #[test] fn test_quantum_d_mass() { let parameters = hydrogen_fh1(); - let temperature = Dual2::from_re(26.7060).derive(); - let r = Dual2::from_re(3.5); + let temperature = 26.7060; + let r = 3.5; let u0 = parameters.qmie_potential_ij(0, 0, r, temperature); let eps = 1.0e-5; let u2 = parameters.qmie_potential_ij(0, 0, r + eps, temperature); let u1 = parameters.qmie_potential_ij(0, 0, r - eps, temperature); - let dudr_num = (u2[0].re() - u1[0].re()) / eps / 2.0; - let d2udr2_num = (u2[1].re() - u1[1].re()) / eps / 2.0; - assert!(((dudr_num - u0[1].re()) / u0[1].re()).abs() < 1.0e-9); - assert!(((d2udr2_num - u0[2].re()) / u0[2].re()).abs() < 1.0e-9); + let dudr_num = (u2[0] - u1[0]) / eps / 2.0; + let d2udr2_num = (u2[1] - u1[1]) / eps / 2.0; + assert!(((dudr_num - u0[1]) / u0[1]).abs() < 1.0e-9); + assert!(((d2udr2_num - u0[2]) / u0[2]).abs() < 1.0e-9); } #[test] fn test_sigma_effective() { let parameters = hydrogen_fh1(); - let temperature = Dual2::from_re(26.7060).derive(); + let temperature = 26.7060; let sigma_eff = parameters.calc_sigma_eff_ij(0, 0, temperature); - println!("{}", sigma_eff.re() - 3.2540054024660556); - assert!((sigma_eff.re() - 3.2540054024660556).abs() < 5.0e-7) + println!("{}", sigma_eff - 3.2540054024660556); + assert!((sigma_eff - 3.2540054024660556).abs() < 5.0e-7) } #[test] fn test_eps_div_k_effective() { let parameters = hydrogen_fh1(); - let temperature = Dual2::from_re(26.7060).derive(); + let temperature = 26.7060; let epsilon_k_eff = parameters.calc_epsilon_k_eff_ij(0, 0, temperature); - println!("{}", epsilon_k_eff.re() - 21.654396207986697); - assert!((epsilon_k_eff.re() - 21.654396207986697).abs() < 1.0e-6) + println!("{}", epsilon_k_eff - 21.654396207986697); + assert!((epsilon_k_eff - 21.654396207986697).abs() < 1.0e-6) } #[test] fn test_zero_integrand() { let parameters = hydrogen_fh1(); - let temperature = Dual2::from_re(26.706).derive(); + let temperature = 26.706; let sigma_eff = parameters.calc_sigma_eff_ij(0, 0, temperature); let r0 = parameters.zero_integrand(0, 0, temperature, sigma_eff); - println!("{}", r0.re() - 2.5265031901173732); - assert!((r0.re() - 2.5265031901173732).abs() < 5.0e-7) + println!("{}", r0 - 2.5265031901173732); + assert!((r0 - 2.5265031901173732).abs() < 5.0e-7) } #[test] fn test_hs_diameter() { let parameters = hydrogen_fh1(); - let temperature = Dual2::from_re(26.7060).derive(); + let temperature = Dual64::from(26.7060).derivative(); let sigma_eff = parameters.calc_sigma_eff_ij(0, 0, temperature); let d_hs = parameters.hs_diameter_ij(0, 0, temperature, sigma_eff); assert!((d_hs.re() - 3.1410453883283341).abs() < 5.0e-8); - assert!((d_hs.v1[0] + 8.4528823966252661e-3).abs() < 1.0e-9); + assert!((d_hs.eps.unwrap() + 8.4528823966252661e-3).abs() < 1.0e-9); } #[test] diff --git a/src/saftvrqmie/eos/non_additive_hs.rs b/src/saftvrqmie/eos/non_additive_hs.rs index 520a5a975..ae8c8cdea 100644 --- a/src/saftvrqmie/eos/non_additive_hs.rs +++ b/src/saftvrqmie/eos/non_additive_hs.rs @@ -11,7 +11,7 @@ pub struct NonAddHardSphere { pub parameters: Arc, } -impl> HelmholtzEnergyDual for NonAddHardSphere { +impl + Copy> HelmholtzEnergyDual for NonAddHardSphere { fn helmholtz_energy(&self, state: &StateHD) -> D { let p = &self.parameters; let n = p.m.len(); @@ -40,7 +40,7 @@ impl fmt::Display for NonAddHardSphere { } } -pub fn reduced_non_additive_hs_energy>( +pub fn reduced_non_additive_hs_energy + Copy>( parameters: &SaftVRQMieParameters, d_hs_ij: &Array2, d_hs_add_ij: &Array2, diff --git a/src/uvtheory/eos/attractive_perturbation_bh.rs b/src/uvtheory/eos/attractive_perturbation_bh.rs index c99580cb5..b4be61452 100644 --- a/src/uvtheory/eos/attractive_perturbation_bh.rs +++ b/src/uvtheory/eos/attractive_perturbation_bh.rs @@ -51,7 +51,7 @@ impl fmt::Display for AttractivePerturbationBH { } } -impl> HelmholtzEnergyDual for AttractivePerturbationBH { +impl + Copy> HelmholtzEnergyDual for AttractivePerturbationBH { /// Helmholtz energy for attractive perturbation, eq. 52 fn helmholtz_energy(&self, state: &StateHD) -> D { let p = &self.parameters; @@ -87,7 +87,7 @@ fn delta_b12u>(t_x: D, mean_field_constant_x: D, weighted_sigma3 -mean_field_constant_x / t_x * 2.0 * PI * weighted_sigma3_ij } -fn residual_virial_coefficient>(p: &UVParameters, x: &Array1, t: D) -> D { +fn residual_virial_coefficient + Copy>(p: &UVParameters, x: &Array1, t: D) -> D { let mut delta_b2bar = D::zero(); for i in 0..p.ncomponents { let xi = x[i]; @@ -101,7 +101,7 @@ fn residual_virial_coefficient>(p: &UVParameters, x: &Array1, delta_b2bar } -fn correlation_integral_bh>( +fn correlation_integral_bh + Copy>( rho_x: D, mean_field_constant_x: D, rep_x: D, @@ -116,7 +116,7 @@ fn correlation_integral_bh>( /// U-fraction according to Barker-Henderson division. /// Eq. 15 -fn u_fraction_bh>(rep_x: D, reduced_density: D, one_fluid_beta: D) -> D { +fn u_fraction_bh + Copy>(rep_x: D, reduced_density: D, one_fluid_beta: D) -> D { let mut c = [D::zero(); 4]; let inv_rep = rep_x.recip(); for i in 0..4 { @@ -130,11 +130,11 @@ fn u_fraction_bh>(rep_x: D, reduced_density: D, one_fluid_beta: /// Activation function used for u-fraction according to Barker-Henderson division. /// Eq. 16 -fn activation>(c: D, one_fluid_beta: D) -> D { +fn activation + Copy>(c: D, one_fluid_beta: D) -> D { one_fluid_beta * c.sqrt() / (one_fluid_beta.powi(2) * c + 1.0).sqrt() } -fn one_fluid_properties>( +fn one_fluid_properties + Copy>( p: &UVParameters, x: &Array1, t: D, @@ -172,7 +172,7 @@ fn one_fluid_properties>( ) } -fn coefficients_bh>(rep: D, att: D, d: D) -> [D; 3] { +fn coefficients_bh + Copy>(rep: D, att: D, d: D) -> [D; 3] { let c11 = d.powd(-rep + 6.0) * ((D::one() * 2.0f64).powd(-rep + 3.0) - d.powd(rep - 3.0)) / (-rep + 3.0) + (-d.powi(3) * 8.0 + 1.0) / 24.0; @@ -190,14 +190,14 @@ fn coefficients_bh>(rep: D, att: D, d: D) -> [D; 3] { [c1, c2, c3] } -fn delta_b2>(reduced_temperature: D, rep: f64, att: f64) -> D { +fn delta_b2 + Copy>(reduced_temperature: D, rep: f64, att: f64) -> D { let rc = 5.0; let alpha = mean_field_constant(rep, att, rc); let yeff = y_eff(reduced_temperature, rep, att); -(yeff * (rc.powi(3) - 1.0) / 3.0 + reduced_temperature.recip() * alpha) * 2.0 * PI } -fn y_eff>(reduced_temperature: D, rep: f64, att: f64) -> D { +fn y_eff + Copy>(reduced_temperature: D, rep: f64, att: f64) -> D { // optimize: move this part to parameter initialization let rc = 5.0; let rs = 1.0; diff --git a/src/uvtheory/eos/attractive_perturbation_uvb3.rs b/src/uvtheory/eos/attractive_perturbation_uvb3.rs index c703dee6a..af8aab27c 100644 --- a/src/uvtheory/eos/attractive_perturbation_uvb3.rs +++ b/src/uvtheory/eos/attractive_perturbation_uvb3.rs @@ -90,7 +90,7 @@ impl fmt::Display for AttractivePerturbationUVB3 { } } -impl> HelmholtzEnergyDual for AttractivePerturbationUVB3 { +impl + Copy> HelmholtzEnergyDual for AttractivePerturbationUVB3 { /// Helmholtz energy for attractive perturbation fn helmholtz_energy(&self, state: &StateHD) -> D { let p = &self.parameters; @@ -131,7 +131,7 @@ impl> HelmholtzEnergyDual for AttractivePerturbationUVB3 { } } -fn delta_b12u>( +fn delta_b12u + Copy>( t_x: D, mean_field_constant_x: D, weighted_sigma3_ij: D, @@ -144,7 +144,7 @@ fn delta_b12u>( * weighted_sigma3_ij } -fn residual_virial_coefficient>(p: &UVParameters, x: &Array1, t: D) -> D { +fn residual_virial_coefficient + Copy>(p: &UVParameters, x: &Array1, t: D) -> D { let mut delta_b2bar = D::zero(); for i in 0..p.ncomponents { @@ -163,7 +163,7 @@ fn residual_virial_coefficient>(p: &UVParameters, x: &Array1, } delta_b2bar } -fn residual_third_virial_coefficient>( +fn residual_third_virial_coefficient + Copy>( p: &UVParameters, x: &Array1, t: D, @@ -192,7 +192,7 @@ fn residual_third_virial_coefficient>( } delta_b3bar } -fn correlation_integral_wca>( +fn correlation_integral_wca + Copy>( rho_x: D, mean_field_constant_x: D, rep_x: D, @@ -209,7 +209,7 @@ fn correlation_integral_wca>( } /// U-fraction with low temperature correction omega -fn u_fraction_wca>(rep_x: D, reduced_density: D, t_x: D) -> D { +fn u_fraction_wca + Copy>(rep_x: D, reduced_density: D, t_x: D) -> D { let omega = if t_x.re() < 175.0 { (-t_x * CU_WCA[5] * (reduced_density - CU_WCA[6]).powi(2)).exp() * ((t_x * CU_WCA[7]).tanh().recip() - 1.0).powi(2) @@ -223,7 +223,7 @@ fn u_fraction_wca>(rep_x: D, reduced_density: D, t_x: D) -> D { } // Coefficients for IWCA -fn coefficients_wca>(rep: D, att: D, d: D) -> [D; 6] { +fn coefficients_wca + Copy>(rep: D, att: D, d: D) -> [D; 6] { let rep_inv = rep.recip(); let rs_x = (rep / att).powd((rep - att).recip()); let tau_x = -d + rs_x; @@ -261,7 +261,7 @@ fn factorial(num: u64) -> u64 { (1..=num).product() } -fn delta_b2>(reduced_temperature: D, rep: f64, att: f64, q: D) -> D { +fn delta_b2 + Copy>(reduced_temperature: D, rep: f64, att: f64, q: D) -> D { let rm = (rep / att).powd((rep - att).recip()); let beta = reduced_temperature.recip(); let b20 = q.powi(3) * 2.0 / 3.0 * PI; @@ -289,13 +289,13 @@ fn delta_b2>(reduced_temperature: D, rep: f64, att: f64, q: D) - for i in 2..16 { let k = factorial(i as u64) as f64 * i as f64; - sum_beta += beta.powi(i as i32) / k + sum_beta += beta.powi(i) / k } (b20 - rm.powi(3) * 2.0 / 3.0 * PI - c1) * y - sum_beta * c2 - beta * c3 - beta.powi(2) * c4 } -fn delta_b31u>( +fn delta_b31u + Copy>( t_x: D, weighted_sigma3_ij: D, rm_x: D, @@ -313,7 +313,14 @@ fn delta_b31u>( t_x.recip() * 4.0 * mie_prefactor(rep_x, att_x) * PI * k1 * weighted_sigma3_ij.powi(2) } -fn delta_b3>(t_x: D, rm_x: f64, rep_x: f64, _att_x: f64, d_x: D, q_x: D) -> D { +fn delta_b3 + Copy>( + t_x: D, + rm_x: f64, + rep_x: f64, + _att_x: f64, + d_x: D, + q_x: D, +) -> D { let beta = t_x.recip(); let b30 = (q_x.powi(3) * PI / 6.0).powi(2) * 10.0; diff --git a/src/uvtheory/eos/attractive_perturbation_wca.rs b/src/uvtheory/eos/attractive_perturbation_wca.rs index 9904e61a4..cc31bab5f 100644 --- a/src/uvtheory/eos/attractive_perturbation_wca.rs +++ b/src/uvtheory/eos/attractive_perturbation_wca.rs @@ -77,7 +77,7 @@ impl fmt::Display for AttractivePerturbationWCA { } } -impl> HelmholtzEnergyDual for AttractivePerturbationWCA { +impl + Copy> HelmholtzEnergyDual for AttractivePerturbationWCA { /// Helmholtz energy for attractive perturbation, eq. 52 fn helmholtz_energy(&self, state: &StateHD) -> D { let p = &self.parameters; @@ -110,7 +110,7 @@ impl> HelmholtzEnergyDual for AttractivePerturbationWCA { } // (S43) & (S53) -fn delta_b12u>( +fn delta_b12u + Copy>( t_x: D, mean_field_constant_x: D, weighted_sigma3_ij: D, @@ -123,7 +123,7 @@ fn delta_b12u>( * weighted_sigma3_ij } -fn residual_virial_coefficient>(p: &UVParameters, x: &Array1, t: D) -> D { +fn residual_virial_coefficient + Copy>(p: &UVParameters, x: &Array1, t: D) -> D { let mut delta_b2bar = D::zero(); for i in 0..p.ncomponents { @@ -145,7 +145,7 @@ fn residual_virial_coefficient>(p: &UVParameters, x: &Array1, delta_b2bar } -fn correlation_integral_wca>( +fn correlation_integral_wca + Copy>( rho_x: D, mean_field_constant_x: D, rep_x: D, @@ -163,13 +163,13 @@ fn correlation_integral_wca>( /// U-fraction according to Barker-Henderson division. /// Eq. 15 -fn u_fraction_wca>(rep_x: D, reduced_density: D) -> D { +fn u_fraction_wca + Copy>(rep_x: D, reduced_density: D) -> D { (reduced_density * CU_WCA[0] + reduced_density.powi(2) * (rep_x.recip() * CU_WCA[2] + CU_WCA[1])) .tanh() } -pub(super) fn one_fluid_properties>( +pub(super) fn one_fluid_properties + Copy>( p: &UVParameters, x: &Array1, t: D, @@ -212,7 +212,7 @@ pub(super) fn one_fluid_properties>( } // Coefficients for IWCA from eq. (S55) -fn coefficients_wca>(rep: D, att: D, d: D) -> [D; 6] { +fn coefficients_wca + Copy>(rep: D, att: D, d: D) -> [D; 6] { let rep_inv = rep.recip(); let rs_x = (rep / att).powd((rep - att).recip()); let tau_x = -d + rs_x; @@ -244,7 +244,7 @@ fn coefficients_wca>(rep: D, att: D, d: D) -> [D; 6] { [c1, c2, c3, c4, c5, c6] } -fn delta_b2>(reduced_temperature: D, rep: f64, att: f64, q: D) -> D { +fn delta_b2 + Copy>(reduced_temperature: D, rep: f64, att: f64, q: D) -> D { let rm = (rep / att).powf(1.0 / (rep - att)); // Check mixing rule!! let rc = 5.0; let alpha = mean_field_constant(rep, att, rc); @@ -256,7 +256,7 @@ fn delta_b2>(reduced_temperature: D, rep: f64, att: f64, q: D) - * PI } -fn y_eff>(reduced_temperature: D, rep: f64, att: f64) -> D { +fn y_eff + Copy>(reduced_temperature: D, rep: f64, att: f64) -> D { // optimize: move this part to parameter initialization let rc = 5.0; let rs = (rep / att).powf(1.0 / (rep - att)); diff --git a/src/uvtheory/eos/hard_sphere_bh.rs b/src/uvtheory/eos/hard_sphere_bh.rs index a3001cb4b..984577b97 100644 --- a/src/uvtheory/eos/hard_sphere_bh.rs +++ b/src/uvtheory/eos/hard_sphere_bh.rs @@ -23,7 +23,7 @@ pub struct HardSphereBH { pub parameters: Arc, } -impl> HelmholtzEnergyDual for HardSphereBH { +impl + Copy> HelmholtzEnergyDual for HardSphereBH { /// Helmholtz energy for hard spheres, eq. 19 (check Volume) fn helmholtz_energy(&self, state: &StateHD) -> D { let d = diameter_bh(&self.parameters, state.temperature); @@ -46,7 +46,10 @@ impl fmt::Display for HardSphereBH { /// Dimensionless Hard-sphere diameter according to Barker-Henderson division. /// Eq. S23 and S24. /// -pub(super) fn diameter_bh>(parameters: &UVParameters, temperature: D) -> Array1 { +pub(super) fn diameter_bh + Copy>( + parameters: &UVParameters, + temperature: D, +) -> Array1 { parameters .cd_bh_pure .iter() @@ -54,14 +57,16 @@ pub(super) fn diameter_bh>(parameters: &UVParameters, temperatur .map(|(i, c)| { let t = temperature / parameters.epsilon_k[i]; let d = t.powf(0.25) * c[1] + t.powf(0.75) * c[2] + t.powf(1.25) * c[3]; - (t * c[0] + d * (t + 1.0).ln() + t.powi(2) * c[4] + 1.0) - .powf(-0.5 / parameters.rep[i] as f64) + (t * c[0] + d * (t + 1.0).ln() + t.powi(2) * c[4] + 1.0).powf(-0.5 / parameters.rep[i]) * parameters.sigma[i] }) .collect() } -pub(super) fn zeta>(partial_density: &Array1, diameter: &Array1) -> [D; 4] { +pub(super) fn zeta + Copy>( + partial_density: &Array1, + diameter: &Array1, +) -> [D; 4] { let mut zeta: [D; 4] = [D::zero(), D::zero(), D::zero(), D::zero()]; for i in 0..partial_density.len() { for k in 0..4 { @@ -72,7 +77,7 @@ pub(super) fn zeta>(partial_density: &Array1, diameter: &Arra zeta } -pub(super) fn packing_fraction>( +pub(super) fn packing_fraction + Copy>( partial_density: &Array1, diameter: &Array1, ) -> D { @@ -81,7 +86,7 @@ pub(super) fn packing_fraction>( }) } -pub(super) fn zeta_23>(molefracs: &Array1, diameter: &Array1) -> D { +pub(super) fn zeta_23 + Copy>(molefracs: &Array1, diameter: &Array1) -> D { let mut zeta: [D; 2] = [D::zero(), D::zero()]; for i in 0..molefracs.len() { for k in 0..2 { @@ -91,7 +96,7 @@ pub(super) fn zeta_23>(molefracs: &Array1, diameter: &Array1< zeta[0] / zeta[1] } -pub(super) fn packing_fraction_b>( +pub(super) fn packing_fraction_b + Copy>( parameters: &UVParameters, diameter: &Array1, eta: D, @@ -111,7 +116,7 @@ pub(super) fn packing_fraction_b>( }) } -pub(super) fn packing_fraction_a>( +pub(super) fn packing_fraction_a + Copy>( parameters: &UVParameters, diameter: &Array1, eta: D, diff --git a/src/uvtheory/eos/hard_sphere_wca.rs b/src/uvtheory/eos/hard_sphere_wca.rs index 8a8e57270..ed636eb6b 100644 --- a/src/uvtheory/eos/hard_sphere_wca.rs +++ b/src/uvtheory/eos/hard_sphere_wca.rs @@ -53,7 +53,7 @@ pub struct HardSphereWCA { pub parameters: Arc, } -impl> HelmholtzEnergyDual for HardSphereWCA { +impl + Copy> HelmholtzEnergyDual for HardSphereWCA { /// Helmholtz energy for hard spheres, eq. 19 (check Volume) fn helmholtz_energy(&self, state: &StateHD) -> D { let d = diameter_wca(&self.parameters, state.temperature); @@ -74,7 +74,7 @@ impl fmt::Display for HardSphereWCA { } /// Dimensionless Hard-sphere diameter according to Weeks-Chandler-Andersen division. -pub(super) fn diameter_wca>( +pub(super) fn diameter_wca + Copy>( parameters: &UVParameters, temperature: D, ) -> Array1 { @@ -96,7 +96,11 @@ pub(super) fn diameter_wca>( .collect() } -pub(super) fn dimensionless_diameter_q_wca>(t_x: D, rep_x: D, att_x: D) -> D { +pub(super) fn dimensionless_diameter_q_wca + Copy>( + t_x: D, + rep_x: D, + att_x: D, +) -> D { let nu = rep_x; let n = att_x; let rs = (nu / n).powd((nu - n).recip()); @@ -122,7 +126,10 @@ pub(super) fn dimensionless_diameter_q_wca>(t_x: D, rep_x: D, at * rs } -pub(super) fn zeta>(partial_density: &Array1, diameter: &Array1) -> [D; 4] { +pub(super) fn zeta + Copy>( + partial_density: &Array1, + diameter: &Array1, +) -> [D; 4] { let mut zeta: [D; 4] = [D::zero(), D::zero(), D::zero(), D::zero()]; for i in 0..partial_density.len() { for k in 0..4 { @@ -133,7 +140,7 @@ pub(super) fn zeta>(partial_density: &Array1, diameter: &Arra zeta } -pub(super) fn packing_fraction>( +pub(super) fn packing_fraction + Copy>( partial_density: &Array1, diameter: &Array1, ) -> D { @@ -142,7 +149,7 @@ pub(super) fn packing_fraction>( }) } -pub(super) fn zeta_23>(molefracs: &Array1, diameter: &Array1) -> D { +pub(super) fn zeta_23 + Copy>(molefracs: &Array1, diameter: &Array1) -> D { let mut zeta: [D; 2] = [D::zero(), D::zero()]; for i in 0..molefracs.len() { for k in 0..2 { @@ -153,7 +160,7 @@ pub(super) fn zeta_23>(molefracs: &Array1, diameter: &Array1< } #[inline] -pub(super) fn dimensionless_length_scale>( +pub(super) fn dimensionless_length_scale + Copy>( parameters: &UVParameters, temperature: D, ) -> Array1 { @@ -172,7 +179,7 @@ pub(super) fn dimensionless_length_scale>( #[inline] -pub(super) fn packing_fraction_b>( +pub(super) fn packing_fraction_b + Copy>( parameters: &UVParameters, eta: D, temperature: D, @@ -194,7 +201,7 @@ pub(super) fn packing_fraction_b>( }) } -pub(super) fn packing_fraction_b_uvb3>( +pub(super) fn packing_fraction_b_uvb3 + Copy>( parameters: &UVParameters, eta: D, temperature: D, @@ -216,7 +223,7 @@ pub(super) fn packing_fraction_b_uvb3>( }) } -pub(super) fn packing_fraction_a>( +pub(super) fn packing_fraction_a + Copy>( parameters: &UVParameters, eta: D, temperature: D, @@ -245,7 +252,7 @@ pub(super) fn packing_fraction_a>( }) } -pub(super) fn packing_fraction_a_uvb3>( +pub(super) fn packing_fraction_a_uvb3 + Copy>( parameters: &UVParameters, eta: D, temperature: D, diff --git a/src/uvtheory/eos/reference_perturbation_bh.rs b/src/uvtheory/eos/reference_perturbation_bh.rs index e54ca0c2e..37d7f67a9 100644 --- a/src/uvtheory/eos/reference_perturbation_bh.rs +++ b/src/uvtheory/eos/reference_perturbation_bh.rs @@ -18,7 +18,7 @@ impl fmt::Display for ReferencePerturbationBH { } } -impl> HelmholtzEnergyDual for ReferencePerturbationBH { +impl + Copy> HelmholtzEnergyDual for ReferencePerturbationBH { /// Helmholtz energy for perturbation reference (Mayer-f), eq. 29 fn helmholtz_energy(&self, state: &StateHD) -> D { let p = &self.parameters; diff --git a/src/uvtheory/eos/reference_perturbation_uvb3.rs b/src/uvtheory/eos/reference_perturbation_uvb3.rs index 30d79d126..8ebfac45e 100644 --- a/src/uvtheory/eos/reference_perturbation_uvb3.rs +++ b/src/uvtheory/eos/reference_perturbation_uvb3.rs @@ -19,7 +19,7 @@ impl fmt::Display for ReferencePerturbationUVB3 { } } -impl> HelmholtzEnergyDual for ReferencePerturbationUVB3 { +impl + Copy> HelmholtzEnergyDual for ReferencePerturbationUVB3 { /// Helmholtz energy for perturbation reference (Mayer-f), eq. 29 fn helmholtz_energy(&self, state: &StateHD) -> D { let p = &self.parameters; diff --git a/src/uvtheory/eos/reference_perturbation_wca.rs b/src/uvtheory/eos/reference_perturbation_wca.rs index 89fb58801..4b709424e 100644 --- a/src/uvtheory/eos/reference_perturbation_wca.rs +++ b/src/uvtheory/eos/reference_perturbation_wca.rs @@ -19,7 +19,7 @@ impl fmt::Display for ReferencePerturbationWCA { } } -impl> HelmholtzEnergyDual for ReferencePerturbationWCA { +impl + Copy> HelmholtzEnergyDual for ReferencePerturbationWCA { /// Helmholtz energy for perturbation reference (Mayer-f), eq. 29 fn helmholtz_energy(&self, state: &StateHD) -> D { let p = &self.parameters; diff --git a/src/uvtheory/parameters.rs b/src/uvtheory/parameters.rs index 4b54c22a7..f602084d6 100644 --- a/src/uvtheory/parameters.rs +++ b/src/uvtheory/parameters.rs @@ -93,12 +93,12 @@ lazy_static! { } #[inline] -pub fn mie_prefactor>(rep: D, att: D) -> D { +pub fn mie_prefactor + Copy>(rep: D, att: D) -> D { rep / (rep - att) * (rep / att).powd(att / (rep - att)) } #[inline] -pub fn mean_field_constant>(rep: D, att: D, x: D) -> D { +pub fn mean_field_constant + Copy>(rep: D, att: D, x: D) -> D { mie_prefactor(rep, att) * (x.powd(-att + 3.0) / (att - 3.0) - x.powd(-rep + 3.0) / (rep - 3.0)) }