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