Minor cleanup of the PC-SAFT implementation (#108)

* Minor cleanup of the PC-SAFT implementation * readd association parameters to string repr of PcSaftParameters
feos-org · Jan 12, 2023 · 74b7f49 · 74b7f49
1 parent dca766e
commit 74b7f49
Show file tree

Hide file tree

Showing 4 changed files with 145 additions and 495 deletions.
diff --git a/benches/state_properties.rs b/benches/state_properties.rs
@@ -77,5 +77,38 @@ fn properties_pcsaft(c: &mut Criterion) {
     });
 }
 
-criterion_group!(bench, properties_pcsaft);
+fn properties_pcsaft_polar(c: &mut Criterion) {
+    let parameters = PcSaftParameters::from_json(
+        vec!["acetone", "butanal", "dimethyl ether"],
+        "./parameters/pcsaft/gross2006.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("state_properties_pcsaft_polar");
+    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, properties_pcsaft_polar);
 criterion_main!(bench);
diff --git a/src/association/mod.rs b/src/association/mod.rs
@@ -80,12 +80,14 @@ impl AssociationParameters {
 
         for (i, record) in records.iter().enumerate() {
             if let Some(record) = record.as_ref() {
-                assoc_comp.push(i);
-                sigma_assoc.push(sigma[i]);
-                kappa_ab.push(record.kappa_ab);
-                epsilon_k_ab.push(record.epsilon_k_ab);
-                na.push(record.na.unwrap_or(1.0));
-                nb.push(record.nb.unwrap_or(1.0));
+                if record.kappa_ab > 0.0 && record.epsilon_k_ab > 0.0 {
+                    assoc_comp.push(i);
+                    sigma_assoc.push(sigma[i]);
+                    kappa_ab.push(record.kappa_ab);
+                    epsilon_k_ab.push(record.epsilon_k_ab);
+                    na.push(record.na.unwrap_or(1.0));
+                    nb.push(record.nb.unwrap_or(1.0));
+                }
             }
         }
 
@@ -238,21 +240,14 @@ impl<P> fmt::Display for Association<P> {
 
 impl<P: HardSphereProperties> Association<P> {
     pub fn assoc_site_frac_ab<D: DualNum<f64>>(deltarho: D, na: f64, nb: f64) -> D {
-        if deltarho.re() > f64::EPSILON.sqrt() {
-            (((deltarho * (na - nb) + 1.0).powi(2) + deltarho * nb * 4.0).sqrt()
-                - (deltarho * (nb - na) + 1.0))
-                / (deltarho * na * 2.0)
-        } else {
-            D::one() + deltarho * nb * (deltarho * (nb + na) - 1.0)
-        }
+        (((deltarho * (na - nb) + 1.0).powi(2) + deltarho * nb * 4.0).sqrt()
+            + (deltarho * (nb - na) + 1.0))
+            .recip()
+            * (2.0 * nb / na)
     }
 
     pub fn assoc_site_frac_a<D: DualNum<f64>>(deltarho: D, na: f64) -> D {
-        if deltarho.re() > f64::EPSILON.sqrt() {
-            ((deltarho * na * 4.0 + 1.0).powi(2) - 1.0).sqrt() / (deltarho * na * 2.0)
-        } else {
-            D::one() + deltarho * na * (deltarho * na * 2.0 - 1.0)
-        }
+        ((deltarho * 4.0 * na + 1.0).sqrt() + 1.0).recip() * 2.0
     }
 
     #[allow(clippy::too_many_arguments)]