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

InterfaceKinetics rewrite -> addition of general BV reactions and affinity formulation.

Browse source
This commit is contained in:
Harry Moffat 2014-08-12 23:02:47 +00:00
parent fcbf41ac73
commit 6741d8f7c6
13 changed files with 449 additions and 107 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
doc/doxygen/Cantera.txt
View file

@ -3,4 +3,7 @@
Use the menu at the top to view detailed documentation of the code.
\ref thermopage
*/

95
doc/doxygen/thermo.txt Normal file
View file

@ -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").
*/

1
include/cantera/base/utilities.h
View file

@ -19,6 +19,7 @@
#include "ct_defs.h"
#include "global.h"
#include <numeric>
#include <algorithm>
namespace Cantera

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

@ -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<size_t> 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<size_t> 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

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

@ -64,8 +64,7 @@ public:
*
* @return Pointer to the new kinetics manager.
*/
virtual Kinetics* newKinetics(XML_Node& phase,
std::vector<ThermoPhase*> th);
virtual Kinetics* newKinetics(XML_Node& phase, std::vector<ThermoPhase*> th);
/**
* Return a new, empty kinetics manager.

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

@ -52,20 +52,40 @@ public:
std::vector<size_t> 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<grouplist_t> rgroups; //!< Optional data used in reaction path diagrams

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

@ -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<size_t, size_t> m_loc;
};

33
include/cantera/kinetics/importKinetics.h
View file

@ -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<size_t>& 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<ThermoPhase*> 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<ThermoPhase*> 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<ThermoPhase*> 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.
/*!

35
include/cantera/kinetics/reaction_defs.h
View file

@ -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;

1
include/cantera/thermo/Phase.h
View file

@ -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"

31
src/kinetics/InterfaceKinetics.cpp
View file

@ -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<bool>(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;

275
src/kinetics/importKinetics.cpp
View file

@ -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<size_t>& 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<string> key, val;
getPairs(rg, key, val);
std::vector<string> 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<string, size_t> 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<XML_Node*> ord;
rxn.getChildren("order",ord);
std::vector<XML_Node*> 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<size_t>& reactants = rdata.reactants;
//const std::vector<doublereal>& rstoich = rdata.rstoich;
const std::vector<size_t>& products = rdata.products;
const std::vector<doublereal>& pstoich = rdata.pstoich;
/*
* Check to see if reaction orders have been specified.
*/
if (rxnNode.hasChild("order")) {
std::vector<XML_Node*> 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<XML_Node*> 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<string> 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<string> 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<int, doublereal>& 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") : "<no equation>";
rdata.equation = (rxnNode.hasChild("equation")) ? rxnNode("equation") : "<no 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<ThermoPhase*> 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<ThermoPhase*> phases(1);
std::vector<ThermoPhase*> 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;
}

5
src/thermo/GibbsExcessVPSSTP.cpp
View file

@ -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));