Cantera  2.3.0
GasKinetics.cpp
Go to the documentation of this file.
1 /**
2  * @file GasKinetics.cpp Homogeneous kinetics in ideal gases
3  */
4 
5 // This file is part of Cantera. See License.txt in the top-level directory or
6 // at http://www.cantera.org/license.txt for license and copyright information.
7 
9 
10 using namespace std;
11 
12 namespace Cantera
13 {
14 GasKinetics::GasKinetics(thermo_t* thermo) :
15  BulkKinetics(thermo),
16  m_logp_ref(0.0),
17  m_logc_ref(0.0),
18  m_logStandConc(0.0),
19  m_pres(0.0)
20 {
21 }
22 
23 Kinetics* GasKinetics::duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const
24 {
25  GasKinetics* gK = new GasKinetics(*this);
26  gK->assignShallowPointers(tpVector);
27  return gK;
28 }
29 
31 {
32  doublereal T = thermo().temperature();
33  doublereal P = thermo().pressure();
34  m_logStandConc = log(thermo().standardConcentration());
35  doublereal logT = log(T);
36 
37  if (T != m_temp) {
38  if (!m_rfn.empty()) {
39  m_rates.update(T, logT, m_rfn.data());
40  }
41 
42  if (!m_rfn_low.empty()) {
43  m_falloff_low_rates.update(T, logT, m_rfn_low.data());
44  m_falloff_high_rates.update(T, logT, m_rfn_high.data());
45  }
46  if (!falloff_work.empty()) {
47  m_falloffn.updateTemp(T, falloff_work.data());
48  }
49  updateKc();
50  m_ROP_ok = false;
51  }
52 
53  if (T != m_temp || P != m_pres) {
54  if (m_plog_rates.nReactions()) {
55  m_plog_rates.update(T, logT, m_rfn.data());
56  m_ROP_ok = false;
57  }
58 
59  if (m_cheb_rates.nReactions()) {
60  m_cheb_rates.update(T, logT, m_rfn.data());
61  m_ROP_ok = false;
62  }
63  }
64  m_pres = P;
65  m_temp = T;
66 }
67 
69 {
70  thermo().getActivityConcentrations(m_conc.data());
71  doublereal ctot = thermo().molarDensity();
72 
73  // 3-body reactions
74  if (!concm_3b_values.empty()) {
75  m_3b_concm.update(m_conc, ctot, concm_3b_values.data());
76  }
77 
78  // Falloff reactions
79  if (!concm_falloff_values.empty()) {
80  m_falloff_concm.update(m_conc, ctot, concm_falloff_values.data());
81  }
82 
83  // P-log reactions
84  if (m_plog_rates.nReactions()) {
85  double logP = log(thermo().pressure());
86  m_plog_rates.update_C(&logP);
87  }
88 
89  // Chebyshev reactions
90  if (m_cheb_rates.nReactions()) {
91  double log10P = log10(thermo().pressure());
92  m_cheb_rates.update_C(&log10P);
93  }
94 
95  m_ROP_ok = false;
96 }
97 
99 {
100  thermo().getStandardChemPotentials(m_grt.data());
101  fill(m_rkcn.begin(), m_rkcn.end(), 0.0);
102 
103  // compute Delta G^0 for all reversible reactions
104  getRevReactionDelta(m_grt.data(), m_rkcn.data());
105 
106  doublereal rrt = 1.0 / thermo().RT();
107  for (size_t i = 0; i < m_revindex.size(); i++) {
108  size_t irxn = m_revindex[i];
109  m_rkcn[irxn] = std::min(exp(m_rkcn[irxn]*rrt - m_dn[irxn]*m_logStandConc),
110  BigNumber);
111  }
112 
113  for (size_t i = 0; i != m_irrev.size(); ++i) {
114  m_rkcn[ m_irrev[i] ] = 0.0;
115  }
116 }
117 
119 {
120  update_rates_T();
121  thermo().getStandardChemPotentials(m_grt.data());
122  fill(m_rkcn.begin(), m_rkcn.end(), 0.0);
123 
124  // compute Delta G^0 for all reactions
125  getReactionDelta(m_grt.data(), m_rkcn.data());
126 
127  doublereal rrt = 1.0 / thermo().RT();
128  for (size_t i = 0; i < nReactions(); i++) {
129  kc[i] = exp(-m_rkcn[i]*rrt + m_dn[i]*m_logStandConc);
130  }
131 
132  // force an update of T-dependent properties, so that m_rkcn will
133  // be updated before it is used next.
134  m_temp = 0.0;
135 }
136 
137 void GasKinetics::processFalloffReactions()
138 {
139  // use m_ropr for temporary storage of reduced pressure
140  vector_fp& pr = m_ropr;
141 
142  for (size_t i = 0; i < m_falloff_low_rates.nReactions(); i++) {
143  pr[i] = concm_falloff_values[i] * m_rfn_low[i] / (m_rfn_high[i] + SmallNumber);
144  AssertFinite(pr[i], "GasKinetics::processFalloffReactions",
145  "pr[{}] is not finite.", i);
146  }
147 
148  m_falloffn.pr_to_falloff(pr.data(), falloff_work.data());
149 
150  for (size_t i = 0; i < m_falloff_low_rates.nReactions(); i++) {
151  if (reactionType(m_fallindx[i]) == FALLOFF_RXN) {
152  pr[i] *= m_rfn_high[i];
153  } else { // CHEMACT_RXN
154  pr[i] *= m_rfn_low[i];
155  }
156  }
157 
158  scatter_copy(pr.begin(), pr.begin() + m_falloff_low_rates.nReactions(),
159  m_ropf.begin(), m_fallindx.begin());
160 }
161 
162 void GasKinetics::updateROP()
163 {
164  update_rates_C();
165  update_rates_T();
166  if (m_ROP_ok) {
167  return;
168  }
169 
170  // copy rate coefficients into ropf
171  m_ropf = m_rfn;
172 
173  // multiply ropf by enhanced 3b conc for all 3b rxns
174  if (!concm_3b_values.empty()) {
175  m_3b_concm.multiply(m_ropf.data(), concm_3b_values.data());
176  }
177 
178  if (m_falloff_high_rates.nReactions()) {
179  processFalloffReactions();
180  }
181 
182  // multiply by perturbation factor
183  multiply_each(m_ropf.begin(), m_ropf.end(), m_perturb.begin());
184 
185  // copy the forward rates to the reverse rates
186  m_ropr = m_ropf;
187 
188  // for reverse rates computed from thermochemistry, multiply the forward
189  // rates copied into m_ropr by the reciprocals of the equilibrium constants
190  multiply_each(m_ropr.begin(), m_ropr.end(), m_rkcn.begin());
191 
192  // multiply ropf by concentration products
193  m_reactantStoich.multiply(m_conc.data(), m_ropf.data());
194 
195  // for reversible reactions, multiply ropr by concentration products
196  m_revProductStoich.multiply(m_conc.data(), m_ropr.data());
197 
198  for (size_t j = 0; j != nReactions(); ++j) {
199  m_ropnet[j] = m_ropf[j] - m_ropr[j];
200  }
201 
202  for (size_t i = 0; i < m_rfn.size(); i++) {
203  AssertFinite(m_rfn[i], "GasKinetics::updateROP",
204  "m_rfn[{}] is not finite.", i);
205  AssertFinite(m_ropf[i], "GasKinetics::updateROP",
206  "m_ropf[{}] is not finite.", i);
207  AssertFinite(m_ropr[i], "GasKinetics::updateROP",
208  "m_ropr[{}] is not finite.", i);
209  }
210  m_ROP_ok = true;
211 }
212 
213 void GasKinetics::getFwdRateConstants(doublereal* kfwd)
214 {
215  update_rates_C();
216  update_rates_T();
217 
218  // copy rate coefficients into ropf
219  m_ropf = m_rfn;
220 
221  // multiply ropf by enhanced 3b conc for all 3b rxns
222  if (!concm_3b_values.empty()) {
223  m_3b_concm.multiply(m_ropf.data(), concm_3b_values.data());
224  }
225 
226  if (m_falloff_high_rates.nReactions()) {
227  processFalloffReactions();
228  }
229 
230  // multiply by perturbation factor
231  multiply_each(m_ropf.begin(), m_ropf.end(), m_perturb.begin());
232 
233  for (size_t i = 0; i < nReactions(); i++) {
234  kfwd[i] = m_ropf[i];
235  }
236 }
237 
238 bool GasKinetics::addReaction(shared_ptr<Reaction> r)
239 {
240  // operations common to all reaction types
241  bool added = BulkKinetics::addReaction(r);
242  if (!added) {
243  return false;
244  }
245 
246  switch (r->reaction_type) {
247  case ELEMENTARY_RXN:
248  addElementaryReaction(dynamic_cast<ElementaryReaction&>(*r));
249  break;
250  case THREE_BODY_RXN:
251  addThreeBodyReaction(dynamic_cast<ThreeBodyReaction&>(*r));
252  break;
253  case FALLOFF_RXN:
254  case CHEMACT_RXN:
255  addFalloffReaction(dynamic_cast<FalloffReaction&>(*r));
256  break;
257  case PLOG_RXN:
258  addPlogReaction(dynamic_cast<PlogReaction&>(*r));
259  break;
260  case CHEBYSHEV_RXN:
261  addChebyshevReaction(dynamic_cast<ChebyshevReaction&>(*r));
262  break;
263  default:
264  throw CanteraError("GasKinetics::addReaction",
265  "Unknown reaction type specified: {}", r->reaction_type);
266  }
267  return true;
268 }
269 
270 void GasKinetics::addFalloffReaction(FalloffReaction& r)
271 {
272  // install high and low rate coeff calculators and extend the high and low
273  // rate coeff value vectors
274  size_t nfall = m_falloff_high_rates.nReactions();
275  m_falloff_high_rates.install(nfall, r.high_rate);
276  m_rfn_high.push_back(0.0);
277  m_falloff_low_rates.install(nfall, r.low_rate);
278  m_rfn_low.push_back(0.0);
279 
280  // add this reaction number to the list of falloff reactions
281  m_fallindx.push_back(nReactions()-1);
282  m_rfallindx[nReactions()-1] = nfall;
283 
284  // install the enhanced third-body concentration calculator
285  map<size_t, double> efficiencies;
286  for (const auto& eff : r.third_body.efficiencies) {
287  size_t k = kineticsSpeciesIndex(eff.first);
288  if (k != npos) {
289  efficiencies[k] = eff.second;
290  } else if (!m_skipUndeclaredThirdBodies) {
291  throw CanteraError("GasKinetics::addFalloffReaction", "Found "
292  "third-body efficiency for undefined species '" + eff.first +
293  "' while adding reaction '" + r.equation() + "'");
294  }
295  }
296  m_falloff_concm.install(nfall, efficiencies,
298  concm_falloff_values.resize(m_falloff_concm.workSize());
299 
300  // install the falloff function calculator for this reaction
301  m_falloffn.install(nfall, r.reaction_type, r.falloff);
302  falloff_work.resize(m_falloffn.workSize());
303 }
304 
305 void GasKinetics::addThreeBodyReaction(ThreeBodyReaction& r)
306 {
307  m_rates.install(nReactions()-1, r.rate);
308  map<size_t, double> efficiencies;
309  for (const auto& eff : r.third_body.efficiencies) {
310  size_t k = kineticsSpeciesIndex(eff.first);
311  if (k != npos) {
312  efficiencies[k] = eff.second;
313  } else if (!m_skipUndeclaredThirdBodies) {
314  throw CanteraError("GasKinetics::addThreeBodyReaction", "Found "
315  "third-body efficiency for undefined species '" + eff.first +
316  "' while adding reaction '" + r.equation() + "'");
317  }
318  }
319  m_3b_concm.install(nReactions()-1, efficiencies,
320  r.third_body.default_efficiency);
321  concm_3b_values.resize(m_3b_concm.workSize());
322 }
323 
324 void GasKinetics::addPlogReaction(PlogReaction& r)
325 {
326  m_plog_rates.install(nReactions()-1, r.rate);
327 }
328 
329 void GasKinetics::addChebyshevReaction(ChebyshevReaction& r)
330 {
331  m_cheb_rates.install(nReactions()-1, r.rate);
332 }
333 
334 void GasKinetics::modifyReaction(size_t i, shared_ptr<Reaction> rNew)
335 {
336  // operations common to all reaction types
338 
339  switch (rNew->reaction_type) {
340  case ELEMENTARY_RXN:
341  modifyElementaryReaction(i, dynamic_cast<ElementaryReaction&>(*rNew));
342  break;
343  case THREE_BODY_RXN:
344  modifyThreeBodyReaction(i, dynamic_cast<ThreeBodyReaction&>(*rNew));
345  break;
346  case FALLOFF_RXN:
347  case CHEMACT_RXN:
348  modifyFalloffReaction(i, dynamic_cast<FalloffReaction&>(*rNew));
349  break;
350  case PLOG_RXN:
351  modifyPlogReaction(i, dynamic_cast<PlogReaction&>(*rNew));
352  break;
353  case CHEBYSHEV_RXN:
354  modifyChebyshevReaction(i, dynamic_cast<ChebyshevReaction&>(*rNew));
355  break;
356  default:
357  throw CanteraError("GasKinetics::modifyReaction",
358  "Unknown reaction type specified: {}", rNew->reaction_type);
359  }
360 
361  // invalidate all cached data
362  m_ROP_ok = false;
363  m_temp += 0.1234;
364  m_pres += 0.1234;
365 }
366 
367 void GasKinetics::modifyThreeBodyReaction(size_t i, ThreeBodyReaction& r)
368 {
369  m_rates.replace(i, r.rate);
370 }
371 
372 void GasKinetics::modifyFalloffReaction(size_t i, FalloffReaction& r)
373 {
374  size_t iFall = m_rfallindx[i];
375  m_falloff_high_rates.replace(iFall, r.high_rate);
376  m_falloff_low_rates.replace(iFall, r.low_rate);
377  m_falloffn.replace(iFall, r.falloff);
378 }
379 
380 void GasKinetics::modifyPlogReaction(size_t i, PlogReaction& r)
381 {
382  m_plog_rates.replace(i, r.rate);
383 }
384 
385 void GasKinetics::modifyChebyshevReaction(size_t i, ChebyshevReaction& r)
386 {
387  m_cheb_rates.replace(i, r.rate);
388 }
389 
391 {
393  m_logp_ref = log(thermo().refPressure()) - log(GasConstant);
394 }
395 
396 void GasKinetics::invalidateCache()
397 {
398  BulkKinetics::invalidateCache();
399  m_pres += 0.13579;
400 }
401 
402 }
virtual void modifyReaction(size_t i, shared_ptr< Reaction > rNew)
Modify the rate expression associated with a reaction.
Definition: Kinetics.cpp:622
const int PLOG_RXN
A pressure-dependent rate expression consisting of several Arrhenius rate expressions evaluated at di...
Definition: reaction_defs.h:51
StoichManagerN m_revProductStoich
Stoichiometry manager for the products of reversible reactions.
Definition: Kinetics.h:913
void pr_to_falloff(doublereal *values, const doublereal *work)
Given a vector of reduced pressures for each falloff reaction, replace each entry by the value of the...
Definition: FalloffMgr.h:79
vector_fp m_ropr
Reverse rate-of-progress for each reaction.
Definition: Kinetics.h:982
int reaction_type
Type of the reaction.
Definition: Reaction.h:45
doublereal temperature() const
Temperature (K).
Definition: Phase.h:601
thermo_t & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism...
Definition: Kinetics.h:276
void replace(size_t rxn, shared_ptr< Falloff > f)
Definition: FalloffMgr.h:54
virtual void getFwdRateConstants(doublereal *kfwd)
Return the forward rate constants.
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
void install(size_t rxn, int reactionType, shared_ptr< Falloff > f)
Install a new falloff function calculator.
Definition: FalloffMgr.h:39
const int CHEBYSHEV_RXN
A general gas-phase pressure-dependent reaction where k(T,P) is defined in terms of a bivariate Cheby...
Definition: reaction_defs.h:57
virtual void init()
Prepare the class for the addition of reactions, after all phases have been added.
virtual void getActivityConcentrations(doublereal *c) const
This method returns an array of generalized concentrations.
Definition: ThermoPhase.h:403
STL namespace.
virtual void getRevReactionDelta(const doublereal *g, doublereal *dg)
Given an array of species properties &#39;g&#39;, return in array &#39;dg&#39; the change in this quantity in the rev...
Definition: Kinetics.cpp:438
vector_fp m_ropnet
Net rate-of-progress for each reaction.
Definition: Kinetics.h:985
const int CHEMACT_RXN
A chemical activation reaction.
Definition: reaction_defs.h:65
virtual int reactionType(size_t i) const
Flag specifying the type of reaction.
Definition: Kinetics.h:651
void updateTemp(doublereal t, doublereal *work)
Update the cached temperature-dependent intermediate results for all installed falloff functions...
Definition: FalloffMgr.h:69
size_t workSize()
Size of the work array required to store intermediate results.
Definition: FalloffMgr.h:59
Kinetics manager for elementary gas-phase chemistry.
Definition: GasKinetics.h:26
void multiply_each(OutputIter x_begin, OutputIter x_end, InputIter y_begin)
Multiply each entry in x by the corresponding entry in y.
Definition: utilities.h:162
doublereal RT() const
Return the Gas Constant multiplied by the current temperature.
Definition: ThermoPhase.h:809
virtual void update_rates_T()
Update temperature-dependent portions of reaction rates and falloff functions.
Definition: GasKinetics.cpp:30
vector_fp m_ropf
Forward rate-of-progress for each reaction.
Definition: Kinetics.h:979
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:93
Partial specialization of Kinetics for chemistry in a single bulk phase.
Definition: BulkKinetics.h:21
const int FALLOFF_RXN
The general form for a gas-phase association or dissociation reaction, with a pressure-dependent rate...
Definition: reaction_defs.h:43
virtual void init()
Prepare the class for the addition of reactions, after all phases have been added.
Definition: Kinetics.h:743
virtual void getEquilibriumConstants(doublereal *kc)
Return a vector of Equilibrium constants.
shared_ptr< Falloff > falloff
Falloff function which determines how low_rate and high_rate are combined to determine the rate const...
Definition: Reaction.h:150
Rate1< Arrhenius > m_falloff_high_rates
Rate expressions for falloff reactions at the high-pressure limit.
Definition: GasKinetics.h:90
std::vector< size_t > m_irrev
Indices of irreversible reactions.
Definition: BulkKinetics.h:52
doublereal molarDensity() const
Molar density (kmol/m^3).
Definition: Phase.cpp:666
const doublereal BigNumber
largest number to compare to inf.
Definition: ct_defs.h:128
vector_fp m_rfn
Forward rate constant for each reaction.
Definition: Kinetics.h:973
virtual bool addReaction(shared_ptr< Reaction > r)
Add a single reaction to the mechanism.
StoichManagerN m_reactantStoich
Stoichiometry manager for the reactants for each reaction.
Definition: Kinetics.h:910
bool m_skipUndeclaredThirdBodies
Definition: Kinetics.h:991
virtual void assignShallowPointers(const std::vector< thermo_t *> &tpVector)
Reassign the pointers within the Kinetics object.
Definition: Kinetics.cpp:129
void updateKc()
Update the equilibrium constants in molar units.
Definition: GasKinetics.cpp:98
vector_fp m_perturb
Vector of perturbation factors for each reaction&#39;s rate of progress vector.
Definition: Kinetics.h:925
A reaction that is first-order in [M] at low pressure, like a third-body reaction, but zeroth-order in [M] as pressure increases.
Definition: Reaction.h:128
Public interface for kinetics managers.
Definition: Kinetics.h:111
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition: Kinetics.h:310
GasKinetics(thermo_t *thermo=0)
Constructor.
Definition: GasKinetics.cpp:14
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
std::vector< size_t > m_revindex
Indices of reversible reactions.
Definition: BulkKinetics.h:51
size_t nReactions() const
Number of reactions in the reaction mechanism.
Definition: Kinetics.h:184
#define AssertFinite(expr, procedure,...)
Throw an exception if the specified exception is not a finite number.
Definition: ctexceptions.h:277
std::vector< size_t > m_fallindx
Reaction index of each falloff reaction.
Definition: GasKinetics.h:80
virtual bool addReaction(shared_ptr< Reaction > r)
Add a single reaction to the mechanism.
virtual doublereal pressure() const
Return the thermodynamic pressure (Pa).
Definition: ThermoPhase.h:278
const int THREE_BODY_RXN
A gas-phase reaction that requires a third-body collision partner.
Definition: reaction_defs.h:37
vector_fp m_rkcn
Reciprocal of the equilibrium constant in concentration units.
Definition: Kinetics.h:976
virtual void update_rates_C()
Update properties that depend on concentrations.
Definition: GasKinetics.cpp:68
virtual Kinetics * duplMyselfAsKinetics(const std::vector< thermo_t *> &tpVector) const
Duplication routine for objects which inherit from Kinetics.
Definition: GasKinetics.cpp:23
const doublereal SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:126
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:157
Composition efficiencies
Map of species to third body efficiency.
Definition: Reaction.h:103
doublereal m_pres
Last pressure at which rates were evaluated.
Definition: GasKinetics.h:108
std::map< size_t, size_t > m_rfallindx
Map of reaction index to falloff reaction index (i.e indices in m_falloff_low_rates and m_falloff_hig...
Definition: GasKinetics.h:84
A reaction with a non-reacting third body "M" that acts to add or remove energy from the reacting spe...
Definition: Reaction.h:112
void scatter_copy(InputIter begin, InputIter end, OutputIter result, IndexIter index)
Copies a contiguous range in a sequence to indexed positions in another sequence. ...
Definition: utilities.h:330
vector_fp m_dn
Difference between the global reactants order and the global products order.
Definition: BulkKinetics.h:57
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:64
const int ELEMENTARY_RXN
A reaction with a rate coefficient that depends only on temperature and voltage that also obeys mass-...
Definition: reaction_defs.h:31
Arrhenius high_rate
The rate constant in the high-pressure limit.
Definition: Reaction.h:143
Arrhenius low_rate
The rate constant in the low-pressure limit.
Definition: Reaction.h:140
virtual void getStandardChemPotentials(doublereal *mu) const
Get the array of chemical potentials at unit activity for the species at their standard states at the...
Definition: ThermoPhase.h:579
virtual void modifyReaction(size_t i, shared_ptr< Reaction > rNew)
Modify the rate expression associated with a reaction.
Namespace for the Cantera kernel.
Definition: application.cpp:29
Rate1< Arrhenius > m_falloff_low_rates
Rate expressions for falloff reactions at the low-pressure limit.
Definition: GasKinetics.h:87
double default_efficiency
The default third body efficiency for species not listed in efficiencies.
Definition: Reaction.h:107
virtual void getReactionDelta(const doublereal *property, doublereal *deltaProperty)
Change in species properties.
Definition: Kinetics.cpp:428
ThirdBody third_body
Relative efficiencies of third-body species in enhancing the reaction rate.
Definition: Reaction.h:146
std::string equation() const
The chemical equation for this reaction.
Definition: Reaction.cpp:89