From ae96e48b4f5f6c072aa7dcaac527bc5e5f750865 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Thu, 18 Apr 2013 22:07:00 +0000
Subject: [PATCH] [Thermo] Adjust NASA polynomials to eliminate discontinuities

We expect the high- and low-temperature NASA polynomials for each species to be
continuous in a number of places, despite many frequently used thermo databases
having data that is inconsistent.

In addition to warning about discontinuities, We now modify the provided
coefficients to ensure continuity at the midpoint temperature. This resolves
numerical issues, e.g. with the algorithm used for setting the state of a
reactor network.
---
 src/thermo/NasaThermo.h             | 17 +++++++++++------
 src/thermo/SpeciesThermoFactory.cpp | 27 +++++++++++++++++++++++++--
 2 files changed, 36 insertions(+), 8 deletions(-)

diff --git a/src/thermo/NasaThermo.h b/src/thermo/NasaThermo.h
index 618421d4b..b07a747fc 100644
--- a/src/thermo/NasaThermo.h
+++ b/src/thermo/NasaThermo.h
@@ -162,7 +162,7 @@ public:
         vector_fp chigh(c+8, c+15);
         vector_fp clow(c+1, c+8);
 
-        checkContinuity(name, tmid, &clow[0], &chigh[0]);
+        ensureContinuity(name, tmid, &clow[0], &chigh[0]);
 
         m_high[igrp-1].push_back(NasaPoly1(index, tmid, thigh,
                                            refPressure, &chigh[0]));
@@ -460,17 +460,22 @@ protected:
     mutable std::map<size_t, std::string> m_name;
 
 private:
-    //! see SpeciesThermoFactory.cpp for the definition
+    //! Adjust polynomials to be continuous at the midpoint temperature.
     /*!
+     * Check to see if the provided coefficients are nearly continuous. Adjust
+     * the values to get more precise contintinuity to avoid convergence
+     * issues with algorithms that expect these quantities to be continuous.
+     * See SpeciesThermoFactory.cpp for the definition.
+     *
      * @param name string name of species
      * @param tmid  Mid temperature, between the two temperature regions
      * @param clow  coefficients for lower temperature region
      * @param chigh coefficients for higher temperature region
      */
-    void checkContinuity(const std::string& name, double tmid,
-                         const doublereal* clow, doublereal* chigh);
+    void ensureContinuity(const std::string& name, double tmid,
+                         doublereal* clow, doublereal* chigh);
 
-    //! for internal use by checkContinuity
+    //! for internal use by ensureContinuity
     /*!
      * @param t temperature
      * @param c coefficient array
@@ -481,7 +486,7 @@ private:
                + c[0]/t;
     }
 
-    //! for internal use by checkContinuity
+    //! for internal use by ensureContinuity
     /*!
      * @param t temperature
      * @param c coefficient array
diff --git a/src/thermo/SpeciesThermoFactory.cpp b/src/thermo/SpeciesThermoFactory.cpp
index 5e6bf3d70..4cb9ee6b0 100644
--- a/src/thermo/SpeciesThermoFactory.cpp
+++ b/src/thermo/SpeciesThermoFactory.cpp
@@ -198,8 +198,8 @@ SpeciesThermo* SpeciesThermoFactory::newSpeciesThermoManager(std::string& stype)
     return (SpeciesThermo*) 0;
 }
 
-void NasaThermo::checkContinuity(const std::string& name, double tmid,
-                                 const doublereal* clow, doublereal* chigh)
+void NasaThermo::ensureContinuity(const std::string& name, double tmid,
+                                 doublereal* clow, doublereal* chigh)
 {
     // heat capacity
     doublereal cplow = poly4(tmid, clow + 2);
@@ -215,6 +215,13 @@ void NasaThermo::checkContinuity(const std::string& name, double tmid,
                  +fp2str(cphigh)+".\n");
     }
 
+    // Adjust coefficients to eliminate any discontinuity
+    chigh[2] += 0.5 * delta;
+    clow[2] -= 0.5 * delta;
+
+    AssertThrowMsg(std::abs(poly4(tmid, clow+2) - poly4(tmid, chigh+2)) < 1e-12,
+                   "NasaThermo::ensureContinuity", "Cp/R does not match");
+
     // enthalpy
     doublereal hrtlow = enthalpy_RT(tmid, clow);
     doublereal hrthigh = enthalpy_RT(tmid, chigh);
@@ -229,6 +236,14 @@ void NasaThermo::checkContinuity(const std::string& name, double tmid,
                  +fp2str(hrthigh)+".\n");
     }
 
+    // Adjust coefficients to eliminate any discontinuity
+    chigh[0] += 0.5 * delta * tmid;
+    clow[0] -= 0.5 * delta * tmid;
+
+    AssertThrowMsg(std::abs(enthalpy_RT(tmid, clow) -
+                            enthalpy_RT(tmid, chigh)) < 1e-12,
+                   "NasaThermo::ensureContinuity", "H/RT does not match");
+
     // entropy
     doublereal srlow = entropy_R(tmid, clow);
     doublereal srhigh = entropy_R(tmid, chigh);
@@ -242,6 +257,14 @@ void NasaThermo::checkContinuity(const std::string& name, double tmid,
         writelog("\tValue computed using high-temperature polynomial: "
                  +fp2str(srhigh)+".\n");
     }
+
+    // Adjust coefficients to eliminate any discontinuity
+    chigh[1] += 0.5 * delta;
+    clow[1] -= 0.5 * delta;
+
+    AssertThrowMsg(std::abs(entropy_R(tmid, clow) -
+                            entropy_R(tmid, chigh)) < 1e-12,
+                   "NasaThermo::ensureContinuity", "S/R does not match");
 }