17#include <boost/math/tools/roots.hpp>
20namespace bmt = boost::math::tools;
34 for (
auto S : m_surfaces) {
35 S->kinetics()->setDerivativeSettings(settings);
53 "Error: reactor is empty.");
78 for (
auto& S : m_surfaces) {
79 S->getCoverages(y + loc);
80 loc += S->thermo()->nSpecies();
86 if (!m_thermo || (m_chem && !
m_kin)) {
87 throw CanteraError(
"Reactor::initialize",
"Reactor contents not set"
88 " for reactor '" + m_name +
"'.");
95 for (
size_t n = 0; n < m_wall.size(); n++) {
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);
109 m_work.resize(maxnt);
114 size_t ns = m_sensParams.size();
115 for (
auto& S : m_surfaces) {
116 ns += S->nSensParams();
139 auto u_err = [
this, U](
double T) {
145 boost::uintmax_t maxiter = 100;
146 pair<double, double> TT;
148 TT = bmt::bracket_and_solve_root(
149 u_err, T, 1.2,
true, bmt::eps_tolerance<double>(48), maxiter);
150 }
catch (std::exception&) {
155 bmt::eps_tolerance<double>(48), maxiter);
156 }
catch (std::exception& err2) {
160 "{}\nat U = {}, rho = {}", err2.what(), U,
m_mass /
m_vol);
163 if (fabs(TT.first - TT.second) > 1e-7*TT.first) {
164 throw CanteraError(
"Reactor::updateState",
"root finding failed");
178 for (
auto& S : m_surfaces) {
179 S->setCoverages(y+loc);
180 loc += S->thermo()->nSpecies();
187 if (updatePressure) {
195 if (
m_net !=
nullptr) {
198 for (
size_t i = 0; i < m_outlet.size(); i++) {
199 m_outlet[i]->setSimTime(time);
200 m_outlet[i]->updateMassFlowRate(time);
202 for (
size_t i = 0; i < m_inlet.size(); i++) {
203 m_inlet[i]->setSimTime(time);
204 m_inlet[i]->updateMassFlowRate(time);
206 for (
size_t i = 0; i < m_wall.size(); i++) {
207 m_wall[i]->setSimTime(time);
213 double& dmdt = RHS[0];
214 double* mdYdt = RHS + 3;
233 for (
size_t k = 0; k <
m_nsp; k++) {
237 mdYdt[k] -= Y[k] * mdot_surf;
252 for (
auto outlet : m_outlet) {
261 for (
auto inlet : m_inlet) {
264 for (
size_t n = 0; n <
m_nsp; n++) {
267 mdYdt[n] += (mdot_spec - mdot * Y[n]);
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();
289 fill(sdot, sdot +
m_nsp, 0.0);
292 for (
auto S : m_surfaces) {
303 for (
size_t k = 1; k < nk; k++) {
304 RHS[loc + k] = m_work[surfloc + k] * rs0 * surf->
size(k);
311 double wallarea = S->area();
312 for (
size_t k = 0; k <
m_nsp; k++) {
313 sdot[k] += m_work[bulkloc + k] * wallarea;
322 "Reactor must be initialized first.");
327 Eigen::ArrayXd yCurrent(m_nv);
331 Eigen::ArrayXd yPerturbed = yCurrent;
332 Eigen::ArrayXd lhsPerturbed(m_nv), lhsCurrent(m_nv);
333 Eigen::ArrayXd rhsPerturbed(m_nv), rhsCurrent(m_nv);
337 eval(time, lhsCurrent.data(), rhsCurrent.data());
339 double rel_perturb = std::sqrt(std::numeric_limits<double>::epsilon());
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;
350 eval(time, lhsPerturbed.data(), rhsPerturbed.data());
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) {
358 static_cast<int>(i),
static_cast<int>(j),
359 (ydotPerturbed - ydotCurrent) / delta_y);
365 Eigen::SparseMatrix<double> jac(m_nv, m_nv);
373 fill(sdot, sdot +
m_nsp, 0.0);
376 for (
auto S : m_surfaces) {
387 for (
size_t k = 1; k < nk; k++) {
388 RHS[loc + k] = m_work[surfloc + k] * rs0 * surf->
size(k);
395 double wallarea = S->area();
396 for (
size_t k = 0; k <
m_nsp; k++) {
397 sdot[k] += m_work[bulkloc + k] * wallarea;
406 "Reaction number out of range ({})", rxn);
411 m_sensParams.emplace_back(
412 SensitivityParameter{rxn, p, 1.0, SensParameterType::reaction});
418 throw CanteraError(
"Reactor::addSensitivitySpeciesEnthalpy",
419 "Species index out of range ({})", k);
425 m_sensParams.emplace_back(
426 SensitivityParameter{k, p, m_thermo->
Hf298SS(k),
427 SensParameterType::enthalpy});
440 for (
auto& S : m_surfaces) {
457 }
else if (nm ==
"mass") {
459 }
else if (nm ==
"volume") {
461 }
else if (nm ==
"int_energy") {
475 }
else if (k >= 3 && k <
neq()) {
477 if (k < m_thermo->nSpecies()) {
482 for (
auto& S : m_surfaces) {
484 if (k < th->nSpecies()) {
491 throw CanteraError(
"Reactor::componentName",
"Index is out of bounds.");
499 for (
auto& p : m_sensParams) {
500 if (p.type == SensParameterType::reaction) {
503 }
else if (p.type == SensParameterType::enthalpy) {
507 for (
auto& S : m_surfaces) {
508 S->setSensitivityParameters(params);
512 m_kin->invalidateCache();
521 for (
auto& p : m_sensParams) {
522 if (p.type == SensParameterType::reaction) {
524 }
else if (p.type == SensParameterType::enthalpy) {
528 for (
auto& S : m_surfaces) {
529 S->resetSensitivityParameters();
533 m_kin->invalidateCache();
541 "Error: reactor is empty.");
547 [](
double val){return val>0;})) {
558 std::fill(limits, limits + m_nv, -1.0);
567 throw CanteraError(
"Reactor::setAdvanceLimit",
"No component named '{}'", nm);
572 "Error: reactor is empty.");
577 "Cannot set limit on a reactor that is not "
578 "assigned to a ReactorNet object.");
582 }
else if (k > m_nv) {
584 "Index out of bounds.");
591 [](
double val){return val>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.
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).
Public interface for kinetics managers.
size_t reactionPhaseIndex() const
Phase where the reactions occur.
double multiplier(size_t i) const
The current value of the multiplier for reaction i.
shared_ptr< Reaction > reaction(size_t i)
Return the Reaction object for reaction i.
size_t nReactions() const
Number of reactions in the reaction mechanism.
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
virtual void setMultiplier(size_t i, double f)
Set the multiplier for reaction i to f.
virtual void getNetProductionRates(double *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
void restoreState(const vector< double > &state)
Restore a state saved on a previous call to saveState.
size_t nSpecies() const
Returns the number of species in the phase.
virtual void setMassFractions_NoNorm(const double *const y)
Set the mass fractions to the specified values without normalizing.
void setState_TD(double t, double rho)
Set the internally stored temperature (K) and density (kg/m^3)
double temperature() const
Temperature (K).
void saveState(vector< double > &state) const
Save the current internal state of the phase.
string speciesName(size_t k) const
Name of the species with index k.
virtual void setDensity(const double density_)
Set the internally stored density (kg/m^3) of the phase.
const double * massFractions() const
Return a const pointer to the mass fraction array.
const vector< double > & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
virtual double density() const
Density (kg/m^3).
void getMassFractions(double *const y) const
Get the species mass fractions.
virtual double pressure() const
Return the thermodynamic pressure (Pa).
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.
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.
virtual string componentName(size_t k)
Return the name of the solution component with index i.
void setKineticsMgr(Kinetics &kin) override
Specify chemical kinetics governing the reactor.
virtual void evalSurfaces(double *LHS, double *RHS, double *sdot)
Evaluate terms related to surface reactions.
virtual size_t componentIndex(const string &nm) const
Return the index in the solution vector for this reactor of the component named nm.
virtual void updateSurfaceState(double *y)
Update the state of SurfPhase objects attached to this reactor.
void insert(G &contents)
Insert something into the reactor.
virtual void applySensitivity(double *params)
Set reaction rate multipliers based on the sensitivity variables in params.
Kinetics * m_kin
Pointer to the homogeneous Kinetics object that handles the reactions.
vector< double > m_wdot
Species net molar production rates.
virtual void eval(double t, double *LHS, double *RHS)
Evaluate the reactor governing equations.
virtual void evalWalls(double t)
Evaluate terms related to Walls.
Eigen::SparseMatrix< double > finiteDifferenceJacobian()
Calculate the reactor-specific Jacobian using a finite difference method.
size_t neq()
Number of equations (state variables) for this reactor.
virtual void addSensitivitySpeciesEnthalpy(size_t k)
Add a sensitivity parameter associated with the enthalpy formation of species k (in the homogeneous p...
double m_Qdot
net heat transfer into the reactor, through walls [W]
vector< double > m_advancelimits
!< Number of variables associated with reactor surfaces
virtual void addSensitivityReaction(size_t rxn)
Add a sensitivity parameter associated with the reaction number rxn (in the homogeneous phase).
virtual size_t nSensParams() const
Number of sensitivity parameters associated with this reactor (including walls)
virtual void setDerivativeSettings(AnyMap &settings)
Use this to set the kinetics objects derivative settings.
virtual void updateState(double *y)
Set the state of the reactor to correspond to the state vector y.
void setChemistry(bool cflag=true) override
Enable or disable changes in reactor composition due to chemical reactions.
void setAdvanceLimits(const double *limits)
Set absolute step size limits during advance.
vector< Eigen::Triplet< double > > m_jac_trips
Vector of triplets representing the jacobian.
void setAdvanceLimit(const string &nm, const double limit)
Set individual step size limit for component name nm
virtual void resetSensitivity(double *params)
Reset the reaction rate multipliers.
vector< double > m_sdot
Production rates of gas phase species on surfaces [kmol/s].
virtual void getState(double *y)
Get the the current state of the reactor.
bool hasAdvanceLimits() const
Check whether Reactor object uses advance limits.
double m_vdot
net rate of volume change from moving walls [m^3/s]
void syncState() override
Set the state of the reactor to correspond to the state of the associated ThermoPhase object.
virtual void getSurfaceInitialConditions(double *y)
Get initial conditions for SurfPhase objects attached to this reactor.
void initialize(double t0=0.0) override
Initialize the reactor.
bool getAdvanceLimits(double *limits) const
Retrieve absolute step size limits during advance.
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...
virtual bool timeIsIndependent() const
Indicates whether the governing equations for this reactor are functions of time or a spatial variabl...
virtual void updateConnected(bool updatePressure)
Update the state information needed by connected reactors, flow devices, and reactor walls.
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.
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.
virtual void initialize()
Called just before the start of integration.
virtual void setDerivativeSettings(const AnyMap &settings)
Set/modify derivative settings.
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"
offset
Offsets of solution components in the 1D solution array.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...