Fix handling of reference pressure in MultiSpeciesThermo

Require all species in a phase to have the same reference pressure. Fixes #1435
Cantera · Jun 14, 2023 · 5973dbd · 5973dbd
1 parent e2cdb1c
commit 5973dbd
Show file tree

Hide file tree

Showing 4 changed files with 30 additions and 5 deletions.
diff --git a/doc/sphinx/yaml/species.rst b/doc/sphinx/yaml/species.rst
@@ -160,6 +160,7 @@ Example::
     thermo:
       model: Shomate
       temperature-ranges: [298, 1300, 6000]
+      reference-pressure: 1 bar
       data:
       - [25.56759, 6.096130, 4.054656, -2.671301, 0.131021,
         -118.0089, 227.3665]

diff --git a/include/cantera/thermo/MultiSpeciesThermo.h b/include/cantera/thermo/MultiSpeciesThermo.h
@@ -126,6 +126,8 @@ class MultiSpeciesThermo
      * for the first species.
      *
      * @param k Species Index
+     * @deprecated The species index parameter is deprecated and will be removed after
+     *     Cantera 3.0. All species in a phase must have the same reference pressure.
      */
     virtual doublereal refPressure(size_t k=npos) const;
 
@@ -226,7 +228,7 @@ class MultiSpeciesThermo
     double m_thigh_min = 1e+30;
 
     //! reference pressure (Pa)
-    double m_p0 = OneAtm;
+    double m_p0 = 0.0;
 
     //! indicates if data for species has been installed
     std::vector<bool> m_installed;

diff --git a/src/thermo/MultiSpeciesThermo.cpp b/src/thermo/MultiSpeciesThermo.cpp
@@ -13,12 +13,13 @@
 #include "cantera/base/stringUtils.h"
 #include "cantera/base/utilities.h"
 #include "cantera/base/ctexceptions.h"
+#include "cantera/base/global.h"
 
 namespace Cantera
 {
 
 void MultiSpeciesThermo::install_STIT(size_t index,
-                                        shared_ptr<SpeciesThermoInterpType> stit_ptr)
+                                      shared_ptr<SpeciesThermoInterpType> stit_ptr)
 {
     if (!stit_ptr) {
         throw CanteraError("MultiSpeciesThermo::install_STIT",
@@ -27,6 +28,15 @@ void MultiSpeciesThermo::install_STIT(size_t index,
     AssertThrowMsg(m_speciesLoc.find(index) == m_speciesLoc.end(),
             "MultiSpeciesThermo::install_STIT",
             "Index position isn't null, duplication of assignment: {}", index);
+    if (m_p0 == 0) {
+        // First species added; use this to set the reference pressure
+        m_p0 = stit_ptr->refPressure();
+    } else if (fabs(m_p0 - stit_ptr->refPressure()) > 1e-6) {
+        throw CanteraError("MultiSpeciesThermo::install_STIT",
+            "Cannot add species {} with reference pressure {}.\n"
+            "Inconsistent with previously-added species with reference pressure {}.",
+            index, stit_ptr->refPressure(), m_p0);
+    }
     int type = stit_ptr->reportType();
     m_speciesLoc[index] = {type, m_sp[type].size()};
     m_sp[type].emplace_back(index, stit_ptr);
@@ -144,6 +154,9 @@ doublereal MultiSpeciesThermo::maxTemp(size_t k) const
 doublereal MultiSpeciesThermo::refPressure(size_t k) const
 {
     if (k != npos) {
+        warn_deprecated("MultiSpeciesThermo::refPressure(size_t k)",
+            "The species index parameter is deprecated and will be removed after"
+            " Cantera 3.0.");
         const SpeciesThermoInterpType* sp = provideSTIT(k);
         if (sp) {
             return sp->refPressure();

diff --git a/test/thermo/thermoParameterizations.cpp b/test/thermo/thermoParameterizations.cpp
@@ -54,10 +54,19 @@ TEST_F(SpeciesThermoInterpTypeTest, install_const_cp)
     EXPECT_DOUBLE_EQ(p2.cp_mass(), p.cp_mass());
 }
 
-TEST_F(SpeciesThermoInterpTypeTest, DISABLED_install_bad_pref)
+TEST_F(SpeciesThermoInterpTypeTest, install_specified_pref)
+{
+    auto sO2 = make_shared<Species>("O2", parseCompString("O:2"));
+    auto sH2 = make_shared<Species>("H2", parseCompString("H:2"));
+    sO2->thermo = make_shared<ConstCpPoly>(200, 5000, 100000, c_o2);
+    sH2->thermo = make_shared<ConstCpPoly>(200, 5000, 100000, c_h2);
+    p.addSpecies(sO2);
+    // Pref does not match
+    EXPECT_DOUBLE_EQ(p.refPressure(), 1e5);
+}
+
+TEST_F(SpeciesThermoInterpTypeTest, install_bad_pref)
 {
-    // Currently broken because MultiSpeciesThermo does not enforce reference
-    // pressure consistency.
     auto sO2 = make_shared<Species>("O2", parseCompString("O:2"));
     auto sH2 = make_shared<Species>("H2", parseCompString("H:2"));
     sO2->thermo = make_shared<ConstCpPoly>(200, 5000, 101325, c_o2);