Cantera  2.1.2
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 // Copyright 2001 California Institute of Technology
6 
10 #include "cantera/kinetics/Group.h"
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("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  size_t nn = m_label.size();
58  if (nn == 0) {
59  return;
60  }
61  doublereal v;
62  map<string, doublereal>::const_iterator i = m_label.begin();
63  for (; i != m_label.end(); ++i) {
64  v = i->second/m_total;
65  if (nn == 1) {
66  s << i->first << "\\l";
67  } else if (v > threshold) {
68  s << i->first;
69  int percent = int(100*v + 0.5);
70  if (percent < 100) {
71  s << " (" << percent << "%)\\l";
72  } else {
73  s << "\\l";
74  }
75  }
76  }
77 }
78 
79 ReactionPathDiagram::ReactionPathDiagram()
80 {
81  name = "reaction_paths";
82  m_flxmax = 0.0;
83  bold_color = "blue";
84  normal_color = "steelblue";
85  dashed_color = "gray";
86  dot_options = "center=1;";
87  m_font = "Helvetica"; // RXNPATH_FONT;
88  bold_min = 0.2;
89  dashed_max = 0.0;
90  label_min = 0.0;
91  threshold = 0.005;
92  flow_type = NetFlow;
93  scale = -1;
94  x_size = -1.0;
95  y_size = -1.0;
96  arrow_width = -5.0;
97  show_details = false;
98  arrow_hue = 0.6666;
99  title = "";
100  m_local = npos;
101 }
102 
104 {
105  // delete the nodes
106  map<size_t, SpeciesNode*>::const_iterator i = m_nodes.begin();
107  for (; i != m_nodes.end(); ++i) {
108  delete i->second;
109  }
110 
111  // delete the paths
112  size_t nn = nPaths();
113  for (size_t n = 0; n < nn; n++) {
114  delete m_pathlist[n];
115  }
116 }
117 
118 vector_int ReactionPathDiagram::reactions()
119 {
120  size_t i, npaths = nPaths();
121  doublereal flmax = 0.0, flxratio;
122  Path* p;
123  for (i = 0; i < npaths; i++) {
124  p = path(i);
125  if (p->flow() > flmax) {
126  flmax = p->flow();
127  }
128  }
129  m_rxns.clear();
130  for (i = 0; i < npaths; i++) {
131  p = path(i);
132  const Path::rxn_path_map& rxns = p->reactionMap();
133  Path::rxn_path_map::const_iterator m = rxns.begin();
134  for (; m != rxns.end(); ++m) {
135  flxratio = m->second/flmax;
136  if (flxratio > threshold) {
137  m_rxns[m->first] = 1;
138  }
139  }
140  }
141  vector_int r;
142  map<size_t, int>::const_iterator begin = m_rxns.begin();
143  for (; begin != m_rxns.end(); ++begin) {
144  r.push_back(int(begin->first));
145  }
146  return r;
147 }
148 
149 void ReactionPathDiagram::add(ReactionPathDiagram& d)
150 {
151  // doublereal f1, f2;
152  // int nnodes = nNodes();
153  // if (nnodes != d.nNodes()) {
154  // throw CanteraError("ReactionPathDiagram::add",
155  // "number of nodes must be the same");
156  // }
157  size_t np = nPaths();
158  size_t n, k1, k2;
159  Path* p = 0;
160  for (n = 0; n < np; n++) {
161  p = path(n);
162  k1 = p->begin()->number;
163  k2 = p->end()->number;
164  p->setFlow(p->flow() + d.flow(k1,k2));
165  }
166 }
167 
168 void ReactionPathDiagram::findMajorPaths(doublereal athreshold, size_t lda,
169  doublereal* a)
170 {
171  size_t nn = nNodes();
172  size_t n, m, k1, k2;
173  doublereal fl, netmax = 0.0;
174  for (n = 0; n < nn; n++) {
175  for (m = n+1; m < nn; m++) {
176  k1 = m_speciesNumber[n];
177  k2 = m_speciesNumber[m];
178  fl = fabs(netFlow(k1,k2));
179  if (fl > netmax) {
180  netmax = fl;
181  }
182  }
183  }
184  for (n = 0; n < nn; n++) {
185  for (m = n+1; m < nn; m++) {
186  k1 = m_speciesNumber[n];
187  k2 = m_speciesNumber[m];
188  fl = fabs(netFlow(k1,k2));
189  if (fl > athreshold*netmax) {
190  a[lda*k1 + k2] = 1;
191  }
192  }
193  }
194 }
195 
196 void ReactionPathDiagram::writeData(ostream& s)
197 {
198  doublereal f1, f2;
199  size_t nnodes = nNodes();
200  size_t i1, i2, k1, k2;
201  s << title << endl;
202  for (i1 = 0; i1 < nnodes; i1++) {
203  k1 = m_speciesNumber[i1];
204  s << m_nodes[k1]->name << " ";
205  }
206  s << endl;
207  for (i1 = 0; i1 < nnodes; i1++) {
208  k1 = m_speciesNumber[i1];
209  for (i2 = i1+1; i2 < nnodes; i2++) {
210  k2 = m_speciesNumber[i2];
211  f1 = flow(k1, k2);
212  f2 = flow(k2, k1);
213  //if (f1 > 0.001 || f2 > 0.001) {
214  s << m_nodes[k1]->name << " " << m_nodes[k2]->name
215  << " " << f1 << " " << -f2 << endl;
216  //}
217  }
218  }
219 }
220 
222 {
223  doublereal flxratio, flmax = 0.0, lwidth;
224  //s.flags(std::ios_base::showpoint+std::ios_base::fixed);
225  s.precision(3);
226 
227  // a directed graph
228  s << "digraph " << name << " {" << endl;
229 
230  // the graph will be no larger than x_size, y_size
231  if (x_size > 0.0) {
232  if (y_size < 0.0) {
233  y_size = x_size;
234  }
235  s << "size = \""
236  << x_size << ","
237  << y_size << "\";"
238  << endl;
239  }
240 
241  //s << "color = white;" << endl;
242  if (dot_options != "") {
243  s << dot_options << endl;
244  }
245 
246  Path* p;
247 
248  size_t kbegin, kend, i1, i2, k1, k2;
249  doublereal flx;
250 
251  // draw paths representing net flows
252  if (flow_type == NetFlow) {
253 
254  // if no scale was specified, normalize
255  // net flows by the maximum net flow
256  if (scale <= 0.0) {
257  for (i1 = 0; i1 < nNodes(); i1++) {
258  k1 = m_speciesNumber[i1];
259  node(k1)->visible = false;
260  for (i2 = i1+1; i2 < nNodes(); i2++) {
261  k2 = m_speciesNumber[i2];
262  flx = netFlow(k1, k2);
263  if (flx < 0.0) {
264  flx = -flx;
265  }
266  if (flx > flmax) {
267  flmax = flx;
268  }
269  }
270  }
271  } else {
272  flmax = scale;
273  }
274 
275  if (flmax < 1.e-10) {
276  flmax = 1.e-10;
277  }
278 
279  // loop over all unique pairs of nodes
280 
281  for (i1 = 0; i1 < nNodes(); i1++) {
282  k1 = m_speciesNumber[i1];
283  for (i2 = i1+1; i2 < nNodes(); i2++) {
284  k2 = m_speciesNumber[i2];
285  flx = netFlow(k1, k2);
286  if (m_local != npos) {
287  if (k1 != m_local && k2 != m_local) {
288  flx = 0.0;
289  }
290  }
291  if (flx != 0.0) {
292  // set beginning and end of the path based on the
293  // sign of the net flow
294 
295  if (flx > 0.0) {
296  kbegin = k1;
297  kend = k2;
298  flxratio = flx/flmax;
299  } else {
300  kbegin = k2;
301  kend = k1;
302  flxratio = -flx/flmax;
303  }
304 
305  // write out path specification if the net flow
306  // is greater than the threshold
307 
308  if (flxratio >= threshold) {
309  // make nodes visible
310  node(kbegin)->visible = true;
311  node(kend)->visible = true;
312 
313  s << "s" << kbegin << " -> s" << kend;
314  s << "[fontname=\""+m_font+"\", style=\"setlinewidth(";
315 
316  if (arrow_width < 0) {
317  lwidth = 1.0 - 4.0
318  * log10(flxratio/threshold)/log10(threshold) + 1.0;
319  s << lwidth << ")\"";
320  s << ", arrowsize="
321  << std::min(6.0, 0.5*lwidth);
322  } else {
323  s << arrow_width << ")\"";
324  s << ", arrowsize=" << flxratio + 1;
325  }
326 
327  doublereal hue = 0.7;
328  doublereal bright = 0.9;
329  s << ", color=" << "\"" << hue << ", "
330  << flxratio + 0.5
331  << ", " << bright << "\"" << endl;
332 
333  if (flxratio > label_min) {
334  s << ", label=\" " << flxratio;
335  if (show_details) {
336  if (flow(kbegin, kend) > 0.0) {
337  s << "\\l fwd: "
338  << flow(kbegin, kend)/flmax << "\\l";
339  path(kbegin, kend)->writeLabel(s);
340  }
341  if (flow(kend, kbegin) > 0.0) {
342  s << " \\l rev: "
343  << flow(kend,kbegin)/flmax << "\\l";
344  path(kend, kbegin)->writeLabel(s);
345  }
346  }
347  s << "\"";
348  }
349  s << "];" << endl;
350  }
351  }
352  }
353  }
354  }
355 
356  else {
357  for (size_t i = 0; i < nPaths(); i++) {
358  p = path(i);
359  if (p->flow() > flmax) {
360  flmax = p->flow();
361  }
362  }
363 
364  for (size_t i = 0; i < nPaths(); i++) {
365  p = path(i);
366  flxratio = p->flow()/flmax;
367  if (m_local != npos) {
368  if (p->begin()->number != m_local
369  && p->end()->number != m_local) {
370  flxratio = 0.0;
371  }
372  }
373  if (flxratio > threshold) {
374  p->begin()->visible = true;
375  p->end()->visible = true;
376  s << "s" << p->begin()->number
377  << " -> s" << p->end()->number;
378 
379  if (arrow_width < 0) {
380  lwidth = 1.0 - 4.0 * log10(flxratio/threshold)/log10(threshold)
381  + 1.0;
382  s << "[fontname=\""+m_font+"\", style=\"setlinewidth("
383  //<< 1.0 - arrow_width*flxratio
384  << lwidth
385  << ")\"";
386  s << ", arrowsize="
387  << std::min(6.0, 0.5*lwidth); // 1 - arrow_width*flxratio;
388  } else {
389  s << ", style=\"setlinewidth("
390  << arrow_width << ")\"";
391  s << ", arrowsize=" << flxratio + 1;
392  }
393  doublereal hue = 0.7; //2.0/(1.0 + pow(log10(flxratio),2)) ;
394  doublereal bright = 0.9;
395  s << ", color=" << "\"" << hue << ", " << flxratio + 0.5
396  << ", " << bright << "\"" << endl;
397 
398  if (flxratio > label_min) {
399  s << ", label = \" " << flxratio;
400  if (show_details) {
401  s << "\\l";
402  p->writeLabel(s);
403  }
404  s << "\"";
405  }
406  s << "];" << endl;
407  }
408  }
409  }
410  s.precision(2);
411  map<size_t, SpeciesNode*>::const_iterator b = m_nodes.begin();
412  for (; b != m_nodes.end(); ++b) {
413  if (b->second->visible) {
414  s << "s" << b->first << " [ fontname=\""+m_font+"\", label=\"" << b->second->name
415  //<< " \\n " << b->second->value
416  << "\"];" << endl;
417  }
418  }
419  s << " label = " << "\"" << "Scale = "
420  << flmax << "\\l " << title << "\";" << endl; //created with Cantera (www.cantera.org)\\l\";"
421  s << " fontname = \""+m_font+"\";" << endl << "}" << endl;
422 }
423 
424 
425 void ReactionPathDiagram::addNode(size_t k, const string& nm, doublereal x)
426 {
427  if (!m_nodes[k]) {
428  m_nodes[k] = new SpeciesNode;
429  m_nodes[k]->number = k;
430  m_nodes[k]->name = nm;
431  m_nodes[k]->value = x;
432  m_speciesNumber.push_back(k);
433  }
434 }
435 
436 void ReactionPathDiagram::linkNodes(size_t k1, size_t k2, size_t rxn,
437  doublereal value, string legend)
438 {
439  SpeciesNode* begin = m_nodes[k1];
440  SpeciesNode* end = m_nodes[k2];
441  Path* ff = m_paths[k1][k2];
442  if (!ff) {
443  ff= new Path(begin, end);
444  m_paths[k1][k2] = ff;
445  m_pathlist.push_back(ff);
446  }
447  ff->addReaction(rxn, value, legend);
448  m_rxns[rxn] = 1;
449  if (ff->flow() > m_flxmax) {
450  m_flxmax = ff->flow();
451  }
452 }
453 
454 std::vector<size_t> ReactionPathDiagram::species()
455 {
456  return m_speciesNumber;
457 }
458 
459 int ReactionPathBuilder::findGroups(ostream& logfile, Kinetics& s)
460 {
461  m_groups.resize(m_nr);
462 
463  for (size_t i = 0; i < m_nr; i++) { // loop over reactions
464  logfile << endl << "Reaction " << i+1 << ": "
465  << s.reactionString(i);
466 
467  size_t nrnet = m_reac[i].size();
468  size_t npnet = m_prod[i].size();
469  const std::vector<size_t>& r = s.reactants(i);
470  const std::vector<size_t>& p = s.products(i);
471 
472  size_t nr = r.size();
473  size_t np = p.size();
474 
475  Group b0, b1, bb;
476 
477  vector<string>& e = m_elementSymbols;
478 
479  const vector<grouplist_t>& rgroups = s.reactantGroups(i);
480  const vector<grouplist_t>& pgroups = s.productGroups(i);
481 
482  if (m_determinate[i]) {
483  logfile << " ... OK." << endl;
484  }
485 
486  else if (rgroups.size() > 0) {
487  logfile << " ... specified groups." << endl;
488  size_t nrg = rgroups.size();
489  size_t npg = pgroups.size();
490  size_t kr, kp, ngrpr, ngrpp;
491  Group gr, gp;
492 
493  if (nrg != nr || npg != np) {
494  return -1;
495  }
496 
497  // loop over reactants
498  for (size_t igr = 0; igr < nrg; igr++) {
499  kr = r[igr];
500  ngrpr = rgroups[igr].size();
501 
502  // loop over products
503  for (size_t igp = 0; igp < npg; igp++) {
504  kp = p[igp];
505  ngrpp = pgroups[igp].size();
506 
507  // loop over pairs of reactant and product groups
508  for (size_t kgr = 0; kgr < ngrpr; kgr++) {
509  gr = Group(rgroups[igr][kgr]);
510  for (size_t kgp = 0; kgp < ngrpp; kgp++) {
511  gp = Group(pgroups[igp][kgp]);
512  if (gr == gp) {
513  m_transfer[i][kr][kp] = gr;
514  }
515  }
516  }
517  }
518  }
519  }
520 
521  else if (nrnet == 2 && npnet == 2) {
522  // indices for the two reactants
523  size_t kr0 = m_reac[i][0];
524  size_t kr1 = m_reac[i][1];
525 
526  // indices for the two products
527  size_t kp0 = m_prod[i][0];
528  size_t kp1 = m_prod[i][1];
529 
530  // references to the Group objects representing the
531  // reactants
532  const Group& r0 = m_sgroup[kr0];
533  const Group& r1 = m_sgroup[kr1];
534  const Group& p0 = m_sgroup[kp0];
535  const Group& p1 = m_sgroup[kp1];
536 
537  const Group* group_a0=0, *group_b0=0, *group_c0=0,
538  *group_a1=0, *group_b1=0, *group_c1=0;
539  b0 = p0 - r0;
540  b1 = p1 - r0;
541  if (b0.valid() && b1.valid()) {
542  logfile << " ... ambiguous." << endl;
543  } else if (!b0.valid() && !b1.valid()) {
544  logfile << " ... cannot express as A + BC = AB + C" << endl;
545  } else {
546  logfile << endl;
547  }
548 
549  if (b0.valid()) {
550  if (b0.sign() > 0) {
551  group_a0 = &r0;
552  group_b0 = &b0;
553  group_c0 = &p1;
554  m_transfer[i][kr0][kp0] = r0;
555  m_transfer[i][kr1][kp0] = b0;
556  m_transfer[i][kr1][kp1] = p1;
557  } else {
558  group_a0 = &r1;
559  group_c0 = &p0;
560  b0 *= -1;
561  group_b0 = &b0;
562  m_transfer[i][kr1][kp1] = r1;
563  m_transfer[i][kr0][kp1] = b0;
564  m_transfer[i][kr0][kp0] = p0;
565  }
566  logfile << " ";
567  group_a0->fmt(logfile,e);
568  logfile << " + ";
569  group_b0->fmt(logfile,e);
570  group_c0->fmt(logfile,e);
571  logfile << " = ";
572  group_a0->fmt(logfile,e);
573  group_b0->fmt(logfile,e);
574  logfile << " + ";
575  group_c0->fmt(logfile,e);
576  if (b1.valid()) {
577  logfile << " [<= default] " << endl;
578  } else {
579  logfile << endl;
580  }
581  }
582 
583 
584  if (b1.valid()) {
585  if (b1.sign() > 0) {
586  group_a1 = &r0;
587  group_b1 = &b1;
588  group_c1 = &p0;
589  if (!b0.valid()) {
590  m_transfer[i][kr0][kp1] = r0;
591  m_transfer[i][kr1][kp1] = b0;
592  m_transfer[i][kr1][kp0] = p0;
593  }
594  } else {
595  group_a1 = &r1;
596  group_c1 = &p1;
597  b1 *= -1;
598  group_b1 = &b1;
599  if (!b0.valid()) {
600  m_transfer[i][kr1][kp0] = r1;
601  m_transfer[i][kr0][kp0] = b0;
602  m_transfer[i][kr0][kp1] = p1;
603  }
604  }
605  logfile << " ";
606  group_a1->fmt(logfile,e);
607  logfile << " + ";
608  group_b1->fmt(logfile,e);
609  group_c1->fmt(logfile,e);
610  logfile << " = ";
611  group_a1->fmt(logfile,e);
612  group_b1->fmt(logfile,e);
613  logfile << " + ";
614  group_c1->fmt(logfile,e);
615  logfile << endl;
616  }
617  } else {
618  logfile << "... cannot parse. [ignored]" << endl;
619  }
620  }
621  return 1;
622 }
623 
624 void ReactionPathBuilder::writeGroup(ostream& out, const Group& g)
625 {
626  g.fmt(out, m_elementSymbols);
627 }
628 
629 void ReactionPathBuilder::findElements(Kinetics& kin)
630 {
631 
632  string ename;
633  m_enamemap.clear();
634  m_nel = 0;
635  size_t np = kin.nPhases();
636  ThermoPhase* p;
637  for (size_t i = 0; i < np; i++) {
638  p = &kin.thermo(i);
639  // iterate over the elements in this phase
640  size_t nel = p->nElements();
641  for (size_t m = 0; m < nel; m++) {
642  ename = p->elementName(m);
643 
644  // if no entry is found for this element name, then
645  // it is a new element. In this case, add the name
646  // to the list of names, increment the element count,
647  // and add an entry to the name->(index+1) map.
648  if (m_enamemap.find(ename) == m_enamemap.end()) {
649  m_enamemap[ename] = m_nel + 1;
650  m_elementSymbols.push_back(ename);
651  m_nel++;
652  }
653  }
654  }
655  m_atoms.resize(kin.nTotalSpecies(), m_nel, 0.0);
656  string sym;
657  // iterate over the elements
658  for (size_t m = 0; m < m_nel; m++) {
659  sym = m_elementSymbols[m];
660  size_t k = 0;
661  // iterate over the phases
662  for (size_t ip = 0; ip < np; ip++) {
663  ThermoPhase* p = &kin.thermo(ip);
664  size_t nsp = p->nSpecies();
665  size_t mlocal = p->elementIndex(sym);
666  for (size_t kp = 0; kp < nsp; kp++) {
667  if (mlocal != npos) {
668  m_atoms(k, m) = p->nAtoms(kp, mlocal);
669  }
670  k++;
671  }
672  }
673  }
674 }
675 
676 int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin)
677 {
678  //m_warn.clear();
679  m_transfer.clear();
680 
681  //const Kinetics::thermo_t& ph = kin.thermo();
682 
683  m_elementSymbols.clear();
684  findElements(kin);
685  //m_nel = ph.nElements();
686  m_ns = kin.nTotalSpecies(); //ph.nSpecies();
687  m_nr = kin.nReactions();
688 
689  //for (m = 0; m < m_nel; m++) {
690  // m_elementSymbols.push_back(ph.elementName(m));
691  //}
692 
693  // all reactants / products, even ones appearing on both sides
694  // of the reaction
695  // mod 8/18/01 dgg
696  vector<vector<size_t> > allProducts;
697  vector<vector<size_t> > allReactants;
698  for (size_t i = 0; i < m_nr; i++) {
699  allReactants.push_back(kin.reactants(i));
700  allProducts.push_back(kin.products(i));
701  }
702 
703  // m_reac and m_prod exclude indices for species that appear on
704  // both sides of the reaction, so that the diagram contains no loops.
705 
706  m_reac.resize(m_nr);
707  m_prod.resize(m_nr);
708 
709  m_ropf.resize(m_nr);
710  m_ropr.resize(m_nr);
711  m_determinate.resize(m_nr);
712 
713  m_x.resize(m_ns); // not currently used ?
714  m_elatoms.resize(m_nel, m_nr);
715 
716  size_t nr, np, n, k;
717  size_t nmol;
718  map<size_t, int> net;
719 
720  for (size_t i = 0; i < m_nr; i++) {
721 
722  // construct the lists of reactant and product indices, not
723  // including molecules that appear on both sides.
724 
725  m_reac[i].clear();
726  m_prod[i].clear();
727  net.clear();
728  nr = allReactants[i].size();
729  np = allProducts[i].size();
730  for (size_t ir = 0; ir < nr; ir++) {
731  net[allReactants[i][ir]]--;
732  }
733  for (size_t ip = 0; ip < np; ip++) {
734  net[allProducts[i][ip]]++;
735  }
736 
737  for (k = 0; k < m_ns; k++) {
738  if (net[k] < 0) {
739  nmol = -net[k];
740  for (size_t jr = 0; jr < nmol; jr++) {
741  m_reac[i].push_back(k);
742  }
743  } else if (net[k] > 0) {
744  nmol = net[k];
745  for (size_t jp = 0; jp < nmol; jp++) {
746  m_prod[i].push_back(k);
747  }
748  }
749  }
750 
751  size_t nrnet = m_reac[i].size();
752  // int npnet = m_prod[i].size();
753 
754  // compute number of atoms of each element in each reaction,
755  // excluding molecules that appear on both sides of the
756  // reaction. We only need to compute this for the reactants,
757  // since the elements are conserved.
758 
759  for (n = 0; n < nrnet; n++) {
760  k = m_reac[i][n];
761  for (size_t m = 0; m < m_nel; m++) {
762  m_elatoms(m,i) += m_atoms(k,m); //ph.nAtoms(k,m);
763  }
764  }
765  }
766 
767  // build species groups
768  vector_int comp(m_nel);
769  m_sgroup.resize(m_ns);
770  for (size_t j = 0; j < m_ns; j++) {
771  for (size_t m = 0; m < m_nel; m++) {
772  comp[m] = int(m_atoms(j,m)); //ph.nAtoms(j,m));
773  }
774  m_sgroup[j] = Group(comp);
775  }
776 
777 
778  // determine whether or not the reaction is "determinate", meaning
779  // that there is no ambiguity about which reactant is the source for
780  // any element in any product. This is false if more than one
781  // reactant contains a given element, *and* more than one product
782  // contains the element. In this case, additional information is
783  // needed to determine the partitioning of the reactant atoms of
784  // that element among the products.
785 
786  int nar, nap;
787  for (size_t i = 0; i < m_nr; i++) {
788  nr = m_reac[i].size();
789  np = m_prod[i].size();
790  m_determinate[i] = true;
791  for (size_t m = 0; m < m_nel; m++) {
792  nar = 0;
793  nap = 0;
794  for (size_t j = 0; j < nr; j++) {
795  // if (ph.nAtoms(m_reac[i][j],m) > 0) nar++;
796  if (m_atoms(m_reac[i][j],m) > 0) {
797  nar++;
798  }
799  }
800  for (size_t j = 0; j < np; j++) {
801  if (m_atoms(m_prod[i][j],m) > 0) {
802  nap++;
803  }
804  }
805  if (nar > 1 && nap > 1) {
806  m_determinate[i] = false;
807  break;
808  }
809  }
810  }
811 
812  findGroups(logfile, kin);
813  return 1;
814 }
815 
816 string reactionLabel(size_t i, size_t kr, size_t nr,
817  const std::vector<size_t>& slist, const Kinetics& s)
818 {
819 
820  //int np = s.nPhases();
821  string label = "";
822  for (size_t l = 0; l < nr; l++) {
823  if (l != kr) {
824  label += " + "+ s.kineticsSpeciesName(slist[l]);
825  }
826  }
827  if (s.reactionType(i) == THREE_BODY_RXN) {
828  label += " + M ";
829  } else if (s.reactionType(i) == FALLOFF_RXN) {
830  label += " (+ M)";
831  }
832  return label;
833 }
834 
835 int ReactionPathBuilder::build(Kinetics& s, const string& element,
836  ostream& output, ReactionPathDiagram& r, bool quiet)
837 {
838  doublereal f, ropf, ropr, fwd, rev;
839  string fwdlabel, revlabel;
840  map<size_t, int> warn;
841 
842  doublereal threshold = 0.0;
843  bool fwd_incl, rev_incl, force_incl;
844 
845  // const Kinetics::thermo_t& ph = s.thermo();
846  size_t m = m_enamemap[element]-1; //ph.elementIndex(element);
847 
848  r.element = element;
849  if (m == npos) {
850  return -1;
851  }
852 
853  s.getFwdRatesOfProgress(DATA_PTR(m_ropf));
854  s.getRevRatesOfProgress(DATA_PTR(m_ropr));
855 
856  //ph.getMoleFractions(m_x.begin());
857 
858  //doublereal sum = 0.0;
859  //for (k = 0; k < kk; k++) {
860  // sum += m_x[k] * ph.nAtoms(k,m);
861  //}
862  //sum *= ph.molarDensity();
863 
864  // species explicitly included or excluded
865  vector<string>& in_nodes = r.included();
866  vector<string>& out_nodes = r.excluded();
867 
868  vector_int status;
869  status.resize(s.nTotalSpecies(), 0);
870  for (size_t ni = 0; ni < in_nodes.size(); ni++) {
871  status[s.kineticsSpeciesIndex(in_nodes[ni])] = 1;
872  }
873  for (size_t ne = 0; ne < out_nodes.size(); ne++) {
874  status[s.kineticsSpeciesIndex(out_nodes[ne])] = -1;
875  }
876 
877  for (size_t i = 0; i < m_nr; i++) {
878  ropf = m_ropf[i];
879  ropr = m_ropr[i];
880 
881  // loop over reactions involving element m
882  if (m_elatoms(m, i) > 0) {
883  size_t nr = m_reac[i].size();
884  size_t np = m_prod[i].size();
885 
886  for (size_t kr = 0; kr < nr; kr++) {
887  size_t kkr = m_reac[i][kr];
888 
889  fwdlabel = reactionLabel(i, kr, nr, m_reac[i], s);
890 
891  for (size_t kp = 0; kp < np; kp++) {
892  size_t kkp = m_prod[i][kp];
893  revlabel = "";
894  for (size_t l = 0; l < np; l++) {
895  if (l != kp) {
896  revlabel += " + "+ s.kineticsSpeciesName(m_prod[i][l]);
897  }
898  }
899  if (s.reactionType(i) == THREE_BODY_RXN) {
900  revlabel += " + M ";
901  } else if (s.reactionType(i) == FALLOFF_RXN) {
902  revlabel += " (+ M)";
903  }
904 
905 
906  // calculate the flow only for pairs that are
907  // not the same species, both contain atoms of
908  // element m, and both are allowed to appear in
909  // the diagram
910 
911  if ((kkr != kkp) && (m_atoms(kkr,m) > 0
912  && m_atoms(kkp,m) > 0)
913  && status[kkr] >= 0 && status[kkp] >= 0) {
914 
915  // if neither species contains the full
916  // number of atoms of element m in the
917  // reaction, then we must consider the
918  // type of reaction to determine which
919  // reactant species was the source of a
920  // given m-atom in the product
921 
922  if ((m_atoms(kkp,m) < m_elatoms(m, i)) &&
923  (m_atoms(kkr,m) < m_elatoms(m, i))) {
924  map<size_t, map<size_t, Group> >& g = m_transfer[i];
925  if (g.empty()) {
926  if (!warn[i]) {
927  if (!quiet) {
928  output << endl;
929  output << "*************** REACTION IGNORED ***************" << endl;
930  output << "Warning: no rule to determine partitioning of " << element
931  << endl << " in reaction " << s.reactionString(i) << "." << endl
932  << "*************** REACTION IGNORED **************" << endl;
933  output << endl;
934  warn[i] = 1;
935  }
936  }
937  f = 0.0;
938  } else {
939  if (!g[kkr][kkp]) {
940  f = 0.0;
941  } else {
942  f = g[kkr][kkp].nAtoms(m);
943  }
944  }
945  }
946 
947  // no ambiguity about where the m-atoms come
948  // from or go to. Either all reactant m atoms
949  // end up in one product, or only one reactant
950  // contains all the m-atoms. In either case,
951  // the number of atoms transferred is given by
952  // the same expression.
953 
954  else {
955  f = m_atoms(kkp,m) * m_atoms(kkr,m) / m_elatoms(m, i);
956  }
957 
958  fwd = ropf*f;
959  rev = ropr*f;
960  force_incl = ((status[kkr] == 1) || (status[kkp] == 1));
961 
962  fwd_incl = ((fwd > threshold) ||
963  (fwd > 0.0 && force_incl));
964  rev_incl = ((rev > threshold) ||
965  (rev > 0.0 && force_incl));
966  if (fwd_incl || rev_incl) {
967  if (!r.hasNode(kkr)) {
968  r.addNode(kkr, s.kineticsSpeciesName(kkr), m_x[kkr]);
969  }
970  if (!r.hasNode(kkp)) {
971  r.addNode(kkp, s.kineticsSpeciesName(kkp), m_x[kkp]);
972  }
973  }
974  if (fwd_incl) {
975  r.linkNodes(kkr, kkp, int(i), fwd, fwdlabel);
976  }
977  if (rev_incl) {
978  r.linkNodes(kkp, kkr, -int(i), rev, revlabel);
979  }
980  }
981  }
982  }
983  }
984  }
985  return 1;
986 }
987 
988 }
bool visible
Visible on graph;.
Definition: ReactionPath.h:42
void exportToDot(std::ostream &s)
Export the reaction path diagram.
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:173
Nodes in reaction path graphs.
Definition: ReactionPath.h:28
const int FALLOFF_RXN
The general form for an association or dissociation reaction, with a pressure-dependent rate...
Definition: reaction_defs.h:38
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:167
doublereal flow(size_t k1, size_t k2)
The one-way flow from node k1 to node k2.
Definition: ReactionPath.h:187
virtual ~ReactionPathDiagram()
Destructor.
This file defines some constants used to specify reaction types.
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:68
Base class for kinetics managers and also contains the kineticsmgr module documentation (see Kinetics...
const int THREE_BODY_RXN
A reaction that requires a third-body collision partner.
Definition: reaction_defs.h:32
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
Classes for reaction path analysis.
doublereal netFlow(size_t k1, size_t k2)
The net flow from node k1 to node k2.
Definition: ReactionPath.h:182