From d86a6bf09799e98342bb5d8653fb04739d0c1752 Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Tue, 3 Feb 2004 03:33:42 +0000 Subject: [PATCH] initial import --- Cantera/src/EdgeKinetics.cpp | 506 +++++++++++++++++++++++++++++++++++ Cantera/src/EdgeKinetics.h | 396 +++++++++++++++++++++++++++ Cantera/src/EdgePhase.h | 40 +++ 3 files changed, 942 insertions(+) create mode 100644 Cantera/src/EdgeKinetics.cpp create mode 100644 Cantera/src/EdgeKinetics.h create mode 100644 Cantera/src/EdgePhase.h diff --git a/Cantera/src/EdgeKinetics.cpp b/Cantera/src/EdgeKinetics.cpp new file mode 100644 index 000000000..c893b1132 --- /dev/null +++ b/Cantera/src/EdgeKinetics.cpp @@ -0,0 +1,506 @@ +/** + * @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 +using namespace std; + + +namespace Cantera { + + ////////////////////////////////////////////////////////////////// + + /** + * Construct an empty EdgeKinetics reaction mechanism. + * @param thermo This is an optional parameter that may be + * used to initialize the inherited Kinetics class with + * one ThermoPhase class object -> in other words it's + * useful for initialization of homogeneous kinetics + * mechanisms. + */ + EdgeKinetics:: + EdgeKinetics(thermo_t* thermo) : + Kinetics(thermo), + m_kk(0), + m_redo_rates(false), + m_nirrev(0), + m_nrev(0), + m_finalized(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, m_kdata->m_rfn.begin()); + applyButlerVolmerCorrection(m_kdata->m_rfn.begin()); + 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(m_conc.begin() + 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); + + 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(m_mu0.begin() + 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_reactantStoich.decrementReactions(m_mu0.begin(), m_rkc.begin()); + m_revProductStoich.incrementReactions(m_mu0.begin(), m_rkc.begin()); + + 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(dmu.begin() + 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_reactantStoich.decrementReactions(dmu.begin(), rmu.begin()); + m_revProductStoich.incrementReactions(dmu.begin(), rmu.begin()); + + 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(m_mu0.begin() + 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_reactantStoich.decrementReactions(m_mu0.begin(), kc); + m_revProductStoich.incrementReactions(m_mu0.begin(), kc); + m_irrevProductStoich.incrementReactions(m_mu0.begin(), 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, ...). 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. + fill(m_rwork.begin(), m_rwork.begin() + m_ii, 0.0); + m_reactantStoich.decrementReactions(m_pot.begin(), m_rwork.begin()); + m_revProductStoich.incrementReactions(m_pot.begin(), m_rwork.begin()); + m_irrevProductStoich.incrementReactions(m_pot.begin(), m_rwork.begin()); + + // 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; + + for (i = 0; i < m_ii; i++) { + eamod = 0.5*m_rwork[i]; + if (eamod != 0.0 && m_E[i] != 0.0) { + ea = GasConstant * m_E[i]; + if (eamod + ea < 0.0) eamod = -ea; + kf[i] *= 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(m_conc.begin(), ropf.begin()); + + // for reversible reactions, multiply ropr by concentration + // products + m_revProductStoich.multiply(m_conc.begin(), ropr.begin()); + + // do global reactions + m_globalReactantStoich.power(m_conc.begin(), 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; + 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(), + rp.begin() ); + // store activation energy + m_E.push_back(r.rateCoeffParameters[2]); + // 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(), + rp.begin() ); + + // 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; + + int rnum = reactionNumber(); + + vector_int rk; + int nr = r.reactants.size(); + for (n = 0; n < nr; n++) { + ns = r.rstoich[n]; + 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++) { + ns = r.pstoich[n]; + 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& r, const vector& 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); + } + +} + + + + + + + + diff --git a/Cantera/src/EdgeKinetics.h b/Cantera/src/EdgeKinetics.h new file mode 100644 index 000000000..288a189a6 --- /dev/null +++ b/Cantera/src/EdgeKinetics.h @@ -0,0 +1,396 @@ +/** + * @file EdgeKinetics.h + * + * $Author$ + * $Revision$ + * $Date$ + */ + +// Copyright 2001 California Institute of Technology + + +#ifndef CT_EDGEKINETICS_H +#define CT_EDGEKINETICS_H + +#include +#include +#include +#include + +#include "mix_defs.h" +#include "Kinetics.h" + +#include "utilities.h" +#include "RateCoeffMgr.h" +#include "StoichManager.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 { + + 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. + */ + EdgeKinetics(thermo_t* thermo = 0); + + /// Destructor. + virtual ~EdgeKinetics(); + + /** + * Identifies the subclass of the Kinetics manager type. + * These are listed in mix_defs.h. + */ + virtual int ID() { return cEdgeKinetics; } + + /** + * Identifies the subclass of the Kinetics manager type. + * These are listed in mix_defs.h. + */ + 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(); + 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(); + 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(); + 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(); + fill(cdot, cdot + m_kk, 0.0); + m_revProductStoich.incrementSpecies( + m_kdata->m_ropf.begin(), cdot); + m_irrevProductStoich.incrementSpecies( + m_kdata->m_ropf.begin(), cdot); + m_reactantStoich.incrementSpecies( + m_kdata->m_ropr.begin(), 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(); + fill(ddot, ddot + m_kk, 0.0); + m_revProductStoich.incrementSpecies( + m_kdata->m_ropr.begin(), ddot); + m_reactantStoich.incrementSpecies( + m_kdata->m_ropf.begin(), 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(); + fill(net, net + m_kk, 0.0); + m_revProductStoich.incrementSpecies( + m_kdata->m_ropnet.begin(), net); + m_irrevProductStoich.incrementSpecies( + m_kdata->m_ropnet.begin(), net); + m_reactantStoich.decrementSpecies( + m_kdata->m_ropnet.begin(), 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 (find(m_revindex.begin(), m_revindex.end(), i) + < m_revindex.end()) return true; + else return false; + } + + /** + * Return a string representing the reaction. + */ + virtual 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. + */ + virtual void finalize(); + virtual bool ready() const; + + + void updateROP(); + + + const vector& reactantGroups(int i) + { return m_rgroups[i]; } + const vector& 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 m_rates; + //Rate1 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 map > m_index; + + vector 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; + + map > m_rgroups; + map > m_pgroups; + + vector m_rxntype; + + mutable vector > m_rrxn; + mutable vector > m_prxn; + + vector_int m_revindex; + vector 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; + + 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 vector& r, + const vector& p); + void updateKc(); + + void registerReaction(int rxnNumber, int type, int loc) { + m_index[rxnNumber] = pair(type, loc); + } + void applyButlerVolmerCorrection(doublereal* kf); + bool m_finalized; + }; +} + +#endif diff --git a/Cantera/src/EdgePhase.h b/Cantera/src/EdgePhase.h new file mode 100644 index 000000000..3cb665a0f --- /dev/null +++ b/Cantera/src/EdgePhase.h @@ -0,0 +1,40 @@ +/** + * + * @file EdgePhase.h + * + */ + +/* $Author$ + * $Date$ + * $Revision$ + * + * Copyright 2002 California Institute of Technology + * + */ + + +#ifndef CT_EDGEPHASE_H +#define CT_EDGEPHASE_H + +#include "mix_defs.h" +#include "ThermoPhase.h" + +namespace Cantera { + + class EdgePhase : public SurfPhase { + + public: + + EdgePhase(doublereal n0 = 0.0); + virtual ~EdgePhase() {} + virtual int eosType() const { return cEdge; } + + }; +} + +#endif + + + + +