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