Cantera  2.3.0
RxnRates.cpp
Go to the documentation of this file.
1 //! @file RxnRates.cpp
2 
3 // This file is part of Cantera. See License.txt in the top-level directory or
4 // at http://www.cantera.org/license.txt for license and copyright information.
5 
7 #include "cantera/base/Array.h"
8 
9 namespace Cantera
10 {
12  : m_logA(-1.0E300)
13  , m_b(0.0)
14  , m_E(0.0)
15  , m_A(0.0)
16 {
17 }
18 
19 Arrhenius::Arrhenius(doublereal A, doublereal b, doublereal E)
20  : m_b(b)
21  , m_E(E)
22  , m_A(A)
23 {
24  if (m_A <= 0.0) {
25  m_logA = -1.0E300;
26  } else {
27  m_logA = std::log(m_A);
28  }
29 }
30 
31 SurfaceArrhenius::SurfaceArrhenius()
32  : m_b(0.0)
33  , m_E(0.0)
34  , m_A(0.0)
35  , m_acov(0.0)
36  , m_ecov(0.0)
37  , m_mcov(0.0)
38 {
39 }
40 
41 SurfaceArrhenius::SurfaceArrhenius(double A, double b, double Ta)
42  : m_b(b)
43  , m_E(Ta)
44  , m_A(A)
45  , m_acov(0.0)
46  , m_ecov(0.0)
47  , m_mcov(0.0)
48 {
49 }
50 
51 void SurfaceArrhenius::addCoverageDependence(size_t k, doublereal a,
52  doublereal m, doublereal e)
53 {
54  m_sp.push_back(k);
55  m_ac.push_back(a);
56  m_ec.push_back(e);
57  if (m != 0.0) {
58  m_msp.push_back(k);
59  m_mc.push_back(m);
60  }
61 }
62 
63 Plog::Plog(const std::multimap<double, Arrhenius>& rates)
64  : logP_(-1000)
65  , logP1_(1000)
66  , logP2_(-1000)
67  , rDeltaP_(-1.0)
68 {
69  size_t j = 0;
70  rates_.reserve(rates.size());
71  // Insert intermediate pressures
72  for (const auto& rate : rates) {
73  double logp = std::log(rate.first);
74  if (pressures_.empty() || pressures_.rbegin()->first != logp) {
75  // starting a new group
76  pressures_[logp] = {j, j+1};
77  } else {
78  // another rate expression at the same pressure
79  pressures_[logp].second = j+1;
80  }
81 
82  j++;
83  rates_.push_back(rate.second);
84  }
85 
86  // Duplicate the first and last groups to handle P < P_0 and P > P_N
87  pressures_.insert({-1000.0, pressures_.begin()->second});
88  pressures_.insert({1000.0, pressures_.rbegin()->second});
89 }
90 
91 void Plog::validate(const std::string& equation)
92 {
93  double T[] = {200.0, 500.0, 1000.0, 2000.0, 5000.0, 10000.0};
94  for (auto iter = pressures_.begin(); iter->first < 1000; iter++) {
95  update_C(&iter->first);
96  for (size_t i=0; i < 6; i++) {
97  double k = updateRC(log(T[i]), 1.0/T[i]);
98  if (!(k >= 0)) {
99  // k is NaN. Increment the iterator so that the error
100  // message will correctly indicate that the problematic rate
101  // expression is at the higher of the adjacent pressures.
102  throw CanteraError("Plog::validate",
103  "Invalid rate coefficient for reaction '{}'\nat P = {}, T = {}",
104  equation, std::exp((++iter)->first), T[i]);
105  }
106  }
107  }
108 }
109 
110 std::vector<std::pair<double, Arrhenius> > Plog::rates() const
111 {
112  std::vector<std::pair<double, Arrhenius> > R;
113  // initial preincrement to skip rate for P --> 0
114  for (auto iter = ++pressures_.begin();
115  iter->first < 1000; // skip rates for (P --> infinity)
116  ++iter) {
117  for (size_t i = iter->second.first;
118  i < iter->second.second;
119  i++) {
120  R.emplace_back(std::exp(iter->first), rates_[i]);
121  }
122  }
123  return R;
124 }
125 
126 
127 ChebyshevRate::ChebyshevRate(double Tmin, double Tmax, double Pmin, double Pmax,
128  const Array2D& coeffs)
129  : Tmin_(Tmin)
130  , Tmax_(Tmax)
131  , Pmin_(Pmin)
132  , Pmax_(Pmax)
133  , nP_(coeffs.nColumns())
134  , nT_(coeffs.nRows())
135  , chebCoeffs_(coeffs.nColumns() * coeffs.nRows(), 0.0)
136  , dotProd_(coeffs.nRows())
137 {
138  double logPmin = std::log10(Pmin);
139  double logPmax = std::log10(Pmax);
140  double TminInv = 1.0 / Tmin;
141  double TmaxInv = 1.0 / Tmax;
142 
143  TrNum_ = - TminInv - TmaxInv;
144  TrDen_ = 1.0 / (TmaxInv - TminInv);
145  PrNum_ = - logPmin - logPmax;
146  PrDen_ = 1.0 / (logPmax - logPmin);
147 
148  for (size_t t = 0; t < nT_; t++) {
149  for (size_t p = 0; p < nP_; p++) {
150  chebCoeffs_[nP_*t + p] = coeffs(t,p);
151  }
152  }
153 }
154 
155 }
doublereal updateRC(doublereal logT, doublereal recipT) const
Update the value the rate constant.
Definition: RxnRates.h:257
Plog()
Default constructor.
Definition: RxnRates.h:220
void addCoverageDependence(size_t k, doublereal a, doublereal m, doublereal e)
Add a coverage dependency for species k, with pre-exponential dependence a, rate constant exponential...
Definition: RxnRates.cpp:51
size_t nP_
number of points in the pressure direction
Definition: RxnRates.h:449
double Tmin() const
Minimum valid temperature [K].
Definition: RxnRates.h:405
double Pmin() const
Minimum valid pressure [Pa].
Definition: RxnRates.h:415
std::vector< std::pair< double, Arrhenius > > rates() const
Return the pressures and Arrhenius expressions which comprise this reaction.
Definition: RxnRates.cpp:110
Arrhenius()
Default constructor.
Definition: RxnRates.cpp:11
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition: Array.h:31
Header file for class Cantera::Array2D.
double TrDen_
terms appearing in the reduced temperature
Definition: RxnRates.h:446
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
void validate(const std::string &equation)
Check to make sure that the rate expression is finite over a range of temperatures at each interpolat...
Definition: RxnRates.cpp:91
const vector_fp & coeffs() const
Access the Chebyshev coefficients.
Definition: RxnRates.h:439
void update_C(const doublereal *c)
Update concentration-dependent parts of the rate coefficient.
Definition: RxnRates.h:227
size_t nT_
number of points in the temperature direction
Definition: RxnRates.h:450
std::map< double, std::pair< size_t, size_t > > pressures_
log(p) to (index range) in the rates_ vector
Definition: RxnRates.h:294
double Tmax() const
Maximum valid temperature [K].
Definition: RxnRates.h:410
ChebyshevRate()
Default constructor.
Definition: RxnRates.h:349
Namespace for the Cantera kernel.
Definition: application.cpp:29
double PrDen_
terms appearing in the reduced pressure
Definition: RxnRates.h:447
vector_fp chebCoeffs_
Chebyshev coefficients, length nP * nT.
Definition: RxnRates.h:451
double Pmax() const
Maximum valid pressure [Pa].
Definition: RxnRates.h:420