Add BinarySolutionTabulatedThermo class

include/cantera/thermo/BinarySolutionTabulatedThermo.h

@@ -0,0 +1,110 @@
+/**
+ * @file BinarySolutionTabulatedThermo.h
+ * Header file for an binary solution model with tabulated standard state
+ * thermodynamic data (see \ref thermoprops and class
+ * \link Cantera::BinarySolutionTabulatedThermo BinarySolutionTabulatedThermo\endlink).
+ */
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://www.cantera.org/license.txt for license and copyright information.
+
+#ifndef CT_BINARYSOLUTIONTABULATEDTHERMO_H
+#define CT_BINARYSOLUTIONTABULATEDTHERMO_H
+
+#include "IdealSolidSolnPhase.h"
+#include "cantera/base/utilities.h"
+
+namespace Cantera
+{
+
+//!   Overloads the virtual methods of class IdealSolidSolnPhase to implement the
+//!   tabulated thermodynamics for one species.
+/**
+ *
+ *
+ * @ingroup thermoprops
+ */
+class BinarySolutionTabulatedThermo : public IdealSolidSolnPhase
+{
+public:
+    /**
+     * Constructor for BinarySolutionTabulatedThermo.
+     * The generalized concentrations can have three different forms
+     * depending on the value of the member attribute #m_formGC, which
+     * is supplied in the constructor or read from the XML data file.
+     *
+     * @param formCG This parameter initializes the #m_formGC variable.
+     */
+    BinarySolutionTabulatedThermo(int formCG=0);
+
+    //! Construct and initialize an BinarySolutionTabulatedThermo ThermoPhase object
+    //! directly from an ASCII input file
+    /*!
+     * This constructor will also fully initialize the object.
+     * The generalized concentrations can have three different forms
+     * depending on the value of the member attribute #m_formGC, which
+     * is supplied in the constructor or read from the XML data file.
+     *
+     * @param infile File name for the XML datafile containing information
+     *               for this phase
+     * @param id     The name of this phase. This is used to look up
+     *               the phase in the XML datafile.
+     * @param formCG This parameter initializes the #m_formGC variable.
+     */
+    BinarySolutionTabulatedThermo(const std::string& infile, const std::string& id="", int formCG=0);
+
+    //! Construct and initialize an BinarySolutionTabulatedThermo ThermoPhase object
+    //! directly from an XML database
+    /*!
+     * The generalized concentrations can have three different forms
+     * depending on the value of the member attribute #m_formGC, which
+     * is supplied in the constructor and/or read from the data file.
+     *
+     * @param root   XML tree containing a description of the phase.
+     *               The tree must be positioned at the XML element
+     *               named phase with id, "id", on input to this routine.
+     * @param id     The name of this phase. This is used to look up
+     *               the phase in the XML datafile.
+     * @param formCG This parameter initializes the #m_formGC variable.
+     */
+    BinarySolutionTabulatedThermo(XML_Node& root, const std::string& id="", int formCG=0);
+
+    virtual std::string type() const {
+        return "BinarySolutionTabulatedThermo";
+    }
+
+    virtual void initThermoXML(XML_Node& phaseNode, const std::string& id_);
+
+protected:
+
+    int m_formGC;
+
+    double m_Pref;
+
+    double m_Pcurrent;
+
+    vector_fp m_speciesMolarVolume;
+
+    //! If the compositions have changed, update the tabulated thermo lookup
+    virtual void compositionChanged();
+
+    //! Species thermodynamics interpolation functions
+    double* interpolate(double x) const;
+
+    //! Current tabulated species index
+    size_t m_kk_tab;
+
+    //! Current tabulated species mole fraction
+    double m_xlast;
+
+    //! Vector for storing tabulated thermo
+    vector_fp m_molefrac_tab;
+    vector_fp m_enthalpy_tab;
+    vector_fp m_entropy_tab;
+
+private:
+    void _updateThermo();
+};
+}
+
+#endif

interfaces/cython/cantera/ctml_writer.py

@@ -22,6 +22,7 @@
 from __future__ import print_function
 
 import sys
+import re
 
 # Python 2/3 compatibility
 try:
@@ -267,13 +268,13 @@ def disable_motz_wise():
     global _motz_wise
     _motz_wise = False
 
-def export_species(filename, fmt = 'CSV'):
+def export_species(filename, fmt='CSV'):
     global _valexport
     global _valfmt
     _valexport = filename
     _valfmt = fmt
 
-def validate(species = 'yes', reactions = 'yes'):
+def validate(species='yes', reactions='yes'):
     """
     Enable or disable validation of species and reactions.
 
@@ -527,6 +528,7 @@ class species(object):
                  note = '',
                  thermo = None,
                  transport = None,
+                 standardState = None,
                  charge = -999,
                  size = 1.0):
         """
@@ -555,6 +557,9 @@ class species(object):
             and must be consistent with the transport model of the phase into which
             the species is imported. To specify parameters for multiple
             transport models, group the entries in parentheses.
+        :param standardState:
+            The species standard state model. Currently used only for IdealSolidSolution and derived
+            classes where it is used to calculate the phase density.
         :param size:
             The species "size". Currently used only for surface species,
             where it represents the number of sites occupied.
@@ -574,6 +579,7 @@ class species(object):
             self._thermo = const_cp()
 
         self._transport = transport
+        self._standardState = standardState
         chrg = 0
         self._charge = charge
         if 'E' in self._atoms:
@@ -590,7 +596,7 @@ class species(object):
         for e in self._atoms.keys():
             _enames[e] = 1
 
-    def export(self, f, fmt = 'CSV'):
+    def export(self, f, fmt='CSV'):
         if fmt == 'CSV':
             s = self._name+','
             for e in _enames:
@@ -639,12 +645,21 @@ class species(object):
                 nt = len(self._transport)
                 for n in range(nt):
                     self._transport[n].build(t)
+        if self._standardState:
+            ss = s.addChild("standardState")
+            ss['model'] = id
+            if isinstance(self._standardState, standardState):
+                self._standardState.build(ss)
+            else:
+                nt = len(self._thermo)
+                for n in range(nt):
+                    self._thermo[n].build(t)
 
 class thermo(object):
     """Base class for species standard-state thermodynamic properties."""
     def _build(self, p):
         return p.addChild("thermo")
-    def export(self, f, fmt = 'CSV'):
+    def export(self, f, fmt='CSV'):
         pass
 
 class Mu0_table(thermo):
@@ -798,6 +813,24 @@ class NASA9(thermo):
         u["size"] = "9"
         u["name"] = "coeffs"
 
+class standardState(object):
+    """Base class for species standard-state properties."""
+    def _build(self, p):
+        return p.addChild("standardState")
+
+class constantIncompressible(standardState):
+    """Constant molar volume."""
+    def __init__(self,
+                 molarVolume = 0.0):
+        """
+        :param molarVolume:
+            Reference-state molar volume. Default: 0.0.
+        """
+        self._mv = molarVolume
+    def build(self, ss):
+        ss['model'] = 'constant_incompressible'
+        mv_units = _ulen+'3/'+_umol
+        addFloat(ss,'molarVolume',self._mv, defunits = mv_units)
 
 class activityCoefficients(object):
     pass
@@ -1268,6 +1301,8 @@ class reaction(object):
             self._kf = [self._kf]
         elif self._type == 'edge':
             self._kf = [self._kf]
+            if self._rateCoeff:
+                kfnode['type'] = self._rateCoeff
         elif self._type == 'threeBody':
             self._kf = [self._kf]
             self.mdim += 1
@@ -1627,11 +1662,13 @@ class edge_reaction(reaction):
                  kf = None,
                  id = '',
                  order = '',
+                 rateCoeff = '',
                  beta = 0.0,
                  options = []):
         reaction.__init__(self, equation, kf, id, order, options)
         self._type = 'edge'
         self._beta = beta
+        self._rateCoeff = rateCoeff
 
 
 #--------------
@@ -2126,6 +2163,161 @@ class incompressible_solid(phase):
         k = ph.addChild("kinetics")
         k['model'] = 'none'
 
+class IdealSolidSolution(phase):
+    """An IdealSolidSolution phase."""
+    def __init__(self,
+                 name = '',
+                 elements = '',
+                 species = '',
+                 note = '',
+                 density = None,
+                 transport = 'None',
+                 initial_state = None,
+                 standard_concentration = None,
+                 options = []):
+
+        phase.__init__(self, name, 3, elements, species, note, 'None',
+                       initial_state, options)
+        self._pure = 0
+        self._stdconc = standard_concentration
+        if self._stdconc is None:
+            raise CTI_Error('In phase ' + name + ': standard_concentration must be specified.')
+        self._tr = transport
+
+    def conc_dim(self):
+        return (1,-3)
+
+    def build(self, p):
+        ph = phase.build(self, p)
+        e = ph.child("thermo")
+        e['model'] = 'IdealSolidSolution'
+        if self._tr:
+            t = ph.addChild('transport')
+            t['model'] = self._tr
+        k = ph.addChild("kinetics")
+        k['model'] = 'none'
+        sc = ph.addChild('standardConc')
+        sc['model'] = self._stdconc
+
+class BinarySolutionTabulatedThermo(phase):
+    """A BinarySolutionTabulatedThermo phase."""
+    def __init__(self,
+                 name = '',
+                 elements = '',
+                 species = '',
+                 note = '',
+                 transport = 'None',
+                 initial_state = None,
+                 standard_concentration = None,
+                 tabulated_species = '',
+                 tabulated_thermo = None,
+                 options = []):
+        phase.__init__(self, name, 3, elements, species, note, 'None',
+                       initial_state, options)
+        self._pure = 0
+        self._tabSpecies = tabulated_species
+        self._tabThermo = tabulated_thermo
+        self._stdconc = standard_concentration
+        self._tr = transport
+        if self._stdconc is None:
+            raise CTI_Error('In phase ' + name
+                + ': standard_concentration must be specified.')
+        if tabulated_species is None:
+            raise CTI_Error('In phase ' + name
+                + ': tabulated_species must be specified.')
+        if tabulated_thermo is None:
+            raise CTI_Error('In phase ' + name
+                + ': Thermo data must be provided for the tabulated_species.')
+
+    def conc_dim(self):
+        return (1,-3)
+    def build(self, p):
+        ph = phase.build(self, p)
+        e = ph.child("thermo")
+        e['model'] = 'BinarySolutionTabulatedThermo'
+        e1 = e.addChild('tabulatedSpecies')
+        e1['name'] = self._tabSpecies
+        t = e.addChild("tabulatedThermo")
+        self._tabThermo.build(t)
+        if self._tr:
+            t = ph.addChild('transport')
+            t['model'] = self._tr
+        k = ph.addChild("kinetics")
+        k['model'] = 'none'
+        sc = ph.addChild('standardConc')
+        sc['model'] = self._stdconc
+
+class table(thermo):
+    """User provided thermo table for BinarySolutionTabulatedThermo"""
+    def __init__(self,
+                 moleFraction = ([],''),
+                 enthalpy = ([],''),
+                 entropy = ([],'')):
+        """
+        :param moleFraction:
+            The mole fraction of the tabulated species. Required parameter.
+        :param enthalpy:
+            The enthalpy of the tabulated species. Required parameter.
+        :param entropy:
+            The entropy of the tabulated species. Required parameter.
+        """
+        self.x = moleFraction
+        self.h = enthalpy
+        self.s = entropy
+    def build(self,t):
+        x = ', '.join('{0:12.5e}'.format(val) for val in self.x[0])
+        u1 = t.addChild("moleFraction", x)
+        u1['units'] = self.x[1]
+        u1['size'] = str(len(self.x[0]))
+        h = ', '.join('{0:12.5e}'.format(val) for val in self.h[0])
+        u2 = t.addChild("enthalpy", h)
+        u2['units'] = self.h[1]
+        u2['size'] = str(len(self.h[0]))
+        s = ', '.join('{0:12.5e}'.format(val) for val in self.s[0])
+        u3 = t.addChild("entropy", s)
+        u3['units'] = self.s[1]
+        u3['size'] = str(len(self.s[0]))
+
+class csvfile(thermo):
+    """User provided CSV file for BinarySolutionTabulatedThermo"""
+    def __init__(self,filename):
+        fh = open(filename)
+        x = []
+        linenumber = 0
+        for line in fh.readlines():
+            linenumber += 1
+            line = line.strip()
+            if not line.startswith("*"):
+                value = re.split(r';|,\s|\s+|\t',line)
+                if len(value) != 3:
+                    raise CTI_Error('In file: ' + filename + ', bad line format at line:' + str(value))
+                else:
+                    y = [float(val) for val in value]
+                x.append(y)
+        fh.close()
+        dat = []
+        for i in range(3):
+            dat.append([row[i] for row in x])
+            self.length = len(dat[0])
+            self.dat = dat
+
+    def build(self,t):
+        energy_units = _uenergy + '/' + 'mol'
+        dat_str = ['','','']
+        nr = 0
+        for rows in self.dat:
+            dat_str[nr] += ', '.join('{0:12.5e}'.format(val) for val in rows)
+            dat_str[nr] += '\n'
+            nr += 1
+        u1 = t.addChild("moleFraction", dat_str[0])
+        u1['units'] = str(1)
+        u1['size'] = str(self.length)
+        u2 = t.addChild("enthalpy", dat_str[1])
+        u2['units'] = energy_units
+        u2['size'] = str(self.length)
+        u3 = t.addChild("entropy", dat_str[2])
+        u3['units'] = energy_units + '/K'
+        u3['size'] = str(self.length)
 
 class lattice(phase):
     def __init__(self,

src/thermo/BinarySolutionTabulatedThermo.cpp

@@ -0,0 +1,226 @@
+/**
+ *  @file BinarySolutionTabulatedThermo.cpp Implementation file for an binary solution model
+ *      with tabulated standard state thermodynamic data (see \ref thermoprops and
+ *      class \link Cantera::BinarySolutionTabulatedThermo BinarySolutionTabulatedThermo\endlink).
+ */
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://www.cantera.org/license.txt for license and copyright information.
+
+#include "cantera/thermo/BinarySolutionTabulatedThermo.h"
+#include "cantera/thermo/PDSS.h"
+#include "cantera/thermo/ThermoFactory.h"
+#include "cantera/base/stringUtils.h"
+#include "cantera/base/ctml.h"
+#include "cantera/thermo/SpeciesThermoFactory.h"
+#include "cantera/thermo/MultiSpeciesThermo.h"
+
+namespace Cantera
+{
+
+BinarySolutionTabulatedThermo::BinarySolutionTabulatedThermo(int formGC) :
+    m_formGC(formGC),
+    m_Pref(OneAtm),
+    m_Pcurrent(OneAtm)
+{
+    if (formGC < 0 || formGC > 2) {
+        throw CanteraError(" BinarySolutionTabulatedThermo Constructor",
+                           " Illegal value of formGC");
+    }
+}
+
+BinarySolutionTabulatedThermo::BinarySolutionTabulatedThermo(const std::string& inputFile,
+        const std::string& id_, int formGC) :
+    m_formGC(formGC),
+    m_Pref(OneAtm),
+    m_Pcurrent(OneAtm)
+{
+    if (formGC < 0 || formGC > 2) {
+        throw CanteraError(" BinarySolutionTabulatedThermo Constructor",
+                           " Illegal value of formGC");
+    }
+    initThermoFile(inputFile, id_);
+}
+
+BinarySolutionTabulatedThermo::BinarySolutionTabulatedThermo(XML_Node& root, const std::string& id_,
+        int formGC) :
+    m_formGC(formGC),
+    m_Pref(OneAtm),
+    m_Pcurrent(OneAtm)
+{
+    if (formGC < 0 || formGC > 2) {
+        throw CanteraError(" BinarySolutionTabulatedThermo Constructor",
+                           " Illegal value of formGC");
+    }
+    importPhase(root, this);
+}
+
+void BinarySolutionTabulatedThermo::compositionChanged()
+{
+    IdealSolidSolnPhase::compositionChanged();
+    _updateThermo();
+}
+
+void BinarySolutionTabulatedThermo::_updateThermo()
+{
+    double tnow = temperature();
+    double xnow = moleFraction(m_kk_tab);
+    double c[4];
+    double *d;
+    double dS_corr = 0.0;
+    double tlow = 0.0, thigh = 0.0;
+    int type = 0;
+    if (m_tlast != tnow || m_xlast != xnow) {
+        c[0] = tnow;
+        d = interpolate(xnow);
+        c[1] = d[0] * 1e3; // 1e3 for conversion J/mol -> J/kmol
+        if (xnow == 0)
+        {
+            dS_corr = -BigNumber;
+        } else if (xnow == 1)
+        {
+            dS_corr = BigNumber;
+        } else
+        {
+            dS_corr = GasConstant*std::log(xnow/(1.0-xnow)) + GasConstant/Faraday*std::log(this->standardConcentration(1-m_kk_tab)/this->standardConcentration(m_kk_tab));
+        }
+        c[2] = d[1] * 1e3 + dS_corr; // 1e3 for conversion J/K/mol -> J/K/kmol
+        c[3] = 0.0;
+        type = m_spthermo.reportType(m_kk_tab);
+        tlow = m_spthermo.minTemp(m_kk_tab);
+        thigh = m_spthermo.maxTemp(m_kk_tab);
+        shared_ptr<SpeciesThermoInterpType> stit(
+                newSpeciesThermoInterpType(type, tlow, thigh, OneAtm, c));
+        m_spthermo.modifySpecies(m_kk_tab, stit);
+        // Update the thermodynamic functions of the reference state.
+        m_spthermo.update(tnow, m_cp0_R.data(), m_h0_RT.data(), m_s0_R.data());
+        doublereal rrt = 1.0 / RT();
+        for (size_t k = 0; k < m_kk; k++) {
+            double deltaE = rrt * m_pe[k];
+            m_h0_RT[k] += deltaE;
+            m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
+        }
+        m_xlast = xnow;
+        m_tlast = tnow;
+    }
+}
+
+void BinarySolutionTabulatedThermo::initThermoXML(XML_Node& phaseNode, const std::string& id_)
+{
+    vector_fp x, h, s;
+    std::vector<std::pair<double,double>> x_h_temp, x_s_temp;
+
+    if (id_.size() > 0) {
+        if (phaseNode.id() != id_) {
+            throw CanteraError("BinarySolutionTabulatedThermo::initThermoXML",
+                    "phasenode and Id are incompatible");
+        }
+    }
+    if (nSpecies()!=2) {
+        throw CanteraError("BinarySolutionTabulatedThermo::initThermoXML",
+                            "No. of species should be equal to 2!");
+    }
+    if (phaseNode.hasChild("thermo")) {
+        XML_Node& thermoNode = phaseNode.child("thermo");
+        std::string mString = thermoNode["model"];
+        if (!caseInsensitiveEquals(mString, "binarysolutiontabulatedthermo")) {
+            throw CanteraError("BinarySolutionTabulatedThermo::initThermoXML",
+                    "Unknown thermo model: " + mString);
+        }
+        if (thermoNode.hasChild("tabulatedSpecies")) {
+            XML_Node& speciesNode = thermoNode.child("tabulatedSpecies");
+            std::string tabulated_species_name = speciesNode["name"];
+            m_kk_tab = speciesIndex(tabulated_species_name);
+            if (m_kk_tab == npos) {
+                throw CanteraError("BinarySolutionTabulatedThermo::initThermoXML",
+                        "Species " + tabulated_species_name + " not found.");
+            }
+            m_xlast = moleFraction(m_kk_tab);
+        }
+        if (thermoNode.hasChild("tabulatedThermo")) {
+            XML_Node& dataNode = thermoNode.child("tabulatedThermo");
+            getFloatArray(dataNode, x, false, "", "moleFraction");
+            getFloatArray(dataNode, h, false, "", "enthalpy");
+            getFloatArray(dataNode, s, false, "", "entropy");
+
+            // Check for data length consistency
+            if ((x.size() != h.size()) || (x.size() != s.size()) || (h.size() != s.size())) {
+                throw CanteraError("BinarySolutionTabulatedThermo::initThermoXML",
+                        "Species tabulated thermo data has different lengths.");
+            }
+            // Sort the x, h, s data in the order of increasing x
+            for(size_t i = 0; i < x.size(); i++){
+                x_h_temp.push_back(std::make_pair(x[i],h[i]));
+                x_s_temp.push_back(std::make_pair(x[i],s[i]));
+            }
+            std::sort(x_h_temp.begin(), x_h_temp.end());
+            std::sort(x_s_temp.begin(), x_s_temp.end());
+
+            // Store the sorted values in different arrays
+            m_molefrac_tab.resize(x_h_temp.size());
+            m_enthalpy_tab.resize(x_h_temp.size());
+            m_entropy_tab.resize(x_h_temp.size());
+            for (size_t i = 0; i < x_h_temp.size(); i++) {
+                m_molefrac_tab[i] = x_h_temp[i].first;
+                m_enthalpy_tab[i] = x_h_temp[i].second;
+                m_entropy_tab[i] = x_s_temp[i].second;
+            }
+        } else {
+            throw CanteraError("BinarySolutionTabulatedThermo::initThermoXML",
+                    "Unspecified tabulated species or thermo");
+        }
+    } else {
+        throw CanteraError("BinarySolutionTabulatedThermo::initThermoXML",
+                "Unspecified thermo model");
+    }
+
+    /*
+     * Form of the standard concentrations. Must have one of:
+     *
+     *     <standardConc model="unity" />
+     *     <standardConc model="molar_volume" />
+     *     <standardConc model="solvent_volume" />
+     */
+    if (phaseNode.hasChild("standardConc")) {
+        XML_Node& scNode = phaseNode.child("standardConc");
+        std::string formString = scNode.attrib("model");
+        if (caseInsensitiveEquals(formString, "unity")) {
+            m_formGC = 0;
+        } else if (caseInsensitiveEquals(formString, "molar_volume")) {
+            m_formGC = 1;
+        } else if (caseInsensitiveEquals(formString, "solvent_volume")) {
+            m_formGC = 2;
+        } else {
+            throw CanteraError("BinarySolutionTabulatedThermo::initThermoXML",
+                    "Unknown standardConc model: " + formString);
+        }
+    } else {
+        throw CanteraError("BinarySolutionTabulatedThermo::initThermoXML",
+                "Unspecified standardConc model");
+    }
+
+    // Call the base initThermo, which handles setting the initial state
+    ThermoPhase::initThermoXML(phaseNode, id_);
+}
+
+double* BinarySolutionTabulatedThermo::interpolate(double x) const
+{
+    static double c[2];
+    // Check if x is out of bound
+    if (x > m_molefrac_tab.back()) {
+        c[0] = m_enthalpy_tab.back();
+        c[1] = m_entropy_tab.back();
+        return c;
+    }
+    if (x < m_molefrac_tab[0]) {
+        c[0] = m_enthalpy_tab[0];
+        c[1] = m_entropy_tab[0];
+        return c;
+    }
+    size_t i = std::distance(m_molefrac_tab.begin(), std::lower_bound(m_molefrac_tab.begin(), m_molefrac_tab.end(), x));
+    c[0] = m_enthalpy_tab[i-1] + (m_enthalpy_tab[i] - m_enthalpy_tab[i-1]) * (x - m_molefrac_tab[i-1])/(m_molefrac_tab[i]- m_molefrac_tab[i-1]);
+    c[1] = m_entropy_tab[i-1] + (m_entropy_tab[i] - m_entropy_tab[i-1]) * (x - m_molefrac_tab[i-1])/(m_molefrac_tab[i]- m_molefrac_tab[i-1]);
+    return c;
+}
+
+}

src/thermo/ThermoFactory.cpp

@@ -36,6 +36,7 @@
 #include "cantera/thermo/IdealMolalSoln.h"
 #include "cantera/thermo/IdealSolnGasVPSS.h"
 #include "cantera/thermo/WaterSSTP.h"
+#include "cantera/thermo/BinarySolutionTabulatedThermo.h"
 #include "cantera/base/stringUtils.h"
 
 using namespace std;
@@ -71,6 +72,7 @@ ThermoFactory::ThermoFactory()
     m_synonyms["RedlichKwongMFTP"] = "RedlichKwong";
     reg("MaskellSolidSolnPhase", []() { return new MaskellSolidSolnPhase(); });
     reg("PureLiquidWater", []() { return new WaterSSTP(); });
+    reg("BinarySolutionTabulatedThermo", []() { return new BinarySolutionTabulatedThermo(); });
 }
 
 ThermoPhase* ThermoFactory::newThermoPhase(const std::string& model)