Skip to content

Fix performance regressions in the Newton solver and isotherm calcula… #404

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
Oct 26, 2026
Merged

Fix performance regressions in the Newton solver and isotherm calcula… #404

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
54 changes: 24 additions & 30 deletions crates/feos-dft/src/adsorption/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -214,16 +214,16 @@ where
_ => unreachable!(),
};
}
let profile = PoreProfile::<D, _>::new(
grid.clone(),
&bulk,
external_potential,
None,
PoreSpecification::ChemicalPotential,
)
.solve(solver)?
.profile;
let mut old_density = Some(&profile.density);
let mut old_profile = Some(
PoreProfile::<D, _>::new(
grid.clone(),
&bulk,
external_potential,
None,
PoreSpecification::ChemicalPotential,
)
.solve(solver)?,
);

for i in 0..pressure.len() {
let mut bulk = State::new_npt(
Expand All @@ -240,27 +240,21 @@ where
.clone();
}

let p = PoreProfile::new(
grid.clone(),
&bulk,
external_potential,
old_density,
PoreSpecification::ChemicalPotential,
);
let p2 = PoreProfile::new(
grid.clone(),
&bulk,
external_potential,
None,
PoreSpecification::ChemicalPotential,
);
profiles.push(p.solve(solver).or_else(|_| p2.solve(solver)));

old_density = if let Some(Ok(l)) = profiles.last() {
Some(&l.profile.density)
} else {
None
// continue from the previous profile (density, fractions of non-bonded sites, and
// convolvers), the profile initialized from the bulk is only created if needed
let new_profile = || {
PoreProfile::new(
grid.clone(),
&bulk,
external_potential,
None,
PoreSpecification::ChemicalPotential,
)
};
let p = old_profile.map_or_else(new_profile, |p| p.update_bulk(&bulk));
profiles.push(p.solve(solver).or_else(|_| new_profile().solve(solver)));

old_profile = profiles.last().and_then(|p| p.as_ref().ok()).cloned();
}

Ok(Adsorption::new(functional, profiles))
Expand Down
26 changes: 19 additions & 7 deletions crates/feos-dft/src/solver/euler_lagrange.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -271,15 +271,27 @@ where
}

/// Norm of the difference between the density and the projected density, together with the
/// mass action residual.
/// mass action residual weighted with the density of the segment that carries the site (so
/// that both parts are densities and the tolerance does not depend on the density level).
fn convergence_norm(&self, x: &Array<f64, D::Larger>, residual: &Array<f64, D::Larger>) -> f64 {
let mass_action_residual = self.association_rows(residual);
let density_residual = Zip::from(&self.density_rows(x))
let rho = self.density_rows(x).mapv(f64::exp);
let density_residual = Zip::from(&rho)
.and(&self.density_rows(residual))
.map_collect(|&x, &r| x.exp() * (1.0 - r.exp()));
((density_residual.mapv(|r| r * r).sum() + mass_action_residual.mapv(|r| r * r).sum())
/ residual.len() as f64)
.sqrt()
.fold(0.0, |acc, &rho, &r| acc + (rho * (1.0 - r.exp())).powi(2));
let mass_action_residual = self.functional.association().map_or(0.0, |a| {
let p = a.sites();
p.a.iter()
.chain(p.b)
.chain(p.c)
.zip(self.association_rows(residual).outer_iter())
.map(|(s, r)| {
Zip::from(&r)
.and(&rho.index_axis(Axis(0), s.assoc_comp))
.fold(0.0, |acc, &r, &rho| acc + (rho * r).powi(2))
})
.sum()
});
((density_residual + mass_action_residual) / residual.len() as f64).sqrt()
}

/// Changes of the logarithms of the density and of X limited to [NEWTON_MAX_STEP].
Expand Down
15 changes: 10 additions & 5 deletions crates/feos-dft/src/solver/newton.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,10 +53,10 @@ pub(crate) fn newton_krylov<D: Dimension, S: NonlinearSystem<D>>(
) -> FeosResult<(bool, usize)> {
let solver = "Newton";

// the residual of the accepted step is reused, so the (more expensive) recording is only
// done if another Newton step is required
let mut res = system.residual(&mut Evaluate, x, solver_log)?;
for k in 0..max_iter {
// residual and linearization
let mut record = Record::new();
let res = system.residual(&mut record, x, solver_log)?;
let norm = system.convergence_norm(x, &res);
if !norm.is_finite() {
return Err(FeosError::IterationFailed(solver.into()));
Expand All @@ -67,6 +67,10 @@ pub(crate) fn newton_krylov<D: Dimension, S: NonlinearSystem<D>>(
}
let norm0 = system.norm(&res);

// linearization
let mut record = Record::new();
system.residual(&mut record, x, solver_log)?;

// inexact Newton: relative tolerance for the linear system, bounded by the absolute one
// (the replays are timed as a whole by GMRES)
let rhs = |delta: &Array<f64, D>| {
Expand All @@ -82,9 +86,10 @@ pub(crate) fn newton_krylov<D: Dimension, S: NonlinearSystem<D>>(
let mut accepted = false;
for _ in 0..NEWTON_MAX_LINE_SEARCH {
*x = &x0 + &(alpha * &delta);
if let Ok(res) = system.residual(&mut Evaluate, x, solver_log) {
let norm = system.norm(&res);
if let Ok(res_new) = system.residual(&mut Evaluate, x, solver_log) {
let norm = system.norm(&res_new);
if norm.is_finite() && norm < (1.0 - 1e-4 * alpha) * norm0 {
res = res_new;
accepted = true;
break;
}
Expand Down
2 changes: 1 addition & 1 deletion crates/feos/tests/gc_pcsaft/dft.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -286,7 +286,7 @@ fn test_association_functional() -> Result<(), Box<dyn Error>> {
),
];
// tight tolerance, so that the comparison is not limited by the convergence of the profiles
let solver = DFTSolver::new(None).newton(None, None, Some(1e-13));
let solver = DFTSolver::new(None).newton(None, None, Some(1e-14));
for (name, parameters, t, gamma_ref, pore_ref) in systems {
for version in [
AssociationVersion::YuWu,
Expand Down
Loading