Cantera  3.2.0a2
Loading...
Searching...
No Matches
Reactor.cpp
Go to the documentation of this file.
1//! @file Reactor.cpp A zero-dimensional reactor
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
16
17#include <boost/math/tools/roots.hpp>
18
19using namespace std;
20namespace bmt = boost::math::tools;
21
22namespace Cantera
23{
24
25Reactor::Reactor(shared_ptr<Solution> sol, const string& name)
26 : ReactorBase(sol, name)
27{
28 m_kin = sol->kinetics().get();
29 if (m_kin->nReactions() == 0) {
30 setChemistry(false);
31 } else {
32 setChemistry(true);
33 }
34 m_vol = 1.0; // By default, the volume is set to 1.0 m^3.
35}
36
37void Reactor::setDerivativeSettings(AnyMap& settings)
38{
39 m_kin->setDerivativeSettings(settings);
40 // translate settings to surfaces
41 for (auto S : m_surfaces) {
42 S->kinetics()->setDerivativeSettings(settings);
43 }
44}
45
46void Reactor::setKinetics(Kinetics& kin)
47{
48 warn_deprecated("Reactor::setKinetics",
49 "After Cantera 3.2, a change of reactor contents after instantiation "
50 "will be disabled.");
51 m_kin = &kin;
52 if (m_kin->nReactions() == 0) {
53 setChemistry(false);
54 } else {
55 setChemistry(true);
56 }
57}
58
59void Reactor::getState(double* y)
60{
61 if (m_thermo == 0) {
62 throw CanteraError("Reactor::getState",
63 "Error: reactor is empty.");
64 }
65 m_thermo->restoreState(m_state);
66
67 // set the first component to the total mass
68 m_mass = m_thermo->density() * m_vol;
69 y[0] = m_mass;
70
71 // set the second component to the total volume
72 y[1] = m_vol;
73
74 // set the third component to the total internal energy
75 y[2] = m_thermo->intEnergy_mass() * m_mass;
76
77 // set components y+3 ... y+K+2 to the mass fractions of each species
78 m_thermo->getMassFractions(y+3);
79
80 // set the remaining components to the surface species
81 // coverages on the walls
82 getSurfaceInitialConditions(y + m_nsp + 3);
83}
84
85void Reactor::getSurfaceInitialConditions(double* y)
86{
87 size_t loc = 0;
88 for (auto& S : m_surfaces) {
89 S->getCoverages(y + loc);
90 loc += S->thermo()->nSpecies();
91 }
92}
93
94void Reactor::initialize(double t0)
95{
96 if (!m_thermo || (m_chem && !m_kin)) {
97 throw CanteraError("Reactor::initialize", "Reactor contents not set"
98 " for reactor '" + m_name + "'.");
99 }
100 m_thermo->restoreState(m_state);
101 m_sdot.resize(m_nsp, 0.0);
102 m_wdot.resize(m_nsp, 0.0);
103 updateConnected(true);
104
105 for (size_t n = 0; n < m_wall.size(); n++) {
106 WallBase* W = m_wall[n];
107 W->initialize();
108 }
109
110 m_nv = m_nsp + 3;
111 m_nv_surf = 0;
112 size_t maxnt = 0;
113 for (auto& S : m_surfaces) {
114 m_nv_surf += S->thermo()->nSpecies();
115 size_t nt = S->kinetics()->nTotalSpecies();
116 maxnt = std::max(maxnt, nt);
117 }
118 m_nv += m_nv_surf;
119 m_work.resize(maxnt);
120}
121
122size_t Reactor::nSensParams() const
123{
124 size_t ns = m_sensParams.size();
125 for (auto& S : m_surfaces) {
126 ns += S->nSensParams();
127 }
128 return ns;
129}
130
131void Reactor::updateState(double* y)
132{
133 // The components of y are [0] the total mass, [1] the total volume,
134 // [2] the total internal energy, [3...K+3] are the mass fractions of each
135 // species, and [K+3...] are the coverages of surface species on each wall.
136 m_mass = y[0];
137 m_vol = y[1];
138 m_thermo->setMassFractions_NoNorm(y+3);
139
140 if (m_energy) {
141 double U = y[2];
142 // Residual function: error in internal energy as a function of T
143 auto u_err = [this, U](double T) {
144 m_thermo->setState_TD(T, m_mass / m_vol);
145 return m_thermo->intEnergy_mass() * m_mass - U;
146 };
147
148 double T = m_thermo->temperature();
149 boost::uintmax_t maxiter = 100;
150 pair<double, double> TT;
151 try {
152 TT = bmt::bracket_and_solve_root(
153 u_err, T, 1.2, true, bmt::eps_tolerance<double>(48), maxiter);
154 } catch (std::exception&) {
155 // Try full-range bisection if bracketing fails (for example, near
156 // temperature limits for the phase's equation of state)
157 try {
158 TT = bmt::bisect(u_err, m_thermo->minTemp(), m_thermo->maxTemp(),
159 bmt::eps_tolerance<double>(48), maxiter);
160 } catch (std::exception& err2) {
161 // Set m_thermo back to a reasonable state if root finding fails
162 m_thermo->setState_TD(T, m_mass / m_vol);
163 throw CanteraError("Reactor::updateState",
164 "{}\nat U = {}, rho = {}", err2.what(), U, m_mass / m_vol);
165 }
166 }
167 if (fabs(TT.first - TT.second) > 1e-7*TT.first) {
168 throw CanteraError("Reactor::updateState", "root finding failed");
169 }
170 m_thermo->setState_TD(TT.second, m_mass / m_vol);
171 } else {
172 m_thermo->setDensity(m_mass/m_vol);
173 }
174
175 updateConnected(true);
176 updateSurfaceState(y + m_nsp + 3);
177}
178
179void Reactor::updateSurfaceState(double* y)
180{
181 size_t loc = 0;
182 for (auto& S : m_surfaces) {
183 S->setCoverages(y+loc);
184 loc += S->thermo()->nSpecies();
185 }
186}
187
188void Reactor::updateConnected(bool updatePressure) {
189 // save parameters needed by other connected reactors
190 m_enthalpy = m_thermo->enthalpy_mass();
191 if (updatePressure) {
192 m_pressure = m_thermo->pressure();
193 }
194 try {
195 m_intEnergy = m_thermo->intEnergy_mass();
196 } catch (NotImplementedError&) {
197 m_intEnergy = NAN;
198 }
199 m_thermo->saveState(m_state);
200
201 // Update the mass flow rate of connected flow devices
202 double time = 0.0;
203 if (m_net != nullptr) {
204 time = (timeIsIndependent()) ? m_net->time() : m_net->distance();
205 }
206 for (size_t i = 0; i < m_outlet.size(); i++) {
207 m_outlet[i]->setSimTime(time);
208 m_outlet[i]->updateMassFlowRate(time);
209 }
210 for (size_t i = 0; i < m_inlet.size(); i++) {
211 m_inlet[i]->setSimTime(time);
212 m_inlet[i]->updateMassFlowRate(time);
213 }
214 for (size_t i = 0; i < m_wall.size(); i++) {
215 m_wall[i]->setSimTime(time);
216 }
217}
218
219void Reactor::eval(double time, double* LHS, double* RHS)
220{
221 double& dmdt = RHS[0];
222 double* mdYdt = RHS + 3; // mass * dY/dt
223
224 evalWalls(time);
225 m_thermo->restoreState(m_state);
226 const vector<double>& mw = m_thermo->molecularWeights();
227 const double* Y = m_thermo->massFractions();
228
229 evalSurfaces(LHS + m_nsp + 3, RHS + m_nsp + 3, m_sdot.data());
230 // mass added to gas phase from surface reactions
231 double mdot_surf = dot(m_sdot.begin(), m_sdot.end(), mw.begin());
232 dmdt = mdot_surf;
233
234 // volume equation
235 RHS[1] = m_vdot;
236
237 if (m_chem) {
238 m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
239 }
240
241 for (size_t k = 0; k < m_nsp; k++) {
242 // production in gas phase and from surfaces
243 mdYdt[k] = (m_wdot[k] * m_vol + m_sdot[k]) * mw[k];
244 // dilution by net surface mass flux
245 mdYdt[k] -= Y[k] * mdot_surf;
246 LHS[k+3] = m_mass;
247 }
248
249 // Energy equation.
250 // @f[
251 // \dot U = -P\dot V + A \dot q + \dot m_{in} h_{in} - \dot m_{out} h.
252 // @f]
253 if (m_energy) {
254 RHS[2] = - m_thermo->pressure() * m_vdot + m_Qdot;
255 } else {
256 RHS[2] = 0.0;
257 }
258
259 // add terms for outlets
260 for (auto outlet : m_outlet) {
261 double mdot = outlet->massFlowRate();
262 dmdt -= mdot; // mass flow out of system
263 if (m_energy) {
264 RHS[2] -= mdot * m_enthalpy;
265 }
266 }
267
268 // add terms for inlets
269 for (auto inlet : m_inlet) {
270 double mdot = inlet->massFlowRate();
271 dmdt += mdot; // mass flow into system
272 for (size_t n = 0; n < m_nsp; n++) {
273 double mdot_spec = inlet->outletSpeciesMassFlowRate(n);
274 // flow of species into system and dilution by other species
275 mdYdt[n] += (mdot_spec - mdot * Y[n]);
276 }
277 if (m_energy) {
278 RHS[2] += mdot * inlet->enthalpy_mass();
279 }
280 }
281}
282
283void Reactor::evalWalls(double t)
284{
285 // time is currently unused
286 m_vdot = 0.0;
287 m_Qdot = 0.0;
288 for (size_t i = 0; i < m_wall.size(); i++) {
289 int f = 2 * m_lr[i] - 1;
290 m_vdot -= f * m_wall[i]->expansionRate();
291 m_Qdot += f * m_wall[i]->heatRate();
292 }
293}
294
295void Reactor::evalSurfaces(double* LHS, double* RHS, double* sdot)
296{
297 fill(sdot, sdot + m_nsp, 0.0);
298 size_t loc = 0; // offset into ydot
299
300 for (auto S : m_surfaces) {
301 Kinetics* kin = S->kinetics();
302 SurfPhase* surf = S->thermo();
303
304 double rs0 = 1.0/surf->siteDensity();
305 size_t nk = surf->nSpecies();
306 double sum = 0.0;
307 S->syncState();
308 kin->getNetProductionRates(&m_work[0]);
309 for (size_t k = 1; k < nk; k++) {
310 RHS[loc + k] = m_work[k] * rs0 * surf->size(k);
311 sum -= RHS[loc + k];
312 }
313 RHS[loc] = sum;
314 loc += nk;
315
316 size_t bulkloc = kin->kineticsSpeciesIndex(m_thermo->speciesName(0));
317 double wallarea = S->area();
318 for (size_t k = 0; k < m_nsp; k++) {
319 sdot[k] += m_work[bulkloc + k] * wallarea;
320 }
321 }
322}
323
324vector<size_t> Reactor::steadyConstraints() const
325{
326 if (!energyEnabled()) {
327 throw CanteraError("Reactor::steadyConstraints", "Steady state solver cannot"
328 " be used with {0} when energy equation is disabled."
329 "\nConsider using IdealGas{0} instead.\n"
330 "See https://github.com/Cantera/enhancements/issues/234", type());
331 }
332 if (nSurfs() != 0) {
333 throw CanteraError("Reactor::steadyConstraints", "Steady state solver cannot"
334 " currently be used when reactor surfaces are present.\n"
335 "See https://github.com/Cantera/enhancements/issues/234.");
336 }
337 return {1}; // volume
338}
339
340Eigen::SparseMatrix<double> Reactor::finiteDifferenceJacobian()
341{
342 if (m_nv == 0) {
343 throw CanteraError("Reactor::finiteDifferenceJacobian",
344 "Reactor must be initialized first.");
345 }
346 // clear former jacobian elements
347 m_jac_trips.clear();
348
349 Eigen::ArrayXd yCurrent(m_nv);
350 getState(yCurrent.data());
351 double time = (m_net != nullptr) ? m_net->time() : 0.0;
352
353 Eigen::ArrayXd yPerturbed = yCurrent;
354 Eigen::ArrayXd lhsPerturbed(m_nv), lhsCurrent(m_nv);
355 Eigen::ArrayXd rhsPerturbed(m_nv), rhsCurrent(m_nv);
356 lhsCurrent = 1.0;
357 rhsCurrent = 0.0;
358 updateState(yCurrent.data());
359 eval(time, lhsCurrent.data(), rhsCurrent.data());
360
361 double rel_perturb = std::sqrt(std::numeric_limits<double>::epsilon());
362 double atol = (m_net != nullptr) ? m_net->atol() : 1e-15;
363
364 for (size_t j = 0; j < m_nv; j++) {
365 yPerturbed = yCurrent;
366 double delta_y = std::max(std::abs(yCurrent[j]), 1000 * atol) * rel_perturb;
367 yPerturbed[j] += delta_y;
368
369 updateState(yPerturbed.data());
370 lhsPerturbed = 1.0;
371 rhsPerturbed = 0.0;
372 eval(time, lhsPerturbed.data(), rhsPerturbed.data());
373
374 // d ydot_i/dy_j
375 for (size_t i = 0; i < m_nv; i++) {
376 double ydotPerturbed = rhsPerturbed[i] / lhsPerturbed[i];
377 double ydotCurrent = rhsCurrent[i] / lhsCurrent[i];
378 if (ydotCurrent != ydotPerturbed) {
379 m_jac_trips.emplace_back(
380 static_cast<int>(i), static_cast<int>(j),
381 (ydotPerturbed - ydotCurrent) / delta_y);
382 }
383 }
384 }
385 updateState(yCurrent.data());
386
387 Eigen::SparseMatrix<double> jac(m_nv, m_nv);
388 jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
389 return jac;
390}
391
392
393void Reactor::evalSurfaces(double* RHS, double* sdot)
394{
395 fill(sdot, sdot + m_nsp, 0.0);
396 size_t loc = 0; // offset into ydot
397
398 for (auto S : m_surfaces) {
399 Kinetics* kin = S->kinetics();
400 SurfPhase* surf = S->thermo();
401
402 double rs0 = 1.0/surf->siteDensity();
403 size_t nk = surf->nSpecies();
404 double sum = 0.0;
405 S->syncState();
406 kin->getNetProductionRates(&m_work[0]);
407 for (size_t k = 1; k < nk; k++) {
408 RHS[loc + k] = m_work[k] * rs0 * surf->size(k);
409 sum -= RHS[loc + k];
410 }
411 RHS[loc] = sum;
412 loc += nk;
413
414 size_t bulkloc = kin->kineticsSpeciesIndex(m_thermo->speciesName(0));
415 double wallarea = S->area();
416 for (size_t k = 0; k < m_nsp; k++) {
417 sdot[k] += m_work[bulkloc + k] * wallarea;
418 }
419 }
420}
421
422void Reactor::addSensitivityReaction(size_t rxn)
423{
424 if (!m_chem || rxn >= m_kin->nReactions()) {
425 throw CanteraError("Reactor::addSensitivityReaction",
426 "Reaction number out of range ({})", rxn);
427 }
428
429 size_t p = network().registerSensitivityParameter(
430 name()+": "+m_kin->reaction(rxn)->equation(), 1.0, 1.0);
431 m_sensParams.emplace_back(
432 SensitivityParameter{rxn, p, 1.0, SensParameterType::reaction});
433}
434
435void Reactor::addSensitivitySpeciesEnthalpy(size_t k)
436{
437 if (k >= m_thermo->nSpecies()) {
438 throw CanteraError("Reactor::addSensitivitySpeciesEnthalpy",
439 "Species index out of range ({})", k);
440 }
441
442 size_t p = network().registerSensitivityParameter(
443 name() + ": " + m_thermo->speciesName(k) + " enthalpy",
444 0.0, GasConstant * 298.15);
445 m_sensParams.emplace_back(
446 SensitivityParameter{k, p, m_thermo->Hf298SS(k),
447 SensParameterType::enthalpy});
448}
449
450size_t Reactor::speciesIndex(const string& nm) const
451{
452 // check for a gas species name
453 size_t k = m_thermo->speciesIndex(nm);
454 if (k != npos) {
455 return k;
456 }
457
458 // check for a wall species
459 size_t offset = m_nsp;
460 for (auto& S : m_surfaces) {
461 ThermoPhase* th = S->thermo();
462 k = th->speciesIndex(nm);
463 if (k != npos) {
464 return k + offset;
465 } else {
466 offset += th->nSpecies();
467 }
468 }
469 return npos;
470}
471
472size_t Reactor::componentIndex(const string& nm) const
473{
474 size_t k = speciesIndex(nm);
475 if (k != npos) {
476 return k + 3;
477 } else if (nm == "mass") {
478 return 0;
479 } else if (nm == "volume") {
480 return 1;
481 } else if (nm == "int_energy") {
482 return 2;
483 } else {
484 return npos;
485 }
486}
487
488string Reactor::componentName(size_t k) {
489 if (k == 0) {
490 return "mass";
491 } else if (k == 1) {
492 return "volume";
493 } else if (k == 2) {
494 return "int_energy";
495 } else if (k >= 3 && k < neq()) {
496 k -= 3;
497 if (k < m_thermo->nSpecies()) {
498 return m_thermo->speciesName(k);
499 } else {
500 k -= m_thermo->nSpecies();
501 }
502 for (auto& S : m_surfaces) {
503 ThermoPhase* th = S->thermo();
504 if (k < th->nSpecies()) {
505 return th->speciesName(k);
506 } else {
507 k -= th->nSpecies();
508 }
509 }
510 }
511 throw CanteraError("Reactor::componentName", "Index is out of bounds.");
512}
513
514double Reactor::upperBound(size_t k) const {
515 if (k == 0) {
516 return BigNumber; // mass
517 } else if (k == 1) {
518 return BigNumber; // volume
519 } else if (k == 2) {
520 return BigNumber; // internal energy
521 } else if (k >= 3 && k < m_nv) {
522 return 1.0; // species mass fraction or surface coverage
523 } else {
524 throw CanteraError("Reactor::upperBound", "Index {} is out of bounds.", k);
525 }
526}
527
528double Reactor::lowerBound(size_t k) const {
529 if (k == 0) {
530 return 0; // mass
531 } else if (k == 1) {
532 return 0; // volume
533 } else if (k == 2) {
534 return -BigNumber; // internal energy
535 } else if (k >= 3 && k < m_nv) {
536 return -Tiny; // species mass fraction or surface coverage
537 } else {
538 throw CanteraError("Reactor::lowerBound", "Index {} is out of bounds.", k);
539 }
540}
541
542void Reactor::resetBadValues(double* y) {
543 for (size_t k = 3; k < m_nv; k++) {
544 y[k] = std::max(y[k], 0.0);
545 }
546}
547
548void Reactor::applySensitivity(double* params)
549{
550 if (!params) {
551 return;
552 }
553 for (auto& p : m_sensParams) {
554 if (p.type == SensParameterType::reaction) {
555 p.value = m_kin->multiplier(p.local);
556 m_kin->setMultiplier(p.local, p.value*params[p.global]);
557 } else if (p.type == SensParameterType::enthalpy) {
558 m_thermo->modifyOneHf298SS(p.local, p.value + params[p.global]);
559 }
560 }
561 for (auto& S : m_surfaces) {
562 S->setSensitivityParameters(params);
563 }
564 m_thermo->invalidateCache();
565 if (m_kin) {
566 m_kin->invalidateCache();
567 }
568}
569
570void Reactor::resetSensitivity(double* params)
571{
572 if (!params) {
573 return;
574 }
575 for (auto& p : m_sensParams) {
576 if (p.type == SensParameterType::reaction) {
577 m_kin->setMultiplier(p.local, p.value);
578 } else if (p.type == SensParameterType::enthalpy) {
579 m_thermo->resetHf298(p.local);
580 }
581 }
582 for (auto& S : m_surfaces) {
583 S->resetSensitivityParameters();
584 }
585 m_thermo->invalidateCache();
586 if (m_kin) {
587 m_kin->invalidateCache();
588 }
589}
590
591void Reactor::setAdvanceLimits(const double *limits)
592{
593 if (m_thermo == 0) {
594 throw CanteraError("Reactor::setAdvanceLimits",
595 "Error: reactor is empty.");
596 }
597 m_advancelimits.assign(limits, limits + m_nv);
598
599 // resize to zero length if no limits are set
600 if (std::none_of(m_advancelimits.begin(), m_advancelimits.end(),
601 [](double val){return val>0;})) {
602 m_advancelimits.resize(0);
603 }
604}
605
606bool Reactor::getAdvanceLimits(double *limits) const
607{
608 bool has_limit = hasAdvanceLimits();
609 if (has_limit) {
610 std::copy(m_advancelimits.begin(), m_advancelimits.end(), limits);
611 } else {
612 std::fill(limits, limits + m_nv, -1.0);
613 }
614 return has_limit;
615}
616
617void Reactor::setAdvanceLimit(const string& nm, const double limit)
618{
619 size_t k = componentIndex(nm);
620 if (k == npos) {
621 throw CanteraError("Reactor::setAdvanceLimit", "No component named '{}'", nm);
622 }
623
624 if (m_thermo == 0) {
625 throw CanteraError("Reactor::setAdvanceLimit",
626 "Error: reactor is empty.");
627 }
628 if (m_nv == 0) {
629 if (m_net == 0) {
630 throw CanteraError("Reactor::setAdvanceLimit",
631 "Cannot set limit on a reactor that is not "
632 "assigned to a ReactorNet object.");
633 } else {
634 m_net->initialize();
635 }
636 } else if (k > m_nv) {
637 throw CanteraError("Reactor::setAdvanceLimit",
638 "Index out of bounds.");
639 }
640 m_advancelimits.resize(m_nv, -1.0);
641 m_advancelimits[k] = limit;
642
643 // resize to zero length if no limits are set
644 if (std::none_of(m_advancelimits.begin(), m_advancelimits.end(),
645 [](double val){return val>0;})) {
646 m_advancelimits.resize(0);
647 }
648}
649
650}
Base class for kinetics managers and also contains the kineticsmgr module documentation (see Kinetics...
Header file for class ReactorSurface.
Header for a simple thermodynamics model of a surface phase derived from ThermoPhase,...
Header file for base class WallBase.
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:431
Base class for exceptions thrown by Cantera classes.
Public interface for kinetics managers.
Definition Kinetics.h:126
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition Kinetics.h:274
virtual void getNetProductionRates(double *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
Definition Kinetics.cpp:411
An error indicating that an unimplemented function has been called.
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:232
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:142
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition Phase.cpp:129
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
Definition SurfPhase.h:98
double size(size_t k) const
Returns the number of sites occupied by one molecule of species k.
Definition SurfPhase.h:221
double siteDensity() const
Returns the site density.
Definition SurfPhase.h:216
Base class for a phase with thermodynamic properties.
Base class for 'walls' (walls, pistons, etc.) connecting reactors.
Definition Wall.h:23
virtual void initialize()
Called just before the start of integration.
Definition Wall.h:68
double dot(InputIter x_begin, InputIter x_end, InputIter2 y_begin)
Function that calculates a templated inner product.
Definition utilities.h:82
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition ct_defs.h:120
Namespace for the Cantera kernel.
Definition AnyMap.cpp:595
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
const double Tiny
Small number to compare differences of mole fractions against.
Definition ct_defs.h:173
offset
Offsets of solution components in the 1D solution array.
Definition Flow1D.h:24
void warn_deprecated(const string &source, const AnyBase &node, const string &message)
A deprecation warning for syntax in an input file.
Definition AnyMap.cpp:1997
const double BigNumber
largest number to compare to inf.
Definition ct_defs.h:160
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...