research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Ignore negative values when setting composition (mass fractions etc.)

Browse source 
This should resolve a number of problems when integrating reactor networks,
as well as some other unexpected behavior.
This commit is contained in:
Ray Speth 2012-03-25 23:25:41 +00:00
parent c56bfab0ed
commit bf1f14c178
1 changed files with 28 additions and 15 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

43
src/thermo/Phase.cpp
View file

@ -316,12 +316,22 @@ void Phase::restoreState(size_t lenstate, const doublereal* state)
void Phase::setMoleFractions(const doublereal* const x)
{
doublereal sum = dot(x, x + m_kk, m_molwts.begin());
doublereal rsum = 1.0/sum;
transform(x, x + m_kk, m_ym.begin(), timesConstant<double>(rsum));
// Use m_y as a temporary work vector for the non-negative mole fractions
doublereal norm = 0.0;
doublereal sum = 0;
for (size_t k = 0; k < m_kk; k++) {
double xk = std::max(x[k], 0.0); // Ignore negative mole fractions
m_y[k] = xk;
norm += xk;
sum += m_molwts[k] * xk;
}
transform(m_y.begin(), m_y.end(), m_ym.begin(),
timesConstant<double>(1.0/sum));
// Now set m_y to the mass fractions
transform(m_ym.begin(), m_ym.begin() + m_kk, m_molwts.begin(),
m_y.begin(), multiplies<double>());
doublereal norm = accumulate(x, x + m_kk, 0.0);
m_mmw = sum/norm;
// Call a routine to determine whether state has changed.
@ -367,15 +377,15 @@ void Phase::setMoleFractionsByName(const std::string& x)
void Phase::setMassFractions(const doublereal* const y)
{
doublereal norm = 0.0, sum = 0.0;
norm = accumulate(y, y + m_kk, 0.0);
copy(y, y + m_kk, m_y.begin());
scale(y, y + m_kk, m_y.begin(), 1.0/norm);
for (size_t k = 0; k < m_kk; k++) {
m_y[k] = std::max(y[k], 0.0); // Ignore negative mass fractions
}
doublereal norm = accumulate(m_y.begin(), m_y.end(), 0.0);
scale(m_y.begin(), m_y.end(), m_y.begin(), 1.0/norm);
transform(m_y.begin(), m_y.begin() + m_kk, m_rmolwts.begin(),
transform(m_y.begin(), m_y.end(), m_rmolwts.begin(),
m_ym.begin(), multiplies<double>());
sum = accumulate(m_ym.begin(), m_ym.begin() + m_kk, 0.0);
m_mmw = 1.0/sum;
m_mmw = 1.0 / accumulate(m_ym.begin(), m_ym.end(), 0.0);
// Call a routine to determine whether state has changed.
stateMFChangeCalc();
@ -585,17 +595,20 @@ void Phase::getConcentrations(doublereal* const c) const
void Phase::setConcentrations(const doublereal* const conc)
{
// Use m_y as temporary storage for non-negative concentrations
doublereal sum = 0.0, norm = 0.0;
for (size_t k = 0; k != m_kk; ++k) {
sum += conc[k]*m_molwts[k];
norm += conc[k];
double ck = std::max(conc[k], 0.0); // Ignore negative concentrations
m_y[k] = ck;
sum += ck * m_molwts[k];
norm += ck;
}
m_mmw = sum/norm;
setDensity(sum);
doublereal rsum = 1.0/sum;
for (size_t k = 0; k != m_kk; ++k) {
m_ym[k] = conc[k] * rsum;
m_y[k] = m_ym[k] * m_molwts[k];
m_ym[k] = m_y[k] * rsum;
m_y[k] = m_ym[k] * m_molwts[k]; // m_y is now the mass fraction
}
// Call a routine to determine whether state has changed.