diff --git a/CMakeLists.txt b/CMakeLists.txt index 22a09a252..2a3d56dcd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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) diff --git a/Cantera/python/Cantera/Kinetics.py b/Cantera/python/Cantera/Kinetics.py index 905b68fa6..7bb706f1d 100755 --- a/Cantera/python/Cantera/Kinetics.py +++ b/Cantera/python/Cantera/Kinetics.py @@ -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) diff --git a/Cantera/python/src/ctkinetics_methods.cpp b/Cantera/python/src/ctkinetics_methods.cpp index a5209b51d..6f0fbdd25 100644 --- a/Cantera/python/src/ctkinetics_methods.cpp +++ b/Cantera/python/src/ctkinetics_methods.cpp @@ -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: ; diff --git a/Cantera/src/kinetics/EdgeKinetics.cpp b/Cantera/src/kinetics/EdgeKinetics.cpp deleted file mode 100644 index 9497b54aa..000000000 --- a/Cantera/src/kinetics/EdgeKinetics.cpp +++ /dev/null @@ -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 -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& 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/kinetics/EdgeKinetics.h b/Cantera/src/kinetics/EdgeKinetics.h index c2f6e1b4d..f1eaf2de8 100644 --- a/Cantera/src/kinetics/EdgeKinetics.h +++ b/Cantera/src/kinetics/EdgeKinetics.h @@ -12,51 +12,11 @@ #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" +#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& reactantGroups(int i) - { return m_rgroups[i]; } - const std::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 std::map > m_index; - - std::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; - - std::map > m_rgroups; - std::map > m_pgroups; - - std::vector m_rxntype; - - mutable std::vector > m_rrxn; - mutable std::vector > m_prxn; - - vector_int m_revindex; - std::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; - 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& r, - const std::vector& p); - void updateKc(); - - void registerReaction(int rxnNumber, int type, int loc) { - m_index[rxnNumber] = std::pair(type, loc); - } - void applyButlerVolmerCorrection(doublereal* kf); - bool m_finalized; - bool m_has_electrochem_rxns; }; } diff --git a/Cantera/src/kinetics/InterfaceKinetics.cpp b/Cantera/src/kinetics/InterfaceKinetics.cpp index a1a3f302b..d2271890a 100644 --- a/Cantera/src/kinetics/InterfaceKinetics.cpp +++ b/Cantera/src/kinetics/InterfaceKinetics.cpp @@ -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; + } } diff --git a/Cantera/src/kinetics/InterfaceKinetics.h b/Cantera/src/kinetics/InterfaceKinetics.h index 90bb62437..18984d8e5 100644 --- a/Cantera/src/kinetics/InterfaceKinetics.h +++ b/Cantera/src/kinetics/InterfaceKinetics.h @@ -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 > 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: + }; } diff --git a/Cantera/src/kinetics/Makefile.in b/Cantera/src/kinetics/Makefile.in index 3ac823355..40b1e4912 100644 --- a/Cantera/src/kinetics/Makefile.in +++ b/Cantera/src/kinetics/Makefile.in @@ -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 diff --git a/Cantera/src/thermo/CMakeLists.txt b/Cantera/src/thermo/CMakeLists.txt index ecbd95d81..b6dc45f02 100644 --- a/Cantera/src/thermo/CMakeLists.txt +++ b/Cantera/src/thermo/CMakeLists.txt @@ -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}) diff --git a/cmake/fortran.cmake b/cmake/fortran.cmake index 413a56f59..21e62eb18 100644 --- a/cmake/fortran.cmake +++ b/cmake/fortran.cmake @@ -3,6 +3,7 @@ #if (NOT BUILD_WITH_F2C) #### Fortran 90 +message("Fortran cmake") if (BUILD_F90_INTERFACE) if (F90 STREQUAL "default") diff --git a/config.cmake b/config.cmake index 521db17c2..b1a3a212e 100755 --- a/config.cmake +++ b/config.cmake @@ -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