Cantera  2.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
StoichManager.h
Go to the documentation of this file.
1 /**
2  * @file StoichManager.h
3  */
4 
5 // Copyright 2001 California Institute of Technology
6 
7 #ifndef CT_STOICH_MGR_H
8 #define CT_STOICH_MGR_H
9 
12 
13 namespace Cantera
14 {
15 
16 /**
17  * @defgroup Stoichiometry Stoichiometry
18  *
19  * Note: these classes are designed for internal use in class
20  * ReactionStoichManager.
21  *
22  * The classes defined here implement simple operations that are used by class
23  * ReactionStoichManager to compute things like rates of progress, species
24  * production rates, etc. In general, a reaction mechanism may involve many
25  * species and many reactions, but any given reaction typically only involves
26  * a few species as reactants, and a few as products. Therefore, the matrix of
27  * stoichiometric coefficients is very sparse. Not only is it sparse, but the
28  * non-zero matrix elements often have the value 1, and in many cases no more
29  * than three coefficients are non-zero for the reactants and/or the products.
30  *
31  * For the present purposes, we will consider each direction of a reversible
32  * reaction to be a separate reaction. We often need to compute quantities
33  * that can formally be written as a matrix product of a stoichiometric
34  * coefficient matrix and a vector of reaction rates. For example, the species
35  * creation rates are given by
36  *
37  * \f[
38  * \dot C_k = \sum_k \nu^{(p)}_{k,i} R_i
39  * \f]
40  *
41  * where \f$ \nu^{(p)_{k,i}} \f$ is the product-side stoichiometric
42  * coefficient of species \a k in reaction \a i. This could be done be
43  * straightforward matrix multiplication, but would be inefficient, since most
44  * of the matrix elements of \f$ \nu^{(p)}_{k,i} \f$ are zero. We could do
45  * better by using sparse-matrix algorithms to compute this product.
46  *
47  * If the reactions are general ones, with non-integral stoichiometric
48  * coefficients, this is about as good as we can do. But we are particularly
49  * concerned here with the performance for very large reaction mechanisms,
50  * which are usually composed of elementary reactions, which have integral
51  * stoichiometric coefficients. Furthermore, very few elementary reactions
52  * involve more than 3 product or reactant molecules.
53  *
54  * But we can do even better if we take account of the special structure of
55  * this matrix for elementary reactions involving three or fewer product
56  * molecules (or reactant molecules).
57  *
58  * To take advantage of this structure, reactions are divided into four
59  * groups. These classes are designed to take advantage of this sparse
60  * structure when computing quantities that can be written as matrix
61  * multiplies.
62  *
63  * They are designed to explicitly unroll loops over species or reactions for
64  * Operations on reactions that require knowing the reaction stoichiometry.
65  *
66  * This module consists of class StoichManager, and classes C1, C2, and C3.
67  * Classes C1, C2, and C3 handle operations involving one, two, or three
68  * species, respectively, in a reaction. Instances are instantiated with a
69  * reaction number, and n species numbers (n = 1 for C1, etc.). All three
70  * classes have the same interface.
71  *
72  * These classes are designed for use by StoichManager, and the operations
73  * implemented are those needed to efficiently compute quantities such as
74  * rates of progress, species production rates, reaction thermochemistry, etc.
75  * The compiler will inline these methods into the body of the corresponding
76  * StoichManager method, and so there is no performance penalty (unless
77  * inlining is turned off).
78  *
79  * To describe the methods, consider class C3 and suppose an instance is
80  * created with reaction number irxn and species numbers k0, k1, and k2.
81  *
82  * - multiply(in, out) : out[irxn] is multiplied by
83  * in[k0] * in[k1] * in[k2]
84  *
85  * - power(in, out) : out[irxn] is multiplied by
86  * (in[k0]^order0) * (in[k1]^order1) * (in[k2]^order2)
87  *
88  * - incrementReaction(in, out) : out[irxn] is incremented by
89  * in[k0] + in[k1] + in[k2]
90  *
91  * - decrementReaction(in, out) : out[irxn] is decremented by
92  * in[k0] + in[k1] + in[k2]
93  *
94  * - incrementSpecies(in, out) : out[k0], out[k1], and out[k2]
95  * are all incremented by in[irxn]
96  *
97  * - decrementSpecies(in, out) : out[k0], out[k1], and out[k2]
98  * are all decremented by in[irxn]
99  *
100  * The function multiply() is usually used when evaluating the forward and
101  * reverse rates of progress of reactions. The rate constants are usually
102  * loaded into out[]. Then multiply() is called to add in the dependence of
103  * the species concentrations to yield a forward and reverse rop.
104  *
105  * The function incrementSpecies() and its cousin decrementSpecies() is used
106  * to translate from rates of progress to species production rates. The vector
107  * in[] is preloaded with the rates of progress of all reactions. Then
108  * incrementSpecies() is called to increment the species production vector,
109  * out[], with the rates of progress.
110  *
111  * The functions incrementReaction() and decrementReaction() are used to find
112  * the standard state equilibrium constant for a reaction. Here, output[] is a
113  * vector of length number of reactions, usually the standard Gibbs free
114  * energies of reaction, while input, usually the standard state Gibbs free
115  * energies of species, is a vector of length number of species.
116  *
117  * Note the stoichiometric coefficient for a species in a reaction is handled
118  * by always assuming it is equal to one and then treating reactants and
119  * products for a reaction separately. Bimolecular reactions involving the
120  * identical species are treated as involving separate species.
121  *
122  * @internal This class should be upgraded to include cases where
123  * real stoichiometric coefficients are used. Shouldn't be that
124  * hard to do, and they occur in engineering simulations with some
125  * regularity.
126  *
127  */
128 
129 static doublereal ppow(doublereal x, doublereal order)
130 {
131  if (x > 0.0) {
132  return std::pow(x, order);
133  } else {
134  return 0.0;
135  }
136 }
137 
138 //! @deprecated To be removed after Cantera 2.2
139 inline static std::string fmt(const std::string& r, size_t n)
140 {
141  return r + "[" + int2str(n) + "]";
142 }
143 
144 /**
145  * Handles one species in a reaction.
146  * See @ref Stoichiometry
147  * @ingroup Stoichiometry
148  * @internal
149  */
150 class C1
151 {
152 public:
153  C1(size_t rxn = 0, size_t ic0 = 0) :
154  m_rxn(rxn),
155  m_ic0(ic0) {
156  }
157 
158  size_t data(std::vector<size_t>& ic) {
159  ic.resize(1);
160  ic[0] = m_ic0;
161  return m_rxn;
162  }
163 
164  void incrementSpecies(const doublereal* R, doublereal* S) const {
165  S[m_ic0] += R[m_rxn];
166  }
167 
168  void decrementSpecies(const doublereal* R, doublereal* S) const {
169  S[m_ic0] -= R[m_rxn];
170  }
171 
172  void multiply(const doublereal* S, doublereal* R) const {
173  R[m_rxn] *= S[m_ic0];
174  }
175 
176  void incrementReaction(const doublereal* S, doublereal* R) const {
177  R[m_rxn] += S[m_ic0];
178  }
179 
180  void decrementReaction(const doublereal* S, doublereal* R) const {
181  R[m_rxn] -= S[m_ic0];
182  }
183 
184  size_t rxnNumber() const {
185  return m_rxn;
186  }
187  size_t speciesIndex(size_t n) const {
188  return m_ic0;
189  }
190  size_t nSpecies() {
191  return 1;
192  }
193 
194  //! @deprecated To be removed after Cantera 2.2
195  void writeMultiply(const std::string& r, std::map<size_t, std::string>& out) {
196  out[m_rxn] = fmt(r, m_ic0);
197  }
198 
199  //! @deprecated To be removed after Cantera 2.2
200  void writeIncrementReaction(const std::string& r, std::map<size_t, std::string>& out) {
201  out[m_rxn] += " + "+fmt(r, m_ic0);
202  }
203 
204  //! @deprecated To be removed after Cantera 2.2
205  void writeDecrementReaction(const std::string& r, std::map<size_t, std::string>& out) {
206  out[m_rxn] += " - "+fmt(r, m_ic0);
207  }
208 
209  //! @deprecated To be removed after Cantera 2.2
210  void writeIncrementSpecies(const std::string& r, std::map<size_t, std::string>& out) {
211  out[m_ic0] += " + "+fmt(r, m_rxn);
212  }
213 
214  //! @deprecated To be removed after Cantera 2.2
215  void writeDecrementSpecies(const std::string& r, std::map<size_t, std::string>& out) {
216  out[m_ic0] += " - "+fmt(r, m_rxn);
217  }
218 
219 private:
220  //! Reaction number
221  size_t m_rxn;
222  //! Species number
223  size_t m_ic0;
224 };
225 
226 /**
227  * Handles two species in a single reaction.
228  * See @ref Stoichiometry
229  * @ingroup Stoichiometry
230  */
231 class C2
232 {
233 public:
234  C2(size_t rxn = 0, size_t ic0 = 0, size_t ic1 = 0)
235  : m_rxn(rxn), m_ic0(ic0), m_ic1(ic1) {}
236 
237  size_t data(std::vector<size_t>& ic) {
238  ic.resize(2);
239  ic[0] = m_ic0;
240  ic[1] = m_ic1;
241  return m_rxn;
242  }
243 
244  void incrementSpecies(const doublereal* R, doublereal* S) const {
245  S[m_ic0] += R[m_rxn];
246  S[m_ic1] += R[m_rxn];
247  }
248 
249  void decrementSpecies(const doublereal* R, doublereal* S) const {
250  S[m_ic0] -= R[m_rxn];
251  S[m_ic1] -= R[m_rxn];
252  }
253 
254  void multiply(const doublereal* S, doublereal* R) const {
255  R[m_rxn] *= S[m_ic0] * S[m_ic1];
256  }
257 
258  void incrementReaction(const doublereal* S, doublereal* R) const {
259  R[m_rxn] += S[m_ic0] + S[m_ic1];
260  }
261 
262  void decrementReaction(const doublereal* S, doublereal* R) const {
263  R[m_rxn] -= (S[m_ic0] + S[m_ic1]);
264  }
265 
266  size_t rxnNumber() const {
267  return m_rxn;
268  }
269  size_t speciesIndex(size_t n) const {
270  return (n == 0 ? m_ic0 : m_ic1);
271  }
272  size_t nSpecies() {
273  return 2;
274  }
275 
276  //! @deprecated To be removed after Cantera 2.2
277  void writeMultiply(const std::string& r, std::map<size_t, std::string>& out) {
278  out[m_rxn] = fmt(r, m_ic0) + " * " + fmt(r, m_ic1);
279  }
280 
281  //! @deprecated To be removed after Cantera 2.2
282  void writeIncrementReaction(const std::string& r, std::map<size_t, std::string>& out) {
283  out[m_rxn] += " + "+fmt(r, m_ic0)+" + "+fmt(r, m_ic1);
284  }
285 
286  //! @deprecated To be removed after Cantera 2.2
287  void writeDecrementReaction(const std::string& r, std::map<size_t, std::string>& out) {
288  out[m_rxn] += " - "+fmt(r, m_ic0)+" - "+fmt(r, m_ic1);
289  }
290 
291  //! @deprecated To be removed after Cantera 2.2
292  void writeIncrementSpecies(const std::string& r, std::map<size_t, std::string>& out) {
293  std::string s = " + "+fmt(r, m_rxn);
294  out[m_ic0] += s;
295  out[m_ic1] += s;
296  }
297 
298  //! @deprecated To be removed after Cantera 2.2
299  void writeDecrementSpecies(const std::string& r, std::map<size_t, std::string>& out) {
300  std::string s = " - "+fmt(r, m_rxn);
301  out[m_ic0] += s;
302  out[m_ic1] += s;
303  }
304 
305 private:
306  //! Reaction index -> index into the ROP vector
307  size_t m_rxn;
308 
309  //! Species index -> index into the species vector for the two species.
310  size_t m_ic0, m_ic1;
311 };
312 
313 /**
314  * Handles three species in a reaction.
315  * See @ref Stoichiometry
316  * @ingroup Stoichiometry
317  */
318 class C3
319 {
320 public:
321  C3(size_t rxn = 0, size_t ic0 = 0, size_t ic1 = 0, size_t ic2 = 0)
322  : m_rxn(rxn), m_ic0(ic0), m_ic1(ic1), m_ic2(ic2) {}
323 
324  size_t data(std::vector<size_t>& ic) {
325  ic.resize(3);
326  ic[0] = m_ic0;
327  ic[1] = m_ic1;
328  ic[2] = m_ic2;
329  return m_rxn;
330  }
331 
332  void incrementSpecies(const doublereal* R, doublereal* S) const {
333  S[m_ic0] += R[m_rxn];
334  S[m_ic1] += R[m_rxn];
335  S[m_ic2] += R[m_rxn];
336  }
337 
338  void decrementSpecies(const doublereal* R, doublereal* S) const {
339  S[m_ic0] -= R[m_rxn];
340  S[m_ic1] -= R[m_rxn];
341  S[m_ic2] -= R[m_rxn];
342  }
343 
344  void multiply(const doublereal* S, doublereal* R) const {
345  R[m_rxn] *= S[m_ic0] * S[m_ic1] * S[m_ic2];
346  }
347 
348  void incrementReaction(const doublereal* S, doublereal* R) const {
349  R[m_rxn] += S[m_ic0] + S[m_ic1] + S[m_ic2];
350  }
351 
352  void decrementReaction(const doublereal* S, doublereal* R) const {
353  R[m_rxn] -= (S[m_ic0] + S[m_ic1] + S[m_ic2]);
354  }
355 
356  size_t rxnNumber() const {
357  return m_rxn;
358  }
359  size_t speciesIndex(size_t n) const {
360  return (n == 0 ? m_ic0 : (n == 1 ? m_ic1 : m_ic2));
361  }
362  size_t nSpecies() {
363  return 3;
364  }
365 
366  //! @deprecated To be removed after Cantera 2.2
367  void writeMultiply(const std::string& r, std::map<size_t, std::string>& out) {
368  out[m_rxn] = fmt(r, m_ic0) + " * " + fmt(r, m_ic1) + " * " + fmt(r, m_ic2);
369  }
370 
371  //! @deprecated To be removed after Cantera 2.2
372  void writeIncrementReaction(const std::string& r, std::map<size_t, std::string>& out) {
373  out[m_rxn] += " + "+fmt(r, m_ic0)+" + "+fmt(r, m_ic1)+" + "+fmt(r, m_ic2);
374  }
375 
376  //! @deprecated To be removed after Cantera 2.2
377  void writeDecrementReaction(const std::string& r, std::map<size_t, std::string>& out) {
378  out[m_rxn] += " - "+fmt(r, m_ic0)+" - "+fmt(r, m_ic1)+" - "+fmt(r, m_ic2);
379  }
380 
381  //! @deprecated To be removed after Cantera 2.2
382  void writeIncrementSpecies(const std::string& r, std::map<size_t, std::string>& out) {
383  std::string s = " + "+fmt(r, m_rxn);
384  out[m_ic0] += s;
385  out[m_ic1] += s;
386  out[m_ic2] += s;
387  }
388 
389  //! @deprecated To be removed after Cantera 2.2
390  void writeDecrementSpecies(const std::string& r, std::map<size_t, std::string>& out) {
391  std::string s = " - "+fmt(r, m_rxn);
392  out[m_ic0] += s;
393  out[m_ic1] += s;
394  out[m_ic2] += s;
395  }
396 private:
397  size_t m_rxn;
398  size_t m_ic0;
399  size_t m_ic1;
400  size_t m_ic2;
401 };
402 
403 /**
404  * Handles any number of species in a reaction, including fractional
405  * stoichiometric coefficients, and arbitrary reaction orders.
406  * See @ref Stoichiometry
407  * @ingroup Stoichiometry
408  */
409 class C_AnyN
410 {
411 public:
412  C_AnyN() :
413  m_n(0),
414  m_rxn(npos) {
415  }
416 
417  C_AnyN(size_t rxn, const std::vector<size_t>& ic, const vector_fp& order_, const vector_fp& stoich_) :
418  m_n(0),
419  m_rxn(rxn) {
420  m_n = ic.size();
421  m_ic.resize(m_n);
422  m_order.resize(m_n);
423  m_stoich.resize(m_n);
424  for (size_t n = 0; n < m_n; n++) {
425  m_ic[n] = ic[n];
426  m_order[n] = order_[n];
427  m_stoich[n] = stoich_[n];
428  }
429  }
430 
431  size_t data(std::vector<size_t>& ic) {
432  ic.resize(m_n);
433  for (size_t n = 0; n < m_n; n++) {
434  ic[n] = m_ic[n];
435  }
436  return m_rxn;
437  }
438 
439  doublereal order(size_t n) const {
440  return m_order[n];
441  }
442  doublereal stoich(size_t n) const {
443  return m_stoich[n];
444  }
445  size_t speciesIndex(size_t n) const {
446  return m_ic[n];
447  }
448 
449  void multiply(const doublereal* input, doublereal* output) const {
450  doublereal oo;
451  for (size_t n = 0; n < m_n; n++) {
452  oo = m_order[n];
453  if (oo != 0.0) {
454  output[m_rxn] *= ppow(input[m_ic[n]], oo);
455  }
456  }
457  }
458 
459  void incrementSpecies(const doublereal* input,
460  doublereal* output) const {
461  doublereal x = input[m_rxn];
462  for (size_t n = 0; n < m_n; n++) {
463  output[m_ic[n]] += m_stoich[n]*x;
464  }
465  }
466 
467  void decrementSpecies(const doublereal* input,
468  doublereal* output) const {
469  doublereal x = input[m_rxn];
470  for (size_t n = 0; n < m_n; n++) {
471  output[m_ic[n]] -= m_stoich[n]*x;
472  }
473  }
474 
475  void incrementReaction(const doublereal* input,
476  doublereal* output) const {
477  for (size_t n = 0; n < m_n; n++) output[m_rxn]
478  += m_stoich[n]*input[m_ic[n]];
479  }
480 
481  void decrementReaction(const doublereal* input,
482  doublereal* output) const {
483  for (size_t n = 0; n < m_n; n++) output[m_rxn]
484  -= m_stoich[n]*input[m_ic[n]];
485  }
486 
487  //! @deprecated To be removed after Cantera 2.2
488  void writeMultiply(const std::string& r, std::map<size_t, std::string>& out) {
489  out[m_rxn] = "";
490  for (size_t n = 0; n < m_n; n++) {
491  if (m_order[n] == 1.0) {
492  out[m_rxn] += fmt(r, m_ic[n]);
493  } else {
494  out[m_rxn] += "pow("+fmt(r, m_ic[n])+","+fp2str(m_order[n])+")";
495  }
496  if (n < m_n-1) {
497  out[m_rxn] += " * ";
498  }
499  }
500  }
501 
502  //! @deprecated To be removed after Cantera 2.2
503  void writeIncrementReaction(const std::string& r, std::map<size_t, std::string>& out) {
504  for (size_t n = 0; n < m_n; n++) {
505  out[m_rxn] += " + "+fp2str(m_stoich[n]) + "*" + fmt(r, m_ic[n]);
506  }
507  }
508 
509  //! @deprecated To be removed after Cantera 2.2
510  void writeDecrementReaction(const std::string& r, std::map<size_t, std::string>& out) {
511  for (size_t n = 0; n < m_n; n++) {
512  out[m_rxn] += " - "+fp2str(m_stoich[n]) + "*" + fmt(r, m_ic[n]);
513  }
514  }
515 
516  //! @deprecated To be removed after Cantera 2.2
517  void writeIncrementSpecies(const std::string& r, std::map<size_t, std::string>& out) {
518  std::string s = fmt(r, m_rxn);
519  for (size_t n = 0; n < m_n; n++) {
520  out[m_ic[n]] += " + "+fp2str(m_stoich[n]) + "*" + s;
521  }
522  }
523 
524  //! @deprecated To be removed after Cantera 2.2
525  void writeDecrementSpecies(const std::string& r, std::map<size_t, std::string>& out) {
526  std::string s = fmt(r, m_rxn);
527  for (size_t n = 0; n < m_n; n++) {
528  out[m_ic[n]] += " - "+fp2str(m_stoich[n]) + "*" + s;
529  }
530  }
531 
532 private:
533 
534  //! Length of the m_ic vector
535  /*!
536  * This is the number of species which participate in the reaction order
537  * and stoichiometric coefficient vectors for the reactant or product description
538  * of the reaction.
539  */
540  size_t m_n;
541 
542  //! ID of the reaction corresponding to this stoichiometric manager
543  /*!
544  * This is used within the interface to select the array position to read and write to
545  * Normally this is associated with the reaction number in an array of quantities indexed
546  * by the reaction number, e.g., ROP[irxn].
547  */
548  size_t m_rxn;
549 
550  //! Vector of species which are involved with this stoichiometric manager calculations
551  /*!
552  * This is an integer list of species which participate in either the
553  * reaction order matrix or the stoichiometric order matrix for this reaction, m_rxn.
554  */
555  std::vector<size_t> m_ic;
556 
557  //! Reaction orders for the reaction
558  /*!
559  * This is either for the reactants or products.
560  * Length = m_n
561  * Species number, m_ic[n], has a reaction order of m_order[n].
562  */
564 
565  //! Stoichiometric coefficients for the reaction, reactant or product side.
566  /*!
567  * This is either for the reactants or products.
568  * Length = m_n
569  * Species number m_ic[m], has a stoichiometric coefficient of m_stoich[n].
570  */
572 };
573 
574 template<class InputIter, class Vec1, class Vec2>
575 inline static void _multiply(InputIter begin, InputIter end,
576  const Vec1& input, Vec2& output)
577 {
578  for (; begin != end; ++begin) {
579  begin->multiply(input, output);
580  }
581 }
582 
583 template<class InputIter, class Vec1, class Vec2>
584 inline static void _incrementSpecies(InputIter begin,
585  InputIter end, const Vec1& input, Vec2& output)
586 {
587  for (; begin != end; ++begin) {
588  begin->incrementSpecies(input, output);
589  }
590 }
591 
592 template<class InputIter, class Vec1, class Vec2>
593 inline static void _decrementSpecies(InputIter begin,
594  InputIter end, const Vec1& input, Vec2& output)
595 {
596  for (; begin != end; ++begin) {
597  begin->decrementSpecies(input, output);
598  }
599 }
600 
601 template<class InputIter, class Vec1, class Vec2>
602 inline static void _incrementReactions(InputIter begin,
603  InputIter end, const Vec1& input, Vec2& output)
604 {
605  for (; begin != end; ++begin) {
606  begin->incrementReaction(input, output);
607  }
608 }
609 
610 template<class InputIter, class Vec1, class Vec2>
611 inline static void _decrementReactions(InputIter begin,
612  InputIter end, const Vec1& input, Vec2& output)
613 {
614  for (; begin != end; ++begin) {
615  begin->decrementReaction(input, output);
616  }
617 }
618 
619 //! @deprecated To be removed after Cantera 2.2
620 template<class InputIter>
621 inline static void _writeIncrementSpecies(InputIter begin, InputIter end,
622  const std::string& r, std::map<size_t, std::string>& out)
623 {
624  for (; begin != end; ++begin) {
625  begin->writeIncrementSpecies(r, out);
626  }
627 }
628 
629 //! @deprecated To be removed after Cantera 2.2
630 template<class InputIter>
631 inline static void _writeDecrementSpecies(InputIter begin, InputIter end,
632  const std::string& r, std::map<size_t, std::string>& out)
633 {
634  for (; begin != end; ++begin) {
635  begin->writeDecrementSpecies(r, out);
636  }
637 }
638 
639 //! @deprecated To be removed after Cantera 2.2
640 template<class InputIter>
641 inline static void _writeIncrementReaction(InputIter begin, InputIter end,
642  const std::string& r, std::map<size_t, std::string>& out)
643 {
644  for (; begin != end; ++begin) {
645  begin->writeIncrementReaction(r, out);
646  }
647 }
648 
649 //! @deprecated To be removed after Cantera 2.2
650 template<class InputIter>
651 inline static void _writeDecrementReaction(InputIter begin, InputIter end,
652  const std::string& r, std::map<size_t, std::string>& out)
653 {
654  for (; begin != end; ++begin) {
655  begin->writeDecrementReaction(r, out);
656  }
657 }
658 
659 //! @deprecated To be removed after Cantera 2.2
660 template<class InputIter>
661 inline static void _writeMultiply(InputIter begin, InputIter end,
662  const std::string& r, std::map<size_t, std::string>& out)
663 {
664  for (; begin != end; ++begin) {
665  begin->writeMultiply(r, out);
666  }
667 }
668 
669 /*
670  * This class handles operations involving the stoichiometric
671  * coefficients on one side of a reaction (reactant or product) for
672  * a set of reactions comprising a reaction mechanism. This class is
673  * used by class ReactionStoichMgr, which contains three instances
674  * of this class (one to handle operations on the reactions, one for
675  * the products of reversible reactions, and one for the products of
676  * irreversible reactions).
677  *
678  * This class is designed for use with elementary reactions, or at
679  * least ones with integral stoichiometric coefficients. Let \f$ M(i) \f$
680  * be the number of molecules on the product or reactant side of
681  * reaction number i.
682  * \f[
683  * r_i = \sum_m^{M_i} s_{k_{m,i}}
684  * \f]
685  * To understand the operations performed by this class, let
686  * \f$ N_{k,i}\f$ denote the stoichiometric coefficient of species k on
687  * one side (reactant or product) in reaction i. Then \b N is a sparse
688  * K by I matrix of stoichiometric coefficients.
689  *
690  * The following matrix operations may be carried out with a vector
691  * S of length K, and a vector R of length I:
692  *
693  * - \f$ S = S + N R\f$ (incrementSpecies)
694  * - \f$ S = S - N R\f$ (decrementSpecies)
695  * - \f$ R = R + N^T S \f$ (incrementReaction)
696  * - \f$ R = R - N^T S \f$ (decrementReaction)
697  *
698  * The actual implementation, however, does not compute these
699  * quantities by matrix multiplication. A faster algorithm is used
700  * that makes use of the fact that the \b integer-valued N matrix is
701  * very sparse, and the non-zero terms are small positive integers.
702  * \f[
703  * S_k = R_{i1} + \dots + R_{iM}
704  * \f]
705  * where M is the number of molecules, and $\f i(m) \f$ is the
706  * See @ref Stoichiometry
707  * @ingroup Stoichiometry
708  */
709 class StoichManagerN
710 {
711 public:
712  /**
713  * Constructor for the StoichManagerN class.
714  *
715  * @internal Consider adding defaulted entries here that supply
716  * the total number of reactions in the mechanism and the total
717  * number of species in the species list. Then, we could use those
718  * numbers to provide error checks during the construction of the
719  * object. Those numbers would also provide some clarity to the
720  * purpose and utility of this class.
721  *
722  * DGG - the problem is that the number of reactions and species
723  * are not known initially.
724  */
725  StoichManagerN() {
726  }
727 
728  /**
729  * Add a single reaction to the list of reactions that this
730  * stoichiometric manager object handles.
731  *
732  * This function is the same as the add() function below. However,
733  * the order of each species in the power list expression is
734  * set to one automatically.
735  */
736  void add(size_t rxn, const std::vector<size_t>& k) {
737  vector_fp order(k.size(), 1.0);
738  vector_fp stoich(k.size(), 1.0);
739  add(rxn, k, order, stoich);
740  }
741 
742  void add(size_t rxn, const std::vector<size_t>& k, const vector_fp& order) {
743  vector_fp stoich(k.size(), 1.0);
744  add(rxn, k, order, stoich);
745  }
746 
747  //! Add a single reaction to the list of reactions that this
748  //! stoichiometric manager object handles.
749  /*!
750  * @param rxn Reaction index of the current reaction. This is used
751  * as an index into vectors which have length n_total_rxn.
752  * @param k This is a vector of integer values specifying the
753  * species indices. The length of this vector species
754  * the number of different species in the description.
755  * The value of the entries are the species indices.
756  * These are used as indexes into vectors which have
757  * length n_total_species.
758  * @param order This is a vector of the same length as vector k.
759  * The order is used for the routine power(), which produces
760  * a power law expression involving the species vector.
761  * @param stoich This is used to handle fractional stoichiometric coefficients
762  * on the product side of irreversible reactions.
763  */
764  void add(size_t rxn, const std::vector<size_t>& k, const vector_fp& order,
765  const vector_fp& stoich) {
766  if (order.size() != k.size()) {
767  throw CanteraError("StoichManagerN::add()", "size of order and species arrays differ");
768  }
769  if (stoich.size() != k.size()) {
770  throw CanteraError("StoichManagerN::add()", "size of stoich and species arrays differ");
771  }
772  bool frac = false;
773  for (size_t n = 0; n < stoich.size(); n++) {
774  if (fmod(stoich[n], 1.0) || stoich[n] != order[n]) {
775  frac = true;
776  break;
777  }
778  }
779  if (frac || k.size() > 3) {
780  m_cn_list.push_back(C_AnyN(rxn, k, order, stoich));
781  } else {
782  // Try to express the reaction with unity stoichiometric
783  // coefficients (by repeating species when necessary) so that the
784  // simpler 'multiply' function can be used to compute the rate
785  // instead of 'power'.
786  std::vector<size_t> kRep;
787  for (size_t n = 0; n < k.size(); n++) {
788  for (size_t i = 0; i < stoich[n]; i++)
789  kRep.push_back(k[n]);
790  }
791 
792  switch (kRep.size()) {
793  case 1:
794  m_c1_list.push_back(C1(rxn, kRep[0]));
795  break;
796  case 2:
797  m_c2_list.push_back(C2(rxn, kRep[0], kRep[1]));
798  break;
799  case 3:
800  m_c3_list.push_back(C3(rxn, kRep[0], kRep[1], kRep[2]));
801  break;
802  default:
803  m_cn_list.push_back(C_AnyN(rxn, k, order, stoich));
804  }
805  }
806  }
807 
808  void multiply(const doublereal* input, doublereal* output) const {
809  _multiply(m_c1_list.begin(), m_c1_list.end(), input, output);
810  _multiply(m_c2_list.begin(), m_c2_list.end(), input, output);
811  _multiply(m_c3_list.begin(), m_c3_list.end(), input, output);
812  _multiply(m_cn_list.begin(), m_cn_list.end(), input, output);
813  }
814 
815  void incrementSpecies(const doublereal* input, doublereal* output) const {
816  _incrementSpecies(m_c1_list.begin(), m_c1_list.end(), input, output);
817  _incrementSpecies(m_c2_list.begin(), m_c2_list.end(), input, output);
818  _incrementSpecies(m_c3_list.begin(), m_c3_list.end(), input, output);
819  _incrementSpecies(m_cn_list.begin(), m_cn_list.end(), input, output);
820  }
821 
822  void decrementSpecies(const doublereal* input, doublereal* output) const {
823  _decrementSpecies(m_c1_list.begin(), m_c1_list.end(), input, output);
824  _decrementSpecies(m_c2_list.begin(), m_c2_list.end(), input, output);
825  _decrementSpecies(m_c3_list.begin(), m_c3_list.end(), input, output);
826  _decrementSpecies(m_cn_list.begin(), m_cn_list.end(), input, output);
827  }
828 
829  void incrementReactions(const doublereal* input, doublereal* output) const {
830  _incrementReactions(m_c1_list.begin(), m_c1_list.end(), input, output);
831  _incrementReactions(m_c2_list.begin(), m_c2_list.end(), input, output);
832  _incrementReactions(m_c3_list.begin(), m_c3_list.end(), input, output);
833  _incrementReactions(m_cn_list.begin(), m_cn_list.end(), input, output);
834  }
835 
836  void decrementReactions(const doublereal* input, doublereal* output) const {
837  _decrementReactions(m_c1_list.begin(), m_c1_list.end(), input, output);
838  _decrementReactions(m_c2_list.begin(), m_c2_list.end(), input, output);
839  _decrementReactions(m_c3_list.begin(), m_c3_list.end(), input, output);
840  _decrementReactions(m_cn_list.begin(), m_cn_list.end(), input, output);
841  }
842 
843  //! @deprecated To be removed after Cantera 2.2
844  void writeIncrementSpecies(const std::string& r, std::map<size_t, std::string>& out) {
845  _writeIncrementSpecies(m_c1_list.begin(), m_c1_list.end(), r, out);
846  _writeIncrementSpecies(m_c2_list.begin(), m_c2_list.end(), r, out);
847  _writeIncrementSpecies(m_c3_list.begin(), m_c3_list.end(), r, out);
848  _writeIncrementSpecies(m_cn_list.begin(), m_cn_list.end(), r, out);
849  }
850 
851  //! @deprecated To be removed after Cantera 2.2
852  void writeDecrementSpecies(const std::string& r, std::map<size_t, std::string>& out) {
853  _writeDecrementSpecies(m_c1_list.begin(), m_c1_list.end(), r, out);
854  _writeDecrementSpecies(m_c2_list.begin(), m_c2_list.end(), r, out);
855  _writeDecrementSpecies(m_c3_list.begin(), m_c3_list.end(), r, out);
856  _writeDecrementSpecies(m_cn_list.begin(), m_cn_list.end(), r, out);
857  }
858 
859  //! @deprecated To be removed after Cantera 2.2
860  void writeIncrementReaction(const std::string& r, std::map<size_t, std::string>& out) {
861  _writeIncrementReaction(m_c1_list.begin(), m_c1_list.end(), r, out);
862  _writeIncrementReaction(m_c2_list.begin(), m_c2_list.end(), r, out);
863  _writeIncrementReaction(m_c3_list.begin(), m_c3_list.end(), r, out);
864  _writeIncrementReaction(m_cn_list.begin(), m_cn_list.end(), r, out);
865  }
866 
867  //! @deprecated To be removed after Cantera 2.2
868  void writeDecrementReaction(const std::string& r, std::map<size_t, std::string>& out) {
869  _writeDecrementReaction(m_c1_list.begin(), m_c1_list.end(), r, out);
870  _writeDecrementReaction(m_c2_list.begin(), m_c2_list.end(), r, out);
871  _writeDecrementReaction(m_c3_list.begin(), m_c3_list.end(), r, out);
872  _writeDecrementReaction(m_cn_list.begin(), m_cn_list.end(), r, out);
873  }
874 
875  //! @deprecated To be removed after Cantera 2.2
876  void writeMultiply(const std::string& r, std::map<size_t, std::string>& out) {
877  _writeMultiply(m_c1_list.begin(), m_c1_list.end(), r, out);
878  _writeMultiply(m_c2_list.begin(), m_c2_list.end(), r, out);
879  _writeMultiply(m_c3_list.begin(), m_c3_list.end(), r, out);
880  _writeMultiply(m_cn_list.begin(), m_cn_list.end(), r, out);
881  }
882 
883 private:
884  std::vector<C1> m_c1_list;
885  std::vector<C2> m_c2_list;
886  std::vector<C3> m_c3_list;
887  std::vector<C_AnyN> m_cn_list;
888 };
889 
890 }
891 
892 #endif
size_t m_n
Length of the m_ic vector.
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:39
void writeIncrementReaction(const std::string &r, std::map< size_t, std::string > &out)
void writeMultiply(const std::string &r, std::map< size_t, std::string > &out)
void writeDecrementReaction(const std::string &r, std::map< size_t, std::string > &out)
static void _writeIncrementSpecies(InputIter begin, InputIter end, const std::string &r, std::map< size_t, std::string > &out)
size_t m_rxn
Reaction index -> index into the ROP vector.
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
size_t m_ic0
Species number.
void writeIncrementReaction(const std::string &r, std::map< size_t, std::string > &out)
static std::string fmt(const std::string &r, size_t n)
void writeMultiply(const std::string &r, std::map< size_t, std::string > &out)
void writeDecrementSpecies(const std::string &r, std::map< size_t, std::string > &out)
void writeIncrementReaction(const std::string &r, std::map< size_t, std::string > &out)
void writeIncrementSpecies(const std::string &r, std::map< size_t, std::string > &out)
static void _writeIncrementReaction(InputIter begin, InputIter end, const std::string &r, std::map< size_t, std::string > &out)
void writeDecrementReaction(const std::string &r, std::map< size_t, std::string > &out)
size_t m_rxn
Reaction number.
void writeMultiply(const std::string &r, std::map< size_t, std::string > &out)
size_t m_ic0
Species index -> index into the species vector for the two species.
static void _writeDecrementReaction(InputIter begin, InputIter end, const std::string &r, std::map< size_t, std::string > &out)
Handles two species in a single reaction.
std::vector< size_t > m_ic
Vector of species which are involved with this stoichiometric manager calculations.
Handles any number of species in a reaction, including fractional stoichiometric coefficients, and arbitrary reaction orders.
std::string fp2str(const double x, const std::string &fmt)
Convert a double into a c++ string.
Definition: stringUtils.cpp:28
void multiply(const DenseMatrix &A, const double *const b, double *const prod)
Multiply A*b and return the result in prod. Uses BLAS routine DGEMV.
void writeDecrementSpecies(const std::string &r, std::map< size_t, std::string > &out)
void writeDecrementReaction(const std::string &r, std::map< size_t, std::string > &out)
void writeIncrementReaction(const std::string &r, std::map< size_t, std::string > &out)
void writeIncrementSpecies(const std::string &r, std::map< size_t, std::string > &out)
void writeMultiply(const std::string &r, std::map< size_t, std::string > &out)
Handles one species in a reaction.
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
void writeDecrementSpecies(const std::string &r, std::map< size_t, std::string > &out)
Contains declarations for string manipulation functions within Cantera.
Handles three species in a reaction.
void writeDecrementReaction(const std::string &r, std::map< size_t, std::string > &out)
static void _writeMultiply(InputIter begin, InputIter end, const std::string &r, std::map< size_t, std::string > &out)
void writeIncrementSpecies(const std::string &r, std::map< size_t, std::string > &out)
size_t m_rxn
ID of the reaction corresponding to this stoichiometric manager.
vector_fp m_stoich
Stoichiometric coefficients for the reaction, reactant or product side.
static void _writeDecrementSpecies(InputIter begin, InputIter end, const std::string &r, std::map< size_t, std::string > &out)
void writeIncrementSpecies(const std::string &r, std::map< size_t, std::string > &out)
vector_fp m_order
Reaction orders for the reaction.
void writeDecrementSpecies(const std::string &r, std::map< size_t, std::string > &out)
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...