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

EdgeKinetics now derives from InterfaceKinetics, and InterfaceKinetics has the beta parameter for electrochemistry

Browse source
This commit is contained in:
Dave Goodwin 2007-05-18 10:54:54 +00:00
parent 5b75d6bae7
commit 39d60dd2d1
11 changed files with 126 additions and 923 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

8
CMakeLists.txt
View file

@ -1,7 +1,11 @@
PROJECT (Cantera)
#----------------------------
# user-configurable settings
#----------------------------
INCLUDE (config.cmake)
#operating system
SET(os ${CMAKE_SYSTEM_NAME})
if (os STREQUAL "Darwin")
@ -15,6 +19,8 @@ else (PYTHON_CMD STREQUAL "default")
SET( PYTHON_EXE ${PYTHON_CMD} )
endif (PYTHON_CMD STREQUAL "default")
INCLUDE(cmake/thermo.cmake)
INCLUDE(cmake/fortran.cmake)
# configuration
@ -22,6 +28,8 @@ CONFIGURE_FILE (
${PROJECT_SOURCE_DIR}/config.h_cmake.in
${PROJECT_BINARY_DIR}/config.h )
include_directories(${PROJECT_BINARY_DIR})
set(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/build/lib/${CMAKE_SYSTEM_PROCESSOR}-${CMAKE_SYSTEM_VERSION})
add_subdirectory(Cantera)

18
Cantera/python/Cantera/Kinetics.py
View file

@ -265,6 +265,24 @@ class Kinetics:
def sourceTerms(self):
return _cantera.kin_getarray(self.ckin,80)
def delta_H(self):
return _cantera.kin_getarray(self.ckin,90)
def delta_G(self):
return _cantera.kin_getarray(self.ckin,91)
def delta_S(self):
return _cantera.kin_getarray(self.ckin,92)
def delta_H0(self):
return _cantera.kin_getarray(self.ckin,93)
def delta_G0(self):
return _cantera.kin_getarray(self.ckin,94)
def delta_S0(self):
return _cantera.kin_getarray(self.ckin,95)
def multiplier(self,i):
return _cantera.kin_multiplier(self.ckin,i)

20
Cantera/python/src/ctkinetics_methods.cpp
View file

@ -135,7 +135,7 @@ kin_getarray(PyObject *self, PyObject *args)
int nrxns = kin_nReactions(kin);
int nsp = kin_nSpecies(kin);
int ix;
if (job < 45) ix = nrxns; else ix = nsp;
if (job < 45 || job >= 90) ix = nrxns; else ix = nsp;
PyArrayObject* x =
(PyArrayObject*)PyArray_FromDims(1, &ix, PyArray_DOUBLE);
@ -177,7 +177,23 @@ kin_getarray(PyObject *self, PyObject *args)
case 80:
iok = kin_getSourceTerms(kin, nsp, xd);
break;
case 90:
iok = kin_getDelta(kin, 0, nrxns, xd);
break;
case 91:
iok = kin_getDelta(kin, 1, nrxns, xd);
break;
case 92:
iok = kin_getDelta(kin, 2, nrxns, xd);
break;
case 93:
iok = kin_getDelta(kin, 3, nrxns, xd);
break;
case 94:
iok = kin_getDelta(kin, 4, nrxns, xd);
break;
case 95:
iok = kin_getDelta(kin, 5, nrxns, xd);
break;
default:
;

531
Cantera/src/kinetics/EdgeKinetics.cpp
View file

@ -1,531 +0,0 @@
/**
* @file EdgeKinetics.cpp
*
*/
// Copyright 2002 California Institute of Technology
// turn off warnings under Windows
#ifdef WIN32
#pragma warning(disable:4786)
#pragma warning(disable:4503)
#endif
#include "EdgeKinetics.h"
#include "SurfPhase.h"
#include "ReactionData.h"
//#include "StoichManager.h"
#include "RateCoeffMgr.h"
#include <iostream>
using namespace std;
namespace Cantera {
//////////////////////////////////////////////////////////////////
/**
* Construct an empty EdgeKinetics reaction mechanism.
*/
EdgeKinetics::
EdgeKinetics() :
Kinetics(),
m_kk(0),
m_redo_rates(false),
m_nirrev(0),
m_nrev(0),
m_finalized(false),
m_has_electrochem_rxns(false)
{
m_kdata = new EdgeKineticsData;
m_kdata->m_temp = 0.0;
}
/**
* Destructor
*/
EdgeKinetics::
~EdgeKinetics(){
delete m_kdata;
}
/**
* Update properties that depend on temperature
*
*/
void EdgeKinetics::
_update_rates_T() {
_update_rates_phi();
doublereal T = thermo(surfacePhaseIndex()).temperature();
if (T != m_kdata->m_temp || m_redo_rates) {
m_kdata->m_logtemp = log(T);
m_rates.update(T, m_kdata->m_logtemp, DATA_PTR(m_kdata->m_rfn));
if (m_has_electrochem_rxns)
applyButlerVolmerCorrection(DATA_PTR(m_kdata->m_rfn));
m_kdata->m_temp = T;
updateKc();
m_kdata->m_ROP_ok = false;
m_redo_rates = false;
}
}
void EdgeKinetics::
_update_rates_phi() {
int np = nPhases();
for (int n = 0; n < np; n++) {
if (thermo(n).electricPotential() != m_phi[n]) {
m_phi[n] = thermo(n).electricPotential();
m_redo_rates = true;
}
}
}
/**
* Update properties that depend on concentrations. This method
* fills out the array of generalized concentrations by calling
* method getActivityConcentrations for each phase, which classes
* representing phases should overload to return the appropriate
* quantities.
*/
void EdgeKinetics::
_update_rates_C() {
int n;
//m_rates.update(m_kdata->m_temp,
// m_kdata->m_logtemp, m_kdata->m_rfn.begin());
int np = nPhases();
for (n = 0; n < np; n++) {
thermo(n).getActivityConcentrations(DATA_PTR(m_conc) + m_start[n]);
}
m_kdata->m_ROP_ok = false;
}
/**
* Update the equilibrium constants in molar units for all
* reversible reactions. Irreversible reactions have their
* equilibrium constant set to zero.
*/
void EdgeKinetics::updateKc() {
int i, irxn;
vector_fp& m_rkc = m_kdata->m_rkcn;
fill(m_rkc.begin(), m_rkc.end(), 0.0);
b
if (m_nrev > 0) {
int n, nsp, k, ik=0;
doublereal rt = GasConstant*thermo(0).temperature();
doublereal rrt = 1.0/rt;
int np = nPhases();
for (n = 0; n < np; n++) {
thermo(n).getStandardChemPotentials(DATA_PTR(m_mu0) + m_start[n]);
nsp = thermo(n).nSpecies();
for (k = 0; k < nsp; k++) {
m_mu0[ik] -= rt*thermo(n).logStandardConc(k);
m_mu0[ik] += Faraday * m_phi[n] * thermo(n).charge(k);
ik++;
}
}
// compute Delta mu^0 for all reversible reactions
m_rxnstoich.getRevReactionDelta(m_ii, DATA_PTR(m_mu0),
DATA_PTR(m_rkc));
for (i = 0; i < m_nrev; i++) {
irxn = m_revindex[i];
m_rkc[irxn] = exp(m_rkc[irxn]*rrt);
}
for (i = 0; i != m_nirrev; ++i) {
m_rkc[ m_irrev[i] ] = 0.0;
}
}
}
void EdgeKinetics::checkPartialEquil() {
int i, irxn;
vector_fp dmu(nTotalSpecies(), 0.0);
vector_fp rmu(nReactions(), 0.0);
if (m_nrev > 0) {
int n, nsp, k, ik=0;
doublereal rt = GasConstant*thermo(0).temperature();
doublereal rrt = 1.0/rt;
int np = nPhases();
for (n = 0; n < np; n++) {
thermo(n).getChemPotentials(DATA_PTR(dmu) + m_start[n]);
nsp = thermo(n).nSpecies();
for (k = 0; k < nsp; k++) {
dmu[ik] += Faraday * m_phi[n] * thermo(n).charge(k);
cout << thermo(n).speciesName(k) << " " << dmu[ik] << endl;
ik++;
}
}
// compute Delta mu^ for all reversible reactions
m_rxnstoich.getRevReactionDelta(m_ii, DATA_PTR(dmu),
DATA_PTR(rmu));
//m_reactantStoich.decrementReactions(DATA_PTR(dmu), DATA_PTR(rmu));
//m_revProductStoich.incrementReactions(DATA_PTR(dmu), DATA_PTR(rmu));
for (i = 0; i < m_nrev; i++) {
irxn = m_revindex[i];
cout << "Reaction " << irxn << " " << exp(rmu[irxn]*rrt) << endl;
}
}
}
/**
* Get the equilibrium constants of all reactions, whether
* reversible or not.
*/
void EdgeKinetics::getEquilibriumConstants(doublereal* kc) {
int i;
int n, nsp, k, ik=0;
doublereal rt = GasConstant*thermo(0).temperature();
doublereal rrt = 1.0/rt;
int np = nPhases();
for (n = 0; n < np; n++) {
thermo(n).getStandardChemPotentials(DATA_PTR(m_mu0) + m_start[n]);
nsp = thermo(n).nSpecies();
for (k = 0; k < nsp; k++) {
m_mu0[ik] -= rt*thermo(n).logStandardConc(k);
m_mu0[ik] += Faraday * m_phi[n] * thermo(n).charge(k);
ik++;
}
}
fill(kc, kc + m_ii, 0.0);
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_mu0), kc);
// m_reactantStoich.decrementReactions(DATA_PTR(m_mu0), kc);
//m_revProductStoich.incrementReactions(DATA_PTR(m_mu0), kc);
//m_irrevProductStoich.incrementReactions(DATA_PTR(m_mu0), kc);
for (i = 0; i < m_ii; i++) {
kc[i] = exp(-kc[i]*rrt);
}
}
/**
* For reactions that transfer charge across a potential difference,
* the activation energies are modified by the potential difference.
* (see, for example, Baird and Falkner, "Electrochemical Methods").
* This method applies this correction.
*/
void EdgeKinetics::applyButlerVolmerCorrection(doublereal* kf) {
int i;
int n, nsp, k, ik=0;
doublereal rt = GasConstant*thermo(0).temperature();
doublereal rrt = 1.0/rt;
int np = nPhases();
// compute the electrical potential energy of each species
for (n = 0; n < np; n++) {
nsp = thermo(n).nSpecies();
for (k = 0; k < nsp; k++) {
m_pot[ik] = Faraday*thermo(n).charge(k)*m_phi[n];
ik++;
}
}
// compute the change in electrical potential energy for each
// reaction. This will only be non-zero if a potential
// difference is present.
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_pot),
DATA_PTR(m_rwork));
// modify the reaction rates. Only modify those with a
// non-zero activation energy, and do not decrease the
// activation energy below zero.
doublereal ea, eamod;
int nct = m_beta.size();
int irxn;
for (i = 0; i < nct; i++) {
irxn = m_ctrxn[i];
eamod = m_beta[i]*m_rwork[irxn];
if (eamod != 0.0 && m_E[i] != 0.0) {
ea = GasConstant * m_E[i];
if (eamod + ea < 0.0) {
writelog("Warning: act energy mod too large!\n");
eamod = -ea;
}
kf[irxn] *= exp(-eamod*rrt);
}
}
}
/**
* Update the rates of progress of the reactions in the reaciton
* mechanism. This routine operates on internal data.
*/
void EdgeKinetics::updateROP() {
_update_rates_T();
_update_rates_C();
if (m_kdata->m_ROP_ok) return;
const vector_fp& rf = m_kdata->m_rfn;
const vector_fp& m_rkc = m_kdata->m_rkcn;
array_fp& ropf = m_kdata->m_ropf;
array_fp& ropr = m_kdata->m_ropr;
array_fp& ropnet = m_kdata->m_ropnet;
// copy rate coefficients into ropf
copy(rf.begin(), rf.end(), ropf.begin());
// multiply by perturbation factor
multiply_each(ropf.begin(), ropf.end(), m_perturb.begin());
// copy the forward rates to the reverse rates
copy(ropf.begin(), ropf.end(), ropr.begin());
// for reverse rates computed from thermochemistry, multiply
// the forward rates copied into m_ropr by the reciprocals of
// the equilibrium constants
multiply_each(ropr.begin(), ropr.end(), m_rkc.begin());
// multiply ropf by concentration products
m_reactantStoich.multiply(DATA_PTR(m_conc), DATA_PTR(ropf));
// for reversible reactions, multiply ropr by concentration
// products
m_revProductStoich.multiply(DATA_PTR(m_conc), DATA_PTR(ropr));
// do global reactions
//m_globalReactantStoich.power(DATA_PTR(m_conc), ropf.begin());
for (int j = 0; j != m_ii; ++j) {
ropnet[j] = ropf[j] - ropr[j];
}
m_kdata->m_ROP_ok = true;
}
/**
* Add a single reaction to the mechanism. This routine
* must be called after init() and before finalize().
* This function branches on the types of reactions allowed
* by the interfaceKinetics manager in order to install
* the reaction correctly in the manager.
* The manager allows the following reaction types
* Elementary
* Surface
* Global
* There is no difference between elementary and surface
* reactions.
*/
void EdgeKinetics::
addReaction(const ReactionData& r) {
int nr = r.reactants.size();
// a global reaction is idnetified as one with
// a reactant stoichiometric coefficient not equal
// to the molecularity for some reactant
bool isglobal = false;
for (int n = 0; n < nr; n++) {
if (r.rstoich[n] != int(r.order[n])) {
isglobal = true; break;
}
}
if (isglobal)
addGlobalReaction(r);
else
addElementaryReaction(r);
installReagents( r );
installGroups(reactionNumber(), r.rgroups, r.pgroups);
incrementRxnCount();
m_rxneqn.push_back(r.equation);
}
void EdgeKinetics::
addElementaryReaction(const ReactionData& r) {
int iloc;
// install rate coeff calculator
vector_fp rp = r.rateCoeffParameters;
// coverage dependence
int ncov = r.cov.size();
for (int m = 0; m < ncov; m++) rp.push_back(r.cov[m]);
iloc = m_rates.install( reactionNumber(), r.rateCoeffType, rp.size(),
DATA_PTR(rp) );
// store activation energy
m_E.push_back(r.rateCoeffParameters[2]);
if (r.beta > 0.0) {
m_has_electrochem_rxns = true;
m_beta.push_back(r.beta);
m_ctrxn.push_back(reactionNumber());
}
// add constant term to rate coeff value vector
m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]);
registerReaction( reactionNumber(), ELEMENTARY_RXN, iloc);
}
void EdgeKinetics::
addGlobalReaction(const ReactionData& r) {
int iloc;
// install rate coeff calculator
vector_fp rp = r.rateCoeffParameters;
int ncov = r.cov.size();
for (int m = 0; m < ncov; m++) rp.push_back(r.cov[m]);
iloc = m_rates.install( reactionNumber(),
r.rateCoeffType, rp.size(),
DATA_PTR(rp) );
// add constant term to rate coeff value vector
m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]);
int nr = r.order.size();
vector_fp ordr(nr);
for (int n = 0; n < nr; n++) {
ordr[n] = r.order[n] - r.rstoich[n];
}
m_globalReactantStoich.add( reactionNumber(),
r.reactants, ordr);
registerReaction( reactionNumber(), GLOBAL_RXN, iloc);
}
void EdgeKinetics::installReagents(const ReactionData& r) {
m_kdata->m_ropf.push_back(0.0); // extend by one for new rxn
m_kdata->m_ropr.push_back(0.0);
m_kdata->m_ropnet.push_back(0.0);
int n, ns, m;
doublereal nsFlt;
int rnum = reactionNumber();
vector_int rk;
int nr = r.reactants.size();
for (n = 0; n < nr; n++) {
nsFlt = r.rstoich[n];
ns = (int) nsFlt;
if ((doublereal) ns != nsFlt) {
if (ns < 1) ns = 1;
}
m_rrxn[r.reactants[n]][rnum] = ns;
for (m = 0; m < ns; m++) {
rk.push_back(r.reactants[n]);
}
}
m_reactants.push_back(rk);
vector_int pk;
int np = r.products.size();
for (n = 0; n < np; n++) {
nsFlt = r.pstoich[n];
ns = (int) nsFlt;
if ((doublereal) ns != nsFlt) {
if (ns < 1) ns = 1;
}
m_prxn[r.products[n]][rnum] = ns;
for (m = 0; m < ns; m++) {
pk.push_back(r.products[n]);
}
}
m_products.push_back(pk);
m_kdata->m_rkcn.push_back(0.0);
m_reactantStoich.add( reactionNumber(), rk);
if (r.reversible) {
m_revProductStoich.add(reactionNumber(), pk);
//m_dn.push_back(pk.size() - rk.size());
m_revindex.push_back(reactionNumber());
m_nrev++;
}
else {
m_irrevProductStoich.add(reactionNumber(), pk);
//m_dn.push_back(pk.size() - rk.size());
m_irrev.push_back( reactionNumber() );
m_nirrev++;
}
}
void EdgeKinetics::installGroups(int irxn,
const vector<grouplist_t>& r, const vector<grouplist_t>& p) {
if (!r.empty()) {
m_rgroups[reactionNumber()] = r;
m_pgroups[reactionNumber()] = p;
}
}
/**
* Prepare the class for the addition of reactions. This function
* must be called after instantiation of the class, but before
* any reactions are actually added to the mechanism.
* This function calculates m_kk the number of species in all
* phases participating in the reaction mechanism. We don't know
* m_kk previously, before all phases have been added.
*/
void EdgeKinetics::init() {
int n;
m_kk = 0;
int np = nPhases();
for (n = 0; n < np; n++) {
m_kk += thermo(n).nSpecies();
}
m_rrxn.resize(m_kk);
m_prxn.resize(m_kk);
m_conc.resize(m_kk);
m_mu0.resize(m_kk);
m_pot.resize(m_kk, 0.0);
m_phi.resize(np, 0.0);
}
/**
* Finish adding reactions and prepare for use. This function
* must be called after all reactions are entered into the mechanism
* and before the mechanism is used to calculate reaction rates.
*
* Here, we resize work arrays based on the number of reactions,
* since we don't know this number up to now.
*/
void EdgeKinetics::finalize() {
m_rwork.resize(nReactions());
m_finalized = true;
}
bool EdgeKinetics::ready() const {
return (m_finalized);
}
}

353
Cantera/src/kinetics/EdgeKinetics.h
View file

@ -12,51 +12,11 @@
#ifndef CT_EDGEKINETICS_H
#define CT_EDGEKINETICS_H
#include <fstream>
#include <math.h>
#include <map>
#include <stdlib.h>
#include "mix_defs.h"
#include "Kinetics.h"
#include "utilities.h"
#include "RateCoeffMgr.h"
#include "StoichManager.h"
#include "InterfaceKinetics.h"
namespace Cantera {
// forward references
class ReactionData;
class EdgeKineticsData;
class ThermoPhase;
class SurfPhase;
class ImplicitSurfChem;
/**
* Holds mechanism-specific data.
*/
class EdgeKineticsData {
public:
EdgeKineticsData() :
m_ROP_ok(false),
m_temp(0.0), m_logtemp(0.0)
{}
virtual ~EdgeKineticsData(){}
doublereal m_logp0, m_logc0;
array_fp m_ropf, m_ropr, m_ropnet;
//array_fp m_rfn_low, m_rfn_high;
bool m_ROP_ok;
doublereal m_temp, m_logtemp;
vector_fp m_rfn;
vector_fp m_rkcn;
};
class EdgeKinetics : public Kinetics {
class EdgeKinetics : public InterfaceKinetics {
public:
@ -64,10 +24,10 @@ namespace Cantera {
* Constructor
*
*/
EdgeKinetics();
EdgeKinetics() : InterfaceKinetics() {}
/// Destructor.
virtual ~EdgeKinetics();
virtual ~EdgeKinetics() {}
/**
* Identifies the subclass of the Kinetics manager type.
@ -81,311 +41,8 @@ namespace Cantera {
*/
virtual int type() { return cEdgeKinetics; }
/**
* Set the electric potential in the nth phase
*
* @param n phase Index in this kinetics object.
* @param V Electric potential (volts)
*/
void setElectricPotential(int n, doublereal V) {
thermo(n).setElectricPotential(V);
m_redo_rates = true;
}
/**
* @name Reaction Rates Of Progress
*/
//@{
/**
* Forward rates of progress.
* Return the forward rates of progress in array fwdROP, which
* must be dimensioned at least as large as the total number
* of reactions.
* Units are kmol/m2/s
*/
virtual void getFwdRatesOfProgress(doublereal* fwdROP) {
updateROP();
std::copy(m_kdata->m_ropf.begin(), m_kdata->m_ropf.end(), fwdROP);
}
/**
* Reverse rates of progress.
* Return the reverse rates of progress in array revROP, which
* must be dimensioned at least as large as the total number
* of reactions.
* Units are kmol/m2/s
*/
virtual void getRevRatesOfProgress(doublereal* revROP) {
updateROP();
std::copy(m_kdata->m_ropr.begin(), m_kdata->m_ropr.end(), revROP);
}
/**
* Net rates of progress. Return the net (forward - reverse)
* rates of progress in array netROP, which must be
* dimensioned at least as large as the total number of
* reactions.
* Units are kmol/m2/s
*/
virtual void getNetRatesOfProgress(doublereal* netROP) {
updateROP();
std::copy(m_kdata->m_ropnet.begin(), m_kdata->m_ropnet.end(), netROP);
}
/**
* Equilibrium constants. Return the equilibrium constants of
* the reactions in concentration units in array kc, which
* must be dimensioned at least as large as the total number
* of reactions.
*/
virtual void getEquilibriumConstants(doublereal* kc);
//@}
/**
* @name Species Production Rates
*/
//@{
/**
* Species creation rates [kmol/m^2/s]. Return the species
* creation rates in array cdot, which must be
* dimensioned at least as large as the total number of
* species in all phases of the kinetics
* model
*
*/
virtual void getCreationRates(doublereal* cdot) {
updateROP();
std::fill(cdot, cdot + m_kk, 0.0);
m_revProductStoich.incrementSpecies(
&m_kdata->m_ropf[0], cdot);
m_irrevProductStoich.incrementSpecies(
&m_kdata->m_ropf[0], cdot);
m_reactantStoich.incrementSpecies(
&m_kdata->m_ropr[0], cdot);
}
/**
* Species destruction rates [kmol/m^2/s]. Return the species
* destruction rates in array ddot, which must be
* dimensioned at least as large as the total number of
* species in all phases of the kinetics
* model
*
*/
virtual void getDestructionRates(doublereal* ddot) {
updateROP();
std::fill(ddot, ddot + m_kk, 0.0);
m_revProductStoich.incrementSpecies(
&m_kdata->m_ropr[0], ddot);
m_reactantStoich.incrementSpecies(
&m_kdata->m_ropf[0], ddot);
}
/**
* Species net production rates [kmol/m^2/s]. Return the species
* net production rates (creation - destruction) in array
* wdot, which must be dimensioned at least as large as the
* total number of species in all phases of the kinetics
* model
*/
virtual void getNetProductionRates(doublereal* net) {
updateROP();
std::fill(net, net + m_kk, 0.0);
m_revProductStoich.incrementSpecies(
&m_kdata->m_ropnet[0], net);
m_irrevProductStoich.incrementSpecies(
&m_kdata->m_ropnet[0], net);
m_reactantStoich.decrementSpecies(
&m_kdata->m_ropnet[0], net);
}
//@}
/**
* @name Reaction Mechanism Informational Query Routines
*/
//@{
/**
* Stoichiometric coefficient of species k as a reactant in
* reaction i.
*/
virtual doublereal reactantStoichCoeff(int k, int i) const {
return m_rrxn[k][i];
}
/**
* Stoichiometric coefficient of species k as a product in
* reaction i.
*/
virtual doublereal productStoichCoeff(int k, int i) const {
return m_prxn[k][i];
}
/**
* Flag specifying the type of reaction. The legal values and
* their meaning are specific to the particular kinetics
* manager.
*/
virtual int reactionType(int i) const {
return m_index[i].first;
}
/**
* True if reaction i has been declared to be reversible. If
* isReversible(i) is false, then the reverse rate of progress
* for reaction i is always zero.
*/
virtual bool isReversible(int i) {
if (std::find(m_revindex.begin(), m_revindex.end(), i)
< m_revindex.end()) return true;
else return false;
}
/**
* Return a string representing the reaction.
*/
virtual std::string reactionString(int i) const {
return m_rxneqn[i];
}
//@}
/**
* @name Reaction Mechanism Construction
*/
//@{
/**
* Prepare the class for the addition of reactions. This function
* must be called after instantiation of the class, but before
* any reactions are actually added to the mechanism.
* This function calculates m_kk the number of species in all
* phases participating in the reaction mechanism. We don't know
* m_kk previously, before all phases have been added.
*/
virtual void init();
/**
* Add a single reaction to the mechanism.
*/
virtual void addReaction(const ReactionData& r);
/**
* Finish adding reactions and prepare for use. This function
* must be called after all reactions are entered into the mechanism
* and before the mechanism is used to calculate reaction rates.
*/
// defined in InterfaceKinetics.cpp
virtual void finalize();
virtual bool ready() const;
void updateROP();
const std::vector<grouplist_t>& reactantGroups(int i)
{ return m_rgroups[i]; }
const std::vector<grouplist_t>& productGroups(int i)
{ return m_pgroups[i]; }
void _update_rates_T();
void _update_rates_phi();
void _update_rates_C();
void checkPartialEquil();
protected:
/**
* m_kk here is the number of species in all of the phases
* that participate in the kinetics mechanism.
*/
int m_kk;
Rate1<SurfaceArrhenius> m_rates;
//Rate1<Arrhenius> m_rates;
bool m_redo_rates;
/**
* Vector of information about reactions in the
* mechanism.
* The key is the reaction index (0 < i < m_ii).
* The first pair is the reactionType of the reaction.
* The second pair is ...
*/
mutable std::map<int, std::pair<int, int> > m_index;
std::vector<int> m_irrev;
StoichManagerN m_reactantStoich;
StoichManagerN m_revProductStoich;
StoichManagerN m_irrevProductStoich;
StoichManagerN m_globalReactantStoich;
int m_nirrev;
/**
* Number of reversible reactions in the mechanism
*/
int m_nrev;
std::map<int, std::vector<grouplist_t> > m_rgroups;
std::map<int, std::vector<grouplist_t> > m_pgroups;
std::vector<int> m_rxntype;
mutable std::vector<std::map<int, doublereal> > m_rrxn;
mutable std::vector<std::map<int, doublereal> > m_prxn;
vector_int m_revindex;
std::vector<std::string> m_rxneqn;
/**
* Temporary data storage used in calculating the rates of
* of reactions.
*/
EdgeKineticsData* m_kdata;
/**
* An array of generalized concentrations
* \f$ C_k \f$ that are defined such that \f$ a_k = C_k /
* C^0_k, \f$ where \f$ C^0_k \f$ is a standard concentration/
* These generalized concentrations are used
* by this kinetics manager class to compute the forward and
* reverse rates of elementary reactions. The "units" for the
* concentrations of each phase depend upon the implementation
* of kinetics within that phase.
* The order of the species within the vector is based on
* the order of listed ThermoPhase objects in the class, and the
* order of the species within each ThermoPhase class.
*/
vector_fp m_conc;
vector_fp m_mu0;
vector_fp m_phi;
vector_fp m_pot;
vector_fp m_rwork;
vector_fp m_E;
vector_fp m_beta;
vector_int m_ctrxn;
private:
int reactionNumber(){ return m_ii;}
void addElementaryReaction(const ReactionData& r);
void addGlobalReaction(const ReactionData& r);
void installReagents(const ReactionData& r);
void installGroups(int irxn, const std::vector<grouplist_t>& r,
const std::vector<grouplist_t>& p);
void updateKc();
void registerReaction(int rxnNumber, int type, int loc) {
m_index[rxnNumber] = std::pair<int, int>(type, loc);
}
void applyButlerVolmerCorrection(doublereal* kf);
bool m_finalized;
bool m_has_electrochem_rxns;
};
}

72
Cantera/src/kinetics/InterfaceKinetics.cpp
View file

@ -13,6 +13,7 @@
#endif
#include "InterfaceKinetics.h"
#include "EdgeKinetics.h"
#include "SurfPhase.h"
#include "ReactionData.h"
@ -44,7 +45,8 @@ namespace Cantera {
m_surf(0),
m_integrator(0),
m_finalized(false),
m_has_coverage_dependence(false)
m_has_coverage_dependence(false),
m_has_electrochem_rxns(false)
{
if (thermo != 0) addPhase(*thermo);
m_kdata = new InterfaceKineticsData;
@ -77,7 +79,8 @@ namespace Cantera {
if (T != m_kdata->m_temp || m_redo_rates) {
m_kdata->m_logtemp = log(T);
m_rates.update(T, m_kdata->m_logtemp, DATA_PTR(m_kdata->m_rfn));
applyButlerVolmerCorrection(DATA_PTR(m_kdata->m_rfn));
if (m_has_electrochem_rxns)
applyButlerVolmerCorrection(DATA_PTR(m_kdata->m_rfn));
m_kdata->m_temp = T;
updateKc();
m_kdata->m_ROP_ok = false;
@ -153,15 +156,14 @@ namespace Cantera {
}
// compute Delta mu^0 for all reversible reactions
//m_reactantStoich.decrementReactions(m_mu0.begin(), m_rkc.begin());
//m_revProductStoich.incrementReactions(m_mu0.begin(), m_rkc.begin());
m_rxnstoich.getRevReactionDelta(m_ii, DATA_PTR(m_mu0),
DATA_PTR(m_rkc));
for (i = 0; i < m_nrev; i++) {
irxn = m_revindex[i];
if (irxn < 0 || irxn >= nReactions()) {
throw CanteraError("InterfaceKinetics","illegal value: irxn = "+int2str(irxn));
throw CanteraError("InterfaceKinetics",
"illegal value: irxn = "+int2str(irxn));
}
m_rkc[irxn] = exp(m_rkc[irxn]*rrt);
}
@ -238,9 +240,6 @@ namespace Cantera {
fill(kc, kc + m_ii, 0.0);
//m_reactantStoich.decrementReactions(m_mu0.begin(), kc);
//m_revProductStoich.incrementReactions(m_mu0.begin(), kc);
//m_irrevProductStoich.incrementReactions(m_mu0.begin(), kc);
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_mu0), kc);
for (i = 0; i < m_ii; i++) {
@ -328,23 +327,18 @@ namespace Cantera {
// activation energy below zero.
doublereal ea, eamod;
for (i = 0; i < m_ii; i++) {
eamod = 0.5*m_rwork[i];
int nct = m_beta.size();
int irxn;
for (i = 0; i < nct; i++) {
irxn = m_ctrxn[i];
eamod = m_beta[i]*m_rwork[irxn];
if (eamod != 0.0 && m_E[i] != 0.0) {
ea = GasConstant * m_E[i];
if (eamod + ea < 0.0) {
writelog("Warning: act energy mod too large!\n");
eamod = -ea;
writelog("warning: modified E < 0.\n");
}
kf[i] *= exp(-eamod*rrt);
// if (kf[i] == 0.0) {
// for (n = 0; n < np; n++) {
// cout << "phi " << n << " " << thermo(n).electricPotential() << " " << m_phi[n] << endl;
// }
// cout << "Zero rate coeff." << endl;
// cout << "eamod = " << eamod << " " << eamod*rrt << endl;
// cout << eamod/Faraday << endl;
// }
kf[irxn] *= exp(-eamod*rrt);
}
}
}
@ -380,10 +374,10 @@ namespace Cantera {
getFwdRateConstants(krev);
if (doIrreversible) {
doublereal *tmpKc = DATA_PTR(m_kdata->m_ropnet);
getEquilibriumConstants(tmpKc);
for (int i = 0; i < m_ii; i++) {
krev[i] /= tmpKc[i];
}
getEquilibriumConstants(tmpKc);
for (int i = 0; i < m_ii; i++) {
krev[i] /= tmpKc[i];
}
}
else {
const vector_fp& rkc = m_kdata->m_rkcn;
@ -661,18 +655,29 @@ namespace Cantera {
void InterfaceKinetics::
addElementaryReaction(const ReactionData& r) {
int iloc;
// install rate coeff calculator
vector_fp rp = r.rateCoeffParameters;
int ncov = r.cov.size();
if (ncov > 3) {
m_has_coverage_dependence = true;
}
for (int m = 0; m < ncov; m++) rp.push_back(r.cov[m]);
iloc = m_rates.install( reactionNumber(),
r.rateCoeffType, rp.size(),
DATA_PTR(rp) );
// store activation energy
m_E.push_back(r.rateCoeffParameters[2]);
if (r.beta > 0.0) {
m_has_electrochem_rxns = true;
m_beta.push_back(r.beta);
m_ctrxn.push_back(reactionNumber());
}
// add constant term to rate coeff value vector
m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]);
registerReaction( reactionNumber(), ELEMENTARY_RXN, iloc);
@ -757,7 +762,6 @@ namespace Cantera {
* of reactants for the rnum'th reaction
*/
m_reactants.push_back(rk);
vector_int pk;
int np = r.products.size();
for (n = 0; n < np; n++) {
@ -835,10 +839,14 @@ namespace Cantera {
*/
void InterfaceKinetics::finalize() {
m_rwork.resize(nReactions());
int ks = surfacePhaseIndex();
int ks = reactionPhaseIndex();
if (ks < 0) throw CanteraError("InterfaceKinetics::finalize",
"no surface phase is present.");
m_surf = (SurfPhase*)&thermo(ks);
if (m_surf->nDim() != 2)
throw CanteraError("InterfaceKinetics::finalize",
"expected interface dimension = 2, but got dimension = "
+int2str(m_surf->nDim()));
m_finalized = true;
}
@ -860,6 +868,18 @@ namespace Cantera {
m_integrator = 0;
}
void EdgeKinetics::finalize() {
m_rwork.resize(nReactions());
int ks = reactionPhaseIndex();
if (ks < 0) throw CanteraError("EdgeKinetics::finalize",
"no edge phase is present.");
m_surf = (SurfPhase*)&thermo(ks);
if (m_surf->nDim() != 1)
throw CanteraError("EdgeKinetics::finalize",
"expected interface dimension = 1, but got dimension = "
+int2str(m_surf->nDim()));
m_finalized = true;
}
}

38
Cantera/src/kinetics/InterfaceKinetics.h
View file

@ -88,10 +88,22 @@ namespace Cantera {
virtual int ID() { return cInterfaceKinetics; }
virtual int type() { return cInterfaceKinetics; }
///
/// @name Reaction Rates Of Progress
///
//@{
/**
* Set the electric potential in the nth phase
*
* @param n phase Index in this kinetics object.
* @param V Electric potential (volts)
*/
void setElectricPotential(int n, doublereal V) {
thermo(n).setElectricPotential(V);
m_redo_rates = true;
}
///
/// @name Reaction Rates Of Progress
///
//@{
//! Return the forward rates of progress for each reaction
/*!
@ -410,7 +422,6 @@ namespace Cantera {
* product stoichiometric coefficient for the species being the value.
*/
mutable std::vector<std::map<int, doublereal> > m_prxn;
//! String expression for each rxn
/*!
* Vector of strings of length m_ii, the number of
@ -489,15 +500,15 @@ namespace Cantera {
//! Pointer to the surface solver
ImplicitSurfChem* m_integrator;
public:
vector_fp m_beta;
vector_int m_ctrxn;
int reactionNumber(){ return m_ii;}
protected:
void addElementaryReaction(const ReactionData& r);
void addGlobalReaction(const ReactionData& r);
void installReagents(const ReactionData& r);
private:
void updateKc();
//! Write values into m_index
@ -512,10 +523,13 @@ namespace Cantera {
void applyButlerVolmerCorrection(doublereal* kf);
protected:
//! boolean indicating whether mechanism has been finalized
bool m_finalized;
bool m_has_coverage_dependence;
//! boolean indicating whether mechanism has been finalized
bool m_finalized;
bool m_has_coverage_dependence;
bool m_has_electrochem_rxns;
private:
};
}

2
Cantera/src/kinetics/Makefile.in
View file

@ -48,7 +48,7 @@ endif
# heterogeneous kinetics
ifeq ($(do_heterokin),1)
HETEROKIN_OBJ=InterfaceKinetics.o ImplicitSurfChem.o EdgeKinetics.o
HETEROKIN_OBJ=InterfaceKinetics.o ImplicitSurfChem.o
HETEROKIN_H =InterfaceKinetics.h ImplicitSurfChem.h EdgeKinetics.h
HETEROKIN = $(HETEROKIN_OBJ)
endif

2
Cantera/src/thermo/CMakeLists.txt
View file

@ -5,6 +5,6 @@ SET (CTTHERMO_SRCS State.cpp Elements.cpp Constituents.cpp Phase.cpp
ThermoFactory.cpp phasereport.cpp StoichSubstance.cpp
PureFluidPhase.cpp LatticeSolidPhase.cpp LatticePhase.cpp)
INCLUDE_DIRECTORIES (../base . )
INCLUDE_DIRECTORIES (../base)
ADD_LIBRARY(ctthermo ${CTTHERMO_SRCS})

1
cmake/fortran.cmake
View file

@ -3,6 +3,7 @@
#if (NOT BUILD_WITH_F2C)
#### Fortran 90
message("Fortran cmake")
if (BUILD_F90_INTERFACE)
if (F90 STREQUAL "default")

4
config.cmake
View file

@ -58,9 +58,9 @@
# build process
# default try to do a full installation, but fall back to a minimal
# one in case of errors
OPTION(BUILD_PYTHON_PACKAGE "Build the Python Package?")
OPTION(CANTERA_BUILD_PYTHON_PACKAGE "Build the Python Package?" ON)
SET(PYTHON_PACKAGE "default" CACHE STRING "full or minimal")
SET(CANTERA_PYTHON_PACKAGE_TYPE "default" CACHE LIST "full or minimal")
# Cantera needs to know where to find the Python interpreter. If
# PYTHON_CMD is set to "default", then cmake will look for the Python