Cantera  2.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ExtraGlobalRxn.cpp
1 /**
2  * @file example2.cpp
3  *
4  */
5 /*
6  * $Id: ExtraGlobalRxn.cpp 571 2013-03-26 16:44:21Z hkmoffa $
7  *
8  */
9 
10 /*
11  * Copyright (2005) Sandia Corporation. Under the terms of
12  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
13  * U.S. Government retains certain rights in this software.
14  */
15 
16 // Example 2
17 //
18 // Read a mechanism, and print to the standard output stream a
19 // well-formatted Chemkin ELEMENT section.
20 //
21 
22 #include "cantera/kinetics/ExtraGlobalRxn.h"
23 
24 
25 
27 
28 // Kinetics includes
29 #include "cantera/kinetics.h"
33 
34 #include <iostream>
35 #include <new>
36 #include <string>
37 #include <vector>
38 #include <typeinfo>
39 
40 
41 
42 using namespace std;
43 using namespace Cantera;
44 
45 namespace Cantera {
46 
47 //============================================================================================================
48 static void erase_vd(std::vector<doublereal>& m_vec, int index)
49 {
50  std::vector<double>::iterator ipos;
51  ipos = m_vec.begin();
52  ipos += index;
53  m_vec.erase(ipos);
54 }
55 //============================================================================================================
56 static void erase_vi(std::vector<int>& m_vec, int index)
57 {
58  std::vector<int>::iterator ipos;
59  ipos = m_vec.begin();
60  ipos += index;
61  m_vec.erase(ipos);
62 }
63 //============================================================================================================
64 //! add the species into the list of products or reactants
65 /*!
66  * Note this function gets called for both the product and reactant sides. However, it's only
67  * called for one side or another.
68  *
69  * @param kkinspec kinetic species index of the product
70  * @param
71  */
72 static void addV(int kkinspec, double ps, std::vector<int>& m_Products,
73  std::vector<doublereal>& m_ProductStoich)
74 {
75  int nsize = static_cast<int>(m_Products.size());
76  for (int i = 0; i < nsize; i++) {
77  if (m_Products[i] == kkinspec) {
78  m_ProductStoich[i] += ps;
79  return;
80  }
81  }
82  m_Products.push_back(kkinspec);
83  m_ProductStoich.push_back(ps);
84 }
85 //============================================================================================================
86 ExtraGlobalRxn::ExtraGlobalRxn(Kinetics* k_ptr) :
87  m_ThisIsASurfaceRxn(false),
88  m_kinetics(k_ptr),
89  m_InterfaceKinetics(0),
90  m_nKinSpecies(0),
91  m_nReactants(0),
92  m_nProducts(0),
93  m_nNetSpecies(0),
94  m_nRxns(0),
95  m_SpecialSpecies(-1),
96  m_SpecialSpeciesProduct(true),
97  iphaseKin(0),
98  m_ok(false),
99  m_reversible(true)
100 {
101  warn_deprecated("class ExtraGlobalRxn",
102  "Unfinished implementation to be removed after Cantera 2.2.");
103  m_InterfaceKinetics = dynamic_cast<InterfaceKinetics*>(k_ptr);
104  if (m_InterfaceKinetics) {
105  m_ThisIsASurfaceRxn = true;
106  }
107  m_nRxns = static_cast<int>(m_kinetics->nReactions());
108  m_ElemRxnVector.resize(m_nRxns,0.0);
109  m_nKinSpecies = static_cast<int>(m_kinetics->nTotalSpecies());
110 }
111 //============================================================================================================
112 void ExtraGlobalRxn::setupElemRxnVector(double* RxnVector,
113  int specialSpecies)
114 {
115  int i;
116  int kkinspec;
117  for (size_t i = 0; i < (size_t) m_nRxns; i++) {
118  m_ElemRxnVector[i] = RxnVector[i];
119  }
120  for (size_t i = 0; i < (size_t) m_nRxns; i++) {
121  if (RxnVector[i] > 0.0) {
122  for (kkinspec = 0; kkinspec < m_nKinSpecies; kkinspec++) {
123  double ps = m_kinetics->productStoichCoeff(kkinspec, i);
124  if (ps > 0.0) {
125  addV(kkinspec, RxnVector[i]* ps, m_Products, m_ProductStoich);
126  addV(kkinspec, RxnVector[i]* ps, m_NetSpecies, m_netStoich);
127  }
128  double rs = m_kinetics->reactantStoichCoeff(kkinspec, i);
129  if (rs > 0.0) {
130  addV(kkinspec, RxnVector[i] * rs, m_Reactants, m_ReactantStoich);
131  addV(kkinspec, -RxnVector[i] * rs, m_NetSpecies, m_netStoich);
132  }
133  }
134  } else if (RxnVector[i] < 0.0) {
135  for (kkinspec = 0; kkinspec < m_nKinSpecies; kkinspec++) {
136  double ps = m_kinetics->productStoichCoeff(kkinspec, i);
137  if (ps > 0.0) {
138  addV(kkinspec,- RxnVector[i]* ps, m_Reactants, m_ReactantStoich);
139  addV(kkinspec, RxnVector[i]* ps, m_NetSpecies, m_netStoich);
140  }
141  double rs = m_kinetics->reactantStoichCoeff(kkinspec, i);
142  if (rs > 0.0) {
143  addV(kkinspec, -RxnVector[i] * rs, m_Products, m_ProductStoich);
144  addV(kkinspec, -RxnVector[i] * rs, m_NetSpecies, m_netStoich);
145  }
146  }
147  }
148  }
149 Recheck:
150  for (i = 0; i < static_cast<int>(m_Products.size()); i++) {
151  if (m_ProductStoich[i] == 0.0) {
152  erase_vi(m_Products, i);
153  erase_vd(m_ProductStoich, i);
154  goto Recheck ;
155  }
156  }
157  for (i = 0; i < static_cast<int>(m_Reactants.size()); i++) {
158  if (m_ReactantStoich[i] == 0.0) {
159  erase_vi(m_Reactants, i);
160  erase_vd(m_ReactantStoich, i);
161  goto Recheck ;
162  }
163  }
164  for (i = 0; i < static_cast<int>(m_NetSpecies.size()); i++) {
165  if (m_netStoich[i] == 0.0) {
166  erase_vi(m_NetSpecies, i);
167  erase_vd(m_netStoich, i);
168  goto Recheck ;
169  }
170  }
171 
172  for (i = 0; i < static_cast<int>(m_Products.size()); i++) {
173  int ik = m_Products[i];
174  for (int j = 0; j < static_cast<int>(m_Reactants.size()); j++) {
175  int jk = m_Reactants[j];
176  if (ik == jk) {
177  if (m_ProductStoich[i] == m_ReactantStoich[j]) {
178  erase_vi(m_Products, i);
179  erase_vd(m_ProductStoich, i);
180  erase_vi(m_Reactants, j);
181  erase_vd(m_ReactantStoich, j);
182  } else if (m_ProductStoich[i] > m_ReactantStoich[j]) {
183  m_ProductStoich[i] -= m_ReactantStoich[j];
184  erase_vi(m_Reactants, j);
185  erase_vd(m_ReactantStoich, j);
186  } else {
187  m_ReactantStoich[j] -= m_ProductStoich[i];
188  erase_vi(m_Products, i);
189  erase_vd(m_ProductStoich, i);
190  }
191  // We just screwed up the indexing -> restart.
192  goto Recheck ;
193  }
194 
195  }
196  }
197  m_nProducts = static_cast<int>(m_Products.size());
198  m_nReactants = static_cast<int>(m_Reactants.size());
199  m_nNetSpecies = static_cast<int>(m_NetSpecies.size());
200 
201  /*
202  * Section to assign the special species
203  */
204  m_SpecialSpecies = specialSpecies;
205  if (specialSpecies == -1) {
206  m_SpecialSpecies = m_Products[0];
207  }
208  bool ifound = false;
209  for (i = 0; i < (int) m_NetSpecies.size(); i++) {
210  int ik = m_NetSpecies[i];
211  if (ik == m_SpecialSpecies) {
212  if (m_netStoich[i] > 0.0) {
213  m_SpecialSpeciesProduct = true;
214  } else {
215  m_SpecialSpeciesProduct = false;
216  }
217  m_SS_index = i;
218  ifound = true;
219  break;
220  }
221  }
222  if (!ifound) {
223  throw CanteraError(":setupElemRxnVector",
224  "Special species not a reactant or product: "
225  + int2str(m_SpecialSpecies));
226  }
227 
228  m_ok = true;
229 }
230 //============================================================================================================
231 std::string ExtraGlobalRxn::reactionString()
232 {
233  string rs;
234  int k, istoich;
235  for (k = 0; k < m_nReactants; k++) {
236  int kkinspecies = m_Reactants[k];
237  double stoich = m_ReactantStoich[k];
238  if (stoich != 1.0) {
239  istoich = (int) stoich;
240  if (fabs((double)istoich - stoich) < 0.00001) {
241  rs += int2str(istoich) + " ";
242  } else {
243  rs += fp2str(stoich) + " ";
244  }
245  }
246  string sName = m_kinetics->kineticsSpeciesName(kkinspecies);
247  rs += sName;
248  if (k < (m_nReactants-1)) {
249  rs += " + ";
250  }
251  }
252  rs += " = ";
253  for (k = 0; k < m_nProducts; k++) {
254  int kkinspecies = m_Products[k];
255  double stoich = m_ProductStoich[k];
256  if (stoich != 1.0) {
257  istoich = (int) stoich;
258  if (fabs((double)istoich - stoich) < 0.00001) {
259  rs += int2str(istoich) + " ";
260  } else {
261  rs += fp2str(stoich) + " ";
262  }
263  }
264  string sName = m_kinetics->kineticsSpeciesName(kkinspecies);
265  rs += sName;
266  if (k < (m_nProducts-1)) {
267  rs += " + ";
268  }
269  }
270  return rs;
271 }
272 //============================================================================================================
273 
274 std::vector<int>& ExtraGlobalRxn::reactants()
275 {
276  return m_Reactants;
277 }
278 //============================================================================================================
279 
280 std::vector<int>& ExtraGlobalRxn::products()
281 {
282  return m_Products;
283 }
284 
285 //============================================================================================================
286 bool ExtraGlobalRxn::isReversible()
287 {
288  return m_reversible;
289 }
290 //============================================================================================================
291 
292 double ExtraGlobalRxn::reactantStoichCoeff(int kKin)
293 {
294  for (int k = 0; k < m_nReactants; k++) {
295  int kkinspec = m_Reactants[k];
296  if (kkinspec == kKin) {
297  return m_ReactantStoich[k];
298  }
299  }
300  return 0.0;
301 }
302 //============================================================================================================
303 double ExtraGlobalRxn::productStoichCoeff(int kKin)
304 {
305  for (int k = 0; k < m_nProducts; k++) {
306  int kkinspec = m_Products[k];
307  if (kkinspec == kKin) {
308  return m_ProductStoich[k];
309  }
310  }
311  return 0.0;
312 }
313 //============================================================================================================
314 
315 double ExtraGlobalRxn::deltaSpecValue(double* speciesVectorProperty)
316 {
317  int k;
318  double val = 0;
319  for (k = 0; k < m_nNetSpecies; k++) {
320  int kkinspec = m_NetSpecies[k];
321  val += speciesVectorProperty[kkinspec] * m_netStoich[k];
322  }
323  return val;
324 }
325 //============================================================================================================
326 
327 double ExtraGlobalRxn::deltaRxnVecValue(double* rxnVectorProperty)
328 {
329  double val = 0;
330  for (int i = 0; i < m_nRxns; i++) {
331  val += m_ElemRxnVector[i] * rxnVectorProperty[i];
332  }
333  return val;
334 }
335 //============================================================================================================
336 double ExtraGlobalRxn::ROPValue(double* ROPElemKinVector)
337 {
338  double val = 0.0;
339  for (int i = 0; i < m_nRxns; i++) {
340  double kstoich = m_kinetics->productStoichCoeff(m_SpecialSpecies, i) - m_kinetics->reactantStoichCoeff(m_SpecialSpecies, i);
341  if (m_ElemRxnVector[i] > 0.0) {
342  val += ROPElemKinVector[i] * kstoich * m_ElemRxnVector[i];
343  } else {
344  val -= ROPElemKinVector[i] * kstoich * m_ElemRxnVector[i];
345  }
346  }
347  if (!m_SpecialSpeciesProduct) {
348  val = -val;
349  }
350  return val;
351 }
352 //============================================================================================================
353 double ExtraGlobalRxn::FwdROPValue(double* FwdROPElemKinVector,
354  double* RevROPElemKinVector)
355 {
356  double val = 0.0;
357  for (int i = 0; i < m_nRxns; i++) {
358  double kstoich = m_kinetics->productStoichCoeff(m_SpecialSpecies, i) - m_kinetics->reactantStoichCoeff(m_SpecialSpecies, i);
359  if (m_ElemRxnVector[i] > 0.0) {
360  val += FwdROPElemKinVector[i] * kstoich * m_ElemRxnVector[i];
361  }
362  if (m_ElemRxnVector[i] < 0.0) {
363  val += RevROPElemKinVector[i] * kstoich * m_ElemRxnVector[i];
364  }
365  }
366  if (!m_SpecialSpeciesProduct) {
367  val = -val;
368  }
369  return val;
370 }
371 //============================================================================================================
372 double ExtraGlobalRxn::RevROPValue(double* FwdROPElemKinVector,
373  double* RevROPElemKinVector)
374 {
375  double val = 0.0;
376  for (int i = 0; i < m_nRxns; i++) {
377  double kstoich = m_kinetics->productStoichCoeff(m_SpecialSpecies, i)- m_kinetics->reactantStoichCoeff(m_SpecialSpecies, i);
378  if (m_ElemRxnVector[i] > 0.0) {
379  val += RevROPElemKinVector[i] * kstoich * m_ElemRxnVector[i];
380  }
381  if (m_ElemRxnVector[i] < 0.0) {
382  val += FwdROPElemKinVector[i] * kstoich * m_ElemRxnVector[i];
383  }
384  }
385  if (!m_SpecialSpeciesProduct) {
386  val = -val;
387  }
388  return val;
389 }
390 //============================================================================================================
391 }
392 
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:39
Header for a simple thermodynamics model of a surface phase derived from ThermoPhase, assuming an ideal solution model (see Thermodynamic Properties and class SurfPhase).
std::string kineticsSpeciesName(size_t k) const
Return the name of the kth species in the kinetics manager.
Definition: Kinetics.cpp:359
int m_nReactants
Number of reactants in the global reaction.
std::vector< doublereal > m_ReactantStoich
Vector of reactant stoichiometries that make up the global reaction.
size_t nTotalSpecies() const
The total number of species in all phases participating in the kinetics mechanism.
Definition: Kinetics.h:297
void warn_deprecated(const std::string &method, const std::string &extra)
Print a warning indicating that method is deprecated.
Definition: global.cpp:78
static void addV(int kkinspec, double ps, std::vector< int > &m_Products, std::vector< doublereal > &m_ProductStoich)
add the species into the list of products or reactants
A kinetics manager for heterogeneous reaction mechanisms.
std::string fp2str(const double x, const std::string &fmt)
Convert a double into a c++ string.
Definition: stringUtils.cpp:28
Public interface for kinetics managers.
Definition: Kinetics.h:128
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:99
Headers for the DenseMatrix object, which deals with dense rectangular matrices and description of th...
Cantera::Kinetics * m_kinetics
This kinetics operator is associated with just one homogeneous phase, associated with tpList[0] phase...
size_t nReactions() const
Number of reactions in the reaction mechanism.
Definition: Kinetics.h:193
virtual double reactantStoichCoeff(size_t k, size_t i) const
Stoichiometric coefficient of species k as a reactant in reaction i.
Definition: Kinetics.cpp:428
std::vector< int > m_Reactants
Vector of reactants that make up the global reaction.
virtual double productStoichCoeff(size_t k, size_t i) const
Stoichiometric coefficient of species k as a product in reaction i.
Definition: Kinetics.cpp:433
Cantera::InterfaceKinetics * m_InterfaceKinetics
This kinetics operator is associated with multiple homogeneous and surface phases.