Cantera  3.1.0a1
Loading...
Searching...
No Matches
ReactionPath.cpp
Go to the documentation of this file.
1/**
2 * @file ReactionPath.cpp
3 * Implementation file for classes used in reaction path analysis.
4 */
5
6// This file is part of Cantera. See License.txt in the top-level directory or
7// at https://cantera.org/license.txt for license and copyright information.
8
12
13#include <boost/algorithm/string.hpp>
14
15namespace ba = boost::algorithm;
16
17using namespace std;
18
19namespace Cantera
20{
21
23{
24 m_paths.push_back(path);
25 if (path->begin() == this) {
26 m_out += path->flow();
27 } else if (path->end() == this) {
28 m_in += path->flow();
29 } else {
30 throw CanteraError("SpeciesNode::addPath", "path added to wrong node");
31 }
32}
33
34void SpeciesNode::printPaths()
35{
36 for (size_t i = 0; i < m_paths.size(); i++) {
37 cout << m_paths[i]->begin()->name << " --> "
38 << m_paths[i]->end()->name << ": "
39 << m_paths[i]->flow() << endl;
40 }
41}
42
44 : m_a(begin), m_b(end)
45{
46 begin->addPath(this);
47 end->addPath(this);
48}
49
50void Path::addReaction(size_t rxnNumber, double value, const string& label)
51{
52 m_rxn[rxnNumber] += value;
53 m_total += value;
54 if (label != "") {
55 m_label[label] += value;
56 }
57}
58
59void Path::writeLabel(ostream& s, double threshold)
60{
61 if (m_label.size() == 0) {
62 return;
63 }
64 double v;
65 for (const auto& [species, value] : m_label) {
66 v = value / m_total;
67 if (m_label.size() == 1) {
68 s << species << "\\l";
69 } else if (v > threshold) {
70 s << species;
71 int percent = int(100*v + 0.5);
72 if (percent < 100) {
73 s << " (" << percent << "%)\\l";
74 } else {
75 s << "\\l";
76 }
77 }
78 }
79}
80
82{
83 // delete the nodes
84 for (const auto& [k, node] : m_nodes) {
85 delete node;
86 }
87
88 // delete the paths
89 size_t nn = nPaths();
90 for (size_t n = 0; n < nn; n++) {
91 delete m_pathlist[n];
92 }
93}
94
95vector<int> ReactionPathDiagram::reactions()
96{
97 double flmax = 0.0;
98 for (size_t i = 0; i < nPaths(); i++) {
99 Path* p = path(i);
100 flmax = std::max(p->flow(), flmax);
101 }
102 m_rxns.clear();
103 for (size_t i = 0; i < nPaths(); i++) {
104 for (const auto& [iRxn, flux] : path(i)->reactionMap()) {
105 double flxratio = flux / flmax;
106 if (flxratio > threshold) {
107 m_rxns.insert(iRxn);
108 }
109 }
110 }
111 vector<int> r;
112 for (const auto& rxn : m_rxns) {
113 r.push_back(int(rxn));
114 }
115 return r;
116}
117
118void ReactionPathDiagram::add(ReactionPathDiagram& d)
119{
120 for (size_t n = 0; n < nPaths(); n++) {
121 Path* p = path(n);
122 size_t k1 = p->begin()->number;
123 size_t k2 = p->end()->number;
124 p->setFlow(p->flow() + d.flow(k1,k2));
125 }
126}
127
128void ReactionPathDiagram::findMajorPaths(double athreshold, size_t lda, double* a)
129{
130 double netmax = 0.0;
131 for (size_t n = 0; n < nNodes(); n++) {
132 for (size_t m = n+1; m < nNodes(); m++) {
133 size_t k1 = m_speciesNumber[n];
134 size_t k2 = m_speciesNumber[m];
135 double fl = fabs(netFlow(k1,k2));
136 netmax = std::max(fl, netmax);
137 }
138 }
139 for (size_t n = 0; n < nNodes(); n++) {
140 for (size_t m = n+1; m < nNodes(); m++) {
141 size_t k1 = m_speciesNumber[n];
142 size_t k2 = m_speciesNumber[m];
143 double fl = fabs(netFlow(k1,k2));
144 if (fl > athreshold*netmax) {
145 a[lda*k1 + k2] = 1;
146 }
147 }
148 }
149}
150
151void ReactionPathDiagram::writeData(ostream& s)
152{
153 s << title << endl;
154 for (size_t i1 = 0; i1 < nNodes(); i1++) {
155 size_t k1 = m_speciesNumber[i1];
156 s << m_nodes[k1]->name << " ";
157 }
158 s << endl;
159 for (size_t i1 = 0; i1 < nNodes(); i1++) {
160 size_t k1 = m_speciesNumber[i1];
161 for (size_t i2 = i1+1; i2 < nNodes(); i2++) {
162 size_t k2 = m_speciesNumber[i2];
163 double f1 = flow(k1, k2);
164 double f2 = flow(k2, k1);
165 s << m_nodes[k1]->name << " " << m_nodes[k2]->name
166 << " " << f1 << " " << -f2 << endl;
167 }
168 }
169}
170
172{
173 double flmax = 0.0;
174 s.precision(3);
175
176 // a directed graph
177 s << "digraph " << name << " {" << endl;
178
179 // the graph will be no larger than x_size, y_size
180 if (x_size > 0.0) {
181 if (y_size < 0.0) {
182 y_size = x_size;
183 }
184 s << "size = \""
185 << x_size << ","
186 << y_size << "\";"
187 << endl;
188 }
189
190 if (dot_options != "") {
191 s << dot_options << endl;
192 }
193
194 // draw paths representing net flows
195 if (flow_type == NetFlow) {
196 // if no scale was specified, normalize net flows by the maximum net
197 // flow
198 if (scale <= 0.0) {
199 for (size_t i1 = 0; i1 < nNodes(); i1++) {
200 size_t k1 = m_speciesNumber[i1];
201 node(k1)->visible = false;
202 for (size_t i2 = i1+1; i2 < nNodes(); i2++) {
203 size_t k2 = m_speciesNumber[i2];
204 double flx = netFlow(k1, k2);
205 if (flx < 0.0) {
206 flx = -flx;
207 }
208 flmax = std::max(flx, flmax);
209 }
210 }
211 } else {
212 flmax = scale;
213 }
214 flmax = std::max(flmax, 1e-10);
215
216 // loop over all unique pairs of nodes
217 for (size_t i1 = 0; i1 < nNodes(); i1++) {
218 size_t k1 = m_speciesNumber[i1];
219 for (size_t i2 = i1+1; i2 < nNodes(); i2++) {
220 size_t k2 = m_speciesNumber[i2];
221 double flx = netFlow(k1, k2);
222 if (m_local != npos && k1 != m_local && k2 != m_local) {
223 flx = 0.0;
224 }
225 if (flx != 0.0) {
226 double flxratio;
227 size_t kbegin, kend;
228 // set beginning and end of the path based on the sign of
229 // the net flow
230 if (flx > 0.0) {
231 kbegin = k1;
232 kend = k2;
233 flxratio = flx/flmax;
234 } else {
235 kbegin = k2;
236 kend = k1;
237 flxratio = -flx/flmax;
238 }
239
240 // write out path specification if the net flow is greater
241 // than the threshold
242 if (flxratio >= threshold) {
243 // make nodes visible
244 node(kbegin)->visible = true;
245 node(kend)->visible = true;
246
247 s << "s" << kbegin << " -> s" << kend;
248 s << "[fontname=\""+m_font+"\", penwidth=";
249
250 if (arrow_width < 0) {
251 double lwidth = 1.0 - 4.0
252 * log10(flxratio/threshold)/log10(threshold) + 1.0;
253 s << lwidth;
254 s << ", arrowsize="
255 << std::min(6.0, 0.5*lwidth);
256 } else {
257 s << arrow_width;
258 s << ", arrowsize=" << flxratio + 1;
259 }
260
261 double hue = 0.7;
262 double bright = 0.9;
263 s << ", color=" << "\"" << hue << ", "
264 << flxratio + 0.5
265 << ", " << bright << "\"" << endl;
266
267 if (flxratio > label_min) {
268 s << ", label=\" " << flxratio;
269 if (show_details) {
270 if (flow(kbegin, kend) > 0.0) {
271 s << "\\l fwd: "
272 << flow(kbegin, kend)/flmax << "\\l";
273 path(kbegin, kend)->writeLabel(s);
274 }
275 if (flow(kend, kbegin) > 0.0) {
276 s << " \\l rev: "
277 << flow(kend,kbegin)/flmax << "\\l";
278 path(kend, kbegin)->writeLabel(s);
279 }
280 }
281 s << "\"";
282 }
283 s << "];" << endl;
284 }
285 }
286 }
287 }
288 } else {
289 if (scale < 0) {
290 for (size_t i = 0; i < nPaths(); i++) {
291 flmax = std::max(path(i)->flow(), flmax);
292 }
293 } else {
294 flmax = scale;
295 }
296
297 for (size_t i = 0; i < nPaths(); i++) {
298 Path* p = path(i);
299 double flxratio = p->flow()/flmax;
300 if (m_local != npos) {
301 if (p->begin()->number != m_local
302 && p->end()->number != m_local) {
303 flxratio = 0.0;
304 }
305 }
306 if (flxratio > threshold) {
307 p->begin()->visible = true;
308 p->end()->visible = true;
309 s << "s" << p->begin()->number
310 << " -> s" << p->end()->number;
311
312 if (arrow_width < 0) {
313 double lwidth = 1.0 - 4.0 * log10(flxratio/threshold)/log10(threshold)
314 + 1.0;
315 s << "[fontname=\""+m_font+"\", penwidth="
316 << lwidth;
317 s << ", arrowsize="
318 << std::min(6.0, 0.5*lwidth);
319 } else {
320 s << ", penwidth="
321 << arrow_width;
322 s << ", arrowsize=" << flxratio + 1;
323 }
324 double hue = 0.7;
325 double bright = 0.9;
326 s << ", color=" << "\"" << hue << ", " << flxratio + 0.5
327 << ", " << bright << "\"" << endl;
328
329 if (flxratio > label_min) {
330 s << ", label = \" " << flxratio;
331 if (show_details) {
332 s << "\\l";
333 p->writeLabel(s);
334 }
335 s << "\"";
336 }
337 s << "];" << endl;
338 }
339 }
340 }
341 s.precision(2);
342 for (const auto& [kSpecies, node] : m_nodes) {
343 if (node->visible) {
344 s << "s" << kSpecies << " [ fontname=\""+m_font+"\", label=\"" << node->name
345 << "\"];" << endl;
346 }
347 }
348 s << " label = " << "\"" << "Scale = "
349 << flmax << "\\l " << title << "\";" << endl;
350 s << " fontname = \""+m_font+"\";" << endl << "}" << endl;
351}
352
353
354void ReactionPathDiagram::addNode(size_t k, const string& nm, double x)
355{
356 if (!m_nodes[k]) {
357 m_nodes[k] = new SpeciesNode;
358 m_nodes[k]->number = k;
359 m_nodes[k]->name = nm;
360 m_nodes[k]->value = x;
361 m_speciesNumber.push_back(k);
362 }
363}
364
365void ReactionPathDiagram::linkNodes(size_t k1, size_t k2, size_t rxn,
366 double value, string legend)
367{
368 Path* ff = m_paths[k1][k2];
369 if (!ff) {
370 ff= new Path(m_nodes[k1], m_nodes[k2]);
371 m_paths[k1][k2] = ff;
372 m_pathlist.push_back(ff);
373 }
374 ff->addReaction(rxn, value, legend);
375 m_rxns.insert(rxn);
376 m_flxmax = std::max(ff->flow(), m_flxmax);
377}
378
379vector<size_t> ReactionPathDiagram::species()
380{
381 return m_speciesNumber;
382}
383
385{
386 m_groups.resize(m_nr);
387 for (size_t i = 0; i < m_nr; i++) { // loop over reactions
388 logfile << endl << "Reaction " << i+1 << ": "
389 << s.reaction(i)->equation();
390
391 if (m_determinate[i]) {
392 logfile << " ... OK." << endl;
393 } else if (m_reac[i].size() == 2 && m_prod[i].size() == 2) {
394 // indices for the two reactants
395 size_t kr0 = m_reac[i][0];
396 size_t kr1 = m_reac[i][1];
397
398 // indices for the two products
399 size_t kp0 = m_prod[i][0];
400 size_t kp1 = m_prod[i][1];
401
402 // references to the Group objects representing the reactants
403 const Group& r0 = m_sgroup[kr0];
404 const Group& r1 = m_sgroup[kr1];
405 const Group& p0 = m_sgroup[kp0];
406 const Group& p1 = m_sgroup[kp1];
407
408 const Group* group_a0=0, *group_b0=0, *group_c0=0,
409 *group_a1=0, *group_b1=0, *group_c1=0;
410 Group b0 = p0 - r0;
411 Group b1 = p1 - r0;
412 if (b0.valid() && b1.valid()) {
413 logfile << " ... ambiguous." << endl;
414 } else if (!b0.valid() && !b1.valid()) {
415 logfile << " ... cannot express as A + BC = AB + C" << endl;
416 } else {
417 logfile << endl;
418 }
419
420 if (b0.valid()) {
421 if (b0.sign() > 0) {
422 group_a0 = &r0;
423 group_b0 = &b0;
424 group_c0 = &p1;
425 m_transfer[i][0][0] = r0;
426 m_transfer[i][1][0] = b0;
427 m_transfer[i][1][1] = p1;
428 } else {
429 group_a0 = &r1;
430 group_c0 = &p0;
431 b0 *= -1;
432 group_b0 = &b0;
433 m_transfer[i][1][1] = r1;
434 m_transfer[i][0][1] = b0;
435 m_transfer[i][0][0] = p0;
436 }
437 logfile << " ";
438 group_a0->fmt(logfile, m_elementSymbols);
439 logfile << " + ";
440 group_b0->fmt(logfile,m_elementSymbols);
441 group_c0->fmt(logfile, m_elementSymbols);
442 logfile << " = ";
443 group_a0->fmt(logfile, m_elementSymbols);
444 group_b0->fmt(logfile, m_elementSymbols);
445 logfile << " + ";
446 group_c0->fmt(logfile, m_elementSymbols);
447 if (b1.valid()) {
448 logfile << " [<= default] " << endl;
449 } else {
450 logfile << endl;
451 }
452 }
453
454 if (b1.valid()) {
455 if (b1.sign() > 0) {
456 group_a1 = &r0;
457 group_b1 = &b1;
458 group_c1 = &p0;
459 if (!b0.valid()) {
460 m_transfer[i][0][1] = r0;
461 m_transfer[i][1][1] = b0;
462 m_transfer[i][1][0] = p0;
463 }
464 } else {
465 group_a1 = &r1;
466 group_c1 = &p1;
467 b1 *= -1;
468 group_b1 = &b1;
469 if (!b0.valid()) {
470 m_transfer[i][1][0] = r1;
471 m_transfer[i][0][0] = b0;
472 m_transfer[i][0][1] = p1;
473 }
474 }
475 logfile << " ";
476 group_a1->fmt(logfile, m_elementSymbols);
477 logfile << " + ";
478 group_b1->fmt(logfile, m_elementSymbols);
479 group_c1->fmt(logfile, m_elementSymbols);
480 logfile << " = ";
481 group_a1->fmt(logfile, m_elementSymbols);
482 group_b1->fmt(logfile, m_elementSymbols);
483 logfile << " + ";
484 group_c1->fmt(logfile, m_elementSymbols);
485 logfile << endl;
486 }
487 } else {
488 logfile << "... cannot parse. [ignored]" << endl;
489 }
490 }
491 return 1;
492}
493
494void ReactionPathBuilder::findElements(Kinetics& kin)
495{
496 m_enamemap.clear();
497 m_nel = 0;
498 for (size_t i = 0; i < kin.nPhases(); i++) {
499 ThermoPhase* p = &kin.thermo(i);
500 // iterate over the elements in this phase
501 for (size_t m = 0; m < p->nElements(); m++) {
502 string ename = p->elementName(m);
503
504 // if no entry is found for this element name, then it is a new
505 // element. In this case, add the name to the list of names,
506 // increment the element count, and add an entry to the
507 // name->(index+1) map.
508 if (m_enamemap.find(ename) == m_enamemap.end()) {
509 m_enamemap[ename] = m_nel + 1;
510 m_elementSymbols.push_back(ename);
511 m_nel++;
512 }
513 }
514 }
515 m_atoms.resize(kin.nTotalSpecies(), m_nel, 0.0);
516 // iterate over the elements
517 for (size_t m = 0; m < m_nel; m++) {
518 size_t k = 0;
519 // iterate over the phases
520 for (size_t ip = 0; ip < kin.nPhases(); ip++) {
521 ThermoPhase* p = &kin.thermo(ip);
522 size_t mlocal = p->elementIndex(m_elementSymbols[m]);
523 for (size_t kp = 0; kp < p->nSpecies(); kp++) {
524 if (mlocal != npos) {
525 m_atoms(k, m) = p->nAtoms(kp, mlocal);
526 }
527 k++;
528 }
529 }
530 }
531}
532
533int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin)
534{
535 m_transfer.clear();
536 m_elementSymbols.clear();
537 findElements(kin);
538 m_ns = kin.nTotalSpecies();
539 m_nr = kin.nReactions();
540
541 // all reactants / products, even ones appearing on both sides of the
542 // reaction
543 vector<vector<size_t>> allProducts(m_nr);
544 vector<vector<size_t>> allReactants(m_nr);
545 for (size_t i = 0; i < m_nr; i++) {
546 for (size_t k = 0; k < m_ns; k++) {
547 for (int n = 0; n < kin.reactantStoichCoeff(k, i); n++) {
548 allReactants[i].push_back(k);
549 }
550 for (int n = 0; n < kin.productStoichCoeff(k, i); n++) {
551 allProducts[i].push_back(k);
552 }
553 }
554 }
555
556 // m_reac and m_prod exclude indices for species that appear on
557 // both sides of the reaction, so that the diagram contains no loops.
558 m_reac.resize(m_nr);
559 m_prod.resize(m_nr);
560 m_ropf.resize(m_nr);
561 m_ropr.resize(m_nr);
562 m_determinate.resize(m_nr);
563 m_x.resize(m_ns); // not currently used ?
564 m_elatoms.resize(m_nel, m_nr);
565
566 for (size_t i = 0; i < m_nr; i++) {
567 // construct the lists of reactant and product indices, not including
568 // molecules that appear on both sides.
569 m_reac[i].clear();
570 m_prod[i].clear();
571 map<size_t, int> net;
572 size_t nr = allReactants[i].size();
573 size_t np = allProducts[i].size();
574 for (size_t ir = 0; ir < nr; ir++) {
575 net[allReactants[i][ir]]--;
576 }
577 for (size_t ip = 0; ip < np; ip++) {
578 net[allProducts[i][ip]]++;
579 }
580
581 for (size_t k = 0; k < m_ns; k++) {
582 if (net[k] < 0) {
583 size_t nmol = -net[k];
584 for (size_t jr = 0; jr < nmol; jr++) {
585 m_reac[i].push_back(k);
586 }
587 } else if (net[k] > 0) {
588 size_t nmol = net[k];
589 for (size_t jp = 0; jp < nmol; jp++) {
590 m_prod[i].push_back(k);
591 }
592 }
593 }
594
595 size_t nrnet = m_reac[i].size();
596
597 // compute number of atoms of each element in each reaction, excluding
598 // molecules that appear on both sides of the reaction. We only need to
599 // compute this for the reactants, since the elements are conserved.
600 for (size_t n = 0; n < nrnet; n++) {
601 size_t k = m_reac[i][n];
602 for (size_t m = 0; m < m_nel; m++) {
603 m_elatoms(m,i) += m_atoms(k,m);
604 }
605 }
606 }
607
608 // build species groups
609 vector<int> comp(m_nel);
610 m_sgroup.resize(m_ns);
611 for (size_t j = 0; j < m_ns; j++) {
612 for (size_t m = 0; m < m_nel; m++) {
613 comp[m] = int(m_atoms(j,m));
614 }
615 m_sgroup[j] = Group(comp);
616 }
617
618 // determine whether or not the reaction is "determinate", meaning that
619 // there is no ambiguity about which reactant is the source for any element
620 // in any product. This is false if more than one reactant contains a given
621 // element, *and* more than one product contains the element. In this case,
622 // additional information is needed to determine the partitioning of the
623 // reactant atoms of that element among the products.
624 for (size_t i = 0; i < m_nr; i++) {
625 size_t nr = m_reac[i].size();
626 size_t np = m_prod[i].size();
627 m_determinate[i] = true;
628 for (size_t m = 0; m < m_nel; m++) {
629 int nar = 0;
630 int nap = 0;
631 for (size_t j = 0; j < nr; j++) {
632 if (m_atoms(m_reac[i][j],m) > 0) {
633 nar++;
634 }
635 }
636 for (size_t j = 0; j < np; j++) {
637 if (m_atoms(m_prod[i][j],m) > 0) {
638 nap++;
639 }
640 }
641 if (nar > 1 && nap > 1) {
642 m_determinate[i] = false;
643 break;
644 }
645 }
646 }
647
648 findGroups(logfile, kin);
649 return 1;
650}
651
652string reactionLabel(size_t i, size_t kr, size_t nr,
653 const vector<size_t>& slist, const Kinetics& s)
654{
655 string label = "";
656 for (size_t j = 0; j < nr; j++) {
657 if (j != kr) {
658 label += " + "+ s.kineticsSpeciesName(slist[j]);
659 }
660 }
661 if (ba::starts_with(s.reaction(i)->type(), "three-body")) {
662 label += " + M ";
663 } else if (ba::starts_with(s.reaction(i)->type(), "falloff")) {
664 label += " (+ M)";
665 }
666 return label;
667}
668
669int ReactionPathBuilder::build(Kinetics& s, const string& element,
670 ostream& output, ReactionPathDiagram& r, bool quiet)
671{
672 map<size_t, int> warn;
673 double threshold = 0.0;
674 size_t m = m_enamemap[element]-1;
675 r.element = element;
676 if (m == npos) {
677 return -1;
678 }
679
680 s.getFwdRatesOfProgress(m_ropf.data());
681 s.getRevRatesOfProgress(m_ropr.data());
682
683 // species explicitly included or excluded
684 vector<string>& in_nodes = r.included();
685 vector<string>& out_nodes = r.excluded();
686
687 vector<int> status(s.nTotalSpecies(), 0);
688 for (size_t ni = 0; ni < in_nodes.size(); ni++) {
689 status[s.kineticsSpeciesIndex(in_nodes[ni])] = 1;
690 }
691 for (size_t ne = 0; ne < out_nodes.size(); ne++) {
692 status[s.kineticsSpeciesIndex(out_nodes[ne])] = -1;
693 }
694
695 for (size_t i = 0; i < m_nr; i++) {
696 double ropf = m_ropf[i];
697 double ropr = m_ropr[i];
698
699 // loop over reactions involving element m
700 if (m_elatoms(m, i) > 0) {
701 size_t nr = m_reac[i].size();
702 size_t np = m_prod[i].size();
703
704 for (size_t kr = 0; kr < nr; kr++) {
705 size_t kkr = m_reac[i][kr];
706 string fwdlabel = reactionLabel(i, kr, nr, m_reac[i], s);
707
708 for (size_t kp = 0; kp < np; kp++) {
709 size_t kkp = m_prod[i][kp];
710 string revlabel = "";
711 for (size_t j = 0; j < np; j++) {
712 if (j != kp) {
713 revlabel += " + "+ s.kineticsSpeciesName(m_prod[i][j]);
714 }
715 }
716 if (ba::starts_with(s.reaction(i)->type(), "three-body")) {
717 revlabel += " + M ";
718 } else if (ba::starts_with(s.reaction(i)->type(), "falloff")) {
719 revlabel += " (+ M)";
720 }
721
722 // calculate the flow only for pairs that are not the same
723 // species, both contain atoms of element m, and both are
724 // allowed to appear in the diagram
725 if ((kkr != kkp) && (m_atoms(kkr,m) > 0
726 && m_atoms(kkp,m) > 0)
727 && status[kkr] >= 0 && status[kkp] >= 0) {
728 // if neither species contains the full number of atoms
729 // of element m in the reaction, then we must consider
730 // the type of reaction to determine which reactant
731 // species was the source of a given m-atom in the
732 // product
733 double f;
734 if ((m_atoms(kkp,m) < m_elatoms(m, i)) &&
735 (m_atoms(kkr,m) < m_elatoms(m, i))) {
736 map<size_t, map<size_t, Group>>& g = m_transfer[i];
737 if (g.empty()) {
738 if (!warn[i] && !quiet) {
739 output << endl;
740 output << "*************** REACTION IGNORED ***************" << endl;
741 output << "Warning: no rule to determine partitioning of " << element
742 << endl << " in reaction " << s.reaction(i)->equation() << "." << endl
743 << "*************** REACTION IGNORED **************" << endl;
744 output << endl;
745 warn[i] = 1;
746 }
747 f = 0.0;
748 } else {
749 if (!g[kr][kp]) {
750 f = 0.0;
751 } else {
752 f = g[kr][kp].nAtoms(m);
753 }
754 }
755 } else {
756 // no ambiguity about where the m-atoms come from or
757 // go to. Either all reactant m atoms end up in one
758 // product, or only one reactant contains all the
759 // m-atoms. In either case, the number of atoms
760 // transferred is given by the same expression.
761 f = m_atoms(kkp,m) * m_atoms(kkr,m) / m_elatoms(m, i);
762 }
763
764 double fwd = ropf*f;
765 double rev = ropr*f;
766 bool force_incl = ((status[kkr] == 1) || (status[kkp] == 1));
767
768 bool fwd_incl = ((fwd > threshold) ||
769 (fwd > 0.0 && force_incl));
770 bool rev_incl = ((rev > threshold) ||
771 (rev > 0.0 && force_incl));
772 if (fwd_incl || rev_incl) {
773 if (!r.hasNode(kkr)) {
774 r.addNode(kkr, s.kineticsSpeciesName(kkr), m_x[kkr]);
775 }
776 if (!r.hasNode(kkp)) {
777 r.addNode(kkp, s.kineticsSpeciesName(kkp), m_x[kkp]);
778 }
779 }
780 if (fwd_incl) {
781 r.linkNodes(kkr, kkp, int(i), fwd, fwdlabel);
782 }
783 if (rev_incl) {
784 r.linkNodes(kkp, kkr, -int(i), rev, revlabel);
785 }
786 }
787 }
788 }
789 }
790 }
791 return 1;
792}
793
794}
Classes for reaction path analysis.
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
virtual void resize(size_t n, size_t m, double v=0.0)
Resize the array, and fill the new entries with 'v'.
Definition Array.cpp:47
Base class for exceptions thrown by Cantera classes.
void resize(size_t n, size_t m, double v=0.0) override
Resize the matrix.
Class Group is an internal class used by class ReactionPath.
Definition Group.h:21
bool valid() const
True if all non-zero atom numbers have the same sign.
Definition Group.h:99
Public interface for kinetics managers.
Definition Kinetics.h:125
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
Definition Kinetics.h:242
shared_ptr< Reaction > reaction(size_t i)
Return the Reaction object for reaction i.
Definition Kinetics.cpp:684
size_t nPhases() const
The number of phases participating in the reaction mechanism.
Definition Kinetics.h:184
size_t nTotalSpecies() const
The total number of species in all phases participating in the kinetics mechanism.
Definition Kinetics.h:254
const SpeciesNode * begin() const
Upstream node.
void addReaction(size_t rxnNumber, double value, const string &label="")
Add a reaction to the path.
void writeLabel(std::ostream &s, double threshold=0.005)
Write the label for a path connecting two species, indicating the percent of the total flow due to ea...
double flow()
The total flow in this path.
Path(SpeciesNode *begin, SpeciesNode *end)
Constructor.
const SpeciesNode * end() const
Downstream node.
size_t elementIndex(const string &name) const
Return the index of element named 'name'.
Definition Phase.cpp:55
size_t nElements() const
Number of elements.
Definition Phase.cpp:30
string elementName(size_t m) const
Name of the element with index m.
Definition Phase.cpp:49
int findGroups(std::ostream &logfile, Kinetics &s)
Analyze a reaction to determine which reactants lead to which products.
map< size_t, map< size_t, map< size_t, Group > > > m_transfer
m_transfer[reaction][reactant number][product number] where "reactant number" means the number of the...
map< size_t, SpeciesNode * > m_nodes
map of species index to SpeciesNode
void exportToDot(std::ostream &s)
Export the reaction path diagram.
set< size_t > m_rxns
Indices of reactions that are included in the diagram.
virtual ~ReactionPathDiagram()
Destructor.
double flow(size_t k1, size_t k2)
The one-way flow from node k1 to node k2.
double netFlow(size_t k1, size_t k2)
The net flow from node k1 to node k2.
Nodes in reaction path graphs.
string name
Label on graph.
bool visible
Visible on graph;.
size_t number
Species number.
void addPath(Path *path)
add a path to or from this node
Base class for a phase with thermodynamic properties.
void warn(const string &warning, const string &method, const string &msg, const Args &... args)
Print a generic warning raised from method.
Definition global.h:251
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