From 6a27a2755ab67d5aff5df8243126fd1772e83843 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 13 Aug 2026 20:52:47 +0200 Subject: [PATCH] Fix calculation of bulk state when moles are specified --- crates/feos-dft/src/adsorption/pore.rs | 1 + crates/feos-dft/src/convolver/mod.rs | 61 +----- .../src/convolver/periodic_convolver.rs | 16 +- crates/feos-dft/src/geometry.rs | 14 +- crates/feos-dft/src/profile/mod.rs | 174 ++++++++++-------- crates/feos-dft/src/profile/properties.rs | 17 +- crates/feos-dft/src/solver.rs | 38 ++-- py-feos/src/dft/adsorption/pore.rs | 7 + py-feos/src/dft/profile.rs | 6 +- 9 files changed, 156 insertions(+), 178 deletions(-) diff --git a/crates/feos-dft/src/adsorption/pore.rs b/crates/feos-dft/src/adsorption/pore.rs index b59474085..3a3cdd46c 100644 --- a/crates/feos-dft/src/adsorption/pore.rs +++ b/crates/feos-dft/src/adsorption/pore.rs @@ -162,6 +162,7 @@ where } pub fn update_bulk(mut self, bulk: &State) -> Self { + self.profile.specification = DFTSpecification::from_state(bulk); self.profile.bulk = bulk.clone(); self.grand_potential = None; self.interfacial_tension = None; diff --git a/crates/feos-dft/src/convolver/mod.rs b/crates/feos-dft/src/convolver/mod.rs index edca1f7e4..7a8b5298f 100644 --- a/crates/feos-dft/src/convolver/mod.rs +++ b/crates/feos-dft/src/convolver/mod.rs @@ -1,10 +1,8 @@ use crate::geometry::{Axis, Geometry, Grid}; use crate::weight_functions::*; -use ndarray::linalg::Dot; use ndarray::prelude::*; use ndarray::{Axis as Axis_nd, RemoveAxis, Slice}; use num_dual::*; -use num_traits::Zero; use rustdct::DctNum; use std::ops::{AddAssign, MulAssign, SubAssign}; use std::sync::Arc; @@ -36,46 +34,6 @@ pub trait Convolver: Send + Sync { ) -> Array; } -pub(crate) struct BulkConvolver { - weight_constants: Vec>, -} - -impl + Copy + Send + Sync> BulkConvolver { - #[expect(clippy::new_ret_no_self)] - pub(crate) fn new(weight_functions: Vec>) -> Arc> { - let weight_constants = weight_functions - .into_iter() - .map(|w| w.weight_constants(Zero::zero(), 0)) - .collect(); - Arc::new(Self { weight_constants }) - } -} - -impl + Copy + Send + Sync> Convolver for BulkConvolver -where - Array2: Dot, Output = Array1>, -{ - fn convolve(&self, _: Array0, _: &WeightFunction) -> Array0 { - unreachable!() - } - - fn weighted_densities(&self, density: &Array1) -> Vec> { - self.weight_constants - .iter() - .map(|w| w.dot(density)) - .collect() - } - - fn functional_derivative(&self, partial_derivatives: &[Array1]) -> Array1 { - self.weight_constants - .iter() - .zip(partial_derivatives.iter()) - .map(|(w, pd)| pd.dot(w)) - .reduce(|a, b| a + b) - .unwrap() - } -} - /// Base structure to hold either information about the weight function through /// `WeightFunctionInfo` or the weight functions themselves via /// `FFTWeightFunctions`. @@ -125,11 +83,10 @@ pub struct ConvolverFFT { cartesian_transforms: Vec>, } -impl ConvolverFFT +impl ConvolverFFT where T: DctNum + DualNum, D::Larger: Dimension, - D::Smaller: Dimension, ::Larger: Dimension, { /// Create the appropriate FFT convolver for the given grid. @@ -139,9 +96,11 @@ where lanczos: Option, ) -> Arc> { match grid { + Grid::Bulk => PeriodicConvolver::new_0d(weight_functions), Grid::Polar(r) => CurvilinearConvolver::new(r, &[], weight_functions, lanczos), Grid::Spherical(r) => CurvilinearConvolver::new(r, &[], weight_functions, lanczos), Grid::Cartesian1(z) => Self::new(Some(z), &[], weight_functions, lanczos), + Grid::Periodical1(z) => PeriodicConvolver::new_1d(z, weight_functions, lanczos), Grid::Cylindrical { r, z } => { CurvilinearConvolver::new(r, &[z], weight_functions, lanczos) } @@ -553,11 +512,10 @@ struct CurvilinearConvolver { convolver_boundary: Arc>, } -impl CurvilinearConvolver +impl CurvilinearConvolver where T: DctNum + DualNum, D::Larger: Dimension, - D::Smaller: Dimension, ::Larger: Dimension, { #[expect(clippy::new_ret_no_self)] @@ -574,10 +532,9 @@ where } } -impl Convolver for CurvilinearConvolver +impl Convolver for CurvilinearConvolver where T: DctNum + DualNum, - D::Smaller: Dimension, D::Larger: Dimension, { fn convolve( @@ -587,9 +544,9 @@ where ) -> Array { // subtract boundary profile from full profile let profile_boundary = profile - .index_axis(Axis(0), profile.shape()[0] - 1) + .slice_axis(Axis(0), Slice::from(profile.shape()[0] - 1..)) .into_owned(); - for mut lane in profile.outer_iter_mut() { + for mut lane in profile.axis_chunks_iter_mut(Axis(0), 1) { lane.sub_assign(&profile_boundary); } @@ -597,14 +554,12 @@ where let mut result = self.convolver.convolve(profile, weight_function); // convolve boundary profile - let profile_boundary = profile_boundary.insert_axis(Axis(0)); let result_boundary = self .convolver_boundary .convolve(profile_boundary, weight_function); // Add boundary result back to full result - let result_boundary = result_boundary.index_axis(Axis(0), 0); - for mut lane in result.outer_iter_mut() { + for mut lane in result.axis_chunks_iter_mut(Axis(0), 1) { lane.add_assign(&result_boundary); } result diff --git a/crates/feos-dft/src/convolver/periodic_convolver.rs b/crates/feos-dft/src/convolver/periodic_convolver.rs index bd594e451..85c24fd18 100644 --- a/crates/feos-dft/src/convolver/periodic_convolver.rs +++ b/crates/feos-dft/src/convolver/periodic_convolver.rs @@ -30,6 +30,18 @@ where D::Larger: Dimension, ::Larger: Dimension, { + pub fn new_0d(weight_functions: &[WeightFunctionInfo]) -> Arc> { + Self::new(&[], |_| {}, weight_functions, None) + } + + pub fn new_1d( + axis: &Axis, + weight_functions: &[WeightFunctionInfo], + lanczos: Option, + ) -> Arc> { + Self::new(&[axis], |_| {}, weight_functions, lanczos) + } + pub fn new_2d( axes: &[&Axis], angle: Angle, @@ -200,7 +212,7 @@ impl PeriodicConvolver { } fn forward_transform(&self, f: ArrayView) -> Array, D2> { - let offset = D2::NDIM.unwrap() - D::NDIM.unwrap(); + let offset = f.ndim() - self.k_abs.ndim(); let mut result = f.mapv(Complex::from); for (i, transform) in self.forward_transforms.iter().enumerate() { for r in result.lanes_mut(Axis_nd(i + offset)).into_iter() { @@ -211,7 +223,7 @@ impl PeriodicConvolver { } fn inverse_transform(&self, mut f: Array, D2>) -> Array { - let offset = D2::NDIM.unwrap() - D::NDIM.unwrap(); + let offset = f.ndim() - self.k_abs.ndim(); for (i, transform) in self.inverse_transforms.iter().enumerate() { for r in f.lanes_mut(Axis_nd(i + offset)).into_iter() { self.transform(transform, r); diff --git a/crates/feos-dft/src/geometry.rs b/crates/feos-dft/src/geometry.rs index e69409e9b..b4a14f5d1 100644 --- a/crates/feos-dft/src/geometry.rs +++ b/crates/feos-dft/src/geometry.rs @@ -6,7 +6,9 @@ use std::f64::consts::{FRAC_PI_3, PI}; /// Grids with up to three dimensions. #[derive(Clone)] pub enum Grid { + Bulk, Cartesian1(Axis), + Periodical1(Axis), Cartesian2(Axis, Axis), Periodical2(Axis, Axis, Angle), Cartesian3(Axis, Axis, Axis), @@ -27,7 +29,8 @@ impl Grid { pub fn axes(&self) -> Vec<&Axis> { match self { - Self::Cartesian1(x) => vec![x], + Self::Bulk => vec![], + Self::Cartesian1(x) | Self::Periodical1(x) => vec![x], Self::Cartesian2(x, y) | Self::Periodical2(x, y, _) => vec![x, y], Self::Cartesian3(x, y, z) | Self::Periodical3(x, y, z, _) => vec![x, y, z], Self::Spherical(r) | Self::Polar(r) => vec![r], @@ -37,7 +40,8 @@ impl Grid { pub fn axes_mut(&mut self) -> Vec<&mut Axis> { match self { - Self::Cartesian1(x) => vec![x], + Self::Bulk => vec![], + Self::Cartesian1(x) | Self::Periodical1(x) => vec![x], Self::Cartesian2(x, y) | Self::Periodical2(x, y, _) => vec![x, y], Self::Cartesian3(x, y, z) | Self::Periodical3(x, y, z, _) => vec![x, y, z], Self::Spherical(r) | Self::Polar(r) => vec![r], @@ -72,7 +76,11 @@ impl Grid { pub fn mesh(&self) -> Vec>> { match self { - Grid::Cartesian1(ax) | Grid::Spherical(ax) | Grid::Polar(ax) => { + Self::Bulk => vec![], + Grid::Cartesian1(ax) + | Self::Periodical1(ax) + | Grid::Spherical(ax) + | Grid::Polar(ax) => { vec![Length::from_reduced(ax.grid.clone()).into_dyn()] } Grid::Cartesian2(u, v) => mesh_2d(u, v, 90.0 * DEGREES), diff --git a/crates/feos-dft/src/profile/mod.rs b/crates/feos-dft/src/profile/mod.rs index 6497faa23..a604cb76f 100644 --- a/crates/feos-dft/src/profile/mod.rs +++ b/crates/feos-dft/src/profile/mod.rs @@ -1,10 +1,10 @@ -use crate::convolver::{BulkConvolver, Convolver, ConvolverFFT}; +use crate::convolver::{Convolver, ConvolverFFT, PeriodicConvolver}; use crate::functional::HelmholtzEnergyFunctional; use crate::geometry::Grid; use crate::solver::{DFTSolver, DFTSolverLog}; use feos_core::{FeosError, FeosResult, ReferenceSystem, State}; use nalgebra::{DVector, Dyn, U1}; -use ndarray::{Array, Array1, ArrayBase, Axis as Axis_nd, Data, Dimension, RemoveAxis}; +use ndarray::{Array, Array1, ArrayBase, Axis as Axis_nd, Data, Dimension, Ix0, arr1}; use num_dual::DualNum; use quantity::{_Volume, Density, Energy, Entropy, Length, Moles, Quantity, Temperature, Volume}; use std::ops::{Add, MulAssign}; @@ -22,7 +22,7 @@ const MAX_POTENTIAL: f64 = 50.0; #[derive(Clone)] pub enum DFTSpecification { /// DFT with specified chemical potential. - ChemicalPotential, + ChemicalPotential(Array1), /// DFT with specified number of particles. /// /// The solution is still a grand canonical density profile, but the chemical @@ -30,21 +30,42 @@ pub enum DFTSpecification { /// with the specified number of particles. Moles(Array1), /// DFT with specified total number of moles. - TotalMoles(f64), + TotalMoles(f64, Array1), } impl DFTSpecification { - fn calculate_bulk_density( - &self, - bulk_density: &Array1, - z: &Array1, - ) -> FeosResult> { + fn calculate_fugacity(&self, z: &Array1) -> FeosResult> { Ok(match self { - Self::ChemicalPotential => bulk_density.clone(), + Self::ChemicalPotential(fugacity) => fugacity.clone(), Self::Moles(moles) => moles / z, - Self::TotalMoles(total_moles) => bulk_density * *total_moles / (bulk_density * z).sum(), + Self::TotalMoles(total_moles, fugacity) => { + fugacity * *total_moles / (fugacity * z).sum() + } }) } + + pub fn from_state(state: &State) -> Self { + let component_index = state.eos.component_index().into_owned(); + let m = arr1(&state.eos.m()); + let partial_density = state.partial_density().into_reduced(); + let temperature = state.temperature.into_reduced(); + let bulk_density = component_index + .iter() + .map(|&i| partial_density[i]) + .collect(); + let bulk_convolver = + PeriodicConvolver::<_, Ix0>::new_0d(&state.eos.weight_functions(temperature)); + let (_, dfdrho_bulk) = state + .eos + .functional_derivative(temperature, &bulk_density, bulk_convolver.as_ref()) + .unwrap(); + let exp_dfdrho = (dfdrho_bulk / m).mapv(f64::exp); + let bonds = state + .eos + .bond_integrals(temperature, &exp_dfdrho, bulk_convolver.as_ref()); + let fugacity = bulk_density * exp_dfdrho * bonds; + Self::ChemicalPotential(fugacity) + } } #[derive(Clone)] @@ -80,10 +101,9 @@ impl DFTProfile { } } -impl DFTProfile +impl DFTProfile where D::Larger: Dimension, - D::Smaller: Dimension, ::Larger: Dimension, { /// Create a new density profile. @@ -147,7 +167,7 @@ where convolver, temperature: bulk.temperature, density, - specification: DFTSpecification::ChemicalPotential, + specification: DFTSpecification::from_state(bulk), external_potential, bulk: bulk.clone(), solver_log: None, @@ -167,7 +187,12 @@ where pub fn fix_total_moles(&mut self) { let rho = self.density.to_reduced(); let moles = self.integrate_reduced_comp(&rho).sum(); - self.specification = DFTSpecification::TotalMoles(moles); + let DFTSpecification::ChemicalPotential(fugacity) = + DFTSpecification::from_state(&self.bulk) + else { + unreachable!() + }; + self.specification = DFTSpecification::TotalMoles(moles, fugacity); } /// Return the external potential in SI units. @@ -290,37 +315,17 @@ where .weighted_densities(&self.density.to_reduced())) } - #[expect(clippy::type_complexity)] - pub fn residual(&self, log: bool) -> FeosResult<(Array, Array1, f64)> { - // Read from profile + pub fn residual(&self, log: bool) -> FeosResult<(Array, f64)> { let density = self.density.to_reduced(); - let partial_density = self.bulk.partial_density().into_reduced(); - let mut bulk_density = self - .bulk - .eos - .component_index() - .iter() - .map(|&i| partial_density[i]) - .collect(); - - let (res, res_bulk, res_norm, _, _) = - self.euler_lagrange_equation(&density, &mut bulk_density, log)?; - Ok((res, res_bulk, res_norm)) + let (res, res_norm, _, _) = self.euler_lagrange_equation(&density, log)?; + Ok((res, res_norm)) } #[expect(clippy::type_complexity)] - pub(crate) fn euler_lagrange_equation( + fn fugacity( &self, density: &Array, - bulk_density: &mut Array1, - log: bool, - ) -> FeosResult<( - Array, - Array1, - f64, - Array, - Array, - )> { + ) -> FeosResult<(Array, Array, Array1)> { // calculate reduced temperature let temperature = self.temperature.to_reduced(); @@ -333,21 +338,10 @@ where // calculate total functional derivative dfdrho += &self.external_potential; - // calculate bulk functional derivative - let bulk_convolver = BulkConvolver::new(self.bulk.eos.weight_functions(temperature)); - let (_, dfdrho_bulk) = self.bulk.eos.functional_derivative( - temperature, - bulk_density, - bulk_convolver.as_ref(), - )?; dfdrho .outer_iter_mut() - .zip(dfdrho_bulk) .zip(self.bulk.eos.m().iter()) - .for_each(|((mut df, df_b), &m)| { - df -= df_b; - df /= m - }); + .for_each(|(mut df, &m)| df /= m); // calculate bond integrals let exp_dfdrho = dfdrho.mapv(|x| (-x).exp()); @@ -355,22 +349,36 @@ where .bulk .eos .bond_integrals(temperature, &exp_dfdrho, self.convolver.as_ref()); - let mut rho_projected = &exp_dfdrho * bonds; + let z = &exp_dfdrho * bonds; - // calculate bulk density based on the given specification - let res_bulk = &*bulk_density - - self.specification.calculate_bulk_density( - bulk_density, - &self.integrate_reduced_comp(&rho_projected), - )?; - *bulk_density -= &res_bulk; + // calculate fugacity based on the given specification + let fugacity = self + .specification + .calculate_fugacity(&self.integrate_reduced_comp(&z))?; - // multiply bulk density + Ok((exp_dfdrho, z, fugacity)) + } + + #[expect(clippy::type_complexity)] + pub(crate) fn euler_lagrange_equation( + &self, + density: &Array, + log: bool, + ) -> FeosResult<( + Array, + f64, + Array, + Array, + )> { + // calculate functional derivatives and fugacity + let (exp_dfdrho, mut rho_projected, fugacity) = self.fugacity(density)?; + + // multiply fugacity rho_projected .outer_iter_mut() - .zip(bulk_density.iter()) - .for_each(|(mut x, &rho_b)| { - x *= rho_b; + .zip(fugacity.iter()) + .for_each(|(mut x, &f)| { + x *= f; }); // calculate residual @@ -391,7 +399,7 @@ where (density - &rho_projected).mapv(|x| x * x).sum().sqrt() / (res.len() as f64).sqrt(); if res_norm.is_finite() { - Ok((res, res_bulk, res_norm, exp_dfdrho, rho_projected)) + Ok((res, res_norm, exp_dfdrho, rho_projected)) } else { Err(FeosError::IterationFailed("Euler-Lagrange equation".into())) } @@ -402,25 +410,35 @@ where let solver = solver.cloned().unwrap_or_default(); // Read from profile - let component_index = self.bulk.eos.component_index().into_owned(); let mut density = self.density.to_reduced(); - let partial_density = self.bulk.partial_density().into_reduced(); - let mut bulk_density = component_index - .iter() - .map(|&i| partial_density[i]) - .collect(); // Call solver(s) - self.call_solver(&mut density, &mut bulk_density, &solver, debug)?; + self.call_solver(&mut density, &solver, debug)?; + + // Update bulk state + if !matches!(self.specification, DFTSpecification::ChemicalPotential(_)) { + // 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)?; + bulk_profile.specification = DFTSpecification::ChemicalPotential(fugacity); + let solver = DFTSolver::new(None).newton(None, None, None, None); + bulk_profile.solve(Some(&solver), false)?; + + // create the state based on the results from the bulk profile + let component_index = self.bulk.eos.component_index(); + let mut partial_density = self.bulk.partial_density(); + bulk_profile + .density + .into_iter() + .enumerate() + .for_each(|(i, r)| partial_density.set(component_index[i], r)); + self.bulk = State::new_density(&self.bulk.eos, self.bulk.temperature, partial_density)?; + } // Update profile self.density = Density::from_reduced(density); - let mut partial_density = self.bulk.partial_density(); - bulk_density - .into_iter() - .enumerate() - .for_each(|(i, r)| partial_density.set(component_index[i], Density::from_reduced(r))); - self.bulk = State::new_density(&self.bulk.eos, self.bulk.temperature, partial_density)?; + Ok(()) } } diff --git a/crates/feos-dft/src/profile/properties.rs b/crates/feos-dft/src/profile/properties.rs index 2fe6a92a9..c506fac24 100644 --- a/crates/feos-dft/src/profile/properties.rs +++ b/crates/feos-dft/src/profile/properties.rs @@ -1,11 +1,11 @@ #![allow(type_alias_bounds)] use super::DFTProfile; -use crate::convolver::{BulkConvolver, Convolver}; +use crate::convolver::{Convolver, PeriodicConvolver}; use crate::functional_contribution::FunctionalContribution; use crate::{ConvolverFFT, DFTSolverLog, HelmholtzEnergyFunctional, WeightFunctionInfo}; use feos_core::{Contributions, FeosResult, ReferenceSystem, Total, Verbosity}; use nalgebra::{DMatrix, DVector}; -use ndarray::{Array, Array1, Axis, Dimension, RemoveAxis}; +use ndarray::{Array, Array1, Axis, Dimension, Ix0, RemoveAxis}; use num_dual::{Dual64, DualNum}; use quantity::{ Density, Energy, Entropy, EntropyDensity, MolarEnergy, Moles, Pressure, Quantity, Temperature, @@ -298,17 +298,8 @@ where { fn density_derivative(&self, lhs: &Array) -> FeosResult> { let rho = self.density.to_reduced(); - let partial_density = self.bulk.partial_density().into_reduced(); - let mut rho_bulk = self - .bulk - .eos - .component_index() - .iter() - .map(|&i| partial_density[i]) - .collect(); - let second_partial_derivatives = self.second_partial_derivatives(&rho)?; - let (_, _, _, exp_dfdrho, _) = self.euler_lagrange_equation(&rho, &mut rho_bulk, false)?; + let (_, _, exp_dfdrho, _) = self.euler_lagrange_equation(&rho, false)?; let rhs = |x: &_| { let delta_functional_derivative = @@ -411,7 +402,7 @@ where .map(|&i| partial_density[i]) .collect(); let rho_bulk_dual = rho_bulk.mapv(Dual64::from); - let bulk_convolver = BulkConvolver::new(weight_functions); + let bulk_convolver = PeriodicConvolver::<_, Ix0>::new_0d(&weight_functions); let (_, dfdrho_bulk) = self.bulk .eos diff --git a/crates/feos-dft/src/solver.rs b/crates/feos-dft/src/solver.rs index 69d2e2757..2baba60fb 100644 --- a/crates/feos-dft/src/solver.rs +++ b/crates/feos-dft/src/solver.rs @@ -216,7 +216,6 @@ where pub(crate) fn call_solver( &mut self, rho: &mut Array, - rho_bulk: &mut Array1, solver: &DFTSolver, debug: bool, ) -> FeosResult<()> { @@ -225,13 +224,11 @@ where let mut log = DFTSolverLog::new(solver.verbosity); for algorithm in &solver.algorithms { let (conv, iter) = match algorithm { - DFTAlgorithm::PicardIteration(picard) => { - self.solve_picard(*picard, rho, rho_bulk, &mut log) - } + DFTAlgorithm::PicardIteration(picard) => self.solve_picard(*picard, rho, &mut log), DFTAlgorithm::AndersonMixing(anderson) => { - self.solve_anderson(*anderson, rho, rho_bulk, &mut log) + self.solve_anderson(*anderson, rho, &mut log) } - DFTAlgorithm::Newton(newton) => self.solve_newton(*newton, rho, rho_bulk, &mut log), + DFTAlgorithm::Newton(newton) => self.solve_newton(*newton, rho, &mut log), }?; converged = conv; iterations += iter; @@ -255,7 +252,6 @@ where &self, picard: PicardIteration, rho: &mut Array, - rho_bulk: &mut Array1, log: &mut DFTSolverLog, ) -> FeosResult<(bool, usize)> { let solver = if picard.log { @@ -266,8 +262,7 @@ where for k in 0..picard.max_iter { // calculate residual - let (res, _, res_norm, _, _) = - self.euler_lagrange_equation(&*rho, rho_bulk, picard.log)?; + let (res, res_norm, _, _) = self.euler_lagrange_equation(&*rho, picard.log)?; log.add_residual(solver, k, res_norm); // check for convergence @@ -276,10 +271,9 @@ where } // apply line search or constant damping - let damping_coefficient = picard.damping_coefficient.map_or_else( - || self.line_search(rho, &res, rho_bulk, res_norm, picard.log), - Ok, - )?; + let damping_coefficient = picard + .damping_coefficient + .map_or_else(|| self.line_search(rho, &res, res_norm, picard.log), Ok)?; // update solution if picard.log { @@ -295,7 +289,6 @@ where &self, rho: &Array, delta_rho: &Array, - rho_bulk: &mut Array1, res0: f64, logarithm: bool, ) -> FeosResult { @@ -311,9 +304,7 @@ where } else { rho + alpha * delta_rho }; - let Ok((_, _, res2, _, _)) = - self.euler_lagrange_equation(&rho_new, rho_bulk, logarithm) - else { + let Ok((_, res2, _, _)) = self.euler_lagrange_equation(&rho_new, logarithm) else { continue; }; if res2 > res0 { @@ -326,9 +317,7 @@ where } else { rho + 0.5 * alpha * delta_rho }; - let Ok((_, _, res1, _, _)) = - self.euler_lagrange_equation(&rho_new, rho_bulk, logarithm) - else { + let Ok((_, res1, _, _)) = self.euler_lagrange_equation(&rho_new, logarithm) else { continue; }; @@ -359,7 +348,6 @@ where &self, anderson: AndersonMixing, rho: &mut Array, - rho_bulk: &mut Array1, log: &mut DFTSolverLog, ) -> FeosResult<(bool, usize)> { let solver = if anderson.log { @@ -382,8 +370,7 @@ where let m = resm.len() + 1; // calculate residual - let (res, _, res_norm, _, _) = - self.euler_lagrange_equation(&*rho, rho_bulk, anderson.log)?; + let (res, res_norm, _, _) = self.euler_lagrange_equation(&*rho, anderson.log)?; log.add_residual(solver, k, res_norm); // check for convergence @@ -434,14 +421,13 @@ where &self, newton: Newton, rho: &mut Array, - rho_bulk: &mut Array1, log: &mut DFTSolverLog, ) -> FeosResult<(bool, usize)> { 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) = - self.euler_lagrange_equation(rho, rho_bulk, newton.log)?; + let (res, res_norm, exp_dfdrho, rho_p) = + self.euler_lagrange_equation(rho, newton.log)?; log.add_residual(solver, k, res_norm); // check convergence diff --git a/py-feos/src/dft/adsorption/pore.rs b/py-feos/src/dft/adsorption/pore.rs index ab66684e2..d07c75b39 100644 --- a/py-feos/src/dft/adsorption/pore.rs +++ b/py-feos/src/dft/adsorption/pore.rs @@ -28,6 +28,13 @@ impl PyGrid { Self(Grid::Cartesian1(x)) } + /// Generate a 1D Cartesian grid with periodic boundary conditions on both sides. + #[staticmethod] + pub fn periodical_1d(n_points: usize, length: Length) -> Self { + let x = AxisDFT::new_cartesian(n_points, length, None); + Self(Grid::Periodical1(x)) + } + /// Generate a polar grid with radial axis. #[staticmethod] pub fn polar(n_points: usize, length: Length) -> Self { diff --git a/py-feos/src/dft/profile.rs b/py-feos/src/dft/profile.rs index 40104de89..e1acef783 100644 --- a/py-feos/src/dft/profile.rs +++ b/py-feos/src/dft/profile.rs @@ -18,9 +18,9 @@ macro_rules! impl_profile { &self, log: bool, py: Python<'py>, - ) -> PyResult<(Bound<'py, PyArrayDyn>, Bound<'py, PyArray1>, f64)> { - let (res_rho, res_mu, res_norm) = self.0.profile.residual(log).map_err(PyFeosError::from)?; - Ok((res_rho.view().into_dyn().to_pyarray(py), res_mu.view().to_pyarray(py), res_norm)) + ) -> PyResult<(Bound<'py, PyArrayDyn>, f64)> { + let (res_rho, res_norm) = self.0.profile.residual(log).map_err(PyFeosError::from)?; + Ok((res_rho.view().into_dyn().to_pyarray(py), res_norm)) } /// Solve the profile in-place. A non-default solver can be provided