Skip to content

Commit 41ba5fa

Browse files
prehnerprehner
committed
Extend Newton solver to DFT specifications other than chemical potentials
1 parent 5af3237 commit 41ba5fa

3 files changed

Lines changed: 55 additions & 29 deletions

File tree

‎crates/feos-dft/src/profile/mod.rs‎

Lines changed: 37 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -34,14 +34,24 @@ pub enum DFTSpecification {
3434
}
3535

3636
impl DFTSpecification {
37-
fn calculate_fugacity(&self, z: &Array1<f64>) -> FeosResult<Array1<f64>> {
38-
Ok(match self {
37+
fn calculate_fugacity(&self, z: &Array1<f64>) -> Array1<f64> {
38+
match self {
3939
Self::ChemicalPotential(fugacity) => fugacity.clone(),
4040
Self::Moles(moles) => moles / z,
4141
Self::TotalMoles(total_moles, fugacity) => {
4242
fugacity * *total_moles / (fugacity * z).sum()
4343
}
44-
})
44+
}
45+
}
46+
47+
pub(crate) fn delta_fugacity(&self, z: &Array1<f64>, delta_z: &Array1<f64>) -> Array1<f64> {
48+
match self {
49+
Self::ChemicalPotential(fugacity) => Array1::zeros(fugacity.len()),
50+
Self::Moles(_) => -delta_z / z,
51+
Self::TotalMoles(_, fugacity) => {
52+
-(fugacity * delta_z).sum() / (fugacity * z).sum() * Array1::ones(fugacity.len())
53+
}
54+
}
4555
}
4656

4757
pub fn from_state<F: HelmholtzEnergyFunctional>(state: &State<F>) -> Self {
@@ -216,7 +226,7 @@ where
216226
profile.sum() * functional_determinant
217227
}
218228

219-
fn integrate_reduced_comp<S: Data<Elem = N>, N: DualNum<Primitive = f64> + Copy>(
229+
pub(crate) fn integrate_reduced_comp<S: Data<Elem = N>, N: DualNum<Primitive = f64> + Copy>(
220230
&self,
221231
profile: &ArrayBase<S, D::Larger>,
222232
) -> Array1<N> {
@@ -320,15 +330,20 @@ where
320330

321331
pub fn residual(&self, log: bool) -> FeosResult<(Array<f64, D::Larger>, f64)> {
322332
let density = self.density.to_reduced();
323-
let (res, res_norm, _, _) = self.euler_lagrange_equation(&density, log)?;
333+
let (res, res_norm, _, _, _) = self.euler_lagrange_equation(&density, log)?;
324334
Ok((res, res_norm))
325335
}
326336

327337
#[expect(clippy::type_complexity)]
328338
fn fugacity(
329339
&self,
330340
density: &Array<f64, D::Larger>,
331-
) -> FeosResult<(Array<f64, D::Larger>, Array<f64, D::Larger>, Array1<f64>)> {
341+
) -> FeosResult<(
342+
Array<f64, D::Larger>,
343+
Array1<f64>,
344+
Array<f64, D::Larger>,
345+
Array1<f64>,
346+
)> {
332347
// calculate reduced temperature
333348
let temperature = self.temperature.to_reduced();
334349

@@ -352,14 +367,21 @@ where
352367
.bulk
353368
.eos
354369
.bond_integrals(temperature, &exp_dfdrho, self.convolver.as_ref());
355-
let z = &exp_dfdrho * bonds;
370+
let mut rho_projected = &exp_dfdrho * bonds;
371+
let z = self.integrate_reduced_comp(&rho_projected);
356372

357373
// calculate fugacity based on the given specification
358-
let fugacity = self
359-
.specification
360-
.calculate_fugacity(&self.integrate_reduced_comp(&z))?;
374+
let fugacity = self.specification.calculate_fugacity(&z);
375+
376+
// multiply fugacity
377+
rho_projected
378+
.outer_iter_mut()
379+
.zip(fugacity.iter())
380+
.for_each(|(mut x, &f)| {
381+
x *= f;
382+
});
361383

362-
Ok((exp_dfdrho, z, fugacity))
384+
Ok((exp_dfdrho, z, rho_projected, fugacity))
363385
}
364386

365387
#[expect(clippy::type_complexity)]
@@ -371,18 +393,11 @@ where
371393
Array<f64, D::Larger>,
372394
f64,
373395
Array<f64, D::Larger>,
396+
Array1<f64>,
374397
Array<f64, D::Larger>,
375398
)> {
376399
// calculate functional derivatives and fugacity
377-
let (exp_dfdrho, mut rho_projected, fugacity) = self.fugacity(density)?;
378-
379-
// multiply fugacity
380-
rho_projected
381-
.outer_iter_mut()
382-
.zip(fugacity.iter())
383-
.for_each(|(mut x, &f)| {
384-
x *= f;
385-
});
400+
let (exp_dfdrho, z, rho_projected, _) = self.fugacity(density)?;
386401

387402
// calculate residual
388403
let mut res = if log {
@@ -402,7 +417,7 @@ where
402417
(density - &rho_projected).mapv(|x| x * x).sum().sqrt() / (res.len() as f64).sqrt();
403418

404419
if res_norm.is_finite() {
405-
Ok((res, res_norm, exp_dfdrho, rho_projected))
420+
Ok((res, res_norm, exp_dfdrho, z, rho_projected))
406421
} else {
407422
Err(FeosError::IterationFailed("Euler-Lagrange equation".into()))
408423
}
@@ -423,7 +438,7 @@ where
423438
// solve a bulk profile with the Newton solver
424439
let mut bulk_profile =
425440
DFTProfile::<Ix0, _>::new(Grid::Bulk, &self.bulk, None, None, None);
426-
let (_, _, fugacity) = self.fugacity(&density)?;
441+
let (_, _, _, fugacity) = self.fugacity(&density)?;
427442
bulk_profile.specification = DFTSpecification::ChemicalPotential(fugacity);
428443
let solver = DFTSolver::new(None).newton(None, None, None, None);
429444
bulk_profile.solve(Some(&solver), false)?;

‎crates/feos-dft/src/profile/properties.rs‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -299,7 +299,7 @@ where
299299
fn density_derivative(&self, lhs: &Array<f64, D::Larger>) -> FeosResult<Array<f64, D::Larger>> {
300300
let rho = self.density.to_reduced();
301301
let second_partial_derivatives = self.second_partial_derivatives(&rho)?;
302-
let (_, _, exp_dfdrho, _) = self.euler_lagrange_equation(&rho, false)?;
302+
let (_, _, exp_dfdrho, _, _) = self.euler_lagrange_equation(&rho, false)?;
303303

304304
let rhs = |x: &_| {
305305
let delta_functional_derivative =

‎crates/feos-dft/src/solver.rs‎

Lines changed: 17 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -262,7 +262,7 @@ where
262262

263263
for k in 0..picard.max_iter {
264264
// calculate residual
265-
let (res, res_norm, _, _) = self.euler_lagrange_equation(&*rho, picard.log)?;
265+
let (res, res_norm, _, _, _) = self.euler_lagrange_equation(&*rho, picard.log)?;
266266
log.add_residual(solver, k, res_norm);
267267

268268
// check for convergence
@@ -304,7 +304,7 @@ where
304304
} else {
305305
rho + alpha * delta_rho
306306
};
307-
let Ok((_, res2, _, _)) = self.euler_lagrange_equation(&rho_new, logarithm) else {
307+
let Ok((_, res2, _, _, _)) = self.euler_lagrange_equation(&rho_new, logarithm) else {
308308
continue;
309309
};
310310
if res2 > res0 {
@@ -317,7 +317,7 @@ where
317317
} else {
318318
rho + 0.5 * alpha * delta_rho
319319
};
320-
let Ok((_, res1, _, _)) = self.euler_lagrange_equation(&rho_new, logarithm) else {
320+
let Ok((_, res1, _, _, _)) = self.euler_lagrange_equation(&rho_new, logarithm) else {
321321
continue;
322322
};
323323

@@ -370,7 +370,7 @@ where
370370
let m = resm.len() + 1;
371371

372372
// calculate residual
373-
let (res, res_norm, _, _) = self.euler_lagrange_equation(&*rho, anderson.log)?;
373+
let (res, res_norm, _, _, _) = self.euler_lagrange_equation(&*rho, anderson.log)?;
374374
log.add_residual(solver, k, res_norm);
375375

376376
// check for convergence
@@ -426,7 +426,7 @@ where
426426
let solver = if newton.log { "Newton (log)" } else { "Newton" };
427427
for k in 0..newton.max_iter {
428428
// calculate initial residual
429-
let (res, res_norm, exp_dfdrho, rho_p) =
429+
let (res, res_norm, exp_dfdrho, z, rho_p) =
430430
self.euler_lagrange_equation(rho, newton.log)?;
431431
log.add_residual(solver, k, res_norm);
432432

@@ -447,8 +447,19 @@ where
447447
.zip(self.bulk.eos.m().iter())
448448
.for_each(|(mut q, &m)| q /= m);
449449
let delta_i = self.delta_bond_integrals(&exp_dfdrho, &delta_functional_derivative);
450+
let mut delta_exp_dfdrho = delta_functional_derivative - delta_i;
451+
let delta_z = -self.integrate_reduced_comp(&(&delta_exp_dfdrho * &exp_dfdrho));
452+
453+
let delta_fugacity = self.specification.delta_fugacity(&z, &delta_z);
454+
delta_exp_dfdrho
455+
.outer_iter_mut()
456+
.zip(delta_fugacity.iter())
457+
.for_each(|(mut z, &f)| {
458+
z -= f;
459+
});
460+
450461
let rho = if newton.log { &*rho } else { &rho_p };
451-
delta_rho + (delta_functional_derivative - delta_i) * rho
462+
delta_rho + delta_exp_dfdrho * rho
452463
};
453464

454465
// update solution

0 commit comments

Comments
 (0)