From 45d50999794d3426670a49e7d27efaf5d54dc4a0 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Wed, 23 Jan 2019 18:38:07 -0500
Subject: [PATCH] Create PDSS_HKFT objects from YAML definitions

---
 src/thermo/PDSS_HKFT.cpp       | 25 ++++++++++++++
 test/data/thermo-models.yaml   | 60 +++++++++++++++++++++++++++++++++-
 test/thermo/thermoFromYaml.cpp | 25 ++++++++++++++
 3 files changed, 109 insertions(+), 1 deletion(-)

diff --git a/src/thermo/PDSS_HKFT.cpp b/src/thermo/PDSS_HKFT.cpp
index 3b2b743f5..31aa9df52 100644
--- a/src/thermo/PDSS_HKFT.cpp
+++ b/src/thermo/PDSS_HKFT.cpp
@@ -236,6 +236,31 @@ void PDSS_HKFT::setState_TP(doublereal temp, doublereal pres)
 void PDSS_HKFT::initThermo()
 {
     PDSS::initThermo();
+    if (m_input.hasKey("h0")) {
+        m_deltaH_formation_tr_pr = m_input.convert("h0", "cal/gmol");
+    }
+    if (m_input.hasKey("g0")) {
+        m_deltaG_formation_tr_pr = m_input.convert("g0", "cal/gmol");
+    }
+    if (m_input.hasKey("s0")) {
+        m_Entrop_tr_pr = m_input.convert("s0", "cal/gmol/K");
+    }
+    auto& units = m_input.units();
+    if (m_input.hasKey("a")) {
+        auto& a = m_input["a"].asVector<AnyValue>(4);
+        m_a1 = units.convert(a[0], "cal/gmol/bar");
+        m_a2 = units.convert(a[1], "cal/gmol");
+        m_a3 = units.convert(a[2], "cal*K/gmol/bar");
+        m_a4 = units.convert(a[3], "cal*K/gmol");
+    }
+    if (m_input.hasKey("c")) {
+        auto& c = m_input["c"].asVector<AnyValue>(2);
+        m_c1 = units.convert(c[0], "cal/gmol/K");
+        m_c2 = units.convert(c[1], "cal*K/gmol");
+    }
+    if (m_input.hasKey("omega")) {
+        m_omega_pr_tr = m_input.convert("omega", "cal/gmol");
+    }
 
     // Ok, if we are missing one, then we construct its value from the other two.
     // This code has been internally verified.
diff --git a/test/data/thermo-models.yaml b/test/data/thermo-models.yaml
index 503e89657..47d4e1e1d 100644
--- a/test/data/thermo-models.yaml
+++ b/test/data/thermo-models.yaml
@@ -129,7 +129,7 @@ phases:
   activity-data:
     temperature-model: complex # "constant" or "linear" are the other options
     A_Debye: 1.175930 kg^0.5/gmol^0.5
-    interactions:
+    interactions: &hmw-ineractions
     - species: [Na+, Cl-]
       beta0: [0.0765, 0.008946, -3.3158E-6, -777.03, -4.4706]
       beta1: [0.2664, 6.1608E-5, 1.0715E-6, 0.0, 0.0]
@@ -154,6 +154,19 @@ phases:
     - {species: [Na+, H+], theta: 0.036}
     - {species: [Cl-, Na+, H+], psi: [-0.004]}
 
+- name: HMW-NaCl-HKFT
+  species:
+  - {HMW-species: [H2O(L)]}
+  - {HKFT-species: [Na+, Cl-, H+, OH-]}
+  thermo: HMW-electrolyte
+  state:
+    T: 323.15
+    P: 101325
+    molalities: {Na+: 6.0954, Cl-: 6.0954, H+: 2.1628e-9, OH-: 1.3977e-6}
+  activity-data:
+    temperature-model: complex # "constant" or "linear" are the other options
+    A_Debye: variable
+    interactions: *hmw-ineractions
 
 species:
 - name: NaCl(s)
@@ -475,3 +488,48 @@ HMW-species:
   equation-of-state:
     model: constant-volume
     molar-volume: 1.3
+
+
+HKFT-species:
+- name: Na+
+  composition: {Na: 1, E: -1}
+  equation-of-state:
+    model: HKFT
+    h0: -57433. cal/gmol
+    s0: 13.96 cal/gmol/K
+    a: [0.1839 cal/gmol/bar, -228.5 cal/gmol, 3.256 cal*K/gmol/bar, -27260. cal*K/gmol]
+    c: [18.18 cal/gmol/K, -29810. cal*K/gmol]
+    omega: 33060 cal/gmol
+
+- name: Cl-
+  composition: {Cl: 1, E: 1}
+  equation-of-state:
+    model: HKFT
+    units: {energy: cal, quantity: mol, pressure: bar}
+    g0: -31379
+    s0: 13.56
+    a: [0.4032, 480.1, 5.563, -28470.]
+    c: [-4.4, -57140.]
+    omega: 145600.
+
+- name: H+
+  composition: {H: 1, E: -1}
+  equation-of-state:
+    model: HKFT
+    g0: 0
+    s0: 0
+    a: [0, 0, 0, 0]
+    c: [0, 0]
+    omega: 0
+
+- name: OH-
+  composition: {O: 1, H: 1, E: 1}
+  equation-of-state:
+    model: HKFT
+    units: {energy: cal, quantity: mol, pressure: bar}
+    g0: -37595.
+    h0: -54977.
+    s0: -2.56
+    a: [0.12527, 7.38, 1.8423, -27821]
+    c: [4.15, -103460.]
+    omega: 172460.
diff --git a/test/thermo/thermoFromYaml.cpp b/test/thermo/thermoFromYaml.cpp
index 5e69af887..a0c140483 100644
--- a/test/thermo/thermoFromYaml.cpp
+++ b/test/thermo/thermoFromYaml.cpp
@@ -302,3 +302,28 @@ TEST(ThermoFromYaml, HMWSoln)
         EXPECT_NEAR(mu0[k]/1e6, mu0Ref[k], 2e-6);
     }
 }
+
+TEST(ThermoFromYaml, HMWSoln_HKFT)
+{
+    AnyMap infile = AnyMap::fromYamlFile("thermo-models.yaml");
+    auto phaseNodes = infile["phases"].asMap("name");
+    auto thermo = newPhase(*phaseNodes.at("HMW-NaCl-HKFT"), infile);
+
+    double mvRef[] = {0.01815224, 0.00157182, 0.01954605, 0.00173137, -0.0020266};
+    double hRef[] = {-2.84097589e+08, -2.38159643e+08, -1.68846908e+08,
+                     3.59728865e+06, -2.29291570e+08};
+    double acoeffRef[] = {0.922402064, 1.21860196, 1.21860175, 5.08172471,
+                          0.59832209};
+
+    // Regression test based on HMWSoln.fromScratch_HKFT
+    size_t N = thermo->nSpecies();
+    vector_fp mv(N), h(N), acoeff(N);
+    thermo->getPartialMolarVolumes(mv.data());
+    thermo->getPartialMolarEnthalpies(h.data());
+    thermo->getActivityCoefficients(acoeff.data());
+    for (size_t k = 0; k < N; k++) {
+        EXPECT_NEAR(mv[k], mvRef[k], 2e-8);
+        EXPECT_NEAR(h[k], hRef[k], 2e0);
+        EXPECT_NEAR(acoeff[k], acoeffRef[k], 2e-8);
+    }
+}