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