Add benchmarks and Cargo profile using link-time optimization

.github/workflows/release.yml

@@ -17,7 +17,7 @@ jobs:
         with:
           manylinux: auto
           command: build
-          args: --release --out dist
+          args: --profile release-lto --out dist
       - name: Upload wheels
         uses: actions/upload-artifact@v2
         with:
@@ -42,11 +42,11 @@ jobs:
         uses: messense/maturin-action@v1
         with:
           target: x86_64
-          args: --release --out dist
+          args: --profile release-lto --out dist
       - name: Build wheels - universal2
         uses: messense/maturin-action@v1
         with:
-          args: --release --universal2 --out dist
+          args: --profile release-lto --universal2 --out dist
       - name: Upload wheels
         uses: actions/upload-artifact@v2
         with:
@@ -74,7 +74,7 @@ jobs:
         uses: messense/maturin-action@v1
         with:
           target: ${{ matrix.target }}
-          args: --release --out dist
+          args: --profile release-lto --out dist
       - name: Upload wheels
         uses: actions/upload-artifact@v2
         with:

CHANGELOG.md

@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 - Added SAFT-VRQ Mie equation of state and Helmholtz energy functional for first order Feynman-Hibbs corrected Mie fluids. [#79](https://github.com/feos-org/feos/pull/79)
 - Added `estimator` module to documentation. [#86](https://github.com/feos-org/feos/pull/86)
+- Added benchmarks for the evaluation of the Helmholtz energy and some properties of the `State` object for PC-SAFT. [#89](https://github.com/feos-org/feos/pull/89)
 
 ### Changed
 - Export `EosVariant` and `FunctionalVariant` directly in the crate root instead of their own modules. [#62](https://github.com/feos-org/feos/pull/62)

Cargo.toml

@@ -46,6 +46,19 @@ optional = true
 
 [dev-dependencies]
 approx = "0.4"
+criterion = "0.4"
+
+[[bench]]
+name = "state_properties"
+harness = false
+
+[[bench]]
+name = "dual_numbers"
+harness = false
+
+[profile.release-lto]
+inherits = "release"
+lto = true
 
 [features]
 default = []

README.md

@@ -96,37 +96,52 @@ If there is no compiled package for your system available from PyPI and you have
 pip install git+https://github.com/feos-org/feos
 ```
 
+This command builds the package without link-time optimization (LTO) that can be used to increase the performance further.
+See the *Building from source* section for information about building the wheel including LTO.
+
 ### Building from source
 
 To compile the code you need the Rust compiler and `maturin` (>=0.13,<0.14) installed.
-To install the package directly into the active environment, use
+To install the package directly into the active environment (virtualenv or conda), use
 
 ```
-maturin develop --release --features python
+maturin develop --release
 ```
 
-and specify the models that you want to include in the python package as additional features, e.g.
+which uses the `python` and `all_models` feature as specified in the `pyproject.toml` file.
+
+Alternatively, you can specify the models or features that you want to include in the python package explicitly, e.g.
 
 ```
 maturin develop --release --features "python pcsaft dft"
 ```
 
-for the PC-SAFT equation of state and Helmholtz energy functional. If you want to include all available models, use
+for the PC-SAFT equation of state and Helmholtz energy functional.
+
+To build wheels including link-time optimization (LTO), use
 
 ```
-maturin develop --release --features "python all_models"
+maturin build --profile="release-lto"
 ```
 
-To build wheels, use
+which will use the `python` and `all_models` features specified in the `pyproject.toml` file.
+Use the following command to build a wheel with specific features:
 
 ```
-maturin build --release --out dist --features "python ..."
+maturin build --profile="release-lto" --features "python ..."
 ```
 
+LTO increases compile times measurably but the resulting wheel is more performant and has a smaller size.
+For development however, we recommend using the `--release` flag.
+
 ## Documentation
 
 For a documentation of the Python API, Python examples, and a guide to the underlying Rust framework check out the [documentation](https://feos-org.github.io/feos/).
 
+## Benchmarks
+
+Check out the [benches](https://github.com/feos-org/feos/tree/main/benches) directory for information about provided Rust benchmarks and how to run them.
+
 ## Developers
 
 This software is currently maintained by members of the groups of

benches/README.md

@@ -0,0 +1,16 @@
+# Benchmarks
+
+This directory contains different benchmarks.
+For best performance, use the `release-lto` profile.
+Depending on the benchmark, you might have to consider different Cargo `features` as denoted in the table below.
+
+For example, to run the benchmarks in `dual_numbers`, which uses PC-SAFT, use
+
+```
+cargo bench --profile=release-lto --features=pcsaft --bench=dual_numbers
+``` 
+
+|Name|Description|Cargo features|
+|--|--|--|
+|`dual_numbers`|Helmholtz energy function evaluated using `StateHD` with different dual number types.|`pcsaft`|
+|`state_properties`|Properties of `State`. Including state creation using the natural variables of the Helmholtz energy (no density iteration).|`pcsaft`|

benches/dual_numbers.rs

@@ -0,0 +1,123 @@
+//! Benchmarks for the evaluation of the Helmholtz energy function
+//! for a given `StateHD` for different types of dual numbers.
+//! These should give an idea about the expected slow-down depending
+//! on the dual number type used without the overhead of the `State`
+//! creation.
+use criterion::{criterion_group, criterion_main, Criterion};
+use feos::pcsaft::{PcSaft, PcSaftParameters};
+use feos_core::{
+    parameter::{IdentifierOption, Parameter},
+    Derivative, EquationOfState, HelmholtzEnergy, HelmholtzEnergyDual, State, StateHD,
+};
+use ndarray::{arr1, Array};
+use num_dual::DualNum;
+use quantity::si::*;
+use std::sync::Arc;
+
+/// Helper function to create a state for given parameters.
+/// - temperature is 80% of critical temperature,
+/// - volume is critical volume,
+/// - molefracs (or moles) for equimolar mixture.
+fn state_pcsaft(parameters: PcSaftParameters) -> State<SIUnit, PcSaft> {
+    let n = parameters.pure_records.len();
+    let eos = Arc::new(PcSaft::new(Arc::new(parameters)));
+    let moles = Array::from_elem(n, 1.0 / n as f64) * 10.0 * MOL;
+    let cp = State::critical_point(&eos, Some(&moles), None, Default::default()).unwrap();
+    let temperature = 0.8 * cp.temperature;
+    State::new_nvt(&eos, temperature, cp.volume, &moles).unwrap()
+}
+
+/// Residual Helmholtz energy given an equation of state and a StateHD.
+fn a_res<D: DualNum<f64>, E: EquationOfState>(inp: (&Arc<E>, &StateHD<D>)) -> D
+where
+    (dyn HelmholtzEnergy + 'static): HelmholtzEnergyDual<D>,
+{
+    inp.0.evaluate_residual(inp.1)
+}
+
+/// Benchmark for evaluation of the Helmholtz energy for different dual number types.
+fn bench_dual_numbers<E: EquationOfState>(
+    c: &mut Criterion,
+    group_name: &str,
+    state: State<SIUnit, E>,
+) {
+    let mut group = c.benchmark_group(group_name);
+    group.bench_function("a_f64", |b| {
+        b.iter(|| a_res((&state.eos, &state.derive0())))
+    });
+    group.bench_function("a_dual", |b| {
+        b.iter(|| a_res((&state.eos, &state.derive1(Derivative::DV))))
+    });
+    group.bench_function("a_hyperdual", |b| {
+        b.iter(|| a_res((&state.eos, &state.derive2(Derivative::DV, Derivative::DV))))
+    });
+    group.bench_function("a_dual3", |b| {
+        b.iter(|| a_res((&state.eos, &state.derive3(Derivative::DV))))
+    });
+}
+
+/// Benchmark for the PC-SAFT equation of state
+fn pcsaft(c: &mut Criterion) {
+    // methane
+    let parameters = PcSaftParameters::from_json(
+        vec!["methane"],
+        "./parameters/pcsaft/gross2001.json",
+        None,
+        IdentifierOption::Name,
+    )
+    .unwrap();
+    bench_dual_numbers(c, "methane", state_pcsaft(parameters));
+
+    // water (4C, polar)
+    let parameters = PcSaftParameters::from_json(
+        vec!["water_4C_polar"],
+        "./parameters/pcsaft/rehner2020.json",
+        None,
+        IdentifierOption::Name,
+    )
+    .unwrap();
+    bench_dual_numbers(c, "water_4c_polar", state_pcsaft(parameters));
+
+    // methane, ethane, propane
+    let parameters = PcSaftParameters::from_json(
+        vec!["methane", "ethane", "propane"],
+        "./parameters/pcsaft/gross2001.json",
+        None,
+        IdentifierOption::Name,
+    )
+    .unwrap();
+    bench_dual_numbers(c, "methane_ethane_propane", state_pcsaft(parameters));
+}
+
+/// Benchmark for the PC-SAFT equation of state.
+/// Binary system of methane and co2 used to model biogas.
+fn methane_co2_pcsaft(c: &mut Criterion) {
+    let parameters = PcSaftParameters::from_multiple_json(
+        &[
+            (vec!["methane"], "./parameters/pcsaft/gross2001.json"),
+            (
+                vec!["carbon dioxide"],
+                "./parameters/pcsaft/gross2005_fit.json",
+            ),
+        ],
+        None,
+        IdentifierOption::Name,
+    )
+    .unwrap();
+    let k_ij = -0.0192211646;
+    let parameters =
+        PcSaftParameters::new_binary(parameters.pure_records.clone(), Some(k_ij.into()));
+    let eos = Arc::new(PcSaft::new(Arc::new(parameters)));
+
+    // 230 K, 50 bar, x0 = 0.15
+    let temperature = 230.0 * KELVIN;
+    let density = 24.16896 * KILO * MOL / METER.powi(3);
+    let volume = 10.0 * MOL / density;
+    let x = arr1(&[0.15, 0.85]);
+    let moles = &x * 10.0 * MOL;
+    let state = State::new_nvt(&eos, temperature, volume, &moles).unwrap();
+    bench_dual_numbers(c, "methane_co2", state);
+}
+
+criterion_group!(bench, pcsaft, methane_co2_pcsaft);
+criterion_main!(bench);

benches/state_properties.rs

@@ -0,0 +1,72 @@
+use criterion::{criterion_group, criterion_main, Criterion};
+use feos::pcsaft::{PcSaft, PcSaftParameters};
+use feos_core::{
+    parameter::{IdentifierOption, Parameter},
+    Contributions, EquationOfState, State,
+};
+use ndarray::arr1;
+use quantity::si::*;
+use std::sync::Arc;
+
+type S = State<SIUnit, PcSaft>;
+
+/// Evaluate a property of a state given the EoS, the property to compute,
+/// temperature, volume, moles, and the contributions to consider.
+fn property<E: EquationOfState, T, F: Fn(&State<SIUnit, E>, Contributions) -> T>(
+    (eos, property, t, v, n, contributions): (
+        &Arc<E>,
+        F,
+        SINumber,
+        SINumber,
+        &SIArray1,
+        Contributions,
+    ),
+) -> T {
+    let state = State::new_nvt(eos, t, v, n).unwrap();
+    property(&state, contributions)
+}
+
+/// Evaluate a property with of a state given the EoS, the property to compute,
+/// temperature, volume, moles.
+fn property_no_contributions<E: EquationOfState, T, F: Fn(&State<SIUnit, E>) -> T>(
+    (eos, property, t, v, n): (&Arc<E>, F, SINumber, SINumber, &SIArray1),
+) -> T {
+    let state = State::new_nvt(eos, t, v, n).unwrap();
+    property(&state)
+}
+
+fn properties_pcsaft(c: &mut Criterion) {
+    let parameters = PcSaftParameters::from_json(
+        vec!["methane", "ethane", "propane"],
+        "./parameters/pcsaft/gross2001.json",
+        None,
+        IdentifierOption::Name,
+    )
+    .unwrap();
+    let eos = Arc::new(PcSaft::new(Arc::new(parameters)));
+    let t = 300.0 * KELVIN;
+    let density = 71.18 * KILO * MOL / METER.powi(3);
+    let v = 100.0 * MOL / density;
+    let x = arr1(&[1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0]);
+    let m = &x * 100.0 * MOL;
+
+    let mut group = c.benchmark_group("methane_ethane_propane");
+    group.bench_function("a", |b| {
+        b.iter(|| property((&eos, S::helmholtz_energy, t, v, &m, Contributions::Total)))
+    });
+    group.bench_function("compressibility", |b| {
+        b.iter(|| property((&eos, S::compressibility, t, v, &m, Contributions::Total)))
+    });
+    group.bench_function("ln_phi", |b| {
+        b.iter(|| property_no_contributions((&eos, S::ln_phi, t, v, &m)))
+    });
+    group.bench_function("c_v", |b| {
+        b.iter(|| property((&eos, S::c_v, t, v, &m, Contributions::ResidualNvt)))
+    });
+    group.bench_function("molar_volume", |b| {
+        b.iter(|| property((&eos, S::molar_volume, t, v, &m, Contributions::ResidualNvt)))
+    });
+}
+
+criterion_group!(bench, properties_pcsaft);
+criterion_main!(bench);

feos-core/src/lib.rs

@@ -42,7 +42,7 @@ pub use errors::{EosError, EosResult};
 pub use phase_equilibria::{
     PhaseDiagram, PhaseDiagramHetero, PhaseEquilibrium, SolverOptions, Verbosity,
 };
-pub use state::{Contributions, DensityInitialization, State, StateBuilder, StateHD, StateVec};
+pub use state::{Contributions, DensityInitialization, State, StateBuilder, StateHD, StateVec, Derivative};
 
 #[cfg(feature = "python")]
 pub mod python;

feos-core/src/state/mod.rs

@@ -188,7 +188,7 @@ where
 /// Derivatives of the helmholtz energy.
 #[derive(Clone, Copy, Eq, Hash, PartialEq, Debug, PartialOrd, Ord)]
 #[allow(non_camel_case_types)]
-pub(crate) enum Derivative {
+pub enum Derivative {
     /// Derivative with respect to system volume.
     DV,
     /// Derivative with respect to temperature.
@@ -670,15 +670,17 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
         Err(EosError::NotConverged("State::update_gibbs_energy".into()))
     }
 
-    fn derive0(&self) -> StateHD<f64> {
+    /// Creates a [StateHD] cloning temperature, volume and moles.
+    pub fn derive0(&self) -> StateHD<f64> {
         StateHD::new(
             self.reduced_temperature,
             self.reduced_volume,
             self.reduced_moles.clone(),
         )
     }
 
-    fn derive1(&self, derivative: Derivative) -> StateHD<Dual64> {
+    /// Creates a [StateHD] taking the first derivative.
+    pub fn derive1(&self, derivative: Derivative) -> StateHD<Dual64> {
         let mut t = Dual64::from(self.reduced_temperature);
         let mut v = Dual64::from(self.reduced_volume);
         let mut n = self.reduced_moles.mapv(Dual64::from);
@@ -690,7 +692,12 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
         StateHD::new(t, v, n)
     }
 
-    fn derive2(&self, derivative1: Derivative, derivative2: Derivative) -> StateHD<HyperDual64> {
+    /// Creates a [StateHD] taking the first and second (partial) derivatives.
+    pub fn derive2(
+        &self,
+        derivative1: Derivative,
+        derivative2: Derivative,
+    ) -> StateHD<HyperDual64> {
         let mut t = HyperDual64::from(self.reduced_temperature);
         let mut v = HyperDual64::from(self.reduced_volume);
         let mut n = self.reduced_moles.mapv(HyperDual64::from);
@@ -707,7 +714,8 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
         StateHD::new(t, v, n)
     }
 
-    fn derive3(&self, derivative: Derivative) -> StateHD<Dual3_64> {
+    /// Creates a [StateHD] taking the first, second, and third derivative with respect to a single property.
+    pub fn derive3(&self, derivative: Derivative) -> StateHD<Dual3_64> {
         let mut t = Dual3_64::from(self.reduced_temperature);
         let mut v = Dual3_64::from(self.reduced_volume);
         let mut n = self.reduced_moles.mapv(Dual3_64::from);