Cantera  3.1.0a1
Loading...
Searching...
No Matches
Falloff.cpp
Go to the documentation of this file.
1/**
2 * @file Falloff.cpp Definitions for member functions of classes derived from
3 * Falloff
4 */
5
6// This file is part of Cantera. See License.txt in the top-level directory or
7// at https://cantera.org/license.txt for license and copyright information.
8
11#include "cantera/base/global.h"
12#include "cantera/base/AnyMap.h"
13#include "cantera/kinetics/Falloff.h"
16
17namespace Cantera
18{
19
20FalloffData::FalloffData()
21{
22 conc_3b.resize(1, NAN);
23 m_conc_3b_buf.resize(1, NAN);
24}
25
26void FalloffData::update(double T)
27{
28 throw CanteraError("FalloffData::update",
29 "Missing state information: 'FalloffData' requires third-body concentration.");
30}
31
32void FalloffData::update(double T, double M)
33{
35 conc_3b[0] = M;
36}
37
38bool FalloffData::update(const ThermoPhase& phase, const Kinetics& kin)
39{
40 double rho_m = phase.molarDensity();
41 int mf = phase.stateMFNumber();
42 double T = phase.temperature();
43 bool changed = false;
44 if (T != temperature) {
46 changed = true;
47 }
48 if (rho_m != molar_density || mf != m_state_mf_number) {
49 molar_density = rho_m;
52 changed = true;
53 }
54 return changed;
55}
56
58{
59 if (m_perturbed) {
60 throw CanteraError("FalloffData::perturbThirdBodies",
61 "Cannot apply another perturbation as state is already perturbed.");
62 }
64 for (auto& c3b : conc_3b) {
65 c3b *= 1. + deltaM;
66 }
67 m_perturbed = true;
68}
69
71{
73 // only restore if there is a valid buffered value
74 if (!m_perturbed) {
75 return;
76 }
78 m_perturbed = false;
79}
80
81FalloffRate::FalloffRate(const AnyMap& node, const UnitStack& rate_units)
82 : FalloffRate()
83{
84 setParameters(node, rate_units);
85}
86
88{
89 ArrheniusRate _low = low;
91 _low.check(m_input.getString("equation", ""));
93 throw CanteraError("FalloffRate::setLowRate",
94 "Detected inconsistent rate definitions;\nhigh and low "
95 "rate pre-exponential factors must have the same sign.");
96 }
97 m_lowRate = std::move(_low);
98}
99
101{
102 ArrheniusRate _high = high;
104 _high.check(m_input.getString("equation", ""));
106 throw CanteraError("FalloffRate::setHighRate",
107 "Detected inconsistent rate definitions;\nhigh and low "
108 "rate pre-exponential factors must have the same sign.");
109 }
110 m_highRate = std::move(_high);
111}
112
113void FalloffRate::setFalloffCoeffs(const vector<double>& c)
114{
115 if (c.size() != 0) {
116 throw InputFileError("FalloffRate::setFalloffCoeffs", m_input,
117 "Incorrect number of parameters. 0 required. Received {}.",
118 c.size());
119 }
120 m_valid = true;
121}
122
123void FalloffRate::getFalloffCoeffs(vector<double>& c) const
124{
125 c.clear();
126}
127
128void FalloffRate::setParameters(const AnyMap& node, const UnitStack& rate_units)
129{
130 ReactionRate::setParameters(node, rate_units);
131 if (node.empty()) {
132 return;
133 }
134
135 m_negativeA_ok = node.getBool("negative-A", false);
136 if (node["type"] == "chemically-activated") {
138 }
139
140 UnitStack low_rate_units = rate_units;
141 UnitStack high_rate_units = rate_units;
142 if (rate_units.size()) {
144 high_rate_units.join(1);
145 } else {
146 low_rate_units.join(-1);
147 }
148 }
149 if (node.hasKey("low-P-rate-constant")) {
151 node["low-P-rate-constant"], node.units(), low_rate_units);
153 }
154 if (node.hasKey("high-P-rate-constant")) {
156 node["high-P-rate-constant"], node.units(), high_rate_units);
158 }
159}
160
162{
163 if (m_negativeA_ok) {
164 node["negative-A"] = true;
165 }
166 AnyMap rNode;
168 if (!rNode.empty()) {
169 node["low-P-rate-constant"] = std::move(rNode);
170 }
171 rNode.clear();
173 if (!rNode.empty()) {
174 node["high-P-rate-constant"] = std::move(rNode);
175 }
176}
177
178void FalloffRate::check(const string& equation)
179{
180 m_lowRate.check(equation);
181 m_highRate.check(equation);
182 if (!m_lowRate.valid() || !m_highRate.valid()) {
183 // arrhenius rates are not initialized
184 return;
185 }
187 throw InputFileError("FalloffRate::check", m_input,
188 "Inconsistent rate definitions found in reaction '{}';\nhigh and low "
189 "rate pre-exponential factors must have the same sign.", equation);
190 }
191}
192
193void FalloffRate::validate(const string& equation, const Kinetics& kin)
194{
195 try {
196 m_lowRate.validate(equation, kin);
197 m_highRate.validate(equation, kin);
198 } catch (CanteraError& err) {
199 throw InputFileError("FalloffRate::validate", m_input, err.getMessage());
200 }
201}
202
203LindemannRate::LindemannRate(const AnyMap& node, const UnitStack& rate_units)
204 : LindemannRate()
205{
206 setParameters(node, rate_units);
207}
208
209LindemannRate::LindemannRate(const ArrheniusRate& low, const ArrheniusRate& high,
210 const vector<double>& c)
211 : LindemannRate()
212{
213 m_lowRate = low;
214 m_highRate = high;
216}
217
218TroeRate::TroeRate(const AnyMap& node, const UnitStack& rate_units)
219 : TroeRate()
220{
221 setParameters(node, rate_units);
222}
223
224TroeRate::TroeRate(const ArrheniusRate& low, const ArrheniusRate& high,
225 const vector<double>& c)
226 : TroeRate()
227{
228 m_lowRate = low;
229 m_highRate = high;
231}
232
233void TroeRate::setFalloffCoeffs(const vector<double>& c)
234{
235 if (c.size() != 3 && c.size() != 4) {
236 throw InputFileError("TroeRate::setFalloffCoeffs", m_input,
237 "Incorrect number of coefficients. 3 or 4 required. Received {}.",
238 c.size());
239 }
240 m_a = c[0];
241 if (std::abs(c[1]) < SmallNumber) {
242 m_rt3 = std::numeric_limits<double>::infinity();
243 } else {
244 m_rt3 = 1.0 / c[1];
245 }
246
247 if (std::abs(c[2]) < SmallNumber) {
248 m_rt1 = std::numeric_limits<double>::infinity();
249 } else {
250 m_rt1 = 1.0 / c[2];
251 }
252
253 if (c.size() == 4) {
254 if (std::abs(c[3]) < SmallNumber) {
255 warn_user("TroeRate::setFalloffCoeffs",
256 "Unexpected parameter value T2=0. Omitting exp(T2/T) term from "
257 "falloff expression. To suppress this warning, remove value "
258 "for T2 from the input file. In the unlikely case that the "
259 "exp(T2/T) term should be included with T2 effectively equal "
260 "to 0, set T2 to a sufficiently small value "
261 "(for example, T2 < 1e-16).");
262 }
263 m_t2 = c[3];
264 } else {
265 m_t2 = 0.;
266 }
267 m_valid = true;
268}
269
270void TroeRate::getFalloffCoeffs(vector<double>& c) const
271{
272 if (std::abs(m_t2) < SmallNumber) {
273 c.resize(3);
274 } else {
275 c.resize(4, 0.);
276 c[3] = m_t2;
277 }
278 c[0] = m_a;
279 c[1] = 1.0 / m_rt3;
280 c[2] = 1.0 / m_rt1;
281}
282
283void TroeRate::updateTemp(double T, double* work) const
284{
285 double Fcent = (1.0 - m_a) * exp(-T*m_rt3) + m_a * exp(-T*m_rt1);
286 if (m_t2) {
287 Fcent += exp(- m_t2 / T);
288 }
289 *work = log10(std::max(Fcent, SmallNumber));
290}
291
292double TroeRate::F(double pr, const double* work) const
293{
294 double lpr = log10(std::max(pr,SmallNumber));
295 double cc = -0.4 - 0.67 * (*work);
296 double nn = 0.75 - 1.27 * (*work);
297 double f1 = (lpr + cc)/ (nn - 0.14 * (lpr + cc));
298 double lgf = (*work) / (1.0 + f1 * f1);
299 return pow(10.0, lgf);
300}
301
302void TroeRate::setParameters(const AnyMap& node, const UnitStack& rate_units)
303{
304 if (node.empty()) {
305 return;
306 }
307
308 FalloffRate::setParameters(node, rate_units);
309 auto& f = node["Troe"].as<AnyMap>();
310 if (f.empty()) {
311 return;
312 }
313 vector<double> params{
314 f["A"].asDouble(),
315 f["T3"].asDouble(),
316 f["T1"].asDouble()
317 };
318 if (f.hasKey("T2")) {
319 params.push_back(f["T2"].asDouble());
320 }
321 setFalloffCoeffs(params);
322}
323
325{
327
328 AnyMap params;
329 if (valid()) {
330 params["A"] = m_a;
331 params["T3"].setQuantity(1.0 / m_rt3, "K");
332 params["T1"].setQuantity(1.0 / m_rt1, "K");
333 if (std::abs(m_t2) > SmallNumber) {
334 params["T2"].setQuantity(m_t2, "K");
335 }
336 }
337 params.setFlowStyle();
338 node["Troe"] = std::move(params);
339}
340
341SriRate::SriRate(const AnyMap& node, const UnitStack& rate_units)
342 : SriRate()
343{
344 setParameters(node, rate_units);
345}
346
347void SriRate::setFalloffCoeffs(const vector<double>& c)
348{
349 if (c.size() != 3 && c.size() != 5) {
350 throw InputFileError("SriRate::setFalloffCoeffs", m_input,
351 "Incorrect number of coefficients. 3 or 5 required. Received {}.",
352 c.size());
353 }
354
355 if (c[2] < 0.0) {
356 throw InputFileError("SriRate::setFalloffCoeffs()", m_input,
357 "m_c parameter is less than zero: {}", c[2]);
358 }
359 m_a = c[0];
360 m_b = c[1];
361 m_c = c[2];
362
363 if (c.size() == 5) {
364 if (c[3] < 0.0) {
365 throw InputFileError("SriRate::setFalloffCoeffs()", m_input,
366 "m_d parameter is less than zero: {}", c[3]);
367 }
368 m_d = c[3];
369 m_e = c[4];
370 } else {
371 m_d = 1.0;
372 m_e = 0.0;
373 }
374 m_valid = true;
375}
376
377void SriRate::getFalloffCoeffs(vector<double>& c) const
378{
379 if (m_e < SmallNumber && std::abs(m_e - 1.) < SmallNumber) {
380 c.resize(3);
381 } else {
382 c.resize(5, 0.);
383 c[3] = m_d;
384 c[4] = m_e;
385 }
386 c[0] = m_a;
387 c[1] = m_b;
388 c[2] = m_c;
389}
390
391void SriRate::updateTemp(double T, double* work) const
392{
393 *work = m_a * exp(- m_b / T);
394 if (m_c != 0.0) {
395 *work += exp(- T/m_c);
396 }
397 work[1] = m_d * pow(T,m_e);
398}
399
400double SriRate::F(double pr, const double* work) const
401{
402 double lpr = log10(std::max(pr,SmallNumber));
403 double xx = 1.0/(1.0 + lpr*lpr);
404 return pow(*work, xx) * work[1];
405}
406
407void SriRate::setParameters(const AnyMap& node, const UnitStack& rate_units)
408{
409 if (node.empty()) {
410 return;
411 }
412
413 FalloffRate::setParameters(node, rate_units);
414 auto& f = node["SRI"].as<AnyMap>();
415 if (f.empty()) {
416 return;
417 }
418 vector<double> params{
419 f["A"].asDouble(),
420 f["B"].asDouble(),
421 f["C"].asDouble()
422 };
423 if (f.hasKey("D")) {
424 params.push_back(f["D"].asDouble());
425 }
426 if (f.hasKey("E")) {
427 params.push_back(f["E"].asDouble());
428 }
429 setFalloffCoeffs(params);
430}
431
433{
435
436 AnyMap params;
437 if (valid()) {
438 params["A"] = m_a;
439 params["B"].setQuantity(m_b, "K");
440 params["C"].setQuantity(m_c, "K");
441 if (m_d != 1.0 || m_e != 0.0) {
442 params["D"] = m_d;
443 params["E"] = m_e;
444 }
445 }
446 params.setFlowStyle();
447 node["SRI"] = std::move(params);
448}
449
450TsangRate::TsangRate(const AnyMap& node, const UnitStack& rate_units)
451 : TsangRate()
452{
453 setParameters(node, rate_units);
454}
455
456void TsangRate::setFalloffCoeffs(const vector<double>& c)
457{
458 if (c.size() != 1 && c.size() != 2) {
459 throw InputFileError("TsangRate::init", m_input,
460 "Incorrect number of coefficients. 1 or 2 required. Received {}.",
461 c.size());
462 }
463 m_a = c[0];
464
465 if (c.size() == 2) {
466 m_b = c[1];
467 }
468 else {
469 m_b = 0.0;
470 }
471 m_valid = true;
472}
473
474void TsangRate::getFalloffCoeffs(vector<double>& c) const
475{
476 if (std::abs(m_b) < SmallNumber) {
477 c.resize(1);
478 } else {
479 c.resize(2, 0.);
480 c[1] = m_b;
481 }
482 c[0] = m_a;
483}
484
485void TsangRate::updateTemp(double T, double* work) const
486{
487 double Fcent = m_a + (m_b * T);
488 *work = log10(std::max(Fcent, SmallNumber));
489}
490
491double TsangRate::F(double pr, const double* work) const
492{ //identical to TroeRate::F
493 double lpr = log10(std::max(pr,SmallNumber));
494 double cc = -0.4 - 0.67 * (*work);
495 double nn = 0.75 - 1.27 * (*work);
496 double f1 = (lpr + cc)/ (nn - 0.14 * (lpr + cc));
497 double lgf = (*work) / (1.0 + f1 * f1);
498 return pow(10.0, lgf);
499}
500
501void TsangRate::setParameters(const AnyMap& node, const UnitStack& rate_units)
502{
503 if (node.empty()) {
504 return;
505 }
506
507 FalloffRate::setParameters(node, rate_units);
508 auto& f = node["Tsang"].as<AnyMap>();
509 if (f.empty()) {
510 return;
511 }
512 vector<double> params{
513 f["A"].asDouble(),
514 f["B"].asDouble()
515 };
516 setFalloffCoeffs(params);
517}
518
520{
522
523 AnyMap params;
524 if (!valid()) {
525 // pass
526 } else {
527 // Parameters do not have unit system (yet)
528 params["A"] = m_a;
529 params["B"] = m_b;
530 }
531 params.setFlowStyle();
532 node["Tsang"] = std::move(params);
533}
534
535}
Base class for kinetics managers and also contains the kineticsmgr module documentation (see Kinetics...
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:427
size_t size() const
Returns the number of elements in this map.
Definition AnyMap.h:622
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition AnyMap.cpp:1423
const UnitSystem & units() const
Return the default units that should be used to convert stored values.
Definition AnyMap.h:630
bool empty() const
Return boolean indicating whether AnyMap is empty.
Definition AnyMap.cpp:1418
void setFlowStyle(bool flow=true)
Use "flow" style when outputting this AnyMap to YAML.
Definition AnyMap.cpp:1726
bool getBool(const string &key, bool default_) const
If key exists, return it as a bool, otherwise return default_.
Definition AnyMap.cpp:1515
void clear()
Erase all items in the mapping.
Definition AnyMap.cpp:1433
const string & getString(const string &key, const string &default_) const
If key exists, return it as a string, otherwise return default_.
Definition AnyMap.cpp:1530
void setAllowNegativePreExponentialFactor(bool value)
Set flag indicating whether negative A values are permitted.
Definition Arrhenius.h:141
void getRateParameters(AnyMap &node) const
Get Arrhenius parameters used to populate the rate-coefficient or equivalent field.
Definition Arrhenius.cpp:76
void validate(const string &equation, const Kinetics &kin) override
Validate the reaction rate expression.
virtual double preExponentialFactor() const
Return the pre-exponential factor A (in m, kmol, s to powers depending on the reaction order)
Definition Arrhenius.h:97
void check(const string &equation) override
Check rate expression.
Arrhenius reaction rate type depends only on temperature.
Definition Arrhenius.h:170
Base class for exceptions thrown by Cantera classes.
virtual string getMessage() const
Method overridden by derived classes to format the error message.
Base class for falloff rate calculators.
Definition Falloff.h:83
void setParameters(const AnyMap &node, const UnitStack &rate_units) override
Set parameters.
Definition Falloff.cpp:128
void setHighRate(const ArrheniusRate &high)
Set reaction rate in the high-pressure limit.
Definition Falloff.cpp:100
virtual void setFalloffCoeffs(const vector< double > &c)
Set coefficients of the falloff parameterization.
Definition Falloff.cpp:113
ArrheniusRate m_highRate
The reaction rate in the high-pressure limit.
Definition Falloff.h:222
void validate(const string &equation, const Kinetics &kin) override
Validate the reaction rate expression.
Definition Falloff.cpp:193
ArrheniusRate m_lowRate
The reaction rate in the low-pressure limit.
Definition Falloff.h:221
bool m_chemicallyActivated
Flag labeling reaction as chemically activated.
Definition Falloff.h:225
void getParameters(AnyMap &node) const override
Get parameters.
Definition Falloff.cpp:161
bool m_negativeA_ok
Flag indicating whether negative A values are permitted.
Definition Falloff.h:227
void check(const string &equation) override
Check basic syntax and settings of reaction rate expression.
Definition Falloff.cpp:178
virtual void getFalloffCoeffs(vector< double > &c) const
Retrieve coefficients of the falloff parameterization.
Definition Falloff.cpp:123
void setLowRate(const ArrheniusRate &low)
Set reaction rate in the low-pressure limit.
Definition Falloff.cpp:87
Error thrown for problems processing information contained in an AnyMap or AnyValue.
Definition AnyMap.h:738
Public interface for kinetics managers.
Definition Kinetics.h:125
virtual const vector< double > & thirdBodyConcentrations() const
Provide direct access to current third-body concentration values.
Definition Kinetics.h:528
The Lindemann falloff parameterization.
Definition Falloff.h:242
virtual double molarDensity() const
Molar density (kmol/m^3).
Definition Phase.cpp:576
double temperature() const
Temperature (K).
Definition Phase.h:562
int stateMFNumber() const
Return the State Mole Fraction Number.
Definition Phase.h:761
virtual void setParameters(const AnyMap &node, const UnitStack &units)
Set parameters.
bool valid() const
Get flag indicating whether reaction rate is set up correctly.
bool m_valid
Flag indicating whether reaction rate is set up correctly.
AnyMap m_input
Input data used for specific models.
The SRI falloff function.
Definition Falloff.h:398
double F(double pr, const double *work) const override
The falloff function.
Definition Falloff.cpp:400
void setParameters(const AnyMap &node, const UnitStack &rate_units) override
Set parameters.
Definition Falloff.cpp:407
void setFalloffCoeffs(const vector< double > &c) override
Set coefficients used by parameterization.
Definition Falloff.cpp:347
void updateTemp(double T, double *work) const override
Update the temperature parameters in the representation.
Definition Falloff.cpp:391
double m_d
parameter d in the 5-parameter SRI falloff function. Dimensionless.
Definition Falloff.h:461
void getParameters(AnyMap &node) const override
Get parameters.
Definition Falloff.cpp:432
double m_a
parameter a in the 5-parameter SRI falloff function. Dimensionless.
Definition Falloff.h:452
void getFalloffCoeffs(vector< double > &c) const override
Retrieve coefficients of the falloff parameterization.
Definition Falloff.cpp:377
double m_c
parameter c in the 5-parameter SRI falloff function. [K]
Definition Falloff.h:458
double m_b
parameter b in the 5-parameter SRI falloff function. [K]
Definition Falloff.h:455
SriRate()
Constructor.
Definition Falloff.h:401
double m_e
parameter d in the 5-parameter SRI falloff function. Dimensionless.
Definition Falloff.h:464
Base class for a phase with thermodynamic properties.
double F(double pr, const double *work) const override
The falloff function.
Definition Falloff.cpp:292
void setParameters(const AnyMap &node, const UnitStack &rate_units) override
Set parameters.
Definition Falloff.cpp:302
void setFalloffCoeffs(const vector< double > &c) override
Set coefficients used by parameterization.
Definition Falloff.cpp:233
void updateTemp(double T, double *work) const override
Update the temperature parameters in the representation.
Definition Falloff.cpp:283
double m_t2
parameter T_2 in the 4-parameter Troe falloff function. [K]
Definition Falloff.h:361
void getParameters(AnyMap &node) const override
Get parameters.
Definition Falloff.cpp:324
double m_rt1
parameter 1/T_1 in the 4-parameter Troe falloff function. [K^-1]
Definition Falloff.h:358
double m_a
parameter a in the 4-parameter Troe falloff function. Dimensionless
Definition Falloff.h:352
void getFalloffCoeffs(vector< double > &c) const override
Retrieve coefficients of the falloff parameterization.
Definition Falloff.cpp:270
TroeRate()
Constructor.
Definition Falloff.h:307
double m_rt3
parameter 1/T_3 in the 4-parameter Troe falloff function. [K^-1]
Definition Falloff.h:355
The 1- or 2-parameter Tsang falloff parameterization.
Definition Falloff.h:492
double F(double pr, const double *work) const override
The falloff function.
Definition Falloff.cpp:491
void setParameters(const AnyMap &node, const UnitStack &rate_units) override
Set parameters.
Definition Falloff.cpp:501
void setFalloffCoeffs(const vector< double > &c) override
Set coefficients used by parameterization.
Definition Falloff.cpp:456
void updateTemp(double T, double *work) const override
Update the temperature parameters in the representation.
Definition Falloff.cpp:485
TsangRate()
Constructor.
Definition Falloff.h:495
void getParameters(AnyMap &node) const override
Get parameters.
Definition Falloff.cpp:519
double m_a
parameter a in the Tsang F_cent formulation. Dimensionless
Definition Falloff.h:546
void getFalloffCoeffs(vector< double > &c) const override
Retrieve coefficients of the falloff parameterization.
Definition Falloff.cpp:474
double m_b
parameter b in the Tsang F_cent formulation. [K^-1]
Definition Falloff.h:549
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
void warn_user(const string &method, const string &msg, const Args &... args)
Print a user warning raised from method as CanteraWarning.
Definition global.h:267
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
const double SmallNumber
smallest number to compare to zero.
Definition ct_defs.h:158
Contains declarations for string manipulation functions within Cantera.
int m_state_mf_number
integer that is incremented when composition changes
Definition Falloff.h:58
vector< double > m_conc_3b_buf
buffered third-body concentrations
Definition Falloff.h:61
void perturbThirdBodies(double deltaM)
Perturb third-body concentration vector of data container.
Definition Falloff.cpp:57
bool update(const ThermoPhase &phase, const Kinetics &kin) override
Update data container based on thermodynamic phase state.
Definition Falloff.cpp:38
vector< double > conc_3b
vector of effective third-body concentrations
Definition Falloff.h:54
double molar_density
used to determine if updates are needed
Definition Falloff.h:53
void restore() override
Restore data container after a perturbation.
Definition Falloff.cpp:70
bool m_perturbed
boolean indicating whether 3-rd body values are perturbed
Definition Falloff.h:60
virtual void update(double T)
Update data container based on temperature T
double temperature
temperature
virtual void restore()
Restore data container after a perturbation.
Unit aggregation utility.
Definition Units.h:105
size_t size() const
Size of UnitStack.
Definition Units.h:118
void join(double exponent)
Join (update) exponent of standard units, where the updated exponent is the sum of the pre-existing e...
Definition Units.cpp:352