From 6741d8f7c6715b2eb66502b3cfc0f5de764fad91 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Tue, 12 Aug 2014 23:02:47 +0000 Subject: [PATCH] InterfaceKinetics rewrite -> addition of general BV reactions and affinity formulation. --- doc/doxygen/Cantera.txt | 3 + doc/doxygen/thermo.txt | 95 +++++++ include/cantera/base/utilities.h | 1 + include/cantera/kinetics/InterfaceKinetics.h | 40 ++- include/cantera/kinetics/KineticsFactory.h | 3 +- include/cantera/kinetics/ReactionData.h | 32 ++- include/cantera/kinetics/StoichManager.h | 2 + include/cantera/kinetics/importKinetics.h | 33 ++- include/cantera/kinetics/reaction_defs.h | 35 ++- include/cantera/thermo/Phase.h | 1 - src/kinetics/InterfaceKinetics.cpp | 31 ++- src/kinetics/importKinetics.cpp | 275 +++++++++++++++---- src/thermo/GibbsExcessVPSSTP.cpp | 5 +- 13 files changed, 449 insertions(+), 107 deletions(-) create mode 100644 doc/doxygen/thermo.txt diff --git a/doc/doxygen/Cantera.txt b/doc/doxygen/Cantera.txt index 63c13a8f8..847b6128a 100644 --- a/doc/doxygen/Cantera.txt +++ b/doc/doxygen/Cantera.txt @@ -3,4 +3,7 @@ Use the menu at the top to view detailed documentation of the code. +\ref thermopage + + */ diff --git a/doc/doxygen/thermo.txt b/doc/doxygen/thermo.txt new file mode 100644 index 000000000..cba0b3e18 --- /dev/null +++ b/doc/doxygen/thermo.txt @@ -0,0 +1,95 @@ +/** + +\page thermopage Thermodynamic Properties + +%Cantera can be used to compute thermodynamic properties of pure +substances, solutions, and mixtures of various types, including ones +containing multiple phases. The first step is to create an object that +represents each phase. + +A simple complete program that creates an object representing a gas mixture and prints its temperature is shown below. +\include ex1.cpp + +Class \link Cantera::ThermoPhase ThermoPhase \endlink +is the base class for %Cantera classes that represent +phases of matter. It defines the public interface for all classes that +represent phases. For example, it specifies that they all have a +method \c temperature() that returns the current temperature, a method +\c setTemperature(double T) that sets the temperature, a method \c +getChemPotentials(double* mu) that writes the species chemical +potentials into array \c mu, and so on. + +Class ThermoPhase can be used to represent the intensive state of any +single-phase solution of multiple species. The phase may be a bulk, +three-dimensional phase (a gas, a liquid, or a solid), or may be a +two-dimensional surface phase, or even a one-dimensional "edge" +phase. The specific attributes of each type of phase are specified by +deriving a class from ThemoPhase and providing implementations for the +virtual methods of ThermoPhase. + +\section The Intensive Thermodynamic State + +Class ThermoPhase and classes derived from it work only with the +intensive thermodynamic state. That is, all extensive properties +(enthalpy, entropy, internal energy, volume, etc.) are computed for a +unit quantity (on a mass or mole basis). For example, there is a +method enthalpy_mole() that returns the molar enthalpy (J/kmol), and a +method enthalpy_mass() that returns the specific enthalpy (J/kg), but +no method enthalpy() that would return the total enthalpy (J). This is +because class ThermoPhase does not store the total amount (mass or +mole) of the phase. + +From thermodynamics, it may be shown that the intensive state of a +single-component phase in equilibrium is fully specified by the values +of any r+1 independent thermodynamic properties, where r is the number +of reversible work modes. If the only reversible work mode is +compression (a "simple compressible substance"), then two properties +suffice to specify the intensive state. + +In principle, any two independent p + +specified, the values of all other intensive properties may be +computed. For example, specifying the pressure and molar entropy + +consisting of a solution of K species +in equilibrium is fully specified by the values of any two independent +thermodynamic properties, in addition to in +Class ThermoPhase stores internally the values of the temperature, the +mass density, and the mass fractions of all species. These values are +sufficient to fix the intensive thermodynamic state of the phase. All +properties for a unit amount (on a mass or mole basis) are determined +once the intensive state is specified. For the extensive properties, class ThermoPhase provides methods that return property values on a molar basis (e.g. enthalpy_mole(), with units J/kmol) or on a mass basis (e.g. enthalpy_mass(), with units J/kg). Since the total mass or total number of moles is not stored, + +Note that the total mass or number of moles is not stored + +Given these values, any other intensive thermodynamic property may + +Note that the total mass or total number of moles is not stored -- therefore the values of all extensive properties (mass, volume, energy) are + +This choice is arbitrary, and for most purposes you can't tell which properties are stored and which are computed. + +The classes that derive from ThermoPhase compute o + +For example, suppose we want to create a class to use to compute the properties of ideal gas mixtures. + +Many of the methods of ThermoPhase are declared virtual, and are meant to be +overloaded in classes derived from ThermoPhase. For example, class \link Cantera::IdealGasPhase IdealGasPhase \endlink +derives from ThermoPhase, and represents ideal gas mixtures. + +Although class ThermoPhase defines the interface for all classes +representing phases, it only provides implementations for a few of the +methods. This is because ThermoPhase does not actually know the +equation of state of any phase -- this information is provided by +classes that derive from ThermoPhase. +The methods implemented by ThermoPhase are ones that apply to all phases, independent of +the equation of state. For example, it implements methods temperature() and setTemperature(), +since the temperature value is stored internally. Also, the mass density is stored internally, so + +There is a list of classes which inherit from the ThermoPhase class (see \ref +thermoprops "Thermodynamic Properties") + +There is a list of classes which handle standard states for species (see +\ref spthermo "Species Standard-State Thermodynamic Properties"). + + +*/ diff --git a/include/cantera/base/utilities.h b/include/cantera/base/utilities.h index c9565558f..005c96731 100644 --- a/include/cantera/base/utilities.h +++ b/include/cantera/base/utilities.h @@ -19,6 +19,7 @@ #include "ct_defs.h" #include "global.h" +#include #include namespace Cantera diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index 0433f1a71..6977e8033 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/include/cantera/kinetics/InterfaceKinetics.h @@ -313,9 +313,9 @@ public: void applyButlerVolmerCorrection(doublereal* const kf); //! When an electrode reaction rate is optionally specified in terms of its - //! exchange current density, adjust to standard reaction rate form. - /** - * For a reaction rate that was given in units of Amps/m2 (exchange current + //! exchange current density, adjust kfwd to the standard reaction rate constant form and units. + /*! + * For a reaction rate constant that was given in units of Amps/m2 (exchange current * density formulation with iECDFormulation == true), convert the rate to * kmoles/m2/s. */ @@ -524,6 +524,12 @@ protected: //! Pointer to the single surface phase SurfPhase* m_surf; + //! Vector of reaction types + /*! + * Length = m_ii the number of reactions in the mechanism. + */ + vector_int reactionTypes_; + //! Pointer to the Implicit surface chemistry object /*! * Note this object is owned by this InterfaceKinetics object. It may only @@ -532,20 +538,42 @@ protected: */ ImplicitSurfChem* m_integrator; + //! Electrochemical transfer coefficient for the forward direction + /*! + * Electrochemical transfer coefficient for all reactions that have transfer reactions + * the reaction is given by m_ctrxn[i] + */ vector_fp m_beta; //! Vector of reaction indexes specifying the id of the current transfer //! reactions in the mechanism /*! * Vector of reaction indices which involve current transfers. This provides - * an index into the m_beta array. + * an index into the m_beta, ctrxn_BVform array. * * irxn = m_ctrxn[i] */ std::vector m_ctrxn; - //! Vector of booleans indicating whether the charge transfer reaction may - //! be described by an exchange current density expression + //! Vector of Reactions which follow the butler volmer methodology for specifying the + //! exchange current density first. Then, the other forms are specified based on this form. + /*! + * Length is equal to the number of reactions with charge transfer coefficients, m_ctrxn[] + * + * m_ctrxn_BVform[i] = 0; This means that the irxn reaction is calculated via the standard forward + * and reverse reaction rates + * m_ctrxn_BVform[i] = 1; This means that the irxn reaction is calculated via the BV format + * directly. + * m_ctrxn_BVform[i] = 2; this means that the irxn reaction is calculated via the BV format + * directly, using concentrations instead of activity concentrations. + */ + std::vector m_ctrxn_BVform; + + //! Vector of booleans indicating whether the charge transfer reaction rate constant + //! is described by an exchange current density rate constant expression + /*! + * Length is equal to the number of reactions with charge transfer coefficients, m_ctrxn[] + */ vector_int m_ctrxn_ecdf; //! Vector of standard concentrations diff --git a/include/cantera/kinetics/KineticsFactory.h b/include/cantera/kinetics/KineticsFactory.h index 0ad0e36fb..f7fef2646 100644 --- a/include/cantera/kinetics/KineticsFactory.h +++ b/include/cantera/kinetics/KineticsFactory.h @@ -64,8 +64,7 @@ public: * * @return Pointer to the new kinetics manager. */ - virtual Kinetics* newKinetics(XML_Node& phase, - std::vector th); + virtual Kinetics* newKinetics(XML_Node& phase, std::vector th); /** * Return a new, empty kinetics manager. diff --git a/include/cantera/kinetics/ReactionData.h b/include/cantera/kinetics/ReactionData.h index 10774eeb9..28587c2d3 100644 --- a/include/cantera/kinetics/ReactionData.h +++ b/include/cantera/kinetics/ReactionData.h @@ -52,20 +52,40 @@ public: std::vector products; //!< Indices of product species //! Reaction order with respect to each reactant species, in the order - //! given by #reactants. Usually the same as the stoichiometric - //! coefficients. + //! given by #reactants. Usually the same as the stoichiometric coefficients. + /*! + * Length is equal to the number of reactants defined in the reaction + * The order of species is given by the reactants vectors. + */ vector_fp rorder; //! Reaction order of the reverse reaction with respect to each product - //! species, in the order given by #products. Usually the same as the - //! stoichiometric coefficients. + //! species, in the order given by #products. Usually the same as the stoichiometric coefficients. + /*! + * Length is equal to the number of products defined in the reaction. + * The order of species is given by the products vectors. + */ vector_fp porder; - //! Reactant stoichiometric coefficients, in the order given by - //! #reactants. + //! Reaction order for the forward direction of the reaction + /*! + * Length is equal to the number of kinetic species defined in the kinetics object + * The order of species is given by kinetics species vector. + */ + vector_fp forwardFullOrder_; + + //! Reactant stoichiometric coefficients, in the order given by #reactants. + /*! + * Length is equal to the number of products defined in the reaction. + * The order of species is given by the products vectors. + */ vector_fp rstoich; //! Product stoichiometric coefficients, in the order given by #products. + /*! + * Length is equal to the number of products defined in the reaction. + * The order of species is given by the products vectors. + */ vector_fp pstoich; std::vector rgroups; //!< Optional data used in reaction path diagrams diff --git a/include/cantera/kinetics/StoichManager.h b/include/cantera/kinetics/StoichManager.h index b93bf644d..dd9832408 100644 --- a/include/cantera/kinetics/StoichManager.h +++ b/include/cantera/kinetics/StoichManager.h @@ -908,6 +908,8 @@ private: /** * Map with the Reaction Number as key and the placement in the * vector of reactions list( i.e., m_c1_list[]) as key + * If for example, m_loc[7], was equal to 5, this means that the 7th overall reaction in the mechanism + * is located in the 5th position of m_c1_list if it unimolecular and only has one reactant/product. */ std::map m_loc; }; diff --git a/include/cantera/kinetics/importKinetics.h b/include/cantera/kinetics/importKinetics.h index 2b66bdf07..6ba288c09 100644 --- a/include/cantera/kinetics/importKinetics.h +++ b/include/cantera/kinetics/importKinetics.h @@ -12,8 +12,8 @@ // Copyright 2002 California Institute of Technology -#ifndef CT_IMPORTCTML_H -#define CT_IMPORTCTML_H +#ifndef CT_IMPORTKINETICS_H +#define CT_IMPORTKINETICS_H #include "cantera/thermo/ThermoPhase.h" #include "Kinetics.h" @@ -83,11 +83,19 @@ void checkRxnElementBalance(Kinetics& kin, * routine to skip this reaction and continue. Otherwise, we * will throw an error. */ -bool getReagents(const XML_Node& rxn, Kinetics& kin, int rp, - std::string default_phase, +bool getReagents(const XML_Node& rxn, Kinetics& kin, int rp, std::string default_phase, std::vector& spnum, vector_fp& stoich, vector_fp& order, const ReactionRules& rules); +//! Get non-mass-action orders for a reaction + + +extern bool getOrders(const XML_Node& rxnNode, Kinetics& kin, + std::string default_phase, const ReactionData& rdata, + vector_fp& order, vector_fp& fullForwardsOrders, + const ReactionRules& rules); + + //! Read the rate coefficient data from the XML file. /*! * Extract the rate coefficient for a reaction from the xml node, kf. @@ -159,7 +167,7 @@ bool installReactionArrays(const XML_Node& p, Kinetics& kin, * * @param kin This is a pointer to a kinetics manager class that will be * initialized with the kinetics mechanism. Inherited Kinetics - * classes may be used here. + * classes should be used here. * * @ingroup kineticsmgr * @@ -167,10 +175,11 @@ bool installReactionArrays(const XML_Node& p, Kinetics& kin, bool importKinetics(const XML_Node& phase, std::vector th, Kinetics* kin); -//!Build a single-phase ThermoPhase object with associated kinetics mechanism. +//! Build a single-phase ThermoPhase object with associated kinetics mechanism. /*! * In a single call, this routine initializes a ThermoPhase object and a - * homogeneous kinetics object for a phase. + * homogeneous kinetics object for a phase. It returns the fully initialized + * ThermoPhase object ptr and kinetics ptr. * * @param root pointer to the XML tree which will be searched to find the * XML phase element. @@ -179,7 +188,7 @@ bool importKinetics(const XML_Node& phase, std::vector th, * @param nm Name of the XML element. Should be "phase" * @param th Pointer to a bare ThermoPhase object, which will be initialized * by this operation. - * @param k Pointer to a bare Kinetics object, which will be initialized + * @param kin Pointer to a bare Kinetics object, which will be initialized * by this operation to a homogeneous kinetics manager * * @return @@ -188,17 +197,17 @@ bool importKinetics(const XML_Node& phase, std::vector th, * For Example * * @code - * ThermoPhase *th = new ThermoPhase(); - * Kinetics *k = new Kinetics(); + * ThermoPhase *th = new ThermoPhase(); + * Kinetics *kin = new Kinetics(); * XML_Node *root = get_XML_File("gri30.xml"); - * ok = buildSolutionFromXML(root, "gri30_mix", "phase", th, k) + * ok = buildSolutionFromXML(root, "gri30_mix", "phase", th, kin) * @endcode * * @ingroup inputfiles * @see importKinetics() */ bool buildSolutionFromXML(XML_Node& root, const std::string& id, - const std::string& nm, ThermoPhase* th, Kinetics* k); + const std::string& nm, ThermoPhase* th, Kinetics* kin); //! Search an XML tree for species data. /*! diff --git a/include/cantera/kinetics/reaction_defs.h b/include/cantera/kinetics/reaction_defs.h index 567d9618c..a2f0d48de 100644 --- a/include/cantera/kinetics/reaction_defs.h +++ b/include/cantera/kinetics/reaction_defs.h @@ -32,7 +32,7 @@ const int ELEMENTARY_RXN = 1; const int THREE_BODY_RXN = 2; /** - * The general form for an association or dissociation reaction, with a + * The general form for a gas-phase association or dissociation reaction, with a * pressure-dependent rate. Example: CH3 + H (+M) <-> CH4 (+M) */ const int FALLOFF_RXN = 4; @@ -46,7 +46,7 @@ const int FALLOFF_RXN = 4; const int PLOG_RXN = 5; /** - * A general pressure-dependent reaction where k(T,P) is defined in terms of + * A general gas-phase pressure-dependent reaction where k(T,P) is defined in terms of * a bivariate Chebyshev polynomial. */ const int CHEBYSHEV_RXN = 6; @@ -61,17 +61,42 @@ const int CHEMACT_RXN = 8; /** * A reaction occurring on a surface. + * NOTE: This is a bit ambiguous, and will be taken out in the future + * The dimensionality of the interface is a separate concept from the type + * of the reaction. */ const int SURFACE_RXN = 20; +//! This is a surface reaction that is formulated using the Butler-Volmer +//! formulation and using concentrations instead of activity concentrations +//! for its exchange current density formulat. +const int BUTLERVOLMER_NOACTIVITYCOEFFS_RXN = 25; + +//! This is a surface reaction that is formulated using the Butler-Volmer +//! formulation. Note the B-V equations can be derived from the forward +//! and reverse rate constants for a single step reaction. However, there +//! are some advantages to using the formulation directly. +const int BUTLERVOLMER_RXN = 26; + +//! This is a surface reaction that is formulated using the affinity +//! representation, common in the geochemistry community. +//! This is generally a global non-mass action reaction with an additional functional +//! form dependence on delta G of reaction. +const int SURFACEAFFINITY_RXN = 27; + + + + /** - * A reaction occurring at a one-dimensional interface between two - * surface phases. + * A reaction occurring at a one-dimensional interface between two surface phases. + * NOTE: This is a bit ambiguous, and will be taken out in the future + * The dimensionality of the interface is a separate concept from the type + * of the reaction. */ const int EDGE_RXN = 22; /** - * A global reaction. These may have non-integral reaction orders, + * A global reaction. These may have non-mass action reaction orders, * and are not allowed to be reversible. */ const int GLOBAL_RXN = 30; diff --git a/include/cantera/thermo/Phase.h b/include/cantera/thermo/Phase.h index c6016a343..bdd1a19d0 100644 --- a/include/cantera/thermo/Phase.h +++ b/include/cantera/thermo/Phase.h @@ -7,7 +7,6 @@ #ifndef CT_PHASE_H #define CT_PHASE_H -#include "cantera/Cantera.h" #include "cantera/base/vec_functions.h" #include "cantera/base/ctml.h" #include "cantera/thermo/Elements.h" diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index 3a520983f..f67e41f24 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -199,6 +199,9 @@ void InterfaceKinetics::_update_rates_T() void InterfaceKinetics::_update_rates_phi() { + // + // Store electric potentials for each phase in the array m_phi[]. + // for (size_t n = 0; n < nPhases(); n++) { if (thermo(n).electricPotential() != m_phi[n]) { m_phi[n] = thermo(n).electricPotential(); @@ -222,7 +225,7 @@ void InterfaceKinetics::_update_rates_C() } m_ROP_ok = false; } - +//============================================================================================================================ void InterfaceKinetics::getActivityConcentrations(doublereal* const conc) { _update_rates_C(); @@ -310,13 +313,13 @@ void InterfaceKinetics::checkPartialEquil() } } } - +//============================================================================================================================ void InterfaceKinetics::getFwdRatesOfProgress(doublereal* fwdROP) { updateROP(); std::copy(m_ropf.begin(), m_ropf.end(), fwdROP); } - +//============================================================================================================================ void InterfaceKinetics::getRevRatesOfProgress(doublereal* revROP) { updateROP(); @@ -399,7 +402,7 @@ void InterfaceKinetics::applyButlerVolmerCorrection(doublereal* const kf) for (size_t n = 0; n < nPhases(); n++) { size_t nsp = thermo(n).nSpecies(); for (size_t k = 0; k < nsp; k++) { - m_pot[ik] = Faraday*thermo(n).charge(k)*m_phi[n]; + m_pot[ik] = Faraday * thermo(n).charge(k) * m_phi[n]; ik++; } } @@ -446,7 +449,7 @@ void InterfaceKinetics::applyButlerVolmerCorrection(doublereal* const kf) } } } - +//================================================================================================================== /* * For a reaction rate that was given in units of Amps/m2 (exchange current * density formulation with iECDFormulation == true), convert the rate to @@ -469,7 +472,7 @@ void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* con } } } - +//================================================================================================================== void InterfaceKinetics::getFwdRateConstants(doublereal* kfwd) { @@ -482,7 +485,7 @@ void InterfaceKinetics::getFwdRateConstants(doublereal* kfwd) multiply_each(kfwd, kfwd + nReactions(), m_perturb.begin()); } - +//================================================================================================================== void InterfaceKinetics::getRevRateConstants(doublereal* krev, bool doIrreversible) { getFwdRateConstants(krev); @@ -520,8 +523,10 @@ void InterfaceKinetics::updateROP() // the forward rates copied into m_ropr by the reciprocals of // the equilibrium constants multiply_each(m_ropr.begin(), m_ropr.end(), m_rkcn.begin()); - - // multiply ropf by concentration products + // + // multiply ropf by the actyivity concentration reaction orders to obtain + // the forward rates of progress. + // m_rxnstoich.multiplyReactants(DATA_PTR(m_conc), DATA_PTR(m_ropf)); // for reversible reactions, multiply ropr by concentration @@ -720,7 +725,7 @@ void InterfaceKinetics::getDeltaSSEntropy(doublereal* deltaS) */ m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaS); } - +//============================================================================================================================ void InterfaceKinetics::addReaction(ReactionData& r) { /* @@ -745,6 +750,7 @@ void InterfaceKinetics::addReaction(ReactionData& r) * increase the size of m_perturb by one as well. */ incrementRxnCount(); + // m_rxneqn.push_back(r.equation); m_rxnPhaseIsReactant.push_back(std::vector(nPhases(), false)); @@ -764,6 +770,7 @@ void InterfaceKinetics::addReaction(ReactionData& r) m_rxnPhaseIsProduct[i][p] = true; } } +//============================================================================================================================ void InterfaceKinetics::addElementaryReaction(ReactionData& r) { @@ -791,7 +798,7 @@ void InterfaceKinetics::addElementaryReaction(ReactionData& r) /* * Install the reaction rate into the vector of reactions handled by this class */ - size_t iloc = m_rates.install(reactionNumber(), r); + size_t iloc = m_rates.install(m_ii, r); /* * Change the reaction rate coefficient type back to its original value @@ -817,7 +824,7 @@ void InterfaceKinetics::addElementaryReaction(ReactionData& r) m_rfn.push_back(r.rateCoeffParameters[0]); registerReaction(reactionNumber(), ELEMENTARY_RXN, iloc); } - +//============================================================================================================================ void InterfaceKinetics::setIOFlag(int ioFlag) { m_ioFlag = ioFlag; diff --git a/src/kinetics/importKinetics.cpp b/src/kinetics/importKinetics.cpp index ed618368f..36026849f 100644 --- a/src/kinetics/importKinetics.cpp +++ b/src/kinetics/importKinetics.cpp @@ -24,6 +24,7 @@ #include "cantera/kinetics/KineticsFactory.h" #include "cantera/kinetics/reaction_defs.h" #include "cantera/kinetics/ReactionData.h" +#include "cantera/kinetics/importKinetics.h" #include "cantera/base/global.h" #include "cantera/base/stringUtils.h" @@ -138,7 +139,7 @@ void checkRxnElementBalance(Kinetics& kin, throw CanteraError("checkRxnElementBalance",msg); } } - +//==================================================================================================================== bool getReagents(const XML_Node& rxn, Kinetics& kin, int rp, std::string default_phase, std::vector& spnum, vector_fp& stoich, vector_fp& order, @@ -165,30 +166,30 @@ bool getReagents(const XML_Node& rxn, Kinetics& kin, int rp, * are stored as a colon separated pair. Get all of these * pairs in the reactions/products object. */ - vector key, val; - getPairs(rg, key, val); + std::vector key, val; + ctml::getPairs(rg, key, val); /* * Loop over each of the pairs and process them */ doublereal ord, stch; - string ph, sp; + string ph, spName; map speciesMap; for (size_t n = 0; n < key.size(); n++) { - sp = key[n]; // sp is the string name for species + spName = key[n]; // sp is the string name for species ph = ""; /* * Search for the species in the kinetics object using the * member function kineticsSpeciesIndex(). We will search * for the species in all phases defined in the kinetics operator. */ - size_t isp = kin.kineticsSpeciesIndex(sp); + size_t isp = kin.kineticsSpeciesIndex(spName); if (isp == npos) { if (rules.skipUndeclaredSpecies) { return false; } else { throw CanteraError("getReagents", - "Undeclared reactant or product species "+sp); + "Undeclared reactant or product species " + spName); return false; } } @@ -209,15 +210,16 @@ bool getReagents(const XML_Node& rxn, Kinetics& kin, int rp, /* * Needed to process reaction orders below. */ - speciesMap[sp] = order.size(); + speciesMap[spName] = order.size(); } /* * Check to see if reaction orders have been specified. */ + if (rp == 1 && rxn.hasChild("order")) { - vector ord; - rxn.getChildren("order",ord); + std::vector ord; + rxn.getChildren("order", ord); doublereal forder; for (size_t nn = 0; nn < ord.size(); nn++) { const XML_Node& oo = *ord[nn]; @@ -240,7 +242,133 @@ bool getReagents(const XML_Node& rxn, Kinetics& kin, int rp, } return true; } +//==================================================================================================================== +// Fill in the fullForwardsOrders array for a specific reaction +/* + * rxnNode XML node for the reaction + */ +bool getOrders(const XML_Node& rxnNode, Kinetics& kin, + std::string default_phase, const ReactionData& rdata, + vector_fp& order, vector_fp& fullForwardsOrders, + const ReactionRules& rules) +{ + // + // Gather the number of species in the kinetics object and resize fullForwardsOrders + // + size_t nsp = kin.nTotalSpecies(); + fullForwardsOrders.resize(nsp, 0.0); + const std::vector& reactants = rdata.reactants; + //const std::vector& rstoich = rdata.rstoich; + const std::vector& products = rdata.products; + const std::vector& pstoich = rdata.pstoich; + + + /* + * Check to see if reaction orders have been specified. + */ + if (rxnNode.hasChild("order")) { + std::vector ord; + rxnNode.getChildren("order", ord); + doublereal forder; + for (size_t nn = 0; nn < ord.size(); nn++) { + const XML_Node& oo = *ord[nn]; + forder = fpValue(oo()); + std::string spName = oo["species"]; + size_t k = kin.kineticsSpeciesIndex(spName); + if (k == npos) { + throw CanteraError("getOrders()", + "Species not in kinetics species list: " + spName); + } + for (size_t n = 0; n < reactants.size(); n++) { + if (reactants[n] == k) { + order[n] = forder; + } + } + } + } + + if (rxnNode.hasChild("orders")) { + std::vector orders; + rxnNode.getChildren("orders", orders); + // + // Doesn't really make sense to have more than one of these blocks + // + if (orders.size() != 1) { + throw CanteraError("getOrders()", " More than one XML orders block"); + } + XML_Node& osNode = *orders[0]; + // + // read the model attribute and figure out how to initialize the full orders vector. + // + string baseHndling = osNode["model"]; + string ss = lowercase(baseHndling); + if (ss == "zeroorders") { + for (size_t k = 0; k < nsp; k++) { + fullForwardsOrders[k] = 0.0; + } + } else if (ss == "reactantorders") { + for (size_t k = 0; k < nsp; k++) { + fullForwardsOrders[k] = 0.0; + } + for (size_t n = 0; n < order.size(); n++) { + size_t k = reactants[n]; + double fac = order[n]; + fullForwardsOrders[k] = fac; + } + } else if (ss == "butlervolmerorders") { + // + // ok first thing to do is get the electrochemical transfer coefficient + // since the order depend on the value. + // Also, if we don't find one, then it's an error + double beta = -10.0; + if (rxnNode.hasChild("rateCoeff")) { + XML_Node& rc = rxnNode.child("rateCoeff"); + if (rc.hasChild("electrochem")) { + XML_Node& eb = rc.child("electrochem"); + string sbeta = eb["beta"]; + beta = fpValueCheck(sbeta); + } + } + if (beta == -10.0) { + throw CanteraError("getOrders()", + "ButlerVolmerOrders model requested but no electrochem beta input"); + } + double betar = 1.0 - beta; + for (size_t k = 0; k < nsp; k++) { + fullForwardsOrders[k] = 0.0; + } + for (size_t n = 0; n < reactants.size(); n++) { + size_t k = reactants[n]; + double fac = order[n]; + fullForwardsOrders[k] += fac * betar; + } + for (size_t n = 0; n < products.size(); n++) { + size_t k = products[n]; + double fac = pstoich[n]; + fullForwardsOrders[k] += fac * beta; + } + } else { + throw CanteraError("getOrders()", "unknown model for orders XML_Node: " + baseHndling); + } + + std::vector key, val; + int numFound = ctml::getPairs(osNode, key, val); + + // + // Fill in the fullForwardsOrders array + // + for (size_t n = 0; n < (size_t) numFound; n++) { + double fac = fpValueCheck(val[n]); + string ss = key[n]; + size_t k = kin.kineticsSpeciesIndex(ss); + fullForwardsOrders[k] = fac; + } + } + + return true; +} +//==================================================================================================================== /** * getArrhenius() parses the xml element called Arrhenius. * The Arrhenius expression is @@ -442,7 +570,7 @@ static void getEfficiencies(const XML_Node& eff, Kinetics& kin, rdata.default_3b_eff = fpValue(eff["default"]); vector key, val; - getPairs(eff, key, val); + ctml::getPairs(eff, key, val); string nm; string phse = kin.thermo(0).id(); for (size_t n = 0; n < key.size(); n++) { @@ -599,37 +727,76 @@ doublereal isDuplicateReaction(std::map& r1, return ratio; } -bool rxninfo::installReaction(int iRxn, const XML_Node& r, Kinetics& kin, +bool rxninfo::installReaction(int iRxn, const XML_Node& rxnNode, Kinetics& kin, string default_phase, ReactionRules& rules, bool validate_rxn) { - // Check to see that we are in fact at a reaction node - if (r.name() != "reaction") { - throw CanteraError(" rxninfo::installReaction", - " expected xml node reaction, got " + r.name()); + // + // Check to see that we are in fact at a reaction node in the XML tree + // + if (rxnNode.name() != "reaction") { + throw CanteraError("rxninfo::installReaction()", + "Expected xml node reaction, got " + rxnNode.name()); } - + // // We use the ReactionData object to store initial values read in from the - // xml data. Then, when we have collected everything we add the reaction to + // xml data. Then, when we have collected everything, we add the reaction to // the kinetics object, kin, at the end of the routine. + // ReactionData& rdata = **m_rdata.insert(m_rdata.end(), new ReactionData()); rdata.validate = validate_rxn; + /* + * Search the reaction element for the attribute "type". + * If found, then branch on the type, to fill in appropriate + * fields in rdata. + */ + rdata.reactionType = ELEMENTARY_RXN; + string typ = rxnNode["type"]; + string ltype = lowercase(typ); + if (typ == "falloff") { + rdata.reactionType = FALLOFF_RXN; + rdata.falloffType = SIMPLE_FALLOFF; + } else if (typ == "chemAct") { + rdata.reactionType = CHEMACT_RXN; + rdata.falloffType = SIMPLE_FALLOFF; + } else if (typ == "threeBody") { + rdata.reactionType = THREE_BODY_RXN; + } else if (typ == "plog") { + rdata.reactionType = PLOG_RXN; + } else if (typ == "chebyshev") { + rdata.reactionType = CHEBYSHEV_RXN; + } else if (typ == "surface") { + rdata.reactionType = SURFACE_RXN; + } else if (typ == "edge") { + rdata.reactionType = EDGE_RXN; + } else if (ltype == "butlervolmer_noactivitycoeffs") { + rdata.reactionType = BUTLERVOLMER_NOACTIVITYCOEFFS_RXN; + } else if (ltype == "butlervolmer") { + rdata.reactionType = BUTLERVOLMER_RXN; + } else if (ltype == "surfaceaffinity") { + rdata.reactionType = SURFACEAFFINITY_RXN; + } else if (ltype == "global") { + rdata.reactionType = GLOBAL_RXN; + } else if (typ != "") { + throw CanteraError("installReaction()", "Unknown reaction type: " + typ); + } + // Check to see if the reaction is specified to be a duplicate of another // reaction. It's an error if the reaction is a duplicate and this is not // set. - rdata.duplicate = (r.hasAttrib("duplicate")) ? 1 : 0; + rdata.duplicate = (rxnNode.hasAttrib("duplicate")) ? 1 : 0; // Check to see if the reaction rate constant can be negative. It's an // error if a negative rate constant is found and this is not set. - rules.allowNegativeA = (r.hasAttrib("negative_A")) ? 1 : 0; + rules.allowNegativeA = (rxnNode.hasAttrib("negative_A")) ? 1 : 0; // Use the contents of the "equation" child element as the reaction's // string representation. Post-process to convert "[" and "]" characters // back into "<" and ">" which cannot easily be stored in an XML file. This // reaction string is used only for display purposes. It is not parsed for // the identities of reactants or products. - rdata.equation = (r.hasChild("equation")) ? r("equation") : ""; + rdata.equation = (rxnNode.hasChild("equation")) ? rxnNode("equation") : ""; static const char* delimiters[] = {" [=] ", " =] ", " = ", "[=]", "=]", "="}; static const char* replacements[] = {" <=> ", " => ", " = ", "<=>", "=>", "="}; for (size_t i = 0; i < 6; i++) { @@ -642,13 +809,14 @@ bool rxninfo::installReaction(int iRxn, const XML_Node& r, Kinetics& kin, break; } } - - // get the reactants - bool ok = getReagents(r, kin, 1, default_phase, rdata.reactants, + // + // get the reactant and their stoichiometries + // + bool ok = getReagents(rxnNode, kin, 1, default_phase, rdata.reactants, rdata.rstoich, rdata.rorder, rules); // Get the products. We store the id of products in rdata.products - ok = ok && getReagents(r, kin, -1, default_phase, rdata.products, + ok = ok && getReagents(rxnNode, kin, -1, default_phase, rdata.products, rdata.pstoich, rdata.porder, rules); // if there was a problem getting either the reactants or the products, @@ -656,26 +824,39 @@ bool rxninfo::installReaction(int iRxn, const XML_Node& r, Kinetics& kin, if (!ok) { return false; } - + // // check whether the reaction is specified to be // reversible. Default is irreversible. - string isrev = r["reversible"]; + // + string isrev = rxnNode["reversible"]; rdata.reversible = (isrev == "yes" || isrev == "true"); + + // HKM this will be removed shortly + // If reaction orders are specified, then this reaction does not follow // mass-action kinetics, and is not an elementary reaction. So check that // it is not reversible, since computing the reverse rate from // thermochemistry only works for elementary reactions. Set the type to // global, so that kinetics managers will know to process the reaction // orders. - if (r.hasChild("order")) { - if (rdata.reversible == true) + if (rxnNode.hasChild("order")) { + if (rdata.reversible == true) { throw CanteraError("installReaction", "reaction orders may only be given for " "irreversible reactions"); + } rdata.global = true; } + // + // Fill in the forwardFullOrder_ array + // + if (rxnNode.hasChild("orders")) { + ok = getOrders(rxnNode, kin, default_phase, rdata, + rdata.rorder, rdata.forwardFullOrder_, rules); + } + // Some reactions can be elementary reactions but have fractional // stoichiometries wrt to some products and reactants. An example of these // are solid reactions involving phase transformations. Species with @@ -712,39 +893,13 @@ bool rxninfo::installReaction(int iRxn, const XML_Node& r, Kinetics& kin, } } - /* - * Search the reaction element for the attribute "type". - * If found, then branch on the type, to fill in appropriate - * fields in rdata. - */ - rdata.reactionType = ELEMENTARY_RXN; - string typ = r["type"]; - if (typ == "falloff") { - rdata.reactionType = FALLOFF_RXN; - rdata.falloffType = SIMPLE_FALLOFF; - } else if (typ == "chemAct") { - rdata.reactionType = CHEMACT_RXN; - rdata.falloffType = SIMPLE_FALLOFF; - } else if (typ == "threeBody") { - rdata.reactionType = THREE_BODY_RXN; - } else if (typ == "plog") { - rdata.reactionType = PLOG_RXN; - } else if (typ == "chebyshev") { - rdata.reactionType = CHEBYSHEV_RXN; - } else if (typ == "surface") { - rdata.reactionType = SURFACE_RXN; - } else if (typ == "edge") { - rdata.reactionType = EDGE_RXN; - } else if (typ != "") { - throw CanteraError("installReaction", "Unknown reaction type: " + typ); - } - + rdata.number = iRxn; rdata.rxn_number = iRxn; // Read the rate coefficient data from the XML file. Trigger an // exception for negative A unless specifically authorized. - getRateCoefficient(r.child("rateCoeff"), kin, rdata, rules); + getRateCoefficient(rxnNode.child("rateCoeff"), kin, rdata, rules); if (validate_rxn) { // Look for undeclared duplicate reactions. @@ -1012,7 +1167,7 @@ bool importKinetics(const XML_Node& phase, std::vector th, } bool buildSolutionFromXML(XML_Node& root, const std::string& id, - const std::string& nm, ThermoPhase* th, Kinetics* k) + const std::string& nm, ThermoPhase* th, Kinetics* kin) { XML_Node* x; x = get_XML_NameID(nm, string("#")+id, &root); @@ -1029,7 +1184,7 @@ bool buildSolutionFromXML(XML_Node& root, const std::string& id, * Create a vector of ThermoPhase pointers of length 1 * having the current th ThermoPhase as the entry. */ - vector phases(1); + std::vector phases(1); phases[0] = th; /* * Fill in the kinetics object k, by querying the @@ -1037,7 +1192,7 @@ bool buildSolutionFromXML(XML_Node& root, const std::string& id, * eventually the source term vector will be constructed * from the list of ThermoPhases in the vector, phases. */ - importKinetics(*x, phases, k); + importKinetics(*x, phases, kin); return true; } diff --git a/src/thermo/GibbsExcessVPSSTP.cpp b/src/thermo/GibbsExcessVPSSTP.cpp index eb3dcaaed..e7edd76e8 100644 --- a/src/thermo/GibbsExcessVPSSTP.cpp +++ b/src/thermo/GibbsExcessVPSSTP.cpp @@ -72,8 +72,7 @@ GibbsExcessVPSSTP& GibbsExcessVPSSTP::operator=(const GibbsExcessVPSSTP& b) return *this; } //========================================================================================================================= -ThermoPhase* -GibbsExcessVPSSTP::duplMyselfAsThermoPhase() const +ThermoPhase* GibbsExcessVPSSTP::duplMyselfAsThermoPhase() const { return new GibbsExcessVPSSTP(*this); } @@ -244,7 +243,7 @@ const vector_fp& GibbsExcessVPSSTP::getPartialMolarVolumesVector() const //========================================================================================================================= double GibbsExcessVPSSTP::checkMFSum(const doublereal* const x) const { - doublereal norm = accumulate(x, x + m_kk, 0.0); + doublereal norm = std::accumulate(x, x + m_kk, 0.0); if (fabs(norm - 1.0) > 1.0E-9) { throw CanteraError("GibbsExcessVPSSTP::checkMFSum", "(MF sum - 1) exceeded tolerance of 1.0E-9:" + fp2str(norm));