8 #ifndef CT_STOICH_MGR_H
9 #define CT_STOICH_MGR_H
11 #include "cantera/numerics/eigen_sparse.h"
134 C1(
size_t rxn = 0,
size_t ic0 = 0) :
139 void incrementSpecies(
const double* R,
double* S)
const {
143 void decrementSpecies(
const double* R,
double* S)
const {
147 void multiply(
const double* S,
double* R)
const {
151 void incrementReaction(
const double* S,
double* R)
const {
155 void decrementReaction(
const double* S,
double* R)
const {
159 void resizeCoeffs(
const map<pair<size_t, size_t>,
size_t>& indices)
164 void derivatives(
const double* S,
const double* R, vector<double>& jac)
const
170 void scale(
const double* R,
double* out,
double factor)
const
192 C2(
size_t rxn = 0,
size_t ic0 = 0,
size_t ic1 = 0)
195 void incrementSpecies(
const double* R,
double* S)
const {
197 S[m_ic1] += R[
m_rxn];
200 void decrementSpecies(
const double* R,
double* S)
const {
202 S[m_ic1] -= R[
m_rxn];
205 void multiply(
const double* S,
double* R)
const {
206 if (S[
m_ic0] < 0 && S[m_ic1] < 0) {
213 void incrementReaction(
const double* S,
double* R)
const {
217 void decrementReaction(
const double* S,
double* R)
const {
221 void resizeCoeffs(
const map<pair<size_t, size_t>,
size_t>& indices)
227 void derivatives(
const double* S,
const double* R, vector<double>& jac)
const
230 jac[m_jc0] += R[
m_rxn] * S[m_ic1];
237 void scale(
const double* R,
double* out,
double factor)
const
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) {}
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];
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];
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)) {
280 R[m_rxn] *= S[m_ic0] * S[m_ic1] * S[m_ic2];
284 void incrementReaction(
const double* S,
double* R)
const {
285 R[m_rxn] += S[m_ic0] + S[m_ic1] + S[m_ic2];
288 void decrementReaction(
const double* S,
double* R)
const {
289 R[m_rxn] -= (S[m_ic0] + S[m_ic1] + S[m_ic2]);
292 void resizeCoeffs(
const map<pair<size_t, size_t>,
size_t>& indices)
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});
299 void derivatives(
const double* S,
const double* R, vector<double>& jac)
const
301 if (S[m_ic1] > 0 && S[m_ic2] > 0) {
302 jac[m_jc0] += R[m_rxn] * S[m_ic1] * S[m_ic2];;
304 if (S[m_ic0] > 0 && S[m_ic2] > 0) {
305 jac[m_jc1] += R[m_rxn] * S[m_ic0] * S[m_ic2];
307 if (S[m_ic0] > 0 && S[m_ic1] > 0) {
308 jac[
m_jc2] += R[m_rxn] * S[m_ic0] * S[m_ic1];
312 void scale(
const double* R,
double* out,
double factor)
const
314 out[m_rxn] = 3 * R[m_rxn] * factor;
337 C_AnyN(
size_t rxn,
const vector<size_t>& ic,
338 const vector<double>& order_,
const vector<double>& stoich_) :
345 for (
size_t n = 0; n <
m_n; n++) {
352 void multiply(
const double* input,
double* output)
const {
353 for (
size_t n = 0; n <
m_n; n++) {
356 double c = input[
m_ic[n]];
358 output[
m_rxn] *= std::pow(c, order);
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++) {
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++) {
380 void incrementReaction(
const double* input,
double* output)
const {
381 for (
size_t n = 0; n <
m_n; n++) {
386 void decrementReaction(
const double* input,
double* output)
const {
387 for (
size_t n = 0; n <
m_n; n++) {
392 void resizeCoeffs(
const map<pair<size_t, size_t>,
size_t>& indices)
394 for (
size_t i = 0; i <
m_n; i++) {
399 for (
size_t n = 0; n <
m_n; ++n) {
404 void derivatives(
const double* S,
const double* R, vector<double>& jac)
const
406 for (
size_t i = 0; i <
m_n; i++) {
408 double prod = R[
m_rxn];
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++) {
414 if (S[
m_ic[j]] > 0) {
425 jac[
m_jc[i]] += prod;
429 void scale(
const double* R,
double* out,
double factor)
const
479 template<
class InputIter,
class Vec1,
class Vec2>
480 inline static void _multiply(InputIter begin, InputIter end,
481 const Vec1& input, Vec2& output)
483 for (; begin != end; ++begin) {
484 begin->multiply(input, output);
488 template<
class InputIter,
class Vec1,
class Vec2>
489 inline static void _incrementSpecies(InputIter begin,
490 InputIter end,
const Vec1& input, Vec2& output)
492 for (; begin != end; ++begin) {
493 begin->incrementSpecies(input, output);
497 template<
class InputIter,
class Vec1,
class Vec2>
498 inline static void _decrementSpecies(InputIter begin,
499 InputIter end,
const Vec1& input, Vec2& output)
501 for (; begin != end; ++begin) {
502 begin->decrementSpecies(input, output);
506 template<
class InputIter,
class Vec1,
class Vec2>
507 inline static void _incrementReactions(InputIter begin,
508 InputIter end,
const Vec1& input, Vec2& output)
510 for (; begin != end; ++begin) {
511 begin->incrementReaction(input, output);
515 template<
class InputIter,
class Vec1,
class Vec2>
516 inline static void _decrementReactions(InputIter begin,
517 InputIter end,
const Vec1& input, Vec2& output)
519 for (; begin != end; ++begin) {
520 begin->decrementReaction(input, output);
524 template<
class InputIter,
class Indices>
525 inline static void _resizeCoeffs(InputIter begin, InputIter end, Indices& ix)
527 for (; begin != end; ++begin) {
528 begin->resizeCoeffs(ix);
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)
536 for (; begin != end; ++begin) {
537 begin->derivatives(S, R, jac);
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)
545 for (; begin != end; ++begin) {
546 begin->scale(in, out, factor);
595 m_stoichCoeffs.setZero();
596 m_stoichCoeffs.resize(0, 0);
605 m_stoichCoeffs.resize(nSpc, nRxn);
606 m_stoichCoeffs.reserve(nCoeffs);
610 Eigen::SparseMatrix<double> tmp = m_stoichCoeffs.transpose();
612 for (
int i = 0; i < tmp.outerSize() + 1; i++) {
615 m_innerIndices.resize(nCoeffs);
616 for (
size_t n = 0; n < nCoeffs; n++) {
617 m_innerIndices[n] = tmp.innerIndexPtr()[n];
619 m_values.resize(nCoeffs, 0.);
622 map<pair<size_t, size_t>,
size_t> indices;
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++;
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);
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);
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);
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()) {
677 "StoichManagerN::add()",
"size of order and species arrays differ");
679 if (stoich.size() != k.size()) {
681 "StoichManagerN::add()",
"size of stoich and species arrays differ");
684 for (
size_t n = 0; n < stoich.size(); n++) {
686 static_cast<int>(k[n]),
static_cast<int>(rxn), stoich[n]);
687 if (fmod(stoich[n], 1.0) || stoich[n] != order[n]) {
691 if (frac || k.size() > 3) {
692 m_cn_list.emplace_back(rxn, k, order, stoich);
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]);
705 switch (kRep.size()) {
707 m_c1_list.emplace_back(rxn, kRep[0]);
710 m_c2_list.emplace_back(rxn, kRep[0], kRep[1]);
713 m_c3_list.emplace_back(rxn, kRep[0], kRep[1], kRep[2]);
716 m_cn_list.emplace_back(rxn, k, order, stoich);
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);
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);
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);
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);
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);
764 throw CanteraError(
"StoichManagerN::stoichCoeffs",
"The object "
765 "is not fully configured; make sure to call resizeCoeffs().");
767 return m_stoichCoeffs;
780 Eigen::SparseMatrix<double>
derivatives(
const double* conc,
const double* rates)
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);
789 return Eigen::Map<Eigen::SparseMatrix<double>>(
790 m_stoichCoeffs.cols(), m_stoichCoeffs.rows(), m_values.size(),
795 void scale(
const double* in,
double* out,
double factor)
const
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);
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;
813 Eigen::SparseMatrix<double> m_stoichCoeffs;
817 vector<int> m_innerIndices;
818 vector<double> m_values;
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.
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.
const size_t npos
index returned by functions to indicate "no position"