Cantera  3.1.0a1
PlogRate.cpp
Go to the documentation of this file.
1 //! @file PlogRate.cpp
2 
3 // This file is part of Cantera. See License.txt in the top-level directory or
4 // at https://cantera.org/license.txt for license and copyright information.
5 
8 
9 namespace Cantera
10 {
11 
12 void PlogData::update(double T)
13 {
14  throw CanteraError("PlogData::update",
15  "Missing state information: 'PlogData' requires pressure.");
16 }
17 
18 bool PlogData::update(const ThermoPhase& phase, const Kinetics& kin)
19 {
20  double T = phase.temperature();
21  double P = phase.pressure();
22  if (P != pressure || T != temperature) {
23  update(T, P);
24  return true;
25  }
26  return false;
27 }
28 
29 void PlogData::perturbPressure(double deltaP)
30 {
31  if (m_pressure_buf > 0.) {
32  throw CanteraError("PlogData::perturbPressure",
33  "Cannot apply another perturbation as state is already perturbed.");
34  }
36  update(temperature, pressure * (1. + deltaP));
37 }
38 
40 {
42  // only restore if there is a valid buffered value
43  if (m_pressure_buf < 0.) {
44  return;
45  }
47  m_pressure_buf = -1.;
48 }
49 
50 // Methods of class PlogRate
51 
52 PlogRate::PlogRate(const std::multimap<double, ArrheniusRate>& rates)
53 {
54  setRates(rates);
55 }
56 
57 PlogRate::PlogRate(const AnyMap& node, const UnitStack& rate_units)
58 {
59  setParameters(node, rate_units);
60 }
61 
62 void PlogRate::setParameters(const AnyMap& node, const UnitStack& rate_units)
63 {
64  ReactionRate::setParameters(node, rate_units);
65  std::multimap<double, ArrheniusRate> multi_rates;
66  if (node.hasKey("rate-constants")) {
67  auto& rates = node["rate-constants"].asVector<AnyMap>();
68  for (const auto& rate : rates) {
69  multi_rates.insert({rate.convert("P", "Pa"),
70  ArrheniusRate(AnyValue(rate), node.units(), rate_units)});
71  }
72  }
73  setRates(multi_rates);
74 }
75 
76 void PlogRate::getParameters(AnyMap& rateNode, const Units& rate_units) const
77 {
78  vector<AnyMap> rateList;
79  if (!valid()) {
80  // object not fully set up
81  return;
82  }
83  for (const auto& [pressure, rate] : getRates()) {
84  AnyMap rateNode_;
85  rateNode_["P"].setQuantity(pressure, "Pa");
86  rate.getRateParameters(rateNode_);
87  rateList.push_back(std::move(rateNode_));
88  }
89  rateNode["rate-constants"] = std::move(rateList);
90 }
91 
92 void PlogRate::setRates(const std::multimap<double, ArrheniusRate>& rates)
93 {
94  size_t j = 0;
95  rates_.clear();
96  pressures_.clear();
97  m_valid = !rates.empty();
98  rates_.reserve(rates.size());
99  // Insert intermediate pressures
100  for (const auto& [pressure, rate] : rates) {
101  double logp = std::log(pressure);
102  if (pressures_.empty() || pressures_.rbegin()->first != logp) {
103  // starting a new group
104  pressures_[logp] = {j, j+1};
105  } else {
106  // another rate expression at the same pressure
107  pressures_[logp].second = j+1;
108  }
109 
110  j++;
111  rates_.push_back(rate);
112  }
113  if (!m_valid) {
114  // ensure that reaction rate can be evaluated (but returns NaN)
115  rates_.reserve(1);
116  pressures_[std::log(OneBar)] = {0, 0};
117  rates_.push_back(ArrheniusRate());
118  }
119 
120  // Duplicate the first and last groups to handle P < P_0 and P > P_N
121  pressures_.insert({-1000.0, pressures_.begin()->second});
122  pressures_.insert({1000.0, pressures_.rbegin()->second});
123 }
124 
125 void PlogRate::validate(const string& equation, const Kinetics& kin)
126 {
127  if (!valid()) {
128  throw InputFileError("PlogRate::validate", m_input,
129  "Rate object for reaction '{}' is not configured.", equation);
130  }
131 
132  fmt::memory_buffer err_reactions;
133  double T[] = {300.0, 500.0, 1000.0, 2000.0, 5000.0, 10000.0};
134  PlogData data;
135 
136  for (auto iter = ++pressures_.begin(); iter->first < 1000; iter++) {
137  data.update(T[0], exp(iter->first)); // iter->first contains log(p)
138  updateFromStruct(data);
139  for (size_t i=0; i < 6; i++) {
140  double k = 0;
141  for (size_t p = ilow1_; p < ilow2_; p++) {
142  k += rates_.at(p).evalRate(log(T[i]), 1.0 / T[i]);
143  }
144  if (!(k > 0)) {
145  fmt_append(err_reactions,
146  "at P = {:.5g}, T = {:.1f}\n", std::exp(iter->first), T[i]);
147  }
148  }
149  }
150  if (err_reactions.size()) {
151  throw InputFileError("PlogRate::validate", m_input,
152  "\nInvalid rate coefficient for reaction '{}'\n{}",
153  equation, to_string(err_reactions));
154  }
155 }
156 
157 std::multimap<double, ArrheniusRate> PlogRate::getRates() const
158 {
159  std::multimap<double, ArrheniusRate> rateMap;
160  // initial preincrement to skip rate for P --> 0
161  for (auto iter = ++pressures_.begin();
162  iter->first < 1000; // skip rates for (P --> infinity)
163  ++iter) {
164  for (size_t i = iter->second.first; i < iter->second.second; i++) {
165  rateMap.insert({std::exp(iter->first), rates_[i]});
166  }
167  }
168  return rateMap;
169 }
170 
171 }
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
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
A wrapper for a variable whose type is determined at runtime.
Definition: AnyMap.h:86
Arrhenius reaction rate type depends only on temperature.
Definition: Arrhenius.h:170
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:66
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
double temperature() const
Temperature (K).
Definition: Phase.h:562
virtual double pressure() const
Return the thermodynamic pressure (Pa).
Definition: Phase.h:580
map< double, pair< size_t, size_t > > pressures_
log(p) to (index range) in the rates_ vector
Definition: PlogRate.h:181
void updateFromStruct(const PlogData &shared_data)
Update information specific to reaction.
Definition: PlogRate.h:110
PlogRate()=default
Default constructor.
std::multimap< double, ArrheniusRate > getRates() const
Return the pressures and Arrhenius expressions which comprise this reaction.
Definition: PlogRate.cpp:157
void setParameters(const AnyMap &node, const UnitStack &rate_units) override
Perform object setup based on AnyMap node information.
Definition: PlogRate.cpp:62
void setRates(const std::multimap< double, ArrheniusRate > &rates)
Set up Plog object.
Definition: PlogRate.cpp:92
void validate(const string &equation, const Kinetics &kin) override
Check to make sure that the rate expression is finite over a range of temperatures at each interpolat...
Definition: PlogRate.cpp:125
size_t ilow1_
Indices to the ranges within rates_ for the lower / upper pressure, such that rates_[ilow1_] through ...
Definition: PlogRate.h:194
virtual void setParameters(const AnyMap &node, const UnitStack &units)
Set parameters.
Definition: ReactionRate.h:101
bool valid() const
Get flag indicating whether reaction rate is set up correctly.
Definition: ReactionRate.h:203
bool m_valid
Flag indicating whether reaction rate is set up correctly.
Definition: ReactionRate.h:236
AnyMap m_input
Input data used for specific models.
Definition: ReactionRate.h:230
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:390
A representation of the units associated with a dimensional quantity.
Definition: Units.h:35
void fmt_append(fmt::memory_buffer &b, Args... args)
Versions 6.2.0 and 6.2.1 of fmtlib do not include this define before they include windows....
Definition: fmt.h:29
const double OneBar
One bar [Pa].
Definition: ct_defs.h:99
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:564
Data container holding shared data specific to PlogRate.
Definition: PlogRate.h:20
void perturbPressure(double deltaP)
Perturb pressure of data container.
Definition: PlogRate.cpp:29
double m_pressure_buf
buffered pressure
Definition: PlogRate.h:53
void update(double T) override
Update data container based on temperature T
Definition: PlogRate.cpp:12
virtual void update(double T)
Update data container based on temperature T
Definition: ReactionData.h:36
void restore() override
Restore data container after a perturbation.
Definition: PlogRate.cpp:39
double pressure
pressure
Definition: PlogRate.h:49
double temperature
temperature
Definition: ReactionData.h:110
virtual void restore()
Restore data container after a perturbation.
Definition: ReactionData.h:90
Unit aggregation utility.
Definition: Units.h:105