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/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)) }); } } 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/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..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)); @@ -101,46 +102,62 @@ 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), - )); - if self.dipole { - v.push(( + Dispersion.helmholtz_energy_density(&self.params, state, zeta[3]), + ); + 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 { - v.push(( + 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 { - v.push(( + 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() { - v.push(( + f( "Association", association.helmholtz_energy_density( &self.params, @@ -148,9 +165,8 @@ impl ResidualDyn for PcSaft { state, &d, ), - )) + ) } - v } } @@ -168,7 +184,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..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,8 +524,9 @@ impl DipoleQuadrupole { #[cfg(test)] mod tests { use super::*; + use crate::hard_sphere::HardSphere; 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, }; @@ -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 }) @@ -564,9 +678,9 @@ 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_qp = Quadrupole.helmholtz_energy_density(&mix, &s) * v; - let a_dp = Dipole.helmholtz_energy_density(&mix, &s) * v; + let a_disp = a_disp(&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); 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 {