From 224ef720e6f70f337bfeacc54ba565288841cae7 Mon Sep 17 00:00:00 2001 From: Manik Mayur Date: Fri, 14 Sep 2018 16:36:56 +0200 Subject: [PATCH] Add BinarySolutionTabulatedThermo class --- .../thermo/BinarySolutionTabulatedThermo.h | 110 +++++++++ interfaces/cython/cantera/ctml_writer.py | 200 +++++++++++++++- src/thermo/BinarySolutionTabulatedThermo.cpp | 226 ++++++++++++++++++ src/thermo/ThermoFactory.cpp | 2 + 4 files changed, 534 insertions(+), 4 deletions(-) create mode 100644 include/cantera/thermo/BinarySolutionTabulatedThermo.h create mode 100644 src/thermo/BinarySolutionTabulatedThermo.cpp diff --git a/include/cantera/thermo/BinarySolutionTabulatedThermo.h b/include/cantera/thermo/BinarySolutionTabulatedThermo.h new file mode 100644 index 000000000..b9393db6c --- /dev/null +++ b/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 diff --git a/interfaces/cython/cantera/ctml_writer.py b/interfaces/cython/cantera/ctml_writer.py index 62b65defe..372cc7eb7 100644 --- a/interfaces/cython/cantera/ctml_writer.py +++ b/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, diff --git a/src/thermo/BinarySolutionTabulatedThermo.cpp b/src/thermo/BinarySolutionTabulatedThermo.cpp new file mode 100644 index 000000000..474fbe17c --- /dev/null +++ b/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 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> 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: + * + * + * + * + */ + 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; +} + +} diff --git a/src/thermo/ThermoFactory.cpp b/src/thermo/ThermoFactory.cpp index 9520647d6..14e1a20c3 100644 --- a/src/thermo/ThermoFactory.cpp +++ b/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)