Cantera  3.0.0
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
25void Reactor::insert(shared_ptr<Solution> sol) {
26 setThermoMgr(*sol->thermo());
27 setKineticsMgr(*sol->kinetics());
28}
29
31{
33 // translate settings to surfaces
34 for (auto S : m_surfaces) {
35 S->kinetics()->setDerivativeSettings(settings);
36 }
37}
38
40{
41 m_kin = &kin;
42 if (m_kin->nReactions() == 0) {
43 setChemistry(false);
44 } else {
45 setChemistry(true);
46 }
47}
48
49void Reactor::getState(double* y)
50{
51 if (m_thermo == 0) {
52 throw CanteraError("Reactor::getState",
53 "Error: reactor is empty.");
54 }
55 m_thermo->restoreState(m_state);
56
57 // set the first component to the total mass
58 m_mass = m_thermo->density() * m_vol;
59 y[0] = m_mass;
60
61 // set the second component to the total volume
62 y[1] = m_vol;
63
64 // set the third component to the total internal energy
65 y[2] = m_thermo->intEnergy_mass() * m_mass;
66
67 // set components y+3 ... y+K+2 to the mass fractions of each species
68 m_thermo->getMassFractions(y+3);
69
70 // set the remaining components to the surface species
71 // coverages on the walls
73}
74
76{
77 size_t loc = 0;
78 for (auto& S : m_surfaces) {
79 S->getCoverages(y + loc);
80 loc += S->thermo()->nSpecies();
81 }
82}
83
84void Reactor::initialize(double t0)
85{
86 if (!m_thermo || (m_chem && !m_kin)) {
87 throw CanteraError("Reactor::initialize", "Reactor contents not set"
88 " for reactor '" + m_name + "'.");
89 }
90 m_thermo->restoreState(m_state);
91 m_sdot.resize(m_nsp, 0.0);
92 m_wdot.resize(m_nsp, 0.0);
93 updateConnected(true);
94
95 for (size_t n = 0; n < m_wall.size(); n++) {
96 WallBase* W = m_wall[n];
97 W->initialize();
98 }
99
100 m_nv = m_nsp + 3;
101 m_nv_surf = 0;
102 size_t maxnt = 0;
103 for (auto& S : m_surfaces) {
104 m_nv_surf += S->thermo()->nSpecies();
105 size_t nt = S->kinetics()->nTotalSpecies();
106 maxnt = std::max(maxnt, nt);
107 }
108 m_nv += m_nv_surf;
109 m_work.resize(maxnt);
110}
111
113{
114 size_t ns = m_sensParams.size();
115 for (auto& S : m_surfaces) {
116 ns += S->nSensParams();
117 }
118 return ns;
119}
120
122{
124 m_mass = m_thermo->density() * m_vol;
125}
126
127void Reactor::updateState(double* y)
128{
129 // The components of y are [0] the total mass, [1] the total volume,
130 // [2] the total internal energy, [3...K+3] are the mass fractions of each
131 // species, and [K+3...] are the coverages of surface species on each wall.
132 m_mass = y[0];
133 m_vol = y[1];
134 m_thermo->setMassFractions_NoNorm(y+3);
135
136 if (m_energy) {
137 double U = y[2];
138 // Residual function: error in internal energy as a function of T
139 auto u_err = [this, U](double T) {
140 m_thermo->setState_TD(T, m_mass / m_vol);
141 return m_thermo->intEnergy_mass() * m_mass - U;
142 };
143
144 double T = m_thermo->temperature();
145 boost::uintmax_t maxiter = 100;
146 pair<double, double> TT;
147 try {
148 TT = bmt::bracket_and_solve_root(
149 u_err, T, 1.2, true, bmt::eps_tolerance<double>(48), maxiter);
150 } catch (std::exception&) {
151 // Try full-range bisection if bracketing fails (for example, near
152 // temperature limits for the phase's equation of state)
153 try {
154 TT = bmt::bisect(u_err, m_thermo->minTemp(), m_thermo->maxTemp(),
155 bmt::eps_tolerance<double>(48), maxiter);
156 } catch (std::exception& err2) {
157 // Set m_thermo back to a reasonable state if root finding fails
158 m_thermo->setState_TD(T, m_mass / m_vol);
159 throw CanteraError("Reactor::updateState",
160 "{}\nat U = {}, rho = {}", err2.what(), U, m_mass / m_vol);
161 }
162 }
163 if (fabs(TT.first - TT.second) > 1e-7*TT.first) {
164 throw CanteraError("Reactor::updateState", "root finding failed");
165 }
166 m_thermo->setState_TD(TT.second, m_mass / m_vol);
167 } else {
168 m_thermo->setDensity(m_mass/m_vol);
169 }
170
171 updateConnected(true);
172 updateSurfaceState(y + m_nsp + 3);
173}
174
176{
177 size_t loc = 0;
178 for (auto& S : m_surfaces) {
179 S->setCoverages(y+loc);
180 loc += S->thermo()->nSpecies();
181 }
182}
183
184void Reactor::updateConnected(bool updatePressure) {
185 // save parameters needed by other connected reactors
186 m_enthalpy = m_thermo->enthalpy_mass();
187 if (updatePressure) {
188 m_pressure = m_thermo->pressure();
189 }
190 m_intEnergy = m_thermo->intEnergy_mass();
191 m_thermo->saveState(m_state);
192
193 // Update the mass flow rate of connected flow devices
194 double time = 0.0;
195 if (m_net != nullptr) {
196 time = (timeIsIndependent()) ? m_net->time() : m_net->distance();
197 }
198 for (size_t i = 0; i < m_outlet.size(); i++) {
199 m_outlet[i]->setSimTime(time);
200 m_outlet[i]->updateMassFlowRate(time);
201 }
202 for (size_t i = 0; i < m_inlet.size(); i++) {
203 m_inlet[i]->setSimTime(time);
204 m_inlet[i]->updateMassFlowRate(time);
205 }
206 for (size_t i = 0; i < m_wall.size(); i++) {
207 m_wall[i]->setSimTime(time);
208 }
209}
210
211void Reactor::eval(double time, double* LHS, double* RHS)
212{
213 double& dmdt = RHS[0];
214 double* mdYdt = RHS + 3; // mass * dY/dt
215
216 evalWalls(time);
217 m_thermo->restoreState(m_state);
218 const vector<double>& mw = m_thermo->molecularWeights();
219 const double* Y = m_thermo->massFractions();
220
221 evalSurfaces(LHS + m_nsp + 3, RHS + m_nsp + 3, m_sdot.data());
222 // mass added to gas phase from surface reactions
223 double mdot_surf = dot(m_sdot.begin(), m_sdot.end(), mw.begin());
224 dmdt = mdot_surf;
225
226 // volume equation
227 RHS[1] = m_vdot;
228
229 if (m_chem) {
230 m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
231 }
232
233 for (size_t k = 0; k < m_nsp; k++) {
234 // production in gas phase and from surfaces
235 mdYdt[k] = (m_wdot[k] * m_vol + m_sdot[k]) * mw[k];
236 // dilution by net surface mass flux
237 mdYdt[k] -= Y[k] * mdot_surf;
238 LHS[k+3] = m_mass;
239 }
240
241 // Energy equation.
242 // @f[
243 // \dot U = -P\dot V + A \dot q + \dot m_{in} h_{in} - \dot m_{out} h.
244 // @f]
245 if (m_energy) {
246 RHS[2] = - m_thermo->pressure() * m_vdot + m_Qdot;
247 } else {
248 RHS[2] = 0.0;
249 }
250
251 // add terms for outlets
252 for (auto outlet : m_outlet) {
253 double mdot = outlet->massFlowRate();
254 dmdt -= mdot; // mass flow out of system
255 if (m_energy) {
256 RHS[2] -= mdot * m_enthalpy;
257 }
258 }
259
260 // add terms for inlets
261 for (auto inlet : m_inlet) {
262 double mdot = inlet->massFlowRate();
263 dmdt += mdot; // mass flow into system
264 for (size_t n = 0; n < m_nsp; n++) {
265 double mdot_spec = inlet->outletSpeciesMassFlowRate(n);
266 // flow of species into system and dilution by other species
267 mdYdt[n] += (mdot_spec - mdot * Y[n]);
268 }
269 if (m_energy) {
270 RHS[2] += mdot * inlet->enthalpy_mass();
271 }
272 }
273}
274
275void Reactor::evalWalls(double t)
276{
277 // time is currently unused
278 m_vdot = 0.0;
279 m_Qdot = 0.0;
280 for (size_t i = 0; i < m_wall.size(); i++) {
281 int f = 2 * m_lr[i] - 1;
282 m_vdot -= f * m_wall[i]->expansionRate();
283 m_Qdot += f * m_wall[i]->heatRate();
284 }
285}
286
287void Reactor::evalSurfaces(double* LHS, double* RHS, double* sdot)
288{
289 fill(sdot, sdot + m_nsp, 0.0);
290 size_t loc = 0; // offset into ydot
291
292 for (auto S : m_surfaces) {
293 Kinetics* kin = S->kinetics();
294 SurfPhase* surf = S->thermo();
295
296 double rs0 = 1.0/surf->siteDensity();
297 size_t nk = surf->nSpecies();
298 double sum = 0.0;
299 S->syncState();
300 kin->getNetProductionRates(&m_work[0]);
301 size_t ns = kin->reactionPhaseIndex();
302 size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
303 for (size_t k = 1; k < nk; k++) {
304 RHS[loc + k] = m_work[surfloc + k] * rs0 * surf->size(k);
305 sum -= RHS[loc + k];
306 }
307 RHS[loc] = sum;
308 loc += nk;
309
310 size_t bulkloc = kin->kineticsSpeciesIndex(m_thermo->speciesName(0));
311 double wallarea = S->area();
312 for (size_t k = 0; k < m_nsp; k++) {
313 sdot[k] += m_work[bulkloc + k] * wallarea;
314 }
315 }
316}
317
318Eigen::SparseMatrix<double> Reactor::finiteDifferenceJacobian()
319{
320 if (m_nv == 0) {
321 throw CanteraError("Reactor::finiteDifferenceJacobian",
322 "Reactor must be initialized first.");
323 }
324 // clear former jacobian elements
325 m_jac_trips.clear();
326
327 Eigen::ArrayXd yCurrent(m_nv);
328 getState(yCurrent.data());
329 double time = (m_net != nullptr) ? m_net->time() : 0.0;
330
331 Eigen::ArrayXd yPerturbed = yCurrent;
332 Eigen::ArrayXd lhsPerturbed(m_nv), lhsCurrent(m_nv);
333 Eigen::ArrayXd rhsPerturbed(m_nv), rhsCurrent(m_nv);
334 lhsCurrent = 1.0;
335 rhsCurrent = 0.0;
336 updateState(yCurrent.data());
337 eval(time, lhsCurrent.data(), rhsCurrent.data());
338
339 double rel_perturb = std::sqrt(std::numeric_limits<double>::epsilon());
340 double atol = (m_net != nullptr) ? m_net->atol() : 1e-15;
341
342 for (size_t j = 0; j < m_nv; j++) {
343 yPerturbed = yCurrent;
344 double delta_y = std::max(std::abs(yCurrent[j]), 1000 * atol) * rel_perturb;
345 yPerturbed[j] += delta_y;
346
347 updateState(yPerturbed.data());
348 lhsPerturbed = 1.0;
349 rhsPerturbed = 0.0;
350 eval(time, lhsPerturbed.data(), rhsPerturbed.data());
351
352 // d ydot_i/dy_j
353 for (size_t i = 0; i < m_nv; i++) {
354 double ydotPerturbed = rhsPerturbed[i] / lhsPerturbed[i];
355 double ydotCurrent = rhsCurrent[i] / lhsCurrent[i];
356 if (ydotCurrent != ydotPerturbed) {
357 m_jac_trips.emplace_back(
358 static_cast<int>(i), static_cast<int>(j),
359 (ydotPerturbed - ydotCurrent) / delta_y);
360 }
361 }
362 }
363 updateState(yCurrent.data());
364
365 Eigen::SparseMatrix<double> jac(m_nv, m_nv);
366 jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
367 return jac;
368}
369
370
371void Reactor::evalSurfaces(double* RHS, double* sdot)
372{
373 fill(sdot, sdot + m_nsp, 0.0);
374 size_t loc = 0; // offset into ydot
375
376 for (auto S : m_surfaces) {
377 Kinetics* kin = S->kinetics();
378 SurfPhase* surf = S->thermo();
379
380 double rs0 = 1.0/surf->siteDensity();
381 size_t nk = surf->nSpecies();
382 double sum = 0.0;
383 S->syncState();
384 kin->getNetProductionRates(&m_work[0]);
385 size_t ns = kin->reactionPhaseIndex();
386 size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
387 for (size_t k = 1; k < nk; k++) {
388 RHS[loc + k] = m_work[surfloc + k] * rs0 * surf->size(k);
389 sum -= RHS[loc + k];
390 }
391 RHS[loc] = sum;
392 loc += nk;
393
394 size_t bulkloc = kin->kineticsSpeciesIndex(m_thermo->speciesName(0));
395 double wallarea = S->area();
396 for (size_t k = 0; k < m_nsp; k++) {
397 sdot[k] += m_work[bulkloc + k] * wallarea;
398 }
399 }
400}
401
403{
404 if (!m_chem || rxn >= m_kin->nReactions()) {
405 throw CanteraError("Reactor::addSensitivityReaction",
406 "Reaction number out of range ({})", rxn);
407 }
408
410 name()+": "+m_kin->reaction(rxn)->equation(), 1.0, 1.0);
411 m_sensParams.emplace_back(
412 SensitivityParameter{rxn, p, 1.0, SensParameterType::reaction});
413}
414
416{
417 if (k >= m_thermo->nSpecies()) {
418 throw CanteraError("Reactor::addSensitivitySpeciesEnthalpy",
419 "Species index out of range ({})", k);
420 }
421
423 name() + ": " + m_thermo->speciesName(k) + " enthalpy",
424 0.0, GasConstant * 298.15);
425 m_sensParams.emplace_back(
426 SensitivityParameter{k, p, m_thermo->Hf298SS(k),
427 SensParameterType::enthalpy});
428}
429
430size_t Reactor::speciesIndex(const string& nm) const
431{
432 // check for a gas species name
433 size_t k = m_thermo->speciesIndex(nm);
434 if (k != npos) {
435 return k;
436 }
437
438 // check for a wall species
439 size_t offset = m_nsp;
440 for (auto& S : m_surfaces) {
441 ThermoPhase* th = S->thermo();
442 k = th->speciesIndex(nm);
443 if (k != npos) {
444 return k + offset;
445 } else {
446 offset += th->nSpecies();
447 }
448 }
449 return npos;
450}
451
452size_t Reactor::componentIndex(const string& nm) const
453{
454 size_t k = speciesIndex(nm);
455 if (k != npos) {
456 return k + 3;
457 } else if (nm == "mass") {
458 return 0;
459 } else if (nm == "volume") {
460 return 1;
461 } else if (nm == "int_energy") {
462 return 2;
463 } else {
464 return npos;
465 }
466}
467
468string Reactor::componentName(size_t k) {
469 if (k == 0) {
470 return "mass";
471 } else if (k == 1) {
472 return "volume";
473 } else if (k == 2) {
474 return "int_energy";
475 } else if (k >= 3 && k < neq()) {
476 k -= 3;
477 if (k < m_thermo->nSpecies()) {
478 return m_thermo->speciesName(k);
479 } else {
480 k -= m_thermo->nSpecies();
481 }
482 for (auto& S : m_surfaces) {
483 ThermoPhase* th = S->thermo();
484 if (k < th->nSpecies()) {
485 return th->speciesName(k);
486 } else {
487 k -= th->nSpecies();
488 }
489 }
490 }
491 throw CanteraError("Reactor::componentName", "Index is out of bounds.");
492}
493
494void Reactor::applySensitivity(double* params)
495{
496 if (!params) {
497 return;
498 }
499 for (auto& p : m_sensParams) {
500 if (p.type == SensParameterType::reaction) {
501 p.value = m_kin->multiplier(p.local);
502 m_kin->setMultiplier(p.local, p.value*params[p.global]);
503 } else if (p.type == SensParameterType::enthalpy) {
504 m_thermo->modifyOneHf298SS(p.local, p.value + params[p.global]);
505 }
506 }
507 for (auto& S : m_surfaces) {
508 S->setSensitivityParameters(params);
509 }
510 m_thermo->invalidateCache();
511 if (m_kin) {
512 m_kin->invalidateCache();
513 }
514}
515
516void Reactor::resetSensitivity(double* params)
517{
518 if (!params) {
519 return;
520 }
521 for (auto& p : m_sensParams) {
522 if (p.type == SensParameterType::reaction) {
523 m_kin->setMultiplier(p.local, p.value);
524 } else if (p.type == SensParameterType::enthalpy) {
525 m_thermo->resetHf298(p.local);
526 }
527 }
528 for (auto& S : m_surfaces) {
529 S->resetSensitivityParameters();
530 }
531 m_thermo->invalidateCache();
532 if (m_kin) {
533 m_kin->invalidateCache();
534 }
535}
536
537void Reactor::setAdvanceLimits(const double *limits)
538{
539 if (m_thermo == 0) {
540 throw CanteraError("Reactor::setAdvanceLimits",
541 "Error: reactor is empty.");
542 }
543 m_advancelimits.assign(limits, limits + m_nv);
544
545 // resize to zero length if no limits are set
546 if (std::none_of(m_advancelimits.begin(), m_advancelimits.end(),
547 [](double val){return val>0;})) {
548 m_advancelimits.resize(0);
549 }
550}
551
552bool Reactor::getAdvanceLimits(double *limits) const
553{
554 bool has_limit = hasAdvanceLimits();
555 if (has_limit) {
556 std::copy(m_advancelimits.begin(), m_advancelimits.end(), limits);
557 } else {
558 std::fill(limits, limits + m_nv, -1.0);
559 }
560 return has_limit;
561}
562
563void Reactor::setAdvanceLimit(const string& nm, const double limit)
564{
565 size_t k = componentIndex(nm);
566 if (k == npos) {
567 throw CanteraError("Reactor::setAdvanceLimit", "No component named '{}'", nm);
568 }
569
570 if (m_thermo == 0) {
571 throw CanteraError("Reactor::setAdvanceLimit",
572 "Error: reactor is empty.");
573 }
574 if (m_nv == 0) {
575 if (m_net == 0) {
576 throw CanteraError("Reactor::setAdvanceLimit",
577 "Cannot set limit on a reactor that is not "
578 "assigned to a ReactorNet object.");
579 } else {
580 m_net->initialize();
581 }
582 } else if (k > m_nv) {
583 throw CanteraError("Reactor::setAdvanceLimit",
584 "Index out of bounds.");
585 }
586 m_advancelimits.resize(m_nv, -1.0);
587 m_advancelimits[k] = limit;
588
589 // resize to zero length if no limits are set
590 if (std::none_of(m_advancelimits.begin(), m_advancelimits.end(),
591 [](double val){return val>0;})) {
592 m_advancelimits.resize(0);
593 }
594}
595
596}
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:427
Base class for exceptions thrown by Cantera classes.
double outletSpeciesMassFlowRate(size_t k)
Mass flow rate (kg/s) of outlet species k.
double enthalpy_mass()
specific enthalpy
double massFlowRate()
Mass flow rate (kg/s).
Definition FlowDevice.h:39
Public interface for kinetics managers.
Definition Kinetics.h:126
size_t reactionPhaseIndex() const
Phase where the reactions occur.
Definition Kinetics.h:235
double multiplier(size_t i) const
The current value of the multiplier for reaction i.
Definition Kinetics.h:1434
shared_ptr< Reaction > reaction(size_t i)
Return the Reaction object for reaction i.
Definition Kinetics.cpp:784
size_t nReactions() const
Number of reactions in the reaction mechanism.
Definition Kinetics.h:153
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition Kinetics.h:288
virtual void setMultiplier(size_t i, double f)
Set the multiplier for reaction i to f.
Definition Kinetics.h:1443
virtual void getNetProductionRates(double *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
Definition Kinetics.cpp:435
void restoreState(const vector< double > &state)
Restore a state saved on a previous call to saveState.
Definition Phase.cpp:275
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:245
virtual void setMassFractions_NoNorm(const double *const y)
Set the mass fractions to the specified values without normalizing.
Definition Phase.cpp:370
void setState_TD(double t, double rho)
Set the internally stored temperature (K) and density (kg/m^3)
Definition Phase.cpp:444
double temperature() const
Temperature (K).
Definition Phase.h:662
void saveState(vector< double > &state) const
Save the current internal state of the phase.
Definition Phase.cpp:251
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:151
virtual void setDensity(const double density_)
Set the internally stored density (kg/m^3) of the phase.
Definition Phase.cpp:707
const double * massFractions() const
Return a const pointer to the mass fraction array.
Definition Phase.h:535
const vector< double > & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition Phase.cpp:501
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition Phase.cpp:138
virtual double density() const
Density (kg/m^3).
Definition Phase.h:687
void getMassFractions(double *const y) const
Get the species mass fractions.
Definition Phase.cpp:584
virtual double pressure() const
Return the thermodynamic pressure (Pa).
Definition Phase.h:680
FlowDevice & outlet(size_t n=0)
Return a reference to the n-th outlet FlowDevice connected to this reactor.
double m_pressure
Current pressure in the reactor [Pa].
ReactorNet * m_net
The ReactorNet that this reactor is part of.
FlowDevice & inlet(size_t n=0)
Return a reference to the n-th inlet FlowDevice connected to this reactor.
vector< int > m_lr
Vector of length nWalls(), indicating whether this reactor is on the left (0) or right (1) of each wa...
double m_vol
Current volume of the reactor [m^3].
virtual void syncState()
Set the state of the reactor to correspond to the state of the associated ThermoPhase object.
double m_intEnergy
Current internal energy of the reactor [J/kg].
size_t m_nsp
Number of homogeneous species in the mixture.
virtual void setThermoMgr(ThermoPhase &thermo)
Specify the mixture contained in the reactor.
ReactorNet & network()
The ReactorNet that this reactor belongs to.
double m_enthalpy
Current specific enthalpy of the reactor [J/kg].
string name() const
Return the name of this reactor.
Definition ReactorBase.h:64
void initialize()
Initialize the reactor network.
double time()
Current value of the simulation time [s], for reactor networks that are solved in the time domain.
size_t registerSensitivityParameter(const string &name, double value, double scale)
Used by Reactor and Wall objects to register the addition of sensitivity parameters so that the React...
double distance()
Current position [m] along the length of the reactor network, for reactors that are solved as a funct...
double atol()
Absolute integration tolerance.
Definition ReactorNet.h:95
virtual string componentName(size_t k)
Return the name of the solution component with index i.
Definition Reactor.cpp:468
void setKineticsMgr(Kinetics &kin) override
Specify chemical kinetics governing the reactor.
Definition Reactor.cpp:39
virtual void evalSurfaces(double *LHS, double *RHS, double *sdot)
Evaluate terms related to surface reactions.
Definition Reactor.cpp:287
virtual size_t componentIndex(const string &nm) const
Return the index in the solution vector for this reactor of the component named nm.
Definition Reactor.cpp:452
virtual void updateSurfaceState(double *y)
Update the state of SurfPhase objects attached to this reactor.
Definition Reactor.cpp:175
void insert(G &contents)
Insert something into the reactor.
Definition Reactor.h:70
virtual void applySensitivity(double *params)
Set reaction rate multipliers based on the sensitivity variables in params.
Definition Reactor.cpp:494
Kinetics * m_kin
Pointer to the homogeneous Kinetics object that handles the reactions.
Definition Reactor.h:277
vector< double > m_wdot
Species net molar production rates.
Definition Reactor.h:289
virtual void eval(double t, double *LHS, double *RHS)
Evaluate the reactor governing equations.
Definition Reactor.cpp:211
virtual void evalWalls(double t)
Evaluate terms related to Walls.
Definition Reactor.cpp:275
Eigen::SparseMatrix< double > finiteDifferenceJacobian()
Calculate the reactor-specific Jacobian using a finite difference method.
Definition Reactor.cpp:318
size_t neq()
Number of equations (state variables) for this reactor.
Definition Reactor.h:103
virtual void addSensitivitySpeciesEnthalpy(size_t k)
Add a sensitivity parameter associated with the enthalpy formation of species k (in the homogeneous p...
Definition Reactor.cpp:415
double m_Qdot
net heat transfer into the reactor, through walls [W]
Definition Reactor.h:281
vector< double > m_advancelimits
!< Number of variables associated with reactor surfaces
Definition Reactor.h:296
virtual void addSensitivityReaction(size_t rxn)
Add a sensitivity parameter associated with the reaction number rxn (in the homogeneous phase).
Definition Reactor.cpp:402
virtual size_t nSensParams() const
Number of sensitivity parameters associated with this reactor (including walls)
Definition Reactor.cpp:112
virtual void setDerivativeSettings(AnyMap &settings)
Use this to set the kinetics objects derivative settings.
Definition Reactor.cpp:30
virtual void updateState(double *y)
Set the state of the reactor to correspond to the state vector y.
Definition Reactor.cpp:127
void setChemistry(bool cflag=true) override
Enable or disable changes in reactor composition due to chemical reactions.
Definition Reactor.h:79
void setAdvanceLimits(const double *limits)
Set absolute step size limits during advance.
Definition Reactor.cpp:537
double m_mass
total mass
Definition Reactor.h:283
vector< Eigen::Triplet< double > > m_jac_trips
Vector of triplets representing the jacobian.
Definition Reactor.h:302
void setAdvanceLimit(const string &nm, const double limit)
Set individual step size limit for component name nm
Definition Reactor.cpp:563
virtual void resetSensitivity(double *params)
Reset the reaction rate multipliers.
Definition Reactor.cpp:516
vector< double > m_sdot
Production rates of gas phase species on surfaces [kmol/s].
Definition Reactor.h:287
virtual void getState(double *y)
Get the the current state of the reactor.
Definition Reactor.cpp:49
bool hasAdvanceLimits() const
Check whether Reactor object uses advance limits.
Definition Reactor.h:186
double m_vdot
net rate of volume change from moving walls [m^3/s]
Definition Reactor.h:279
void syncState() override
Set the state of the reactor to correspond to the state of the associated ThermoPhase object.
Definition Reactor.cpp:121
virtual void getSurfaceInitialConditions(double *y)
Get initial conditions for SurfPhase objects attached to this reactor.
Definition Reactor.cpp:75
void initialize(double t0=0.0) override
Initialize the reactor.
Definition Reactor.cpp:84
bool getAdvanceLimits(double *limits) const
Retrieve absolute step size limits during advance.
Definition Reactor.cpp:552
virtual size_t speciesIndex(const string &nm) const
Return the index in the solution vector for this reactor of the species named nm, in either the homog...
Definition Reactor.cpp:430
virtual bool timeIsIndependent() const
Indicates whether the governing equations for this reactor are functions of time or a spatial variabl...
Definition Reactor.h:61
virtual void updateConnected(bool updatePressure)
Update the state information needed by connected reactors, flow devices, and reactor walls.
Definition Reactor.cpp:184
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:226
double siteDensity() const
Returns the site density.
Definition SurfPhase.h:221
Base class for a phase with thermodynamic properties.
virtual double minTemp(size_t k=npos) const
Minimum temperature for which the thermodynamic data for the species or phase are valid.
virtual double maxTemp(size_t k=npos) const
Maximum temperature for which the thermodynamic data for the species are valid.
virtual void modifyOneHf298SS(const size_t k, const double Hf298New)
Modify the value of the 298 K Heat of Formation of one species in the phase (J kmol-1)
void invalidateCache() override
Invalidate any cached values which are normally updated only when a change in state is detected.
virtual void resetHf298(const size_t k=npos)
Restore the original heat of formation of one or more species.
double intEnergy_mass() const
Specific internal energy. Units: J/kg.
double Hf298SS(const size_t k) const
Report the 298 K Heat of Formation of the standard state of one species (J kmol-1)
double enthalpy_mass() const
Specific enthalpy. Units: J/kg.
Base class for 'walls' (walls, pistons, etc.) connecting reactors.
Definition Wall.h:22
virtual void initialize()
Called just before the start of integration.
Definition Wall.h:93
virtual void setDerivativeSettings(const AnyMap &settings)
Set/modify derivative settings.
Definition Kinetics.h:732
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:564
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:195
offset
Offsets of solution components in the 1D solution array.
Definition StFlow.h:24
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...