Cantera  2.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 // Copyright 2001-2004 California Institute of Technology
9 
10 #ifndef CT_KINETICS_H
11 #define CT_KINETICS_H
12 
15 #include "cantera/thermo/mix_defs.h"
16 
17 namespace Cantera
18 {
19 
20 // forward references
21 class ReactionData;
22 
23 /**
24  * @defgroup chemkinetics Chemical Kinetics
25  */
26 
27 /// @defgroup kineticsmgr Kinetics Managers
28 /// @section kinmodman Models and Managers
29 ///
30 /// A kinetics manager is a C++ class that implements a kinetics
31 /// model; a kinetics model is a set of mathematical equation
32 /// describing how various kinetic quantities are to be computed --
33 /// reaction rates, species production rates, etc. Many different
34 /// kinetics models might be defined to handle different types of
35 /// kinetic processes. For example, one kinetics model might use
36 /// expressions valid for elementary reactions in ideal gas
37 /// mixtures. It might, for example, require the reaction orders
38 /// to be integral and equal to the forward stoichiometric
39 /// coefficients, require that each reaction be reversible with a
40 /// reverse rate satisfying detailed balance, include
41 /// pressure-dependent unimolecular reactions, etc. Another
42 /// kinetics model might be designed for heterogeneous chemistry
43 /// at interfaces, and might allow empirical reaction orders,
44 /// coverage-dependent activation energies, irreversible
45 /// 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
49 /// model equations may be complex and expensive to evaluate, a
50 /// kinetics manager may adopt various strategies to 'manage' the
51 /// computation and evaluate the expressions efficiently. For
52 /// example, if there are rate coefficients or other quantities
53 /// that depend only on temperature, a manager class may choose to
54 /// store these quantities internally, and re-evaluate them only
55 /// when the temperature has actually changed. Or a manager
56 /// designed for use with reaction mechanisms with a few repeated
57 /// activation energies might precompute the terms \f$ exp(-E/RT)
58 /// \f$, instead of evaluating the exponential repeatedly for each
59 /// reaction. There are many other possible 'management styles',
60 /// each of which might be better suited to some reaction
61 /// mechanisms than others.
62 ///
63 /// But however a manager structures the internal computation, the
64 /// tasks the manager class must perform are, for the most part,
65 /// the same. It must be able to compute reaction rates, species
66 /// production rates, equilibrium constants, etc. Therefore, all
67 /// kinetics manager classes should have a common set of public
68 /// methods, but differ in how they implement these methods.
69 ///
70 /// A kinetics manager computes reaction rates of progress,
71 /// species production rates, equilibrium constants, and similar
72 /// quantities for a reaction mechanism. All kinetics manager
73 /// classes derive from class Kinetics, which defines a common
74 /// public interface for all kinetics managers. Each derived class
75 /// overloads the virtual methods of Kinetics to implement a
76 /// particular kinetics model.
77 ///
78 /// For example, class GasKinetics implements reaction rate
79 /// expressions appropriate for homogeneous reactions in ideal gas
80 /// mixtures, and class InterfaceKinetics implements expressions
81 /// appropriate for heterogeneous mechanisms at interfaces,
82 /// including how to handle reactions involving charged species of
83 /// phases with different electric potentials --- something that
84 /// class GasKinetics doesn't deal with at all.
85 ///
86 /// Kinetics managers may be also created that hard-wire a
87 /// particular reaction mechanism in C++ code. This can often
88 /// result in faster performance. An example of this is the
89 /// kinetics manager GRI30_Kinetics that hard-wires the rate
90 /// expressions for the natural gas combustion mechanism GRI-3.0.
91 ///
92 /// Many of the methods of class Kinetics write into arrays the
93 /// values of some quantity for each species, for example the net
94 /// production rate. These methods always write the results into
95 /// flat arrays, ordered by phase in the order the phase was
96 /// added, and within a phase in the order the species were added
97 /// to the phase (which is the same ordering as in the input
98 /// file). Example: suppose a heterogeneous mechanism involves
99 /// three phases -- a bulk phase 'a', another bulk phase 'b', and
100 /// the surface phase 'a:b' at the a/b interface. Phase 'a'
101 /// contains 12 species, phase 'b' contains 3, and at the
102 /// interface there are 5 adsorbed species defined in phase
103 /// 'a:b'. Then methods like getNetProductionRates(doublereal* net)
104 /// will write and output array of length 20, beginning at the location
105 /// pointed to by 'net'. The first 12 values will be the net production
106 /// rates for all 12 species of phase 'a' (even if some do not participate
107 /// in the reactions), the next 3 will be for phase 'b', and finally the
108 /// net production rates for the surface species will occupy the last
109 /// 5 locations.
110 /// @ingroup chemkinetics
111 
112 
113 //! Public interface for kinetics managers.
114 /*!
115  * This class serves as a base class to derive 'kinetics
116  * managers', which are classes that manage homogeneous chemistry
117  * within one phase, or heterogeneous chemistry at one
118  * interface. The virtual methods of this class are meant to be
119  * overloaded in subclasses. The non-virtual methods perform
120  * generic functions and are implemented in Kinetics. They should
121  * not be overloaded. Only those methods required by a subclass
122  * need to be overloaded; the rest will throw exceptions if
123  * called.
124  *
125  * When the nomenclature "kinetics species index" is used below,
126  * this means that the species index ranges over all species in
127  * all phases handled by the kinetics manager.
128  *
129  * @ingroup kineticsmgr
130  */
131 class Kinetics
132 {
133 
134 public:
135  /**
136  * @name Constructors and General Information about Mechanism
137  */
138  //@{
139 
140  /// Default constructor.
141  Kinetics();
142 
143  /// Destructor.
144  virtual ~Kinetics();
145 
146  //!Copy Constructor for the %Kinetics object.
147  /*!
148  * Currently, this is not fully implemented. If called it will
149  * throw an exception.
150  */
151  Kinetics(const Kinetics&);
152 
153  //! Assignment operator
154  /*!
155  * This is NOT a virtual function.
156  *
157  * @param right Reference to %Kinetics object to be copied into the
158  * current one.
159  */
160  Kinetics& operator=(const Kinetics& right);
161 
162 
163  //! Duplication routine for objects which inherit from Kinetics
164  /*!
165  * This virtual routine can be used to duplicate %Kinetics objects
166  * inherited from %Kinetics even if the application only has
167  * a pointer to %Kinetics to work with.
168  *
169  * These routines are basically wrappers around the derived copy constructor.
170  *
171  * @param tpVector Vector of shallow pointers to ThermoPhase objects. this is the
172  * m_thermo vector within this object
173  */
174  virtual Kinetics* duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const;
175 
176  //! Reassign the shallow pointers within the %FKinetics object
177  /*!
178  * This type or routine is absolute necessary because the Kinetics object doesn't
179  * own the ThermoPhase objects. After a duplication, we need to point to different
180  * ThermoPhase objects.
181  *
182  * We check that the ThermoPhase objects are aligned in the same order and have
183  * the following identical properties to the ones that they are replacing.
184  * id()
185  * eosType()
186  * nSpecies()
187  *
188  * @param tpVector Vector of shallow pointers to ThermoPhase objects. this is the
189  * m_thermo vector within this object
190  */
191  virtual void assignShallowPointers(const std::vector<thermo_t*> & tpVector);
192 
193  //! Identifies the subclass of the Kinetics manager type.
194  /*!
195  * These are listed in mix_defs.h.
196  * @deprecated use type() instead
197  */
198  DEPRECATED(virtual int ID() const);
199 
200  //! Identifies the kinetics manager type.
201  /*!
202  * Each class derived from Kinetics should overload this method to
203  * return a unique integer. Standard values are defined in file
204  * mix_defs.h.
205  */
206  virtual int type() const;
207 
208  //! Number of reactions in the reaction mechanism.
209  size_t nReactions() const {
210  return m_ii;
211  }
212 
213  //! Check that the specified reaction index is in range
214  //! Throws an exception if i is greater than nReactions()
215  void checkReactionIndex(size_t m) const;
216 
217  //! Check that an array size is at least nReactions()
218  //! Throws an exception if ii is less than nReactions(). Used before calls
219  //! which take an array pointer.
220  void checkReactionArraySize(size_t ii) const;
221 
222  //! Check that the specified species index is in range
223  //! Throws an exception if k is greater than nSpecies()-1
224  void checkSpeciesIndex(size_t k) const;
225 
226  //! Check that an array size is at least nSpecies()
227  //! Throws an exception if kk is less than nSpecies(). Used before calls
228  //! which take an array pointer.
229  void checkSpeciesArraySize(size_t mm) const;
230  //@}
231 
232 
233  /**
234  * @name Information/Lookup Functions about Phases and Species
235  */
236  //@{
237 
238  /**
239  * The number of phases participating in the reaction
240  * mechanism. For a homogeneous reaction mechanism, this will
241  * always return 1, but for a heterogeneous mechanism it will
242  * return the total number of phases in the mechanism.
243  */
244  size_t nPhases() const {
245  return m_thermo.size();
246  }
247 
248  //! Check that the specified phase index is in range
249  //! Throws an exception if m is greater than nPhases()
250  void checkPhaseIndex(size_t m) const;
251 
252  //! Check that an array size is at least nPhases()
253  //! Throws an exception if mm is less than nPhases(). Used before calls
254  //! which take an array pointer.
255  void checkPhaseArraySize(size_t mm) const;
256 
257  /**
258  * Return the phase index of a phase in the list of phases
259  * defined within the object.
260  *
261  * @param ph std::string name of the phase
262  *
263  * If a -1 is returned, then the phase is not defined in
264  * the Kinetics object.
265  */
266  size_t phaseIndex(std::string ph) {
267  if (m_phaseindex.find(ph) == m_phaseindex.end()) {
268  return npos;
269  } else {
270  return m_phaseindex[ph] - 1;
271  }
272  }
273 
274  /**
275  * This returns the integer index of the phase which has
276  * ThermoPhase type cSurf. For heterogeneous mechanisms, this
277  * identifies the one surface phase. For homogeneous
278  * mechanisms, this returns -1.
279  */
280  size_t surfacePhaseIndex() {
281  return m_surfphase;
282  }
283 
284  /**
285  * Phase where the reactions occur. For heterogeneous
286  * mechanisms, one of the phases in the list of phases
287  * represents the 2D interface or 1D edge at which the
288  * reactions take place. This method returns the index of the
289  * phase with the smallest spatial dimension (1, 2, or 3)
290  * among the list of phases. If there is more than one, the
291  * index of the first one is returned. For homogeneous
292  * mechanisms, the value 0 is returned.
293  */
295  return m_rxnphase;
296  }
297 
298 
299  /**
300  * This method returns a reference to the nth ThermoPhase
301  * object defined in this kinetics mechanism. It is typically
302  * used so that member functions of the ThermoPhase object may
303  * be called. For homogeneous mechanisms, there is only one
304  * object, and this method can be called without an argument
305  * to access it.
306  *
307  * @param n Index of the ThermoPhase being sought.
308  */
309  thermo_t& thermo(size_t n=0) {
310  return *m_thermo[n];
311  }
312  const thermo_t& thermo(size_t n=0) const {
313  return *m_thermo[n];
314  }
315 
316  /**
317  * This method returns a reference to the nth ThermoPhase
318  * defined in this kinetics mechanism.
319  * It is typically used so that member functions of the
320  * ThermoPhase may be called. @deprecated This method is redundant.
321  *
322  * @param n Index of the ThermoPhase being sought.
323  */
324  DEPRECATED(thermo_t& phase(size_t n=0)) {
325  deprecatedMethod("Kinetics","phase","thermo");
326  return *m_thermo[n];
327  }
328  /**
329  * This method returns a reference to the nth ThermoPhase
330  * defined in this kinetics mechanism.
331  * It is typically used so that member functions of the
332  * ThermoPhase may be called. @deprecated This method is redundant.
333  *
334  * @param n Index of the ThermoPhase being sought.
335  */
336  DEPRECATED(const thermo_t& phase(size_t n=0) const) {
337  deprecatedMethod("Kinetics","phase","thermo");
338  return *m_thermo[n];
339  }
340 
341  /**
342  * The total number of species in all phases participating in
343  * the kinetics mechanism. This is useful to dimension arrays
344  * for use in calls to methods that return the species
345  * production rates, for example.
346  */
347  size_t nTotalSpecies() const {
348  size_t n=0, np;
349  np = nPhases();
350  for (size_t p = 0; p < np; p++) {
351  n += thermo(p).nSpecies();
352  }
353  return n;
354  }
355 
356  /**
357  * Returns the starting index of the species in the nth phase
358  * associated with the reaction mechanism.
359  *
360  * @param n Return the index of first species in the nth phase
361  * associated with the reaction mechanism.
362  */
363  size_t start(size_t n) {
364  deprecatedMethod("Kinetics","start","kineticsSpeciesIndex(0,n)");
365  return m_start[n];
366  }
367 
368 
369  /**
370  * The location of species k of phase n in species arrays.
371  * Kinetics manager classes return species production rates in
372  * flat arrays, with the species of each phases following one
373  * another, in the order the phases were added. This method
374  * is useful to find the value for a particular species of a
375  * particular phase in arrays returned from methods like
376  * getCreationRates that return an array of species-specific
377  * quantities.
378  *
379  * Example: suppose a heterogeneous mechanism involves three
380  * phases. The first contains 12 species, the second 26, and
381  * the third 3. Then species arrays must have size at least
382  * 41, and positions 0 - 11 are the values for the species in
383  * the first phase, positions 12 - 37 are the values for the
384  * species in the second phase, etc. Then
385  * kineticsSpeciesIndex(7, 0) = 7, kineticsSpeciesIndex(4, 1)
386  * = 16, and kineticsSpeciesIndex(2, 2) = 40.
387  *
388  * @param k species index
389  * @param n phase index for the species
390  */
391  size_t kineticsSpeciesIndex(size_t k, size_t n) const {
392  return m_start[n] + k;
393  }
394 
395 
396  //! Return the std::string name of the kth species in the kinetics
397  //! manager.
398  /*!
399  * k is an integer from 0 to ktot - 1, where ktot is
400  * the number of species in the kinetics manager, which is the
401  * sum of the number of species in all phases participating in
402  * the kinetics manager. If k is out of bounds, the std::string
403  * "<unknown>" is returned.
404  *
405  * @param k species index
406  */
407  std::string kineticsSpeciesName(size_t k) const;
408 
409  /**
410  * This routine will look up a species number based on the input
411  * std::string nm. The lookup of species will occur for all phases
412  * listed in the kinetics object.
413  *
414  * return
415  * - If a match is found, the position in the species list is returned.
416  * - If no match is found, the value -1 is returned.
417  *
418  * @param nm Input string name of the species
419  */
420  size_t kineticsSpeciesIndex(const std::string& nm) const;
421 
422  /**
423  * This routine will look up a species number based on the input
424  * std::string nm. The lookup of species will occur in the specified
425  * phase of the object, or all phases if ph is "<any>".
426  *
427  * return
428  * - If a match is found, the position in the species list is returned.
429  * - If no match is found, the value npos (-1) is returned.
430  *
431  * @param nm Input string name of the species
432  * @param ph Input string name of the phase.
433  */
434  size_t kineticsSpeciesIndex(const std::string& nm,
435  const std::string& ph) const;
436 
437  /**
438  * This function looks up the std::string name of a species and
439  * returns a reference to the ThermoPhase object of the
440  * phase where the species resides.
441  * Will throw an error if the species std::string doesn't match.
442  *
443  * @param nm String containing the name of the species.
444  */
445  thermo_t& speciesPhase(std::string nm);
446 
447  /**
448  * This function takes as an argument the kineticsSpecies index
449  * (i.e., the list index in the list of species in the kinetics
450  * manager) and returns the species' owning ThermoPhase object.
451  *
452  * @param k Species index
453  */
454  thermo_t& speciesPhase(size_t k) {
455  return thermo(speciesPhaseIndex(k));
456  }
457 
458  /**
459  * This function takes as an argument the kineticsSpecies index
460  * (i.e., the list index in the list of species in the kinetics
461  * manager) and returns the index of the phase owning the
462  * species.
463  *
464  * @param k Species index
465  */
466  size_t speciesPhaseIndex(size_t k);
467 
468  //@}
469 
470 
471 
472  /**
473  * @name Reaction Rates Of Progress
474  */
475  //@{
476 
477  //! Return the forward rates of progress of the reactions
478  /*!
479  * Forward rates of progress. Return the forward rates of
480  * progress in array fwdROP, which must be dimensioned at
481  * least as large as the total number of reactions.
482  *
483  * @param fwdROP Output vector containing forward rates
484  * of progress of the reactions. Length: m_ii.
485  */
486  virtual void getFwdRatesOfProgress(doublereal* fwdROP) {
487  err("getFwdRatesOfProgress");
488  }
489 
490  //! Return the Reverse rates of progress of the reactions
491  /*!
492  * Return the reverse rates of
493  * progress in array revROP, which must be dimensioned at
494  * least as large as the total number of reactions.
495  *
496  * @param revROP Output vector containing reverse rates
497  * of progress of the reactions. Length: m_ii.
498  */
499  virtual void getRevRatesOfProgress(doublereal* revROP) {
500  err("getRevRatesOfProgress");
501  }
502 
503  /**
504  * Net rates of progress. Return the net (forward - reverse)
505  * rates of progress in array netROP, which must be
506  * dimensioned at least as large as the total number of
507  * reactions.
508  *
509  * @param netROP Output vector of the net ROP. Length: m_ii.
510  */
511  virtual void getNetRatesOfProgress(doublereal* netROP) {
512  err("getNetRatesOfProgress");
513  }
514 
515 
516 
517  //! Return a vector of Equilibrium constants.
518  /*!
519  * Return the equilibrium constants of
520  * the reactions in concentration units in array kc, which
521  * must be dimensioned at least as large as the total number
522  * of reactions.
523  *
524  * @param kc Output vector containing the equilibrium constants.
525  * Length: m_ii.
526  */
527  virtual void getEquilibriumConstants(doublereal* kc) {
528  err("getEquilibriumConstants");
529  }
530 
531  /**
532  * Change in species properties. Given an array of molar species
533  * property values \f$ z_k, k = 1, \dots, K \f$, return the
534  * array of reaction values
535  * \f[
536  * \Delta Z_i = \sum_k \nu_{k,i} z_k, i = 1, \dots, I.
537  * \f]
538  * For example, if this method is called with the array of
539  * standard-state molar Gibbs free energies for the species,
540  * then the values returned in array \c deltaProperty would be
541  * the standard-state Gibbs free energies of reaction for each
542  * reaction.
543  *
544  * @param property Input vector of property value. Length: m_kk.
545  * @param deltaProperty Output vector of deltaRxn. Length: m_ii.
546  */
547  virtual void getReactionDelta(const doublereal* property,
548  doublereal* deltaProperty) {
549  err("getReactionDelta");
550  }
551 
552  //! Return the vector of values for the reaction gibbs free energy change.
553  /*!
554  * These values depend upon the concentration of the solution.
555  *
556  * units = J kmol-1
557  *
558  * @param deltaG Output vector of deltaG's for reactions
559  * Length: m_ii.
560  */
561  virtual void getDeltaGibbs(doublereal* deltaG) {
562  err("getDeltaGibbs");
563  }
564 
565  //! Return the vector of values for the reaction electrochemical free energy change.
566  /*!
567  * These values depend upon the concentration of the solution and
568  * the voltage of the phases
569  *
570  * units = J kmol-1
571  *
572  * @param deltaM Output vector of deltaM's for reactions
573  * Length: m_ii.
574  */
575  virtual void getDeltaElectrochemPotentials(doublereal* deltaM) {
576  err("getDeltaElectrochemPotentials");
577  }
578 
579  /**
580  * Return the vector of values for the reactions change in
581  * enthalpy. These values depend upon the concentration of
582  * the solution.
583  *
584  * units = J kmol-1
585  *
586  * @param deltaH Output vector of deltaH's for reactions
587  * Length: m_ii.
588  */
589  virtual void getDeltaEnthalpy(doublereal* deltaH) {
590  err("getDeltaEnthalpy");
591  }
592 
593  /**
594  * Return the vector of values for the reactions change in
595  * entropy. These values depend upon the concentration of the
596  * solution.
597  *
598  * units = J kmol-1 Kelvin-1
599  *
600  * @param deltaS Output vector of deltaS's for reactions
601  * Length: m_ii.
602  */
603  virtual void getDeltaEntropy(doublereal* deltaS) {
604  err("getDeltaEntropy");
605  }
606 
607  /**
608  * Return the vector of values for the reaction standard state
609  * gibbs free energy change. These values don't depend upon
610  * the concentration of the solution.
611  *
612  * units = J kmol-1
613  *
614  * @param deltaG Output vector of ss deltaG's for reactions
615  * Length: m_ii.
616  */
617  virtual void getDeltaSSGibbs(doublereal* deltaG) {
618  err("getDeltaSSGibbs");
619  }
620 
621  /**
622  * Return the vector of values for the change in the standard
623  * state enthalpies of reaction. These values don't depend
624  * upon the concentration of the solution.
625  *
626  * units = J kmol-1
627  *
628  * @param deltaH Output vector of ss deltaH's for reactions
629  * Length: m_ii.
630  */
631  virtual void getDeltaSSEnthalpy(doublereal* deltaH) {
632  err("getDeltaSSEnthalpy");
633  }
634 
635  /**
636  * Return the vector of values for the change in the standard
637  * state entropies for each reaction. These values don't
638  * depend upon the concentration of the solution.
639  *
640  * units = J kmol-1 Kelvin-1
641  *
642  * @param deltaS Output vector of ss deltaS's for reactions
643  * Length: m_ii.
644  */
645  virtual void getDeltaSSEntropy(doublereal* deltaS) {
646  err("getDeltaSSEntropy");
647  }
648 
649 
650  //@}
651  /**
652  * @name Species Production Rates
653  */
654  //@{
655 
656  /**
657  * Species creation rates [kmol/m^3/s or kmol/m^2/s]. Return the
658  * species creation rates in array cdot, which must be
659  * dimensioned at least as large as the total number of
660  * species in all phases. @see nTotalSpecies.
661  *
662  * @param cdot Output vector of creation rates.
663  * Length: m_kk.
664  */
665  virtual void getCreationRates(doublereal* cdot) {
666  err("getCreationRates");
667  }
668 
669  /**
670  * Species destruction rates [kmol/m^3/s or kmol/m^2/s]. Return
671  * the species destruction rates in array ddot, which must be
672  * dimensioned at least as large as the total number of
673  * species. @see nTotalSpecies.
674  *
675  * @param ddot Output vector of destruction rates.
676  * Length: m_kk.
677  */
678  virtual void getDestructionRates(doublereal* ddot) {
679  err("getDestructionRates");
680  }
681 
682  /**
683  * Species net production rates [kmol/m^3/s or kmol/m^2/s]. Return
684  * the species net production rates (creation - destruction)
685  * in array wdot, which must be dimensioned at least as large
686  * as the total number of species. @see nTotalSpecies.
687  *
688  * @param wdot Output vector of net production rates.
689  * Length: m_kk.
690  */
691  virtual void getNetProductionRates(doublereal* wdot) {
692  err("getNetProductionRates");
693  }
694 
695  //@}
696 
697 
698  /**
699  * @name Reaction Mechanism Informational Query Routines
700  */
701  //@{
702 
703  /**
704  * Stoichiometric coefficient of species k as a reactant in
705  * reaction i.
706  *
707  * @param k kinetic species index
708  * @param i reaction index
709  */
710  virtual doublereal reactantStoichCoeff(size_t k, size_t i) const {
711  err("reactantStoichCoeff");
712  return -1.0;
713  }
714 
715  /**
716  * Stoichiometric coefficient of species k as a product in
717  * reaction i.
718  *
719  * @param k kinetic species index
720  * @param i reaction index
721  */
722  virtual doublereal productStoichCoeff(size_t k, size_t i) const {
723  err("productStoichCoeff");
724  return -1.0;
725  }
726 
727  //! Reactant order of species k in reaction i.
728  /*!
729  * This is the nominal order of the activity concentration in
730  * determining the forward rate of progress of the reaction
731  *
732  * @param k kinetic species index
733  * @param i reaction index
734  */
735  virtual doublereal reactantOrder(size_t k, size_t i) const {
736  err("reactantOrder");
737  return -1.0;
738  }
739 
740  //! product Order of species k in reaction i.
741  /*!
742  * This is the nominal order of the activity concentration of species k in
743  * determining the reverse rate of progress of the reaction i
744  *
745  * For irreversible reactions, this will all be zero.
746  *
747  * @param k kinetic species index
748  * @param i reaction index
749  */
750  virtual doublereal productOrder(int k, int i) const {
751  err("productOrder");
752  return -1.0;
753  }
754 
755  //! Get the vector of activity concentrations used in the kinetics object
756  /*!
757  * @param conc (output) Vector of activity concentrations. Length is
758  * equal to the number of species in the kinetics object
759  */
760  virtual void getActivityConcentrations(doublereal* const conc) {
761  err("getActivityConcentrations");
762  }
763 
764  /**
765  * Returns a read-only reference to the vector of reactant
766  * index numbers for reaction i.
767  *
768  * @param i reaction index
769  */
770  virtual const std::vector<size_t>& reactants(size_t i) const {
771  return m_reactants[i];
772  }
773 
774  /**
775  * Returns a read-only reference to the vector of product
776  * index numbers for reaction i.
777  *
778  * @param i reaction index
779  */
780  virtual const std::vector<size_t>& products(size_t i) const {
781  return m_products[i];
782  }
783 
784  /**
785  * Flag specifying the type of reaction. The legal values and
786  * their meaning are specific to the particular kinetics
787  * manager.
788  *
789  * @param i reaction index
790  */
791  virtual int reactionType(size_t i) const {
792  err("reactionType");
793  return -1;
794  }
795 
796  /**
797  * True if reaction i has been declared to be reversible. If
798  * isReversible(i) is false, then the reverse rate of progress
799  * for reaction i is always zero.
800  *
801  * @param i reaction index
802  */
803  virtual bool isReversible(size_t i) {
804  err("isReversible");
805  return false;
806  }
807 
808  /**
809  * Return a std::string representing the reaction.
810  *
811  * @param i reaction index
812  */
813  virtual std::string reactionString(size_t i) const {
814  err("reactionStd::String");
815  return "<null>";
816  }
817 
818  /**
819  * Return the forward rate constants
820  *
821  * length is the number of reactions. units depends
822  * on many issues. @todo DGG: recommend changing name to
823  * getFwdRateCoefficients.
824  *
825  * @param kfwd Output vector containing the forward reaction rate constants.
826  * Length: m_ii.
827  */
828  virtual void getFwdRateConstants(doublereal* kfwd) {
829  err("getFwdRateConstants");
830  }
831 
832  /**
833  * Return the reverse rate constants.
834  *
835  * length is the number of reactions. units depends
836  * on many issues. Note, this routine will return rate constants
837  * for irreversible reactions if the default for
838  * doIrreversible is overridden. @todo DGG: recommend changing name to
839  * getRevRateCoefficients.
840  *
841  * @param krev Output vector of reverse rate constants.
842  * @param doIrreversible boolean indicating whether irreversible reactions
843  * should be included.
844  */
845  virtual void getRevRateConstants(doublereal* krev,
846  bool doIrreversible = false) {
847  err("getFwdRateConstants");
848  }
849 
850 
851  /**
852  * Return the activation energies in Kelvin.
853  *
854  * length is the number of reactions
855  *
856  * @param E Ouptut vector of activation energies.
857  * Length: m_ii.
858  */
859  virtual void getActivationEnergies(doublereal* E) {
860  err("getActivationEnergies");
861  }
862 
863 
864  //@}
865  /**
866  * @name Reaction Mechanism Construction
867  */
868  //@{
869 
870  //! Add a phase to the kinetics manager object.
871  /*!
872  * This must be done before the function init() is called or
873  * before any reactions are input.
874  * The following fields are updated:
875  * m_start -> vector of integers, containing the
876  * starting position of the species for
877  * each phase in the kinetics mechanism.
878  * m_surfphase -> index of the surface phase.
879  * m_thermo -> vector of pointers to ThermoPhase phases
880  * that participate in the kinetics
881  * mechanism.
882  * m_phaseindex -> map containing the std::string id of each
883  * ThermoPhase phase as a key and the
884  * index of the phase within the kinetics
885  * manager object as the value.
886  *
887  * @param thermo Reference to the ThermoPhase to be added.
888  */
889  virtual void addPhase(thermo_t& thermo);
890 
891  /**
892  * Prepare the class for the addition of reactions. This
893  * method is called by function importKinetics after all
894  * phases have been added but before any reactions have
895  * been. The base class method does nothing, but derived
896  * classes may use this to perform any initialization
897  * (allocating arrays, etc.) that requires knowing the phases
898  * and species, but before any reactions are added.
899  */
900  virtual void init() {}
901 
902  /**
903  * Finish adding reactions and prepare for use. This method is
904  * called by function importKinetics after all reactions have
905  * been entered into the mechanism and before the mechanism is
906  * used to calculate reaction rates. The base class method
907  * does nothing, but derived classes may use this to perform
908  * any initialization (allocating arrays, etc.) that must be
909  * done after the reactions are entered.
910  */
911  virtual void finalize();
912 
913  /**
914  * Add a single reaction to the mechanism. This routine
915  * must be called after init() and before finalize().
916  *
917  * @param r Reference to the ReactionData object for the reaction
918  * to be added.
919  */
920  virtual void addReaction(ReactionData& r) {
921  err("addReaction");
922  }
923 
924  virtual const std::vector<grouplist_t>& reactantGroups(size_t i) {
925  //err("reactantGroups");
926  return m_dummygroups;
927  }
928 
929  virtual const std::vector<grouplist_t>& productGroups(size_t i) {
930  //err("productGroups");
931  return m_dummygroups;
932  }
933 
934 
935  //@}
936  /**
937  * @name Altering Reaction Rates
938  *
939  * These methods alter reaction rates. They are designed
940  * primarily for carrying out sensitivity analysis, but may be
941  * used for any purpose requiring dynamic alteration of rate
942  * constants. For each reaction, a real-valued multiplier may
943  * be defined that multiplies the reaction rate
944  * coefficient. The multiplier may be set to zero to
945  * completely remove a reaction from the mechanism.
946  */
947  //@{
948 
949  //! The current value of the multiplier for reaction i.
950  /*!
951  * @param i index of the reaction
952  */
953  doublereal multiplier(size_t i) const {
954  return m_perturb[i];
955  }
956 
957  //! Set the multiplier for reaction i to f.
958  /*!
959  * @param i index of the reaction
960  * @param f value of the multiplier.
961  */
962  void setMultiplier(size_t i, doublereal f) {
963  m_perturb[i] = f;
964  }
965 
966  //@}
967 
968  /**
969  * Increment the number of reactions in the mechanism by one.
970  * @todo Should be protected?
971  */
973  m_ii++;
974  m_perturb.push_back(1.0);
975  }
976 
977  /**
978  * Returns true if the kinetics manager has been properly
979  * initialized and finalized.
980  */
981  virtual bool ready() const {
982  return false;
983  }
984 
985 
986  /**
987  * Extract from array \c data the portion pertaining to phase \c phase.
988  *
989  * @param data data
990  * @param phase phase
991  * @param phase_data phase_data
992  */
993  void selectPhase(const doublereal* data, const thermo_t* phase,
994  doublereal* phase_data);
995 
996 protected:
997 
998 
999  //! Number of reactions in the mechanism
1000  size_t m_ii;
1001 
1002  //! The number of species in all of the phases
1003  //! that participate in this kinetics mechanism.
1004  size_t m_kk;
1005 
1006  /// Vector of perturbation factors for each reaction's rate of
1007  /// progress vector. It is initialized to one.
1008  ///
1010 
1011  /**
1012  * This is a vector of vectors containing the reactants for
1013  * each reaction. The outer vector is over the number of
1014  * reactions, m_ii. The inner vector is a list of species
1015  * indices. If the stoichiometric coefficient for a reactant
1016  * is greater than one, then the reactant is listed
1017  * contiguously in the vector a number of times equal to its
1018  * stoichiometric coefficient.
1019  * NOTE: These vectors will be wrong if there are real
1020  * stoichiometric coefficients in the expression.
1021  */
1022  std::vector<std::vector<size_t> > m_reactants;
1023 
1024  /**
1025  * This is a vector of vectors containing the products for
1026  * each reaction. The outer vector is over the number of
1027  * reactions, m_ii. The inner vector is a list of species
1028  * indices. If the stoichiometric coefficient for a product is
1029  * greater than one, then the reactant is listed contiguously
1030  * in the vector a number of times equal to its stoichiometric
1031  * coefficient.
1032  * NOTE: These vectors will be wrong if there are real
1033  * stoichiometric coefficients in the expression.
1034  */
1035  std::vector<std::vector<size_t> > m_products;
1036 
1037  //! m_thermo is a vector of pointers to ThermoPhase objects that are involved with this kinetics operator
1038  /*!
1039  * For homogeneous kinetics applications, this vector
1040  * will only have one entry. For interfacial reactions, this
1041  * vector will consist of multiple entries; some of them will
1042  * be surface phases, and the other ones will be bulk phases.
1043  * The order that the objects are listed determines the order
1044  * in which the species comprising each phase are listed in
1045  * the source term vector, originating from the reaction
1046  * mechanism.
1047  *
1048  * Note that this kinetics object doesn't own these ThermoPhase
1049  * objects and is not responsible for creating or deleting
1050  * them.
1051  */
1052  std::vector<thermo_t*> m_thermo;
1053 
1054  /**
1055  * m_start is a vector of integers specifying the beginning position
1056  * for the species vector for the n'th phase in the kinetics
1057  * class.
1058  */
1059  std::vector<size_t> m_start;
1060 
1061  /**
1062  * Mapping of the phase id, i.e., the id attribute in the xml
1063  * phase element to the position of the phase within the
1064  * kinetics object. Positions start with the value of 1. The
1065  * member function, phaseIndex() decrements by one before
1066  * returning the index value, so that missing phases return
1067  * -1.
1068  */
1069  std::map<std::string, size_t> m_phaseindex;
1070 
1071  //! Index in the list of phases of the one surface phase.
1072  /*!
1073  *
1074  */
1075  size_t m_surfphase;
1076 
1077  //! Phase Index where reactions are assumed to be taking place
1078  /*!
1079  * We calculate this by assuming that the phase with the lowest dimensionality is the phase where reactions
1080  * are taking place
1081  * @deprecated
1082  */
1083  size_t m_rxnphase;
1084 
1085  //! number of spatial dimensions of lowest-dimensional phase.
1086  size_t m_mindim;
1087 
1088 private:
1089 
1090  //! Vector of group lists
1091  std::vector<grouplist_t> m_dummygroups;
1092 
1093 
1094  //! Private function of the class Kinetics, indicating that a function
1095  //! inherited from the base class hasn't had a definition assigned to it
1096  /*!
1097  * @param m String message
1098  */
1099  void err(std::string m) const;
1100 
1101 };
1102 
1103 }
1104 
1105 #endif