Cantera  2.4.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 
24 {
25  doublereal T = thermo().temperature();
26  doublereal P = thermo().pressure();
27  m_logStandConc = log(thermo().standardConcentration());
28  doublereal logT = log(T);
29 
30  if (T != m_temp) {
31  if (!m_rfn.empty()) {
32  m_rates.update(T, logT, m_rfn.data());
33  }
34 
35  if (!m_rfn_low.empty()) {
36  m_falloff_low_rates.update(T, logT, m_rfn_low.data());
37  m_falloff_high_rates.update(T, logT, m_rfn_high.data());
38  }
39  if (!falloff_work.empty()) {
40  m_falloffn.updateTemp(T, falloff_work.data());
41  }
42  updateKc();
43  m_ROP_ok = false;
44  }
45 
46  if (T != m_temp || P != m_pres) {
47  if (m_plog_rates.nReactions()) {
48  m_plog_rates.update(T, logT, m_rfn.data());
49  m_ROP_ok = false;
50  }
51 
52  if (m_cheb_rates.nReactions()) {
53  m_cheb_rates.update(T, logT, m_rfn.data());
54  m_ROP_ok = false;
55  }
56  }
57  m_pres = P;
58  m_temp = T;
59 }
60 
62 {
63  thermo().getActivityConcentrations(m_conc.data());
64  doublereal ctot = thermo().molarDensity();
65 
66  // 3-body reactions
67  if (!concm_3b_values.empty()) {
68  m_3b_concm.update(m_conc, ctot, concm_3b_values.data());
69  }
70 
71  // Falloff reactions
72  if (!concm_falloff_values.empty()) {
73  m_falloff_concm.update(m_conc, ctot, concm_falloff_values.data());
74  }
75 
76  // P-log reactions
77  if (m_plog_rates.nReactions()) {
78  double logP = log(thermo().pressure());
79  m_plog_rates.update_C(&logP);
80  }
81 
82  // Chebyshev reactions
83  if (m_cheb_rates.nReactions()) {
84  double log10P = log10(thermo().pressure());
85  m_cheb_rates.update_C(&log10P);
86  }
87 
88  m_ROP_ok = false;
89 }
90 
92 {
93  thermo().getStandardChemPotentials(m_grt.data());
94  fill(m_rkcn.begin(), m_rkcn.end(), 0.0);
95 
96  // compute Delta G^0 for all reversible reactions
97  getRevReactionDelta(m_grt.data(), m_rkcn.data());
98 
99  doublereal rrt = 1.0 / thermo().RT();
100  for (size_t i = 0; i < m_revindex.size(); i++) {
101  size_t irxn = m_revindex[i];
102  m_rkcn[irxn] = std::min(exp(m_rkcn[irxn]*rrt - m_dn[irxn]*m_logStandConc),
103  BigNumber);
104  }
105 
106  for (size_t i = 0; i != m_irrev.size(); ++i) {
107  m_rkcn[ m_irrev[i] ] = 0.0;
108  }
109 }
110 
112 {
113  update_rates_T();
114  thermo().getStandardChemPotentials(m_grt.data());
115  fill(m_rkcn.begin(), m_rkcn.end(), 0.0);
116 
117  // compute Delta G^0 for all reactions
118  getReactionDelta(m_grt.data(), m_rkcn.data());
119 
120  doublereal rrt = 1.0 / thermo().RT();
121  for (size_t i = 0; i < nReactions(); i++) {
122  kc[i] = exp(-m_rkcn[i]*rrt + m_dn[i]*m_logStandConc);
123  }
124 
125  // force an update of T-dependent properties, so that m_rkcn will
126  // be updated before it is used next.
127  m_temp = 0.0;
128 }
129 
130 void GasKinetics::processFalloffReactions()
131 {
132  // use m_ropr for temporary storage of reduced pressure
133  vector_fp& pr = m_ropr;
134 
135  for (size_t i = 0; i < m_falloff_low_rates.nReactions(); i++) {
136  pr[i] = concm_falloff_values[i] * m_rfn_low[i] / (m_rfn_high[i] + SmallNumber);
137  AssertFinite(pr[i], "GasKinetics::processFalloffReactions",
138  "pr[{}] is not finite.", i);
139  }
140 
141  m_falloffn.pr_to_falloff(pr.data(), falloff_work.data());
142 
143  for (size_t i = 0; i < m_falloff_low_rates.nReactions(); i++) {
144  if (reactionType(m_fallindx[i]) == FALLOFF_RXN) {
145  pr[i] *= m_rfn_high[i];
146  } else { // CHEMACT_RXN
147  pr[i] *= m_rfn_low[i];
148  }
149  }
150 
151  scatter_copy(pr.begin(), pr.begin() + m_falloff_low_rates.nReactions(),
152  m_ropf.begin(), m_fallindx.begin());
153 }
154 
155 void GasKinetics::updateROP()
156 {
157  update_rates_C();
158  update_rates_T();
159  if (m_ROP_ok) {
160  return;
161  }
162 
163  // copy rate coefficients into ropf
164  m_ropf = m_rfn;
165 
166  // multiply ropf by enhanced 3b conc for all 3b rxns
167  if (!concm_3b_values.empty()) {
168  m_3b_concm.multiply(m_ropf.data(), concm_3b_values.data());
169  }
170 
171  if (m_falloff_high_rates.nReactions()) {
172  processFalloffReactions();
173  }
174 
175  // multiply by perturbation factor
176  multiply_each(m_ropf.begin(), m_ropf.end(), m_perturb.begin());
177 
178  // copy the forward rates to the reverse rates
179  m_ropr = m_ropf;
180 
181  // for reverse rates computed from thermochemistry, multiply the forward
182  // rates copied into m_ropr by the reciprocals of the equilibrium constants
183  multiply_each(m_ropr.begin(), m_ropr.end(), m_rkcn.begin());
184 
185  // multiply ropf by concentration products
186  m_reactantStoich.multiply(m_conc.data(), m_ropf.data());
187 
188  // for reversible reactions, multiply ropr by concentration products
189  m_revProductStoich.multiply(m_conc.data(), m_ropr.data());
190 
191  for (size_t j = 0; j != nReactions(); ++j) {
192  m_ropnet[j] = m_ropf[j] - m_ropr[j];
193  }
194 
195  for (size_t i = 0; i < m_rfn.size(); i++) {
196  AssertFinite(m_rfn[i], "GasKinetics::updateROP",
197  "m_rfn[{}] is not finite.", i);
198  AssertFinite(m_ropf[i], "GasKinetics::updateROP",
199  "m_ropf[{}] is not finite.", i);
200  AssertFinite(m_ropr[i], "GasKinetics::updateROP",
201  "m_ropr[{}] is not finite.", i);
202  }
203  m_ROP_ok = true;
204 }
205 
206 void GasKinetics::getFwdRateConstants(doublereal* kfwd)
207 {
208  update_rates_C();
209  update_rates_T();
210 
211  // copy rate coefficients into ropf
212  m_ropf = m_rfn;
213 
214  // multiply ropf by enhanced 3b conc for all 3b rxns
215  if (!concm_3b_values.empty()) {
216  m_3b_concm.multiply(m_ropf.data(), concm_3b_values.data());
217  }
218 
219  if (m_falloff_high_rates.nReactions()) {
220  processFalloffReactions();
221  }
222 
223  // multiply by perturbation factor
224  multiply_each(m_ropf.begin(), m_ropf.end(), m_perturb.begin());
225 
226  for (size_t i = 0; i < nReactions(); i++) {
227  kfwd[i] = m_ropf[i];
228  }
229 }
230 
231 bool GasKinetics::addReaction(shared_ptr<Reaction> r)
232 {
233  // operations common to all reaction types
234  bool added = BulkKinetics::addReaction(r);
235  if (!added) {
236  return false;
237  }
238 
239  switch (r->reaction_type) {
240  case ELEMENTARY_RXN:
241  addElementaryReaction(dynamic_cast<ElementaryReaction&>(*r));
242  break;
243  case THREE_BODY_RXN:
244  addThreeBodyReaction(dynamic_cast<ThreeBodyReaction&>(*r));
245  break;
246  case FALLOFF_RXN:
247  case CHEMACT_RXN:
248  addFalloffReaction(dynamic_cast<FalloffReaction&>(*r));
249  break;
250  case PLOG_RXN:
251  addPlogReaction(dynamic_cast<PlogReaction&>(*r));
252  break;
253  case CHEBYSHEV_RXN:
254  addChebyshevReaction(dynamic_cast<ChebyshevReaction&>(*r));
255  break;
256  default:
257  throw CanteraError("GasKinetics::addReaction",
258  "Unknown reaction type specified: {}", r->reaction_type);
259  }
260  return true;
261 }
262 
263 void GasKinetics::addFalloffReaction(FalloffReaction& r)
264 {
265  // install high and low rate coeff calculators and extend the high and low
266  // rate coeff value vectors
267  size_t nfall = m_falloff_high_rates.nReactions();
268  m_falloff_high_rates.install(nfall, r.high_rate);
269  m_rfn_high.push_back(0.0);
270  m_falloff_low_rates.install(nfall, r.low_rate);
271  m_rfn_low.push_back(0.0);
272 
273  // add this reaction number to the list of falloff reactions
274  m_fallindx.push_back(nReactions()-1);
275  m_rfallindx[nReactions()-1] = nfall;
276 
277  // install the enhanced third-body concentration calculator
278  map<size_t, double> efficiencies;
279  for (const auto& eff : r.third_body.efficiencies) {
280  size_t k = kineticsSpeciesIndex(eff.first);
281  if (k != npos) {
282  efficiencies[k] = eff.second;
283  } else if (!m_skipUndeclaredThirdBodies) {
284  throw CanteraError("GasKinetics::addFalloffReaction", "Found "
285  "third-body efficiency for undefined species '" + eff.first +
286  "' while adding reaction '" + r.equation() + "'");
287  }
288  }
289  m_falloff_concm.install(nfall, efficiencies,
291  concm_falloff_values.resize(m_falloff_concm.workSize());
292 
293  // install the falloff function calculator for this reaction
294  m_falloffn.install(nfall, r.reaction_type, r.falloff);
295  falloff_work.resize(m_falloffn.workSize());
296 }
297 
298 void GasKinetics::addThreeBodyReaction(ThreeBodyReaction& r)
299 {
300  m_rates.install(nReactions()-1, r.rate);
301  map<size_t, double> efficiencies;
302  for (const auto& eff : r.third_body.efficiencies) {
303  size_t k = kineticsSpeciesIndex(eff.first);
304  if (k != npos) {
305  efficiencies[k] = eff.second;
306  } else if (!m_skipUndeclaredThirdBodies) {
307  throw CanteraError("GasKinetics::addThreeBodyReaction", "Found "
308  "third-body efficiency for undefined species '" + eff.first +
309  "' while adding reaction '" + r.equation() + "'");
310  }
311  }
312  m_3b_concm.install(nReactions()-1, efficiencies,
313  r.third_body.default_efficiency);
314  concm_3b_values.resize(m_3b_concm.workSize());
315 }
316 
317 void GasKinetics::addPlogReaction(PlogReaction& r)
318 {
319  m_plog_rates.install(nReactions()-1, r.rate);
320 }
321 
322 void GasKinetics::addChebyshevReaction(ChebyshevReaction& r)
323 {
324  m_cheb_rates.install(nReactions()-1, r.rate);
325 }
326 
327 void GasKinetics::modifyReaction(size_t i, shared_ptr<Reaction> rNew)
328 {
329  // operations common to all reaction types
331 
332  switch (rNew->reaction_type) {
333  case ELEMENTARY_RXN:
334  modifyElementaryReaction(i, dynamic_cast<ElementaryReaction&>(*rNew));
335  break;
336  case THREE_BODY_RXN:
337  modifyThreeBodyReaction(i, dynamic_cast<ThreeBodyReaction&>(*rNew));
338  break;
339  case FALLOFF_RXN:
340  case CHEMACT_RXN:
341  modifyFalloffReaction(i, dynamic_cast<FalloffReaction&>(*rNew));
342  break;
343  case PLOG_RXN:
344  modifyPlogReaction(i, dynamic_cast<PlogReaction&>(*rNew));
345  break;
346  case CHEBYSHEV_RXN:
347  modifyChebyshevReaction(i, dynamic_cast<ChebyshevReaction&>(*rNew));
348  break;
349  default:
350  throw CanteraError("GasKinetics::modifyReaction",
351  "Unknown reaction type specified: {}", rNew->reaction_type);
352  }
353 
354  // invalidate all cached data
355  m_ROP_ok = false;
356  m_temp += 0.1234;
357  m_pres += 0.1234;
358 }
359 
360 void GasKinetics::modifyThreeBodyReaction(size_t i, ThreeBodyReaction& r)
361 {
362  m_rates.replace(i, r.rate);
363 }
364 
365 void GasKinetics::modifyFalloffReaction(size_t i, FalloffReaction& r)
366 {
367  size_t iFall = m_rfallindx[i];
368  m_falloff_high_rates.replace(iFall, r.high_rate);
369  m_falloff_low_rates.replace(iFall, r.low_rate);
370  m_falloffn.replace(iFall, r.falloff);
371 }
372 
373 void GasKinetics::modifyPlogReaction(size_t i, PlogReaction& r)
374 {
375  m_plog_rates.replace(i, r.rate);
376 }
377 
378 void GasKinetics::modifyChebyshevReaction(size_t i, ChebyshevReaction& r)
379 {
380  m_cheb_rates.replace(i, r.rate);
381 }
382 
384 {
386  m_logp_ref = log(thermo().refPressure()) - log(GasConstant);
387 }
388 
389 void GasKinetics::invalidateCache()
390 {
391  BulkKinetics::invalidateCache();
392  m_pres += 0.13579;
393 }
394 
395 }
virtual void modifyReaction(size_t i, shared_ptr< Reaction > rNew)
Modify the rate expression associated with a reaction.
Definition: Kinetics.cpp:572
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:852
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:921
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:227
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:370
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:383
vector_fp m_ropnet
Net rate-of-progress for each reaction.
Definition: Kinetics.h:924
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:602
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
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:748
virtual void update_rates_T()
Update temperature-dependent portions of reaction rates and falloff functions.
Definition: GasKinetics.cpp:23
vector_fp m_ropf
Forward rate-of-progress for each reaction.
Definition: Kinetics.h:918
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:699
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:82
std::vector< size_t > m_irrev
Indices of irreversible reactions.
Definition: BulkKinetics.h:51
doublereal molarDensity() const
Molar density (kmol/m^3).
Definition: Phase.cpp:590
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:912
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:849
bool m_skipUndeclaredThirdBodies
Definition: Kinetics.h:930
void updateKc()
Update the equilibrium constants in molar units.
Definition: GasKinetics.cpp:91
vector_fp m_perturb
Vector of perturbation factors for each reaction&#39;s rate of progress vector.
Definition: Kinetics.h:864
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
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition: Kinetics.h:261
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:50
size_t nReactions() const
Number of reactions in the reaction mechanism.
Definition: Kinetics.h:135
#define AssertFinite(expr, procedure,...)
Throw an exception if the specified exception is not a finite number.
Definition: ctexceptions.h:271
std::vector< size_t > m_fallindx
Reaction index of each falloff reaction.
Definition: GasKinetics.h:72
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:245
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:915
virtual void update_rates_C()
Update properties that depend on concentrations.
Definition: GasKinetics.cpp:61
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:100
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:76
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:56
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:546
virtual void modifyReaction(size_t i, shared_ptr< Reaction > rNew)
Modify the rate expression associated with a reaction.
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
Rate1< Arrhenius > m_falloff_low_rates
Rate expressions for falloff reactions at the low-pressure limit.
Definition: GasKinetics.h:79
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:373
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