diff --git a/Cantera/src/ConstDensityThermo.h b/Cantera/src/ConstDensityThermo.h index 7b1ca82c3..b02aabdf7 100755 --- a/Cantera/src/ConstDensityThermo.h +++ b/Cantera/src/ConstDensityThermo.h @@ -19,6 +19,7 @@ #include "mix_defs.h" #include "ThermoPhase.h" #include "SpeciesThermo.h" +#include "utilities.h" namespace Cantera { diff --git a/Cantera/src/Elements.cpp b/Cantera/src/Elements.cpp index 011e85fd4..c1b10c334 100644 --- a/Cantera/src/Elements.cpp +++ b/Cantera/src/Elements.cpp @@ -22,12 +22,14 @@ #include "Elements.h" #include "xml.h" +#include "ctexceptions.h" #ifdef USE_DGG_CODE #include #endif namespace Cantera { + /* awData structure */ /** diff --git a/Cantera/src/Elements.h b/Cantera/src/Elements.h index 60ee12fdc..f5a0daf46 100644 --- a/Cantera/src/Elements.h +++ b/Cantera/src/Elements.h @@ -18,30 +18,12 @@ #undef USE_DGG_CODE #include "ct_defs.h" -#include "ctexceptions.h" +//#include "ctexceptions.h" namespace Cantera { class XML_Node; - -#ifdef INCL_DEPRECATED_METHODS - /** - * Holds element name and atomic weight. Used for output of - * element properties and initial initialization only. - */ - struct ElementData { - string name; - doublereal atomicWeight; - }; -#endif - - class ElementRangeError : public CanteraError { - public: - ElementRangeError(string func, int m, int mmax) : - CanteraError(func, "Element index " + int2str(m) + - " outside valid range of 0 to " + int2str(mmax-1)) {} - }; - + class ElementRangeError; /** Elements Class: Object contains the elements that make up species. * @@ -78,11 +60,6 @@ namespace Cantera { */ int nElements() const { return m_mm; } -#ifdef INCL_DEPRECATED_METHODS - /// Returns an ElementData struct that contains the parameters - /// for element m. - ElementData element(int m) const; -#endif /** Function that returns the index of an element. * * Index of element named \c name. The index is an integer @@ -179,14 +156,7 @@ namespace Cantera { * Note, a string search is the primary way to identify elements. */ vector m_elementNames; -#ifdef USE_DGG_CODE - /** - * Map of elements to indecises - * - * NOTE: this is redundent - */ - map m_definedElements; -#endif + /** * Number of Constituents Objects that use this object * diff --git a/Cantera/src/GasKinetics.cpp b/Cantera/src/GasKinetics.cpp index ba263552d..22bb6d077 100755 --- a/Cantera/src/GasKinetics.cpp +++ b/Cantera/src/GasKinetics.cpp @@ -677,11 +677,3 @@ namespace Cantera { } } - - - - - - - - diff --git a/Cantera/src/IdealGasPhase.h b/Cantera/src/IdealGasPhase.h index e0b23d741..3df19b1b0 100644 --- a/Cantera/src/IdealGasPhase.h +++ b/Cantera/src/IdealGasPhase.h @@ -20,6 +20,7 @@ #include "mix_defs.h" #include "ThermoPhase.h" #include "SpeciesThermo.h" +#include "utilities.h" namespace Cantera { diff --git a/Cantera/src/InterfaceKinetics.cpp b/Cantera/src/InterfaceKinetics.cpp index 4e8925bfb..97bb8e40e 100644 --- a/Cantera/src/InterfaceKinetics.cpp +++ b/Cantera/src/InterfaceKinetics.cpp @@ -46,7 +46,8 @@ namespace Cantera { m_nrev(0), m_surf(0), m_integrator(0), - m_finalized(false) + m_finalized(false), + m_has_coverage_dependence(false) { m_kdata = new InterfaceKineticsData; m_kdata->m_temp = 0.0; @@ -69,6 +70,11 @@ namespace Cantera { void InterfaceKinetics:: _update_rates_T() { _update_rates_phi(); + if (m_has_coverage_dependence) { + m_surf->getCoverages(m_conc.begin()); + m_rates.update_C(m_conc.begin()); + m_redo_rates = true; + } doublereal T = thermo(surfacePhaseIndex()).temperature(); if (T != m_kdata->m_temp || m_redo_rates) { m_kdata->m_logtemp = log(T); @@ -104,16 +110,6 @@ namespace Cantera { _update_rates_C() { int n; - /** - * First evaluate the coverage-dependent terms in the reaction - * rates. - */ - // UNCOMMENT and fix - //m_surf->getCoverages(m_conc.begin()); - //m_rates.update_C(m_conc.begin()); - //m_rates.update(m_kdata->m_temp, - // m_kdata->m_logtemp, m_kdata->m_rfn.begin()); - int np = nPhases(); for (n = 0; n < np; n++) { /* @@ -159,8 +155,9 @@ namespace Cantera { } // compute Delta mu^0 for all reversible reactions - m_reactantStoich.decrementReactions(m_mu0.begin(), m_rkc.begin()); - m_revProductStoich.incrementReactions(m_mu0.begin(), m_rkc.begin()); + //m_reactantStoich.decrementReactions(m_mu0.begin(), m_rkc.begin()); + //m_revProductStoich.incrementReactions(m_mu0.begin(), m_rkc.begin()); + m_rxnstoich.getRevReactionDelta(m_ii, m_mu0.begin(), m_rkc.begin()); for (i = 0; i < m_nrev; i++) { irxn = m_revindex[i]; @@ -183,8 +180,8 @@ namespace Cantera { if (m_nrev > 0) { int n, nsp, k, ik=0; - doublereal rt = GasConstant*thermo(0).temperature(); - doublereal rrt = 1.0/rt; + //doublereal rt = GasConstant*thermo(0).temperature(); + // doublereal rrt = 1.0/rt; int np = nPhases(); for (n = 0; n < np; n++) { thermo(n).getChemPotentials(dmu.begin() + m_start[n]); @@ -197,8 +194,9 @@ namespace Cantera { } // compute Delta mu^ for all reversible reactions - m_reactantStoich.decrementReactions(dmu.begin(), rmu.begin()); - m_revProductStoich.incrementReactions(dmu.begin(), rmu.begin()); + //m_reactantStoich.decrementReactions(dmu.begin(), rmu.begin()); + //m_revProductStoich.incrementReactions(dmu.begin(), rmu.begin()); + m_rxnstoich.getRevReactionDelta(m_ii, dmu.begin(), rmu.begin()); for (i = 0; i < m_nrev; i++) { irxn = m_revindex[i]; @@ -231,9 +229,10 @@ namespace Cantera { fill(kc, kc + m_ii, 0.0); - m_reactantStoich.decrementReactions(m_mu0.begin(), kc); - m_revProductStoich.incrementReactions(m_mu0.begin(), kc); - m_irrevProductStoich.incrementReactions(m_mu0.begin(), kc); + //m_reactantStoich.decrementReactions(m_mu0.begin(), kc); + //m_revProductStoich.incrementReactions(m_mu0.begin(), kc); + //m_irrevProductStoich.incrementReactions(m_mu0.begin(), kc); + m_rxnstoich.getReactionDelta(m_ii, m_mu0.begin(), kc); for (i = 0; i < m_ii; i++) { kc[i] = exp(-kc[i]*rrt); @@ -266,10 +265,11 @@ namespace Cantera { // compute the change in electrical potential energy for each // reaction. This will only be non-zero if a potential // difference is present. - fill(m_rwork.begin(), m_rwork.begin() + m_ii, 0.0); - m_reactantStoich.decrementReactions(m_pot.begin(), m_rwork.begin()); - m_revProductStoich.incrementReactions(m_pot.begin(), m_rwork.begin()); - m_irrevProductStoich.incrementReactions(m_pot.begin(), m_rwork.begin()); + //fill(m_rwork.begin(), m_rwork.begin() + m_ii, 0.0); + //m_reactantStoich.decrementReactions(m_pot.begin(), m_rwork.begin()); + //m_revProductStoich.incrementReactions(m_pot.begin(), m_rwork.begin()); + //m_irrevProductStoich.incrementReactions(m_pot.begin(), m_rwork.begin()); + m_rxnstoich.getReactionDelta(m_ii, m_pot.begin(), m_rwork.begin()); // modify the reaction rates. Only modify those with a // non-zero activation energy, and do not decrease the @@ -355,14 +355,16 @@ namespace Cantera { multiply_each(ropr.begin(), ropr.end(), m_rkc.begin()); // multiply ropf by concentration products - m_reactantStoich.multiply(m_conc.begin(), ropf.begin()); + m_rxnstoich.multiplyReactants(m_conc.begin(), ropf.begin()); + //m_reactantStoich.multiply(m_conc.begin(), ropf.begin()); // for reversible reactions, multiply ropr by concentration // products - m_revProductStoich.multiply(m_conc.begin(), ropr.begin()); + m_rxnstoich.multiplyRevProducts(m_conc.begin(), ropr.begin()); + //m_revProductStoich.multiply(m_conc.begin(), ropr.begin()); // do global reactions - m_globalReactantStoich.power(m_conc.begin(), ropf.begin()); + //m_globalReactantStoich.power(m_conc.begin(), ropf.begin()); for (int j = 0; j != m_ii; ++j) { ropnet[j] = ropf[j] - ropr[j]; @@ -388,21 +390,37 @@ namespace Cantera { void InterfaceKinetics:: addReaction(const ReactionData& r) { - int nr = r.reactants.size(); + // int nr = r.reactants.size(); - // a global reaction is idnetified as one with + // a global reaction is identified as one with // a reactant stoichiometric coefficient not equal // to the molecularity for some reactant - bool isglobal = false; - for (int n = 0; n < nr; n++) { - if (r.rstoich[n] != int(r.order[n])) { - isglobal = true; break; - } - } - if (isglobal) - addGlobalReaction(r); - else - addElementaryReaction(r); +// bool isglobal = false; +// for (int n = 0; n < nr; n++) { +// if (r.rstoich[n] != int(r.order[n])) { +// isglobal = true; break; +// } +// } + // if (isglobal) + // addGlobalReaction(r); + //else + + // if (r.global) + // cout << r.equation << " is global " << endl; + + addElementaryReaction(r); + + //if (r.global) { + // int nr = r.order.size(); + // vector_fp ordr(nr); + // for (int n = 0; n < nr; n++) { + // ordr[n] = r.order[n] - r.rstoich[n]; + // cout << r.reactants[n] << " " << r.order[n] << " " << ordr[n] << endl; + // } + //m_globalReactantStoich.add( reactionNumber(), + // r.reactants, ordr); + //} + // if (r.reactionType == ELEMENTARY_RXN) // addElementaryReaction(r); @@ -425,6 +443,9 @@ namespace Cantera { // install rate coeff calculator vector_fp rp = r.rateCoeffParameters; int ncov = r.cov.size(); + if (ncov > 3) { + m_has_coverage_dependence = true; + } for (int m = 0; m < ncov; m++) rp.push_back(r.cov[m]); iloc = m_rates.install( reactionNumber(), r.rateCoeffType, rp.size(), @@ -437,31 +458,32 @@ namespace Cantera { } - void InterfaceKinetics:: - addGlobalReaction(const ReactionData& r) { +// void InterfaceKinetics:: +// addGlobalReaction(const ReactionData& r) { - int iloc; - // install rate coeff calculator - vector_fp rp = r.rateCoeffParameters; - int ncov = r.cov.size(); - for (int m = 0; m < ncov; m++) rp.push_back(r.cov[m]); - iloc = m_rates.install( reactionNumber(), - r.rateCoeffType, rp.size(), - rp.begin() ); +// int iloc; +// // install rate coeff calculator +// vector_fp rp = r.rateCoeffParameters; +// int ncov = r.cov.size(); +// for (int m = 0; m < ncov; m++) rp.push_back(r.cov[m]); +// iloc = m_rates.install( reactionNumber(), +// r.rateCoeffType, rp.size(), +// rp.begin() ); +// // store activation energy +// m_E.push_back(r.rateCoeffParameters[2]); +// // add constant term to rate coeff value vector +// m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]); - // add constant term to rate coeff value vector - m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]); +// int nr = r.order.size(); +// vector_fp ordr(nr); +// for (int n = 0; n < nr; n++) { +// ordr[n] = r.order[n] - r.rstoich[n]; +// } +// m_globalReactantStoich.add( reactionNumber(), +// r.reactants, ordr); - int nr = r.order.size(); - vector_fp ordr(nr); - for (int n = 0; n < nr; n++) { - ordr[n] = r.order[n] - r.rstoich[n]; - } - m_globalReactantStoich.add( reactionNumber(), - r.reactants, ordr); - - registerReaction( reactionNumber(), GLOBAL_RXN, iloc); - } +// registerReaction( reactionNumber(), GLOBAL_RXN, iloc); +// } void InterfaceKinetics::installReagents(const ReactionData& r) { @@ -473,6 +495,16 @@ namespace Cantera { int rnum = reactionNumber(); + // vectors rk and pk are lists of species numbers, with + // repeated entries for species with stoichiometric + // coefficients > 1. This allows the reaction to be defined + // with unity reaction order for each reactant, and so the + // faster method 'multiply' can be used to compute the rate of + // progress instead of 'power'. + + // Note that this procedure is used for global reactions also. + // The + vector_int rk; int nr = r.reactants.size(); for (n = 0; n < nr; n++) { @@ -497,15 +529,17 @@ namespace Cantera { m_kdata->m_rkcn.push_back(0.0); - m_reactantStoich.add( reactionNumber(), rk); + m_rxnstoich.add( reactionNumber(), r); + + //m_reactantStoich.add( reactionNumber(), rk); if (r.reversible) { - m_revProductStoich.add(reactionNumber(), pk); + // m_revProductStoich.add(reactionNumber(), pk); m_revindex.push_back(reactionNumber()); m_nrev++; } else { - m_irrevProductStoich.add(reactionNumber(), pk); + //m_irrevProductStoich.add(reactionNumber(), pk); m_irrev.push_back( reactionNumber() ); m_nirrev++; } diff --git a/Cantera/src/InterfaceKinetics.h b/Cantera/src/InterfaceKinetics.h index 5f310b5ad..c8df81f45 100644 --- a/Cantera/src/InterfaceKinetics.h +++ b/Cantera/src/InterfaceKinetics.h @@ -22,7 +22,8 @@ #include "utilities.h" #include "RateCoeffMgr.h" -#include "StoichManager.h" +#include "ReactionStoichMgr.h" +//#include "StoichManager.h" namespace Cantera { @@ -165,13 +166,15 @@ namespace Cantera { */ virtual void getCreationRates(doublereal* cdot) { updateROP(); - fill(cdot, cdot + m_kk, 0.0); - m_revProductStoich.incrementSpecies( - m_kdata->m_ropf.begin(), cdot); - m_irrevProductStoich.incrementSpecies( - m_kdata->m_ropf.begin(), cdot); - m_reactantStoich.incrementSpecies( - m_kdata->m_ropr.begin(), cdot); + m_rxnstoich.getCreationRates(m_kk, m_kdata->m_ropf.begin(), + m_kdata->m_ropr.begin(), cdot); + //fill(cdot, cdot + m_kk, 0.0); + //m_revProductStoich.incrementSpecies( + // m_kdata->m_ropf.begin(), cdot); + //m_irrevProductStoich.incrementSpecies( + // m_kdata->m_ropf.begin(), cdot); + //m_reactantStoich.incrementSpecies( + // m_kdata->m_ropr.begin(), cdot); } /** @@ -184,11 +187,13 @@ namespace Cantera { */ virtual void getDestructionRates(doublereal* ddot) { updateROP(); - fill(ddot, ddot + m_kk, 0.0); - m_revProductStoich.incrementSpecies( - m_kdata->m_ropr.begin(), ddot); - m_reactantStoich.incrementSpecies( - m_kdata->m_ropf.begin(), ddot); + m_rxnstoich.getDestructionRates(m_kk, m_kdata->m_ropf.begin(), + m_kdata->m_ropr.begin(), ddot); + // fill(ddot, ddot + m_kk, 0.0); + //m_revProductStoich.incrementSpecies( + // m_kdata->m_ropr.begin(), ddot); + //m_reactantStoich.incrementSpecies( + // m_kdata->m_ropf.begin(), ddot); } /** @@ -200,13 +205,14 @@ namespace Cantera { */ virtual void getNetProductionRates(doublereal* net) { updateROP(); - fill(net, net + m_kk, 0.0); - m_revProductStoich.incrementSpecies( - m_kdata->m_ropnet.begin(), net); - m_irrevProductStoich.incrementSpecies( - m_kdata->m_ropnet.begin(), net); - m_reactantStoich.decrementSpecies( - m_kdata->m_ropnet.begin(), net); + //fill(net, net + m_kk, 0.0); + //m_revProductStoich.incrementSpecies( + // m_kdata->m_ropnet.begin(), net); + //m_irrevProductStoich.incrementSpecies( + // m_kdata->m_ropnet.begin(), net); + //m_reactantStoich.decrementSpecies( + // m_kdata->m_ropnet.begin(), net); + m_rxnstoich.getNetProductionRates(m_kk, m_kdata->m_ropnet.begin(), net); } //@} @@ -331,11 +337,12 @@ namespace Cantera { vector m_irrev; - StoichManagerN m_reactantStoich; - StoichManagerN m_revProductStoich; - StoichManagerN m_irrevProductStoich; + // StoichManagerN m_reactantStoich; + //StoichManagerN m_revProductStoich; + //StoichManagerN m_irrevProductStoich; - StoichManagerN m_globalReactantStoich; + //StoichManagerN m_globalReactantStoich; + ReactionStoichMgr m_rxnstoich; int m_nirrev; @@ -401,6 +408,7 @@ namespace Cantera { } void applyButlerVolmerCorrection(doublereal* kf); bool m_finalized; + bool m_has_coverage_dependence; }; } diff --git a/Cantera/src/Phase.cpp b/Cantera/src/Phase.cpp index 32f8efda5..3ab16468b 100755 --- a/Cantera/src/Phase.cpp +++ b/Cantera/src/Phase.cpp @@ -23,7 +23,7 @@ namespace Cantera { void Phase::saveState(int lenstate, doublereal* state) const { state[0] = temperature(); state[1] = density(); - getMassFractions(lenstate - 2, state + 2); + getMassFractions(state + 2); } void Phase::restoreState(vector_fp& state) { diff --git a/Cantera/src/ReactionData.h b/Cantera/src/ReactionData.h index 7c2d8d0cb..457c1f9c9 100755 --- a/Cantera/src/ReactionData.h +++ b/Cantera/src/ReactionData.h @@ -33,6 +33,7 @@ namespace Cantera { error = 0; equation = ""; default_3b_eff = 1.0; + global = false; } ~ReactionData(){} @@ -56,6 +57,7 @@ namespace Cantera { string equation; doublereal default_3b_eff; vector_fp cov; + bool global; }; } diff --git a/Cantera/src/ReactionStoichMgr.cpp b/Cantera/src/ReactionStoichMgr.cpp index 37adee654..40c81f580 100644 --- a/Cantera/src/ReactionStoichMgr.cpp +++ b/Cantera/src/ReactionStoichMgr.cpp @@ -1,3 +1,14 @@ +//------------------------------------------------ +/// +/// @file ReactionStoichMgr.cpp +/// +/// +//------------------------------------------------ + +// $Author$ +// $Revision$ +// $Date$ + // turn off warnings under Windows #ifdef WIN32 #pragma warning(disable:4786) @@ -7,6 +18,9 @@ #include "ReactionStoichMgr.h" #include "StoichManager.h" +#include "ctexceptions.h" +#include "diagnostics.h" +#include "ReactionData.h" namespace Cantera { @@ -18,6 +32,7 @@ namespace Cantera { m_reactants = new StoichManagerN; m_revproducts = new StoichManagerN; m_irrevproducts = new StoichManagerN; + m_global = new StoichManagerN; } // delete the three stoichiometry managers @@ -25,37 +40,115 @@ namespace Cantera { delete m_reactants; delete m_revproducts; delete m_irrevproducts; + delete m_global; } void ReactionStoichMgr:: add(int rxn, const vector_int& reactants, const vector_int& products, bool reversible) { - vector_fp forder(reactants.size(), 1.0); - add(rxn, reactants, products, reversible, forder); - } + m_reactants->add(rxn, reactants); - void ReactionStoichMgr:: - add(int rxn, const vector_int& reactants, const vector_int& products, - bool reversible, const vector_fp& fwdOrder) { - - // add the reactants with the specified forward order - m_reactants->add(rxn, reactants, fwdOrder); - - // depending on whether the reversible flag is set or not, add the - // products either to the reversible or irreversible product - // stoichiometry manager. if (reversible) m_revproducts->add(rxn, products); else m_irrevproducts->add(rxn, products); } + void ReactionStoichMgr:: + add(int rxn, const ReactionData& r) { + + vector_int rk; + int n, ns, m, nr = r.reactants.size(); + for (n = 0; n < nr; n++) { + ns = r.rstoich[n]; + // m_rrxn[r.reactants[n]][rnum] = ns; + for (m = 0; m < ns; m++) { + rk.push_back(r.reactants[n]); + } + } + + vector_int pk; + int np = r.products.size(); + for (n = 0; n < np; n++) { + ns = r.pstoich[n]; + // m_prxn[r.products[n]][rnum] = ns; + for (m = 0; m < ns; m++) { + pk.push_back(r.products[n]); + } + } + + m_reactants->add( rxn, rk); + + if (r.reversible) { + m_revproducts->add(rxn, pk); + } + else { + m_irrevproducts->add(rxn, pk); + } + + if (r.global) { + vector_fp delta_order(nr,0.0); + for (n = 0; n < nr; n++) { + delta_order[n] = r.order[n] - r.rstoich[n]; + cout << "rxn stoich " << r.reactants[n] << " " << r.order[n] << " " << delta_order[n] << endl; + } + m_global->add(rxn, r.reactants, delta_order); + } + } + + +// void ReactionStoichMgr:: +// add(int rxn, const vector_int& reactants, const vector_int& products, +// bool reversible, const vector_fp& fwdOrder) { + +// #ifdef DIAGNOSE_RXNSTOICHMGR +// printf("ReactionStoichMgr::add adding reaction number %d\n", rxn); +// #endif + +// // add the reactants with the specified forward order +// m_reactants->add(rxn, reactants, fwdOrder); + +// // check whether any orders are not equal to 1.0 +// // if so, add an entry to the global stoich manager +// // with orders decremented by one +// int nr = reactants.size(); +// int n; +// bool global = false; +// for (n = 0; n < nr; n++) { +// if (fwdOrder[n] != 1.0) { + +// #ifdef DIAGNOSE_RXNSTOICHMGR +// printf(".... global reaction: fwdOrder[%d] = %f.\n", n, fwdOrder[n]); +// #endif +// global = true; +// if (reversible) { +// throw CanteraError("ReactionStoichMgr::add", +// "reversible global reactions not allowed"); +// } +// } +// } +// if (global) { +// vector_fp fwdOrder_minus_one; +// for (n = 0; n < nr; n++) +// fwdOrder_minus_one[n] = fwdOrder[n] - 1.0; +// m_global->add(rxn, reactants, fwdOrder_minus_one); +// } + + +// // depending on whether the reversible flag is set or not, add the +// // products either to the reversible or irreversible product +// // stoichiometry manager. +// if (reversible) +// m_revproducts->add(rxn, products); +// else +// m_irrevproducts->add(rxn, products); +// } + void ReactionStoichMgr:: getCreationRates(int nsp, const doublereal* ropf, const doublereal* ropr, doublereal* c) { - // zero out the target array fill(c, c + nsp, 0.0); m_revproducts->incrementSpecies(ropf, c); m_irrevproducts->incrementSpecies(ropf, c); @@ -95,6 +188,7 @@ namespace Cantera { void ReactionStoichMgr:: multiplyReactants(const doublereal* c, doublereal* r) { m_reactants->multiply(c, r); + m_global->power(c, r); } void ReactionStoichMgr:: diff --git a/Cantera/src/ReactionStoichMgr.h b/Cantera/src/ReactionStoichMgr.h index 9b9e23a0c..7eb667691 100644 --- a/Cantera/src/ReactionStoichMgr.h +++ b/Cantera/src/ReactionStoichMgr.h @@ -18,6 +18,7 @@ namespace Cantera { class StoichManagerN; + class ReactionData; /** * Reaction mechanism stoichiometry manager. This is an internal class used @@ -110,11 +111,11 @@ namespace Cantera { * species with index in the corresponding location in 'reactants.' * */ - void add(int rxn, const vector_int& reactants, const vector_int& products, - bool reversible, const vector_fp& fwdOrder); - + // void add(int rxn, const vector_int& reactants, const vector_int& products, + // bool reversible, const vector_fp& fwdOrder); + void add(int rxn, const ReactionData& r); /** * Species creation rates. @@ -219,9 +220,10 @@ namespace Cantera { protected: - StoichManagerN* m_reactants; - StoichManagerN* m_revproducts; - StoichManagerN* m_irrevproducts; + StoichManagerN* m_reactants; + StoichManagerN* m_revproducts; + StoichManagerN* m_irrevproducts; + StoichManagerN* m_global; }; } diff --git a/Cantera/src/State.cpp b/Cantera/src/State.cpp index 6a89a3b5c..fee3d5765 100644 --- a/Cantera/src/State.cpp +++ b/Cantera/src/State.cpp @@ -116,6 +116,26 @@ namespace Cantera { } } + void State::getConcentrations(doublereal* c) const { + scale(m_ym.begin(), m_ym.end(), c, m_dens); + } + + doublereal State::mean_X(const doublereal* Q) const { + return m_mmw*dot(m_ym.begin(), m_ym.end(), Q); + } + + doublereal State::mean_Y(const doublereal* Q) const { + return dot(m_ym.begin(), m_ym.end(), Q); + } + + void State::getMoleFractions(doublereal* x) const { + scale(m_ym.begin(), m_ym.end(), x, m_mmw); + } + + void State::getMassFractions(doublereal* y) const { + copy(m_y.begin(), m_y.end(), y); + } + void State::init(const array_fp& mw) { m_kk = mw.size(); m_molwts.resize(m_kk); diff --git a/Cantera/src/State.h b/Cantera/src/State.h index 9d13956aa..0563417fa 100755 --- a/Cantera/src/State.h +++ b/Cantera/src/State.h @@ -19,8 +19,9 @@ #ifndef CT_STATE2_H #define CT_STATE2_H -#include "utilities.h" -#include "ctexceptions.h" +#include "ct_defs.h" +//#include "utilities.h" +//#include "ctexceptions.h" namespace Cantera { @@ -62,9 +63,7 @@ namespace Cantera { * @param x On return, x contains the mole fractions. Must have a * length greater than or equal to the number of species. */ - void getMoleFractions(doublereal* x) const { - scale(m_ym.begin(), m_ym.end(), x, m_mmw); - } + void getMoleFractions(doublereal* x) const; /// The mole fraction of species k. doublereal moleFraction(int k) const; @@ -89,18 +88,16 @@ namespace Cantera { * @param y On return, y contains the mass fractions. Array y * must have a length at least as large as the number of species. */ - void getMassFractions(size_t leny, doublereal* y) const { - copy(m_y.begin(), m_y.end(), y); - } + //void getMassFractions(size_t leny, doublereal* y) const { + // copy(m_y.begin(), m_y.end(), y); + //} /** * Get the species mass fractions. @param y On return, y * contains the mass fractions. Array \i y must have a length * greater than or equal to the number of species. */ - void getMassFractions(doublereal* y) const { - copy(m_y.begin(), m_y.end(), y); - } + void getMassFractions(doublereal* y) const; /// Mass fraction of species k. doublereal massFraction(int k) const; @@ -127,10 +124,7 @@ namespace Cantera { * Array \i c must have a length greater than or equal to * the number of species. */ - void getConcentrations(doublereal* c) const { - doublereal f = m_dens; - scale(m_ym.begin(), m_ym.end(), c, f); - } + void getConcentrations(doublereal* c) const; /** * Evaluate the mole-fraction-weighted mean of Q: @@ -138,10 +132,7 @@ namespace Cantera { * Array Q should contain pure-species molar property * values. */ - doublereal mean_X(const doublereal* Q) const { - return m_mmw*dot(m_ym.begin(), m_ym.end(), Q); - } - + doublereal mean_X(const doublereal* Q) const; /** * Evaluate the mass-fraction-weighted mean of Q: @@ -149,9 +140,7 @@ namespace Cantera { * Array Q should contain pure-species property * values in mass units. */ - doublereal mean_Y(const doublereal* Q) const { - return dot(m_y.begin(), m_y.end(), Q); - } + doublereal mean_Y(const doublereal* Q) const; /** * The mean molecular weight. Units: (kg/kmol) @@ -179,10 +168,10 @@ namespace Cantera { /// Set the density to value rho (kg/m^3). void setDensity(doublereal rho) { - if (rho != m_dens) { + // if (rho != m_dens) { m_dens = rho; //m_C_updater.need_update(); - } + //} } /// Set the molar density to value n (kmol/m^3). diff --git a/Cantera/src/SurfPhase.cpp b/Cantera/src/SurfPhase.cpp index c2a7f0666..5d50faaa6 100644 --- a/Cantera/src/SurfPhase.cpp +++ b/Cantera/src/SurfPhase.cpp @@ -14,6 +14,7 @@ #include "SurfPhase.h" #include "EdgePhase.h" +#include "utilities.h" #include using namespace std; @@ -225,19 +226,18 @@ namespace Cantera { m_tlast = tnow; doublereal rt = GasConstant * tnow; int k; - //doublereal deltaE; for (k = 0; k < m_kk; k++) { m_h0[k] *= rt; m_s0[k] *= GasConstant; m_cp0[k] *= GasConstant; - //deltaE = m_pe[k]; - //m_h0[k] += deltaE; m_mu0[k] = m_h0[k] - tnow*m_s0[k]; } m_tlast = tnow; } } + + EdgePhase::EdgePhase(doublereal n0) : SurfPhase(n0) { setNDim(1); } diff --git a/Cantera/src/ThermoPhase.cpp b/Cantera/src/ThermoPhase.cpp index 811849aab..5f1ebf895 100644 --- a/Cantera/src/ThermoPhase.cpp +++ b/Cantera/src/ThermoPhase.cpp @@ -23,6 +23,13 @@ namespace Cantera { + void ThermoPhase::getActivities(doublereal* a) { + getActivityConcentrations(a); + int nsp = nSpecies(); + int k; + for (k = 0; k < nsp; k++) a[k] /= standardConcentration(k); + } + void ThermoPhase::setState_TPX(doublereal t, doublereal p, const doublereal* x) { setMoleFractions(x); setTemperature(t); setPressure(p); diff --git a/Cantera/src/ThermoPhase.h b/Cantera/src/ThermoPhase.h index 6dbcfbb1c..d822c8625 100755 --- a/Cantera/src/ThermoPhase.h +++ b/Cantera/src/ThermoPhase.h @@ -550,12 +550,7 @@ namespace Cantera { } //@} - void getActivities(doublereal* a) { - getActivityConcentrations(a); - int nsp = nSpecies(); - doublereal rc = 1.0/standardConcentration(); - scale(a, a + nsp, a, rc); - } + void getActivities(doublereal* a); /** diff --git a/Cantera/src/ctexceptions.h b/Cantera/src/ctexceptions.h index d0331b6be..23d985e28 100755 --- a/Cantera/src/ctexceptions.h +++ b/Cantera/src/ctexceptions.h @@ -24,17 +24,17 @@ namespace Cantera { CanteraError() {} CanteraError(string proc, string msg) { setError(proc, msg); - m_msg = msg; + //m_msg = msg; } virtual ~CanteraError(){} - string errorMessage() { return m_msg; } - void append(string msg) { m_msg += msg; } - void saveError(string procedure) { - setError(procedure, m_msg); - m_msg = ""; - } + //string errorMessage() { return m_msg; } + //void append(string msg) { m_msg += msg; } + //void saveError(string procedure) { + // setError(procedure, m_msg); + // m_msg = ""; + //} protected: - string m_msg; + //string m_msg; }; class ArraySizeError : public CanteraError { @@ -42,6 +42,15 @@ namespace Cantera { ArraySizeError(string proc, int sz, int reqd) : CanteraError(proc, "Array size ("+int2str(sz)+") too small. Must be at least "+int2str(reqd)) {} }; + + class ElementRangeError : public CanteraError { + public: + ElementRangeError(string func, int m, int mmax) : + CanteraError(func, "Element index " + int2str(m) + + " outside valid range of 0 to " + int2str(mmax-1)) {} + }; + + } #endif diff --git a/Cantera/src/diagnostics.h b/Cantera/src/diagnostics.h new file mode 100644 index 000000000..f898af8f4 --- /dev/null +++ b/Cantera/src/diagnostics.h @@ -0,0 +1,10 @@ +#ifndef CT_DIAGNOSTICS_H +#define CT_DIAGNOSTICS_H + +#ifdef DIAGNOSE_ALL + +#define DIAGNOSE_RXNSTOICHMGR // ReactionStoichMgr + +#endif + +#endif diff --git a/Cantera/src/importCTML.cpp b/Cantera/src/importCTML.cpp index 39310256b..734be8ff1 100755 --- a/Cantera/src/importCTML.cpp +++ b/Cantera/src/importCTML.cpp @@ -1146,7 +1146,7 @@ next: throw CanteraError("installReaction", "reaction orders may only be given for " "irreversible reactions"); - //typ = "global"; + rdata.global = true; } diff --git a/Cantera/src/oneD/StFlow.cpp b/Cantera/src/oneD/StFlow.cpp index c62f1a6cd..05227acc7 100644 --- a/Cantera/src/oneD/StFlow.cpp +++ b/Cantera/src/oneD/StFlow.cpp @@ -64,7 +64,7 @@ namespace Cantera { // copy mass fractions int nsp0 = oldmech.nSpecies(); - int nsp1 = newmech.nSpecies(); + //int nsp1 = newmech.nSpecies(); // loop over the species in the old mechanism for (int k = 0; k < nsp0; k++) { @@ -86,17 +86,17 @@ namespace Cantera { // normalize mass fractions for (j = 0; j < points; j++) { newmech.setMassFractions(&newSoln[nv_new*j + 4]); - newmech.getMassFractions(nsp1,&newSoln[nv_new*j + 4]); + newmech.getMassFractions(&newSoln[nv_new*j + 4]); } } //---------------------- drawline ---------------------------------- - static void drawline(ostream& s) { - s << "\n-------------------------------------" - << "------------------------------------------"; - } + //static void drawline(ostream& s) { + // s << "\n-------------------------------------" + // << "------------------------------------------"; + //} static void drawline() { writelog("\n-------------------------------------" diff --git a/Cantera/src/oneD/boundaries1D.cpp b/Cantera/src/oneD/boundaries1D.cpp index 8908aa6b8..ffc0269c3 100644 --- a/Cantera/src/oneD/boundaries1D.cpp +++ b/Cantera/src/oneD/boundaries1D.cpp @@ -638,7 +638,7 @@ namespace Cantera { //scale(m_work.begin(), m_work.end(), m_work.begin(), m_mult[0]); - bool enabled = true; + // bool enabled = true; int ioffset = m_kin->kineticsSpeciesIndex(0, m_surfindex); if (m_enabled) { diff --git a/Cantera/src/zeroD/Reactor.cpp b/Cantera/src/zeroD/Reactor.cpp index 28d67c52f..ae5e09324 100644 --- a/Cantera/src/zeroD/Reactor.cpp +++ b/Cantera/src/zeroD/Reactor.cpp @@ -68,7 +68,7 @@ namespace Cantera { // set components y + 2 ... y + K + 1 to the // mass M_k of each species - m_mix->getMassFractions(leny-2, y+2); + m_mix->getMassFractions(y+2); scale(y + 2, y + m_nsp + 2, y + 2, mass); // set the first component to the total internal