diff --git a/doc/sphinx/cti/phases.rst b/doc/sphinx/cti/phases.rst index fd8d4d7c4..d08dfc36f 100644 --- a/doc/sphinx/cti/phases.rst +++ b/doc/sphinx/cti/phases.rst @@ -269,6 +269,9 @@ Option String Meaning ``'skip_undeclared_third_bodies'`` When importing reactions with third body efficiencies, ignore any efficiencies for undeclared species, rather than flagging them as an error. +``'allow_discontinuous_thermo'`` Disable the automatic adjustment of NASA polynomials to + eliminate discontinuities in enthalpy and entropy at the + midpoint temperature. ================================== ================ Using the ``options`` field, it is possible to extract a sub-mechanism from a large diff --git a/include/cantera/thermo/SpeciesThermo.h b/include/cantera/thermo/SpeciesThermo.h index 10f8f7fb7..e64f4da0b 100644 --- a/include/cantera/thermo/SpeciesThermo.h +++ b/include/cantera/thermo/SpeciesThermo.h @@ -156,7 +156,7 @@ class SpeciesThermo public: //! Constructor - SpeciesThermo() {} + SpeciesThermo() : m_allow_discontinuities(false) {} //! Destructor virtual ~SpeciesThermo() {} @@ -332,6 +332,8 @@ public: */ virtual void modifyOneHf298(const int k, const doublereal Hf298New) = 0; #endif + + bool m_allow_discontinuities; }; //@} } diff --git a/interfaces/python/ctml_writer.py b/interfaces/python/ctml_writer.py index e7f168c42..df569cd2f 100644 --- a/interfaces/python/ctml_writer.py +++ b/interfaces/python/ctml_writer.py @@ -1810,6 +1810,10 @@ class phase(object): if self._comment: ph.addChild('note',self._comment) + thermo = ph.addChild('thermo') + if 'allow_discontinuous_thermo' in self._options: + thermo['allow_discontinuities'] = 'true' + return ph @@ -1854,8 +1858,7 @@ class ideal_gas(phase): def build(self, p): ph = phase.build(self, p) - e = ph.addChild("thermo") - e['model'] = 'IdealGas' + ph.child('thermo')['model'] = 'IdealGas' k = ph.addChild("kinetics") k['model'] = self._kin t = ph.addChild('transport') @@ -1900,7 +1903,7 @@ class stoichiometric_solid(phase): def build(self, p): ph = phase.build(self, p) - e = ph.addChild("thermo") + e = ph.child('thermo') e['model'] = 'StoichSubstance' addFloat(e, 'density', self._dens, defunits = _umass+'/'+_ulen+'3') if self._tr: @@ -1955,7 +1958,7 @@ class metal(phase): def build(self, p): ph = phase.build(self, p) - e = ph.addChild("thermo") + e = ph.child("thermo") e['model'] = 'Metal' addFloat(e, 'density', self._dens, defunits = _umass+'/'+_ulen+'3') if self._tr: @@ -1993,7 +1996,7 @@ class semiconductor(phase): def build(self, p): ph = phase.build(self, p) - e = ph.addChild("thermo") + e = ph.child("thermo") e['model'] = 'Semiconductor' addFloat(e, 'density', self._dens, defunits = _umass+'/'+_ulen+'3') addFloat(e, 'effectiveMass_e', self._emass, defunits = _umass) @@ -2031,7 +2034,7 @@ class incompressible_solid(phase): def build(self, p): ph = phase.build(self, p) - e = ph.addChild("thermo") + e = ph.child("thermo") e['model'] = 'Incompressible' addFloat(e, 'density', self._dens, defunits = _umass+'/'+_ulen+'3') if self._tr: @@ -2071,7 +2074,7 @@ class lattice(phase): #if visible == 0: # return ph = phase.build(self, p) - e = ph.addChild('thermo') + e = ph.child('thermo') e['model'] = 'Lattice' addFloat(e, 'site_density', self._n, defunits = _umol+'/'+_ulen+'3') if self._vac: @@ -2128,7 +2131,7 @@ class lattice_solid(phase): def build(self, p): ph = phase.build(self, p) - e = ph.addChild("thermo") + e = ph.child("thermo") e['model'] = 'LatticeSolid' if self._lattices: @@ -2171,7 +2174,7 @@ class liquid_vapor(phase): def build(self, p): ph = phase.build(self, p) - e = ph.addChild("thermo") + e = ph.child("thermo") e['model'] = 'PureFluid' e['fluid_type'] = repr(self._subflag) k = ph.addChild("kinetics") @@ -2202,7 +2205,7 @@ class RedlichKwongMFTP(phase): def build(self, p): ph = phase.build(self,p) - e = ph.addChild("thermo") + e = ph.child("thermo") e['model'] = 'RedlichKwongMFTP' if self._activityCoefficients: a = e.addChild("activityCoefficients") @@ -2250,7 +2253,7 @@ class redlich_kwong(phase): def build(self, p): ph = phase.build(self, p) - e = ph.addChild("thermo") + e = ph.child("thermo") e['model'] = 'PureFluid' e['fluid_type'] = repr(self._subflag) addFloat(e, 'Tc', self._tc, defunits = "K") @@ -2300,7 +2303,7 @@ class ideal_interface(phase): def build(self, p): ph = phase.build(self, p) - e = ph.addChild("thermo") + e = ph.child("thermo") e['model'] = 'Surface' addFloat(e, 'site_density', self._sitedens, defunits = _umol+'/'+_ulen+'2') k = ph.addChild("kinetics") @@ -2340,7 +2343,7 @@ class edge(phase): def build(self, p): ph = phase.build(self, p) - e = ph.addChild("thermo") + e = ph.child("thermo") e['model'] = 'Edge' addFloat(e, 'site_density', self._sitedens, defunits = _umol+'/'+_ulen) k = ph.addChild("kinetics") @@ -2476,7 +2479,7 @@ class edge(phase): ## def build(self, p): ## ph = phase.build(self, p) -## e = ph.addChild("thermo") +## e = ph.child("thermo") ## sc = e.addChild("standardConc") ## sc['model'] = self._stdconc ## e['model'] = 'HMW' diff --git a/src/thermo/NasaThermo.cpp b/src/thermo/NasaThermo.cpp index ca5f22806..9916d2ad5 100644 --- a/src/thermo/NasaThermo.cpp +++ b/src/thermo/NasaThermo.cpp @@ -82,12 +82,14 @@ void NasaThermo::install(const std::string& name, size_t index, int type, vector_fp chigh(c+8, c+15); vector_fp clow(c+1, c+8); - checkContinuity(name, tmid, &clow[0], &chigh[0]); - //if (maxError > 1e-6) { - // fixDiscontinuities(tlow, tmid, thigh, &clow[0], &chigh[0]); - // AssertThrowMsg(checkContinuity(name, tmid, &clow[0], &chigh[0]) < 1e-12, - // "NasaThermo::install", "Polynomials still not continuous"); - //} + if (!m_allow_discontinuities) { + doublereal maxError = checkContinuity(name, tmid, &clow[0], &chigh[0]); + if (maxError > 1e-6) { + fixDiscontinuities(tlow, tmid, thigh, &clow[0], &chigh[0]); + AssertThrowMsg(checkContinuity(name, tmid, &clow[0], &chigh[0]) < 1e-12, + "NasaThermo::install", "Polynomials still not continuous"); + } + } m_high[igrp-1].push_back(NasaPoly1(index, tmid, thigh, ref_pressure, &chigh[0])); diff --git a/src/thermo/ThermoFactory.cpp b/src/thermo/ThermoFactory.cpp index a62d02a3e..36c7c0dc9 100644 --- a/src/thermo/ThermoFactory.cpp +++ b/src/thermo/ThermoFactory.cpp @@ -434,8 +434,8 @@ bool importPhase(XML_Node& phase, ThermoPhase* th, // Set equation of state parameters. The parameters are // specific to each subclass of ThermoPhase, so this is done // by method setParametersFromXML in each subclass. + const XML_Node& eos = phase.child("thermo"); if (phase.hasChild("thermo")) { - const XML_Node& eos = phase.child("thermo"); th->setParametersFromXML(eos); } else { throw CanteraError("importPhase", @@ -561,6 +561,10 @@ bool importPhase(XML_Node& phase, ThermoPhase* th, // used, and selects a class that can handle the // parameterizations found. spth = newSpeciesThermoMgr(spDataNodeList); + if (eos["allow_discontinuities"] == "true") { + std::cout << "ALLOWING DISCONTINUOUS THERMO!" << std::endl; + spth->m_allow_discontinuities = true; + } // install it in the phase object th->setSpeciesThermo(spth);