added more systems (methane and water) for dual number evaluations, s…

…maller description for README
feos-org · g-bauer · Dec 15, 2022 · Dec 7, 2022 · Dec 7, 2022 · Dec 13, 2022
commit bb93730fa9cb26bab90297441a7770f021426662
diff --git a/benches/README.md b/benches/README.md
@@ -12,5 +12,5 @@ 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. Substances: methane + carbon dioxide. Model: PC-SAFT incl. dipole-dipole contributions.|`pcsaft`|
-|`state_properties`|Properties of `State`. Including state creation using the natural variables of the Helmholtz energy (no density iteration). Substances: methane + ethane + propane. Model: PC-SAFT (hard-sphere, hard-chain, dispersion).|`pcsaft`|
+|`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`|
diff --git a/benches/dual_numbers.rs b/benches/dual_numbers.rs
@@ -3,16 +3,11 @@
 //! These should give an idea about the expected slow-down depending
 //! on the dual number type used without the overhead of the `State`
 //! creation.
-//!
-//! The example system is the binary mixture of CH4/CO2 modelled
-//! with the PCP-SAFT equation of state. The considered Helmholtz
-//! energy contributions are: hard-sphere, hard-chain, dispersion,
-//! and polar (dipolar)..
 use criterion::{criterion_group, criterion_main, Criterion};
 use feos::pcsaft::{PcSaft, PcSaftParameters};
 use feos_core::{
     parameter::{IdentifierOption, Parameter},
-    EquationOfState, HelmholtzEnergy, HelmholtzEnergyDual, StateHD,
+    EquationOfState, HelmholtzEnergy, HelmholtzEnergyDual, State, StateHD,
 };
 use ndarray::arr1;
 use num_dual::{Dual3, Dual64, DualNum, HyperDual64};
@@ -27,7 +22,137 @@ where
     inp.0.evaluate_residual(inp.1)
 }
 
-fn benchmark_dual_numbers(c: &mut Criterion) {
+fn methane_pcsaft(c: &mut Criterion) {
+    let parameters = PcSaftParameters::from_json(
+        vec!["methane"],
+        "./parameters/pcsaft/gross2001.json",
+        None,
+        IdentifierOption::Name,
+    )
+    .unwrap();
+    let eos = Arc::new(PcSaft::new(Arc::new(parameters)));
+
+    let cp = State::critical_point(&eos, None, None, Default::default()).unwrap();
+
+    let t = 0.8 * cp.temperature;
+    let density = cp.density;
+    let volume = 10.0 * MOL / density;
+    let x = arr1(&[1.0]);
+    let m = &x * 10.0 * MOL;
+
+    let mut group = c.benchmark_group("PC-SAFT methane");
+
+    // real valued evaluation
+    let s = StateHD::new(
+        t.to_reduced(KELVIN).unwrap(),
+        volume.to_reduced(ANGSTROM.powi(3)).unwrap(),
+        m.to_reduced(MOL).unwrap(),
+    );
+    group.bench_function("a_f64", |b| b.iter(|| a_res((&eos, &s))));
+
+    // da_dv - dual number
+    let s = StateHD::new(
+        Dual64::from_re(t.to_reduced(KELVIN).unwrap()),
+        Dual64::from_re(volume.to_reduced(ANGSTROM.powi(3)).unwrap()).derive(),
+        m.to_reduced(MOL)
+            .unwrap()
+            .iter()
+            .map(|&mi| Dual64::from_re(mi))
+            .collect(),
+    );
+    group.bench_function("a_dual", |b| b.iter(|| a_res((&eos, &s))));
+
+    // d2a_dv2 - hyperdual number
+    let s = StateHD::new(
+        HyperDual64::from_re(t.to_reduced(KELVIN).unwrap()).derive1(),
+        HyperDual64::from_re(volume.to_reduced(ANGSTROM.powi(3)).unwrap()).derive2(),
+        m.to_reduced(MOL)
+            .unwrap()
+            .iter()
+            .map(|&mi| HyperDual64::from_re(mi))
+            .collect(),
+    );
+    group.bench_function("a_hyperdual", |b| b.iter(|| a_res((&eos, &s))));
+
+    // d3a_dv3 - dual3 number
+    let s = StateHD::new(
+        Dual3::from_re(t.to_reduced(KELVIN).unwrap()),
+        Dual3::from_re(volume.to_reduced(ANGSTROM.powi(3)).unwrap()).derive(),
+        m.to_reduced(MOL)
+            .unwrap()
+            .iter()
+            .map(|&mi| Dual3::from_re(mi))
+            .collect(),
+    );
+    group.bench_function("a_dual3", |b| b.iter(|| a_res((&eos, &s))));
+}
+
+fn water_4c_polar_pcsaft(c: &mut Criterion) {
+    let parameters = PcSaftParameters::from_json(
+        vec!["water_4C_polar"],
+        "./parameters/pcsaft/rehner2020.json",
+        None,
+        IdentifierOption::Name,
+    )
+    .unwrap();
+    let eos = Arc::new(PcSaft::new(Arc::new(parameters)));
+
+    let cp = State::critical_point(&eos, None, None, Default::default()).unwrap();
+
+    let t = 0.8 * cp.temperature;
+    let density = cp.density;
+    let volume = 10.0 * MOL / density;
+    let x = arr1(&[1.0]);
+    let m = &x * 10.0 * MOL;
+
+    let mut group = c.benchmark_group("PC-SAFT water (4C, polar)");
+
+    // real valued evaluation
+    let s = StateHD::new(
+        t.to_reduced(KELVIN).unwrap(),
+        volume.to_reduced(ANGSTROM.powi(3)).unwrap(),
+        m.to_reduced(MOL).unwrap(),
+    );
+    group.bench_function("a_f64", |b| b.iter(|| a_res((&eos, &s))));
+
+    // da_dv - dual number
+    let s = StateHD::new(
+        Dual64::from_re(t.to_reduced(KELVIN).unwrap()),
+        Dual64::from_re(volume.to_reduced(ANGSTROM.powi(3)).unwrap()).derive(),
+        m.to_reduced(MOL)
+            .unwrap()
+            .iter()
+            .map(|&mi| Dual64::from_re(mi))
+            .collect(),
+    );
+    group.bench_function("a_dual", |b| b.iter(|| a_res((&eos, &s))));
+
+    // d2a_dv2 - hyperdual number
+    let s = StateHD::new(
+        HyperDual64::from_re(t.to_reduced(KELVIN).unwrap()).derive1(),
+        HyperDual64::from_re(volume.to_reduced(ANGSTROM.powi(3)).unwrap()).derive2(),
+        m.to_reduced(MOL)
+            .unwrap()
+            .iter()
+            .map(|&mi| HyperDual64::from_re(mi))
+            .collect(),
+    );
+    group.bench_function("a_hyperdual", |b| b.iter(|| a_res((&eos, &s))));
+
+    // d3a_dv3 - dual3 number
+    let s = StateHD::new(
+        Dual3::from_re(t.to_reduced(KELVIN).unwrap()),
+        Dual3::from_re(volume.to_reduced(ANGSTROM.powi(3)).unwrap()).derive(),
+        m.to_reduced(MOL)
+            .unwrap()
+            .iter()
+            .map(|&mi| Dual3::from_re(mi))
+            .collect(),
+    );
+    group.bench_function("a_dual3", |b| b.iter(|| a_res((&eos, &s))));
+}
+
+fn methane_co2_pcsaft(c: &mut Criterion) {
     let parameters = PcSaftParameters::from_multiple_json(
         &[
             (vec!["methane"], "./parameters/pcsaft/gross2001.json"),
@@ -52,7 +177,7 @@ fn benchmark_dual_numbers(c: &mut Criterion) {
     let x = arr1(&[0.15, 0.85]);
     let m = &x * 10.0 * MOL;
 
-    let mut group = c.benchmark_group("dual_numbers");
+    let mut group = c.benchmark_group("PC-SAFT methane + co2");
 
     // real valued evaluation
     let s = StateHD::new(
@@ -99,5 +224,10 @@ fn benchmark_dual_numbers(c: &mut Criterion) {
     group.bench_function("a_dual3", |b| b.iter(|| a_res((&eos, &s))));
 }
 
-criterion_group!(bench, benchmark_dual_numbers);
+criterion_group!(
+    bench,
+    methane_pcsaft,
+    water_4c_polar_pcsaft,
+    methane_co2_pcsaft
+);
 criterion_main!(bench);
diff --git a/benches/state_properties.rs b/benches/state_properties.rs
@@ -113,7 +113,7 @@ fn molar_volume(inp: (&Arc<PcSaft>, SINumber, SINumber, &Array1<f64>)) -> SIArra
     .molar_volume(Contributions::ResidualNvt)
 }
 
-fn benchmark_properties(c: &mut Criterion) {
+fn properties_pcsaft(c: &mut Criterion) {
     let parameters = PcSaftParameters::from_json(
         vec!["methane", "ethane", "propane"],
         "./parameters/pcsaft/gross2001.json",
@@ -128,7 +128,7 @@ fn benchmark_properties(c: &mut Criterion) {
     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("properties");
+    let mut group = c.benchmark_group("PC-SAFT methane + ethane + propane");
     group.bench_function("a", |b| b.iter(|| helmholtz_energy((&eos, t, volume, &m))));
     group.bench_function("compressibility", |b| {
         b.iter(|| compressibility((&eos, t, density, &x)))
@@ -140,5 +140,5 @@ fn benchmark_properties(c: &mut Criterion) {
     });
 }
 
-criterion_group!(bench, benchmark_properties,);
+criterion_group!(bench, properties_pcsaft);
 criterion_main!(bench);