From b82c32f3f7c79bbd156a0c5f0ce345d73306a6b3 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Fri, 27 Jan 2023 19:10:20 +0100 Subject: [PATCH 1/3] Fix association contribution for some association schemes and improve performance --- src/association/mod.rs | 54 ++++++++++++++++++++++++++++-------------- 1 file changed, 36 insertions(+), 18 deletions(-) diff --git a/src/association/mod.rs b/src/association/mod.rs index cc53cc8ba..aa1b74bb7 100644 --- a/src/association/mod.rs +++ b/src/association/mod.rs @@ -7,6 +7,7 @@ use num_dual::linalg::{norm, LU}; use num_dual::*; use serde::{Deserialize, Serialize}; use std::fmt; +use std::ops::SubAssign; use std::sync::Arc; #[cfg(feature = "dft")] @@ -243,7 +244,7 @@ impl Association

{ (((deltarho * (na - nb) + 1.0).powi(2) + deltarho * nb * 4.0).sqrt() + (deltarho * (nb - na) + 1.0)) .recip() - * (2.0 * nb / na) + * 2.0 } pub fn assoc_site_frac_a>(deltarho: D, na: f64) -> D { @@ -336,37 +337,32 @@ impl Association

{ tol: f64, ) -> Result { // gradient - let mut g: Array1 = Array::zeros(2 * nassoc); + let mut g = x.map(D::recip); // Hessian let mut h: Array2 = Array::zeros((2 * nassoc, 2 * nassoc)); - // slice arrays - let (xa, xb) = x.multi_slice_mut((s![..nassoc], s![nassoc..])); - let (mut ga, mut gb) = g.multi_slice_mut((s![..nassoc], s![nassoc..])); - let (mut haa, mut hab, mut hba, mut hbb) = h.multi_slice_mut(( - s![..nassoc, ..nassoc], - s![..nassoc, nassoc..], - s![nassoc.., ..nassoc], - s![nassoc.., nassoc..], - )); + // split x array + let (xa, xb) = x.view().split_at(Axis(0), nassoc); // calculate gradients and approximate Hessian for i in 0..nassoc { let d = &delta.index_axis(Axis(0), i) * rho; let dnx = (&xb * nb * &d).sum() + 1.0; - ga[i] = xa[i].recip() - dnx; - hab.index_axis_mut(Axis(0), i).assign(&(&d * &(-nb))); - haa[(i, i)] = -dnx / xa[i]; + g[i] -= dnx; + for j in 0..nassoc { + h[(i, nassoc + j)] = -d[j] * nb[j]; + h[(nassoc + i, j)] = -d[j] * na[j]; + } + h[(i, i)] = -dnx / xa[i]; let dnx = (&xa * na * &d).sum() + 1.0; - gb[i] = xb[i].recip() - dnx; - hba.index_axis_mut(Axis(0), i).assign(&(&d * &(-na))); - hbb[(i, i)] = -dnx / xb[i]; + g[nassoc + i] -= dnx; + h[(nassoc + i, nassoc + i)] = -dnx / xb[i]; } // Newton step - x.assign(&(&*x - &LU::new(h)?.solve(&g))); + x.sub_assign(&LU::new(h)?.solve(&g)); // check convergence Ok(norm(&g.map(D::re)) < tol) @@ -378,7 +374,9 @@ impl Association

{ mod tests_pcsaft { use super::*; use crate::pcsaft::parameters::utils::water_parameters; + use crate::pcsaft::PcSaftParameters; use approx::assert_relative_eq; + use feos_core::parameter::Parameter; #[test] fn helmholtz_energy() { @@ -403,6 +401,26 @@ mod tests_pcsaft { let a_rust = assoc.helmholtz_energy(&s) / n; assert_relative_eq!(a_rust, -4.229878997054543, epsilon = 1e-10); } + + #[test] + fn helmholtz_energy_cross_3b() { + let mut params = water_parameters(); + let mut record = params.pure_records.pop().unwrap(); + let mut association_record = record.model_record.association_record.unwrap(); + association_record.na = Some(2.0); + record.model_record.association_record = Some(association_record); + let params = Arc::new(PcSaftParameters::new_pure(record)); + let assoc = Association::new(¶ms, ¶ms.association, 50, 1e-10); + let cross_assoc = + Association::new_cross_association(¶ms, ¶ms.association, 50, 1e-10); + let t = 350.0; + let v = 41.248289328513216; + let n = 1.23; + let s = StateHD::new(t, v, arr1(&[n])); + let a_assoc = assoc.helmholtz_energy(&s) / n; + let a_cross_assoc = cross_assoc.helmholtz_energy(&s) / n; + assert_relative_eq!(a_assoc, a_cross_assoc, epsilon = 1e-10); + } } #[cfg(test)] From 6c518545115c47e66f7d7bd13b776cbc9ff99610 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Fri, 27 Jan 2023 19:13:19 +0100 Subject: [PATCH 2/3] add changelog entry --- CHANGELOG.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 75d7d0507..8f9a99f18 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,11 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [Unreleased] +### Changed +- Replaced some slow array operations to make calculations with multiple associating molecules significantly faster. [#129](https://github.com/feos-org/feos/pull/129) + +### Fixed +- Fixed a regression introduced in [#108](https://github.com/feos-org/feos/pull/108) that lead to incorrect results for the 3B association scheme. [#129](https://github.com/feos-org/feos/pull/129) ## [0.4.0] - 2023-01-27 ### Added From 090d5a30cee2b4e5878aeb9bb0f271d6c4a7d3fb Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Sat, 28 Jan 2023 14:54:11 +0100 Subject: [PATCH 3/3] update version --- CHANGELOG.md | 2 ++ Cargo.toml | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8f9a99f18..980ac7246 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,8 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [Unreleased] + +## [0.4.1] - 2023-01-28 ### Changed - Replaced some slow array operations to make calculations with multiple associating molecules significantly faster. [#129](https://github.com/feos-org/feos/pull/129) diff --git a/Cargo.toml b/Cargo.toml index fcb9f38c1..92b7a8cf0 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "feos" -version = "0.4.0" +version = "0.4.1" authors = ["Gernot Bauer ", "Philipp Rehner "] edition = "2018" readme = "README.md"