17#include <boost/math/tools/roots.hpp>
20namespace bmt = boost::math::tools;
25Reactor::Reactor(shared_ptr<Solution> sol,
const string& name)
26 : ReactorBase(sol, name)
28 m_kin = sol->kinetics().get();
29 if (m_kin->nReactions() == 0) {
37void Reactor::setDerivativeSettings(
AnyMap& settings)
39 m_kin->setDerivativeSettings(settings);
41 for (
auto S : m_surfaces) {
42 S->kinetics()->setDerivativeSettings(settings);
49 "After Cantera 3.2, a change of reactor contents after instantiation "
52 if (m_kin->nReactions() == 0) {
59void Reactor::getState(
double* y)
63 "Error: reactor is empty.");
65 m_thermo->restoreState(m_state);
68 m_mass = m_thermo->density() * m_vol;
75 y[2] = m_thermo->intEnergy_mass() * m_mass;
78 m_thermo->getMassFractions(y+3);
82 getSurfaceInitialConditions(y + m_nsp + 3);
85void Reactor::getSurfaceInitialConditions(
double* y)
88 for (
auto& S : m_surfaces) {
89 S->getCoverages(y + loc);
90 loc += S->thermo()->nSpecies();
94void Reactor::initialize(
double t0)
96 if (!m_thermo || (m_chem && !m_kin)) {
97 throw CanteraError(
"Reactor::initialize",
"Reactor contents not set"
98 " for reactor '" + m_name +
"'.");
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);
105 for (
size_t n = 0; n < m_wall.size(); n++) {
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);
119 m_work.resize(maxnt);
122size_t Reactor::nSensParams()
const
124 size_t ns = m_sensParams.size();
125 for (
auto& S : m_surfaces) {
126 ns += S->nSensParams();
131void Reactor::updateState(
double* y)
138 m_thermo->setMassFractions_NoNorm(y+3);
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;
148 double T = m_thermo->temperature();
149 boost::uintmax_t maxiter = 100;
150 pair<double, double> TT;
152 TT = bmt::bracket_and_solve_root(
153 u_err, T, 1.2,
true, bmt::eps_tolerance<double>(48), maxiter);
154 }
catch (std::exception&) {
158 TT = bmt::bisect(u_err, m_thermo->minTemp(), m_thermo->maxTemp(),
159 bmt::eps_tolerance<double>(48), maxiter);
160 }
catch (std::exception& err2) {
162 m_thermo->setState_TD(T, m_mass / m_vol);
164 "{}\nat U = {}, rho = {}", err2.what(), U, m_mass / m_vol);
167 if (fabs(TT.first - TT.second) > 1e-7*TT.first) {
168 throw CanteraError(
"Reactor::updateState",
"root finding failed");
170 m_thermo->setState_TD(TT.second, m_mass / m_vol);
172 m_thermo->setDensity(m_mass/m_vol);
175 updateConnected(
true);
176 updateSurfaceState(y + m_nsp + 3);
179void Reactor::updateSurfaceState(
double* y)
182 for (
auto& S : m_surfaces) {
183 S->setCoverages(y+loc);
184 loc += S->thermo()->nSpecies();
188void Reactor::updateConnected(
bool updatePressure) {
190 m_enthalpy = m_thermo->enthalpy_mass();
191 if (updatePressure) {
192 m_pressure = m_thermo->pressure();
195 m_intEnergy = m_thermo->intEnergy_mass();
199 m_thermo->saveState(m_state);
203 if (m_net !=
nullptr) {
204 time = (timeIsIndependent()) ? m_net->time() : m_net->distance();
206 for (
size_t i = 0; i < m_outlet.size(); i++) {
207 m_outlet[i]->setSimTime(time);
208 m_outlet[i]->updateMassFlowRate(time);
210 for (
size_t i = 0; i < m_inlet.size(); i++) {
211 m_inlet[i]->setSimTime(time);
212 m_inlet[i]->updateMassFlowRate(time);
214 for (
size_t i = 0; i < m_wall.size(); i++) {
215 m_wall[i]->setSimTime(time);
219void Reactor::eval(
double time,
double* LHS,
double* RHS)
221 double& dmdt = RHS[0];
222 double* mdYdt = RHS + 3;
225 m_thermo->restoreState(m_state);
226 const vector<double>& mw = m_thermo->molecularWeights();
227 const double* Y = m_thermo->massFractions();
229 evalSurfaces(LHS + m_nsp + 3, RHS + m_nsp + 3, m_sdot.data());
231 double mdot_surf =
dot(m_sdot.begin(), m_sdot.end(), mw.begin());
238 m_kin->getNetProductionRates(&m_wdot[0]);
241 for (
size_t k = 0; k < m_nsp; k++) {
243 mdYdt[k] = (m_wdot[k] * m_vol + m_sdot[k]) * mw[k];
245 mdYdt[k] -= Y[k] * mdot_surf;
254 RHS[2] = - m_thermo->pressure() * m_vdot + m_Qdot;
260 for (
auto outlet : m_outlet) {
261 double mdot = outlet->massFlowRate();
264 RHS[2] -= mdot * m_enthalpy;
269 for (
auto inlet : m_inlet) {
270 double mdot = inlet->massFlowRate();
272 for (
size_t n = 0; n < m_nsp; n++) {
273 double mdot_spec = inlet->outletSpeciesMassFlowRate(n);
275 mdYdt[n] += (mdot_spec - mdot * Y[n]);
278 RHS[2] += mdot * inlet->enthalpy_mass();
283void Reactor::evalWalls(
double t)
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();
295void Reactor::evalSurfaces(
double* LHS,
double* RHS,
double* sdot)
297 fill(sdot, sdot + m_nsp, 0.0);
300 for (
auto S : m_surfaces) {
309 for (
size_t k = 1; k < nk; k++) {
310 RHS[loc + k] = m_work[k] * rs0 * surf->
size(k);
317 double wallarea = S->area();
318 for (
size_t k = 0; k < m_nsp; k++) {
319 sdot[k] += m_work[bulkloc + k] * wallarea;
324vector<size_t> Reactor::steadyConstraints()
const
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());
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.");
340Eigen::SparseMatrix<double> Reactor::finiteDifferenceJacobian()
344 "Reactor must be initialized first.");
349 Eigen::ArrayXd yCurrent(m_nv);
350 getState(yCurrent.data());
351 double time = (m_net !=
nullptr) ? m_net->time() : 0.0;
353 Eigen::ArrayXd yPerturbed = yCurrent;
354 Eigen::ArrayXd lhsPerturbed(m_nv), lhsCurrent(m_nv);
355 Eigen::ArrayXd rhsPerturbed(m_nv), rhsCurrent(m_nv);
358 updateState(yCurrent.data());
359 eval(time, lhsCurrent.data(), rhsCurrent.data());
361 double rel_perturb = std::sqrt(std::numeric_limits<double>::epsilon());
362 double atol = (m_net !=
nullptr) ? m_net->atol() : 1e-15;
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;
369 updateState(yPerturbed.data());
372 eval(time, lhsPerturbed.data(), rhsPerturbed.data());
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);
385 updateState(yCurrent.data());
387 Eigen::SparseMatrix<double> jac(m_nv, m_nv);
388 jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
393void Reactor::evalSurfaces(
double* RHS,
double* sdot)
395 fill(sdot, sdot + m_nsp, 0.0);
398 for (
auto S : m_surfaces) {
407 for (
size_t k = 1; k < nk; k++) {
408 RHS[loc + k] = m_work[k] * rs0 * surf->
size(k);
415 double wallarea = S->area();
416 for (
size_t k = 0; k < m_nsp; k++) {
417 sdot[k] += m_work[bulkloc + k] * wallarea;
422void Reactor::addSensitivityReaction(
size_t rxn)
424 if (!m_chem || rxn >= m_kin->nReactions()) {
426 "Reaction number out of range ({})", rxn);
429 size_t p = network().registerSensitivityParameter(
430 name()+
": "+m_kin->reaction(rxn)->equation(), 1.0, 1.0);
431 m_sensParams.emplace_back(
435void Reactor::addSensitivitySpeciesEnthalpy(
size_t k)
437 if (k >= m_thermo->nSpecies()) {
438 throw CanteraError(
"Reactor::addSensitivitySpeciesEnthalpy",
439 "Species index out of range ({})", k);
442 size_t p = network().registerSensitivityParameter(
443 name() +
": " + m_thermo->speciesName(k) +
" enthalpy",
445 m_sensParams.emplace_back(
447 SensParameterType::enthalpy});
450size_t Reactor::speciesIndex(
const string& nm)
const
453 size_t k = m_thermo->speciesIndex(nm);
460 for (
auto& S : m_surfaces) {
472size_t Reactor::componentIndex(
const string& nm)
const
474 size_t k = speciesIndex(nm);
477 }
else if (nm ==
"mass") {
479 }
else if (nm ==
"volume") {
481 }
else if (nm ==
"int_energy") {
488string Reactor::componentName(
size_t k) {
495 }
else if (k >= 3 && k < neq()) {
497 if (k < m_thermo->nSpecies()) {
498 return m_thermo->speciesName(k);
500 k -= m_thermo->nSpecies();
502 for (
auto& S : m_surfaces) {
504 if (k < th->nSpecies()) {
511 throw CanteraError(
"Reactor::componentName",
"Index is out of bounds.");
514double Reactor::upperBound(
size_t k)
const {
521 }
else if (k >= 3 && k < m_nv) {
524 throw CanteraError(
"Reactor::upperBound",
"Index {} is out of bounds.", k);
528double Reactor::lowerBound(
size_t k)
const {
535 }
else if (k >= 3 && k < m_nv) {
538 throw CanteraError(
"Reactor::lowerBound",
"Index {} is out of bounds.", k);
542void Reactor::resetBadValues(
double* y) {
543 for (
size_t k = 3; k < m_nv; k++) {
544 y[k] = std::max(y[k], 0.0);
548void Reactor::applySensitivity(
double* params)
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]);
561 for (
auto& S : m_surfaces) {
562 S->setSensitivityParameters(params);
564 m_thermo->invalidateCache();
566 m_kin->invalidateCache();
570void Reactor::resetSensitivity(
double* params)
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);
582 for (
auto& S : m_surfaces) {
583 S->resetSensitivityParameters();
585 m_thermo->invalidateCache();
587 m_kin->invalidateCache();
591void Reactor::setAdvanceLimits(
const double *limits)
595 "Error: reactor is empty.");
597 m_advancelimits.assign(limits, limits + m_nv);
600 if (std::none_of(m_advancelimits.begin(), m_advancelimits.end(),
601 [](
double val){return val>0;})) {
602 m_advancelimits.resize(0);
606bool Reactor::getAdvanceLimits(
double *limits)
const
608 bool has_limit = hasAdvanceLimits();
610 std::copy(m_advancelimits.begin(), m_advancelimits.end(), limits);
612 std::fill(limits, limits + m_nv, -1.0);
617void Reactor::setAdvanceLimit(
const string& nm,
const double limit)
619 size_t k = componentIndex(nm);
621 throw CanteraError(
"Reactor::setAdvanceLimit",
"No component named '{}'", nm);
626 "Error: reactor is empty.");
631 "Cannot set limit on a reactor that is not "
632 "assigned to a ReactorNet object.");
636 }
else if (k > m_nv) {
638 "Index out of bounds.");
640 m_advancelimits.resize(m_nv, -1.0);
641 m_advancelimits[k] = limit;
644 if (std::none_of(m_advancelimits.begin(), m_advancelimits.end(),
645 [](
double val){return val>0;})) {
646 m_advancelimits.resize(0);
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.
Base class for exceptions thrown by Cantera classes.
Public interface for kinetics managers.
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
virtual void getNetProductionRates(double *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
An error indicating that an unimplemented function has been called.
size_t nSpecies() const
Returns the number of species in the phase.
string speciesName(size_t k) const
Name of the species with index k.
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
double size(size_t k) const
Returns the number of sites occupied by one molecule of species k.
double siteDensity() const
Returns the site density.
Base class for a phase with thermodynamic properties.
Base class for 'walls' (walls, pistons, etc.) connecting reactors.
virtual void initialize()
Called just before the start of integration.
double dot(InputIter x_begin, InputIter x_end, InputIter2 y_begin)
Function that calculates a templated inner product.
const double GasConstant
Universal Gas Constant [J/kmol/K].
Namespace for the Cantera kernel.
const size_t npos
index returned by functions to indicate "no position"
const double Tiny
Small number to compare differences of mole fractions against.
offset
Offsets of solution components in the 1D solution array.
void warn_deprecated(const string &source, const AnyBase &node, const string &message)
A deprecation warning for syntax in an input file.
const double BigNumber
largest number to compare to inf.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...