Skip to content

Fix calculation of bulk state when moles are specified #382

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
prehner merged 1 commit into from
Aug 14, 2026
Merged

Fix calculation of bulk state when moles are specified #382

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions crates/feos-dft/src/adsorption/pore.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,6 +162,7 @@ where
}

pub fn update_bulk(mut self, bulk: &State<F>) -> Self {
self.profile.specification = DFTSpecification::from_state(bulk);
self.profile.bulk = bulk.clone();
self.grand_potential = None;
self.interfacial_tension = None;
Expand Down
61 changes: 8 additions & 53 deletions crates/feos-dft/src/convolver/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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;
Expand Down Expand Up @@ -36,46 +34,6 @@ pub trait Convolver<T, D: Dimension>: Send + Sync {
) -> Array<T, D::Larger>;
}

pub(crate) struct BulkConvolver<T> {
weight_constants: Vec<Array2<T>>,
}

impl<T: DualNum<f64> + Copy + Send + Sync> BulkConvolver<T> {
#[expect(clippy::new_ret_no_self)]
pub(crate) fn new(weight_functions: Vec<WeightFunctionInfo<T>>) -> Arc<dyn Convolver<T, Ix0>> {
let weight_constants = weight_functions
.into_iter()
.map(|w| w.weight_constants(Zero::zero(), 0))
.collect();
Arc::new(Self { weight_constants })
}
}

impl<T: DualNum<f64> + Copy + Send + Sync> Convolver<T, Ix0> for BulkConvolver<T>
where
Array2<T>: Dot<Array1<T>, Output = Array1<T>>,
{
fn convolve(&self, _: Array0<T>, _: &WeightFunction<T>) -> Array0<T> {
unreachable!()
}

fn weighted_densities(&self, density: &Array1<T>) -> Vec<Array1<T>> {
self.weight_constants
.iter()
.map(|w| w.dot(density))
.collect()
}

fn functional_derivative(&self, partial_derivatives: &[Array1<T>]) -> Array1<T> {
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`.
Expand Down Expand Up @@ -125,11 +83,10 @@ pub struct ConvolverFFT<T, D: Dimension> {
cartesian_transforms: Vec<CartesianTransform<T>>,
}

impl<T, D: Dimension + RemoveAxis + 'static> ConvolverFFT<T, D>
impl<T, D: Dimension + 'static> ConvolverFFT<T, D>
where
T: DctNum + DualNum<f64>,
D::Larger: Dimension<Smaller = D>,
D::Smaller: Dimension<Larger = D>,
<D::Larger as Dimension>::Larger: Dimension<Smaller = D::Larger>,
{
/// Create the appropriate FFT convolver for the given grid.
Expand All @@ -139,9 +96,11 @@ where
lanczos: Option<i32>,
) -> Arc<dyn Convolver<T, D>> {
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)
}
Expand Down Expand Up @@ -553,11 +512,10 @@ struct CurvilinearConvolver<T, D> {
convolver_boundary: Arc<dyn Convolver<T, D>>,
}

impl<T, D: Dimension + RemoveAxis + 'static> CurvilinearConvolver<T, D>
impl<T, D: Dimension + 'static> CurvilinearConvolver<T, D>
where
T: DctNum + DualNum<f64>,
D::Larger: Dimension<Smaller = D>,
D::Smaller: Dimension<Larger = D>,
<D::Larger as Dimension>::Larger: Dimension<Smaller = D::Larger>,
{
#[expect(clippy::new_ret_no_self)]
Expand All @@ -574,10 +532,9 @@ where
}
}

impl<T, D: Dimension + RemoveAxis> Convolver<T, D> for CurvilinearConvolver<T, D>
impl<T, D: Dimension> Convolver<T, D> for CurvilinearConvolver<T, D>
where
T: DctNum + DualNum<f64>,
D::Smaller: Dimension<Larger = D>,
D::Larger: Dimension<Smaller = D>,
{
fn convolve(
Expand All @@ -587,24 +544,22 @@ where
) -> Array<T, D> {
// 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);
}

// convolve full profile
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
Expand Down
16 changes: 14 additions & 2 deletions crates/feos-dft/src/convolver/periodic_convolver.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,18 @@ where
D::Larger: Dimension<Smaller = D>,
<D::Larger as Dimension>::Larger: Dimension<Smaller = D::Larger>,
{
pub fn new_0d(weight_functions: &[WeightFunctionInfo<T>]) -> Arc<dyn Convolver<T, D>> {
Self::new(&[], |_| {}, weight_functions, None)
}

pub fn new_1d(
axis: &Axis,
weight_functions: &[WeightFunctionInfo<T>],
lanczos: Option<i32>,
) -> Arc<dyn Convolver<T, D>> {
Self::new(&[axis], |_| {}, weight_functions, lanczos)
}

pub fn new_2d(
axes: &[&Axis],
angle: Angle,
Expand Down Expand Up @@ -200,7 +212,7 @@ impl<T: FftNum, D: Dimension> PeriodicConvolver<T, D> {
}

fn forward_transform<D2: Dimension>(&self, f: ArrayView<T, D2>) -> Array<Complex<T>, 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() {
Expand All @@ -211,7 +223,7 @@ impl<T: FftNum, D: Dimension> PeriodicConvolver<T, D> {
}

fn inverse_transform<D2: Dimension>(&self, mut f: Array<Complex<T>, D2>) -> Array<T, D2> {
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);
Expand Down
14 changes: 11 additions & 3 deletions crates/feos-dft/src/geometry.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand All @@ -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],
Expand All @@ -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],
Expand Down Expand Up @@ -72,7 +76,11 @@ impl Grid {

pub fn mesh(&self) -> Vec<Length<ArrayD<f64>>> {
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),
Expand Down
Loading