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

Interfacekinetics rewrite:

Browse source 
Decided to break out an ElectrodeKinetics object for now. We'll see how it goes.
This commit is contained in:
Harry Moffat 2014-08-21 00:49:41 +00:00
parent 0d9b7d7868
commit 69c09709d1
15 changed files with 1259 additions and 72 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

87
include/cantera/kinetics/ElectrodeKinetics.h Normal file
View file

@ -0,0 +1,87 @@
/**
* @file ElectrodeKinetics.h
*
* @ingroup chemkinetics
*/
/*
* Copywrite (2005) Sandia Corporation. Under the terms of
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
* U.S. Government retains certain rights in this software.
*/
#ifndef CT_ELECTRODEKINETICS_H
#define CT_ELECTRODEKINETICS_H
#include "InterfaceKinetics.h"
namespace Cantera
{
//! A kinetics manager for heterogeneous reaction mechanisms. The
//! reactions are assumed to occur at a 2D interface between two 3D phases.
/*!
* This class is a slight addition to the InterfaceKinetics class, adding
* several concepts. First we explicity identify the electrode and solution
* phases. We will also assume that there is an electron phase.
*
* @ingroup chemkinetics
*/
class ElectrodeKinetics : public InterfaceKinetics
{
public:
//! Constructor
/*!
* @param thermo The optional parameter may be used to initialize
* the object with one ThermoPhase object.
* HKM Note -> Since the interface kinetics
* object will probably require multiple thermophase
* objects, this is probably not a good idea
* to have this parameter.
*/
ElectrodeKinetics(thermo_t* thermo = 0);
/// Destructor.
virtual ~ElectrodeKinetics();
//! Copy Constructor for the %Kinetics object.
ElectrodeKinetics(const ElectrodeKinetics& right);
//! Assignment operator
ElectrodeKinetics& operator=(const ElectrodeKinetics& right);
virtual Kinetics* duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const;
virtual int type() const;
//void addGlobalReaction(ReactionData& r);
protected:
//! index of the metal phase in the list of phases for this surface
size_t metalPhaseRS_;
size_t electronPhaseRS_;
//! Index of the solution phase in the list of phases for this surface
size_t solnPhaseRS_;
//! Index of the electrons species in the list of species for this surface kinetics, if none set it to -1
size_t kElectronRS_;
};
}
#endif

120
include/cantera/kinetics/ExtraGlobalRxn.h Normal file
View file

@ -0,0 +1,120 @@
/**
* @file ReactingVolDomain.h
*
*/
/*
* Copywrite (2005) Sandia Corporation. Under the terms of
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
* U.S. Government retains certain rights in this software.
*/
#ifndef EXTRAGLOBALRXN_H
#define EXTRAGLOBALRXN_H
#include "cantera/kinetics/InterfaceKinetics.h"
#include <string>
#include <vector>
namespace Cantera
{
//! Class describing an extra global reaction, which is defined as
//! a linear combination of actuals reactions, global or mass-action, creating a global stoichiometric result
/*!
* This is useful for defining thermodynamics of global processes that occur
* on a surface or in a homogeneous phase.
*
* The class is set up via the function setupElemRxnVector(RxnVector, specialSpecies) which defines
* the vector of stoichiometric coefficients representing the base reaction to combine in order to
* achieve the global result that's to be calculated. specialSpecies is the index of the species
* within the kinetics object that is used to identify the global reaction. Rates of progress
* are defined in terms of the production rate of the special species.
*
*/
class ExtraGlobalRxn
{
public:
//! Constructor takes a default kinetics pointer
/*!
* @param[in] k_ptr Pointer to a Kinetics class that will be used as the basis
* for constructing this class.
*/
ExtraGlobalRxn(Kinetics* k_ptr);
//! Destructor
virtual ~ExtraGlobalRxn();
void setupElemRxnVector(double* RxnVector,
int specialSpecies = -1);
std::string reactionString();
double deltaSpecValue(double* speciesVectorProperty);
std::vector<int>& reactants();
std::vector<int>& products();
bool isReversible();
double ROPValue(double* ROPKinVector);
double FwdROPValue(double* FwdROPElemKinVector, double* RevROPElemKinVector);
double RevROPValue(double* FwdROPElemKinVector, double* RevROPElemKinVector);
double reactantStoichCoeff(int kKin);
double productStoichCoeff(int kKin);
bool m_ThisIsASurfaceRxn;
double deltaRxnVecValue(double* rxnVectorProperty);
//! This kinetics operator is associated with just one
//! homogeneous phase, associated with tpList[0] phase
/*!
* Kinetics object pointer
*/
Cantera::Kinetics* m_kinetics;
//! This kinetics operator is associated with multiple
//! homogeneous and surface phases.
/*!
* This object owns the Kinetics object
*/
Cantera::InterfaceKinetics* m_InterfaceKinetics;
int m_nKinSpecies;
//! Number of reactants in the global reaction
int m_nReactants;
//! Vector of reactants that make up the global reaction
/*!
* This is a list of reactants using the kinetic species index
*/
std::vector<int> m_Reactants;
//! Vector of reactant stoichiometries that make up the global reaction
/*!
* This is a list of reactant stoichiometries. The species index is given in
* the member m_Reactants using the kinetic species index.
*/
std::vector<doublereal> m_ReactantStoich;
int m_nProducts;
std::vector<int> m_Products;
std::vector<doublereal> m_ProductStoich;
int m_nNetSpecies;
std::vector<int> m_NetSpecies;
std::vector<doublereal> m_netStoich;
int m_nRxns;
std::vector<doublereal> m_ElemRxnVector;
int m_SpecialSpecies;
bool m_SpecialSpeciesProduct;
int m_SS_index;
int iphaseKin;
bool m_ok;
bool m_reversible;
};
}
#endif

43
include/cantera/kinetics/InterfaceKinetics.h
View file

@ -105,6 +105,16 @@ public:
*/
void updateExchangeCurrentQuantities();
//! Return the vector of values for the reaction gibbs free energy change.
/*!
* (virtual from Kinetics.h)
* These values depend upon the concentration of the solution.
*
* units = J kmol-1
*
* @param deltaG Output vector of deltaG's for reactions Length: m_ii.
* If 0, this updates the internally stored values only.
*/
virtual void getDeltaGibbs(doublereal* deltaG);
virtual void getDeltaElectrochemPotentials(doublereal* deltaM);
@ -306,11 +316,13 @@ public:
/*!
* @param rxnNumber reaction number
* @param type reaction type
* @param loc location ??
* @param loc location location in the reaction rate coefficient calculator
*
* Right now we only use one reaction rate coefficient calculated named ELEMENTARY_RXN
* Therefore, this is not used within the code)
* (type, loc) is stored as a std::pair
*/
void registerReaction(size_t rxnNumber, int type, size_t loc) {
m_index[rxnNumber] = std::pair<int, size_t>(type, loc);
}
void registerReaction(size_t rxnNumber, int type, size_t loc);
//! Apply modifications for the fowward reaction rate for interfacial charge transfer reactions
/*!
@ -396,6 +408,8 @@ protected:
//! Temporary work vector of length m_kk
vector_fp m_grt;
//! List of reactions numbers which are reversible reactions
/*!
* This is a vector of reaction numbers. Each reaction in the list is
@ -515,6 +529,14 @@ protected:
*/
vector_fp m_mu0;
//! Vector of chemical potentials for all species
/*!
* This vector contains a vector of chemical potentials for all of the species in the kinetics object
*
* Length = m_kk. Units = J/kmol.
*/
vector_fp m_mu;
//! Vector of standard state electrochemical potentials modified by
//! a standard concentration term.
/*!
@ -636,6 +658,19 @@ protected:
* units are Joule kmol-1
*/
vector_fp m_deltaG0;
//! Vector of deltaG[] of reaction, the delta gibbs free energies for each reaction
/*!
* Length is the number of reactions
* units are Joule kmol-1
*/
vector_fp m_deltaG;
//! Vector of the products of the standard concentrations of the reactants
/*!
* Units vary wrt what the units of the standard concentrations are
* Length = number of reactions.
*/
vector_fp m_ProdStanConcReac;
doublereal m_logp0;

1
include/cantera/kinetics/Kinetics.h
View file

@ -465,6 +465,7 @@ public:
//! Return the vector of values for the reaction gibbs free energy change.
/*!
* (virtual from Kinetics.h)
* These values depend upon the concentration of the solution.
*
* units = J kmol-1

3
include/cantera/kinetics/ReactionData.h
View file

@ -13,6 +13,9 @@ namespace Cantera
//! Intermediate class which stores data about a reaction and its rate
//! parameterization before adding the reaction to a Kinetics object.
/*!
* All data in this class is public.
*/
class ReactionData
{
public:

15
include/cantera/kinetics/ReactionStoichMgr.h
View file

@ -66,8 +66,9 @@ public:
ReactionStoichMgr& operator=(const ReactionStoichMgr& right);
/**
* Add a reaction with mass-action kinetics. Vectors
//! Add a reaction with mass-action kinetics.
/*!Vectors
* 'reactants' and 'products' contain the integer species
* indices of the reactants and products, respectively. Note
* that if more than one molecule of a given species is
@ -86,11 +87,11 @@ public:
* - reactants: (2, 2) [ note repeated index ]
* - products: (1)
*
* @param rxn Reaction number. This number will be used as the index
* into the rate of progress vector in the methods below.
* @param reactants vector of integer reactant indices
* @param products vector of integer product indices
* @param reversible true if the reaction is reversible, false otherwise
* @param rxn Reaction number. This number will be used as the index
* into the rate of progress vector in the methods below.
* @param reactants Vector of integer reactant indices
* @param products Vector of integer product indices
* @param reversible True if the reaction is reversible, false otherwise
*/
virtual void add(size_t rxn, const std::vector<size_t>& reactants,
const std::vector<size_t>& products, bool reversible);

117
include/cantera/kinetics/RxnMolChange.h Normal file
View file

@ -0,0 +1,117 @@
/**
* @file RxnMolChange.cpp
*
*/
/*
* Copywrite (2005) Sandia Corporation. Under the terms of
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
* U.S. Government retains certain rights in this software.
*/
#ifndef RXNMOLCHANGE_H
#define RXNMOLCHANGE_H
#include <vector>
namespace Cantera
{
class ExtraGlobalRxn;
class Kinetics;
//! Class that includes some bookeeping entries for a reaction or a global reaction defined on a surface
/*!
* Note that all indexes refer to a specific interfacial or homogeneous kinetics object. It does not
* refer to the Phase list indexes.
*/
class RxnMolChange
{
public:
//! Main constructor for the class
/*!
* @param kinPtr Pointer to the kinetics base class
* @param irxn Specific reaction index.
*/
RxnMolChange(Cantera::Kinetics* kinPtr, int irxn);
//! Destructor
~RxnMolChange();
//! Constructor for the object if the object refers to a global reaction
/*!
* @param kinPtr Pointer to the kinetics base class
* @param egr Specific reaction index.
*/
RxnMolChange(Cantera::Kinetics* kinPtr, Cantera::ExtraGlobalRxn* egr);
//! Vector of mole changes for each phase in the Kinetics object due to the current reaction
/*!
* This is the sum of the product stoichiometric coefficient minum the reactant stoichioemtric coefficient
* for all the species in a phase.
* The index is over all the phases listed in the Kinetics object.
*/
std::vector<double> m_phaseMoleChange;
std::vector<double> m_phaseReactantMoles;
std::vector<double> m_phaseProductMoles;
//! Vector of mass changes for each phase in the Kinetics object due to the current reaction
/*!
* This is the sum of the product stoichiometric coefficient minum the reactant stoichioemtric coefficient
* index multiplied by the molecular weight for all species in a phase.
* The index is over all of the phases listed in the Kinetics object.
*/
std::vector<double> m_phaseMassChange;
//! Vector of mass changes for each phase in the Kinetics object due to the current reaction
/*!
* This is the sum of the product stoichiometric coefficient minum the reactant stoichioemtric coefficient
* index multiplied by the charg for all species in a phase.
* The index is over all of the phases listed in the Kinetics object.
*/
std::vector<double> m_phaseChargeChange;
//! Vector of phase types in the reaction
/*!
* Collection of eosTypes for all phases in the kinetics object
* The index is over all of the phases listed in the Kinetics object.
*/
std::vector<int> m_phaseTypes;
//! Vector of phase dimensions for the reaction
/*!
* Collection of nDims for all phases in the kinetics object
* The index is over all of the phases listed in the Kinetics object.
*/
std::vector<int> m_phaseDims;
//! Number of phases in the kientics object
int m_nPhases;
//! Shallow pointer pointing to the kinetics object
Cantera::Kinetics* m_kinBase;
//! Reaction number within the kinetics object
/*!
* If this is neg 1, then this reaction refers to a global reaction
* specified by the m_egr pointer.
*/
int m_iRxn;
//! Maximum change in charge of any phase due to this reaction
double m_ChargeTransferInRxn;
//! Electrochemical beta parameter for the reaction
double m_beta;
//! Pointer to the specification of the global reaction
/*!
* This is 0, if the class refers to a single reaction in the kinetics object
*/
Cantera::ExtraGlobalRxn* m_egr;
};
}
#endif

52
include/cantera/kinetics/StoichManager.h
View file

@ -8,6 +8,7 @@
#define CT_STOICH_MGR_H
#include "cantera/base/stringUtils.h"
#include "cantera/base/ctexceptions.h"
namespace Cantera
{
@ -452,7 +453,8 @@ public:
m_rxn(right.m_rxn),
m_ic(right.m_ic),
m_order(right.m_order),
m_stoich(right.m_stoich) {
m_stoich(right.m_stoich)
{
}
C_AnyN& operator=(const C_AnyN& right) {
@ -560,27 +562,44 @@ public:
}
private:
//! Length of the m_ic vector
/*!
* This is the number of species which have non-zero entries in either the
* reaction order matrix or the stoichiometric order matrix for this reaction.
* This is the number of species which participate in the reaction order
* and stoichiometric coefficient vectors for the reactant or product description
* of the reaction.
*/
size_t m_n;
//! ID of the reaction corresponding to this stoichiometric manager
/*!
* This is used within the interface to select the
* This is used within the interface to select the array position to read and write to
* Normally this is associated with the reaction number in an array of quantities indexed
* by the reaction number, e.g., ROP[irxn].
*/
size_t m_rxn;
//! Vector of species which are involved with this stoichiometric manager calculations
/*!
* This is an integer list of species which have non-zero entries in either the
* This is an integer list of species which participate in either the
* reaction order matrix or the stoichiometric order matrix for this reaction, m_rxn.
* It's used as the index into the arrays m_order[] and m_stoich[].
*/
std::vector<size_t> m_ic;
//! Reaction orders for the reaction
/*!
* This is either for the reactants or products.
* Length = m_n
* Species number, m_ic[n], has a reaction order of m_order[n].
*/
vector_fp m_order;
//! Stoichiometric coefficients for the reaction, reactant or product side.
/*!
* This is either for the reactants or products.
* Length = m_n
* Species number m_ic[m], has a stoichiometric coefficient of m_stoich[n].
*/
vector_fp m_stoich;
};
@ -689,7 +708,7 @@ inline static void _writeMultiply(InputIter begin, InputIter end,
* be the number of molecules on the product or reactant side of
* reaction number i.
* \f[
* r_i = \sum_m^{M_i} s_{k_{m,i}}
* r_i = \sum_m^{M_i} s_{k_{m,i}}
* \f]
* To understand the operations performed by this class, let
* \f$ N_{k,i}\f$ denote the stoichiometric coefficient of species k on
@ -702,7 +721,7 @@ inline static void _writeMultiply(InputIter begin, InputIter end,
* - \f$ S = S + N R\f$ (incrementSpecies)
* - \f$ S = S - N R\f$ (decrementSpecies)
* - \f$ R = R + N^T S \f$ (incrementReaction)
* - \f$ R = R - N^T S \f$ (deccrementReaction)
* - \f$ R = R - N^T S \f$ (decrementReaction)
*
* The actual implementation, however, does not compute these
* quantities by matrix multiplication. A faster algorithm is used
@ -794,6 +813,12 @@ public:
void add(size_t rxn, const std::vector<size_t>& k, const vector_fp& order,
const vector_fp& stoich) {
m_n[rxn] = k.size();
if (order.size() != k.size()) {
throw CanteraError("StoichManagerN::add()", "size of order and species arrays differ");
}
if (stoich.size() != k.size()) {
throw CanteraError("StoichManagerN::add()", "size of stoich and species arrays differ");
}
bool frac = false;
for (size_t n = 0; n < stoich.size(); n++) {
if (stoich[n] != 1.0 || order[n] != 1.0) {
@ -900,11 +925,16 @@ private:
std::vector<C2> m_c2_list;
std::vector<C3> m_c3_list;
std::vector<C_AnyN> m_cn_list;
/**
* Map with the Reaction Number as key and the Number of species
* as the value.
//! Map with the Reaction Number as key and the Number of species
//! as the value.
/*!
* An example of this is given below:
*
* m_n[irxn] = nSpecies
*/
std::map<size_t, size_t> m_n;
/**
* Map with the Reaction Number as key and the placement in the
* vector of reactions list( i.e., m_c1_list[]) as key

74
src/kinetics/ElectrodeKinetics.cpp Normal file
View file

@ -0,0 +1,74 @@
/**
* @file ElectrodeKinetics.cpp
*/
#include "cantera/kinetics/ElectrodeKinetics.h"
using namespace std;
namespace Cantera
{
//============================================================================================================================
ElectrodeKinetics::ElectrodeKinetics(thermo_t* thermo) :
InterfaceKinetics(thermo),
metalPhaseRS_(npos),
electronPhaseRS_(npos),
solnPhaseRS_(npos),
kElectronRS_(npos)
{
}
//============================================================================================================================
ElectrodeKinetics::~ElectrodeKinetics()
{
}
//============================================================================================================================
ElectrodeKinetics::ElectrodeKinetics(const ElectrodeKinetics& right) :
InterfaceKinetics()
{
/*
* Call the assignment operator
*/
operator=(right);
}
//============================================================================================================================
ElectrodeKinetics& ElectrodeKinetics::operator=(const ElectrodeKinetics& right)
{
/*
* Check for self assignment.
*/
if (this == &right) {
return *this;
}
InterfaceKinetics::operator=(right);
metalPhaseRS_ = right.metalPhaseRS_;
electronPhaseRS_ = right.electronPhaseRS_;
solnPhaseRS_ = right.solnPhaseRS_;
kElectronRS_ = right.kElectronRS_;
return *this;
}
//============================================================================================================================
int ElectrodeKinetics::type() const
{
return cInterfaceKinetics;
}
//============================================================================================================================
Kinetics* ElectrodeKinetics::duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const
{
ElectrodeKinetics* iK = new ElectrodeKinetics(*this);
iK->assignShallowPointers(tpVector);
return iK;
}
//============================================================================================================================
//==================================================================================================================
}

394
src/kinetics/ExtraGlobalRxn.cpp Normal file
View file

@ -0,0 +1,394 @@
/**
* @file example2.cpp
*
*/
/*
* $Id: ExtraGlobalRxn.cpp 571 2013-03-26 16:44:21Z hkmoffa $
*
*/
/*
* Copywrite (2005) Sandia Corporation. Under the terms of
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
* U.S. Government retains certain rights in this software.
*/
// Example 2
//
// Read a mechanism, and print to the standard output stream a
// well-formatted Chemkin ELEMENT section.
//
#include "cantera/kinetics/ExtraGlobalRxn.h"
#include "cantera/numerics/DenseMatrix.h"
// Kinetics includes
#include "cantera/kinetics.h"
#include "cantera/kinetics/InterfaceKinetics.h"
#include "cantera/thermo/SurfPhase.h"
#include "cantera/kinetics/KineticsFactory.h"
#include <iostream>
#include <new>
#include <string>
#include <vector>
#include <typeinfo>
using namespace std;
using namespace Cantera;
namespace Cantera {
//============================================================================================================
static void erase_vd(std::vector<doublereal>& m_vec, int index)
{
std::vector<double>::iterator ipos;
ipos = m_vec.begin();
ipos += index;
m_vec.erase(ipos);
}
//============================================================================================================
static void erase_vi(std::vector<int>& m_vec, int index)
{
std::vector<int>::iterator ipos;
ipos = m_vec.begin();
ipos += index;
m_vec.erase(ipos);
}
//============================================================================================================
//! add the species into the list of products or reactants
/*!
* Note this function gets called for both the product and reactant sides. However, it's only
* called for one side or another.
*
* @param kkinspec kinetic species index of the product
* @param
*/
static void addV(int kkinspec, double ps, std::vector<int>& m_Products,
std::vector<doublereal>& m_ProductStoich)
{
int nsize = m_Products.size();
for (int i = 0; i < nsize; i++) {
if (m_Products[i] == kkinspec) {
m_ProductStoich[i] += ps;
return;
}
}
m_Products.push_back(kkinspec);
m_ProductStoich.push_back(ps);
}
//============================================================================================================
ExtraGlobalRxn::ExtraGlobalRxn(Kinetics* k_ptr) :
m_ThisIsASurfaceRxn(false),
m_kinetics(k_ptr),
m_InterfaceKinetics(0),
m_nKinSpecies(0),
m_nReactants(0),
m_nProducts(0),
m_nNetSpecies(0),
m_nRxns(0),
m_SpecialSpecies(-1),
m_SpecialSpeciesProduct(true),
iphaseKin(0),
m_ok(false),
m_reversible(true)
{
m_InterfaceKinetics = dynamic_cast<InterfaceKinetics*>(k_ptr);
if (m_InterfaceKinetics) {
m_ThisIsASurfaceRxn = true;
}
m_nRxns = m_kinetics->nReactions();
m_ElemRxnVector.resize(m_nRxns,0.0);
m_nKinSpecies = m_kinetics->nTotalSpecies();
}
//============================================================================================================
ExtraGlobalRxn::~ExtraGlobalRxn()
{
}
//============================================================================================================
void ExtraGlobalRxn::setupElemRxnVector(double* RxnVector,
int specialSpecies)
{
int i;
int kkinspec;
for (size_t i = 0; i < (size_t) m_nRxns; i++) {
m_ElemRxnVector[i] = RxnVector[i];
}
for (size_t i = 0; i < (size_t) m_nRxns; i++) {
if (RxnVector[i] > 0.0) {
for (kkinspec = 0; kkinspec < m_nKinSpecies; kkinspec++) {
double ps = m_kinetics->productStoichCoeff(kkinspec, i);
if (ps > 0.0) {
addV(kkinspec, RxnVector[i]* ps, m_Products, m_ProductStoich);
addV(kkinspec, RxnVector[i]* ps, m_NetSpecies, m_netStoich);
}
double rs = m_kinetics->reactantStoichCoeff(kkinspec, i);
if (rs > 0.0) {
addV(kkinspec, RxnVector[i] * rs, m_Reactants, m_ReactantStoich);
addV(kkinspec, -RxnVector[i] * rs, m_NetSpecies, m_netStoich);
}
}
} else if (RxnVector[i] < 0.0) {
for (kkinspec = 0; kkinspec < m_nKinSpecies; kkinspec++) {
double ps = m_kinetics->productStoichCoeff(kkinspec, i);
if (ps > 0.0) {
addV(kkinspec,- RxnVector[i]* ps, m_Reactants, m_ReactantStoich);
addV(kkinspec, RxnVector[i]* ps, m_NetSpecies, m_netStoich);
}
double rs = m_kinetics->reactantStoichCoeff(kkinspec, i);
if (rs > 0.0) {
addV(kkinspec, -RxnVector[i] * rs, m_Products, m_ProductStoich);
addV(kkinspec, -RxnVector[i] * rs, m_NetSpecies, m_netStoich);
}
}
}
}
Recheck:
for (i = 0; i < static_cast<int>(m_Products.size()); i++) {
if (m_ProductStoich[i] == 0.0) {
erase_vi(m_Products, i);
erase_vd(m_ProductStoich, i);
goto Recheck ;
}
}
for (i = 0; i < static_cast<int>(m_Reactants.size()); i++) {
if (m_ReactantStoich[i] == 0.0) {
erase_vi(m_Reactants, i);
erase_vd(m_ReactantStoich, i);
goto Recheck ;
}
}
for (i = 0; i < static_cast<int>(m_NetSpecies.size()); i++) {
if (m_netStoich[i] == 0.0) {
erase_vi(m_NetSpecies, i);
erase_vd(m_netStoich, i);
goto Recheck ;
}
}
for (i = 0; i < static_cast<int>(m_Products.size()); i++) {
int ik = m_Products[i];
for (int j = 0; j < static_cast<int>(m_Reactants.size()); j++) {
int jk = m_Reactants[j];
if (ik == jk) {
if (m_ProductStoich[i] == m_ReactantStoich[j]) {
erase_vi(m_Products, i);
erase_vd(m_ProductStoich, i);
erase_vi(m_Reactants, j);
erase_vd(m_ReactantStoich, j);
} else if (m_ProductStoich[i] > m_ReactantStoich[j]) {
m_ProductStoich[i] -= m_ReactantStoich[j];
erase_vi(m_Reactants, j);
erase_vd(m_ReactantStoich, j);
} else {
m_ReactantStoich[j] -= m_ProductStoich[i];
erase_vi(m_Products, i);
erase_vd(m_ProductStoich, i);
}
// We just screwed up the indexing -> restart.
goto Recheck ;
}
}
}
m_nProducts = m_Products.size();
m_nReactants = m_Reactants.size();
m_nNetSpecies = m_NetSpecies.size();
/*
* Section to assign the special species
*/
m_SpecialSpecies = specialSpecies;
if (specialSpecies == -1) {
m_SpecialSpecies = m_Products[0];
}
bool ifound = false;
for (i = 0; i < (int) m_NetSpecies.size(); i++) {
int ik = m_NetSpecies[i];
if (ik == m_SpecialSpecies) {
if (m_netStoich[i] > 0.0) {
m_SpecialSpeciesProduct = true;
} else {
m_SpecialSpeciesProduct = false;
}
m_SS_index = i;
ifound = true;
break;
}
}
if (!ifound) {
throw CanteraError(":setupElemRxnVector",
"Special species not a reactant or product: "
+ int2str(m_SpecialSpecies));
}
m_ok = true;
}
//============================================================================================================
std::string ExtraGlobalRxn::reactionString()
{
string rs;
int k, istoich;
for (k = 0; k < m_nReactants; k++) {
int kkinspecies = m_Reactants[k];
double stoich = m_ReactantStoich[k];
if (stoich != 1.0) {
istoich = (int) stoich;
if (fabs((double)istoich - stoich) < 0.00001) {
rs += int2str(istoich) + " ";
} else {
rs += fp2str(stoich) + " ";
}
}
string sName = m_kinetics->kineticsSpeciesName(kkinspecies);
rs += sName;
if (k < (m_nReactants-1)) {
rs += " + ";
}
}
rs += " = ";
for (k = 0; k < m_nProducts; k++) {
int kkinspecies = m_Products[k];
double stoich = m_ProductStoich[k];
if (stoich != 1.0) {
istoich = (int) stoich;
if (fabs((double)istoich - stoich) < 0.00001) {
rs += int2str(istoich) + " ";
} else {
rs += fp2str(stoich) + " ";
}
}
string sName = m_kinetics->kineticsSpeciesName(kkinspecies);
rs += sName;
if (k < (m_nProducts-1)) {
rs += " + ";
}
}
return rs;
}
//============================================================================================================
std::vector<int>& ExtraGlobalRxn::reactants()
{
return m_Reactants;
}
//============================================================================================================
std::vector<int>& ExtraGlobalRxn::products()
{
return m_Products;
}
//============================================================================================================
bool ExtraGlobalRxn::isReversible()
{
return m_reversible;
}
//============================================================================================================
double ExtraGlobalRxn::reactantStoichCoeff(int kKin)
{
for (int k = 0; k < m_nReactants; k++) {
int kkinspec = m_Reactants[k];
if (kkinspec == kKin) {
return m_ReactantStoich[k];
}
}
return 0.0;
}
//============================================================================================================
double ExtraGlobalRxn::productStoichCoeff(int kKin)
{
for (int k = 0; k < m_nProducts; k++) {
int kkinspec = m_Products[k];
if (kkinspec == kKin) {
return m_ProductStoich[k];
}
}
return 0.0;
}
//============================================================================================================
double ExtraGlobalRxn::deltaSpecValue(double* speciesVectorProperty)
{
int k;
double val = 0;
for (k = 0; k < m_nNetSpecies; k++) {
int kkinspec = m_NetSpecies[k];
val += speciesVectorProperty[kkinspec] * m_netStoich[k];
}
return val;
}
//============================================================================================================
double ExtraGlobalRxn::deltaRxnVecValue(double* rxnVectorProperty)
{
double val = 0;
for (int i = 0; i < m_nRxns; i++) {
val += m_ElemRxnVector[i] * rxnVectorProperty[i];
}
return val;
}
//============================================================================================================
double ExtraGlobalRxn::ROPValue(double* ROPElemKinVector)
{
double val = 0.0;
for (int i = 0; i < m_nRxns; i++) {
double kstoich = m_kinetics->productStoichCoeff(m_SpecialSpecies, i) - m_kinetics->reactantStoichCoeff(m_SpecialSpecies, i);
if (m_ElemRxnVector[i] > 0.0) {
val += ROPElemKinVector[i] * kstoich * m_ElemRxnVector[i];
} else {
val -= ROPElemKinVector[i] * kstoich * m_ElemRxnVector[i];
}
}
if (!m_SpecialSpeciesProduct) {
val = -val;
}
return val;
}
//============================================================================================================
double ExtraGlobalRxn::FwdROPValue(double* FwdROPElemKinVector,
double* RevROPElemKinVector)
{
double val = 0.0;
for (int i = 0; i < m_nRxns; i++) {
double kstoich = m_kinetics->productStoichCoeff(m_SpecialSpecies, i) - m_kinetics->reactantStoichCoeff(m_SpecialSpecies, i);
if (m_ElemRxnVector[i] > 0.0) {
val += FwdROPElemKinVector[i] * kstoich * m_ElemRxnVector[i];
}
if (m_ElemRxnVector[i] < 0.0) {
val += RevROPElemKinVector[i] * kstoich * m_ElemRxnVector[i];
}
}
if (!m_SpecialSpeciesProduct) {
val = -val;
}
return val;
}
//============================================================================================================
double ExtraGlobalRxn::RevROPValue(double* FwdROPElemKinVector,
double* RevROPElemKinVector)
{
double val = 0.0;
for (int i = 0; i < m_nRxns; i++) {
double kstoich = m_kinetics->productStoichCoeff(m_SpecialSpecies, i)- m_kinetics->reactantStoichCoeff(m_SpecialSpecies, i);
if (m_ElemRxnVector[i] > 0.0) {
val += RevROPElemKinVector[i] * kstoich * m_ElemRxnVector[i];
}
if (m_ElemRxnVector[i] < 0.0) {
val += FwdROPElemKinVector[i] * kstoich * m_ElemRxnVector[i];
}
}
if (!m_SpecialSpeciesProduct) {
val = -val;
}
return val;
}
//============================================================================================================
}

157
src/kinetics/InterfaceKinetics.cpp
View file

@ -31,6 +31,7 @@ InterfaceKinetics::InterfaceKinetics(thermo_t* thermo) :
m_ctrxn_ecdf(0),
m_StandardConc(0),
m_deltaG0(0),
m_deltaG(0),
m_ProdStanConcReac(0),
m_logp0(0.0),
m_logc0(0.0),
@ -71,6 +72,7 @@ InterfaceKinetics::InterfaceKinetics(const InterfaceKinetics& right) :
m_ctrxn_ecdf(0),
m_StandardConc(0),
m_deltaG0(0),
m_deltaG(0),
m_ProdStanConcReac(0),
m_logp0(0.0),
m_logc0(0.0),
@ -121,6 +123,7 @@ InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right)
m_conc = right.m_conc;
m_actConc = right.m_actConc;
m_mu0 = right.m_mu0;
m_mu = right.m_mu;
m_mu0_Kc = right.m_mu0_Kc;
m_phi = right.m_phi;
m_pot = right.m_pot;
@ -134,6 +137,7 @@ InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right)
m_ctrxn_ecdf = right.m_ctrxn_ecdf;
m_StandardConc = right.m_StandardConc;
m_deltaG0 = right.m_deltaG0;
m_deltaG = right.m_deltaG;
m_ProdStanConcReac = right.m_ProdStanConcReac;
m_logp0 = right.m_logp0;
m_logc0 = right.m_logc0;
@ -347,13 +351,13 @@ void InterfaceKinetics::getRevRatesOfProgress(doublereal* revROP)
updateROP();
std::copy(m_ropr.begin(), m_ropr.end(), revROP);
}
//===========================================================================================================
void InterfaceKinetics::getNetRatesOfProgress(doublereal* netROP)
{
updateROP();
std::copy(m_ropnet.begin(), m_ropnet.end(), netROP);
}
//===========================================================================================================
void InterfaceKinetics::getEquilibriumConstants(doublereal* kc)
{
updateMu0();
@ -367,8 +371,9 @@ void InterfaceKinetics::getEquilibriumConstants(doublereal* kc)
kc[i] = exp(-kc[i]*rrt);
}
}
/** values needed to convert from exchange current density to surface reaction rate.
//===========================================================================================================
/*
* values needed to convert from exchange current density to surface reaction rate.
*/
void InterfaceKinetics::updateExchangeCurrentQuantities()
{
@ -391,26 +396,28 @@ void InterfaceKinetics::updateExchangeCurrentQuantities()
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_mu0), DATA_PTR(m_deltaG0));
//
// Calculate the product of the standard concentrations of the reactants
//
for (size_t i = 0; i < m_ii; i++) {
m_ProdStanConcReac[i] = 1.0;
}
m_rxnstoich.multiplyReactants(DATA_PTR(m_StandardConc), DATA_PTR(m_ProdStanConcReac));
}
//===========================================================================================================
void InterfaceKinetics::getCreationRates(doublereal* cdot)
{
updateROP();
m_rxnstoich.getCreationRates(m_kk, &m_ropf[0], &m_ropr[0], cdot);
}
//===========================================================================================================
void InterfaceKinetics::getDestructionRates(doublereal* ddot)
{
updateROP();
m_rxnstoich.getDestructionRates(m_kk, &m_ropf[0], &m_ropr[0], ddot);
}
//===========================================================================================================
void InterfaceKinetics::getNetProductionRates(doublereal* net)
{
updateROP();
@ -566,34 +573,87 @@ void InterfaceKinetics::updateROP()
_update_rates_T();
// get updated activities (rates updated below)
_update_rates_C();
double TT = m_surf->temperature();
double rtdf = GasConstant * TT / Faraday;
if (m_ROP_ok) {
return;
}
// copy rate coefficients into ropf
//
// Copy the reaction rate coefficients, m_rfn, into m_ropf
//
copy(m_rfn.begin(), m_rfn.end(), m_ropf.begin());
// multiply by perturbation factor
//
// Multiply by the perturbation factor
//
multiply_each(m_ropf.begin(), m_ropf.end(), m_perturb.begin());
// copy the forward rates to the reverse rates
//
// Copy the forward rate constants to the reverse rate constants
//
copy(m_ropf.begin(), m_ropf.end(), m_ropr.begin());
// for reverse rates computed from thermochemistry, multiply
//
// 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 the actyivity concentration reaction orders to obtain
// the forward rates of progress.
//
m_rxnstoich.multiplyReactants(DATA_PTR(m_actConc), DATA_PTR(m_ropf));
//
// For reversible reactions, multiply ropr by the activity concentration products
//
m_rxnstoich.multiplyRevProducts(DATA_PTR(m_actConc), DATA_PTR(m_ropr));
//
// Fix up these calculations for cases where the above formalism doesn't hold
//
double OCV = 0.0;
for (size_t jrxn = 0; jrxn != m_ii; ++jrxn) {
int reactionType = reactionTypes_[jrxn];
if (reactionType == BUTLERVOLMER_RXN) {
//
// OK, the reaction rate constant contains the current density rate constant calculation
// the rxnstoich calculation contained the dependence of the current density on the activity concentrations
// We finish up with the ROP calculation
//
// Calculate the overpotential of the reaction
//
// double nStoichElectrons = - rmc->m_phaseChargeChange[metalPhaseRS_];
double nStoichElectrons=1;
//*nStoich = nStoichElectrons;
getDeltaGibbs(0);
if (nStoichElectrons != 0.0) {
OCV = m_deltaG[jrxn]/Faraday/ nStoichElectrons;
}
/*
double exp1 = nu * nStoich * beta / rtdf
double exp2 = -nu * nStoich * Faraday * (1.0 - beta) / (GasConstant * temp);
double val = io * (exp(exp1) - exp(exp2));
doublereal BVterm = exp(exp1 ) - exp(exp2);
m_ropnet[j] = m_ropf[j] * BVterm
m_ropf[j] =
//
m_ropr[j] = m_ropnet[j] - m_ropf[j];
*/
}
}
// for reversible reactions, multiply ropr by the activity concentration
// products
m_rxnstoich.multiplyRevProducts(DATA_PTR(m_actConc),
DATA_PTR(m_ropr));
for (size_t j = 0; j != m_ii; ++j) {
m_ropnet[j] = m_ropf[j] - m_ropr[j];
@ -670,13 +730,18 @@ void InterfaceKinetics::getDeltaGibbs(doublereal* deltaG)
* kinetics mechanism
*/
for (size_t n = 0; n < nPhases(); n++) {
thermo(n).getChemPotentials(DATA_PTR(m_grt) + m_start[n]);
m_thermo[n]->getChemPotentials(DATA_PTR(m_mu) + m_start[n]);
}
//
// Use the stoichiometric manager to find deltaG for each
// reaction.
//
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_mu), DATA_PTR(m_deltaG));
if (deltaG != 0 && (DATA_PTR(m_deltaG) != deltaG)) {
for (size_t j = 0; j < m_ii; ++j) {
deltaG[j] = m_deltaG[j];
}
}
/*
* Use the stoichiometric manager to find deltaG for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaG);
}
//==================================================================================================================
void InterfaceKinetics::getDeltaElectrochemPotentials(doublereal* deltaM)
@ -709,7 +774,7 @@ void InterfaceKinetics::getDeltaEnthalpy(doublereal* deltaH)
*/
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaH);
}
//===========================================================================================================
void InterfaceKinetics::getDeltaEntropy(doublereal* deltaS)
{
/*
@ -725,7 +790,7 @@ void InterfaceKinetics::getDeltaEntropy(doublereal* deltaS)
*/
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaS);
}
//===========================================================================================================
void InterfaceKinetics::getDeltaSSGibbs(doublereal* deltaGSS)
{
/*
@ -743,7 +808,7 @@ void InterfaceKinetics::getDeltaSSGibbs(doublereal* deltaGSS)
*/
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_mu0), deltaGSS);
}
//===========================================================================================================
void InterfaceKinetics::getDeltaSSEnthalpy(doublereal* deltaH)
{
/*
@ -755,7 +820,7 @@ void InterfaceKinetics::getDeltaSSEnthalpy(doublereal* deltaH)
for (size_t n = 0; n < nPhases(); n++) {
thermo(n).getEnthalpy_RT(DATA_PTR(m_grt) + m_start[n]);
}
doublereal RT = thermo().temperature() * GasConstant;
doublereal RT = thermo(0).temperature() * GasConstant;
for (size_t k = 0; k < m_kk; k++) {
m_grt[k] *= RT;
}
@ -765,7 +830,7 @@ void InterfaceKinetics::getDeltaSSEnthalpy(doublereal* deltaH)
*/
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaH);
}
//===========================================================================================================
void InterfaceKinetics::getDeltaSSEntropy(doublereal* deltaS)
{
/*
@ -859,6 +924,12 @@ void InterfaceKinetics::addElementaryReaction(ReactionData& rdata)
rp.push_back(rdata.cov[m]);
}
//
// Find out the reaction type
//
int reactionType = rdata.reactionType;
reactionTypes_.push_back(reactionType);
/*
* Temporarily change the reaction rate coefficient type to surface arrhenius.
* This is what is expected. We'll handle exchange current types below by hand.
@ -919,7 +990,8 @@ void InterfaceKinetics::addGlobalReaction(ReactionData& rdata)
//
// Find out the reaction type
//
int reactionType = rdata.reactionType;
int reactionType = rdata.reactionType;
reactionTypes_.push_back(reactionType);
/*
* Temporarily change the reaction rate coefficient type to surface arrhenius.
@ -945,10 +1017,15 @@ void InterfaceKinetics::addGlobalReaction(ReactionData& rdata)
// store activation energy
m_E.push_back(rdata.rateCoeffParameters[2]);
//
// Add the reaction into the list of electrochemical extras
//
if (rdata.beta > 0.0) {
m_has_electrochem_rxns = true;
m_beta.push_back(rdata.beta);
// Push back the id of the reaction
m_ctrxn.push_back(m_ii);
// set the default to be the normal forward / reverse calculation method
m_ctrxn_BVform.push_back(0);
if (rdata.rateCoeffType == EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE) {
m_has_exchange_current_density_formulation = true;
@ -1083,6 +1160,7 @@ void InterfaceKinetics::init()
m_actConc.resize(m_kk);
m_conc.resize(m_kk);
m_mu0.resize(m_kk);
m_mu.resize(m_kk);
m_mu0_Kc.resize(m_kk);
m_grt.resize(m_kk);
m_pot.resize(m_kk, 0.0);
@ -1105,6 +1183,8 @@ void InterfaceKinetics::finalize()
m_StandardConc.resize(m_kk, 0.0);
m_deltaG0.resize(safe_reaction_size, 0.0);
m_deltaG.resize(safe_reaction_size, 0.0);
m_ProdStanConcReac.resize(safe_reaction_size, 0.0);
if (m_thermo.size() != m_phaseExists.size()) {
@ -1219,6 +1299,20 @@ void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStabl
}
}
//==================================================================================================================
// Write values into m_index
/*
* @param rxnNumber reaction number
* @param type reaction type
* @param loc location location in the reaction rate coefficient calculator
*/
void InterfaceKinetics::registerReaction(size_t rxnNumber, int type, size_t loc)
{
//
// type and loc is storred as a pair of values.
//
m_index[rxnNumber] = std::pair<int, size_t>(type, loc);
}
//==================================================================================================================
void EdgeKinetics::finalize()
{
deltaElectricEnergy_.resize(std::max<size_t>(m_ii, 1));
@ -1243,5 +1337,6 @@ void EdgeKinetics::finalize()
m_finalized = true;
}
//==================================================================================================================
}

3
src/kinetics/Kinetics.cpp
View file

@ -1,7 +1,6 @@
/**
* @file Kinetics.cpp Declarations for the base class for kinetics managers
* (see \ref kineticsmgr and class \link Cantera::Kinetics
* Kinetics\endlink).
* (see \ref kineticsmgr and class \link Cantera::Kinetics Kinetics \endlink).
*
* Kinetics managers calculate rates of progress of species due to
* homogeneous or heterogeneous kinetics.

99
src/kinetics/ReactionStoichMgr.cpp
View file

@ -45,7 +45,7 @@ ReactionStoichMgr& ReactionStoichMgr::operator=(const ReactionStoichMgr& right)
}
return *this;
}
//=========================================================================================================
void ReactionStoichMgr::add(size_t rxn, const std::vector<size_t>& reactants,
const std::vector<size_t>& products,
bool reversible)
@ -59,12 +59,65 @@ void ReactionStoichMgr::add(size_t rxn, const std::vector<size_t>& reactants,
m_irrevproducts.add(rxn, products);
}
}
//=========================================================================================================
// Add the reaction into the stoichiometric manager
void ReactionStoichMgr::add(size_t rxn, const ReactionData& r)
{
size_t k;
std::vector<size_t> rk;
doublereal frac;
doublereal oo, os, of;
bool doGlobal = false;
std::vector<size_t> extReactants = r.reactants;
vector_fp extRStoich = r.rstoich;
vector_fp extROrder = r.rorder;
//
// If we have a complete global reaction then we need to do something more complete
// than the previous treatment. Basically we will use the reactant manager to calculate the
// global forward reaction rate of progress.
//
if (r.forwardFullOrder_.size() > 0) {
//
// Trigger a treatment where the order of the reaction and the stoichiometry
// are treated as different.
//
doGlobal = true;
size_t nsp = r.forwardFullOrder_.size();
//
// Set up a signal vector to indicate whether the species has been added into
// the input vectors for the stoich manager
//
vector_int kHandled(nsp, 0);
//
// Loop over the reactants which are also nonzero stoichioemtric entries
// making sure the forwardFullOrder_ entries take precedence over rorder entries
//
for (size_t kk = 0; kk < r.reactants.size(); kk++) {
k = r.reactants[kk];
os = r.rstoich[kk];
oo = r.rorder[kk];
of = r.forwardFullOrder_[k];
if (of != oo) {
extROrder[kk] = of;
}
kHandled[k] = 1;
}
for (k = 0; k < nsp; k++) {
of = r.forwardFullOrder_[k];
if (of != 0.0) {
if (kHandled[k] == 0) {
//
// Add extra entries to reactant inputs. Set their reactant stoichiometric entries to zero.
//
extReactants.push_back(k);
extROrder.push_back(of);
extRStoich.push_back(0.0);
}
}
}
}
bool isfrac = false;
for (size_t n = 0; n < r.reactants.size(); n++) {
size_t ns = size_t(r.rstoich[n]);
@ -77,12 +130,22 @@ void ReactionStoichMgr::add(size_t rxn, const ReactionData& r)
}
}
// if the reaction has fractional stoichiometric coefficients
// or specified reaction orders, then add it in a general reaction
if (isfrac || r.global || rk.size() > 3) {
m_reactants.add(rxn, r.reactants, r.rorder, r.rstoich);
//
// If the reaction is non-mass action add it in in a general way
// Reactants get extra terms for the forward reaction rate of progress
// that may have zero stoichiometries.
//
if (doGlobal) {
m_reactants.add(rxn, extReactants, extROrder, extRStoich);
} else {
m_reactants.add(rxn, rk);
//
// this is confusing. The only issue should be whether rorder is different than rstoich!
//
if (isfrac || r.global || rk.size() > 3) {
m_reactants.add(rxn, r.reactants, r.rorder, r.rstoich);
} else {
m_reactants.add(rxn, rk);
}
}
std::vector<size_t> pk;
@ -99,22 +162,26 @@ void ReactionStoichMgr::add(size_t rxn, const ReactionData& r)
}
if (r.reversible) {
if (isfrac && !r.isReversibleWithFrac) {
throw CanteraError("ReactionStoichMgr::add",
"Fractional product stoichiometric coefficients only allowed "
"\nfor irreversible reactions and most reversible reactions");
}
//
// this is confusing. The only issue should be whether porder is different than pstoich!
//
if (pk.size() > 3 || r.isReversibleWithFrac) {
m_revproducts.add(rxn, r.products, r.porder, r.pstoich);
} else {
m_revproducts.add(rxn, pk);
}
} else if (isfrac || pk.size() > 3) {
m_irrevproducts.add(rxn, r.products, r.porder, r.pstoich);
} else {
m_irrevproducts.add(rxn, pk);
//
// this is confusing. The only issue should be whether porder is different than pstoich!
//
if (isfrac || pk.size() > 3) {
m_irrevproducts.add(rxn, r.products, r.porder, r.pstoich);
} else {
m_irrevproducts.add(rxn, pk);
}
}
}
//=========================================================================================================
void ReactionStoichMgr::getCreationRates(size_t nsp, const doublereal* ropf,
const doublereal* ropr, doublereal* c)

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

@ -0,0 +1,164 @@
/**
* @file example2.cpp
*
* $Id: RxnMolChange.cpp 571 2013-03-26 16:44:21Z hkmoffa $
*
*/
/*
* Copywrite (2005) Sandia Corporation. Under the terms of
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
* U.S. Government retains certain rights in this software.
*/
#include "cantera/kinetics/RxnMolChange.h"
#include "cantera/thermo.h"
#include "cantera/kinetics.h"
#include "cantera/kinetics/InterfaceKinetics.h"
#include "cantera/kinetics/ExtraGlobalRxn.h"
#include <iostream>
#include <new>
using namespace Cantera;
using namespace std;
namespace Cantera {
RxnMolChange::RxnMolChange(Cantera::Kinetics* kinPtr, int irxn) :
m_nPhases(0),
m_kinBase(kinPtr),
m_iRxn(irxn),
m_ChargeTransferInRxn(0.0),
m_beta(0.0),
m_egr(0)
{
int nReac = kinPtr->nReactions();
int iph;
AssertTrace(irxn >= 0);
AssertTrace(irxn < nReac);
m_nPhases = kinPtr->nPhases();
m_phaseMoleChange.resize(m_nPhases, 0.0);
m_phaseReactantMoles.resize(m_nPhases, 0.0);
m_phaseProductMoles.resize(m_nPhases, 0.0);
m_phaseMassChange.resize(m_nPhases, 0.0);
m_phaseChargeChange.resize(m_nPhases, 0.0);
m_phaseTypes.resize(m_nPhases, 0);
m_phaseDims.resize(m_nPhases, 0);
int m_kk = kinPtr->nTotalSpecies();
for (int kKin = 0; kKin < m_kk; kKin++) {
iph = m_kinBase->speciesPhaseIndex(kKin);
Cantera::ThermoPhase& tpRef = m_kinBase->thermo(iph);
int kLoc = kKin - m_kinBase->kineticsSpeciesIndex(0, iph);
double rsc = m_kinBase->reactantStoichCoeff(kKin, irxn);
double psc = m_kinBase->productStoichCoeff(kKin, irxn);
double nsc = psc - rsc;
m_phaseMoleChange[iph] += (nsc);
m_phaseReactantMoles[iph] += rsc;
m_phaseProductMoles[iph] += psc;
double mw = tpRef.molecularWeight(kLoc);
m_phaseMassChange[iph] += (nsc) * mw;
double chg = tpRef.charge(kLoc);
m_phaseChargeChange[iph] += nsc * chg;
}
for (iph = 0; iph < m_nPhases; iph++) {
Cantera::ThermoPhase& tpRef = m_kinBase->thermo(iph);
m_phaseDims[iph] = tpRef.nDim();
m_phaseTypes[iph] = tpRef.eosType();
if (m_phaseChargeChange[iph] != 0.0) {
double tmp = fabs(m_phaseChargeChange[iph]);
if (tmp > m_ChargeTransferInRxn) {
m_ChargeTransferInRxn = tmp;
}
}
}
if (m_ChargeTransferInRxn) {
Cantera::InterfaceKinetics* iK = dynamic_cast<Cantera::InterfaceKinetics*>(kinPtr);
if (iK) {
m_beta = iK->electrochem_beta(irxn);
} else {
throw Cantera::CanteraError("RxnMolChange", "unknown condition on charge");
}
}
}
RxnMolChange::RxnMolChange(Cantera::Kinetics* kinPtr, Cantera::ExtraGlobalRxn* egr) :
m_nPhases(0),
m_kinBase(kinPtr),
m_iRxn(-1),
m_ChargeTransferInRxn(0.0),
m_beta(0.0),
m_egr(egr)
{
int iph;
AssertTrace(egr != 0);
m_nPhases = kinPtr->nPhases();
m_phaseMoleChange.resize(m_nPhases, 0.0);
m_phaseReactantMoles.resize(m_nPhases, 0.0);
m_phaseProductMoles.resize(m_nPhases, 0.0);
m_phaseMassChange.resize(m_nPhases, 0.0);
m_phaseChargeChange.resize(m_nPhases, 0.0);
m_phaseTypes.resize(m_nPhases, 0);
m_phaseDims.resize(m_nPhases, 0);
int m_kk = kinPtr->nTotalSpecies();
for (int kKin = 0; kKin < m_kk; kKin++) {
iph = m_kinBase->speciesPhaseIndex(kKin);
ThermoPhase& tpRef = m_kinBase->thermo(iph);
int kLoc = kKin - m_kinBase->kineticsSpeciesIndex(0, iph);
double rsc = egr->reactantStoichCoeff(kKin);
double psc = egr->productStoichCoeff(kKin);
double nsc = psc - rsc;
m_phaseMoleChange[iph] += (nsc);
m_phaseReactantMoles[iph] += rsc;
m_phaseProductMoles[iph] += psc;
double mw = tpRef.molecularWeight(kLoc);
m_phaseMassChange[iph] += (nsc) * mw;
double chg = tpRef.charge(kLoc);
m_phaseChargeChange[iph] += nsc * chg;
}
for (iph = 0; iph < m_nPhases; iph++) {
ThermoPhase& tpRef = m_kinBase->thermo(iph);
m_phaseDims[iph] = tpRef.nDim();
m_phaseTypes[iph] = tpRef.eosType();
if (m_phaseChargeChange[iph] != 0.0) {
double tmp = fabs(m_phaseChargeChange[iph]);
if (tmp > m_ChargeTransferInRxn) {
m_ChargeTransferInRxn = tmp;
}
}
}
if (m_ChargeTransferInRxn) {
InterfaceKinetics* iK = dynamic_cast<InterfaceKinetics*>(kinPtr);
if (iK) {
m_beta = 0.0;
} else {
throw CanteraError("RxnMolChange", "unknown condition on charge");
}
}
}
RxnMolChange::~RxnMolChange()
{
}
}

2
src/kinetics/importKinetics.cpp
View file

@ -82,7 +82,7 @@ public:
}
}
};
//=======================================================================================================
void checkRxnElementBalance(Kinetics& kin,
const ReactionData& rdata, doublereal errorTolerance)
{