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