Skip to content

Extend Newton solver to DFT specifications other than chemical potentials #384

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 31, 2026
Merged

Extend Newton solver to DFT specifications other than chemical potentials #384

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
59 changes: 37 additions & 22 deletions crates/feos-dft/src/profile/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,14 +34,24 @@ pub enum DFTSpecification {
}

impl DFTSpecification {
fn calculate_fugacity(&self, z: &Array1<f64>) -> FeosResult<Array1<f64>> {
Ok(match self {
fn calculate_fugacity(&self, z: &Array1<f64>) -> Array1<f64> {
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<f64>, delta_z: &Array1<f64>) -> Array1<f64> {
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<F: HelmholtzEnergyFunctional>(state: &State<F>) -> Self {
Expand Down Expand Up @@ -216,7 +226,7 @@ where
profile.sum() * functional_determinant
}

fn integrate_reduced_comp<S: Data<Elem = N>, N: DualNum<Primitive = f64> + Copy>(
pub(crate) fn integrate_reduced_comp<S: Data<Elem = N>, N: DualNum<Primitive = f64> + Copy>(
&self,
profile: &ArrayBase<S, D::Larger>,
) -> Array1<N> {
Expand Down Expand Up @@ -320,15 +330,20 @@ where

pub fn residual(&self, log: bool) -> FeosResult<(Array<f64, D::Larger>, 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))
}

#[expect(clippy::type_complexity)]
fn fugacity(
&self,
density: &Array<f64, D::Larger>,
) -> FeosResult<(Array<f64, D::Larger>, Array<f64, D::Larger>, Array1<f64>)> {
) -> FeosResult<(
Array<f64, D::Larger>,
Array1<f64>,
Array<f64, D::Larger>,
Array1<f64>,
)> {
// calculate reduced temperature
let temperature = self.temperature.to_reduced();

Expand All @@ -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)]
Expand All @@ -371,18 +393,11 @@ where
Array<f64, D::Larger>,
f64,
Array<f64, D::Larger>,
Array1<f64>,
Array<f64, D::Larger>,
)> {
// 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 {
Expand All @@ -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()))
}
Expand All @@ -423,7 +438,7 @@ where
// solve a bulk profile with the Newton solver
let mut bulk_profile =
DFTProfile::<Ix0, _>::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)?;
Expand Down
2 changes: 1 addition & 1 deletion crates/feos-dft/src/profile/properties.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -299,7 +299,7 @@ where
fn density_derivative(&self, lhs: &Array<f64, D::Larger>) -> FeosResult<Array<f64, D::Larger>> {
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 =
Expand Down
23 changes: 17 additions & 6 deletions crates/feos-dft/src/solver.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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 {
Expand All @@ -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;
};

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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);

Expand All @@ -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
Expand Down