feos-org / feos

FeOs - A Framework for Equations of State and Classical Density Functional Theory
Other
120 stars 23 forks source link

Add non-mixed second order derivatives to `State` property evaluation #94

Closed g-bauer closed 1 year ago

g-bauer commented 1 year ago

This PR introduces the new variant PartialDerivative::Second and renames the existing variant to PartialDerivative::SecondMixed. PartialDerivative::Second can be used to calculate second order derivatives using Dual2 instead of HyperDual, which is sufficient for non-mixed derivatives and computationally less expensive.

The new PartialDerivative looks like so

#[derive(Clone, Copy, Eq, Hash, PartialEq, Debug)]
pub(crate) enum PartialDerivative {
    Zeroth,
    First(Derivative),
    Second(Derivative),
    SecondMixed(Derivative, Derivative),
    Third(Derivative),
}

We can use this variant to calculate non-mixed second order derivatives like so:

fn dp_dv_(&self, evaluate: Evaluate) -> QuantityScalar<U> {
    -self.get_or_compute_derivative(PartialDerivative::Second(DV), evaluate) // this PR
//  -self.get_or_compute_derivative(PartialDerivative::Second(DV, DV), evaluate) // current main
}

fn dp_dt_(&self, evaluate: Evaluate) -> QuantityScalar<U> {
    -self.get_or_compute_derivative(PartialDerivative::SecondMixed(DV, DT), evaluate) // this PR
//  -self.get_or_compute_derivative(PartialDerivative::Second(DV, DT), evaluate) // current main
}

The Cache, in which partial derivatives are stored, puts results into PartialDerivative::SecondMixed so that there are no cache misses when the derivatives are calculated with the (less efficient) variant that uses HyperDuals.

// state/cache.rs, method of Cache
pub fn get_or_insert_with_d2_64<F: FnOnce() -> Dual2_64>(
    &mut self,
    derivative: Derivative,
    f: F,
) -> f64 {
    if let Some(&value) = self.map.get(&PartialDerivative::SecondMixed(derivative, derivative)) {
        self.hit += 1;
        value
    } else {
        self.miss += 1;
        let value = f();
        self.map.insert(PartialDerivative::Zeroth, value.re);
        self.map
            .insert(PartialDerivative::First(derivative), value.v1[0]);
        self.map
            .insert(PartialDerivative::SecondMixed(derivative, derivative), value.v2[0]);
        value.v2[0]
    }
}

This change benefits the density_iteration function in which we currently repeatedly call

// method of State
pub(crate) fn p_dpdrho(&self) -> (QuantityScalar<U>, QuantityScalar<U>) {
    let dp_dv = self.dp_dv(Contributions::Total); // <- currently uses HyperDual
    (
        self.pressure(Contributions::Total),
        (-self.volume * dp_dv / self.density),
    )
}

In the added benchmarks for different ways to create States, this change shows an ~8-9% improvement (i.e. execution speed) for the State::new_npt and an ~2-4% improvement for PhaseEquilibrium::tp_flash although more systems should be added to get a proper idea of the impact.

g-bauer commented 1 year ago

Here are some benchmarks for PC-SAFT comparing current main to this PR. Reported are %-changes is execution time. For new_npt, liq. and vap. means that liquid or vapor densities were chosen as initial densities, respectively, so no stability analysis is performed. The temperatures and pressures are chosen so that those phases are stable. For the bubble and dew point benchmarks, the temperature is used as input.

methane methane + ethane methane + ethane + propane
State::new_npt (liq) -9.7020% -13.284% -13.934%
State::new_npt (vap) -9.4964% -14.938% -14.187%
PhaseEquilibrium::tp_flash -1.7909% -4.2929%
PhaseEquilibrium::bubble_point -3.7754% -4.8505%
PhaseEquilibrium::dew_point -5.5984% -9.8986%

Note that these benchmarks fluctuate quite heavily (around +-1.5%) when ran multiple times. Not sure how to get more stable results.