Skip to content

Commit

Permalink
Fix mole fractions near pure-species limit in BinarySolutionTabulated…
Browse files Browse the repository at this point in the history
…Thermo

Fixes behavior when the tabulated species mole fraction has become 1.0 to within
machine precision, while the other species mole fraction is not quite zero.
  • Loading branch information
speth committed Feb 26, 2023
1 parent 0cc44de commit ac60a59
Show file tree
Hide file tree
Showing 2 changed files with 10 additions and 13 deletions.
3 changes: 0 additions & 3 deletions include/cantera/thermo/BinarySolutionTabulatedThermo.h
Original file line number Diff line number Diff line change
Expand Up @@ -240,9 +240,6 @@ class BinarySolutionTabulatedThermo : public IdealSolidSolnPhase
//! Current tabulated species index
size_t m_kk_tab;

//! Current tabulated species mole fraction
mutable double m_xlast;

//! Tabulated contribution to h0[m_kk_tab] at the current composition
mutable double m_h0_tab;

Expand Down
20 changes: 10 additions & 10 deletions src/thermo/BinarySolutionTabulatedThermo.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -22,8 +22,6 @@ namespace Cantera
BinarySolutionTabulatedThermo::BinarySolutionTabulatedThermo(const std::string& inputFile,
const std::string& id_)
: m_kk_tab(npos)
, m_xlast(-1)

{
initThermoFile(inputFile, id_);
}
Expand All @@ -36,22 +34,24 @@ void BinarySolutionTabulatedThermo::compositionChanged()

void BinarySolutionTabulatedThermo::_updateThermo() const
{
double xnow = moleFraction(m_kk_tab);
bool x_changed = (m_xlast != xnow);
static const int cacheId = m_cache.getId();
CachedScalar cached = m_cache.getScalar(cacheId);
bool x_changed = !cached.validate(stateMFNumber());

if (x_changed) {
m_h0_tab = interpolate(xnow, m_enthalpy_tab);
m_s0_tab = interpolate(xnow, m_entropy_tab);
if (xnow == 0) {
double x_tab = moleFraction(m_kk_tab);
double x_other = moleFraction(1 - m_kk_tab);
m_h0_tab = interpolate(x_tab, m_enthalpy_tab);
m_s0_tab = interpolate(x_tab, m_entropy_tab);
if (x_tab == 0) {
m_s0_tab = -BigNumber;
} else if (xnow == 1) {
} else if (x_other == 0) {
m_s0_tab = BigNumber;
} else {
m_s0_tab += GasConstant*std::log(xnow/(1.0-xnow)) +
m_s0_tab += GasConstant*std::log(x_tab / x_other) +
GasConstant/Faraday*std::log(standardConcentration(1-m_kk_tab)
/standardConcentration(m_kk_tab));
}
m_xlast = xnow;
}

double tnow = temperature();
Expand Down

0 comments on commit ac60a59

Please sign in to comment.