diff --git a/Cantera/src/KineticsFactory.cpp b/Cantera/src/KineticsFactory.cpp index e09c931f0..9b9ae7b17 100644 --- a/Cantera/src/KineticsFactory.cpp +++ b/Cantera/src/KineticsFactory.cpp @@ -44,7 +44,7 @@ namespace Cantera { * Input * ------ * phaseData = This is an XML_Node that contains the xml data - * describing the phase. Of particular note to ths + * describing the phase. Of particular note to this * routine is the child xml element called "kinetics". * The element has one attribute called "model", * with a string value. The value of this string diff --git a/Cantera/src/ThermoPhase.h b/Cantera/src/ThermoPhase.h index cfdb8a3b7..3692d10d7 100755 --- a/Cantera/src/ThermoPhase.h +++ b/Cantera/src/ThermoPhase.h @@ -39,6 +39,9 @@ namespace Cantera { * * These classes are used to compute thermodynamic properties of * phases of matter. + * + * @see newPhase(std::string file, std::string id) Description for how to read ThermoPhases from XML files. + * @see newPhase(XML_Node &phase) How to call the Factory routine to create and initialize ThermoPhase objects. */ /** diff --git a/Cantera/src/global.h b/Cantera/src/global.h index 00a821bc3..efb23a344 100755 --- a/Cantera/src/global.h +++ b/Cantera/src/global.h @@ -121,6 +121,25 @@ namespace Cantera { * %Cantera was built is added to the search path. * * Additional directories may be added by calling function addDirectory. + * + * There are two different types of input files within %Cantera: + * ctml: This is an xml file layed out in such a way that %Cantera can + * interpret the contents. + * cti: A human-readable ascii format for information that %Cantera + * will read. + * + * %Cantera can take its input from both types of files. However, given + * a file in cti format, the initial operation that %Cantera will perform + * is to translate the cti file into a ctml file. + * The translation is carried out via a system call to a python interpretor + * program that actually carries out the translation. In general, a new + * ctml file is created by the translation that is written to the current + * local directory. + * The ctml file is then read back into %Cantera as the input. + * + * Other input routines in other modules: + * @see importKinetics() + * * @{ */ diff --git a/Cantera/src/importCTML.cpp b/Cantera/src/importCTML.cpp index f96451612..c6506c8b7 100755 --- a/Cantera/src/importCTML.cpp +++ b/Cantera/src/importCTML.cpp @@ -43,427 +43,438 @@ using namespace ctml; using namespace std; -// these are all used to check for duplicate reactions +//! these are all used to check for duplicate reactions class rxninfo { public: - vector< map > rdata; - vector eqn; - vector dup, nr, typ; - vector rev; + //! rdata + std::vector< std::map > rdata; + //! string name + std::vector eqn; + //! string vector of ints + std::vector dup,; + //! string vector of ints + std::vector nr; + //! string vector of ints + std::vector typ; + //! vector of bools. + std::vector rev; }; +//! Temporary storage of rxninfo rxninfo* _rxns = 0; +//! @name utilitydefines. +//@{ #define _reactiondata _rxns->rdata #define _eqn _rxns->eqn #define _dup _rxns->dup #define _nr _rxns->nr #define _typ _rxns->typ #define _rev _rxns->rev +//@} namespace Cantera { - /* - * First we define a couple of typedefs that will - * be used throught this file - */ - typedef const vector nodeset_t; - typedef XML_Node node_t; - - const doublereal DefaultPref = 1.01325e5; // one atm + /* + * First we define a couple of typedefs that will + * be used throught this file + */ + //! typedef for a pointer to an XML_Node + typedef const vector nodeset_t; + //! typedef for an XML_Node + typedef XML_Node node_t; - /// split a string at a '#' sign. Used to separate a file name - /// from an id string. - static void split(const string& src, string& file, string& id) { - string::size_type ipound = src.find('#'); - if (ipound != string::npos) { - id = src.substr(ipound+1,src.size()); - file = src.substr(0,ipound); - } - else { - id = ""; - file = src; - } + /// split a string at a '#' sign. Used to separate a file name + /// from an id string. + static void split(const std::string& src, std::string& file, std::string& id) { + string::size_type ipound = src.find('#'); + if (ipound != string::npos) { + id = src.substr(ipound+1,src.size()); + file = src.substr(0,ipound); } - - - /** - * This routine will locate an XML node in either the input - * XML tree or in another input file specified by the file - * part of the file_ID string. Searches are based on the - * ID attribute of the XML element only. - * - * @param file_ID This is a concatenation of two strings seperated - * by the "#" character. The string before the - * pound character is the file name of an xml - * file to carry out the search. The string after - * the # character is the ID attribute - * of the xml element to search for. - * The string is interpreted as a file string if - * no # character is in the string. - * - * @param root If the file string is empty, searches for the - * xml element with matching ID attribute are - * carried out from this XML node. - */ - XML_Node* get_XML_Node(const string& file_ID, XML_Node* root) { - string fname, idstr; - XML_Node *db, *doc; - split(file_ID, fname, idstr); - if (fname == "") { - if (!root) throw CanteraError("get_XML_Node", - "no file name given. file_ID = "+file_ID); - db = root->findID(idstr, 3); - } - else { - doc = get_XML_File(fname); - if (!doc) throw CanteraError("get_XML_Node", - "get_XML_File failed trying to open "+fname); - db = doc->findID(idstr, 3); - } - if (!db) { - throw CanteraError("get_XML_Node", - "id tag '"+idstr+"' not found."); - } - return db; + else { + id = ""; + file = src; } + } - /** - * This routine will locate an XML node in either the input - * XML tree or in another input file specified by the file - * part of the file_ID string. Searches are based on the - * XML element name and the ID attribute of the XML element. - * An exact match of both is usually required. However, the - * ID attribute may be set to "", in which case the first - * xml element with the correct element name will be returned. - * - * @param nameTarget This is the XML element name to look for. - * - * @param file_ID This is a concatenation of two strings seperated - * by the "#" character. The string before the - * pound character is the file name of an xml - * file to carry out the search. The string after - * the # character is the ID attribute - * of the xml element to search for. - * The string is interpreted as a file string if - * no # character is in the string. - * - * @param root If the file string is empty, searches for the - * xml element with matching ID attribute are - * carried out from this XML node. - */ - XML_Node* get_XML_NameID(const string& nameTarget, - const string& file_ID, - XML_Node* root) { - string fname, idTarget; - XML_Node *db, *doc; - split(file_ID, fname, idTarget); - if (fname == "") { - if (!root) return 0; - db = root->findNameID(nameTarget, idTarget); - } else { - doc = get_XML_File(fname); - if (!doc) return 0; - db = doc->findNameID(nameTarget, idTarget); - } - return db; + /* + * This routine will locate an XML node in either the input + * XML tree or in another input file specified by the file + * part of the file_ID string. Searches are based on the + * ID attribute of the XML element only. + * + * @param file_ID This is a concatenation of two strings seperated + * by the "#" character. The string before the + * pound character is the file name of an xml + * file to carry out the search. The string after + * the # character is the ID attribute + * of the xml element to search for. + * The string is interpreted as a file string if + * no # character is in the string. + * + * @param root If the file string is empty, searches for the + * xml element with matching ID attribute are + * carried out from this XML node. + */ + XML_Node* get_XML_Node(const std::string& file_ID, XML_Node* root) { + string fname, idstr; + XML_Node *db, *doc; + split(file_ID, fname, idstr); + if (fname == "") { + if (!root) throw CanteraError("get_XML_Node", + "no file name given. file_ID = "+file_ID); + db = root->findID(idstr, 3); + } + else { + doc = get_XML_File(fname); + if (!doc) throw CanteraError("get_XML_Node", + "get_XML_File failed trying to open "+fname); + db = doc->findID(idstr, 3); } + if (!db) { + throw CanteraError("get_XML_Node", + "id tag '"+idstr+"' not found."); + } + return db; + } + + /* + * This routine will locate an XML node in either the input + * XML tree or in another input file specified by the file + * part of the file_ID string. Searches are based on the + * XML element name and the ID attribute of the XML element. + * An exact match of both is usually required. However, the + * ID attribute may be set to "", in which case the first + * xml element with the correct element name will be returned. + * + * @param nameTarget This is the XML element name to look for. + * + * @param file_ID This is a concatenation of two strings seperated + * by the "#" character. The string before the + * pound character is the file name of an xml + * file to carry out the search. The string after + * the # character is the ID attribute + * of the xml element to search for. + * The string is interpreted as a file string if + * no # character is in the string. + * + * @param root If the file string is empty, searches for the + * xml element with matching ID attribute are + * carried out from this XML node. + */ + XML_Node* get_XML_NameID(const std::string& nameTarget, + const std::string& file_ID, + XML_Node* root) { + string fname, idTarget; + XML_Node *db, *doc; + split(file_ID, fname, idTarget); + if (fname == "") { + if (!root) return 0; + db = root->findNameID(nameTarget, idTarget); + } else { + doc = get_XML_File(fname); + if (!doc) return 0; + db = doc->findNameID(nameTarget, idTarget); + } + return db; + } - /** - * Install a species into a ThermoPhase object, which defines - * the phase thermodynamics and speciation. - * - * This routine first gathers the information from the Species XML - * tree and calls addUniqueSpecies() to add it to the - * ThermoPhase object, p. - * This information consists of: - * ecomp[] = element composition of species. - * chgr = electric charge of species - * name = string name of species - * sz = size of the species - * (option double used a lot in thermo) - * - * Then, the routine processes the "thermo" XML element and - * calls underlying utility routines to read the XML elements - * containing the thermodynamic information for the reference - * state of the species. Failures or lack of information trigger - * an "UnknownSpeciesThermoModel" exception being thrown. - */ - bool installSpecies(int k, const XML_Node& s, thermo_t& p, - SpeciesThermo& spthermo, int rule, SpeciesThermoFactory* factory) { + /* + * Install a species into a ThermoPhase object, which defines + * the phase thermodynamics and speciation. + * + * This routine first gathers the information from the Species XML + * tree and calls addUniqueSpecies() to add it to the + * ThermoPhase object, p. + * This information consists of: + * ecomp[] = element composition of species. + * chgr = electric charge of species + * name = string name of species + * sz = size of the species + * (option double used a lot in thermo) + * + * Then, the routine processes the "thermo" XML element and + * calls underlying utility routines to read the XML elements + * containing the thermodynamic information for the reference + * state of the species. Failures or lack of information trigger + * an "UnknownSpeciesThermoModel" exception being thrown. + */ + bool installSpecies(int k, const XML_Node& s, thermo_t& p, + SpeciesThermo& spthermo, int rule, + SpeciesThermoFactory* factory) { - // get the composition of the species - const XML_Node& a = s.child("atomArray"); - map comp; - getMap(a, comp); + // get the composition of the species + const XML_Node& a = s.child("atomArray"); + map comp; + getMap(a, comp); - // check that all elements in the species - // exist in 'p'. If rule != 0, quietly skip - // this species if some elements are undeclared; - // otherwise, throw an exception - map::const_iterator _b = comp.begin(); - for (; _b != comp.end(); ++_b) { - if (p.elementIndex(_b->first) < 0) { - if (rule == 0) { - throw CanteraError("installSpecies", - "Species " + s["name"] + - " contains undeclared element " + _b->first); - } - else - return false; - } + // check that all elements in the species + // exist in 'p'. If rule != 0, quietly skip + // this species if some elements are undeclared; + // otherwise, throw an exception + map::const_iterator _b = comp.begin(); + for (; _b != comp.end(); ++_b) { + if (p.elementIndex(_b->first) < 0) { + if (rule == 0) { + throw CanteraError("installSpecies", + "Species " + s["name"] + + " contains undeclared element " + _b->first); } + else + return false; + } + } - // construct a vector of atom numbers for each - // element in phase p. Elements not declared in the - // species (i.e., not in map comp) will have zero - // entries in the vector. - int m, nel = p.nElements(); - vector_fp ecomp(nel, 0.0); - for (m = 0; m < nel; m++) { - ecomp[m] = atoi(comp[p.elementName(m)].c_str()); - } + // construct a vector of atom numbers for each + // element in phase p. Elements not declared in the + // species (i.e., not in map comp) will have zero + // entries in the vector. + int m, nel = p.nElements(); + vector_fp ecomp(nel, 0.0); + for (m = 0; m < nel; m++) { + ecomp[m] = atoi(comp[p.elementName(m)].c_str()); + } - // get the species charge, if any. Note that the charge need - // not be explicitly specified if special element 'E' - // (electron) is one of the elements. - doublereal chrg = 0.0; - if (s.hasChild("charge")) chrg = getFloat(s, "charge"); + // get the species charge, if any. Note that the charge need + // not be explicitly specified if special element 'E' + // (electron) is one of the elements. + doublereal chrg = 0.0; + if (s.hasChild("charge")) chrg = getFloat(s, "charge"); - // get the species size, if any. (This is used by surface - // phases to represent how many sites a species occupies.) - doublereal sz = 1.0; - if (s.hasChild("size")) sz = getFloat(s, "size"); + // get the species size, if any. (This is used by surface + // phases to represent how many sites a species occupies.) + doublereal sz = 1.0; + if (s.hasChild("size")) sz = getFloat(s, "size"); - // add the species to phase p. - p.addUniqueSpecies(s["name"], &ecomp[0], chrg, sz); + // add the species to phase p. + p.addUniqueSpecies(s["name"], &ecomp[0], chrg, sz); - // install the thermo parameterization for this species into - // the species thermo manager for phase p. - factory->installThermoForSpecies(k, s, spthermo); + // install the thermo parameterization for this species into + // the species thermo manager for phase p. + factory->installThermoForSpecies(k, s, spthermo); - return true; + return true; + } + + + + /* + * Check a reaction to see if the elements balance. + */ + void checkRxnElementBalance(Kinetics& kin, + const ReactionData &rdata, doublereal errorTolerance) { + int index, klocal, n, kp, kr, m, nel; + double kstoich; + + map bal, balr, balp; + bal.clear(); + balp.clear(); + balr.clear(); + + int np = rdata.products.size(); + + // iterate over the products + for (index = 0; index < np; index++) { + kp = rdata.products[index]; // index of the product in 'kin' + n = kin.speciesPhaseIndex(kp); // phase this product belongs to + klocal = kp - kin.kineticsSpeciesIndex(0,n); // index within this phase + kstoich = rdata.pstoich[index]; // product stoichiometric coeff + const ThermoPhase& ph = kin.speciesPhase(kp); + nel = ph.nElements(); + for (m = 0; m < nel; m++) { + bal[ph.elementName(m)] += kstoich*ph.nAtoms(klocal,m); + balp[ph.elementName(m)] += kstoich*ph.nAtoms(klocal,m); + } + } + int nr = rdata.reactants.size(); + for (index = 0; index < nr; index++) { + kr = rdata.reactants[index]; + n = kin.speciesPhaseIndex(kr); + //klocal = kr - kin.start(n); + klocal = kr - kin.kineticsSpeciesIndex(0,n); + kstoich = rdata.rstoich[index]; + const ThermoPhase& ph = kin.speciesPhase(kr); + nel = ph.nElements(); + for (m = 0; m < nel; m++) { + bal[ph.elementName(m)] -= kstoich*ph.nAtoms(klocal,m); + balr[ph.elementName(m)] += kstoich*ph.nAtoms(klocal,m); + } } + map::iterator b = bal.begin(); + string msg = "\n\tElement Reactants Products"; + bool ok = true; + doublereal err; + for (; b != bal.end(); ++b) { + err = fabs(b->second/(balr[b->first] + balp[b->first])); + if (err > errorTolerance) { + ok = false; + msg += "\n\t"+b->first+" "+ fp2str(balr[b->first]) + +" "+ fp2str(balp[b->first]); + } + } + if (!ok) { + msg = "The following reaction is unbalanced:\n\t" + + rdata.equation + "\n" + msg + "\n"; + throw CanteraError("checkRxnElementBalance",msg); + } + } - /** - * Check a reaction to see if the elements balance. + /** + * Get the reactants or products of a reaction. The information + * is returned in the spnum, stoich, and order vectors. The + * length of the vectors is the number of different types of + * reactants or products found for the reaction. + * + * Input + * -------- + * rxn -> xml node pointing to the reaction element + * in the xml tree. + * kin -> Reference to the kinetics object to install + * the information into. + * rp = 1 -> Go get the reactants for a reaction + * -1 -> Go get the products for a reaction + * default_phase = String name for the default phase + * to loop up species in. + * Output + * ----------- + * spnum = vector of species numbers found. + * Length is number of reactants or products. + * stoich = stoichiometric coefficient of the reactant or product + * Length is number of reactants or products. + * order = Order of the reactant and product in the reaction + * rate expression + * rule = If we fail to find a species, we will throw an error + * if rule != 1. If rule = 1, we simply return false, + * allowing the calling routine to skip this reaction + * and continue. + */ + static bool getReagents(const XML_Node& rxn, kinetics_t& kin, int rp, + string default_phase, + vector_int& spnum, vector_fp& stoich, vector_fp& order, + int rule) { + + string rptype; + + /* + * The id of reactants and products are kept in child elements + * of reaction, named "reactants" and "products". We search + * the xml tree for these children based on the value of rp, + * and store the xml element pointer here. */ - void checkRxnElementBalance(Kinetics& kin, - const ReactionData &rdata, doublereal errorTolerance) { - int index, klocal, n, kp, kr, m, nel; - double kstoich; + if (rp == 1) rptype = "reactants"; + else rptype = "products"; + const XML_Node& rg = rxn.child(rptype); - map bal, balr, balp; - bal.clear(); - balp.clear(); - balr.clear(); - - int np = rdata.products.size(); - - // iterate over the products - for (index = 0; index < np; index++) { - kp = rdata.products[index]; // index of the product in 'kin' - n = kin.speciesPhaseIndex(kp); // phase this product belongs to - klocal = kp - kin.kineticsSpeciesIndex(0,n); // index within this phase - kstoich = rdata.pstoich[index]; // product stoichiometric coeff - const ThermoPhase& ph = kin.speciesPhase(kp); - nel = ph.nElements(); - for (m = 0; m < nel; m++) { - bal[ph.elementName(m)] += kstoich*ph.nAtoms(klocal,m); - balp[ph.elementName(m)] += kstoich*ph.nAtoms(klocal,m); - } - } - int nr = rdata.reactants.size(); - for (index = 0; index < nr; index++) { - kr = rdata.reactants[index]; - n = kin.speciesPhaseIndex(kr); - //klocal = kr - kin.start(n); - klocal = kr - kin.kineticsSpeciesIndex(0,n); - kstoich = rdata.rstoich[index]; - const ThermoPhase& ph = kin.speciesPhase(kr); - nel = ph.nElements(); - for (m = 0; m < nel; m++) { - bal[ph.elementName(m)] -= kstoich*ph.nAtoms(klocal,m); - balr[ph.elementName(m)] += kstoich*ph.nAtoms(klocal,m); - } - } - - map::iterator b = bal.begin(); - string msg = "\n\tElement Reactants Products"; - bool ok = true; - doublereal err; - for (; b != bal.end(); ++b) { - err = fabs(b->second/(balr[b->first] + balp[b->first])); - if (err > errorTolerance) { - ok = false; - msg += "\n\t"+b->first+" "+ fp2str(balr[b->first]) - +" "+ fp2str(balp[b->first]); - } - } - if (!ok) { - msg = "The following reaction is unbalanced:\n\t" - + rdata.equation + "\n" + msg + "\n"; - throw CanteraError("checkRxnElementBalance",msg); - } - } - - - - /** - * Get the reactants or products of a reaction. The information - * is returned in the spnum, stoich, and order vectors. The - * length of the vectors is the number of different types of - * reactants or products found for the reaction. - * - * Input - * -------- - * rxn -> xml node pointing to the reaction element - * in the xml tree. - * kin -> Reference to the kinetics object to install - * the information into. - * rp = 1 -> Go get the reactants for a reaction - * -1 -> Go get the products for a reaction - * default_phase = String name for the default phase - * to loop up species in. - * Output - * ----------- - * spnum = vector of species numbers found. - * Length is number of reactants or products. - * stoich = stoichiometric coefficient of the reactant or product - * Length is number of reactants or products. - * order = Order of the reactant and product in the reaction - * rate expression - * rule = If we fail to find a species, we will throw an error - * if rule != 1. If rule = 1, we simply return false, - * allowing the calling routine to skip this reaction - * and continue. + /* + * The species and stoichiometric coefficient for the species + * are stored as a colon seperated pair. Get all of these + * pairs in the reactions/products object. */ - static bool getReagents(const XML_Node& rxn, kinetics_t& kin, int rp, - string default_phase, - vector_int& spnum, vector_fp& stoich, vector_fp& order, - int rule) { + vector key, val; + getPairs(rg, key, val); - string rptype; + int ns = static_cast(key.size()); - /* - * The id of reactants and products are kept in child elements - * of reaction, named "reactants" and "products". We search - * the xml tree for these children based on the value of rp, - * and store the xml element pointer here. - */ - if (rp == 1) rptype = "reactants"; - else rptype = "products"; - const XML_Node& rg = rxn.child(rptype); + /* + * Loop over each of the pairs and process them + */ + int isp; + doublereal ord, stch; + string ph, sp; + map speciesMap; + for (int n = 0; n < ns; n++) { + sp = 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. + */ + isp = kin.kineticsSpeciesIndex(sp,""); + if (isp < 0) { + if (rule == 1) + return false; + else { + throw CanteraError("getReagents", + "Undeclared reactant or product species "+sp); + return false; + } + } - /* - * The species and stoichiometric coefficient for the species - * are stored as a colon seperated pair. Get all of these - * pairs in the reactions/products object. - */ - vector key, val; - getPairs(rg, key, val); - - int ns = static_cast(key.size()); - - /* - * Loop over each of the pairs and process them - */ - int isp; - doublereal ord, stch; - string ph, sp; - map speciesMap; - for (int n = 0; n < ns; n++) { - sp = 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. - */ - isp = kin.kineticsSpeciesIndex(sp,""); - if (isp < 0) { - if (rule == 1) - return false; - else { - throw CanteraError("getReagents", - "Undeclared reactant or product species "+sp); - return false; - } - } - - /* - * For each reagent, we store the the species number, isp - * the stoichiometric coefficient, val[n], and the order - * species in the reaction rate expression. We assume mass - * action kinetics here, but will modify this below for - * specified species. - */ - spnum.push_back(isp); - stch = atof(val[n].c_str()); - stoich.push_back(stch); - ord = doublereal(stch); - order.push_back(ord); - /* - * Needed to process reaction orders below. - */ - speciesMap[sp] = order.size(); - } - - /* - * Check to see if reactant reaction orders have been specified. - */ - if (rp == 1 && rxn.hasChild("order")) { - vector ord; - rxn.getChildren("order",ord); - int norder = static_cast(ord.size()); - int loc; - doublereal forder; - for (int nn = 0; nn < norder; nn++) { - const XML_Node& oo = *ord[nn]; - string sp = oo["species"]; - loc = speciesMap[sp]; - if (loc == 0) - throw CanteraError("getReagents", - "reaction order specified for non-reactant: " - +sp); - forder = fpValue(oo()); - if (forder < 0.0) { - throw CanteraError("getReagents", - "reaction order must be non-negative"); - } - // replace the forward stoichiometric coefficient - // stored above in 'order' with the specified - // reaction order - order[loc-1] = forder; - } - } - return true; + /* + * For each reagent, we store the the species number, isp + * the stoichiometric coefficient, val[n], and the order + * species in the reaction rate expression. We assume mass + * action kinetics here, but will modify this below for + * specified species. + */ + spnum.push_back(isp); + stch = atof(val[n].c_str()); + stoich.push_back(stch); + ord = doublereal(stch); + order.push_back(ord); + /* + * Needed to process reaction orders below. + */ + speciesMap[sp] = order.size(); } + + /* + * Check to see if reactant reaction orders have been specified. + */ + if (rp == 1 && rxn.hasChild("order")) { + vector ord; + rxn.getChildren("order",ord); + int norder = static_cast(ord.size()); + int loc; + doublereal forder; + for (int nn = 0; nn < norder; nn++) { + const XML_Node& oo = *ord[nn]; + string sp = oo["species"]; + loc = speciesMap[sp]; + if (loc == 0) + throw CanteraError("getReagents", + "reaction order specified for non-reactant: " + +sp); + forder = fpValue(oo()); + if (forder < 0.0) { + throw CanteraError("getReagents", + "reaction order must be non-negative"); + } + // replace the forward stoichiometric coefficient + // stored above in 'order' with the specified + // reaction order + order[loc-1] = forder; + } + } + return true; + } - /** - * getArrhenius() parses the xml element called Arrhenius. - * The Arrhenius expression is - * \f[ k = A T^(b) exp (-E_a / RT). \f] - */ - static void getArrhenius(const XML_Node& node, int& highlow, - doublereal& A, doublereal& b, doublereal& E) { + /** + * getArrhenius() parses the xml element called Arrhenius. + * The Arrhenius expression is + * \f[ k = A T^(b) exp (-E_a / RT). \f] + */ + static void getArrhenius(const XML_Node& node, int& highlow, + doublereal& A, doublereal& b, doublereal& E) { - if (node["name"] == "k0") - highlow = 0; - else highlow = 1; - /* - * We parse the children for the A, b, and E conponents. - */ - A = getFloat(node, "A", "-"); - b = getFloat(node, "b"); - E = getFloat(node, "E", "actEnergy"); - E /= GasConstant; - } + if (node["name"] == "k0") + highlow = 0; + else highlow = 1; + /* + * We parse the children for the A, b, and E conponents. + */ + A = getFloat(node, "A", "-"); + b = getFloat(node, "b"); + E = getFloat(node, "E", "actEnergy"); + E /= GasConstant; + } /** * getStick() processes the XML element called Stick that specifies @@ -600,915 +611,911 @@ namespace Cantera { rdata.falloffParameters = c; } - /** - * Get the enhanced collision efficiencies. It is assumed that the - * reaction mechanism is homogeneous, so that all species belong - * to phase(0) of 'kin'. + /** + * Get the enhanced collision efficiencies. It is assumed that the + * reaction mechanism is homogeneous, so that all species belong + * to phase(0) of 'kin'. + */ + static void getEfficiencies(const node_t& eff, kinetics_t& kin, ReactionData& rdata) { + + // set the default collision efficiency + rdata.default_3b_eff = fpValue(eff["default"]); + + vector key, val; + getPairs(eff, key, val); + int ne = static_cast(key.size()); + string nm; + string phse = kin.thermo(0).id(); + int n, k; + for (n = 0; n < ne; n++) { // ; bb != ee; ++bb) { + nm = key[n];// bb->first; + k = kin.kineticsSpeciesIndex(nm, phse); + rdata.thirdBodyEfficiencies[k] = fpValue(val[n]); // bb->second; + } + } + + /* + * Extract the rate coefficient for a reaction from the xml node, kf. + * kf should point to a XML element named "rateCoeff". + * rdata is the partially filled ReactionData object for the reaction. + * This function will fill in more fields in the ReactionData object. + */ + void getRateCoefficient(const node_t& kf, kinetics_t& kin, + ReactionData& rdata, int negA) { + + int nc = kf.nChildren(); + nodeset_t& kf_children = kf.children(); + vector_fp clow(3,0.0), chigh(3,0.0); + // int nr = nReacMolecules(rdata); + for (int m = 0; m < nc; m++) { + const node_t& c = *kf_children[m]; + string nm = c.name(); + int highlow=0; + + if (nm == "Arrhenius") { + vector_fp coeff(3); + if (c["type"] == "stick") { + getStick(c, kin, rdata, coeff[0], coeff[1], coeff[2]); + chigh = coeff; + } + else { + getArrhenius(c, highlow, coeff[0], coeff[1], coeff[2]); + if (highlow == 1 || rdata.reactionType == THREE_BODY_RXN + || rdata.reactionType == ELEMENTARY_RXN) + chigh = coeff; + else clow = coeff; + } + if (rdata.reactionType == SURFACE_RXN) { + getCoverageDependence(c, + kin.thermo(kin.surfacePhaseIndex()), rdata); + } + + if (coeff[0] <= 0.0 && negA == 0) { + throw CanteraError("getRateCoefficient", + "negative or zero A coefficient for reaction "+int2str(rdata.number)); + } + } + else if (nm == "falloff") { + getFalloff(c, rdata); + } + else if (nm == "efficiencies") { + getEfficiencies(c, kin, rdata); + } + else if (nm == "electrochem") { + rdata.beta = fpValue(c["beta"]); + } + } + /* + * Store the coefficients in the ReactionData object for return + * from this function. */ - static void getEfficiencies(const node_t& eff, kinetics_t& kin, ReactionData& rdata) { + if (rdata.reactionType == CHEMACT_RXN) + rdata.rateCoeffParameters = clow; + else + rdata.rateCoeffParameters = chigh; - // set the default collision efficiency - rdata.default_3b_eff = fpValue(eff["default"]); - - vector key, val; - getPairs(eff, key, val); - int ne = static_cast(key.size()); - string nm; - string phse = kin.thermo(0).id(); - int n, k; - for (n = 0; n < ne; n++) { // ; bb != ee; ++bb) { - nm = key[n];// bb->first; - k = kin.kineticsSpeciesIndex(nm, phse); - rdata.thirdBodyEfficiencies[k] = fpValue(val[n]); // bb->second; - } - } - - /** - * Extract the rate coefficient for a reaction from the xml node, kf. - * kf should point to a XML element named "rateCoeff". - * rdata is the partially filled ReactionData object for the reaction. - * This function will fill in more fields in the ReactionData object. - * - */ - void getRateCoefficient(const node_t& kf, kinetics_t& kin, - ReactionData& rdata, int negA) { - - int nc = kf.nChildren(); - nodeset_t& kf_children = kf.children(); - vector_fp clow(3,0.0), chigh(3,0.0); - // int nr = nReacMolecules(rdata); - for (int m = 0; m < nc; m++) { - const node_t& c = *kf_children[m]; - string nm = c.name(); - int highlow=0; - - if (nm == "Arrhenius") { - vector_fp coeff(3); - if (c["type"] == "stick") { - getStick(c, kin, rdata, coeff[0], coeff[1], coeff[2]); - chigh = coeff; - } - else { - getArrhenius(c, highlow, coeff[0], coeff[1], coeff[2]); - if (highlow == 1 || rdata.reactionType == THREE_BODY_RXN - || rdata.reactionType == ELEMENTARY_RXN) - chigh = coeff; - else clow = coeff; - } - if (rdata.reactionType == SURFACE_RXN) { - getCoverageDependence(c, - kin.thermo(kin.surfacePhaseIndex()), rdata); - } - - if (coeff[0] <= 0.0 && negA == 0) { - throw CanteraError("getRateCoefficient", - "negative or zero A coefficient for reaction "+int2str(rdata.number)); - } - } - else if (nm == "falloff") { - getFalloff(c, rdata); - } - else if (nm == "efficiencies") { - getEfficiencies(c, kin, rdata); - } - else if (nm == "electrochem") { - rdata.beta = fpValue(c["beta"]); - } - } - /* - * Store the coefficients in the ReactionData object for return - * from this function. - */ - if (rdata.reactionType == CHEMACT_RXN) - rdata.rateCoeffParameters = clow; - else - rdata.rateCoeffParameters = chigh; - - if (rdata.reactionType == FALLOFF_RXN) - rdata.auxRateCoeffParameters = clow; - else if (rdata.reactionType == CHEMACT_RXN) - rdata.auxRateCoeffParameters = chigh; - } + if (rdata.reactionType == FALLOFF_RXN) + rdata.auxRateCoeffParameters = clow; + else if (rdata.reactionType == CHEMACT_RXN) + rdata.auxRateCoeffParameters = chigh; + } - /** - * Create a new ThermoPhase object and initializes it according to - * the XML tree database. This routine first looks up the - * identity of the model for the solution thermodynamics in the - * model attribute of the thermo child of the xml phase - * node. Then, it does a string lookup on the model to figure out - * what ThermoPhase derived class is assigned. It creates a new - * instance of that class, and then calls importPhase() to - * populate that class with the correct parameters from the XML - * tree. - */ - ThermoPhase* newPhase(XML_Node& xmlphase) { - const XML_Node& th = xmlphase.child("thermo"); - string model = th["model"]; - ThermoPhase* t = newThermoPhase(model); - importPhase(xmlphase, t); - return t; - } + /* + * Create a new ThermoPhase object and initializes it according to + * the XML tree database. This routine first looks up the + * identity of the model for the solution thermodynamics in the + * model attribute of the thermo child of the xml phase + * node. Then, it does a string lookup on the model to figure out + * what ThermoPhase derived class is assigned. It creates a new + * instance of that class, and then calls importPhase() to + * populate that class with the correct parameters from the XML + * tree. + */ + ThermoPhase* newPhase(XML_Node& xmlphase) { + const XML_Node& th = xmlphase.child("thermo"); + string model = th["model"]; + ThermoPhase* t = newThermoPhase(model); + importPhase(xmlphase, t); + return t; + } - ThermoPhase* newPhase(string infile, string id) { - XML_Node* root = get_XML_File(infile); - if (id == "-") id = ""; - XML_Node* x = get_XML_Node(string("#")+id, root); - if (x) - return newPhase(*x); - else - return 0; - } + ThermoPhase* newPhase(std::string infile, std::string id) { + XML_Node* root = get_XML_File(infile); + if (id == "-") id = ""; + XML_Node* x = get_XML_Node(string("#")+id, root); + if (x) + return newPhase(*x); + else + return 0; + } - /** - * Import a phase specification. - * Here we read an XML description of the phase. - * We import descriptions of the elements that make up the - * species in a phase. - * We import information about the species, including their - * reference state thermodynamic polynomials. We then freeze - * the state of the species, and finally call initThermo() - * a member function of the ThermoPhase object to "finish" - * the description. - * - * - * @param phase This object must be the phase node of a - * complete XML tree - * description of the phase, including all of the - * species data. In other words while "phase" must - * point to an XML phase object, it must have - * sibling nodes "speciesData" that describe - * the species in the phase. - * @param th Pointer to the ThermoPhase object which will - * handle the thermodynamics for this phase. - * We initialize part of the Thermophase object - * here, especially for those objects which are - * part of the Cantera Kernel. - */ - bool importPhase(XML_Node& phase, ThermoPhase* th, - SpeciesThermoFactory* spfactory) { + /* + * Import a phase specification. + * Here we read an XML description of the phase. + * We import descriptions of the elements that make up the + * species in a phase. + * We import information about the species, including their + * reference state thermodynamic polynomials. We then freeze + * the state of the species, and finally call initThermo() + * a member function of the ThermoPhase object to "finish" + * the description. + * + * + * @param phase This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * @param th Pointer to the ThermoPhase object which will + * handle the thermodynamics for this phase. + * We initialize part of the Thermophase object + * here, especially for those objects which are + * part of the Cantera Kernel. + */ + bool importPhase(XML_Node& phase, ThermoPhase* th, + SpeciesThermoFactory* spfactory) { - // Check the the supplied XML node in fact represents a - // phase. - if (phase.name() != "phase") - throw CanteraError("importPhase", - "Current const XML_Node is not a phase element."); + // Check the the supplied XML node in fact represents a + // phase. + if (phase.name() != "phase") + throw CanteraError("importPhase", + "Current const XML_Node is not a phase element."); - // if no species thermo factory was supplied, - // use the default one. - if (!spfactory) - spfactory = SpeciesThermoFactory::factory(); + // if no species thermo factory was supplied, + // use the default one. + if (!spfactory) + spfactory = SpeciesThermoFactory::factory(); - // set the id attribute of the phase to the 'id' attribute - // in the XML tree. - th->setID(phase.id()); - th->setName(phase.id()); + // set the id attribute of the phase to the 'id' attribute + // in the XML tree. + th->setID(phase.id()); + th->setName(phase.id()); - // Number of spatial dimensions. Defaults to 3 (bulk phase) - if (phase.hasAttrib("dim")) { - int idim = intValue(phase["dim"]); - if (idim < 1 || idim > 3) - throw CanteraError("importPhase", - "unphysical number of dimensions: "+phase["dim"]); - th->setNDim(idim); - } - else - th->setNDim(3); // default + // Number of spatial dimensions. Defaults to 3 (bulk phase) + if (phase.hasAttrib("dim")) { + int idim = intValue(phase["dim"]); + if (idim < 1 || idim > 3) + throw CanteraError("importPhase", + "unphysical number of dimensions: "+phase["dim"]); + th->setNDim(idim); + } + else + th->setNDim(3); // default - // Set equation of state parameters. The parameters are - // specific to each subclass of ThermoPhase, so this is done - // by method setParametersFromXML in each subclass. - if (phase.hasChild("thermo")) { - const XML_Node& eos = phase.child("thermo"); - th->setParametersFromXML(eos); - } - - - /*************************************************************** - * Add the elements. - ***************************************************************/ - th->addElementsFromXML(phase); - - - /*************************************************************** - * Add the species. - * - * Species definitions may be imported from multiple - * sources. For each one, a speciesArray element must be - * present. - ***************************************************************/ - XML_Node* db = 0; - vector sparrays; - phase.getChildren("speciesArray", sparrays); - int jsp, nspa = static_cast(sparrays.size()); - vector dbases; - vector_int sprule(nspa,0); - - // loop over the speciesArray elements - for (jsp = 0; jsp < nspa; jsp++) { - - const XML_Node& species = *sparrays[jsp]; - - // If the speciesArray element has a child element - // - // then set sprule[jsp] to 1, so - // that any species with an undeclared element will be - // quietly skipped when importing species. - if (species.hasChild("skip")) { - const XML_Node& sk = species.child("skip"); - string eskip = sk["element"]; - if (eskip == "undeclared") { - sprule[jsp] = 1; - } - string dskip = sk["species"]; - if (dskip == "duplicate") { - sprule[jsp] += 10; - } - } - - string fname, idstr; - - // get a pointer to the node containing the species - // definitions for the species declared in this - // speciesArray element. This may be in the local file - // containing the phase element, or may be in another - // file. - db = get_XML_Node(species["datasrc"], &phase.root()); - - // add this node to the list of species database nodes. - dbases.push_back(db); - } - - - // if the phase has a species thermo manager already installed, - // delete it since we are adding new species. - delete &th->speciesThermo(); - - // create a new species thermo manager. Function - // 'newSpeciesThermoMgr' looks at the species in the database - // to see what thermodynamic property parameterizations are - // used, and selects a class that can handle the - // parameterizations found. - SpeciesThermo* spth = newSpeciesThermoMgr(dbases); - - // install it in the phase object - th->setSpeciesThermo(spth); - SpeciesThermo& spthermo = th->speciesThermo(); - - // used to check that each species is declared only once - map declared; - - int i, k = 0; - - // loop over the species arrays - for (jsp = 0; jsp < nspa; jsp++) { - - const XML_Node& species = *sparrays[jsp]; - db = dbases[jsp]; - - // Get the array of species name strings. - vector spnames; - getStringArray(species, spnames); - int nsp = static_cast(spnames.size()); - - // if 'all' is specified, then add all species - // defined in this database to the phase - if (nsp == 1 && spnames[0] == "all") { - vector allsp; - db->getChildren("species",allsp); - nsp = static_cast(allsp.size()); - spnames.resize(nsp); - for (int nn = 0; nn < nsp; nn++) { - spnames[nn] = (*allsp[nn])["name"]; - } - } - else if (nsp == 1 && spnames[0] == "unique") { - vector uniquesp; - db->getChildren("species",uniquesp); - nsp = static_cast(uniquesp.size()); - spnames.clear(); - spnames.resize(nsp); - string spnm; - for (int nn = 0; nn < nsp; nn++) { - spnm = (*uniquesp[nn])["name"]; - if (!declared[spnm]) spnames[nn] = spnm; - } - } - - string name; - bool skip; - for (i = 0; i < nsp; i++) { - name = spnames[i]; - skip = false; - if (name == "") skip = true; - // Check that every species is only declared once - if (declared[name]) { - if (sprule[jsp] >= 10) - skip = true; - else - throw CanteraError("importPhase", - "duplicate species: "+name); - } - if (!skip) { - declared[name] = true; - - // Find the species in the database by name. - XML_Node* s = db->findByAttr("name",spnames[i]); - if (s) { - if (installSpecies(k, *s, *th, spthermo, sprule[jsp], - spfactory)) - ++k; - } - else { - throw CanteraError("importPhase","no data for species " - +name); - } - } - } - } - - // done adding species. - th->freezeSpecies(); - - th->saveSpeciesData(db); - - // Perform any required subclass-specific initialization. - string id = ""; - th->initThermoXML(phase, id); - - return true; - } - - - - /** - * This function returns true if two reactions are duplicates of - * one another, and false otherwise. The input arguments are two - * maps from species number to stoichiometric coefficient, one for - * each reaction. The reactions are considered duplicates if their - * stoichiometric coefficients have the same ratio for all - * species. - */ - doublereal isDuplicateReaction(map& r1, - map& r2) { - - map::const_iterator b = r1.begin(), e = r1.end(); - int k1 = b->first; - doublereal ratio = 0.0; - if (r1[k1] == 0.0 || r2[k1] == 0.0) goto next; - ratio = r2[k1]/r1[k1]; - ++b; - for (; b != e; ++b) { - k1 = b->first; - if (r1[k1] == 0.0 || r2[k1] == 0.0) goto next; - if (fabs(r2[k1]/r1[k1] - ratio) > 1.e-8) - goto next; - } - return ratio; -next: - ratio = 0.0; - b = r1.begin(); - k1 = b->first; - if (r1[k1] == 0.0 || r2[-k1] == 0.0) return 0.0; - ratio = r2[-k1]/r1[k1]; - ++b; - for (; b != e; ++b) { - k1 = b->first; - if (r1[k1] == 0.0 || r2[-k1] == 0.0) return 0.0; - if (fabs(r2[-k1]/r1[k1] - ratio) > 1.e-8) - return 0.0; - } - return ratio; + // Set equation of state parameters. The parameters are + // specific to each subclass of ThermoPhase, so this is done + // by method setParametersFromXML in each subclass. + if (phase.hasChild("thermo")) { + const XML_Node& eos = phase.child("thermo"); + th->setParametersFromXML(eos); } - /** - * Install an individual reaction into a kinetics manager. The - * data for the reaction is in the xml_node r. In other words, r - * points directly to a ctml element named "reaction". i refers - * to the number id of the reaction in the kinetics object. - * - * @param i Reaction number. - * @param r XML_Node containing reaction data. - * @param k Kinetics manager to which reaction will be added. - * @param default_phase ... - * @param rule Rule for handling reactions with missing species - * (skip or flag as error) - * @param validate_rxn If true, check that this reaction is not a - * duplicate of one already entered, and check that the reaction - * balances. - */ - static bool installReaction(int i, const XML_Node& r, Kinetics* k, - string default_phase, int rule, - bool validate_rxn) { - - Kinetics& kin = *k; - /* - * 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 the kinetics object, k, at the end - * of the routine. (Someday this may be rewritten to skip building - * the ReactionData object). - */ - ReactionData rdata; - rdata.reactionType = ELEMENTARY_RXN; // default - vector_int reac, prod; - string eqn, type; - int nn, eqlen; - vector_fp dummy; - - // check to see if the reaction is specified to be a duplicate - // of another reaction, or to allow a negative pre-exponential. - int dup = 0; - if (r.hasAttrib("duplicate")) dup = 1; - int negA = 0; - if (r.hasAttrib("negative_A")) negA = 1; + /*************************************************************** + * Add the elements. + ***************************************************************/ + th->addElementsFromXML(phase); - /* - * This seemingly simple expression goes and finds the child element, - * "equation". Then it treats all of the contents of the "equation" - * as a string, and returns it the variable eqn. We post process - * the string to convert [ and ] characters into < and >, which - * cannot be stored in an XML file. - */ - if (r.hasChild("equation")) - eqn = r("equation"); - else - eqn = ""; + /*************************************************************** + * Add the species. + * + * Species definitions may be imported from multiple + * sources. For each one, a speciesArray element must be + * present. + ***************************************************************/ + XML_Node* db = 0; + vector sparrays; + phase.getChildren("speciesArray", sparrays); + int jsp, nspa = static_cast(sparrays.size()); + vector dbases; + vector_int sprule(nspa,0); - eqlen = static_cast(eqn.size()); - for (nn = 0; nn < eqlen; nn++) { - if (eqn[nn] == '[') eqn[nn] = '<'; - if (eqn[nn] == ']') eqn[nn] = '>'; - } + // loop over the speciesArray elements + for (jsp = 0; jsp < nspa; jsp++) { + const XML_Node& species = *sparrays[jsp]; - bool ok; - // get the reactants - ok = getReagents(r, kin, 1, default_phase, rdata.reactants, - rdata.rstoich, rdata.order, rule); - - /* - * Get the products. We store the id of products in rdata.products - */ - ok = ok && getReagents(r, kin, -1, default_phase, rdata.products, - rdata.pstoich, dummy, rule); - - // if there was a problem getting either the reactants or the products, - // then abort. - if (!ok) return false; - - // check whether the reaction is specified to be - // reversible. Default is irreversible. - rdata.reversible = false; - string isrev = r["reversible"]; - if (isrev == "yes" || isrev == "true") - rdata.reversible = true; - - - string typ = r["type"]; - - /* - * 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) - throw CanteraError("installReaction", - "reaction orders may only be given for " - "irreversible reactions"); - rdata.global = true; - } - - - /* - * Seaarch the reaction element for the attribute "type". - * If found, then branch on the type, to fill in appropriate - * fields in rdata. - */ - - 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 == "surface") { - rdata.reactionType = SURFACE_RXN; - } - else if (typ == "edge") { - rdata.reactionType = EDGE_RXN; - } - else if (typ != "") - throw CanteraError("installReaction", - "Unknown reaction type: " + typ); - - /* - * Look for undeclared duplicate reactions. - */ - if (validate_rxn) { - doublereal c = 0.0; - - map rxnstoich; - rxnstoich.clear(); - int nr = rdata.reactants.size(); - for (nn = 0; nn < nr; nn++) { - rxnstoich[-1 - rdata.reactants[nn]] -= rdata.rstoich[nn]; - } - int np = rdata.products.size(); - for (nn = 0; nn < np; nn++) { - rxnstoich[rdata.products[nn]+1] += rdata.pstoich[nn]; - } - int nrxns = static_cast(_reactiondata.size()); - for (nn = 0; nn < nrxns; nn++) { - if ((int(rdata.reactants.size()) == _nr[nn]) - && (rdata.reactionType == _typ[nn])) { - c = isDuplicateReaction(rxnstoich, _reactiondata[nn]); - if (c > 0.0 - || (c < 0.0 && rdata.reversible) - || (c < 0.0 && _rev[nn])) { - if ((!dup || !_dup[nn])) { - string msg = string("Undeclared duplicate reactions detected: \n") - +"Reaction "+int2str(nn+1)+": "+_eqn[nn] - +"\nReaction "+int2str(i+1)+": "+eqn+"\n"; - _reactiondata.clear(); - _eqn.clear(); - _rev.clear(); - _nr.clear(); - _typ.clear(); - _dup.clear(); - throw CanteraError("installReaction",msg); - } - } - } - } - _dup.push_back(dup); - _rev.push_back(rdata.reversible); - _eqn.push_back(eqn); - _nr.push_back(rdata.reactants.size()); - _typ.push_back(rdata.reactionType); - _reactiondata.push_back(rxnstoich); - } - - rdata.equation = eqn; - rdata.number = i; - rdata.rxn_number = i; - - /* - * Read the rate coefficient data from the XML file. Trigger an - * exception for negative A unless specifically authorized. - */ - getRateCoefficient(r.child("rateCoeff"), kin, rdata, negA); - - /* - * Check to see that the elements balance in the reaction. - * Throw an error if they don't - */ - if (validate_rxn) { - checkRxnElementBalance(kin, rdata); + // If the speciesArray element has a child element + // + // then set sprule[jsp] to 1, so + // that any species with an undeclared element will be + // quietly skipped when importing species. + if (species.hasChild("skip")) { + const XML_Node& sk = species.child("skip"); + string eskip = sk["element"]; + if (eskip == "undeclared") { + sprule[jsp] = 1; } + string dskip = sk["species"]; + if (dskip == "duplicate") { + sprule[jsp] += 10; + } + } - /* - * Ok we have read everything in about the reaction. Add it - * to the kinetics object by calling the Kinetics member function, - * addReaction() - */ - kin.addReaction(rdata); - return true; + string fname, idstr; + + // get a pointer to the node containing the species + // definitions for the species declared in this + // speciesArray element. This may be in the local file + // containing the phase element, or may be in another + // file. + db = get_XML_Node(species["datasrc"], &phase.root()); + + // add this node to the list of species database nodes. + dbases.push_back(db); } + // if the phase has a species thermo manager already installed, + // delete it since we are adding new species. + delete &th->speciesThermo(); - /** - * Take information from the XML tree, p, about reactions - * and install them into the kinetics object, kin. - * default_phase is the default phase to assume when - * looking up species. - * - * At this point, p usually refers to the phase xml element. - * One of the children of this element is reactionArray, - * the element which determines where in the xml file to - * look up the reaction rate data pertaining to the phase. - * - * On return, if reaction instantiation goes correctly, return true. - * If there is a problem, return false. - */ - bool installReactionArrays(const XML_Node& p, Kinetics& kin, - string default_phase, bool check_for_duplicates) { + // create a new species thermo manager. Function + // 'newSpeciesThermoMgr' looks at the species in the database + // to see what thermodynamic property parameterizations are + // used, and selects a class that can handle the + // parameterizations found. + SpeciesThermo* spth = newSpeciesThermoMgr(dbases); - if (_rxns == 0) { - _rxns = new rxninfo; - } - _eqn.clear(); - _dup.clear(); - _nr.clear(); - _typ.clear(); - _reactiondata.clear(); - _rev.clear(); + // install it in the phase object + th->setSpeciesThermo(spth); + SpeciesThermo& spthermo = th->speciesThermo(); - vector rarrays; - int itot = 0; - /* - * Search the children of the phase element for the - * xml element named reactionArray. If we can't find it, - * then return signaling having not found any reactions. - * Apparently, we allow multiple reactionArray elements here - * Each one will be processed sequentially, with the - * end result being purely additive. - */ - p.getChildren("reactionArray",rarrays); - int na = static_cast(rarrays.size()); - if (na == 0) return false; - for (int n = 0; n < na; n++) { - /* - * Go get a reference to the current xml element, - * reactionArray. We will process this element now. - */ - const XML_Node& rxns = *rarrays[n]; - /* - * The reactionArray element has an attribute called, - * datasrc. The value of the attribute is the xml - * element comprising the top of the - * tree of reactions for the phase. - * Find this datasrc element starting with the root - * of the current xml node. - */ - const XML_Node* rdata = get_XML_Node(rxns["datasrc"], &rxns.root()); - //const XML_Node* rdata = find_XML(rxns["datasrc"],&rxns.root(), - // "","","reactionData"); - /* - * If the reactionArray element has a child element named - * "skip", and if the attribute of skip called "species" has - * a value of "undeclared", we will set rxnrule = 1. - * rxnrule is passed to the routine that parses each individual - * reaction. I believe what this means is that the parser will - * skip all reactions containing an undefined species without - * throwing an error condition. - */ - int rxnrule = 0; - if (rxns.hasChild("skip")) { - const XML_Node& sk = rxns.child("skip"); - string sskip = sk["species"]; - if (sskip == "undeclared") { - rxnrule = 1; - } - } - int i, nrxns = 0; - /* - * Search for child elements called include. We only include - * a reaction if it's tagged by one of the include fields. - * Or, we include all reactions if there are no include fields. - */ - vector incl; - rxns.getChildren("include",incl); - int ninc = static_cast(incl.size()); + // used to check that each species is declared only once + map declared; - vector allrxns; - rdata->getChildren("reaction",allrxns); - nrxns = static_cast(allrxns.size()); - // if no 'include' directive, then include all reactions - if (ninc == 0) { - for (i = 0; i < nrxns; i++) { - const XML_Node* r = allrxns[i]; - if (r) { - if (installReaction(itot, *r, &kin, - default_phase, rxnrule, check_for_duplicates)) ++itot; - } - } + int i, k = 0; + + // loop over the species arrays + for (jsp = 0; jsp < nspa; jsp++) { + + const XML_Node& species = *sparrays[jsp]; + db = dbases[jsp]; + + // Get the array of species name strings. + vector spnames; + getStringArray(species, spnames); + int nsp = static_cast(spnames.size()); + + // if 'all' is specified, then add all species + // defined in this database to the phase + if (nsp == 1 && spnames[0] == "all") { + vector allsp; + db->getChildren("species",allsp); + nsp = static_cast(allsp.size()); + spnames.resize(nsp); + for (int nn = 0; nn < nsp; nn++) { + spnames[nn] = (*allsp[nn])["name"]; + } + } + else if (nsp == 1 && spnames[0] == "unique") { + vector uniquesp; + db->getChildren("species",uniquesp); + nsp = static_cast(uniquesp.size()); + spnames.clear(); + spnames.resize(nsp); + string spnm; + for (int nn = 0; nn < nsp; nn++) { + spnm = (*uniquesp[nn])["name"]; + if (!declared[spnm]) spnames[nn] = spnm; + } + } + + string name; + bool skip; + for (i = 0; i < nsp; i++) { + name = spnames[i]; + skip = false; + if (name == "") skip = true; + // Check that every species is only declared once + if (declared[name]) { + if (sprule[jsp] >= 10) + skip = true; + else + throw CanteraError("importPhase", + "duplicate species: "+name); + } + if (!skip) { + declared[name] = true; + + // Find the species in the database by name. + XML_Node* s = db->findByAttr("name",spnames[i]); + if (s) { + if (installSpecies(k, *s, *th, spthermo, sprule[jsp], + spfactory)) + ++k; } else { - for (int nii = 0; nii < ninc; nii++) { - const XML_Node& ii = *incl[nii]; - string imin = ii["min"]; - string imax = ii["max"]; + throw CanteraError("importPhase","no data for species " + +name); + } + } + } + } - string::size_type iwild = string::npos; - if (imax == imin) { - iwild = imin.find("*"); - if (iwild != string::npos) { - imin = imin.substr(0,iwild); - imax = imin; - } - } + // done adding species. + th->freezeSpecies(); - for (i = 0; i < nrxns; i++) { - const XML_Node* r = allrxns[i]; - string rxid; - if (r) { - rxid = (*r)["id"]; - if (iwild != string::npos) { - rxid = rxid.substr(0,iwild); - } - /* - * To decide whether the reaction is included or not - * we do a lexical min max and operation. This - * sometimes has surprising results. - */ - if ((rxid >= imin) && (rxid <= imax)) { - if (installReaction(itot, *r, &kin, - default_phase, rxnrule, check_for_duplicates)) ++itot; - } - } + th->saveSpeciesData(db); + + // Perform any required subclass-specific initialization. + string id = ""; + th->initThermoXML(phase, id); + + return true; + } + + + /* + * This function returns true if two reactions are duplicates of + * one another, and false otherwise. The input arguments are two + * maps from species number to stoichiometric coefficient, one for + * each reaction. The reactions are considered duplicates if their + * stoichiometric coefficients have the same ratio for all + * species. + */ + doublereal isDuplicateReaction(std::map& r1, + std::map& r2) { + + map::const_iterator b = r1.begin(), e = r1.end(); + int k1 = b->first; + doublereal ratio = 0.0; + if (r1[k1] == 0.0 || r2[k1] == 0.0) goto next; + ratio = r2[k1]/r1[k1]; + ++b; + for (; b != e; ++b) { + k1 = b->first; + if (r1[k1] == 0.0 || r2[k1] == 0.0) goto next; + if (fabs(r2[k1]/r1[k1] - ratio) > 1.e-8) + goto next; + } + return ratio; + next: + ratio = 0.0; + b = r1.begin(); + k1 = b->first; + if (r1[k1] == 0.0 || r2[-k1] == 0.0) return 0.0; + ratio = r2[-k1]/r1[k1]; + ++b; + for (; b != e; ++b) { + k1 = b->first; + if (r1[k1] == 0.0 || r2[-k1] == 0.0) return 0.0; + if (fabs(r2[-k1]/r1[k1] - ratio) > 1.e-8) + return 0.0; + } + return ratio; + } + + + /** + * Install an individual reaction into a kinetics manager. The + * data for the reaction is in the xml_node r. In other words, r + * points directly to a ctml element named "reaction". i refers + * to the number id of the reaction in the kinetics object. + * + * @param i Reaction number. + * @param r XML_Node containing reaction data. + * @param k Kinetics manager to which reaction will be added. + * @param default_phase Default phase for locating a species + * @param rule Rule for handling reactions with missing species + * (skip or flag as error) + * @param validate_rxn If true, check that this reaction is not a + * duplicate of one already entered, and check that the reaction + * balances. + * + * @ingroup kineticsmgr + */ + static bool installReaction(int i, const XML_Node& r, Kinetics* k, + string default_phase, int rule, + bool validate_rxn) { + + Kinetics& kin = *k; + /* + * 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 the kinetics object, k, at the end + * of the routine. (Someday this may be rewritten to skip building + * the ReactionData object). + */ + ReactionData rdata; + rdata.reactionType = ELEMENTARY_RXN; // default + vector_int reac, prod; + string eqn, type; + int nn, eqlen; + vector_fp dummy; + + // check to see if the reaction is specified to be a duplicate + // of another reaction, or to allow a negative pre-exponential. + int dup = 0; + if (r.hasAttrib("duplicate")) dup = 1; + int negA = 0; + if (r.hasAttrib("negative_A")) negA = 1; + + + /* + * This seemingly simple expression goes and finds the child element, + * "equation". Then it treats all of the contents of the "equation" + * as a string, and returns it the variable eqn. We post process + * the string to convert [ and ] characters into < and >, which + * cannot be stored in an XML file. + */ + if (r.hasChild("equation")) + eqn = r("equation"); + else + eqn = ""; + + eqlen = static_cast(eqn.size()); + for (nn = 0; nn < eqlen; nn++) { + if (eqn[nn] == '[') eqn[nn] = '<'; + if (eqn[nn] == ']') eqn[nn] = '>'; + } + + + bool ok; + // get the reactants + ok = getReagents(r, kin, 1, default_phase, rdata.reactants, + rdata.rstoich, rdata.order, rule); + + /* + * Get the products. We store the id of products in rdata.products + */ + ok = ok && getReagents(r, kin, -1, default_phase, rdata.products, + rdata.pstoich, dummy, rule); + + // if there was a problem getting either the reactants or the products, + // then abort. + if (!ok) return false; + + // check whether the reaction is specified to be + // reversible. Default is irreversible. + rdata.reversible = false; + string isrev = r["reversible"]; + if (isrev == "yes" || isrev == "true") + rdata.reversible = true; + + + string typ = r["type"]; + + /* + * 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) + throw CanteraError("installReaction", + "reaction orders may only be given for " + "irreversible reactions"); + rdata.global = true; + } + + + /* + * Seaarch the reaction element for the attribute "type". + * If found, then branch on the type, to fill in appropriate + * fields in rdata. + */ + + 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 == "surface") { + rdata.reactionType = SURFACE_RXN; + } + else if (typ == "edge") { + rdata.reactionType = EDGE_RXN; + } + else if (typ != "") + throw CanteraError("installReaction", + "Unknown reaction type: " + typ); + + /* + * Look for undeclared duplicate reactions. + */ + if (validate_rxn) { + doublereal c = 0.0; + + map rxnstoich; + rxnstoich.clear(); + int nr = rdata.reactants.size(); + for (nn = 0; nn < nr; nn++) { + rxnstoich[-1 - rdata.reactants[nn]] -= rdata.rstoich[nn]; + } + int np = rdata.products.size(); + for (nn = 0; nn < np; nn++) { + rxnstoich[rdata.products[nn]+1] += rdata.pstoich[nn]; + } + int nrxns = static_cast(_reactiondata.size()); + for (nn = 0; nn < nrxns; nn++) { + if ((int(rdata.reactants.size()) == _nr[nn]) + && (rdata.reactionType == _typ[nn])) { + c = isDuplicateReaction(rxnstoich, _reactiondata[nn]); + if (c > 0.0 + || (c < 0.0 && rdata.reversible) + || (c < 0.0 && _rev[nn])) { + if ((!dup || !_dup[nn])) { + string msg = string("Undeclared duplicate reactions detected: \n") + +"Reaction "+int2str(nn+1)+": "+_eqn[nn] + +"\nReaction "+int2str(i+1)+": "+eqn+"\n"; + _reactiondata.clear(); + _eqn.clear(); + _rev.clear(); + _nr.clear(); + _typ.clear(); + _dup.clear(); + throw CanteraError("installReaction",msg); + } + } + } + } + _dup.push_back(dup); + _rev.push_back(rdata.reversible); + _eqn.push_back(eqn); + _nr.push_back(rdata.reactants.size()); + _typ.push_back(rdata.reactionType); + _reactiondata.push_back(rxnstoich); + } + + rdata.equation = eqn; + rdata.number = i; + rdata.rxn_number = i; + + /* + * Read the rate coefficient data from the XML file. Trigger an + * exception for negative A unless specifically authorized. + */ + getRateCoefficient(r.child("rateCoeff"), kin, rdata, negA); + + /* + * Check to see that the elements balance in the reaction. + * Throw an error if they don't + */ + if (validate_rxn) { + checkRxnElementBalance(kin, rdata); + } + + /* + * Ok we have read everything in about the reaction. Add it + * to the kinetics object by calling the Kinetics member function, + * addReaction() + */ + kin.addReaction(rdata); + return true; + } + + /* + * Take information from the XML tree, p, about reactions + * and install them into the kinetics object, kin. + * default_phase is the default phase to assume when + * looking up species. + * + * At this point, p usually refers to the phase xml element. + * One of the children of this element is reactionArray, + * the element which determines where in the xml file to + * look up the reaction rate data pertaining to the phase. + * + * On return, if reaction instantiation goes correctly, return true. + * If there is a problem, return false. + */ + bool installReactionArrays(const XML_Node& p, Kinetics& kin, + std::string default_phase, bool check_for_duplicates) { + + if (_rxns == 0) { + _rxns = new rxninfo; + } + _eqn.clear(); + _dup.clear(); + _nr.clear(); + _typ.clear(); + _reactiondata.clear(); + _rev.clear(); + + vector rarrays; + int itot = 0; + /* + * Search the children of the phase element for the + * xml element named reactionArray. If we can't find it, + * then return signaling having not found any reactions. + * Apparently, we allow multiple reactionArray elements here + * Each one will be processed sequentially, with the + * end result being purely additive. + */ + p.getChildren("reactionArray",rarrays); + int na = static_cast(rarrays.size()); + if (na == 0) return false; + for (int n = 0; n < na; n++) { + /* + * Go get a reference to the current xml element, + * reactionArray. We will process this element now. + */ + const XML_Node& rxns = *rarrays[n]; + /* + * The reactionArray element has an attribute called, + * datasrc. The value of the attribute is the xml + * element comprising the top of the + * tree of reactions for the phase. + * Find this datasrc element starting with the root + * of the current xml node. + */ + const XML_Node* rdata = get_XML_Node(rxns["datasrc"], &rxns.root()); + //const XML_Node* rdata = find_XML(rxns["datasrc"],&rxns.root(), + // "","","reactionData"); + /* + * If the reactionArray element has a child element named + * "skip", and if the attribute of skip called "species" has + * a value of "undeclared", we will set rxnrule = 1. + * rxnrule is passed to the routine that parses each individual + * reaction. I believe what this means is that the parser will + * skip all reactions containing an undefined species without + * throwing an error condition. + */ + int rxnrule = 0; + if (rxns.hasChild("skip")) { + const XML_Node& sk = rxns.child("skip"); + string sskip = sk["species"]; + if (sskip == "undeclared") { + rxnrule = 1; + } + } + int i, nrxns = 0; + /* + * Search for child elements called include. We only include + * a reaction if it's tagged by one of the include fields. + * Or, we include all reactions if there are no include fields. + */ + vector incl; + rxns.getChildren("include",incl); + int ninc = static_cast(incl.size()); + + vector allrxns; + rdata->getChildren("reaction",allrxns); + nrxns = static_cast(allrxns.size()); + // if no 'include' directive, then include all reactions + if (ninc == 0) { + for (i = 0; i < nrxns; i++) { + const XML_Node* r = allrxns[i]; + if (r) { + if (installReaction(itot, *r, &kin, + default_phase, rxnrule, check_for_duplicates)) ++itot; + } + } + } + else { + for (int nii = 0; nii < ninc; nii++) { + const XML_Node& ii = *incl[nii]; + string imin = ii["min"]; + string imax = ii["max"]; + + string::size_type iwild = string::npos; + if (imax == imin) { + iwild = imin.find("*"); + if (iwild != string::npos) { + imin = imin.substr(0,iwild); + imax = imin; + } + } + + for (i = 0; i < nrxns; i++) { + const XML_Node* r = allrxns[i]; + string rxid; + if (r) { + rxid = (*r)["id"]; + if (iwild != string::npos) { + rxid = rxid.substr(0,iwild); + } + /* + * To decide whether the reaction is included or not + * we do a lexical min max and operation. This + * sometimes has surprising results. + */ + if ((rxid >= imin) && (rxid <= imax)) { + if (installReaction(itot, *r, &kin, + default_phase, rxnrule, check_for_duplicates)) ++itot; } } } - } - - /* - * Finalize the installation of the kinetics, now that we know - * the true number of reactions in the mechanism, itot. - */ - kin.finalize(); - //writer = 0; - _eqn.clear(); - _dup.clear(); - _nr.clear(); - _typ.clear(); - _reactiondata.clear(); - delete _rxns; - _rxns = 0; - return true; - } - - - /** - * Import a reaction mechanism for a phase or an interface. - * - * @param phase This is an xml node containing a description - * of a phase. Within the phase is a XML element - * called reactionArray containing the location - * of the description of the reactions that make - * up the kinetics object. - * Also within the phase is an XML element called - * phaseArray containing a listing of other phases - * that participate in the kinetics mechanism. - * - * @param th This is a list of ThermoPhase pointers containing - * the phases that participate in the kinetics - * reactions. All of the phases must have already - * been initialized and formed within Cantera. - * However, their pointers should not have been - * added to the Kinetics object; this addition - * is carried out here. - * - * @param k This is a pointer to the kinetics manager class - * that will be initialized with a kinetics - * mechanism. - */ - bool importKinetics(const XML_Node& phase, vector th, - Kinetics* k) { - - if (k == 0) return false; - - Kinetics& kin = *k; - - // This phase will be the default one - string default_phase = phase["id"]; - - bool check_for_duplicates = false; - if (phase.parent()->hasChild("validate")) { - const XML_Node& d = phase.parent()->child("validate"); - if (d["reactions"] == "yes") check_for_duplicates = true; - } - - // if other phases are involved in the reaction mechanism, - // they must be listed in a 'phaseArray' child - // element. Homogeneous mechanisms do not need to include a - // phaseArray element. - - vector phase_ids; - if (phase.hasChild("phaseArray")) { - const XML_Node& pa = phase.child("phaseArray"); - getStringArray(pa, phase_ids); - } - phase_ids.push_back(default_phase); - - int np = static_cast(phase_ids.size()); - int nt = static_cast(th.size()); - - // for each referenced phase, attempt to find its id among those - // phases specified. - bool phase_ok; - - string phase_id; - string msg = ""; - for (int n = 0; n < np; n++) { - phase_id = phase_ids[n]; - phase_ok = false; - - // loop over the supplied 'ThermoPhase' objects representing - // phases, to find an object with the same id. - for (int m = 0; m < nt; m++) { - if (th[m]->id() == phase_id) { - phase_ok = true; - - // if no phase with this id has been added to - //the kinetics manager yet, then add this one - if (kin.phaseIndex(phase_id) < 0) { - kin.addPhase(*th[m]); - } - } - msg += " "+th[m]->id(); - } - if (!phase_ok) { - throw CanteraError("importKinetics", - "phase "+phase_id+" not found. Supplied phases are:"+msg); - } - } - - // allocates arrays, etc. Must be called after the phases have - // been added to 'kin', so that the number of species in each - // phase is known. - kin.init(); - - // Install the reactions. - return installReactionArrays(phase, kin, default_phase, check_for_duplicates); - } - - /** - * Build a single-phase ThermoPhase object with associated kinetics - * mechanism. - */ - bool buildSolutionFromXML(XML_Node& root, string id, string nm, - ThermoPhase* th, Kinetics* k) { - XML_Node* x; - try { - - x = get_XML_NameID(nm, string("#")+id, &root); - // x = get_XML_Node(string("#")+id, &root); - - /* - * Fill in the ThermoPhase object by querying the - * const XML_Node tree located at x. - */ - importPhase(*x, th); - /* - * Create a vector of ThermoPhase pointers of length 1 - * having the current th ThermoPhase as the entry. - */ - vector phases(1); - phases[0] = th; - /* - * Fill in the kinetics object k, by querying the - * const XML_Node tree located by x. The source terms and - * eventually the source term vector will be constructed - * from the list of ThermoPhases in the vector, phases. - */ - importKinetics(*x, phases, k); - - return true; - } - catch (CanteraError) { - throw CanteraError("buildSolutionFromXML","error encountered"); - return false; - } - } - - /** - * Search an XML tree for species data. - * - * This utility routine will search the XML tree for the species - * named by the string, kname. It will return the XML_Node - * pointer. - * Failures of any kind return the null pointer. - */ - const XML_Node *speciesXML_Node(string kname, - const XML_Node *phaseSpecies) { - /* - * First look at the species database. - * -> Look for the subelement "stoichIsMods" - * in each of the species SS databases. - */ - if (!phaseSpecies) return ((const XML_Node *) 0); - string jname; - vector xspecies; - phaseSpecies->getChildren("species", xspecies); - int jj = xspecies.size(); - for (int j = 0; j < jj; j++) { - const XML_Node& sp = *xspecies[j]; - jname = sp["name"]; - if (jname == kname) { - return &sp; - } + } } - return ((const XML_Node *) 0); } + /* + * Finalize the installation of the kinetics, now that we know + * the true number of reactions in the mechanism, itot. + */ + kin.finalize(); + //writer = 0; + _eqn.clear(); + _dup.clear(); + _nr.clear(); + _typ.clear(); + _reactiondata.clear(); + delete _rxns; + _rxns = 0; + return true; + } + /* + * Import a reaction mechanism for a phase or an interface. + * + * @param phase This is an xml node containing a description + * of a phase. Within the phase is a XML element + * called reactionArray containing the location + * of the description of the reactions that make + * up the kinetics object. + * Also within the phase is an XML element called + * phaseArray containing a listing of other phases + * that participate in the kinetics mechanism. + * + * @param th This is a list of ThermoPhase pointers containing + * the phases that participate in the kinetics + * reactions. All of the phases must have already + * been initialized and formed within Cantera. + * However, their pointers should not have been + * added to the Kinetics object; this addition + * is carried out here. + * + * @param k This is a pointer to the kinetics manager class + * that will be initialized with a kinetics + * mechanism. + */ + bool importKinetics(const XML_Node& phase, std::vector th, + Kinetics* k) { + + if (k == 0) return false; + + Kinetics& kin = *k; + + // This phase will be the default one + string default_phase = phase["id"]; + + bool check_for_duplicates = false; + if (phase.parent()->hasChild("validate")) { + const XML_Node& d = phase.parent()->child("validate"); + if (d["reactions"] == "yes") check_for_duplicates = true; + } + + // if other phases are involved in the reaction mechanism, + // they must be listed in a 'phaseArray' child + // element. Homogeneous mechanisms do not need to include a + // phaseArray element. + + vector phase_ids; + if (phase.hasChild("phaseArray")) { + const XML_Node& pa = phase.child("phaseArray"); + getStringArray(pa, phase_ids); + } + phase_ids.push_back(default_phase); + + int np = static_cast(phase_ids.size()); + int nt = static_cast(th.size()); + + // for each referenced phase, attempt to find its id among those + // phases specified. + bool phase_ok; + + string phase_id; + string msg = ""; + for (int n = 0; n < np; n++) { + phase_id = phase_ids[n]; + phase_ok = false; + + // loop over the supplied 'ThermoPhase' objects representing + // phases, to find an object with the same id. + for (int m = 0; m < nt; m++) { + if (th[m]->id() == phase_id) { + phase_ok = true; + + // if no phase with this id has been added to + //the kinetics manager yet, then add this one + if (kin.phaseIndex(phase_id) < 0) { + kin.addPhase(*th[m]); + } + } + msg += " "+th[m]->id(); + } + if (!phase_ok) { + throw CanteraError("importKinetics", + "phase "+phase_id+" not found. Supplied phases are:"+msg); + } + } + + // allocates arrays, etc. Must be called after the phases have + // been added to 'kin', so that the number of species in each + // phase is known. + kin.init(); + + // Install the reactions. + return installReactionArrays(phase, kin, default_phase, check_for_duplicates); + } + + /* + * Build a single-phase ThermoPhase object with associated kinetics + * mechanism. + */ + bool buildSolutionFromXML(XML_Node& root, std::string id, std::string nm, + ThermoPhase* th, Kinetics* k) { + XML_Node* x; + try { + + x = get_XML_NameID(nm, string("#")+id, &root); + // x = get_XML_Node(string("#")+id, &root); + + /* + * Fill in the ThermoPhase object by querying the + * const XML_Node tree located at x. + */ + importPhase(*x, th); + /* + * Create a vector of ThermoPhase pointers of length 1 + * having the current th ThermoPhase as the entry. + */ + vector phases(1); + phases[0] = th; + /* + * Fill in the kinetics object k, by querying the + * const XML_Node tree located by x. The source terms and + * eventually the source term vector will be constructed + * from the list of ThermoPhases in the vector, phases. + */ + importKinetics(*x, phases, k); + + return true; + } + catch (CanteraError) { + throw CanteraError("buildSolutionFromXML","error encountered"); + return false; + } + } + + /* + * Search an XML tree for species data. + * + * This utility routine will search the XML tree for the species + * named by the string, kname. It will return the XML_Node + * pointer. + * Failures of any kind return the null pointer. + */ + const XML_Node *speciesXML_Node(std::string kname, + const XML_Node *phaseSpeciesData) { + /* + * First look at the species database. + * -> Look for the subelement "stoichIsMods" + * in each of the species SS databases. + */ + if (!phaseSpeciesData) return ((const XML_Node *) 0); + string jname; + vector xspecies; + phaseSpeciesData->getChildren("species", xspecies); + int jj = xspecies.size(); + for (int j = 0; j < jj; j++) { + const XML_Node& sp = *xspecies[j]; + jname = sp["name"]; + if (jname == kname) { + return &sp; + } + } + return ((const XML_Node *) 0); + } + } diff --git a/Cantera/src/importCTML.h b/Cantera/src/importCTML.h index a3a1e2d1e..9194bccb6 100755 --- a/Cantera/src/importCTML.h +++ b/Cantera/src/importCTML.h @@ -14,57 +14,439 @@ #define CT_IMPORTCTML_H #include -//using namespace std; #include "ThermoPhase.h" #include "Kinetics.h" namespace Cantera { - class Kinetics; - class SpeciesThermoFactory; - //class ThermoPhase; - class XML_Node; + class Kinetics; + class SpeciesThermoFactory; + class XML_Node; - bool isCTMLFile(std::string infile); + //! Search for an XML_Node either wiithin an existing XML tree structure, or in another file, + //! based on the file name or the XML id attribute. + /*! + * This routine will locate an XML node in either the input + * XML tree or in another input file specified by the file + * part of the file_ID string. Searches are based on the + * ID attribute of the XML element only. + * + * @param file_ID This is a concatenation of two strings seperated + * by the "#" character. The string before the + * pound character is the file name of an xml + * file to carry out the search. The string after + * the # character is the ID attribute + * of the xml element to search for. + * The string is interpreted as a file string if + * no # character is in the string. + * + * @param root If the file string is empty, searches for the + * xml element with matching ID attribute are + * carried out from this XML node. + * + * @return + * This routine will process the XML file, creating an XML + * tree structure. It returns a pointer to the top of the tree. + * + * + * For example, + * @code + * + * XML_Node* xn = get_XML_Node("phase", "gri30.xml#gri30_mix", 0); + * + * @endcode + * + * will search in the file gri30.xml for an XML element of the following form, where + * the XML element name, phase, is an optional hit: + * @code + * & r1, + std::map& r2); + + + //! This function will check a specific reaction to see if the elements balance. + /*! + * @param kin Kinetics object + * @param rdata Object containing the information about one reaction + * @param errorTolerance double containing the error tolerance. + */ + void checkRxnElementBalance(Kinetics& kin, + const ReactionData &rdata, + doublereal errorTolerance = 1.0e-3); - XML_Node* get_XML_NameID(const std::string& nameTarget, - const std::string& file_ID, XML_Node* root); + //! Read the rate coefficient data from the XML file. + /*! + * Extract the rate coefficient for a reaction from the xml node, kf. + * kf should point to a XML element named "rateCoeff". + * rdata is the partially filled ReactionData object for the reaction. + * This function will fill in more fields in the ReactionData object. + * + * @param kf XML_Node containing information about the rate coefficients. + * @param kin kinetics manager + * @param rdata ReactionData referece + * @param negA Boolean indicating whether negative A's are ok. + * + * Trigger anexception for negative A unless specifically authorized. + */ + void getRateCoefficient(const XML_Node& kf, kinetics_t& kin, + ReactionData& rdata, int negA); - bool installSpecies(int k, const XML_Node& s, thermo_t& p, - SpeciesThermo& spthermo, int rule, - SpeciesThermoFactory* spfactory); + + //! Create a new ThermoPhase object and initializes it according to the XML tree database. + /*! + * This routine first looks up the + * identity of the model for the solution thermodynamics in the + * model attribute of the thermo child of the xml phase + * node. Then, it does a string lookup using Cantera's internal ThermoPhase Factory routines + * on the model to figure out + * what ThermoPhase derived class should be assigned. It creates a new + * instance of that class, and then calls importPhase() to + * populate that class with the correct parameters from the XML + * tree. + * + * @param phase XML_Node reference pointing to the phase XML element. + * + * @return + * Returns a pointer to the completed and initialized ThermoPhase object. + * + * @ingroup inputfiles + */ + ThermoPhase* newPhase(XML_Node& phase); - bool importPhase(XML_Node& phase, ThermoPhase* th, - SpeciesThermoFactory* spfactory = 0); + //! Create a new ThermoPhase object and initializes it according to a file name and phase id + /*! + * This routine will first find a matching XML phase description given a file name + * and phase id string. Then, this routine will looks up the + * identity of the model for the solution thermodynamics in the + * model attribute of the thermo child of the xml phase + * node. Then, it does a string lookup using Cantera's internal ThermoPhase Factory routines + * on the model to figure out + * what ThermoPhase derived class should be assigned. It creates a new + * instance of that class, and then calls importPhase() to + * populate that class with the correct parameters from the XML + * tree. + * + * @param file String file name to look up the phase + * @param id XML id attribute of the phase. + * + * @return + * Returns a pointer to the completed and initialized ThermoPhase object. + * + * @ingroup inputfiles + */ + ThermoPhase* newPhase(std::string file, std::string id); - doublereal isDuplicateReaction(std::map& r1, - std::map& r2); - - bool importKinetics(const XML_Node& phase, std::vector th, - Kinetics* kin); - bool installReactionArrays(const XML_Node& parent, Kinetics& kin, - std::string default_phase, bool check_for_duplicates = false); - - /*** - * This function will check a specific reaction to see if it the - * elements balance. - */ - void checkRxnElementBalance(Kinetics& kin, - const ReactionData &rdata, doublereal errorTolerance = 1.0e-3); - - void getRateCoefficient(const XML_Node& kf, kinetics_t& kin, - ReactionData& rdata, int negA); - - ThermoPhase* newPhase(XML_Node& phase); - ThermoPhase* newPhase(std::string file, std::string id); - bool buildSolutionFromXML(XML_Node& root, std::string id, std::string nm, - ThermoPhase* th, Kinetics* k); - - const XML_Node *speciesXML_Node(std::string kname, - const XML_Node *phaseSpecies); + //! Install information about reactions into the kinetics object, kin. + /*! + * At this point, parent usually refers to the phase xml element. + * One of the children of this element is reactionArray, + * the element which determines where in the xml file to + * look up the reaction rate data. + * + * This is a wrapper routine around the static function installReaction() + * + * @param p parent XML phase element + * @param kin Kinetics object to install reactions into + * @param default_phase The default_phase is the default phase to assume when + * looking up species. + * @param check_for_duplicates Check for reactions with exactly the same + * reactants and products. + * + * @return + * On return, if reaction instantiation goes correctly, return true. + * If there is a problem, return false. + * + * @ingroup kineticsmgr + */ + bool installReactionArrays(const XML_Node& p, Kinetics& kin, + std::string default_phase, + bool check_for_duplicates = false); + + //! Import a reaction mechanism for a phase or an interface. + /*! + * This routine will import a reaction mechanism into a + * kinetics object. The reaction + * mechanism may either be homogeneous or heterogeneous, + * involving multiple ThermoPhase objects. + * The hosting phase should be included as the first argument. + * For example, if phase I is an interface phase between bulk + * phases A and B. Then, the XML_Node for phase I should be + * the first argument. + * The vector of %ThermoPhase objects should be consist of pointers + * to phases I, A, and B. + * + * @param phase This is an xml node containing a description + * of a phase. Within the phase is a XML element + * called reactionArray containing the location + * of the description of the reactions that make + * up the kinetics object. + * Also within the phase is an XML element called + * phaseArray containing a listing of other phases + * that participate in the kinetics mechanism. + * + * @param th This is a list of ThermoPhase pointers containing + * the phases that participate in the kinetics + * reactions. All of the phases must have already + * been initialized and formed within Cantera. + * However, their pointers should not have been + * added to the Kinetics object; this addition + * is carried out here. + * + * @param kin This is a pointer to a bare kinetics manager class + * that will be initialized with the kinetics + * mechanism. + * + * @ingroup kineticsmgr + * + */ + bool importKinetics(const XML_Node& phase, std::vector th, + Kinetics* kin); + + //!Build a single-phase ThermoPhase object with associated kinetics mechanism. + /*! + * In a single call, this routine initializes a ThermoPhase object and a + * homogenous kinetics object for a phase. + * + * @param root pointer to the XML tree which will be searched to find the + * XML phase element. + * + * @param id Name of the phase to be searched for. + * @param nm Name of the XML element. Should be "phase" + * @param th Pointer to a bare ThermoPhase object, which will be initialized + * by this operaton. + * @param k Pointer to a bare Kinetics object, which will be initialized + * by this operation to a homogeneous kinetics manager + * + * @return + * Returns true if all went well. If there are errors, it will return false. + * + * For Example + * + * @code + * ThermoPhase *th = new ThermoPhase(); + * Kinetics *k = new Kinetics(); + * XML_Node *root = get_XML_File("gri30.xml"); + * ok = buildSolutionFromXML(root, "gri30_mix", "phase", th, k) + * @endcode + * + * @ingroup inputfiles + * @see importKinetics() + */ + bool buildSolutionFromXML(XML_Node& root, std::string id, std::string nm, + ThermoPhase* th, Kinetics* k); + + //! Search an XML tree for species data. + /*! + * + * This utility routine will search the XML tree for the species + * named by the string, kname. It will return the XML_Node + * pointer. + * Failures of any kind return the null pointer. + * + * @param kname species Name + * @param phaseSpeciesData Pointer to the phase XML node pertaining to the + * species database for the phase to be found + * + * @return + * Returns a pointer to teh XML node containing the species data. + * + * @ingroup inputfiles + */ + const XML_Node *speciesXML_Node(std::string kname, + const XML_Node *phaseSpeciesData); }