From 71ef5ebfcced8bcbca8b10adb40b0fc01bcfa126 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Mon, 3 Jul 2006 19:57:29 +0000 Subject: [PATCH] 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. --- Cantera/src/RateCoeffMgr.h | 5 +- Cantera/src/RxnRates.h | 461 +++++++++++++++++++++++-------------- 2 files changed, 291 insertions(+), 175 deletions(-) diff --git a/Cantera/src/RateCoeffMgr.h b/Cantera/src/RateCoeffMgr.h index c14b6d409..ae8f4e1c3 100755 --- a/Cantera/src/RateCoeffMgr.h +++ b/Cantera/src/RateCoeffMgr.h @@ -22,7 +22,7 @@ namespace Cantera { /** * This rate coefficient manager supports one parameterization of - * any type. + * the rate constant of any type. */ template 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); } } diff --git a/Cantera/src/RxnRates.h b/Cantera/src/RxnRates.h index 1ae3cd393..b34040c5e 100755 --- a/Cantera/src/RxnRates.h +++ b/Cantera/src/RxnRates.h @@ -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 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 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