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