From 0713dbeaf14b1e12a34e4d9131eb0e28051efd2f Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Mon, 31 Aug 2026 17:24:31 +0200 Subject: [PATCH] Extend Newton solver to DFT specifications other than chemical potentials --- crates/feos-dft/src/profile/mod.rs | 59 ++++++++++++++--------- crates/feos-dft/src/profile/properties.rs | 2 +- crates/feos-dft/src/solver.rs | 23 ++++++--- 3 files changed, 55 insertions(+), 29 deletions(-) diff --git a/crates/feos-dft/src/profile/mod.rs b/crates/feos-dft/src/profile/mod.rs index 6def8ded2..a9052174a 100644 --- a/crates/feos-dft/src/profile/mod.rs +++ b/crates/feos-dft/src/profile/mod.rs @@ -34,14 +34,24 @@ pub enum DFTSpecification { } impl DFTSpecification { - fn calculate_fugacity(&self, z: &Array1) -> FeosResult> { - Ok(match self { + fn calculate_fugacity(&self, z: &Array1) -> Array1 { + match self { Self::ChemicalPotential(fugacity) => fugacity.clone(), Self::Moles(moles) => moles / z, Self::TotalMoles(total_moles, fugacity) => { fugacity * *total_moles / (fugacity * z).sum() } - }) + } + } + + pub(crate) fn delta_fugacity(&self, z: &Array1, delta_z: &Array1) -> Array1 { + match self { + Self::ChemicalPotential(fugacity) => Array1::zeros(fugacity.len()), + Self::Moles(_) => -delta_z / z, + Self::TotalMoles(_, fugacity) => { + -(fugacity * delta_z).sum() / (fugacity * z).sum() * Array1::ones(fugacity.len()) + } + } } pub fn from_state(state: &State) -> Self { @@ -216,7 +226,7 @@ where profile.sum() * functional_determinant } - fn integrate_reduced_comp, N: DualNum + Copy>( + pub(crate) fn integrate_reduced_comp, N: DualNum + Copy>( &self, profile: &ArrayBase, ) -> Array1 { @@ -320,7 +330,7 @@ where pub fn residual(&self, log: bool) -> FeosResult<(Array, f64)> { let density = self.density.to_reduced(); - let (res, res_norm, _, _) = self.euler_lagrange_equation(&density, log)?; + let (res, res_norm, _, _, _) = self.euler_lagrange_equation(&density, log)?; Ok((res, res_norm)) } @@ -328,7 +338,12 @@ where fn fugacity( &self, density: &Array, - ) -> FeosResult<(Array, Array, Array1)> { + ) -> FeosResult<( + Array, + Array1, + Array, + Array1, + )> { // calculate reduced temperature let temperature = self.temperature.to_reduced(); @@ -352,14 +367,21 @@ where .bulk .eos .bond_integrals(temperature, &exp_dfdrho, self.convolver.as_ref()); - let z = &exp_dfdrho * bonds; + let mut rho_projected = &exp_dfdrho * bonds; + let z = self.integrate_reduced_comp(&rho_projected); // calculate fugacity based on the given specification - let fugacity = self - .specification - .calculate_fugacity(&self.integrate_reduced_comp(&z))?; + let fugacity = self.specification.calculate_fugacity(&z); + + // multiply fugacity + rho_projected + .outer_iter_mut() + .zip(fugacity.iter()) + .for_each(|(mut x, &f)| { + x *= f; + }); - Ok((exp_dfdrho, z, fugacity)) + Ok((exp_dfdrho, z, rho_projected, fugacity)) } #[expect(clippy::type_complexity)] @@ -371,18 +393,11 @@ where Array, f64, Array, + Array1, Array, )> { // calculate functional derivatives and fugacity - let (exp_dfdrho, mut rho_projected, fugacity) = self.fugacity(density)?; - - // multiply fugacity - rho_projected - .outer_iter_mut() - .zip(fugacity.iter()) - .for_each(|(mut x, &f)| { - x *= f; - }); + let (exp_dfdrho, z, rho_projected, _) = self.fugacity(density)?; // calculate residual let mut res = if log { @@ -402,7 +417,7 @@ where (density - &rho_projected).mapv(|x| x * x).sum().sqrt() / (res.len() as f64).sqrt(); if res_norm.is_finite() { - Ok((res, res_norm, exp_dfdrho, rho_projected)) + Ok((res, res_norm, exp_dfdrho, z, rho_projected)) } else { Err(FeosError::IterationFailed("Euler-Lagrange equation".into())) } @@ -423,7 +438,7 @@ where // solve a bulk profile with the Newton solver let mut bulk_profile = DFTProfile::::new(Grid::Bulk, &self.bulk, None, None, None); - let (_, _, fugacity) = self.fugacity(&density)?; + let (_, _, _, fugacity) = self.fugacity(&density)?; bulk_profile.specification = DFTSpecification::ChemicalPotential(fugacity); let solver = DFTSolver::new(None).newton(None, None, None, None); bulk_profile.solve(Some(&solver), false)?; diff --git a/crates/feos-dft/src/profile/properties.rs b/crates/feos-dft/src/profile/properties.rs index 8e602b34b..4b7516853 100644 --- a/crates/feos-dft/src/profile/properties.rs +++ b/crates/feos-dft/src/profile/properties.rs @@ -299,7 +299,7 @@ where fn density_derivative(&self, lhs: &Array) -> FeosResult> { let rho = self.density.to_reduced(); let second_partial_derivatives = self.second_partial_derivatives(&rho)?; - let (_, _, exp_dfdrho, _) = self.euler_lagrange_equation(&rho, false)?; + let (_, _, exp_dfdrho, _, _) = self.euler_lagrange_equation(&rho, false)?; let rhs = |x: &_| { let delta_functional_derivative = diff --git a/crates/feos-dft/src/solver.rs b/crates/feos-dft/src/solver.rs index 2baba60fb..15baf742f 100644 --- a/crates/feos-dft/src/solver.rs +++ b/crates/feos-dft/src/solver.rs @@ -262,7 +262,7 @@ where for k in 0..picard.max_iter { // calculate residual - let (res, res_norm, _, _) = self.euler_lagrange_equation(&*rho, picard.log)?; + let (res, res_norm, _, _, _) = self.euler_lagrange_equation(&*rho, picard.log)?; log.add_residual(solver, k, res_norm); // check for convergence @@ -304,7 +304,7 @@ where } else { rho + alpha * delta_rho }; - let Ok((_, res2, _, _)) = self.euler_lagrange_equation(&rho_new, logarithm) else { + let Ok((_, res2, _, _, _)) = self.euler_lagrange_equation(&rho_new, logarithm) else { continue; }; if res2 > res0 { @@ -317,7 +317,7 @@ where } else { rho + 0.5 * alpha * delta_rho }; - let Ok((_, res1, _, _)) = self.euler_lagrange_equation(&rho_new, logarithm) else { + let Ok((_, res1, _, _, _)) = self.euler_lagrange_equation(&rho_new, logarithm) else { continue; }; @@ -370,7 +370,7 @@ where let m = resm.len() + 1; // calculate residual - let (res, res_norm, _, _) = self.euler_lagrange_equation(&*rho, anderson.log)?; + let (res, res_norm, _, _, _) = self.euler_lagrange_equation(&*rho, anderson.log)?; log.add_residual(solver, k, res_norm); // check for convergence @@ -426,7 +426,7 @@ where let solver = if newton.log { "Newton (log)" } else { "Newton" }; for k in 0..newton.max_iter { // calculate initial residual - let (res, res_norm, exp_dfdrho, rho_p) = + let (res, res_norm, exp_dfdrho, z, rho_p) = self.euler_lagrange_equation(rho, newton.log)?; log.add_residual(solver, k, res_norm); @@ -447,8 +447,19 @@ where .zip(self.bulk.eos.m().iter()) .for_each(|(mut q, &m)| q /= m); let delta_i = self.delta_bond_integrals(&exp_dfdrho, &delta_functional_derivative); + let mut delta_exp_dfdrho = delta_functional_derivative - delta_i; + let delta_z = -self.integrate_reduced_comp(&(&delta_exp_dfdrho * &exp_dfdrho)); + + let delta_fugacity = self.specification.delta_fugacity(&z, &delta_z); + delta_exp_dfdrho + .outer_iter_mut() + .zip(delta_fugacity.iter()) + .for_each(|(mut z, &f)| { + z -= f; + }); + let rho = if newton.log { &*rho } else { &rho_p }; - delta_rho + (delta_functional_derivative - delta_i) * rho + delta_rho + delta_exp_dfdrho * rho }; // update solution