Fix for negative pre-exponential factors producing NaN rates
of progress. This bug was reported in the Cantera newsgroup.
The fix consists of replacing
Rc = exp ( mlogA + blogT - Ea / RT)
with
Rc = A * exp (b logT - Ea/RT)
Therefore, the log of A never has to be taken.
Note, this also allows for a zero rate of progress for a reaction.
Consequences:
The permissible range of the rate constant values may be altered.
Numerical roundoff differences will occur, since the order of
operations in key rate constant evalulations has been changed.
Speed of the calculation may be affected. I have no idea which way.
The fix should be considered as provisional. There are different pathways
for fixing this. This is merely the simplest.
This commit is contained in:
Harry Moffat
parent
228ebc78b1
commit
71ef5ebfcc
2 changed files with 291 additions and 175 deletions
|
|
@ -22,7 +22,7 @@ namespace Cantera {
|
|||
|
||||
/**
|
||||
* This rate coefficient manager supports one parameterization of
|
||||
* any type.
|
||||
* the rate constant of any type.
|
||||
*/
|
||||
template<class R>
|
||||
class Rate1 {
|
||||
|
|
@ -104,7 +104,8 @@ namespace Cantera {
|
|||
doublereal recipT = 1.0/T;
|
||||
int i = 0;
|
||||
for (; b != e; ++b, ++i) {
|
||||
values[m_rxn[i]] = exp(b->update(logT, recipT));
|
||||
// values[m_rxn[i]] = exp(b->update(logT, recipT));
|
||||
values[m_rxn[i]] = b->updateRC(logT, recipT);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -20,200 +20,315 @@
|
|||
namespace Cantera {
|
||||
|
||||
|
||||
class Arrhenius {
|
||||
class Arrhenius {
|
||||
|
||||
public:
|
||||
static int type(){ return ARRHENIUS; }
|
||||
Arrhenius() : m_b (0.0), m_E (0.0) {}
|
||||
Arrhenius( int csize, const doublereal* c )
|
||||
: m_b (c[1]), m_E (c[2]) { m_logA = log(c[0]);}
|
||||
Arrhenius( doublereal A, doublereal b, doublereal E)
|
||||
: m_b (b), m_E (E) { m_logA = log(A);}
|
||||
public:
|
||||
static int type(){ return ARRHENIUS; }
|
||||
Arrhenius() :
|
||||
m_logA(-1.0E300),
|
||||
m_b (0.0),
|
||||
m_E (0.0),
|
||||
m_A(0.0)
|
||||
{
|
||||
}
|
||||
Arrhenius(int csize, const doublereal* c) :
|
||||
m_b (c[1]),
|
||||
m_E (c[2]),
|
||||
m_A (c[0])
|
||||
{
|
||||
if (m_A <= 0.0) {
|
||||
m_logA = -1.0E300;
|
||||
} else {
|
||||
m_logA = log(m_A);
|
||||
}
|
||||
}
|
||||
|
||||
void update_C(const doublereal* c) {}
|
||||
Arrhenius(doublereal A, doublereal b, doublereal E) :
|
||||
m_b (b),
|
||||
m_E (E),
|
||||
m_A (A)
|
||||
{
|
||||
if (m_A <= 0.0) {
|
||||
m_logA = -1.0E300;
|
||||
} else {
|
||||
m_logA = log(m_A);
|
||||
}
|
||||
}
|
||||
|
||||
void update_C(const doublereal* c) {}
|
||||
|
||||
doublereal update(doublereal logT, doublereal recipT) const {
|
||||
return m_logA + m_b*logT - m_E*recipT;
|
||||
}
|
||||
|
||||
/// no longer used
|
||||
//doublereal update_dT(doublereal logT, doublereal recipT) const {
|
||||
// return recipT*(m_b + m_E*recipT);
|
||||
//}
|
||||
|
||||
void writeUpdateRHS(ostream& s) const {
|
||||
s << " exp(" << m_logA;
|
||||
if (m_b != 0.0) s << " + " << m_b << " * tlog";
|
||||
if (m_E != 0.0) s << " - " << m_E << " * rt";
|
||||
s << ");" << endl;
|
||||
}
|
||||
|
||||
doublereal activationEnergy_R() const {
|
||||
return m_E;
|
||||
}
|
||||
|
||||
static bool alwaysComputeRate() { return false;}
|
||||
|
||||
protected:
|
||||
doublereal m_logA, m_b, m_E;
|
||||
};
|
||||
|
||||
|
||||
class ArrheniusSum {
|
||||
|
||||
public:
|
||||
static int type(){ return ARRHENIUS_SUM; }
|
||||
ArrheniusSum() : m_nterms(0) {}
|
||||
ArrheniusSum( int csize, const doublereal* c ) {
|
||||
m_nterms = 0;
|
||||
addArrheniusTerm(c[0], c[1], c[2]);
|
||||
}
|
||||
void addArrheniusTerm(doublereal A, doublereal b, doublereal E) {
|
||||
m_terms.push_back(Arrhenius(A, b, E));
|
||||
m_nterms++;
|
||||
}
|
||||
|
||||
void update_C(const doublereal* c) {}
|
||||
|
||||
doublereal update(doublereal logT, doublereal recipT) const {
|
||||
int n;
|
||||
doublereal f, fexp = 0.0;
|
||||
for (n = 0; n < m_nterms; n++) {
|
||||
f = m_terms[n].update(logT, recipT);
|
||||
fexp += exp(f);
|
||||
}
|
||||
return log(fexp);
|
||||
}
|
||||
|
||||
// doublereal update_dT(doublereal logT, doublereal recipT) const {
|
||||
// throw CanteraError("ArrheniusSum::update_dT","not implemented.");
|
||||
//}
|
||||
|
||||
void writeUpdateRHS(ostream& s) const {
|
||||
;
|
||||
}
|
||||
|
||||
//doublereal activationEnergy_R() const {
|
||||
// return m_E;
|
||||
//}
|
||||
|
||||
static bool alwaysComputeRate() { return true;}
|
||||
|
||||
protected:
|
||||
vector<Arrhenius> m_terms;
|
||||
int m_nterms;
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
* Update the value of the logarithm of the rate constant.
|
||||
*
|
||||
* Note, this function should never be called for negative A values.
|
||||
* If it does then it will produce a negative overflow result, and
|
||||
* a zero net forwards reaction rate, instead of a negative reaction
|
||||
* rate constant that is the expected result.
|
||||
*/
|
||||
doublereal update(doublereal logT, doublereal recipT) const {
|
||||
return m_logA + m_b*logT - m_E*recipT;
|
||||
}
|
||||
|
||||
/**
|
||||
* An Arrhenius rate with coverage-dependent terms.
|
||||
* Update the value the rate constant.
|
||||
*
|
||||
* This function returns the actual value of the rate constant.
|
||||
* It can be safely called for negative values of the pre-exponential
|
||||
* factor.
|
||||
*/
|
||||
class SurfaceArrhenius {
|
||||
doublereal updateRC(doublereal logT, doublereal recipT) const {
|
||||
return m_A * exp(m_b*logT - m_E*recipT);
|
||||
}
|
||||
|
||||
public:
|
||||
static int type(){ return ARRHENIUS; }
|
||||
SurfaceArrhenius() : m_b (0.0), m_E (0.0),
|
||||
m_acov(0.0), m_ecov(0.0),
|
||||
m_mcov(0.0), m_ncov(0), m_nmcov(0)
|
||||
{}
|
||||
SurfaceArrhenius( int csize, const doublereal* c )
|
||||
: m_b (c[1]), m_E (c[2]),
|
||||
m_acov(0.0), m_ecov(0.0), m_mcov(0.0), m_ncov(0), m_nmcov(0)
|
||||
{ m_logA = log(c[0]);
|
||||
if (csize >= 7) {
|
||||
for (int n = 3; n < csize-3; n += 4) {
|
||||
addCoverageDependence(int(c[n]),
|
||||
c[n+1], c[n+2], c[n+3]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void addCoverageDependence(int k, doublereal a,
|
||||
doublereal m, doublereal e) {
|
||||
m_ncov++;
|
||||
m_sp.push_back(k);
|
||||
m_ac.push_back(a);
|
||||
m_ec.push_back(e);
|
||||
if (m != 0.0) {
|
||||
m_msp.push_back(k);
|
||||
m_mc.push_back(m);
|
||||
m_nmcov++;
|
||||
}
|
||||
}
|
||||
|
||||
void update_C(const doublereal* theta) {
|
||||
m_acov = 0.0;
|
||||
m_ecov = 0.0;
|
||||
m_mcov = 0.0;
|
||||
int n, k;
|
||||
doublereal th;
|
||||
for (n = 0; n < m_ncov; n++) {
|
||||
k = m_sp[n];
|
||||
m_acov += m_ac[n] * theta[k];
|
||||
m_ecov += m_ec[n] * theta[k];
|
||||
}
|
||||
for (n = 0; n < m_nmcov; n++) {
|
||||
k = m_msp[n];
|
||||
// changed n to k, dgg 1/22/04
|
||||
th = fmaxx(theta[k], Tiny);
|
||||
// th = fmaxx(theta[n], Tiny);
|
||||
m_mcov += m_mc[n]*log(th);
|
||||
}
|
||||
}
|
||||
/// no longer used
|
||||
//doublereal update_dT(doublereal logT, doublereal recipT) const {
|
||||
// return recipT*(m_b + m_E*recipT);
|
||||
//}
|
||||
|
||||
doublereal update(doublereal logT, doublereal recipT) const {
|
||||
return m_logA + m_acov + m_b*logT
|
||||
- (m_E + m_ecov)*recipT + m_mcov;
|
||||
}
|
||||
void writeUpdateRHS(ostream& s) const {
|
||||
s << " exp(" << m_logA;
|
||||
if (m_b != 0.0) s << " + " << m_b << " * tlog";
|
||||
if (m_E != 0.0) s << " - " << m_E << " * rt";
|
||||
s << ");" << endl;
|
||||
}
|
||||
|
||||
doublereal activationEnergy_R() const {
|
||||
return m_E + m_ecov;
|
||||
}
|
||||
doublereal activationEnergy_R() const {
|
||||
return m_E;
|
||||
}
|
||||
|
||||
static bool alwaysComputeRate() { return true;}
|
||||
static bool alwaysComputeRate() { return false;}
|
||||
|
||||
protected:
|
||||
doublereal m_logA, m_b, m_E;
|
||||
doublereal m_acov, m_ecov, m_mcov;
|
||||
vector_int m_sp, m_msp;
|
||||
vector_fp m_ac, m_ec, m_mc;
|
||||
int m_ncov, m_nmcov;
|
||||
};
|
||||
protected:
|
||||
doublereal m_logA, m_b, m_E, m_A;
|
||||
};
|
||||
|
||||
|
||||
class ArrheniusSum {
|
||||
|
||||
public:
|
||||
static int type(){ return ARRHENIUS_SUM; }
|
||||
ArrheniusSum() : m_nterms(0) {}
|
||||
ArrheniusSum( int csize, const doublereal* c ) {
|
||||
m_nterms = 0;
|
||||
addArrheniusTerm(c[0], c[1], c[2]);
|
||||
}
|
||||
void addArrheniusTerm(doublereal A, doublereal b, doublereal E) {
|
||||
m_terms.push_back(Arrhenius(A, b, E));
|
||||
m_nterms++;
|
||||
}
|
||||
|
||||
void update_C(const doublereal* c) {}
|
||||
|
||||
/**
|
||||
* Update the value of the logarithm of the rate constant.
|
||||
*
|
||||
*/
|
||||
doublereal update(doublereal logT, doublereal recipT) const {
|
||||
int n;
|
||||
doublereal f, fsum = 0.0;
|
||||
for (n = 0; n < m_nterms; n++) {
|
||||
f = m_terms[n].updateRC(logT, recipT);
|
||||
fsum += f;
|
||||
}
|
||||
//if (fsum <= 0.0) {
|
||||
// throw CanteraError("ArrheniusSum::update",
|
||||
// "Error: negative total reaction rate");
|
||||
//}
|
||||
return log(fsum);
|
||||
}
|
||||
|
||||
/**
|
||||
* Update the value the rate constant.
|
||||
*
|
||||
* This function returns the actual value of the rate constant.
|
||||
* It can be safely called for negative values of the pre-exponential
|
||||
* factor.
|
||||
*/
|
||||
doublereal updateRC(doublereal logT, doublereal recipT) const {
|
||||
int n;
|
||||
doublereal f, fsum = 0.0;
|
||||
for (n = 0; n < m_nterms; n++) {
|
||||
f = m_terms[n].updateRC(logT, recipT);
|
||||
fsum += f;
|
||||
}
|
||||
//if (fsum <= 0.0) {
|
||||
// throw CanteraError("ArrheniusSum::update",
|
||||
// "Error: negative total reaction rate");
|
||||
//}
|
||||
return fsum;
|
||||
}
|
||||
|
||||
|
||||
// doublereal update_dT(doublereal logT, doublereal recipT) const {
|
||||
// throw CanteraError("ArrheniusSum::update_dT","not implemented.");
|
||||
//}
|
||||
|
||||
void writeUpdateRHS(ostream& s) const {
|
||||
;
|
||||
}
|
||||
|
||||
//doublereal activationEnergy_R() const {
|
||||
// return m_E;
|
||||
//}
|
||||
|
||||
static bool alwaysComputeRate() { return true;}
|
||||
|
||||
protected:
|
||||
vector<Arrhenius> m_terms;
|
||||
int m_nterms;
|
||||
};
|
||||
|
||||
/**
|
||||
* An Arrhenius rate with coverage-dependent terms.
|
||||
*/
|
||||
class SurfaceArrhenius {
|
||||
|
||||
public:
|
||||
static int type(){ return ARRHENIUS; }
|
||||
SurfaceArrhenius() :
|
||||
m_logA(-1.0E300),
|
||||
m_b (0.0),
|
||||
m_E (0.0),
|
||||
m_A(0.0),
|
||||
m_acov(0.0),
|
||||
m_ecov(0.0),
|
||||
m_mcov(0.0),
|
||||
m_ncov(0),
|
||||
m_nmcov(0)
|
||||
{
|
||||
}
|
||||
|
||||
SurfaceArrhenius( int csize, const doublereal* c ) :
|
||||
m_b (c[1]),
|
||||
m_E (c[2]),
|
||||
m_A (c[0]),
|
||||
m_acov(0.0),
|
||||
m_ecov(0.0),
|
||||
m_mcov(0.0),
|
||||
m_ncov(0),
|
||||
m_nmcov(0)
|
||||
{
|
||||
if (m_A <= 0.0) {
|
||||
m_logA = -1.0E300;
|
||||
} else {
|
||||
m_logA = log(c[0]);
|
||||
}
|
||||
if (csize >= 7) {
|
||||
for (int n = 3; n < csize-3; n += 4) {
|
||||
addCoverageDependence(int(c[n]),
|
||||
c[n+1], c[n+2], c[n+3]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void addCoverageDependence(int k, doublereal a,
|
||||
doublereal m, doublereal e) {
|
||||
m_ncov++;
|
||||
m_sp.push_back(k);
|
||||
m_ac.push_back(a);
|
||||
m_ec.push_back(e);
|
||||
if (m != 0.0) {
|
||||
m_msp.push_back(k);
|
||||
m_mc.push_back(m);
|
||||
m_nmcov++;
|
||||
}
|
||||
}
|
||||
|
||||
void update_C(const doublereal* theta) {
|
||||
m_acov = 0.0;
|
||||
m_ecov = 0.0;
|
||||
m_mcov = 0.0;
|
||||
int n, k;
|
||||
doublereal th;
|
||||
for (n = 0; n < m_ncov; n++) {
|
||||
k = m_sp[n];
|
||||
m_acov += m_ac[n] * theta[k];
|
||||
m_ecov += m_ec[n] * theta[k];
|
||||
}
|
||||
for (n = 0; n < m_nmcov; n++) {
|
||||
k = m_msp[n];
|
||||
// changed n to k, dgg 1/22/04
|
||||
th = fmaxx(theta[k], Tiny);
|
||||
// th = fmaxx(theta[n], Tiny);
|
||||
m_mcov += m_mc[n]*log(th);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Update the value of the logarithm of the rate constant.
|
||||
*
|
||||
* This calculation is not safe for negative values of
|
||||
* the preexponential.
|
||||
*/
|
||||
doublereal update(doublereal logT, doublereal recipT) const {
|
||||
return m_logA + m_acov + m_b*logT
|
||||
- (m_E + m_ecov)*recipT + m_mcov;
|
||||
}
|
||||
|
||||
/**
|
||||
* Update the value the rate constant.
|
||||
*
|
||||
* This function returns the actual value of the rate constant.
|
||||
* It can be safely called for negative values of the pre-exponential
|
||||
* factor.
|
||||
*/
|
||||
doublereal updateRC(doublereal logT, doublereal recipT) const {
|
||||
return m_A * exp(m_acov + m_b*logT - (m_E + m_ecov)*recipT + m_mcov);
|
||||
}
|
||||
|
||||
doublereal activationEnergy_R() const {
|
||||
return m_E + m_ecov;
|
||||
}
|
||||
|
||||
static bool alwaysComputeRate() { return true;}
|
||||
|
||||
protected:
|
||||
doublereal m_logA, m_b, m_E, m_A;
|
||||
doublereal m_acov, m_ecov, m_mcov;
|
||||
vector_int m_sp, m_msp;
|
||||
vector_fp m_ac, m_ec, m_mc;
|
||||
int m_ncov, m_nmcov;
|
||||
};
|
||||
|
||||
|
||||
#ifdef INCL_TST
|
||||
|
||||
class TST {
|
||||
class TST {
|
||||
|
||||
public:
|
||||
static int type(){ return TSTRATE; }
|
||||
TST() {}
|
||||
TST( const vector_fp& c ) {
|
||||
m_b.resize(10);
|
||||
copy(c.begin(), c.begin() + 10, m_b.begin());
|
||||
m_k = int(c[10]);
|
||||
}
|
||||
public:
|
||||
static int type(){ return TSTRATE; }
|
||||
TST() {}
|
||||
TST( const vector_fp& c ) {
|
||||
m_b.resize(10);
|
||||
copy(c.begin(), c.begin() + 10, m_b.begin());
|
||||
m_k = int(c[10]);
|
||||
}
|
||||
|
||||
void update_C(const vector_fp& c) {
|
||||
doublereal ck = c[m_k];
|
||||
delta_s0 = m_b[0] + m_b[1]*ck + m_b[2]*ck*ck;
|
||||
delta_e0 = m_b[5] + m_b[6]*ck + m_b[7]*ck*ck;
|
||||
}
|
||||
void update_C(const vector_fp& c) {
|
||||
doublereal ck = c[m_k];
|
||||
delta_s0 = m_b[0] + m_b[1]*ck + m_b[2]*ck*ck;
|
||||
delta_e0 = m_b[5] + m_b[6]*ck + m_b[7]*ck*ck;
|
||||
}
|
||||
|
||||
doublereal update(doublereal logT, doublereal recipT) const {
|
||||
doublereal delta_s = delta_s0*(1.0 + m_b[3]*logT + m_b[4]*recipT);
|
||||
doublereal delta_E = delta_e0*(1.0 + m_b[8]*logT + m_b[9]*recipT);
|
||||
return logBoltz_Planck + logT + delta_s - delta_E*recipT;
|
||||
}
|
||||
doublereal update(doublereal logT, doublereal recipT) const {
|
||||
doublereal delta_s = delta_s0*(1.0 + m_b[3]*logT + m_b[4]*recipT);
|
||||
doublereal delta_E = delta_e0*(1.0 + m_b[8]*logT + m_b[9]*recipT);
|
||||
return logBoltz_Planck + logT + delta_s - delta_E*recipT;
|
||||
}
|
||||
|
||||
doublereal updateRC(doublereal logT, doublereal recipT) const {
|
||||
double lres = update(logT, recipT);
|
||||
return exp(lres);
|
||||
}
|
||||
|
||||
void writeUpdateRHS(ostream& s) const {}
|
||||
void writeUpdateRHS(ostream& s) const {}
|
||||
|
||||
protected:
|
||||
doublereal delta_s0, delta_e0;
|
||||
int m_k;
|
||||
vector_fp m_b;
|
||||
};
|
||||
protected:
|
||||
doublereal delta_s0, delta_e0;
|
||||
int m_k;
|
||||
vector_fp m_b;
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue