1 // This file is part of Cantera. See License.txt in the top-level directory or
2 // at https://cantera.org/license.txt for license and copyright information.
3
4 #include "cantera/kinetics/BulkKinetics.h"
5 #include "cantera/kinetics/Reaction.h"
6 #include "cantera/thermo/ThermoPhase.h"
7
8 namespace Cantera
9 {
10
BulkKinetics(ThermoPhase * thermo)11 BulkKinetics::BulkKinetics(ThermoPhase* thermo) :
12 m_ROP_ok(false),
13 m_temp(0.0)
14 {
15 if (thermo) {
16 addPhase(*thermo);
17 }
18 }
19
isReversible(size_t i)20 bool BulkKinetics::isReversible(size_t i) {
21 return std::find(m_revindex.begin(), m_revindex.end(), i) < m_revindex.end();
22 }
23
getDeltaGibbs(doublereal * deltaG)24 void BulkKinetics::getDeltaGibbs(doublereal* deltaG)
25 {
26 // Get the chemical potentials of the species in the ideal gas solution.
27 thermo().getChemPotentials(m_grt.data());
28 // Use the stoichiometric manager to find deltaG for each reaction.
29 getReactionDelta(m_grt.data(), deltaG);
30 }
31
getDeltaEnthalpy(doublereal * deltaH)32 void BulkKinetics::getDeltaEnthalpy(doublereal* deltaH)
33 {
34 // Get the partial molar enthalpy of all species in the ideal gas.
35 thermo().getPartialMolarEnthalpies(m_grt.data());
36 // Use the stoichiometric manager to find deltaH for each reaction.
37 getReactionDelta(m_grt.data(), deltaH);
38 }
39
getDeltaEntropy(doublereal * deltaS)40 void BulkKinetics::getDeltaEntropy(doublereal* deltaS)
41 {
42 // Get the partial molar entropy of all species in the solid solution.
43 thermo().getPartialMolarEntropies(m_grt.data());
44 // Use the stoichiometric manager to find deltaS for each reaction.
45 getReactionDelta(m_grt.data(), deltaS);
46 }
47
getDeltaSSGibbs(doublereal * deltaG)48 void BulkKinetics::getDeltaSSGibbs(doublereal* deltaG)
49 {
50 // Get the standard state chemical potentials of the species. This is the
51 // array of chemical potentials at unit activity. We define these here as
52 // the chemical potentials of the pure species at the temperature and
53 // pressure of the solution.
54 thermo().getStandardChemPotentials(m_grt.data());
55 // Use the stoichiometric manager to find deltaG for each reaction.
56 getReactionDelta(m_grt.data(), deltaG);
57 }
58
getDeltaSSEnthalpy(doublereal * deltaH)59 void BulkKinetics::getDeltaSSEnthalpy(doublereal* deltaH)
60 {
61 // Get the standard state enthalpies of the species.
62 thermo().getEnthalpy_RT(m_grt.data());
63 for (size_t k = 0; k < m_kk; k++) {
64 m_grt[k] *= thermo().RT();
65 }
66 // Use the stoichiometric manager to find deltaH for each reaction.
67 getReactionDelta(m_grt.data(), deltaH);
68 }
69
getDeltaSSEntropy(doublereal * deltaS)70 void BulkKinetics::getDeltaSSEntropy(doublereal* deltaS)
71 {
72 // Get the standard state entropy of the species. We define these here as
73 // the entropies of the pure species at the temperature and pressure of the
74 // solution.
75 thermo().getEntropy_R(m_grt.data());
76 for (size_t k = 0; k < m_kk; k++) {
77 m_grt[k] *= GasConstant;
78 }
79 // Use the stoichiometric manager to find deltaS for each reaction.
80 getReactionDelta(m_grt.data(), deltaS);
81 }
82
getRevRateConstants(double * krev,bool doIrreversible)83 void BulkKinetics::getRevRateConstants(double* krev, bool doIrreversible)
84 {
85 // go get the forward rate constants. -> note, we don't really care about
86 // speed or redundancy in these informational routines.
87 getFwdRateConstants(krev);
88
89 if (doIrreversible) {
90 getEquilibriumConstants(m_ropnet.data());
91 for (size_t i = 0; i < nReactions(); i++) {
92 krev[i] /= m_ropnet[i];
93 }
94 } else {
95 // m_rkcn[] is zero for irreversible reactions
96 for (size_t i = 0; i < nReactions(); i++) {
97 krev[i] *= m_rkcn[i];
98 }
99 }
100 }
101
addReaction(shared_ptr<Reaction> r,bool resize)102 bool BulkKinetics::addReaction(shared_ptr<Reaction> r, bool resize)
103 {
104 bool added = Kinetics::addReaction(r, resize);
105 if (!added) {
106 // undeclared species, etc.
107 return false;
108 }
109 double dn = 0.0;
110 for (const auto& sp : r->products) {
111 dn += sp.second;
112 }
113 for (const auto& sp : r->reactants) {
114 dn -= sp.second;
115 }
116
117 m_dn.push_back(dn);
118
119 if (r->reversible) {
120 m_revindex.push_back(nReactions()-1);
121 } else {
122 m_irrev.push_back(nReactions()-1);
123 }
124
125 if (!(r->usesLegacy())) {
126 shared_ptr<ReactionRateBase> rate = r->rate();
127 // If neccessary, add new MultiBulkRate evaluator
128 if (m_bulk_types.find(rate->type()) == m_bulk_types.end()) {
129 m_bulk_types[rate->type()] = m_bulk_rates.size();
130 m_bulk_rates.push_back(rate->newMultiRate());
131 m_bulk_rates.back()->resizeSpecies(m_kk);
132 }
133
134 // Add reaction rate to evaluator
135 size_t index = m_bulk_types[rate->type()];
136 m_bulk_rates[index]->add(nReactions() - 1, *rate);
137
138 // Add reaction to third-body evaluator
139 if (r->thirdBody() != nullptr) {
140 addThirdBody(r);
141 }
142 }
143
144 m_concm.push_back(0.0);
145
146 return true;
147 }
148
addThirdBody(shared_ptr<Reaction> r)149 void BulkKinetics::addThirdBody(shared_ptr<Reaction> r)
150 {
151 std::map<size_t, double> efficiencies;
152 for (const auto& eff : r->thirdBody()->efficiencies) {
153 size_t k = kineticsSpeciesIndex(eff.first);
154 if (k != npos) {
155 efficiencies[k] = eff.second;
156 } else if (!m_skipUndeclaredThirdBodies) {
157 throw CanteraError("BulkKinetics::addThirdBody", "Found "
158 "third-body efficiency for undefined species '" + eff.first +
159 "' while adding reaction '" + r->equation() + "'");
160 }
161 }
162 m_multi_concm.install(nReactions() - 1, efficiencies,
163 r->thirdBody()->default_efficiency);
164 concm_multi_values.resize(m_multi_concm.workSize());
165 m_multi_indices.push_back(nReactions() - 1);
166 }
167
addElementaryReaction(ElementaryReaction2 & r)168 void BulkKinetics::addElementaryReaction(ElementaryReaction2& r)
169 {
170 m_rates.install(nReactions()-1, r.rate);
171 }
172
modifyReaction(size_t i,shared_ptr<Reaction> rNew)173 void BulkKinetics::modifyReaction(size_t i, shared_ptr<Reaction> rNew)
174 {
175 // operations common to all reaction types
176 Kinetics::modifyReaction(i, rNew);
177
178 if (!(rNew->usesLegacy())) {
179 shared_ptr<ReactionRateBase> rate = rNew->rate();
180 // Ensure that MultiBulkRate evaluator is available
181 if (m_bulk_types.find(rate->type()) == m_bulk_types.end()) {
182 throw CanteraError("BulkKinetics::modifyReaction",
183 "Evaluator not available for type '{}'.", rate->type());
184 }
185
186 // Replace reaction rate to evaluator
187 size_t index = m_bulk_types[rate->type()];
188 m_bulk_rates[index]->replace(i, *rate);
189 }
190 }
191
modifyElementaryReaction(size_t i,ElementaryReaction2 & rNew)192 void BulkKinetics::modifyElementaryReaction(size_t i, ElementaryReaction2& rNew)
193 {
194 m_rates.replace(i, rNew.rate);
195 }
196
resizeSpecies()197 void BulkKinetics::resizeSpecies()
198 {
199 Kinetics::resizeSpecies();
200 m_act_conc.resize(m_kk);
201 m_phys_conc.resize(m_kk);
202 m_grt.resize(m_kk);
203 for (auto& rates : m_bulk_rates) {
204 rates->resizeSpecies(m_kk);
205 }
206 }
207
setMultiplier(size_t i,double f)208 void BulkKinetics::setMultiplier(size_t i, double f) {
209 Kinetics::setMultiplier(i, f);
210 m_ROP_ok = false;
211 }
212
invalidateCache()213 void BulkKinetics::invalidateCache()
214 {
215 Kinetics::invalidateCache();
216 m_ROP_ok = false;
217 m_temp += 0.13579;
218 }
219
220 }
221