From f9cd39e080bc3e2cde2e3b4d5f399953ed124b5b Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Fri, 16 Oct 2026 14:07:17 +0200 Subject: [PATCH 1/3] Removed materialization of geometry coefficients. Lifted zeta computation out of individual contributions. Additional refactoring to remove allocations in other places. --- crates/feos-core/src/equation_of_state/mod.rs | 8 ++ .../src/equation_of_state/residual.rs | 16 ++++ .../src/phase_equilibria/bubble_dew.rs | 3 + crates/feos-core/src/state/mod.rs | 9 ++- crates/feos-derive/src/residual.rs | 11 +++ crates/feos/src/epcsaft/parameters.rs | 2 +- crates/feos/src/hard_sphere/mod.rs | 80 +++++++++++++------ crates/feos/src/pcsaft/eos/dispersion.rs | 31 ++++--- crates/feos/src/pcsaft/eos/hard_chain.rs | 36 ++++++--- crates/feos/src/pcsaft/eos/mod.rs | 61 +++++++++----- crates/feos/src/pcsaft/eos/polar.rs | 4 +- crates/feos/src/saftvrmie/parameters.rs | 2 +- 12 files changed, 187 insertions(+), 76 deletions(-) diff --git a/crates/feos-core/src/equation_of_state/mod.rs b/crates/feos-core/src/equation_of_state/mod.rs index 14e45c07a..f2746ed86 100644 --- a/crates/feos-core/src/equation_of_state/mod.rs +++ b/crates/feos-core/src/equation_of_state/mod.rs @@ -64,6 +64,14 @@ impl ResidualDyn for EquationOfState, R> { self.residual .reduced_helmholtz_energy_density_contributions(state) } + + fn reduced_residual_helmholtz_energy_density + Copy>( + &self, + state: &StateHD, + ) -> D { + self.residual + .reduced_residual_helmholtz_energy_density(state) + } } impl Subset for EquationOfState, R> { diff --git a/crates/feos-core/src/equation_of_state/residual.rs b/crates/feos-core/src/equation_of_state/residual.rs index af5fc88f6..270db60c0 100644 --- a/crates/feos-core/src/equation_of_state/residual.rs +++ b/crates/feos-core/src/equation_of_state/residual.rs @@ -70,6 +70,19 @@ pub trait ResidualDyn { &self, state: &StateHD, ) -> Vec<(&'static str, D)>; + + /// Evaluate the residual reduced Helmholtz energy density $\beta f^\mathrm{res}$. + /// + /// The default sums [ResidualDyn::reduced_helmholtz_energy_density_contributions]. + /// Override it to avoid allocating the contributions in performance critical code. + fn reduced_residual_helmholtz_energy_density + Copy>( + &self, + state: &StateHD, + ) -> D { + self.reduced_helmholtz_energy_density_contributions(state) + .iter() + .fold(D::zero(), |acc, (_, a)| acc + a) + } } impl + Clone, T: ResidualDyn, D: DualNum + Copy> Residual for C { @@ -93,6 +106,9 @@ impl + Clone, T: ResidualDyn, D: DualNum + Copy> Resid ) -> Vec<(&'static str, D)> { ResidualDyn::reduced_helmholtz_energy_density_contributions(self.deref(), state) } + fn reduced_residual_helmholtz_energy_density(&self, state: &StateHD) -> D { + ResidualDyn::reduced_residual_helmholtz_energy_density(self.deref(), state) + } } /// A residual Helmholtz energy model. diff --git a/crates/feos-core/src/phase_equilibria/bubble_dew.rs b/crates/feos-core/src/phase_equilibria/bubble_dew.rs index f956d4d7d..030881ede 100644 --- a/crates/feos-core/src/phase_equilibria/bubble_dew.rs +++ b/crates/feos-core/src/phase_equilibria/bubble_dew.rs @@ -949,6 +949,9 @@ fn log_iteration( x2: &[f64], newton: bool, ) { + if verbosity < Verbosity::Iter { + return; + } let error = error.map_or_else(|| format!("{:14}", ""), |e| format!("{:<14.8e}", e.abs())); log_iter!( verbosity, diff --git a/crates/feos-core/src/state/mod.rs b/crates/feos-core/src/state/mod.rs index daa7a7ff4..ef63dacd2 100644 --- a/crates/feos-core/src/state/mod.rs +++ b/crates/feos-core/src/state/mod.rs @@ -82,17 +82,20 @@ where { /// Create a new `StateHD` for given temperature, volume and composition. pub fn new(temperature: D, volume: D, moles: &OVector) -> Self { - Self::new_density(temperature, &(moles / volume)) + Self::from_partial_density(temperature, moles / volume) } /// Create a new `StateHD` for given temperature and partial densities pub fn new_density(temperature: D, partial_density: &OVector) -> Self { - let molefracs = partial_density / partial_density.sum(); + Self::from_partial_density(temperature, partial_density.clone()) + } + fn from_partial_density(temperature: D, partial_density: OVector) -> Self { + let molefracs = &partial_density / partial_density.sum(); Self { temperature, molefracs, - partial_density: partial_density.clone(), + partial_density, } } diff --git a/crates/feos-derive/src/residual.rs b/crates/feos-derive/src/residual.rs index df3f5367f..f478bd91b 100644 --- a/crates/feos-derive/src/residual.rs +++ b/crates/feos-derive/src/residual.rs @@ -42,6 +42,12 @@ fn impl_residual( Self::#name(residual) => residual.reduced_helmholtz_energy_density_contributions(state) } }); + let reduced_residual_helmholtz_energy_density = variants.iter().map(|v| { + let name = &v.ident; + quote! { + Self::#name(residual) => residual.reduced_residual_helmholtz_energy_density(state) + } + }); // let subset = variants.iter().map(|v| { // let name = &v.ident; // quote! { @@ -66,6 +72,11 @@ fn impl_residual( #(#reduced_helmholtz_energy_density_contributions,)* } } + fn reduced_residual_helmholtz_energy_density + Copy>(&self, state: &StateHD) -> D { + match self { + #(#reduced_residual_helmholtz_energy_density,)* + } + } // fn subset(&self, component_list: &[usize]) -> Self { // match self { // #(#subset,)* diff --git a/crates/feos/src/epcsaft/parameters.rs b/crates/feos/src/epcsaft/parameters.rs index 364de81e1..94da71e4d 100644 --- a/crates/feos/src/epcsaft/parameters.rs +++ b/crates/feos/src/epcsaft/parameters.rs @@ -358,7 +358,7 @@ impl ElectrolytePcSaftPars { impl HardSphereProperties for ElectrolytePcSaftPars { fn monomer_shape>(&self, _: N) -> MonomerShape<'_, N> { - MonomerShape::NonSpherical(self.m.map(N::from)) + MonomerShape::NonSpherical(&self.m) } fn hs_diameter + Copy>(&self, temperature: D) -> DVector { diff --git a/crates/feos/src/hard_sphere/mod.rs b/crates/feos/src/hard_sphere/mod.rs index 40d4936b4..0eb94586e 100644 --- a/crates/feos/src/hard_sphere/mod.rs +++ b/crates/feos/src/hard_sphere/mod.rs @@ -17,7 +17,7 @@ pub enum MonomerShape<'a, D> { Spherical(usize), /// For non-spherical molecules in a homosegmented approach, the /// chain length parameter $m$. - NonSpherical(DVector), + NonSpherical(&'a DVector), /// For non-spherical molecules in a heterosegmented approach, /// the geometry factors for every segment and the component /// index for every segment. @@ -50,7 +50,10 @@ pub trait HardSphereProperties { let m = DVector::from_element(n, D::from(1.0)); [m.clone(), m.clone(), m.clone(), m] } - MonomerShape::NonSpherical(m) => [m.clone(), m.clone(), m.clone(), m], + MonomerShape::NonSpherical(m) => { + let m = m.map(D::from); + [m.clone(), m.clone(), m.clone(), m] + } MonomerShape::Heterosegmented(g, _) => g, } } @@ -62,20 +65,56 @@ pub trait HardSphereProperties { partial_density: &DVector, k: [i32; N], ) -> [D; N] { - let component_index = self.component_index(); - let geometry_coefficients = self.geometry_coefficients(temperature); let diameter = self.hs_diameter(temperature); - let mut zeta = [D::zero(); N]; - for i in 0..diameter.len() { - for (z, &k) in zeta.iter_mut().zip(k.iter()) { - *z += partial_density[component_index[i]] - * diameter[i].powi(k) - * (geometry_coefficients[k as usize][i] * FRAC_PI_6); + weighted_zeta( + self.monomer_shape(temperature), + &diameter, + partial_density, + k, + ) + } +} + +/// $\frac{\pi}{6}\sum_\alpha C_{k,\alpha} w_{i(\alpha)}d_\alpha^k$ for every requested $k$, +/// where the weights $w_i$ are either partial densities or mole fractions. +/// +/// Matches on the [MonomerShape] directly instead of materializing +/// [HardSphereProperties::component_index] and +/// [HardSphereProperties::geometry_coefficients], which would allocate. +#[inline] +fn weighted_zeta + Copy, const N: usize>( + monomer_shape: MonomerShape<'_, D>, + diameter: &DVector, + weights: &DVector, + k: [i32; N], +) -> [D; N] { + let mut zeta = [D::zero(); N]; + match monomer_shape { + MonomerShape::Spherical(_) => { + for i in 0..diameter.len() { + for (z, &k) in zeta.iter_mut().zip(k.iter()) { + *z += weights[i] * diameter[i].powi(k) * FRAC_PI_6; + } + } + } + MonomerShape::NonSpherical(m) => { + for i in 0..diameter.len() { + for (z, &k) in zeta.iter_mut().zip(k.iter()) { + *z += weights[i] * diameter[i].powi(k) * (m[i] * FRAC_PI_6); + } + } + } + MonomerShape::Heterosegmented(g, component_index) => { + for i in 0..diameter.len() { + for (z, &k) in zeta.iter_mut().zip(k.iter()) { + *z += weights[component_index[i]] + * diameter[i].powi(k) + * (g[k as usize][i] * FRAC_PI_6); + } } } - - zeta } + zeta } /// Implementation of the BMCSL equation of state for hard-sphere mixtures. @@ -103,17 +142,12 @@ impl HardSphere { ) -> (D, [D; 4], DVector) { let p = parameters; let diameter = p.hs_diameter(state.temperature); - - let component_index = p.component_index(); - let geometry_coefficients = p.geometry_coefficients(state.temperature); - let mut zeta = [D::zero(); 4]; - for i in 0..diameter.len() { - for (z, &k) in zeta.iter_mut().zip([0, 1, 2, 3].iter()) { - *z += state.molefracs[component_index[i]] - * diameter[i].powi(k) - * (geometry_coefficients[k as usize][i] * FRAC_PI_6); - } - } + let mut zeta = weighted_zeta( + p.monomer_shape(state.temperature), + &diameter, + &state.molefracs, + [0, 1, 2, 3], + ); let zeta_23 = zeta[2] / zeta[3]; let density = state.partial_density.sum(); zeta.iter_mut().for_each(|z| *z *= density); diff --git a/crates/feos/src/pcsaft/eos/dispersion.rs b/crates/feos/src/pcsaft/eos/dispersion.rs index 50fb3c1fc..32f37980a 100644 --- a/crates/feos/src/pcsaft/eos/dispersion.rs +++ b/crates/feos/src/pcsaft/eos/dispersion.rs @@ -1,5 +1,4 @@ use super::PcSaftPars; -use crate::hard_sphere::HardSphereProperties; use feos_core::StateHD; use num_dual::DualNum; use std::f64::consts::PI; @@ -62,22 +61,24 @@ pub const B2: [f64; 7] = [ pub struct Dispersion; impl Dispersion { + /// `eta` is the packing fraction $\zeta_3$, shared with the other contributions. pub fn helmholtz_energy_density + Copy>( &self, parameters: &PcSaftPars, state: &StateHD, + eta: D, ) -> D { // auxiliary variables let n = parameters.m.len(); let p = parameters; let rho = &state.partial_density; - // packing fraction - let m = p.m.map(D::from); - let [eta] = parameters.zeta(state.temperature, &state.partial_density, [3]); - // mean segment number - let m = state.molefracs.dot(&m); + let m = state + .molefracs + .iter() + .zip(&p.m) + .fold(D::zero(), |acc, (&x, &m)| acc + x * m); // inverse temperature let t_inv = state.temperature.recip(); @@ -118,14 +119,20 @@ impl Dispersion { } #[cfg(test)] -mod tests { +pub(crate) mod tests { use super::*; + use crate::hard_sphere::HardSphere; use crate::pcsaft::parameters::utils::{ butane_parameters, propane_butane_parameters, propane_parameters, }; use approx::assert_relative_eq; use nalgebra::dvector; + pub(crate) fn a_disp(p: &PcSaftPars, s: &StateHD) -> f64 { + let (_, zeta, _) = HardSphere.helmholtz_energy_density_and_properties(p, s); + Dispersion.helmholtz_energy_density(p, s, zeta[3]) + } + #[test] fn helmholtz_energy() { let params = &propane_parameters().params; @@ -133,7 +140,7 @@ mod tests { let v = 1000.0; let n = 1.0; let s = StateHD::new(t, v, &dvector![n]); - let a_rust = Dispersion.helmholtz_energy_density(params, &s) * v; + let a_rust = a_disp(params, &s) * v; assert_relative_eq!(a_rust, -1.0622531100351962, epsilon = 1e-10); } @@ -146,12 +153,12 @@ mod tests { let v = 1000.0; let n = 1.0; let s = StateHD::new(t, v, &dvector![n]); - let a1 = Dispersion.helmholtz_energy_density(propane, &s); - let a2 = Dispersion.helmholtz_energy_density(butane, &s); + let a1 = a_disp(propane, &s); + let a2 = a_disp(butane, &s); let s1m = StateHD::new(t, v, &dvector![n, 0.0]); - let a1m = Dispersion.helmholtz_energy_density(mix, &s1m); + let a1m = a_disp(mix, &s1m); let s2m = StateHD::new(t, v, &dvector![0.0, n]); - let a2m = Dispersion.helmholtz_energy_density(mix, &s2m); + let a2m = a_disp(mix, &s2m); assert_relative_eq!(a1, a1m, epsilon = 1e-14); assert_relative_eq!(a2, a2m, epsilon = 1e-14); } diff --git a/crates/feos/src/pcsaft/eos/hard_chain.rs b/crates/feos/src/pcsaft/eos/hard_chain.rs index e894f695a..04c175cfb 100644 --- a/crates/feos/src/pcsaft/eos/hard_chain.rs +++ b/crates/feos/src/pcsaft/eos/hard_chain.rs @@ -1,45 +1,55 @@ -use crate::hard_sphere::HardSphereProperties; use crate::pcsaft::parameters::PcSaftPars; use feos_core::StateHD; +use nalgebra::DVector; use num_dual::*; pub struct HardChain; impl HardChain { + /// `diameter` and `zeta` are the hard-sphere diameters and packing + /// fractions, shared with the other contributions. #[inline] pub fn helmholtz_energy_density + Copy>( &self, parameters: &PcSaftPars, state: &StateHD, + diameter: &DVector, + zeta: [D; 4], ) -> D { - let d = parameters.hs_diameter(state.temperature); - let [zeta2, zeta3] = parameters.zeta(state.temperature, &state.partial_density, [2, 3]); + let [_, _, zeta2, zeta3] = zeta; let frac_1mz3 = -(zeta3 - 1.0).recip(); let c = zeta2 * frac_1mz3 * frac_1mz3; - let g_hs = d.map(|d| frac_1mz3 + d * c * 1.5 - d.powi(2) * c.powi(2) * (zeta3 - 1.0) * 0.5); - state - .partial_density - .component_mul(&g_hs.map(|g| g.ln())) - .dot(&(-parameters.m.map(|m| D::from(m - 1.0)))) + let mut a = D::zero(); + for (i, &d) in diameter.iter().enumerate() { + let g_hs = frac_1mz3 + d * c * 1.5 - d.powi(2) * c.powi(2) * (zeta3 - 1.0) * 0.5; + a -= state.partial_density[i] * g_hs.ln() * (parameters.m[i] - 1.0); + } + a } } #[cfg(test)] mod tests { use super::*; + use crate::hard_sphere::HardSphere; use crate::pcsaft::parameters::utils::{ butane_parameters, propane_butane_parameters, propane_parameters, }; use approx::assert_relative_eq; use nalgebra::dvector; + fn a_hc(p: &PcSaftPars, s: &StateHD) -> f64 { + let (_, zeta, d) = HardSphere.helmholtz_energy_density_and_properties(p, s); + HardChain.helmholtz_energy_density(p, s, &d, zeta) + } + #[test] fn helmholtz_energy() { let t = 250.0; let v = 1000.0; let n = 1.0; let s = StateHD::new(t, v, &dvector![n]); - let a_rust = HardChain.helmholtz_energy_density(&propane_parameters().params, &s) * v; + let a_rust = a_hc(&propane_parameters().params, &s) * v; assert_relative_eq!(a_rust, -0.12402626171926148, epsilon = 1e-10); } @@ -49,12 +59,12 @@ mod tests { let v = 2.5; let n = 1.0; let s = StateHD::new(t, v, &dvector![n]); - let a1 = HardChain.helmholtz_energy_density(&propane_parameters().params, &s); - let a2 = HardChain.helmholtz_energy_density(&butane_parameters().params, &s); + let a1 = a_hc(&propane_parameters().params, &s); + let a2 = a_hc(&butane_parameters().params, &s); let s1m = StateHD::new(t, v, &dvector![n, 0.0]); - let a1m = HardChain.helmholtz_energy_density(&propane_butane_parameters().params, &s1m); + let a1m = a_hc(&propane_butane_parameters().params, &s1m); let s2m = StateHD::new(t, v, &dvector![0.0, n]); - let a2m = HardChain.helmholtz_energy_density(&propane_butane_parameters().params, &s2m); + let a2m = a_hc(&propane_butane_parameters().params, &s2m); assert_relative_eq!(a1, a1m, epsilon = 1e-14); assert_relative_eq!(a2, a2m, epsilon = 1e-14); } diff --git a/crates/feos/src/pcsaft/eos/mod.rs b/crates/feos/src/pcsaft/eos/mod.rs index 6a68f67f7..02e5c0162 100644 --- a/crates/feos/src/pcsaft/eos/mod.rs +++ b/crates/feos/src/pcsaft/eos/mod.rs @@ -101,46 +101,66 @@ impl ResidualDyn for PcSaft { state: &StateHD, ) -> Vec<(&'static str, D)> { let mut v = Vec::with_capacity(7); - let d = self.params.hs_diameter(state.temperature); + self.for_each_contribution(state, |name, a| v.push((name, a))); + v + } - v.push(( - "Hard Sphere", - HardSphere.helmholtz_energy_density(&self.params, state), - )); + fn reduced_residual_helmholtz_energy_density + Copy>( + &self, + state: &StateHD, + ) -> D { + let mut a = D::zero(); + self.for_each_contribution(state, |_, a_i| a += a_i); + a + } +} + +impl PcSaft { + /// Evaluate every Helmholtz energy contribution and pass it to `f` + /// together with its name. + #[inline] + fn for_each_contribution + Copy>( + &self, + state: &StateHD, + mut f: impl FnMut(&'static str, D), + ) { + let (a_hs, zeta, d) = + HardSphere.helmholtz_energy_density_and_properties(&self.params, state); + f("Hard Sphere", a_hs); if self.hard_chain { - v.push(( + f( "Hard Chain", - HardChain.helmholtz_energy_density(&self.params, state), - )) + HardChain.helmholtz_energy_density(&self.params, state, &d, zeta), + ) } - v.push(( + f( "Dispersion", - Dispersion.helmholtz_energy_density(&self.params, state), - )); + Dispersion.helmholtz_energy_density(&self.params, state, zeta[3]), + ); if self.dipole { - v.push(( + f( "Dipole", Dipole.helmholtz_energy_density(&self.params, state), - )) + ) } if self.quadrupole { - v.push(( + f( "Quadrupole", Quadrupole.helmholtz_energy_density(&self.params, state), - )) + ) } if self.dipole_quadrupole { - v.push(( + f( "DipoleQuadrupole", DipoleQuadrupole.helmholtz_energy_density( &self.params, state, self.options.dq_variant, ), - )) + ) } if let Some(association) = self.association.as_ref() { - v.push(( + f( "Association", association.helmholtz_energy_density( &self.params, @@ -148,9 +168,8 @@ impl ResidualDyn for PcSaft { state, &d, ), - )) + ) } - v } } @@ -168,7 +187,7 @@ impl Molarweight for PcSaft { impl HardSphereProperties for PcSaftPars { fn monomer_shape>(&self, _: N) -> MonomerShape<'_, N> { - MonomerShape::NonSpherical(self.m.map(N::from)) + MonomerShape::NonSpherical(&self.m) } fn hs_diameter + Copy>(&self, temperature: D) -> DVector { diff --git a/crates/feos/src/pcsaft/eos/polar.rs b/crates/feos/src/pcsaft/eos/polar.rs index 096bdc1c2..6f23f9b61 100644 --- a/crates/feos/src/pcsaft/eos/polar.rs +++ b/crates/feos/src/pcsaft/eos/polar.rs @@ -424,7 +424,7 @@ impl DipoleQuadrupole { mod tests { use super::*; use crate::pcsaft::PcSaftParameters; - use crate::pcsaft::eos::dispersion::Dispersion; + use crate::pcsaft::eos::dispersion::tests::a_disp; use crate::pcsaft::parameters::utils::{ carbon_dioxide_parameters, dme_co2_parameters, dme_parameters, }; @@ -564,7 +564,7 @@ mod tests { let n = n_dme + n_co2; let s = StateHD::new(t, v / n, &(dvector![n_dme, n_co2] / n)); // let a_dpqp = dpqp.helmholtz_energy_density(&s)*v; - let a_disp = Dispersion.helmholtz_energy_density(&mix, &s) * v; + let a_disp = a_disp(&mix, &s) * v; let a_qp = Quadrupole.helmholtz_energy_density(&mix, &s) * v; let a_dp = Dipole.helmholtz_energy_density(&mix, &s) * v; assert_relative_eq!(a_disp, -1.6283622072860044, epsilon = 1e-6); diff --git a/crates/feos/src/saftvrmie/parameters.rs b/crates/feos/src/saftvrmie/parameters.rs index b40cd524b..a9bc98f32 100644 --- a/crates/feos/src/saftvrmie/parameters.rs +++ b/crates/feos/src/saftvrmie/parameters.rs @@ -244,7 +244,7 @@ fn beta_u_mie + Copy>(r: D, la: f64, lr: f64, sigma: f64, c_eps_ impl HardSphereProperties for SaftVRMiePars { fn monomer_shape>(&self, _: N) -> MonomerShape<'_, N> { - MonomerShape::NonSpherical(self.m.map(N::from)) + MonomerShape::NonSpherical(&self.m) } fn hs_diameter + Copy>(&self, temperature: D) -> DVector { From b29b5dd72d498be1f5705cd640cbd9b2bb06aa86 Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Fri, 16 Oct 2026 14:27:45 +0200 Subject: [PATCH 2/3] Use mutable buffers allocated once for cross association and port LU from num-dual (for now). Remove some nalgebra methods that allocate in the hot loops in favor of iterator methods that do not allocate. Lift zeta out of polar contributions. Refactor polar contributions to eliminate repeated computations that can be done once. --- crates/feos/src/association/mod.rs | 126 +++++++--- crates/feos/src/pcsaft/eos/mod.rs | 29 +-- crates/feos/src/pcsaft/eos/polar.rs | 372 ++++++++++++++++++---------- 3 files changed, 351 insertions(+), 176 deletions(-) diff --git a/crates/feos/src/association/mod.rs b/crates/feos/src/association/mod.rs index 161cae031..063ee2782 100644 --- a/crates/feos/src/association/mod.rs +++ b/crates/feos/src/association/mod.rs @@ -4,7 +4,7 @@ use crate::hard_sphere::HardSphereProperties; use feos_core::parameter::{AssociationParameters, AssociationSite, BinaryParameters}; use feos_core::{FeosError, FeosResult, StateHD}; use nalgebra::{DMatrix, DVector}; -use num_dual::linalg::LU; +use num_dual::linalg::LinAlgError; use num_dual::*; #[cfg(feature = "dft")] @@ -208,9 +208,12 @@ impl Association { let delta_ab_re = delta_ab.map(|d| d.re()); let delta_cc_re = delta_cc.map(|d| d.re()); let rho_re = rho.map(|r| r.re()); + let n = rho.len(); + let (mut g, mut h) = (DVector::zeros(n), DMatrix::zeros(n, n)); for k in 0..self.max_iter { if Self::newton_step_cross_association( &mut x, + (&mut g, &mut h), &delta_ab_re, &delta_cc_re, &rho_re, @@ -225,8 +228,16 @@ impl Association { // calculate derivatives let mut x_dual = x.map(D::from); + let (mut g, mut h) = (DVector::zeros(n), DMatrix::zeros(n, n)); for _ in 0..D::NDERIV { - Self::newton_step_cross_association(&mut x_dual, delta_ab, delta_cc, rho, self.tol)?; + Self::newton_step_cross_association( + &mut x_dual, + (&mut g, &mut h), + delta_ab, + delta_cc, + rho, + self.tol, + )?; } // save monomer fraction @@ -236,65 +247,71 @@ impl Association { // Helmholtz energy density let f = |x: D| x.ln() - x * 0.5 + 0.5; - Ok(rho.dot(&x_dual.map(f))) + Ok(rho + .iter() + .zip(&x_dual) + .fold(D::zero(), |acc, (&rho, &x)| acc + rho * f(x))) } + /// One Newton step for the monomer fractions `x`. `g` and `h` are scratch + /// buffers for the gradient and the Hessian, reused across steps. fn newton_step_cross_association + Copy>( x: &mut DVector, + (g, h): (&mut DVector, &mut DMatrix), delta_ab: &DMatrix, delta_cc: &DMatrix, rho: &DVector, tol: f64, ) -> FeosResult { let nassoc = x.len(); - // gradient - let mut g = x.map(|x| x.recip()); - // Hessian - let mut h: DMatrix = DMatrix::zeros(nassoc, nassoc); + h.fill(D::zero()); - // split arrays + // sites are ordered [A sites, B sites, C sites] let (a, b) = delta_ab.shape(); let (c, _) = delta_cc.shape(); - let (xa, xc) = x.rows_range_pair(..a + b, a + b..); - let (xa, xb) = xa.rows_range_pair(..a, a..); - let (rhoa, rhoc) = rho.rows_range_pair(..a + b, a + b..); - let (rhoa, rhob) = rhoa.rows_range_pair(..a, a..); + let (b0, c0) = (a, a + b); for i in 0..nassoc { // calculate gradients - let dnx = if i < a { - let d = delta_ab.row(i).transpose(); - xb.component_mul(&rhob).dot(&d) + 1.0 - } else if i < a + b { - let d = delta_ab.column(i - a); - xa.component_mul(&rhoa).dot(&d) + 1.0 + let dnx = if i < b0 { + (0..b).fold(D::zero(), |acc, j| { + acc + x[b0 + j] * rho[b0 + j] * delta_ab[(i, j)] + }) + } else if i < c0 { + (0..a).fold(D::zero(), |acc, j| { + acc + x[j] * rho[j] * delta_ab[(j, i - b0)] + }) } else { - let d = delta_cc.column(i - a - b); - xc.component_mul(&rhoc).dot(&d) + 1.0 - }; - g[i] -= dnx; + (0..c).fold(D::zero(), |acc, j| { + acc + x[c0 + j] * rho[c0 + j] * delta_cc[(j, i - c0)] + }) + } + 1.0; + g[i] = x[i].recip() - dnx; // approximate hessian h[(i, i)] = -dnx / x[i]; - if i < a { + if i < b0 { for j in 0..b { - h[(i, a + j)] = -delta_ab[(i, j)] * rhob[j]; + h[(i, b0 + j)] = -delta_ab[(i, j)] * rho[b0 + j]; } - } else if i < a + b { + } else if i < c0 { for j in 0..a { - h[(i, j)] = -delta_ab[(j, i - a)] * rhoa[j]; + h[(i, j)] = -delta_ab[(j, i - b0)] * rho[j]; } } else { for j in 0..c { - h[(i, a + b + j)] -= delta_cc[(i - a - b, j)] * rhoc[j]; + h[(i, c0 + j)] -= delta_cc[(i - c0, j)] * rho[c0 + j]; } } } + // check convergence (before g is overwritten by the Newton step) + let converged = g.iter().map(|g| g.re().powi(2)).sum::().sqrt() < tol; + // Newton step // avoid stepping to negative values for x (see Michelsen 2006) - let delta_x = LU::new(h)?.solve(&g); - x.iter_mut().zip(&delta_x).for_each(|(x, &delta_x)| { + solve_in_place(h, g)?; + x.iter_mut().zip(g.iter()).for_each(|(x, &delta_x)| { if delta_x.re() < x.re() * 0.8 { *x -= delta_x } else { @@ -302,9 +319,56 @@ impl Association { } }); - // check convergence - Ok(g.map(|g| g.re()).norm() < tol) + Ok(converged) + } +} + +/// Solve $Ax=b$ in place: `b` is overwritten with $x$ and `a` with its LU factors. +/// +/// Gaussian elimination with partial pivoting, performing the same operations +/// as [num_dual::linalg::LU] (row swaps are applied to `b` directly instead of +/// being recorded in a permutation vector), but without allocating. +fn solve_in_place + Copy>( + a: &mut DMatrix, + b: &mut DVector, +) -> Result<(), LinAlgError> { + let n = b.len(); + for i in 0..n { + let mut max_a = 0.0; + let mut imax = i; + for k in i..n { + let abs_a = a[(k, i)].re().abs(); + if abs_a > max_a { + max_a = abs_a; + imax = k; + } + } + if max_a == 0.0 { + return Err(LinAlgError()); + } + if imax != i { + a.swap_rows(i, imax); + b.swap_rows(i, imax); + } + for j in i + 1..n { + a[(j, i)] = a[(j, i)] / a[(i, i)]; + for k in i + 1..n { + a[(j, k)] = a[(j, k)] - a[(j, i)] * a[(i, k)]; + } + } + } + for i in 0..n { + for k in 0..i { + b[i] = b[i] - a[(i, k)] * b[k]; + } + } + for i in (0..n).rev() { + for k in i + 1..n { + b[i] = b[i] - a[(i, k)] * b[k]; + } + b[i] /= a[(i, i)]; } + Ok(()) } #[cfg(test)] diff --git a/crates/feos/src/pcsaft/eos/mod.rs b/crates/feos/src/pcsaft/eos/mod.rs index 02e5c0162..e82669fac 100644 --- a/crates/feos/src/pcsaft/eos/mod.rs +++ b/crates/feos/src/pcsaft/eos/mod.rs @@ -48,9 +48,9 @@ pub struct PcSaft { params: PcSaftPars, options: PcSaftOptions, hard_chain: bool, - dipole: bool, - quadrupole: bool, - dipole_quadrupole: bool, + dipole: Option, + quadrupole: Option, + dipole_quadrupole: Option, association: Option, } @@ -63,9 +63,10 @@ impl PcSaft { let params = PcSaftPars::new(¶meters); let hard_chain = params.m.iter().any(|m| (m - 1.0).abs() > 1e-15); - let dipole = params.ndipole > 0; - let quadrupole = params.nquadpole > 0; - let dipole_quadrupole = params.ndipole > 0 && params.nquadpole > 0; + let dipole = (params.ndipole > 0).then(|| Dipole::new(¶ms)); + let quadrupole = (params.nquadpole > 0).then(|| Quadrupole::new(¶ms)); + let dipole_quadrupole = (params.ndipole > 0 && params.nquadpole > 0) + .then(|| DipoleQuadrupole::new(¶ms, options.dq_variant)); let association = (!parameters.association.is_empty()) .then(|| Association::new(options.max_iter_cross_assoc, options.tol_cross_assoc)); @@ -137,26 +138,22 @@ impl PcSaft { "Dispersion", Dispersion.helmholtz_energy_density(&self.params, state, zeta[3]), ); - if self.dipole { + if let Some(dipole) = &self.dipole { f( "Dipole", - Dipole.helmholtz_energy_density(&self.params, state), + dipole.helmholtz_energy_density(&self.params, state, zeta[3]), ) } - if self.quadrupole { + if let Some(quadrupole) = &self.quadrupole { f( "Quadrupole", - Quadrupole.helmholtz_energy_density(&self.params, state), + quadrupole.helmholtz_energy_density(&self.params, state, zeta[3]), ) } - if self.dipole_quadrupole { + if let Some(dipole_quadrupole) = &self.dipole_quadrupole { f( "DipoleQuadrupole", - DipoleQuadrupole.helmholtz_energy_density( - &self.params, - state, - self.options.dq_variant, - ), + dipole_quadrupole.helmholtz_energy_density(&self.params, state, zeta[3]), ) } if let Some(association) = self.association.as_ref() { diff --git a/crates/feos/src/pcsaft/eos/polar.rs b/crates/feos/src/pcsaft/eos/polar.rs index 6f23f9b61..161080a9c 100644 --- a/crates/feos/src/pcsaft/eos/polar.rs +++ b/crates/feos/src/pcsaft/eos/polar.rs @@ -1,5 +1,4 @@ use super::PcSaftPars; -use crate::hard_sphere::HardSphereProperties; use feos_core::StateHD; use nalgebra::DMatrix; use num_dual::DualNum; @@ -159,29 +158,73 @@ fn triplet_integral_ijk_dq + Copy>(mijk: f64, etas: &[D], c: &[[ .sum() } -pub struct Dipole; +/// Combinatorial factor of a pair term in the sum over $i \le j$. +fn pair_factor(i: usize, j: usize) -> f64 { + if i == j { 1.0 } else { 2.0 } +} + +/// Combinatorial factor of a triplet term in the sum over $i \le j \le k$. +fn triplet_factor(i: usize, j: usize, k: usize) -> f64 { + if i == k { + 1.0 + } else if i == j || j == k { + 3.0 + } else { + 6.0 + } +} + +/// Dipole contribution. Everything that only depends on the parameters is +/// precomputed in [Dipole::new]. +pub struct Dipole { + m: MeanSegmentNumbers, + /// $c_{ij}\tilde\mu_i^2\tilde\mu_j^2/\sigma_{ij}^3$ for $i \le j$ (times $1/T^2$ at runtime) + f2: DMatrix, + /// $c_{ijk}\tilde\mu_i^2\tilde\mu_j^2\tilde\mu_k^2/(\sigma_{ij}\sigma_{ik}\sigma_{jk})$ for $i \le j \le k$ (times $1/T^3$ at runtime) + f3: Vec>, +} impl Dipole { + pub fn new(parameters: &PcSaftPars) -> Self { + let p = parameters; + let n = p.ndipole; + let mu2: Vec = p + .dipole_comp + .iter() + .map(|&i| p.sigma[i].powi(3) * p.epsilon_k[i] * p.mu2[i]) + .collect(); + let mut f2 = DMatrix::zeros(n, n); + let mut f3 = vec![DMatrix::zeros(n, n); n]; + for i in 0..n { + let di = p.dipole_comp[i]; + for j in i..n { + let dj = p.dipole_comp[j]; + f2[(i, j)] = pair_factor(i, j) * mu2[i] * mu2[j] / p.sigma_ij[(di, dj)].powi(3); + for k in j..n { + let dk = p.dipole_comp[k]; + f3[i][(j, k)] = triplet_factor(i, j, k) * mu2[i] * mu2[j] * mu2[k] + / (p.sigma_ij[(di, dj)] * p.sigma_ij[(di, dk)] * p.sigma_ij[(dj, dk)]); + } + } + } + Self { + m: MeanSegmentNumbers::new(p, Multipole::Dipole), + f2, + f3, + } + } + + /// `eta` is the packing fraction $\zeta_3$, shared with the other contributions. #[inline] pub fn helmholtz_energy_density + Copy>( &self, parameters: &PcSaftPars, state: &StateHD, + eta: D, ) -> D { - let m = MeanSegmentNumbers::new(parameters, Multipole::Dipole); let p = parameters; - let t_inv = state.temperature.inv(); - let eps_ij_t = p.e_k_ij.map(|v| t_inv * v); - let sig_ij_3 = p.sigma_ij.map(|v| v.powi(3)); - let mu2_term: Vec = p - .dipole_comp - .iter() - .map(|&i| t_inv * sig_ij_3[(i, i)] * p.epsilon_k[i] * p.mu2[i]) - .collect(); - let rho = &state.partial_density; - let [eta] = parameters.zeta(state.temperature, rho, [3]); let eta2 = eta * eta; let etas = [D::one(), eta, eta2, eta2 * eta, eta2 * eta2]; @@ -191,39 +234,35 @@ impl Dipole { let di = p.dipole_comp[i]; for j in i..p.ndipole { let dj = p.dipole_comp[j]; - let c = if i == j { 1.0 } else { 2.0 }; phi2 -= rho[di] * rho[dj] - * mu2_term[i] - * mu2_term[j] + * self.f2[(i, j)] * pair_integral_ij( - m.mij1[(i, j)], - m.mij2[(i, j)], + self.m.mij1[(i, j)], + self.m.mij2[(i, j)], &etas, &AD, &BD, - eps_ij_t[(di, dj)], - ) - / sig_ij_3[(di, dj)] - * c; + t_inv * p.e_k_ij[(di, dj)], + ); for k in j..p.ndipole { let dk = p.dipole_comp[k]; - let c = if i == k { - 1.0 - } else if i == j || j == k { - 3.0 - } else { - 6.0 - }; - phi3 -= rho[di] * rho[dj] * rho[dk] * mu2_term[i] * mu2_term[j] * mu2_term[k] - / (p.sigma_ij[(di, dj)] * p.sigma_ij[(di, dk)] * p.sigma_ij[(dj, dk)]) - * triplet_integral_ijk(m.mijk1[i][(j, k)], m.mijk2[i][(j, k)], &etas, &CD) - * c; + phi3 -= rho[di] + * rho[dj] + * rho[dk] + * self.f3[i][(j, k)] + * triplet_integral_ijk( + self.m.mijk1[i][(j, k)], + self.m.mijk2[i][(j, k)], + &etas, + &CD, + ); } } } - phi2 *= PI; - phi3 *= PI_SQ_43; + let t_inv2 = t_inv * t_inv; + phi2 *= t_inv2 * PI; + phi3 *= t_inv2 * t_inv * PI_SQ_43; let mut result = phi2 * phi2 / (phi2 - phi3); if result.re().is_nan() { result = phi2 @@ -232,29 +271,58 @@ impl Dipole { } } -pub struct Quadrupole; +/// Quadrupole contribution. Everything that only depends on the parameters is +/// precomputed in [Quadrupole::new]. +pub struct Quadrupole { + m: MeanSegmentNumbers, + /// $c_{ij}\tilde Q_i^2\tilde Q_j^2/\sigma_{ij}^7$ for $i \le j$ (times $1/T^2$ at runtime) + f2: DMatrix, + /// $c_{ijk}\tilde Q_i^2\tilde Q_j^2\tilde Q_k^2/(\sigma_{ij}\sigma_{ik}\sigma_{jk})^3$ for $i \le j \le k$ (times $1/T^3$ at runtime) + f3: Vec>, +} impl Quadrupole { + pub fn new(parameters: &PcSaftPars) -> Self { + let p = parameters; + let n = p.nquadpole; + let q2: Vec = p + .quadpole_comp + .iter() + .map(|&i| p.sigma[i].powi(5) * p.epsilon_k[i] * p.q2[i]) + .collect(); + let mut f2 = DMatrix::zeros(n, n); + let mut f3 = vec![DMatrix::zeros(n, n); n]; + for i in 0..n { + let di = p.quadpole_comp[i]; + for j in i..n { + let dj = p.quadpole_comp[j]; + f2[(i, j)] = pair_factor(i, j) * q2[i] * q2[j] / p.sigma_ij[(di, dj)].powi(7); + for k in j..n { + let dk = p.quadpole_comp[k]; + f3[i][(j, k)] = triplet_factor(i, j, k) * q2[i] * q2[j] * q2[k] + / (p.sigma_ij[(di, dj)] * p.sigma_ij[(di, dk)] * p.sigma_ij[(dj, dk)]) + .powi(3); + } + } + } + Self { + m: MeanSegmentNumbers::new(p, Multipole::Quadrupole), + f2, + f3, + } + } + + /// `eta` is the packing fraction $\zeta_3$, shared with the other contributions. #[inline] pub fn helmholtz_energy_density + Copy>( &self, parameters: &PcSaftPars, state: &StateHD, + eta: D, ) -> D { - let m = MeanSegmentNumbers::new(parameters, Multipole::Quadrupole); let p = parameters; - let t_inv = state.temperature.inv(); - let eps_ij_t = p.e_k_ij.map(|v| t_inv * v); - let sig_ij_3 = p.sigma_ij.map(|v| v.powi(3)); - let q2_term: Vec = p - .quadpole_comp - .iter() - .map(|&i| t_inv * p.sigma[i].powi(5) * p.epsilon_k[i] * p.q2[i]) - .collect(); - let rho = &state.partial_density; - let [eta] = parameters.zeta(state.temperature, rho, [3]); let eta2 = eta * eta; let etas = [D::one(), eta, eta2, eta2 * eta, eta2 * eta2]; @@ -264,40 +332,35 @@ impl Quadrupole { let di = p.quadpole_comp[i]; for j in i..p.nquadpole { let dj = p.quadpole_comp[j]; - let c = if i == j { 1.0 } else { 2.0 }; - phi2 -= (rho[di] + phi2 -= rho[di] * rho[dj] - * q2_term[i] - * q2_term[j] + * self.f2[(i, j)] * pair_integral_ij( - m.mij1[(i, j)], - m.mij2[(i, j)], + self.m.mij1[(i, j)], + self.m.mij2[(i, j)], &etas, &AQ, &BQ, - eps_ij_t[(di, dj)], - )) - / p.sigma_ij[(di, dj)].powi(7) - * c; + t_inv * p.e_k_ij[(di, dj)], + ); for k in j..p.nquadpole { let dk = p.quadpole_comp[k]; - let c = if i == k { - 1.0 - } else if i == j || j == k { - 3.0 - } else { - 6.0 - }; - phi3 += rho[di] * rho[dj] * rho[dk] * q2_term[i] * q2_term[j] * q2_term[k] - / (sig_ij_3[(di, dj)] * sig_ij_3[(di, dk)] * sig_ij_3[(dj, dk)]) - * triplet_integral_ijk(m.mijk1[i][(j, k)], m.mijk2[i][(j, k)], &etas, &CQ) - * c; + phi3 += rho[di] + * rho[dj] + * rho[dk] + * self.f3[i][(j, k)] + * triplet_integral_ijk( + self.m.mijk1[i][(j, k)], + self.m.mijk2[i][(j, k)], + &etas, + &CQ, + ); } } } - - phi2 *= PI * 0.5625; - phi3 *= PI * PI * 0.5625; + let t_inv2 = t_inv * t_inv; + phi2 *= t_inv2 * (PI * 0.5625); + phi3 *= t_inv2 * t_inv * (PI * PI * 0.5625); let mut result = phi2 * phi2 / (phi2 - phi3); if result.re().is_nan() { result = phi2 @@ -313,105 +376,143 @@ pub enum DQVariants { DQ44, } -pub struct DipoleQuadrupole; +/// Dipole-quadrupole cross contribution. Everything that only depends on the +/// parameters (and the combination rule) is precomputed in [DipoleQuadrupole::new]. +pub struct DipoleQuadrupole { + mdq1: DMatrix, + mdq2: DMatrix, + mdqd: Vec>, + mdqq: Vec>, + /// pair prefactor for dipole $i$ and quadrupole $j$ (times $1/T^2$ at runtime) + f2: DMatrix, + /// triplet prefactor for dipole $i$, quadrupole $j$, dipole $k$ (times $1/T^3$ at runtime) + f3d: Vec>, + /// triplet prefactor for dipole $i$, quadrupoles $j$ and $k$, including $\alpha$ (times $1/T^3$ at runtime) + f3q: Vec>, +} impl DipoleQuadrupole { - #[inline] - pub fn helmholtz_energy_density + Copy>( - &self, - parameters: &PcSaftPars, - state: &StateHD, - variant: DQVariants, - ) -> D { + pub fn new(parameters: &PcSaftPars, variant: DQVariants) -> Self { let p = parameters; - - let t_inv = state.temperature.inv(); - let eps_ij_t = p.e_k_ij.map(|v| t_inv * v); - - let q2_term: Vec = p + let (nd, nq) = (p.ndipole, p.nquadpole); + let q2: Vec = p .quadpole_comp .iter() - .map(|&i| t_inv * p.sigma[i].powi(4) * p.epsilon_k[i] * p.q2[i]) + .map(|&i| p.sigma[i].powi(4) * p.epsilon_k[i] * p.q2[i]) .collect(); - let mu2_term: Vec = p + let mu2: Vec = p .dipole_comp .iter() - .map(|&i| t_inv * p.sigma[i].powi(4) * p.epsilon_k[i] * p.mu2[i]) + .map(|&i| p.sigma[i].powi(4) * p.epsilon_k[i] * p.mu2[i]) .collect(); - let rho = &state.partial_density; - let [eta] = parameters.zeta(state.temperature, rho, [3]); - let eta2 = eta * eta; - let etas = [D::one(), eta, eta2, eta2 * eta, eta2 * eta2]; - - // mean segment number - let mut mdq1 = DMatrix::zeros(p.ndipole, p.nquadpole); - let mut mdq2 = DMatrix::zeros(p.ndipole, p.nquadpole); - let mut mdqd = vec![DMatrix::zeros(p.nquadpole, p.ndipole); p.ndipole]; - let mut mdqq = vec![DMatrix::zeros(p.nquadpole, p.nquadpole); p.ndipole]; - for i in 0..p.ndipole { + // mean segment numbers and prefactors + let mut mdq1 = DMatrix::zeros(nd, nq); + let mut mdq2 = DMatrix::zeros(nd, nq); + let mut mdqd = vec![DMatrix::zeros(nq, nd); nd]; + let mut mdqq = vec![DMatrix::zeros(nq, nq); nd]; + let mut f2 = DMatrix::zeros(nd, nq); + let mut f3d = vec![DMatrix::zeros(nq, nd); nd]; + let mut f3q = vec![DMatrix::zeros(nq, nq); nd]; + for i in 0..nd { let di = p.dipole_comp[i]; let mi = p.m[di].min(2.0); - for j in 0..p.nquadpole { + for j in 0..nq { let qj = p.quadpole_comp[j]; let mj = p.m[qj].min(2.0); let m = (mi * mj).sqrt(); mdq1[(i, j)] = (m - 1.0) / m; mdq2[(i, j)] = mdq1[(i, j)] * (m - 2.0) / m; - for k in 0..p.ndipole { + let mu2_q2 = match variant { + DQVariants::DQ35 => mu2[i] / p.sigma[di] * q2[j] * p.sigma[qj], + DQVariants::DQ44 => mu2[i] * q2[j], + }; + f2[(i, j)] = mu2_q2 / p.sigma_ij[(di, qj)].powi(5); + for k in 0..nd { let dk = p.dipole_comp[k]; let mk = p.m[dk].min(2.0); let m = (mi * mj * mk).cbrt(); mdqd[i][(j, k)] = (m - 1.0) / m; + f3d[i][(j, k)] = mu2[i] * q2[j] * mu2[k] + / (p.sigma_ij[(di, qj)] * p.sigma_ij[(di, dk)] * p.sigma_ij[(qj, dk)]) + .powi(2); } - for k in 0..p.nquadpole { + for k in 0..nq { let qk = p.quadpole_comp[k]; let mk = p.m[qk].min(2.0); let m = (mi * mj * mk).cbrt(); mdqq[i][(j, k)] = (m - 1.0) / m; + f3q[i][(j, k)] = mu2[i] * q2[j] * q2[k] + / (p.sigma_ij[(di, qj)] * p.sigma_ij[(di, qk)] * p.sigma_ij[(qj, qk)]) + .powi(2) + * ALPHA; } } } + Self { + mdq1, + mdq2, + mdqd, + mdqq, + f2, + f3d, + f3q, + } + } + + /// `eta` is the packing fraction $\zeta_3$, shared with the other contributions. + #[inline] + pub fn helmholtz_energy_density + Copy>( + &self, + parameters: &PcSaftPars, + state: &StateHD, + eta: D, + ) -> D { + let p = parameters; + let t_inv = state.temperature.inv(); + let rho = &state.partial_density; + let eta2 = eta * eta; + let etas = [D::one(), eta, eta2, eta2 * eta, eta2 * eta2]; let mut phi2 = D::zero(); let mut phi3 = D::zero(); for i in 0..p.ndipole { + let di = p.dipole_comp[i]; for j in 0..p.nquadpole { - let di = p.dipole_comp[i]; let qj = p.quadpole_comp[j]; - let mu2_q2_term = match variant { - DQVariants::DQ35 => mu2_term[i] / p.sigma[di] * q2_term[j] * p.sigma[qj], - DQVariants::DQ44 => mu2_term[i] * q2_term[j], - }; - phi2 -= rho[di] * rho[qj] * mu2_q2_term / p.sigma_ij[(di, qj)].powi(5) + phi2 -= rho[di] + * rho[qj] + * self.f2[(i, j)] * pair_integral_ij( - mdq1[(i, j)], - mdq2[(i, j)], + self.mdq1[(i, j)], + self.mdq2[(i, j)], &etas, &ADQ, &BDQ, - eps_ij_t[(di, qj)], + t_inv * p.e_k_ij[(di, qj)], ); for k in 0..p.ndipole { let dk = p.dipole_comp[k]; - phi3 += rho[di] * rho[qj] * rho[dk] * mu2_term[i] * q2_term[j] * mu2_term[k] - / (p.sigma_ij[(di, qj)] * p.sigma_ij[(di, dk)] * p.sigma_ij[(qj, dk)]) - .powi(2) - * triplet_integral_ijk_dq(mdqd[i][(j, k)], &etas, &CDQ); + phi3 += rho[di] + * rho[qj] + * rho[dk] + * self.f3d[i][(j, k)] + * triplet_integral_ijk_dq(self.mdqd[i][(j, k)], &etas, &CDQ); } for k in 0..p.nquadpole { let qk = p.quadpole_comp[k]; - phi3 += rho[di] * rho[qj] * rho[qk] * mu2_term[i] * q2_term[j] * q2_term[k] - / (p.sigma_ij[(di, qj)] * p.sigma_ij[(di, qk)] * p.sigma_ij[(qj, qk)]) - .powi(2) - * ALPHA - * triplet_integral_ijk_dq(mdqq[i][(j, k)], &etas, &CDQ); + phi3 += rho[di] + * rho[qj] + * rho[qk] + * self.f3q[i][(j, k)] + * triplet_integral_ijk_dq(self.mdqq[i][(j, k)], &etas, &CDQ); } } } - phi2 *= PI * 2.25; - phi3 *= PI * PI; + let t_inv2 = t_inv * t_inv; + phi2 *= t_inv2 * (PI * 2.25); + phi3 *= t_inv2 * t_inv * (PI * PI); let mut result = phi2 * phi2 / (phi2 - phi3); if result.re().is_nan() { result = phi2 @@ -423,6 +524,7 @@ impl DipoleQuadrupole { #[cfg(test)] mod tests { use super::*; + use crate::hard_sphere::HardSphere; use crate::pcsaft::PcSaftParameters; use crate::pcsaft::eos::dispersion::tests::a_disp; use crate::pcsaft::parameters::utils::{ @@ -434,6 +536,18 @@ mod tests { use itertools::Itertools; use nalgebra::{DVector, dvector}; + fn eta(p: &PcSaftPars, s: &StateHD) -> f64 { + HardSphere.helmholtz_energy_density_and_properties(p, s).1[3] + } + + fn a_dipole(p: &PcSaftPars, s: &StateHD) -> f64 { + Dipole::new(p).helmholtz_energy_density(p, s, eta(p, s)) + } + + fn a_quadrupole(p: &PcSaftPars, s: &StateHD) -> f64 { + Quadrupole::new(p).helmholtz_energy_density(p, s, eta(p, s)) + } + #[test] fn test_dipolar_contribution() { let dme = dme_parameters(); @@ -441,7 +555,7 @@ mod tests { let v = 1000.0; let n = 1.0; let s = StateHD::new(t, v, &dvector![n]); - let a = Dipole.helmholtz_energy_density(&dme, &s) * v; + let a = a_dipole(&dme, &s) * v; assert_relative_eq!(a, -1.40501033595417E-002, epsilon = 1e-6); } @@ -459,7 +573,7 @@ mod tests { let t = 350.0; let v = 1000.0; let s = StateHD::new(t, v, &dvector![1.0, 2.0, 3.0]); - let a = Dipole.helmholtz_energy_density(¶meters, &s) * v; + let a = a_dipole(¶meters, &s) * v; assert_relative_eq!(a, -1.4126308106201688, epsilon = 1e-10); } @@ -482,7 +596,7 @@ mod tests { let t = 350.0; let v = 1000.0; let s = StateHD::new(t, v, &DVector::from(n)); - let a = Dipole.helmholtz_energy_density(¶meters, &s) * v; + let a = a_dipole(¶meters, &s) * v; println!("{components:?}: {a}"); a }) @@ -499,7 +613,7 @@ mod tests { let v = 1000.0; let n = 1.0; let s = StateHD::new(t, v, &dvector![n]); - let a = Quadrupole.helmholtz_energy_density(&co2, &s) * v; + let a = a_quadrupole(&co2, &s) * v; assert_relative_eq!(a, -4.38559558854186E-002, epsilon = 1e-6); } @@ -517,7 +631,7 @@ mod tests { let t = 350.0; let v = 1000.0; let s = StateHD::new(t, v, &dvector![1.0, 2.0, 3.0]); - let a = Quadrupole.helmholtz_energy_density(¶meters, &s) * v; + let a = a_quadrupole(¶meters, &s) * v; assert_relative_eq!(a, -0.2748017468669352, epsilon = 1e-10); } @@ -540,7 +654,7 @@ mod tests { let t = 350.0; let v = 1000.0; let s = StateHD::new(t, v, &DVector::from(n)); - let a = Quadrupole.helmholtz_energy_density(¶meters, &s) * v; + let a = a_quadrupole(¶meters, &s) * v; println!("{components:?}: {a}"); a }) @@ -565,8 +679,8 @@ mod tests { let s = StateHD::new(t, v / n, &(dvector![n_dme, n_co2] / n)); // let a_dpqp = dpqp.helmholtz_energy_density(&s)*v; let a_disp = a_disp(&mix, &s) * v; - let a_qp = Quadrupole.helmholtz_energy_density(&mix, &s) * v; - let a_dp = Dipole.helmholtz_energy_density(&mix, &s) * v; + let a_qp = a_quadrupole(&mix, &s) * v; + let a_dp = a_dipole(&mix, &s) * v; assert_relative_eq!(a_disp, -1.6283622072860044, epsilon = 1e-6); assert_relative_eq!(a_dp, -1.35361827881345E-002, epsilon = 1e-6); assert_relative_eq!(a_qp, -4.20168059082731E-002, epsilon = 1e-6); From a82c070ccb6352c2b1555ec39d8951afa099ab62 Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Fri, 16 Oct 2026 15:22:44 +0200 Subject: [PATCH 3/3] Fix benchmark bug. --- crates/feos/benches/contributions.rs | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/crates/feos/benches/contributions.rs b/crates/feos/benches/contributions.rs index 7dc8332c2..26f0b1a5a 100644 --- a/crates/feos/benches/contributions.rs +++ b/crates/feos/benches/contributions.rs @@ -74,13 +74,14 @@ fn pcsaft(c: &mut Criterion) { State::new_npt(&&eos, t, p, &moles, Some(DensityInitialization::Liquid)).unwrap(); let temperature = Dual64::from(state.temperature.into_reduced()).derivative(); let molar_volume = Dual::from(1.0 / state.density.into_reduced()); - let moles = state.moles().unwrap().to_reduced().map(Dual::from); + // one particle in the volume per particle, i.e., at the state's density + let molefracs = state.molefracs.map(Dual::from); // let state_hd = state.derive1(Derivative::DT); let name1 = comp1.identifier.name.as_deref().unwrap(); let name2 = comp2.identifier.name.as_deref().unwrap(); let mix = format!("{name1}_{name2}"); group.bench_function(mix, |b| { - b.iter(|| (&eos).residual_helmholtz_energy(temperature, molar_volume, &moles)) + b.iter(|| (&eos).residual_helmholtz_energy(temperature, molar_volume, &molefracs)) }); } }