1 /** 2 * @file Kinetics.h 3 * Base class for kinetics managers and also contains the kineticsmgr 4 * module documentation (see \ref kineticsmgr and class 5 * \link Cantera::Kinetics Kinetics\endlink). 6 */ 7 8 // This file is part of Cantera. See License.txt in the top-level directory or 9 // at https://cantera.org/license.txt for license and copyright information. 10 11 #ifndef CT_KINETICS_H 12 #define CT_KINETICS_H 13 14 #include "StoichManager.h" 15 #include "cantera/base/ValueCache.h" 16 #include "cantera/kinetics/ReactionFactory.h" 17 18 namespace Cantera 19 { 20 21 class ThermoPhase; 22 class Reaction; 23 class Solution; 24 class AnyMap; 25 26 /** 27 * @defgroup chemkinetics Chemical Kinetics 28 */ 29 30 /// @defgroup kineticsmgr Kinetics Managers 31 /// @section kinmodman Models and Managers 32 /// 33 /// A kinetics manager is a C++ class that implements a kinetics model; a 34 /// kinetics model is a set of mathematical equation describing how various 35 /// kinetic quantities are to be computed -- reaction rates, species production 36 /// rates, etc. Many different kinetics models might be defined to handle 37 /// different types of kinetic processes. For example, one kinetics model might 38 /// use expressions valid for elementary reactions in ideal gas mixtures. It 39 /// might, for example, require the reaction orders to be integral and equal to 40 /// the forward stoichiometric coefficients, require that each reaction be 41 /// reversible with a reverse rate satisfying detailed balance, include 42 /// pressure-dependent unimolecular reactions, etc. Another kinetics model might 43 /// be designed for heterogeneous chemistry at interfaces, and might allow 44 /// empirical reaction orders, coverage-dependent activation energies, 45 /// irreversible reactions, and include effects of potential differences across 46 /// the interface on reaction rates. 47 /// 48 /// A kinetics manager implements a kinetics model. Since the model equations 49 /// may be complex and expensive to evaluate, a kinetics manager may adopt 50 /// various strategies to 'manage' the computation and evaluate the expressions 51 /// efficiently. For example, if there are rate coefficients or other quantities 52 /// that depend only on temperature, a manager class may choose to store these 53 /// quantities internally, and re-evaluate them only when the temperature has 54 /// actually changed. Or a manager designed for use with reaction mechanisms 55 /// with a few repeated activation energies might precompute the terms \f$ 56 /// exp(-E/RT) \f$, instead of evaluating the exponential repeatedly for each 57 /// reaction. There are many other possible 'management styles', each of which 58 /// might be better suited to some reaction mechanisms than others. 59 /// 60 /// But however a manager structures the internal computation, the tasks the 61 /// manager class must perform are, for the most part, the same. It must be able 62 /// to compute reaction rates, species production rates, equilibrium constants, 63 /// etc. Therefore, all kinetics manager classes should have a common set of 64 /// public methods, but differ in how they implement these methods. 65 /// 66 /// A kinetics manager computes reaction rates of progress, species production 67 /// rates, equilibrium constants, and similar quantities for a reaction 68 /// mechanism. All kinetics manager classes derive from class Kinetics, which 69 /// defines a common public interface for all kinetics managers. Each derived 70 /// class overloads the virtual methods of Kinetics to implement a particular 71 /// kinetics model. 72 /// 73 /// For example, class GasKinetics implements reaction rate expressions 74 /// appropriate for homogeneous reactions in ideal gas mixtures, and class 75 /// InterfaceKinetics implements expressions appropriate for heterogeneous 76 /// mechanisms at interfaces, including how to handle reactions involving 77 /// charged species of phases with different electric potentials --- something 78 /// that class GasKinetics doesn't deal with at all. 79 /// 80 /// Many of the methods of class Kinetics write into arrays the values of some 81 /// quantity for each species, for example the net production rate. These 82 /// methods always write the results into flat arrays, ordered by phase in the 83 /// order the phase was added, and within a phase in the order the species were 84 /// added to the phase (which is the same ordering as in the input file). 85 /// Example: suppose a heterogeneous mechanism involves three phases -- a bulk 86 /// phase 'a', another bulk phase 'b', and the surface phase 'a:b' at the a/b 87 /// interface. Phase 'a' contains 12 species, phase 'b' contains 3, and at the 88 /// interface there are 5 adsorbed species defined in phase 'a:b'. Then methods 89 /// like getNetProductionRates(doublereal* net) will write and output array of 90 /// length 20, beginning at the location pointed to by 'net'. The first 12 91 /// values will be the net production rates for all 12 species of phase 'a' 92 /// (even if some do not participate in the reactions), the next 3 will be for 93 /// phase 'b', and finally the net production rates for the surface species will 94 /// occupy the last 5 locations. 95 /// @ingroup chemkinetics 96 97 98 //! Public interface for kinetics managers. 99 /*! 100 * This class serves as a base class to derive 'kinetics managers', which are 101 * classes that manage homogeneous chemistry within one phase, or heterogeneous 102 * chemistry at one interface. The virtual methods of this class are meant to be 103 * overloaded in subclasses. The non-virtual methods perform generic functions 104 * and are implemented in Kinetics. They should not be overloaded. Only those 105 * methods required by a subclass need to be overloaded; the rest will throw 106 * exceptions if called. 107 * 108 * When the nomenclature "kinetics species index" is used below, this means that 109 * the species index ranges over all species in all phases handled by the 110 * kinetics manager. 111 * 112 * @ingroup kineticsmgr 113 */ 114 class Kinetics 115 { 116 public: 117 /** 118 * @name Constructors and General Information about Mechanism 119 */ 120 //@{ 121 122 /// Default constructor. 123 Kinetics(); 124 125 virtual ~Kinetics(); 126 127 //! Kinetics objects are not copyable or assignable 128 Kinetics(const Kinetics&) = delete; 129 Kinetics& operator=(const Kinetics&)= delete; 130 131 //! Identifies the Kinetics manager type. 132 //! Each class derived from Kinetics should override this method to return 133 //! a meaningful identifier. kineticsType()134 virtual std::string kineticsType() const { 135 return "Kinetics"; 136 } 137 138 //! Finalize Kinetics object and associated objects 139 virtual void resizeReactions(); 140 141 //! Number of reactions in the reaction mechanism. nReactions()142 size_t nReactions() const { 143 return m_reactions.size(); 144 } 145 146 //! Check that the specified reaction index is in range 147 //! Throws an exception if i is greater than nReactions() 148 void checkReactionIndex(size_t m) const; 149 150 //! Check that an array size is at least nReactions() 151 //! Throws an exception if ii is less than nReactions(). Used before calls 152 //! which take an array pointer. 153 void checkReactionArraySize(size_t ii) const; 154 155 //! Check that the specified species index is in range 156 //! Throws an exception if k is greater than nSpecies()-1 157 void checkSpeciesIndex(size_t k) const; 158 159 //! Check that an array size is at least nSpecies() 160 //! Throws an exception if kk is less than nSpecies(). Used before calls 161 //! which take an array pointer. 162 void checkSpeciesArraySize(size_t mm) const; 163 164 //@} 165 //! @name Information/Lookup Functions about Phases and Species 166 //@{ 167 168 /** 169 * The number of phases participating in the reaction mechanism. For a 170 * homogeneous reaction mechanism, this will always return 1, but for a 171 * heterogeneous mechanism it will return the total number of phases in the 172 * mechanism. 173 */ nPhases()174 size_t nPhases() const { 175 return m_thermo.size(); 176 } 177 178 //! Check that the specified phase index is in range 179 //! Throws an exception if m is greater than nPhases() 180 void checkPhaseIndex(size_t m) const; 181 182 //! Check that an array size is at least nPhases() 183 //! Throws an exception if mm is less than nPhases(). Used before calls 184 //! which take an array pointer. 185 void checkPhaseArraySize(size_t mm) const; 186 187 /** 188 * Return the phase index of a phase in the list of phases defined within 189 * the object. 190 * 191 * @param ph std::string name of the phase 192 * 193 * If a -1 is returned, then the phase is not defined in the Kinetics 194 * object. 195 */ phaseIndex(const std::string & ph)196 size_t phaseIndex(const std::string& ph) const { 197 if (m_phaseindex.find(ph) == m_phaseindex.end()) { 198 return npos; 199 } else { 200 return m_phaseindex.at(ph) - 1; 201 } 202 } 203 204 /** 205 * This returns the integer index of the phase which has ThermoPhase type 206 * cSurf. For heterogeneous mechanisms, this identifies the one surface 207 * phase. For homogeneous mechanisms, this returns -1. 208 */ surfacePhaseIndex()209 size_t surfacePhaseIndex() const { 210 return m_surfphase; 211 } 212 213 /** 214 * Phase where the reactions occur. For heterogeneous mechanisms, one of 215 * the phases in the list of phases represents the 2D interface or 1D edge 216 * at which the reactions take place. This method returns the index of the 217 * phase with the smallest spatial dimension (1, 2, or 3) among the list 218 * of phases. If there is more than one, the index of the first one is 219 * returned. For homogeneous mechanisms, the value 0 is returned. 220 */ reactionPhaseIndex()221 size_t reactionPhaseIndex() const { 222 return m_rxnphase; 223 } 224 225 /** 226 * This method returns a reference to the nth ThermoPhase object defined 227 * in this kinetics mechanism. It is typically used so that member 228 * functions of the ThermoPhase object may be called. For homogeneous 229 * mechanisms, there is only one object, and this method can be called 230 * without an argument to access it. 231 * 232 * @param n Index of the ThermoPhase being sought. 233 */ 234 ThermoPhase& thermo(size_t n=0) { 235 return *m_thermo[n]; 236 } 237 const ThermoPhase& thermo(size_t n=0) const { 238 return *m_thermo[n]; 239 } 240 241 /** 242 * The total number of species in all phases participating in the kinetics 243 * mechanism. This is useful to dimension arrays for use in calls to 244 * methods that return the species production rates, for example. 245 */ nTotalSpecies()246 size_t nTotalSpecies() const { 247 return m_kk; 248 } 249 250 /** 251 * The location of species k of phase n in species arrays. Kinetics manager 252 * classes return species production rates in flat arrays, with the species 253 * of each phases following one another, in the order the phases were added. 254 * This method is useful to find the value for a particular species of a 255 * particular phase in arrays returned from methods like getCreationRates 256 * that return an array of species-specific quantities. 257 * 258 * Example: suppose a heterogeneous mechanism involves three phases. The 259 * first contains 12 species, the second 26, and the third 3. Then species 260 * arrays must have size at least 41, and positions 0 - 11 are the values 261 * for the species in the first phase, positions 12 - 37 are the values for 262 * the species in the second phase, etc. Then kineticsSpeciesIndex(7, 0) = 263 * 7, kineticsSpeciesIndex(4, 1) = 16, and kineticsSpeciesIndex(2, 2) = 40. 264 * 265 * @param k species index 266 * @param n phase index for the species 267 */ kineticsSpeciesIndex(size_t k,size_t n)268 size_t kineticsSpeciesIndex(size_t k, size_t n) const { 269 return m_start[n] + k; 270 } 271 272 //! Return the name of the kth species in the kinetics manager. 273 /*! 274 * k is an integer from 0 to ktot - 1, where ktot is the number of 275 * species in the kinetics manager, which is the sum of the number of 276 * species in all phases participating in the kinetics manager. If k is 277 * out of bounds, the string "<unknown>" is returned. 278 * 279 * @param k species index 280 */ 281 std::string kineticsSpeciesName(size_t k) const; 282 283 /** 284 * This routine will look up a species number based on the input 285 * std::string nm. The lookup of species will occur for all phases 286 * listed in the kinetics object. 287 * 288 * return 289 * - If a match is found, the position in the species list is returned. 290 * - If no match is found, the value -1 is returned. 291 * 292 * @param nm Input string name of the species 293 */ 294 size_t kineticsSpeciesIndex(const std::string& nm) const; 295 296 /** 297 * This routine will look up a species number based on the input 298 * std::string nm. The lookup of species will occur in the specified 299 * phase of the object, or all phases if ph is "<any>". 300 * 301 * return 302 * - If a match is found, the position in the species list is returned. 303 * - If no match is found, the value npos (-1) is returned. 304 * 305 * @param nm Input string name of the species 306 * @param ph Input string name of the phase. 307 */ 308 size_t kineticsSpeciesIndex(const std::string& nm, 309 const std::string& ph) const; 310 311 /** 312 * This function looks up the name of a species and returns a 313 * reference to the ThermoPhase object of the phase where the species 314 * resides. Will throw an error if the species doesn't match. 315 * 316 * @param nm String containing the name of the species. 317 */ 318 ThermoPhase& speciesPhase(const std::string& nm); 319 const ThermoPhase& speciesPhase(const std::string& nm) const; 320 321 /** 322 * This function takes as an argument the kineticsSpecies index 323 * (i.e., the list index in the list of species in the kinetics 324 * manager) and returns the species' owning ThermoPhase object. 325 * 326 * @param k Species index 327 */ speciesPhase(size_t k)328 ThermoPhase& speciesPhase(size_t k) { 329 return thermo(speciesPhaseIndex(k)); 330 } 331 332 /** 333 * This function takes as an argument the kineticsSpecies index (i.e., the 334 * list index in the list of species in the kinetics manager) and returns 335 * the index of the phase owning the species. 336 * 337 * @param k Species index 338 */ 339 size_t speciesPhaseIndex(size_t k) const; 340 341 //! @} 342 //! @name Reaction Rates Of Progress 343 //! @{ 344 345 //! Return the forward rates of progress of the reactions 346 /*! 347 * Forward rates of progress. Return the forward rates of 348 * progress in array fwdROP, which must be dimensioned at 349 * least as large as the total number of reactions. 350 * 351 * @param fwdROP Output vector containing forward rates 352 * of progress of the reactions. Length: nReactions(). 353 */ 354 virtual void getFwdRatesOfProgress(doublereal* fwdROP); 355 356 //! Return the Reverse rates of progress of the reactions 357 /*! 358 * Return the reverse rates of progress in array revROP, which must be 359 * dimensioned at least as large as the total number of reactions. 360 * 361 * @param revROP Output vector containing reverse rates 362 * of progress of the reactions. Length: nReactions(). 363 */ 364 virtual void getRevRatesOfProgress(doublereal* revROP); 365 366 /** 367 * Net rates of progress. Return the net (forward - reverse) rates of 368 * progress in array netROP, which must be dimensioned at least as large 369 * as the total number of reactions. 370 * 371 * @param netROP Output vector of the net ROP. Length: nReactions(). 372 */ 373 virtual void getNetRatesOfProgress(doublereal* netROP); 374 375 //! Return a vector of Equilibrium constants. 376 /*! 377 * Return the equilibrium constants of the reactions in concentration 378 * units in array kc, which must be dimensioned at least as large as the 379 * total number of reactions. 380 * 381 * \f[ 382 * Kc_i = exp [ \Delta G_{ss,i} ] prod(Cs_k) exp(\sum_k \nu_{k,i} F \phi_n) ] 383 * \f] 384 * 385 * @param kc Output vector containing the equilibrium constants. 386 * Length: nReactions(). 387 */ getEquilibriumConstants(doublereal * kc)388 virtual void getEquilibriumConstants(doublereal* kc) { 389 throw NotImplementedError("Kinetics::getEquilibriumConstants"); 390 } 391 392 /** 393 * Change in species properties. Given an array of molar species property 394 * values \f$ z_k, k = 1, \dots, K \f$, return the array of reaction values 395 * \f[ 396 * \Delta Z_i = \sum_k \nu_{k,i} z_k, i = 1, \dots, I. 397 * \f] 398 * For example, if this method is called with the array of standard-state 399 * molar Gibbs free energies for the species, then the values returned in 400 * array \c deltaProperty would be the standard-state Gibbs free energies of 401 * reaction for each reaction. 402 * 403 * @param property Input vector of property value. Length: m_kk. 404 * @param deltaProperty Output vector of deltaRxn. Length: nReactions(). 405 */ 406 virtual void getReactionDelta(const doublereal* property, 407 doublereal* deltaProperty); 408 409 /** 410 * Given an array of species properties 'g', return in array 'dg' the 411 * change in this quantity in the reversible reactions. Array 'g' must 412 * have a length at least as great as the number of species, and array 413 * 'dg' must have a length as great as the total number of reactions. 414 * This method only computes 'dg' for the reversible reactions, and the 415 * entries of 'dg' for the irreversible reactions are unaltered. This is 416 * primarily designed for use in calculating reverse rate coefficients 417 * from thermochemistry for reversible reactions. 418 */ 419 virtual void getRevReactionDelta(const doublereal* g, doublereal* dg); 420 421 //! Return the vector of values for the reaction Gibbs free energy change. 422 /*! 423 * (virtual from Kinetics.h) 424 * These values depend upon the concentration of the solution. 425 * 426 * units = J kmol-1 427 * 428 * @param deltaG Output vector of deltaG's for reactions Length: 429 * nReactions(). 430 */ getDeltaGibbs(doublereal * deltaG)431 virtual void getDeltaGibbs(doublereal* deltaG) { 432 throw NotImplementedError("Kinetics::getDeltaGibbs"); 433 } 434 435 //! Return the vector of values for the reaction electrochemical free 436 //! energy change. 437 /*! 438 * These values depend upon the concentration of the solution and the 439 * voltage of the phases 440 * 441 * units = J kmol-1 442 * 443 * @param deltaM Output vector of deltaM's for reactions Length: 444 * nReactions(). 445 */ getDeltaElectrochemPotentials(doublereal * deltaM)446 virtual void getDeltaElectrochemPotentials(doublereal* deltaM) { 447 throw NotImplementedError("Kinetics::getDeltaElectrochemPotentials"); 448 } 449 450 /** 451 * Return the vector of values for the reactions change in enthalpy. 452 * These values depend upon the concentration of the solution. 453 * 454 * units = J kmol-1 455 * 456 * @param deltaH Output vector of deltaH's for reactions Length: 457 * nReactions(). 458 */ getDeltaEnthalpy(doublereal * deltaH)459 virtual void getDeltaEnthalpy(doublereal* deltaH) { 460 throw NotImplementedError("Kinetics::getDeltaEnthalpy"); 461 } 462 463 /** 464 * Return the vector of values for the reactions change in entropy. These 465 * values depend upon the concentration of the solution. 466 * 467 * units = J kmol-1 Kelvin-1 468 * 469 * @param deltaS Output vector of deltaS's for reactions Length: 470 * nReactions(). 471 */ getDeltaEntropy(doublereal * deltaS)472 virtual void getDeltaEntropy(doublereal* deltaS) { 473 throw NotImplementedError("Kinetics::getDeltaEntropy"); 474 } 475 476 /** 477 * Return the vector of values for the reaction standard state Gibbs free 478 * energy change. These values don't depend upon the concentration of the 479 * solution. 480 * 481 * units = J kmol-1 482 * 483 * @param deltaG Output vector of ss deltaG's for reactions Length: 484 * nReactions(). 485 */ getDeltaSSGibbs(doublereal * deltaG)486 virtual void getDeltaSSGibbs(doublereal* deltaG) { 487 throw NotImplementedError("Kinetics::getDeltaSSGibbs"); 488 } 489 490 /** 491 * Return the vector of values for the change in the standard state 492 * enthalpies of reaction. These values don't depend upon the concentration 493 * of the solution. 494 * 495 * units = J kmol-1 496 * 497 * @param deltaH Output vector of ss deltaH's for reactions Length: 498 * nReactions(). 499 */ getDeltaSSEnthalpy(doublereal * deltaH)500 virtual void getDeltaSSEnthalpy(doublereal* deltaH) { 501 throw NotImplementedError("Kinetics::getDeltaSSEnthalpy"); 502 } 503 504 /** 505 * Return the vector of values for the change in the standard state 506 * entropies for each reaction. These values don't depend upon the 507 * concentration of the solution. 508 * 509 * units = J kmol-1 Kelvin-1 510 * 511 * @param deltaS Output vector of ss deltaS's for reactions Length: 512 * nReactions(). 513 */ getDeltaSSEntropy(doublereal * deltaS)514 virtual void getDeltaSSEntropy(doublereal* deltaS) { 515 throw NotImplementedError("Kinetics::getDeltaSSEntropy"); 516 } 517 518 //! @} 519 //! @name Species Production Rates 520 //! @{ 521 522 /** 523 * Species creation rates [kmol/m^3/s or kmol/m^2/s]. Return the species 524 * creation rates in array cdot, which must be dimensioned at least as 525 * large as the total number of species in all phases. @see nTotalSpecies. 526 * 527 * @param cdot Output vector of creation rates. Length: m_kk. 528 */ 529 virtual void getCreationRates(doublereal* cdot); 530 531 /** 532 * Species destruction rates [kmol/m^3/s or kmol/m^2/s]. Return the species 533 * destruction rates in array ddot, which must be dimensioned at least as 534 * large as the total number of species. @see nTotalSpecies. 535 * 536 * @param ddot Output vector of destruction rates. Length: m_kk. 537 */ 538 virtual void getDestructionRates(doublereal* ddot); 539 540 /** 541 * Species net production rates [kmol/m^3/s or kmol/m^2/s]. Return the 542 * species net production rates (creation - destruction) in array wdot, 543 * which must be dimensioned at least as large as the total number of 544 * species. @see nTotalSpecies. 545 * 546 * @param wdot Output vector of net production rates. Length: m_kk. 547 */ 548 virtual void getNetProductionRates(doublereal* wdot); 549 550 //! @} 551 //! @name Reaction Mechanism Informational Query Routines 552 //! @{ 553 554 /** 555 * Stoichiometric coefficient of species k as a reactant in reaction i. 556 * 557 * @param k kinetic species index 558 * @param i reaction index 559 */ 560 virtual double reactantStoichCoeff(size_t k, size_t i) const; 561 562 /** 563 * Stoichiometric coefficient matrix for reactants. 564 */ reactantStoichCoeffs()565 Eigen::SparseMatrix<double> reactantStoichCoeffs() const { 566 return m_reactantStoich.stoichCoeffs(); 567 } 568 569 /** 570 * Stoichiometric coefficient of species k as a product in reaction i. 571 * 572 * @param k kinetic species index 573 * @param i reaction index 574 */ 575 virtual double productStoichCoeff(size_t k, size_t i) const; 576 577 /** 578 * Stoichiometric coefficient matrix for products. 579 */ productStoichCoeffs()580 Eigen::SparseMatrix<double> productStoichCoeffs() const { 581 return m_productStoich.stoichCoeffs(); 582 } 583 584 /** 585 * Stoichiometric coefficient matrix for products of reversible reactions. 586 */ revProductStoichCoeffs()587 Eigen::SparseMatrix<double> revProductStoichCoeffs() const { 588 return m_revProductStoich.stoichCoeffs(); 589 } 590 591 //! Reactant order of species k in reaction i. 592 /*! 593 * This is the nominal order of the activity concentration in 594 * determining the forward rate of progress of the reaction 595 * 596 * @param k kinetic species index 597 * @param i reaction index 598 */ reactantOrder(size_t k,size_t i)599 virtual doublereal reactantOrder(size_t k, size_t i) const { 600 throw NotImplementedError("Kinetics::reactantOrder"); 601 } 602 603 //! product Order of species k in reaction i. 604 /*! 605 * This is the nominal order of the activity concentration of species k in 606 * determining the reverse rate of progress of the reaction i 607 * 608 * For irreversible reactions, this will all be zero. 609 * 610 * @param k kinetic species index 611 * @param i reaction index 612 */ productOrder(int k,int i)613 virtual doublereal productOrder(int k, int i) const { 614 throw NotImplementedError("Kinetics::productOrder"); 615 } 616 617 //! Get the vector of activity concentrations used in the kinetics object 618 /*! 619 * @param[out] conc Vector of activity concentrations. Length is equal 620 * to the number of species in the kinetics object 621 */ getActivityConcentrations(doublereal * const conc)622 virtual void getActivityConcentrations(doublereal* const conc) { 623 throw NotImplementedError("Kinetics::getActivityConcentrations"); 624 } 625 626 /** 627 * Flag specifying the type of reaction. The legal values and their meaning 628 * are specific to the particular kinetics manager. 629 * 630 * @param i reaction index 631 * 632 * @deprecated To be changed after Cantera 2.6. 633 */ 634 virtual int reactionType(size_t i) const; 635 636 /** 637 * String specifying the type of reaction. 638 * 639 * @param i reaction index 640 */ 641 virtual std::string reactionTypeStr(size_t i) const; 642 643 /** 644 * True if reaction i has been declared to be reversible. If isReversible(i) 645 * is false, then the reverse rate of progress for reaction i is always 646 * zero. 647 * 648 * @param i reaction index 649 */ isReversible(size_t i)650 virtual bool isReversible(size_t i) { 651 throw NotImplementedError("Kinetics::isReversible"); 652 } 653 654 /** 655 * Return a string representing the reaction. 656 * 657 * @param i reaction index 658 */ 659 std::string reactionString(size_t i) const; 660 661 //! Returns a string containing the reactants side of the reaction equation. 662 std::string reactantString(size_t i) const; 663 664 //! Returns a string containing the products side of the reaction equation. 665 std::string productString(size_t i) const; 666 667 /** 668 * Return the forward rate constants 669 * 670 * The computed values include all temperature-dependent, pressure-dependent, 671 * and third body contributions. Length is the number of reactions. Units are 672 * a combination of kmol, m^3 and s, that depend on the rate expression for 673 * the reaction. 674 * 675 * @param kfwd Output vector containing the forward reaction rate 676 * constants. Length: nReactions(). 677 * 678 * @deprecated Behavior to change after Cantera 2.6; for Cantera 2.6, rate 679 * constants of three-body reactions are multiplied with third-body 680 * concentrations (no change to legacy behavior). After Cantera 2.6, 681 * results will no longer include third-body concentrations and be 682 * consistent with conventional definitions (see Eq. 9.75 in 683 * Kee, Coltrin and Glarborg, 'Chemically eacting Flow', Wiley 684 * Interscience, 2003). 685 * For new behavior, set 'Cantera::use_legacy_rate_constants(false)'. 686 */ getFwdRateConstants(double * kfwd)687 virtual void getFwdRateConstants(double* kfwd) { 688 throw NotImplementedError("Kinetics::getFwdRateConstants"); 689 } 690 691 /** 692 * Return the reverse rate constants. 693 * 694 * The computed values include all temperature-dependent, pressure-dependent, 695 * and third body contributions. Length is the number of reactions. Units are 696 * a combination of kmol, m^3 and s, that depend on the rate expression for 697 * the reaction. Note, this routine will return rate constants for 698 * irreversible reactions if the default for `doIrreversible` is overridden. 699 * 700 * @param krev Output vector of reverse rate constants 701 * @param doIrreversible boolean indicating whether irreversible reactions 702 * should be included. 703 * 704 * @deprecated Behavior to change after Cantera 2.6; for Cantera 2.6, rate 705 * constants of three-body reactions are multiplied with third-body 706 * concentrations (no change to legacy behavior). After Cantera 2.6, 707 * results will no longer include third-body concentrations and be 708 * consistent with conventional definitions (see Eq. 9.75 in 709 * Kee, Coltrin and Glarborg, 'Chemically eacting Flow', Wiley 710 * Interscience, 2003). 711 * For new behavior, set 'Cantera::use_legacy_rate_constants(false)'. 712 */ 713 virtual void getRevRateConstants(double* krev, 714 bool doIrreversible = false) { 715 throw NotImplementedError("Kinetics::getRevRateConstants"); 716 } 717 718 //! @} 719 //! @name Reaction Mechanism Construction 720 //! @{ 721 722 //! Add a phase to the kinetics manager object. 723 /*! 724 * This must be done before the function init() is called or before any 725 * reactions are input. The following fields are updated: 726 * 727 * - #m_start -> vector of integers, containing the starting position of 728 * the species for each phase in the kinetics mechanism. 729 * - #m_surfphase -> index of the surface phase. 730 * - #m_thermo -> vector of pointers to ThermoPhase phases that 731 * participate in the kinetics mechanism. 732 * - #m_phaseindex -> map containing the std::string id of each 733 * ThermoPhase phase as a key and the index of the phase within the 734 * kinetics manager object as the value. 735 * 736 * @param thermo Reference to the ThermoPhase to be added. 737 */ 738 virtual void addPhase(ThermoPhase& thermo); 739 740 /** 741 * Prepare the class for the addition of reactions, after all phases have 742 * been added. This method is called automatically when the first reaction 743 * is added. It needs to be called directly only in the degenerate case 744 * where there are no reactions. The base class method does nothing, but 745 * derived classes may use this to perform any initialization (allocating 746 * arrays, etc.) that requires knowing the phases. 747 */ init()748 virtual void init() {} 749 750 //! Return the parameters for a phase definition which are needed to 751 //! reconstruct an identical object using the newKinetics function. This 752 //! excludes the reaction definitions, which are handled separately. 753 AnyMap parameters(); 754 755 /** 756 * Resize arrays with sizes that depend on the total number of species. 757 * Automatically called before adding each Reaction and Phase. 758 */ 759 virtual void resizeSpecies(); 760 761 /** 762 * Add a single reaction to the mechanism. Derived classes should call the 763 * base class method in addition to handling their own specialized behavior. 764 * 765 * @param r Pointer to the Reaction object to be added. 766 * @param resize If `true`, resizeReactions is called after reaction is added. 767 * @return `true` if the reaction is added or `false` if it was skipped 768 */ 769 virtual bool addReaction(shared_ptr<Reaction> r, bool resize=true); 770 771 /** 772 * Modify the rate expression associated with a reaction. The 773 * stoichiometric equation, type of the reaction, reaction orders, third 774 * body efficiencies, reversibility, etc. must be unchanged. 775 * 776 * @param i Index of the reaction to be modified 777 * @param rNew Reaction with the new rate expressions 778 */ 779 virtual void modifyReaction(size_t i, shared_ptr<Reaction> rNew); 780 781 /** 782 * Return the Reaction object for reaction *i*. Changes to this object do 783 * not affect the Kinetics object until the #modifyReaction function is 784 * called. 785 */ 786 shared_ptr<Reaction> reaction(size_t i); 787 788 shared_ptr<const Reaction> reaction(size_t i) const; 789 790 //! Determine behavior when adding a new reaction that contains species not 791 //! defined in any of the phases associated with this kinetics manager. If 792 //! set to true, the reaction will silently be ignored. If false, (the 793 //! default) an exception will be raised. skipUndeclaredSpecies(bool skip)794 void skipUndeclaredSpecies(bool skip) { 795 m_skipUndeclaredSpecies = skip; 796 } skipUndeclaredSpecies()797 bool skipUndeclaredSpecies() const { 798 return m_skipUndeclaredSpecies; 799 } 800 801 //! Determine behavior when adding a new reaction that contains third-body 802 //! efficiencies for species not defined in any of the phases associated 803 //! with this kinetics manager. If set to true, the given third-body 804 //! efficiency will be ignored. If false, (the default) an exception will be 805 //! raised. skipUndeclaredThirdBodies(bool skip)806 void skipUndeclaredThirdBodies(bool skip) { 807 m_skipUndeclaredThirdBodies = skip; 808 } skipUndeclaredThirdBodies()809 bool skipUndeclaredThirdBodies() const { 810 return m_skipUndeclaredThirdBodies; 811 } 812 813 //@} 814 //! @name Altering Reaction Rates 815 /*! 816 * These methods alter reaction rates. They are designed primarily for 817 * carrying out sensitivity analysis, but may be used for any purpose 818 * requiring dynamic alteration of rate constants. For each reaction, a 819 * real-valued multiplier may be defined that multiplies the reaction rate 820 * coefficient. The multiplier may be set to zero to completely remove a 821 * reaction from the mechanism. 822 */ 823 //@{ 824 825 //! The current value of the multiplier for reaction i. 826 /*! 827 * @param i index of the reaction 828 */ multiplier(size_t i)829 doublereal multiplier(size_t i) const { 830 return m_perturb[i]; 831 } 832 833 //! Set the multiplier for reaction i to f. 834 /*! 835 * @param i index of the reaction 836 * @param f value of the multiplier. 837 */ setMultiplier(size_t i,doublereal f)838 virtual void setMultiplier(size_t i, doublereal f) { 839 m_perturb[i] = f; 840 } 841 invalidateCache()842 virtual void invalidateCache() {}; 843 844 //@} 845 846 //! Check for unmarked duplicate reactions and unmatched marked duplicates 847 /** 848 * If `throw_err` is true, then an exception will be thrown if either any 849 * unmarked duplicate reactions are found, or if any marked duplicate 850 * reactions do not have a matching duplicate reaction. If `throw_err` is 851 * false, the indices of the first pair of duplicate reactions found will be 852 * returned, or the index of the unmatched duplicate will be returned as 853 * both elements of the pair. If no unmarked duplicates or unmatched marked 854 * duplicate reactions are found, returns `(npos, npos)`. 855 */ 856 virtual std::pair<size_t, size_t> checkDuplicates(bool throw_err=true) const; 857 858 /*! 859 * Takes as input an array of properties for all species in the mechanism 860 * and copies those values belonging to a particular phase to the output 861 * array. 862 * @param data Input data array. 863 * @param phase Pointer to one of the phase objects participating in this 864 * reaction mechanism 865 * @param phase_data Output array where the values for the the specified 866 * phase are to be written. 867 */ 868 void selectPhase(const double* data, const ThermoPhase* phase, 869 double* phase_data); 870 871 //! Set root Solution holding all phase information setRoot(std::shared_ptr<Solution> root)872 virtual void setRoot(std::shared_ptr<Solution> root) { 873 m_root = root; 874 } 875 876 protected: 877 //! Cache for saved calculations within each Kinetics object. 878 ValueCache m_cache; 879 880 // Update internal rate-of-progress variables #m_ropf and #m_ropr. updateROP()881 virtual void updateROP() { 882 throw NotImplementedError("Kinetics::updateROP"); 883 } 884 885 //! Check whether `r1` and `r2` represent duplicate stoichiometries 886 //! This function returns a ratio if two reactions are duplicates of 887 //! one another, and 0.0 otherwise. 888 /*! 889 * `r1` and `r2` are maps of species key to stoichiometric coefficient, one 890 * for each reaction, where the species key is `1+k` for reactants and 891 * `-1-k` for products and `k` is the species index. 892 * 893 * @return 0.0 if the stoichiometries are not multiples of one another 894 * Otherwise, it returns the ratio of the stoichiometric coefficients. 895 * 896 * @ingroup kineticsmgr 897 */ 898 double checkDuplicateStoich(std::map<int, double>& r1, 899 std::map<int, double>& r2) const; 900 901 //! Check that the specified reaction is balanced (same number of atoms for 902 //! each element in the reactants and products). Raises an exception if the 903 //! reaction is not balanced. 904 /*! 905 * @deprecated To be removed in Cantera 2.6. 906 * Replaceable by Reaction::checkBalance. 907 */ 908 void checkReactionBalance(const Reaction& R); 909 910 //! @name Stoichiometry management 911 /*! 912 * These objects and functions handle turning reaction extents into species 913 * production rates and also handle turning thermo properties into reaction 914 * thermo properties. 915 */ 916 //@{ 917 918 //! Stoichiometry manager for the reactants for each reaction 919 StoichManagerN m_reactantStoich; 920 921 //! Stoichiometry manager for the products for each reaction 922 StoichManagerN m_productStoich; 923 924 //! Stoichiometry manager for the products of reversible reactions 925 StoichManagerN m_revProductStoich; 926 //@} 927 928 //! Boolean indicating whether Kinetics object is fully configured 929 bool m_ready; 930 931 //! The number of species in all of the phases 932 //! that participate in this kinetics mechanism. 933 size_t m_kk; 934 935 /// Vector of perturbation factors for each reaction's rate of 936 /// progress vector. It is initialized to one. 937 vector_fp m_perturb; 938 939 //! Vector of Reaction objects represented by this Kinetics manager 940 std::vector<shared_ptr<Reaction> > m_reactions; 941 942 //! m_thermo is a vector of pointers to ThermoPhase objects that are 943 //! involved with this kinetics operator 944 /*! 945 * For homogeneous kinetics applications, this vector will only have one 946 * entry. For interfacial reactions, this vector will consist of multiple 947 * entries; some of them will be surface phases, and the other ones will be 948 * bulk phases. The order that the objects are listed determines the order 949 * in which the species comprising each phase are listed in the source term 950 * vector, originating from the reaction mechanism. 951 * 952 * Note that this kinetics object doesn't own these ThermoPhase objects 953 * and is not responsible for creating or deleting them. 954 */ 955 std::vector<ThermoPhase*> m_thermo; 956 957 /** 958 * m_start is a vector of integers specifying the beginning position for the 959 * species vector for the n'th phase in the kinetics class. 960 */ 961 std::vector<size_t> m_start; 962 963 /** 964 * Mapping of the phase name to the position of the phase within the 965 * kinetics object. Positions start with the value of 1. The member 966 * function, phaseIndex() decrements by one before returning the index 967 * value, so that missing phases return -1. 968 */ 969 std::map<std::string, size_t> m_phaseindex; 970 971 //! Index in the list of phases of the one surface phase. 972 size_t m_surfphase; 973 974 //! Phase Index where reactions are assumed to be taking place 975 /*! 976 * We calculate this by assuming that the phase with the lowest 977 * dimensionality is the phase where reactions are taking place. 978 */ 979 size_t m_rxnphase; 980 981 //! number of spatial dimensions of lowest-dimensional phase. 982 size_t m_mindim; 983 984 //! Forward rate constant for each reaction 985 vector_fp m_rfn; 986 987 //! Reciprocal of the equilibrium constant in concentration units 988 vector_fp m_rkcn; 989 990 //! Forward rate-of-progress for each reaction 991 vector_fp m_ropf; 992 993 //! Reverse rate-of-progress for each reaction 994 vector_fp m_ropr; 995 996 //! Net rate-of-progress for each reaction 997 vector_fp m_ropnet; 998 999 //! The enthalpy change for each reaction to calculate Blowers-Masel rates 1000 vector_fp m_dH; 1001 1002 //! @see skipUndeclaredSpecies() 1003 bool m_skipUndeclaredSpecies; 1004 1005 //! @see skipUndeclaredThirdBodies() 1006 bool m_skipUndeclaredThirdBodies; 1007 1008 //! reference to Solution 1009 std::weak_ptr<Solution> m_root; 1010 }; 1011 1012 } 1013 1014 #endif 1015