Cantera 2.6.0
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
7
8#ifndef CT_STOICH_MGR_H
9#define CT_STOICH_MGR_H
10
11#include "cantera/numerics/eigen_sparse.h"
13
14namespace 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 \a k in reaction \a 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
124/**
125 * Handles one species in a reaction.
126 * See @ref Stoichiometry
127 * @ingroup Stoichiometry
128 * @internal
129 */
130class C1
131{
132public:
133 C1(size_t rxn = 0, size_t ic0 = 0) :
134 m_rxn(rxn),
135 m_ic0(ic0) {
136 }
137
138 void incrementSpecies(const doublereal* R, doublereal* S) const {
139 S[m_ic0] += R[m_rxn];
140 }
141
142 void decrementSpecies(const doublereal* R, doublereal* S) const {
143 S[m_ic0] -= R[m_rxn];
144 }
145
146 void multiply(const doublereal* S, doublereal* R) const {
147 R[m_rxn] *= S[m_ic0];
148 }
149
150 void incrementReaction(const doublereal* S, doublereal* R) const {
151 R[m_rxn] += S[m_ic0];
152 }
153
154 void decrementReaction(const doublereal* S, doublereal* R) const {
155 R[m_rxn] -= S[m_ic0];
156 }
157
158 void resizeCoeffs(const std::map<std::pair<size_t, size_t>, size_t>& indices)
159 {
160 m_jc0 = indices.at({m_rxn, m_ic0});
161 }
162
163 void derivatives(const double* S, const double* R, vector_fp& jac) const
164 {
165 jac[m_jc0] += R[m_rxn]; // index (m_ic0, m_rxn)
166 }
167
168
169 void scale(const double* R, double* out, double factor) const
170 {
171 out[m_rxn] = R[m_rxn] * factor;
172 }
173
174private:
175 //! Reaction number
176 size_t m_rxn;
177 //! Species number
178 size_t m_ic0;
179
180 size_t m_jc0; //!< Index in derivative triplet vector
181};
182
183/**
184 * Handles two species in a single reaction.
185 * See @ref Stoichiometry
186 * @ingroup Stoichiometry
187 */
188class C2
189{
190public:
191 C2(size_t rxn = 0, size_t ic0 = 0, size_t ic1 = 0)
192 : m_rxn(rxn), m_ic0(ic0), m_ic1(ic1) {}
193
194 void incrementSpecies(const doublereal* R, doublereal* S) const {
195 S[m_ic0] += R[m_rxn];
196 S[m_ic1] += R[m_rxn];
197 }
198
199 void decrementSpecies(const doublereal* R, doublereal* S) const {
200 S[m_ic0] -= R[m_rxn];
201 S[m_ic1] -= R[m_rxn];
202 }
203
204 void multiply(const doublereal* S, doublereal* R) const {
205 if (S[m_ic0] < 0 && S[m_ic1] < 0) {
206 R[m_rxn] = 0;
207 } else {
208 R[m_rxn] *= S[m_ic0] * S[m_ic1];
209 }
210 }
211
212 void incrementReaction(const doublereal* S, doublereal* R) const {
213 R[m_rxn] += S[m_ic0] + S[m_ic1];
214 }
215
216 void decrementReaction(const doublereal* S, doublereal* R) const {
217 R[m_rxn] -= (S[m_ic0] + S[m_ic1]);
218 }
219
220 void resizeCoeffs(const std::map<std::pair<size_t, size_t>, size_t>& indices)
221 {
222 m_jc0 = indices.at({m_rxn, m_ic0});
223 m_jc1 = indices.at({m_rxn, m_ic1});
224 }
225
226 void derivatives(const double* S, const double* R, vector_fp& jac) const
227 {
228 if (S[m_ic1] > 0) {
229 jac[m_jc0] += R[m_rxn] * S[m_ic1]; // index (m_ic0, m_rxn)
230 }
231 if (S[m_ic0] > 0) {
232 jac[m_jc1] += R[m_rxn] * S[m_ic0]; // index (m_ic1, m_rxn)
233 }
234 }
235
236 void scale(const double* R, double* out, double factor) const
237 {
238 out[m_rxn] = 2 * R[m_rxn] * factor;
239 }
240
241private:
242 //! Reaction index -> index into the ROP vector
243 size_t m_rxn;
244
245 //! Species index -> index into the species vector for the two species.
246 size_t m_ic0, m_ic1;
247
248 size_t m_jc0, m_jc1; //!< Indices in derivative triplet vector
249};
250
251/**
252 * Handles three species in a reaction.
253 * See @ref Stoichiometry
254 * @ingroup Stoichiometry
255 */
256class C3
257{
258public:
259 C3(size_t rxn = 0, size_t ic0 = 0, size_t ic1 = 0, size_t ic2 = 0)
260 : m_rxn(rxn), m_ic0(ic0), m_ic1(ic1), m_ic2(ic2) {}
261
262 void incrementSpecies(const doublereal* R, doublereal* S) const {
263 S[m_ic0] += R[m_rxn];
264 S[m_ic1] += R[m_rxn];
265 S[m_ic2] += 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 S[m_ic2] -= R[m_rxn];
272 }
273
274 void multiply(const doublereal* S, doublereal* R) const {
275 if ((S[m_ic0] < 0 && (S[m_ic1] < 0 || S[m_ic2] < 0)) ||
276 (S[m_ic1] < 0 && S[m_ic2] < 0)) {
277 R[m_rxn] = 0;
278 } else {
279 R[m_rxn] *= S[m_ic0] * S[m_ic1] * S[m_ic2];
280 }
281 }
282
283 void incrementReaction(const doublereal* S, doublereal* R) const {
284 R[m_rxn] += S[m_ic0] + S[m_ic1] + S[m_ic2];
285 }
286
287 void decrementReaction(const doublereal* S, doublereal* R) const {
288 R[m_rxn] -= (S[m_ic0] + S[m_ic1] + S[m_ic2]);
289 }
290
291 void resizeCoeffs(const std::map<std::pair<size_t, size_t>, size_t>& indices)
292 {
293 m_jc0 = indices.at({m_rxn, m_ic0});
294 m_jc1 = indices.at({m_rxn, m_ic1});
295 m_jc2 = indices.at({m_rxn, m_ic2});
296 }
297
298 void derivatives(const double* S, const double* R, vector_fp& jac) const
299 {
300 if (S[m_ic1] > 0 && S[m_ic2] > 0) {
301 jac[m_jc0] += R[m_rxn] * S[m_ic1] * S[m_ic2];; // index (m_ic0, m_rxn)
302 }
303 if (S[m_ic0] > 0 && S[m_ic2] > 0) {
304 jac[m_jc1] += R[m_rxn] * S[m_ic0] * S[m_ic2]; // index (m_ic1, m_ic1)
305 }
306 if (S[m_ic0] > 0 && S[m_ic1] > 0) {
307 jac[m_jc2] += R[m_rxn] * S[m_ic0] * S[m_ic1]; // index (m_ic2, m_ic2)
308 }
309 }
310
311 void scale(const double* R, double* out, double factor) const
312 {
313 out[m_rxn] = 3 * R[m_rxn] * factor;
314 }
315
316private:
317 size_t m_rxn;
318 size_t m_ic0;
319 size_t m_ic1;
320 size_t m_ic2;
321
322 size_t m_jc0, m_jc1, m_jc2; //!< Indices in derivative triplet vector
323};
324
325/**
326 * Handles any number of species in a reaction, including fractional
327 * stoichiometric coefficients, and arbitrary reaction orders.
328 * See @ref Stoichiometry
329 * @ingroup Stoichiometry
330 */
332{
333public:
334 C_AnyN() :
335 m_n(0),
336 m_rxn(npos) {
337 }
338
339 C_AnyN(size_t rxn, const std::vector<size_t>& ic, const vector_fp& order_, const vector_fp& stoich_) :
340 m_n(0),
341 m_rxn(rxn) {
342 m_n = ic.size();
343 m_ic.resize(m_n);
344 m_jc.resize(m_n, 0);
345 m_order.resize(m_n);
346 m_stoich.resize(m_n);
347 for (size_t n = 0; n < m_n; n++) {
348 m_ic[n] = ic[n];
349 m_order[n] = order_[n];
350 m_stoich[n] = stoich_[n];
351 }
352 }
353
354 void multiply(const doublereal* input, doublereal* output) const {
355 for (size_t n = 0; n < m_n; n++) {
356 double order = m_order[n];
357 if (order != 0.0) {
358 double c = input[m_ic[n]];
359 if (c > 0.0) {
360 output[m_rxn] *= std::pow(c, order);
361 } else {
362 output[m_rxn] = 0.0;
363 }
364 }
365 }
366 }
367
368 void incrementSpecies(const doublereal* input,
369 doublereal* output) const {
370 doublereal x = input[m_rxn];
371 for (size_t n = 0; n < m_n; n++) {
372 output[m_ic[n]] += m_stoich[n]*x;
373 }
374 }
375
376 void decrementSpecies(const doublereal* input,
377 doublereal* output) const {
378 doublereal x = input[m_rxn];
379 for (size_t n = 0; n < m_n; n++) {
380 output[m_ic[n]] -= m_stoich[n]*x;
381 }
382 }
383
384 void incrementReaction(const doublereal* input,
385 doublereal* output) const {
386 for (size_t n = 0; n < m_n; n++) {
387 output[m_rxn] += m_stoich[n]*input[m_ic[n]];
388 }
389 }
390
391 void decrementReaction(const doublereal* input,
392 doublereal* output) const {
393 for (size_t n = 0; n < m_n; n++) {
394 output[m_rxn] -= m_stoich[n]*input[m_ic[n]];
395 }
396 }
397
398 void resizeCoeffs(const std::map<std::pair<size_t, size_t>, size_t>& indices)
399 {
400 for (size_t i = 0; i < m_n; i++) {
401 m_jc[i] = indices.at({m_rxn, m_ic[i]});
402 }
403
404 m_sum_order = 0.;
405 for (size_t n = 0; n < m_n; ++n) {
406 m_sum_order += m_order[n];
407 }
408 }
409
410 void derivatives(const double* S, const double* R, vector_fp& jac) const
411 {
412 for (size_t i = 0; i < m_n; i++) {
413 // calculate derivative
414 double prod = R[m_rxn];
415 double order_i = m_order[i];
416 if (S[m_ic[i]] > 0. && order_i != 0.) {
417 prod *= order_i * std::pow(S[m_ic[i]], order_i - 1);
418 for (size_t j = 0; j < m_n; j++) {
419 if (i != j) {
420 if (S[m_ic[j]] > 0) {
421 prod *= std::pow(S[m_ic[j]], m_order[j]);
422 } else {
423 prod = 0.;
424 break;
425 }
426 }
427 }
428 } else {
429 prod = 0.;
430 }
431 jac[m_jc[i]] += prod;
432 }
433 }
434
435 void scale(const double* R, double* out, double factor) const
436 {
437 out[m_rxn] = m_sum_order * R[m_rxn] * factor;
438 }
439
440private:
441 //! Length of the m_ic vector
442 /*!
443 * This is the number of species which participate in the reaction order
444 * and stoichiometric coefficient vectors for the reactant or product
445 * description of the reaction.
446 */
447 size_t m_n;
448
449 //! ID of the reaction corresponding to this stoichiometric manager
450 /*!
451 * This is used within the interface to select the array position to read
452 * and write to Normally this is associated with the reaction number in an
453 * array of quantities indexed by the reaction number, for example, ROP[irxn].
454 */
455 size_t m_rxn;
456
457 //! Vector of species which are involved with this stoichiometric manager
458 //! calculations
459 /*!
460 * This is an integer list of species which participate in either the
461 * reaction order matrix or the stoichiometric order matrix for this
462 * reaction, m_rxn.
463 */
464 std::vector<size_t> m_ic;
465
466 //! Reaction orders for the reaction
467 /*!
468 * This is either for the reactants or products. Length = m_n. Species
469 * number, m_ic[n], has a reaction order of m_order[n].
470 */
472
473 //! Stoichiometric coefficients for the reaction, reactant or product side.
474 /*!
475 * This is either for the reactants or products. Length = m_n. Species
476 * number m_ic[m], has a stoichiometric coefficient of m_stoich[n].
477 */
479
480 std::vector<size_t> m_jc; //!< Indices in derivative triplet vector
481
482 double m_sum_order; //!< Sum of reaction order vector
483};
484
485template<class InputIter, class Vec1, class Vec2>
486inline static void _multiply(InputIter begin, InputIter end,
487 const Vec1& input, Vec2& output)
488{
489 for (; begin != end; ++begin) {
490 begin->multiply(input, output);
491 }
492}
493
494template<class InputIter, class Vec1, class Vec2>
495inline static void _incrementSpecies(InputIter begin,
496 InputIter end, const Vec1& input, Vec2& output)
497{
498 for (; begin != end; ++begin) {
499 begin->incrementSpecies(input, output);
500 }
501}
502
503template<class InputIter, class Vec1, class Vec2>
504inline static void _decrementSpecies(InputIter begin,
505 InputIter end, const Vec1& input, Vec2& output)
506{
507 for (; begin != end; ++begin) {
508 begin->decrementSpecies(input, output);
509 }
510}
511
512template<class InputIter, class Vec1, class Vec2>
513inline static void _incrementReactions(InputIter begin,
514 InputIter end, const Vec1& input, Vec2& output)
515{
516 for (; begin != end; ++begin) {
517 begin->incrementReaction(input, output);
518 }
519}
520
521template<class InputIter, class Vec1, class Vec2>
522inline static void _decrementReactions(InputIter begin,
523 InputIter end, const Vec1& input, Vec2& output)
524{
525 for (; begin != end; ++begin) {
526 begin->decrementReaction(input, output);
527 }
528}
529
530template<class InputIter, class Indices>
531inline static void _resizeCoeffs(InputIter begin, InputIter end, Indices& ix)
532{
533 for (; begin != end; ++begin) {
534 begin->resizeCoeffs(ix);
535 }
536}
537
538template<class InputIter, class Vec1, class Vec2, class Vec3>
539inline static void _derivatives(InputIter begin, InputIter end,
540 const Vec1& S, const Vec2& R, Vec3& jac)
541{
542 for (; begin != end; ++begin) {
543 begin->derivatives(S, R, jac);
544 }
545}
546
547template<class InputIter, class Vec1, class Vec2>
548inline static void _scale(InputIter begin, InputIter end,
549 const Vec1& in, Vec2& out, double factor)
550{
551 for (; begin != end; ++begin) {
552 begin->scale(in, out, factor);
553 }
554}
555
556/*!
557 * This class handles operations involving the stoichiometric coefficients on
558 * one side of a reaction (reactant or product) for a set of reactions
559 * comprising a reaction mechanism. This class is used by class Kinetics, which
560 * contains three instances of this class (one to handle operations on the
561 * reactions, one for the products of reversible reactions, and one for the
562 * products of irreversible reactions).
563 *
564 * This class is designed for use with elementary reactions, or at least ones
565 * with integral stoichiometric coefficients. Let \f$M(i) \f$ be the number of
566 * molecules on the product or reactant side of reaction number i.
567 * \f[
568 * r_i = \sum_m^{M_i} s_{k_{m,i}}
569 * \f]
570 * To understand the operations performed by this class, let \f$N_{k,i}\f$
571 * denote the stoichiometric coefficient of species k on one side (reactant or
572 * product) in reaction i. Then \b N is a sparse K by I matrix of stoichiometric
573 * coefficients.
574 *
575 * The following matrix operations may be carried out with a vector S of length
576 * K, and a vector R of length I:
577 *
578 * - \f$S = S + N R\f$ (incrementSpecies)
579 * - \f$S = S - N R\f$ (decrementSpecies)
580 * - \f$R = R + N^T S \f$ (incrementReaction)
581 * - \f$R = R - N^T S \f$ (decrementReaction)
582 *
583 * The actual implementation, however, does not compute these quantities by
584 * matrix multiplication. A faster algorithm is used that makes use of the fact
585 * that the \b integer-valued N matrix is very sparse, and the non-zero terms
586 * are small positive integers.
587 * \f[
588 * S_k = R_{i1} + \dots + R_{iM}
589 * \f]
590 * where M is the number of molecules, and \f$i(m) \f$ is the
591 * See @ref Stoichiometry
592 * @ingroup Stoichiometry
593 */
595{
596public:
597 /**
598 * Constructor for the StoichManagerN class.
599 *
600 * @internal Consider adding defaulted entries here that supply the total
601 * number of reactions in the mechanism and the total number of species
602 * in the species list. Then, we could use those numbers to provide
603 * error checks during the construction of the object. Those numbers
604 * would also provide some clarity to the purpose and utility of this
605 * class.
606 *
607 * DGG - the problem is that the number of reactions and species are not
608 * known initially.
609 */
611 m_stoichCoeffs.setZero();
612 m_stoichCoeffs.resize(0, 0);
613 }
614
615 //! Resize the sparse coefficient matrix
616 void resizeCoeffs(size_t nSpc, size_t nRxn)
617 {
618 size_t nCoeffs = m_coeffList.size();
619
620 // Stoichiometric coefficient matrix
621 m_stoichCoeffs.resize(nSpc, nRxn);
622 m_stoichCoeffs.reserve(nCoeffs);
623 m_stoichCoeffs.setFromTriplets(m_coeffList.begin(), m_coeffList.end());
624
625 // Set up outer/inner indices for mapped derivative output
626 Eigen::SparseMatrix<double> tmp = m_stoichCoeffs.transpose();
627 m_outerIndices.resize(nSpc + 1); // number of columns + 1
628 for (int i = 0; i < tmp.outerSize() + 1; i++) {
629 m_outerIndices[i] = tmp.outerIndexPtr()[i];
630 }
631 m_innerIndices.resize(nCoeffs);
632 for (size_t n = 0; n < nCoeffs; n++) {
633 m_innerIndices[n] = tmp.innerIndexPtr()[n];
634 }
635 m_values.resize(nCoeffs, 0.);
636
637 // Set up index pairs for derivatives
638 std::map<std::pair<size_t, size_t>, size_t> indices;
639 size_t n = 0;
640 for (int i = 0; i < tmp.outerSize(); i++) {
641 for (Eigen::SparseMatrix<double>::InnerIterator it(tmp, i); it; ++it) {
642 indices[{static_cast<size_t>(it.row()),
643 static_cast<size_t>(it.col())}] = n++;
644 }
645 }
646 // update reaction setup
647 _resizeCoeffs(m_c1_list.begin(), m_c1_list.end(), indices);
648 _resizeCoeffs(m_c2_list.begin(), m_c2_list.end(), indices);
649 _resizeCoeffs(m_c3_list.begin(), m_c3_list.end(), indices);
650 _resizeCoeffs(m_cn_list.begin(), m_cn_list.end(), indices);
651
653 }
654
655 /**
656 * Add a single reaction to the list of reactions that this stoichiometric
657 * manager object handles.
658 *
659 * This function is the same as the add() function below. However, the order
660 * of each species in the power list expression is set to one automatically.
661 */
662 void add(size_t rxn, const std::vector<size_t>& k) {
663 vector_fp order(k.size(), 1.0);
664 vector_fp stoich(k.size(), 1.0);
666 }
667
668 void add(size_t rxn, const std::vector<size_t>& k, const vector_fp& order) {
669 vector_fp stoich(k.size(), 1.0);
671 }
672
673 //! Add a single reaction to the list of reactions that this
674 //! stoichiometric manager object handles.
675 /*!
676 * @param rxn Reaction index of the current reaction. This is used as an
677 * index into vectors which have length n_total_rxn.
678 * @param k This is a vector of integer values specifying the species
679 * indices. The length of this vector species the number of different
680 * species in the description. The value of the entries are the species
681 * indices. These are used as indexes into vectors which have length
682 * n_total_species.
683 * @param order This is a vector of the same length as vector k. The order
684 * is used for the routine power(), which produces a power law
685 * expression involving the species vector.
686 * @param stoich This is used to handle fractional stoichiometric
687 * coefficients on the product side of irreversible reactions.
688 */
689 void add(size_t rxn, const std::vector<size_t>& k, const vector_fp& order,
690 const vector_fp& stoich) {
691 if (order.size() != k.size()) {
692 throw CanteraError("StoichManagerN::add()", "size of order and species arrays differ");
693 }
694 if (stoich.size() != k.size()) {
695 throw CanteraError("StoichManagerN::add()", "size of stoich and species arrays differ");
696 }
697 bool frac = false;
698 for (size_t n = 0; n < stoich.size(); n++) {
699 m_coeffList.emplace_back(
700 static_cast<int>(k[n]), static_cast<int>(rxn), stoich[n]);
701 if (fmod(stoich[n], 1.0) || stoich[n] != order[n]) {
702 frac = true;
703 }
704 }
705 if (frac || k.size() > 3) {
706 m_cn_list.emplace_back(rxn, k, order, stoich);
707 } else {
708 // Try to express the reaction with unity stoichiometric
709 // coefficients (by repeating species when necessary) so that the
710 // simpler 'multiply' function can be used to compute the rate
712 std::vector<size_t> kRep;
713 for (size_t n = 0; n < k.size(); n++) {
714 for (size_t i = 0; i < stoich[n]; i++) {
715 kRep.push_back(k[n]);
716 }
717 }
718
719 switch (kRep.size()) {
720 case 1:
721 m_c1_list.emplace_back(rxn, kRep[0]);
722 break;
723 case 2:
724 m_c2_list.emplace_back(rxn, kRep[0], kRep[1]);
725 break;
726 case 3:
727 m_c3_list.emplace_back(rxn, kRep[0], kRep[1], kRep[2]);
728 break;
729 default:
730 m_cn_list.emplace_back(rxn, k, order, stoich);
731 }
732 }
734 }
735
736 void multiply(const doublereal* input, doublereal* output) const {
737 _multiply(m_c1_list.begin(), m_c1_list.end(), input, output);
738 _multiply(m_c2_list.begin(), m_c2_list.end(), input, output);
739 _multiply(m_c3_list.begin(), m_c3_list.end(), input, output);
740 _multiply(m_cn_list.begin(), m_cn_list.end(), input, output);
741 }
742
743 void incrementSpecies(const doublereal* input, doublereal* output) const {
744 _incrementSpecies(m_c1_list.begin(), m_c1_list.end(), input, output);
745 _incrementSpecies(m_c2_list.begin(), m_c2_list.end(), input, output);
746 _incrementSpecies(m_c3_list.begin(), m_c3_list.end(), input, output);
747 _incrementSpecies(m_cn_list.begin(), m_cn_list.end(), input, output);
748 }
749
750 void decrementSpecies(const doublereal* input, doublereal* output) const {
751 _decrementSpecies(m_c1_list.begin(), m_c1_list.end(), input, output);
752 _decrementSpecies(m_c2_list.begin(), m_c2_list.end(), input, output);
753 _decrementSpecies(m_c3_list.begin(), m_c3_list.end(), input, output);
754 _decrementSpecies(m_cn_list.begin(), m_cn_list.end(), input, output);
755 }
756
757 void incrementReactions(const doublereal* input, doublereal* output) const {
758 _incrementReactions(m_c1_list.begin(), m_c1_list.end(), input, output);
759 _incrementReactions(m_c2_list.begin(), m_c2_list.end(), input, output);
760 _incrementReactions(m_c3_list.begin(), m_c3_list.end(), input, output);
761 _incrementReactions(m_cn_list.begin(), m_cn_list.end(), input, output);
762 }
763
764 void decrementReactions(const doublereal* input, doublereal* output) const {
765 _decrementReactions(m_c1_list.begin(), m_c1_list.end(), input, output);
766 _decrementReactions(m_c2_list.begin(), m_c2_list.end(), input, output);
767 _decrementReactions(m_c3_list.begin(), m_c3_list.end(), input, output);
768 _decrementReactions(m_cn_list.begin(), m_cn_list.end(), input, output);
769 }
770
771 //! Return matrix containing stoichiometric coefficients
772 const Eigen::SparseMatrix<double>& stoichCoeffs() const
773 {
775 // This can happen if a user overrides default behavior:
776 // Kinetics::resizeReactions is not called after adding reactions via
777 // Kinetics::addReaction with the 'resize' flag set to 'false'
778 throw CanteraError("StoichManagerN::stoichCoeffs", "The object "
779 "is not fully configured; make sure to call resizeCoeffs().");
780 }
781 return m_stoichCoeffs;
782 }
783
784 //! Calculate derivatives with respect to species concentrations.
785 /*!
786 * The species derivative is the term of the Jacobian that accounts for
787 * species products in the law of mass action. As StoichManagerN does not account
788 * for third bodies or rate constants that depend on species concentrations,
789 * corresponding derivatives are not included.
790 *
791 * @param conc Species concentration.
792 * @param rates Rates-of-progress.
793 */
794 Eigen::SparseMatrix<double> derivatives(const double* conc, const double* rates)
795 {
796 // calculate derivative entries using known sparse storage order
797 std::fill(m_values.begin(), m_values.end(), 0.);
798 _derivatives(m_c1_list.begin(), m_c1_list.end(), conc, rates, m_values);
799 _derivatives(m_c2_list.begin(), m_c2_list.end(), conc, rates, m_values);
800 _derivatives(m_c3_list.begin(), m_c3_list.end(), conc, rates, m_values);
801 _derivatives(m_cn_list.begin(), m_cn_list.end(), conc, rates, m_values);
802
803 return Eigen::Map<Eigen::SparseMatrix<double>>(
804 m_stoichCoeffs.cols(), m_stoichCoeffs.rows(), m_values.size(),
805 m_outerIndices.data(), m_innerIndices.data(), m_values.data());
806 }
807
808 //! Scale input by reaction order and factor
809 void scale(const double* in, double* out, double factor) const
810 {
811 _scale(m_c1_list.begin(), m_c1_list.end(), in, out, factor);
812 _scale(m_c2_list.begin(), m_c2_list.end(), in, out, factor);
813 _scale(m_c3_list.begin(), m_c3_list.end(), in, out, factor);
814 _scale(m_cn_list.begin(), m_cn_list.end(), in, out, factor);
815 }
816
817private:
818 bool m_ready; //!< Boolean flag indicating whether object is fully configured
819
820 std::vector<C1> m_c1_list;
821 std::vector<C2> m_c2_list;
822 std::vector<C3> m_c3_list;
823 std::vector<C_AnyN> m_cn_list;
824
825 //! Sparse matrices for stoichiometric coefficients
826 SparseTriplets m_coeffList;
827 Eigen::SparseMatrix<double> m_stoichCoeffs;
828
829 //! Storage indicies used to build derivatives
830 std::vector<int> m_outerIndices;
831 std::vector<int> m_innerIndices;
832 vector_fp m_values;
833};
834
835}
836
837#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.
std::vector< size_t > m_ic
Vector of species which are involved with this stoichiometric manager calculations.
double m_sum_order
Sum of reaction order vector.
std::vector< size_t > m_jc
Indices in derivative triplet vector.
vector_fp m_order
Reaction orders for the reaction.
size_t m_n
Length of the m_ic vector.
vector_fp m_stoich
Stoichiometric coefficients for the reaction, reactant or product side.
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:61
void add(size_t rxn, const std::vector< size_t > &k)
Add a single reaction to the list of reactions that this stoichiometric manager object handles.
Boolean flag indicating whether object is fully configured.
Eigen::SparseMatrix< double > derivatives(const double *conc, const double *rates)
Calculate derivatives with respect to species concentrations.
SparseTriplets m_coeffList
Sparse matrices for stoichiometric coefficients.
std::vector< int > m_outerIndices
Storage indicies used to build derivatives.
void add(size_t rxn, const std::vector< size_t > &k, const vector_fp &order, const vector_fp &stoich)
Add a single reaction to the list of reactions that this stoichiometric manager object handles.
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.
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.h:29
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:192
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:184