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

[Kinetics] Add class BulkKinetics

Browse source 
This class serves as a common base for GasKinetics and AqueousKinetics,
eliminating most of the redundancy between the two classes.
This commit is contained in:
Ray Speth 2014-10-30 21:09:59 +00:00
parent 6f0fb5cf8d
commit 3b71d75ada
6 changed files with 262 additions and 662 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

130
include/cantera/kinetics/AqueousKinetics.h
View file

@ -1,6 +1,5 @@
/**
* @file AqueousKinetics.h
*
* @ingroup chemkinetics
*/
@ -9,44 +8,23 @@
#ifndef CT_AQUEOUSKINETICS_H
#define CT_AQUEOUSKINETICS_H
#include "Kinetics.h"
#include "ReactionStoichMgr.h"
#include "RateCoeffMgr.h"
#include "cantera/base/utilities.h"
#include "BulkKinetics.h"
namespace Cantera
{
// forward references
class ReactionData;
/**
* Kinetics manager for elementary aqueous-phase chemistry. This
* kinetics manager implements standard mass-action reaction rate
* expressions for liquids
*
*
* Concentration
* Kinetics manager for elementary aqueous-phase chemistry. This kinetics
* manager implements standard mass-action reaction rate expressions for liquids
*
* @ingroup kinetics
* @deprecated Not actually implemented
*/
class AqueousKinetics : public Kinetics
class AqueousKinetics : public BulkKinetics
{
public:
//! @name Constructors
//! @{
/// Constructor. Creates an empty reaction mechanism.
AqueousKinetics(thermo_t* thermo = 0);
AqueousKinetics(const AqueousKinetics& right);
AqueousKinetics& operator=(const AqueousKinetics& right);
//! Duplication routine for objects which inherit from Kinetics
/*!
* This virtual routine can be used to duplicate %Kinetics objects
@ -59,117 +37,25 @@ public:
* m_thermo vector within this object
*/
virtual Kinetics* duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const;
//@}
virtual int type() const {
return cAqueousKinetics;
}
virtual doublereal reactantStoichCoeff(size_t k, size_t i) const {
return getValue(m_rrxn[k], i, 0.0);
}
virtual doublereal productStoichCoeff(size_t k, size_t i) const {
return getValue(m_prxn[k], i, 0.0);
}
//! @name Reaction Rates Of Progress
//@{
virtual void getEquilibriumConstants(doublereal* kc);
virtual void getDeltaGibbs(doublereal* deltaG);
virtual void getDeltaEnthalpy(doublereal* deltaH);
virtual void getDeltaEntropy(doublereal* deltaS);
virtual void getDeltaSSGibbs(doublereal* deltaG);
virtual void getDeltaSSEnthalpy(doublereal* deltaH);
virtual void getDeltaSSEntropy(doublereal* deltaS);
//! @}
//! @name Reaction Mechanism Informational Query Routines
//! @{
virtual bool isReversible(size_t i) {
if (std::find(m_revindex.begin(), m_revindex.end(), i)
< m_revindex.end()) {
return true;
} else {
return false;
}
}
virtual void getFwdRateConstants(doublereal* kfwd);
virtual void getRevRateConstants(doublereal* krev,
bool doIrreversible = false);
//! @}
//! @name Reaction Mechanism Setup Routines
//! @{
virtual void init();
virtual void addReaction(ReactionData& r);
virtual void finalize();
virtual bool ready() const;
virtual void update_T();
virtual void update_C();
void updateROP();
/*!
* Update temperature-dependent portions of reaction rates and
* falloff functions.
*/
//! Update temperature-dependent portions of reaction rates
void _update_rates_T();
/*!
* Update properties that depend on concentrations. Currently only
* the enhanced collision partner concentrations are updated here.
*/
//! Update properties that depend on concentrations.
void _update_rates_C();
//@}
protected:
size_t m_nfall;
Rate1<Arrhenius> m_rates;
std::vector<size_t> m_irrev;
size_t m_nirrev;
size_t m_nrev;
/**
* Difference between the input global reactants order
* and the input global products order. Changed to a double
* to account for the fact that we can have real-valued
* stoichiometries.
*/
vector_fp m_dn;
std::vector<size_t> m_revindex;
vector_fp m_conc;
vector_fp m_grt;
//! @name Aqueous kinetics data
//!@{
bool m_ROP_ok;
doublereal m_temp;
//!@}
private:
void addElementaryReaction(ReactionData& r);
/**
* Update the equilibrium constants in molar units.
*/
//! Update the equilibrium constants in molar units.
void updateKc();
bool m_finalized;
virtual void addReaction(ReactionData& r);
};
}

63
include/cantera/kinetics/BulkKinetics.h Normal file
View file

@ -0,0 +1,63 @@
/**
* @file BulkKinetics.h
* @ingroup chemkinetics
*/
#ifndef CT_BULKKINETICS_H
#define CT_BULKKINETICS_H
#include "Kinetics.h"
#include "RateCoeffMgr.h"
namespace Cantera
{
//! Partial specialization of Kinetics for chemistry in a single bulk phase
class BulkKinetics : public Kinetics
{
public:
BulkKinetics(thermo_t* thermo = 0);
virtual Kinetics* duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const;
virtual bool isReversible(size_t i);
virtual void getDeltaGibbs(doublereal* deltaG);
virtual void getDeltaEnthalpy(doublereal* deltaH);
virtual void getDeltaEntropy(doublereal* deltaS);
virtual void getDeltaSSGibbs(doublereal* deltaG);
virtual void getDeltaSSEnthalpy(doublereal* deltaH);
virtual void getDeltaSSEntropy(doublereal* deltaS);
virtual void getRevRateConstants(doublereal* krev,
bool doIrreversible = false);
virtual void addReaction(ReactionData& r);
virtual void init();
virtual void finalize();
virtual bool ready() const;
protected:
virtual void addElementaryReaction(ReactionData& r);
Rate1<Arrhenius> m_rates;
std::vector<size_t> m_revindex; //!< Indices of reversible reactions
std::vector<size_t> m_irrev; //!< Indices of irreversible reactions
//! Difference between the global reactants order and the global products
//! order. Of type "double" to account for the fact that we can have real-
//! valued stoichiometries.
vector_fp m_dn;
vector_fp m_conc;
vector_fp m_grt;
bool m_ROP_ok;
doublereal m_temp;
bool m_finalized;
};
}
#endif

61
include/cantera/kinetics/GasKinetics.h
View file

@ -9,26 +9,20 @@
#ifndef CT_GASKINETICS_H
#define CT_GASKINETICS_H
#include "Kinetics.h"
#include "ReactionStoichMgr.h"
#include "BulkKinetics.h"
#include "ThirdBodyMgr.h"
#include "FalloffMgr.h"
#include "RateCoeffMgr.h"
namespace Cantera
{
// forward references
class Enhanced3BConc;
class ReactionData;
/**
* Kinetics manager for elementary gas-phase chemistry. This
* kinetics manager implements standard mass-action reaction rate
* expressions for low-density gases.
* @ingroup kinetics
*/
class GasKinetics : public Kinetics
class GasKinetics : public BulkKinetics
{
public:
//! @name Constructors and General Information
@ -40,12 +34,6 @@ public:
*/
GasKinetics(thermo_t* thermo = 0);
//! Copy Constructor
GasKinetics(const GasKinetics& right);
//! Assignment operator
GasKinetics& operator=(const GasKinetics& right);
virtual Kinetics* duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const;
virtual int type() const {
@ -57,32 +45,8 @@ public:
//! @{
virtual void getEquilibriumConstants(doublereal* kc);
virtual void getDeltaGibbs(doublereal* deltaG);
virtual void getDeltaEnthalpy(doublereal* deltaH);
virtual void getDeltaEntropy(doublereal* deltaS);
virtual void getDeltaSSGibbs(doublereal* deltaG);
virtual void getDeltaSSEnthalpy(doublereal* deltaH);
virtual void getDeltaSSEntropy(doublereal* deltaS);
//! @}
//! @name Reaction Mechanism Informational Query Routines
//! @{
virtual bool isReversible(size_t i) {
if (std::find(m_revindex.begin(), m_revindex.end(), i)
< m_revindex.end()) {
return true;
} else {
return false;
}
}
virtual void getFwdRateConstants(doublereal* kfwd);
virtual void getRevRateConstants(doublereal* krev,
bool doIrreversible = false);
//! @}
//! @name Reaction Mechanism Setup Routines
//! @{
@ -111,30 +75,15 @@ protected:
Rate1<Arrhenius> m_falloff_low_rates;
Rate1<Arrhenius> m_falloff_high_rates;
Rate1<Arrhenius> m_rates;
FalloffMgr m_falloffn;
ThirdBodyMgr<Enhanced3BConc> m_3b_concm;
ThirdBodyMgr<Enhanced3BConc> m_falloff_concm;
std::vector<size_t> m_irrev;
Rate1<Plog> m_plog_rates;
Rate1<ChebyshevRate> m_cheb_rates;
size_t m_nirrev;
size_t m_nrev;
/**
* Difference between the input global reactants order
* and the input global products order. Changed to a double
* to account for the fact that we can have real-valued
* stoichiometries.
*/
vector_fp m_dn;
std::vector<size_t> m_revindex;
//! @name Reaction rate data
//!@{
doublereal m_logp_ref;
@ -142,21 +91,15 @@ protected:
doublereal m_logStandConc;
vector_fp m_rfn_low;
vector_fp m_rfn_high;
bool m_ROP_ok;
doublereal m_temp;
doublereal m_pres; //!< Last pressure at which rates were evaluated
vector_fp falloff_work;
vector_fp concm_3b_values;
vector_fp concm_falloff_values;
//!@}
vector_fp m_conc;
void processFalloffReactions();
vector_fp m_grt;
private:
void addElementaryReaction(ReactionData& r);
void addThreeBodyReaction(ReactionData& r);
void addFalloffReaction(ReactionData& r);
void addPlogReaction(ReactionData& r);

249
src/kinetics/AqueousKinetics.cpp
View file

@ -3,7 +3,6 @@
*
* Homogeneous kinetics in an aqueous phase, either condensed
* or dilute in salts
*
*/
/*
* Copyright (2006) Sandia Corporation. Under the terms of
@ -20,58 +19,8 @@ namespace Cantera
{
AqueousKinetics::AqueousKinetics(thermo_t* thermo) :
m_nfall(0),
m_nirrev(0),
m_nrev(0),
m_ROP_ok(false),
m_temp(0.0),
m_finalized(false)
BulkKinetics(thermo)
{
if (thermo != 0) {
addPhase(*thermo);
}
}
AqueousKinetics::AqueousKinetics(const AqueousKinetics& right) :
m_nfall(0),
m_nirrev(0),
m_nrev(0),
m_ROP_ok(false),
m_temp(0.0),
m_finalized(false)
{
*this = right;
}
AqueousKinetics& AqueousKinetics::operator=(const AqueousKinetics& right)
{
if (this == &right) {
return *this;
}
Kinetics::operator=(right);
m_nfall = right.m_nfall;
m_rates = right.m_rates;
m_irrev = right.m_irrev;
m_nirrev = right.m_nirrev;
m_nrev = right.m_nrev;
m_rxntype = right.m_rxntype;
m_dn = right.m_dn;
m_revindex = right.m_revindex;
m_ROP_ok = right.m_ROP_ok;
m_temp = right.m_temp;
m_conc = right.m_conc;
m_grt = right.m_grt;
m_finalized = right.m_finalized;
throw CanteraError("GasKinetics::operator=()",
"Unfinished implementation");
return *this;
}
Kinetics* AqueousKinetics::duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const
@ -81,19 +30,10 @@ Kinetics* AqueousKinetics::duplMyselfAsKinetics(const std::vector<thermo_t*> & t
return gK;
}
void AqueousKinetics::update_T()
{
}
void AqueousKinetics::update_C()
{
}
void AqueousKinetics::_update_rates_T()
{
doublereal T = thermo().temperature();
doublereal logT = log(T);
m_rates.update(T, logT, &m_rfn[0]);
m_rates.update(T, log(T), &m_rfn[0]);
m_temp = T;
updateKc();
@ -103,7 +43,6 @@ void AqueousKinetics::_update_rates_T()
void AqueousKinetics::_update_rates_C()
{
thermo().getActivityConcentrations(&m_conc[0]);
m_ROP_ok = false;
}
@ -122,12 +61,12 @@ void AqueousKinetics::updateKc()
m_rxnstoich.getRevReactionDelta(m_ii, &m_grt[0], &m_rkcn[0]);
doublereal rrt = 1.0/(GasConstant * thermo().temperature());
for (size_t i = 0; i < m_nrev; i++) {
for (size_t i = 0; i < m_revindex.size(); i++) {
size_t irxn = m_revindex[i];
m_rkcn[irxn] = exp(m_rkcn[irxn]*rrt);
}
for (size_t i = 0; i != m_nirrev; ++i) {
for (size_t i = 0; i != m_irrev.size(); ++i) {
m_rkcn[ m_irrev[i] ] = 0.0;
}
}
@ -157,103 +96,6 @@ void AqueousKinetics::getEquilibriumConstants(doublereal* kc)
m_temp = 0.0;
}
void AqueousKinetics::getDeltaGibbs(doublereal* deltaG)
{
/*
* Get the chemical potentials of the species in the
* ideal gas solution.
*/
thermo().getChemPotentials(&m_grt[0]);
/*
* Use the stoichiometric manager to find deltaG for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaG);
}
void AqueousKinetics::getDeltaEnthalpy(doublereal* deltaH)
{
/*
* Get the partial molar enthalpy of all species in the
* ideal gas.
*/
thermo().getPartialMolarEnthalpies(&m_grt[0]);
/*
* Use the stoichiometric manager to find deltaG for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaH);
}
void AqueousKinetics::getDeltaEntropy(doublereal* deltaS)
{
/*
* Get the partial molar entropy of all species in the
* solid solution.
*/
thermo().getPartialMolarEntropies(&m_grt[0]);
/*
* Use the stoichiometric manager to find deltaS for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaS);
}
void AqueousKinetics::getDeltaSSGibbs(doublereal* deltaG)
{
/*
* Get the standard state chemical potentials of the species.
* This is the array of chemical potentials at unit activity
* We define these here as the chemical potentials of the pure
* species at the temperature and pressure of the solution.
*/
thermo().getStandardChemPotentials(&m_grt[0]);
/*
* Use the stoichiometric manager to find deltaG for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaG);
}
void AqueousKinetics::getDeltaSSEnthalpy(doublereal* deltaH)
{
/*
* Get the standard state enthalpies of the species.
* This is the array of chemical potentials at unit activity
* We define these here as the enthalpies of the pure
* species at the temperature and pressure of the solution.
*/
thermo().getEnthalpy_RT(&m_grt[0]);
doublereal RT = thermo().temperature() * GasConstant;
for (size_t k = 0; k < m_kk; k++) {
m_grt[k] *= RT;
}
/*
* Use the stoichiometric manager to find deltaG for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaH);
}
void AqueousKinetics::getDeltaSSEntropy(doublereal* deltaS)
{
/*
* Get the standard state entropy of the species.
* We define these here as the entropies of the pure
* species at the temperature and pressure of the solution.
*/
thermo().getEntropy_R(&m_grt[0]);
doublereal R = GasConstant;
for (size_t k = 0; k < m_kk; k++) {
m_grt[k] *= R;
}
/*
* Use the stoichiometric manager to find deltaS for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaS);
}
void AqueousKinetics::updateROP()
{
_update_rates_T();
@ -272,19 +114,15 @@ void AqueousKinetics::updateROP()
// copy the forward rates to the reverse rates
copy(m_ropf.begin(), m_ropf.end(), m_ropr.begin());
// for reverse rates computed from thermochemistry, multiply
// the forward rates copied into m_ropr by the reciprocals of
// the equilibrium constants
// for reverse rates computed from thermochemistry, multiply the forward
// rates copied into m_ropr by the reciprocals of the equilibrium constants
multiply_each(m_ropr.begin(), m_ropr.end(), m_rkcn.begin());
// multiply ropf by concentration products
m_rxnstoich.multiplyReactants(&m_conc[0], &m_ropf[0]);
//m_reactantStoich.multiply(m_conc.begin(), ropf.begin());
// for reversible reactions, multiply ropr by concentration
// products
// for reversible reactions, multiply ropr by concentration products
m_rxnstoich.multiplyRevProducts(&m_conc[0], &m_ropr[0]);
//m_revProductStoich.multiply(m_conc.begin(), ropr.begin());
for (size_t j = 0; j != m_ii; ++j) {
m_ropnet[j] = m_ropf[j] - m_ropr[j];
@ -309,84 +147,13 @@ void AqueousKinetics::getFwdRateConstants(doublereal* kfwd)
}
}
void AqueousKinetics::getRevRateConstants(doublereal* krev,
bool doIrreversible)
{
/*
* go get the forward rate constants. -> note, we don't
* really care about speed or redundancy in these
* informational routines.
*/
getFwdRateConstants(krev);
if (doIrreversible) {
getEquilibriumConstants(&m_ropnet[0]);
for (size_t i = 0; i < m_ii; i++) {
krev[i] /= m_ropnet[i];
}
} else {
/*
* m_rkcn[] is zero for irreversible reactions
*/
for (size_t i = 0; i < m_ii; i++) {
krev[i] *= m_rkcn[i];
}
}
}
void AqueousKinetics::addReaction(ReactionData& r)
{
if (r.reactionType == ELEMENTARY_RXN) {
addElementaryReaction(r);
}
m_dn.push_back(accumulate(r.pstoich.begin(), r.pstoich.end(), 0.0) -
accumulate(r.rstoich.begin(), r.rstoich.end(), 0.0));
if (r.reversible) {
m_revindex.push_back(nReactions());
m_nrev++;
} else {
m_irrev.push_back(nReactions());
m_nirrev++;
}
Kinetics::addReaction(r);
}
void AqueousKinetics::addElementaryReaction(ReactionData& r)
{
// install rate coeff calculator
m_rates.install(nReactions(), r);
}
void AqueousKinetics::init()
{
m_kk = thermo().nSpecies();
m_rrxn.resize(m_kk);
m_prxn.resize(m_kk);
m_conc.resize(m_kk);
m_grt.resize(m_kk);
}
void AqueousKinetics::finalize()
{
if (!m_finalized) {
m_finalized = true;
// Guarantee that these arrays can be converted to double* even in the
// special case where there are no reactions defined.
if (!m_ii) {
m_perturb.resize(1, 1.0);
m_ropf.resize(1, 0.0);
m_ropr.resize(1, 0.0);
m_ropnet.resize(1, 0.0);
m_rkcn.resize(1, 0.0);
}
}
}
bool AqueousKinetics::ready() const
{
return m_finalized;
BulkKinetics::addReaction(r);
}
}

164
src/kinetics/BulkKinetics.cpp Normal file
View file

@ -0,0 +1,164 @@
#include "cantera/kinetics/BulkKinetics.h"
namespace Cantera
{
BulkKinetics::BulkKinetics(thermo_t* thermo) :
m_ROP_ok(false),
m_temp(0.0),
m_finalized(false)
{
if (thermo) {
addPhase(*thermo);
}
}
Kinetics* BulkKinetics::duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const
{
BulkKinetics* kin = new BulkKinetics(*this);
kin->assignShallowPointers(tpVector);
return kin;
}
bool BulkKinetics::isReversible(size_t i) {
if (std::find(m_revindex.begin(), m_revindex.end(), i)
< m_revindex.end()) {
return true;
} else {
return false;
}
}
void BulkKinetics::getDeltaGibbs(doublereal* deltaG)
{
// Get the chemical potentials of the species in the ideal gas solution.
thermo().getChemPotentials(&m_grt[0]);
// Use the stoichiometric manager to find deltaG for each reaction.
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaG);
}
void BulkKinetics::getDeltaEnthalpy(doublereal* deltaH)
{
// Get the partial molar enthalpy of all species in the ideal gas.
thermo().getPartialMolarEnthalpies(&m_grt[0]);
// Use the stoichiometric manager to find deltaH for each reaction.
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaH);
}
void BulkKinetics::getDeltaEntropy(doublereal* deltaS)
{
// Get the partial molar entropy of all species in the solid solution.
thermo().getPartialMolarEntropies(&m_grt[0]);
// Use the stoichiometric manager to find deltaS for each reaction.
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaS);
}
void BulkKinetics::getDeltaSSGibbs(doublereal* deltaG)
{
// Get the standard state chemical potentials of the species. This is the
// array of chemical potentials at unit activity. We define these here as
// the chemical potentials of the pure species at the temperature and
// pressure of the solution.
thermo().getStandardChemPotentials(&m_grt[0]);
// Use the stoichiometric manager to find deltaG for each reaction.
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaG);
}
void BulkKinetics::getDeltaSSEnthalpy(doublereal* deltaH)
{
// Get the standard state enthalpies of the species.
thermo().getEnthalpy_RT(&m_grt[0]);
doublereal RT = thermo().temperature() * GasConstant;
for (size_t k = 0; k < m_kk; k++) {
m_grt[k] *= RT;
}
// Use the stoichiometric manager to find deltaH for each reaction.
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaH);
}
void BulkKinetics::getDeltaSSEntropy(doublereal* deltaS)
{
// Get the standard state entropy of the species. We define these here as
// the entropies of the pure species at the temperature and pressure of the
// solution.
thermo().getEntropy_R(&m_grt[0]);
doublereal R = GasConstant;
for (size_t k = 0; k < m_kk; k++) {
m_grt[k] *= R;
}
// Use the stoichiometric manager to find deltaS for each reaction.
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaS);
}
void BulkKinetics::getRevRateConstants(doublereal* krev, bool doIrreversible)
{
/*
* go get the forward rate constants. -> note, we don't
* really care about speed or redundancy in these
* informational routines.
*/
getFwdRateConstants(krev);
if (doIrreversible) {
getEquilibriumConstants(&m_ropnet[0]);
for (size_t i = 0; i < m_ii; i++) {
krev[i] /= m_ropnet[i];
}
} else {
// m_rkcn[] is zero for irreversible reactions
for (size_t i = 0; i < m_ii; i++) {
krev[i] *= m_rkcn[i];
}
}
}
void BulkKinetics::addReaction(ReactionData& r)
{
m_dn.push_back(accumulate(r.pstoich.begin(), r.pstoich.end(), 0.0) -
accumulate(r.rstoich.begin(), r.rstoich.end(), 0.0));
if (r.reversible) {
m_revindex.push_back(nReactions());
} else {
m_irrev.push_back(nReactions());
}
Kinetics::addReaction(r);
}
void BulkKinetics::addElementaryReaction(ReactionData& r)
{
m_rates.install(nReactions(), r);
}
void BulkKinetics::init()
{
m_kk = thermo().nSpecies();
m_rrxn.resize(m_kk);
m_prxn.resize(m_kk);
m_conc.resize(m_kk);
m_grt.resize(m_kk);
}
void BulkKinetics::finalize()
{
m_finalized = true;
// Guarantee that these arrays can be converted to double* even in the
// special case where there are no reactions defined.
if (!m_ii) {
m_perturb.resize(1, 1.0);
m_ropf.resize(1, 0.0);
m_ropr.resize(1, 0.0);
m_ropnet.resize(1, 0.0);
m_rkcn.resize(1, 0.0);
}
}
bool BulkKinetics::ready() const
{
return m_finalized;
}
}

257
src/kinetics/GasKinetics.cpp
View file

@ -13,73 +13,13 @@ using namespace std;
namespace Cantera
{
GasKinetics::GasKinetics(thermo_t* thermo) :
BulkKinetics(thermo),
m_nfall(0),
m_nirrev(0),
m_nrev(0),
m_logp_ref(0.0),
m_logc_ref(0.0),
m_logStandConc(0.0),
m_ROP_ok(false),
m_temp(0.0),
m_pres(0.0),
m_finalized(false)
m_pres(0.0)
{
if (thermo != 0) {
addPhase(*thermo);
}
}
GasKinetics::GasKinetics(const GasKinetics& right)
{
*this = right;
}
GasKinetics& GasKinetics::operator=(const GasKinetics& right)
{
if (this == &right) {
return *this;
}
Kinetics::operator=(right);
m_nfall = right.m_nfall;
m_fallindx = right.m_fallindx;
m_falloff_low_rates = right.m_falloff_low_rates;
m_falloff_high_rates = right.m_falloff_high_rates;
m_rates = right.m_rates;
m_falloffn = right.m_falloffn;
m_3b_concm = right.m_3b_concm;
m_falloff_concm = right.m_falloff_concm;
m_irrev = right.m_irrev;
m_plog_rates = right.m_plog_rates;
m_cheb_rates = right.m_cheb_rates;
m_nirrev = right.m_nirrev;
m_nrev = right.m_nrev;
m_rrxn = right.m_rrxn;
m_prxn = right.m_prxn;
m_dn = right.m_dn;
m_revindex = right.m_revindex;
m_logp_ref = right.m_logp_ref;
m_logc_ref = right.m_logc_ref;
m_logStandConc = right.m_logStandConc;
m_rfn_low = right.m_rfn_low;
m_rfn_high = right.m_rfn_high;
m_ROP_ok = right.m_ROP_ok;
m_temp = right.m_temp;
falloff_work = right.falloff_work;
concm_3b_values = right.concm_3b_values;
concm_falloff_values = right.concm_falloff_values;
m_conc = right.m_conc;
m_grt = right.m_grt;
m_finalized = right.m_finalized;
throw CanteraError("GasKinetics::operator=()",
"Unfinished implementation");
return *this;
}
Kinetics* GasKinetics::duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const
@ -166,13 +106,13 @@ void GasKinetics::updateKc()
m_rxnstoich.getRevReactionDelta(m_ii, &m_grt[0], &m_rkcn[0]);
doublereal rrt = 1.0/(GasConstant * thermo().temperature());
for (size_t i = 0; i < m_nrev; i++) {
for (size_t i = 0; i < m_revindex.size(); i++) {
size_t irxn = m_revindex[i];
m_rkcn[irxn] = std::min(exp(m_rkcn[irxn]*rrt - m_dn[irxn]*m_logStandConc),
BigNumber);
}
for (size_t i = 0; i != m_nirrev; ++i) {
for (size_t i = 0; i != m_irrev.size(); ++i) {
m_rkcn[ m_irrev[i] ] = 0.0;
}
}
@ -196,103 +136,6 @@ void GasKinetics::getEquilibriumConstants(doublereal* kc)
m_temp = 0.0;
}
void GasKinetics::getDeltaGibbs(doublereal* deltaG)
{
/*
* Get the chemical potentials of the species in the
* ideal gas solution.
*/
thermo().getChemPotentials(&m_grt[0]);
/*
* Use the stoichiometric manager to find deltaG for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaG);
}
void GasKinetics::getDeltaEnthalpy(doublereal* deltaH)
{
/*
* Get the partial molar enthalpy of all species in the
* ideal gas.
*/
thermo().getPartialMolarEnthalpies(&m_grt[0]);
/*
* Use the stoichiometric manager to find deltaG for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaH);
}
void GasKinetics::getDeltaEntropy(doublereal* deltaS)
{
/*
* Get the partial molar entropy of all species in the
* solid solution.
*/
thermo().getPartialMolarEntropies(&m_grt[0]);
/*
* Use the stoichiometric manager to find deltaS for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaS);
}
void GasKinetics::getDeltaSSGibbs(doublereal* deltaG)
{
/*
* Get the standard state chemical potentials of the species.
* This is the array of chemical potentials at unit activity
* We define these here as the chemical potentials of the pure
* species at the temperature and pressure of the solution.
*/
thermo().getStandardChemPotentials(&m_grt[0]);
/*
* Use the stoichiometric manager to find deltaG for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaG);
}
void GasKinetics::getDeltaSSEnthalpy(doublereal* deltaH)
{
/*
* Get the standard state enthalpies of the species.
* This is the array of chemical potentials at unit activity
* We define these here as the enthalpies of the pure
* species at the temperature and pressure of the solution.
*/
thermo().getEnthalpy_RT(&m_grt[0]);
doublereal RT = thermo().temperature() * GasConstant;
for (size_t k = 0; k < m_kk; k++) {
m_grt[k] *= RT;
}
/*
* Use the stoichiometric manager to find deltaG for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaH);
}
void GasKinetics::getDeltaSSEntropy(doublereal* deltaS)
{
/*
* Get the standard state entropy of the species.
* We define these here as the entropies of the pure
* species at the temperature and pressure of the solution.
*/
thermo().getEntropy_R(&m_grt[0]);
doublereal R = GasConstant;
for (size_t k = 0; k < m_kk; k++) {
m_grt[k] *= R;
}
/*
* Use the stoichiometric manager to find deltaS for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], deltaS);
}
void GasKinetics::processFalloffReactions()
{
// use m_ropr for temporary storage of reduced pressure
@ -346,19 +189,15 @@ void GasKinetics::updateROP()
// copy the forward rates to the reverse rates
copy(m_ropf.begin(), m_ropf.end(), m_ropr.begin());
// for reverse rates computed from thermochemistry, multiply
// the forward rates copied into m_ropr by the reciprocals of
// the equilibrium constants
// for reverse rates computed from thermochemistry, multiply the forward
// rates copied into m_ropr by the reciprocals of the equilibrium constants
multiply_each(m_ropr.begin(), m_ropr.end(), m_rkcn.begin());
// multiply ropf by concentration products
m_rxnstoich.multiplyReactants(&m_conc[0], &m_ropf[0]);
//m_reactantStoich.multiply(m_conc.begin(), ropf.begin());
// for reversible reactions, multiply ropr by concentration
// products
// for reversible reactions, multiply ropr by concentration products
m_rxnstoich.multiplyRevProducts(&m_conc[0], &m_ropr[0]);
//m_revProductStoich.multiply(m_conc.begin(), ropr.begin());
for (size_t j = 0; j != m_ii; ++j) {
m_ropnet[j] = m_ropf[j] - m_ropr[j];
@ -389,10 +228,6 @@ void GasKinetics::getFwdRateConstants(doublereal* kfwd)
m_3b_concm.multiply(&m_ropf[0], &concm_3b_values[0]);
}
/*
* This routine is hardcoded to replace some of the values
* of the ropf vector.
*/
if (m_nfall) {
processFalloffReactions();
}
@ -405,28 +240,6 @@ void GasKinetics::getFwdRateConstants(doublereal* kfwd)
}
}
void GasKinetics::getRevRateConstants(doublereal* krev, bool doIrreversible)
{
/*
* go get the forward rate constants. -> note, we don't
* really care about speed or redundancy in these
* informational routines.
*/
getFwdRateConstants(krev);
if (doIrreversible) {
getEquilibriumConstants(&m_ropnet[0]);
for (size_t i = 0; i < m_ii; i++) {
krev[i] /= m_ropnet[i];
}
} else {
// m_rkcn[] is zero for irreversible reactions
for (size_t i = 0; i < m_ii; i++) {
krev[i] *= m_rkcn[i];
}
}
}
void GasKinetics::addReaction(ReactionData& r)
{
switch (r.reactionType) {
@ -451,17 +264,7 @@ void GasKinetics::addReaction(ReactionData& r)
}
// operations common to all reaction types
m_dn.push_back(accumulate(r.pstoich.begin(), r.pstoich.end(), 0.0) -
accumulate(r.rstoich.begin(), r.rstoich.end(), 0.0));
if (r.reversible) {
m_revindex.push_back(nReactions());
m_nrev++;
} else {
m_irrev.push_back(nReactions());
m_nirrev++;
}
Kinetics::addReaction(r);
BulkKinetics::addReaction(r);
}
void GasKinetics::addFalloffReaction(ReactionData& r)
@ -474,17 +277,15 @@ void GasKinetics::addFalloffReaction(ReactionData& r)
m_falloff_low_rates.install(m_nfall, r);
m_rfn_low.push_back(r.rateCoeffParameters[0]);
// add this reaction number to the list of
// falloff reactions
// add this reaction number to the list of falloff reactions
m_fallindx.push_back(nReactions());
// install the enhanced third-body concentration
// calculator for this reaction
// install the enhanced third-body concentration calculator for this
// reaction
m_falloff_concm.install(m_nfall, r.thirdBodyEfficiencies,
r.default_3b_eff);
// install the falloff function calculator for
// this reaction
// install the falloff function calculator for this reaction
m_falloffn.install(m_nfall, r.falloffType, r.reactionType,
r.falloffParameters);
@ -492,14 +293,8 @@ void GasKinetics::addFalloffReaction(ReactionData& r)
++m_nfall;
}
void GasKinetics::addElementaryReaction(ReactionData& r)
{
m_rates.install(nReactions(), r);
}
void GasKinetics::addThreeBodyReaction(ReactionData& r)
{
// install rate coeff calculator
m_rates.install(nReactions(), r);
m_3b_concm.install(nReactions(), r.thirdBodyEfficiencies,
r.default_3b_eff);
@ -507,44 +302,26 @@ void GasKinetics::addThreeBodyReaction(ReactionData& r)
void GasKinetics::addPlogReaction(ReactionData& r)
{
// install rate coefficient calculator
m_plog_rates.install(nReactions(), r);
}
void GasKinetics::addChebyshevReaction(ReactionData& r)
{
// install rate coefficient calculator
m_cheb_rates.install(nReactions(), r);
}
void GasKinetics::init()
{
m_kk = thermo().nSpecies();
m_rrxn.resize(m_kk);
m_prxn.resize(m_kk);
m_conc.resize(m_kk);
m_grt.resize(m_kk);
BulkKinetics::init();
m_logp_ref = log(thermo().refPressure()) - log(GasConstant);
}
void GasKinetics::finalize()
{
if (!m_finalized) {
falloff_work.resize(m_falloffn.workSize());
concm_3b_values.resize(m_3b_concm.workSize());
concm_falloff_values.resize(m_falloff_concm.workSize());
m_finalized = true;
// Guarantee that these arrays can be converted to double* even in the
// special case where there are no reactions defined.
if (!m_ii) {
m_perturb.resize(1, 1.0);
m_ropf.resize(1, 0.0);
m_ropr.resize(1, 0.0);
m_ropnet.resize(1, 0.0);
m_rkcn.resize(1, 0.0);
}
}
BulkKinetics::finalize();
falloff_work.resize(m_falloffn.workSize());
concm_3b_values.resize(m_3b_concm.workSize());
concm_falloff_values.resize(m_falloff_concm.workSize());
}
bool GasKinetics::ready() const