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