Cantera  3.1.0a1
Loading...
Searching...
No Matches
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
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 @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 */
131class C1
132{
133public:
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
175private:
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 */
189class C2
190{
191public:
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
242private:
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 */
257class C3
258{
259public:
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
317private:
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 */
333{
334public:
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
434private:
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
479template<class InputIter, class Vec1, class Vec2>
480inline 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
488template<class InputIter, class Vec1, class Vec2>
489inline 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
497template<class InputIter, class Vec1, class Vec2>
498inline 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
506template<class InputIter, class Vec1, class Vec2>
507inline 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
515template<class InputIter, class Vec1, class Vec2>
516inline 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
524template<class InputIter, class Indices>
525inline static void _resizeCoeffs(InputIter begin, InputIter end, Indices& ix)
526{
527 for (; begin != end; ++begin) {
528 begin->resizeCoeffs(ix);
529 }
530}
531
532template<class InputIter, class Vec1, class Vec2, class Vec3>
533inline 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
541template<class InputIter, class Vec1, class Vec2>
542inline 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{
590public:
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
803private:
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.
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.
bool m_ready
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.
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