22 if (thermo.
type() !=
"ideal-gas") {
23 throw CanteraError(
"IdealGasConstPressureMoleReactor::setThermoMgr",
24 "Incompatible phase type provided");
32 throw CanteraError(
"IdealGasConstPressureMoleReactor::getState",
33 "Error: reactor is empty.");
67 double& mcpdTdt = RHS[0];
68 double* dndt = RHS + m_sidx;
87 for (
size_t n = 0; n <
m_nsp; n++) {
96 for (
auto outlet : m_outlet) {
97 for (
size_t n = 0; n <
m_nsp; n++) {
104 for (
auto inlet : m_inlet) {
107 for (
size_t n = 0; n <
m_nsp; n++) {
111 mcpdTdt -=
m_hk[n] * imw[n] * mdot_spec;
125 throw CanteraError(
"IdealGasConstPressureMoleReactor::jacobian",
126 "Reactor must be initialized first.");
134 size_t ssize = m_nv - m_sidx;
137 if (!m_surfaces.empty()) {
138 vector<Eigen::Triplet<double>> species_trips(dnk_dnj.nonZeros());
139 for (
int k = 0; k < dnk_dnj.outerSize(); k++) {
140 for (Eigen::SparseMatrix<double>::InnerIterator it(dnk_dnj, k); it; ++it) {
141 species_trips.emplace_back(
static_cast<int>(it.row()),
142 static_cast<int>(it.col()), it.value());
146 dnk_dnj.resize(ssize, ssize);
147 dnk_dnj.setFromTriplets(species_trips.begin(), species_trips.end());
150 Eigen::VectorXd netProductionRates = Eigen::VectorXd::Zero(ssize);
154 for (
auto &S: m_surfaces) {
155 auto curr_kin = S->kinetics();
156 vector<double> prod_rates(curr_kin->nTotalSpecies());
157 curr_kin->getNetProductionRates(prod_rates.data());
158 for (
size_t i = 0; i < curr_kin->nTotalSpecies(); i++) {
159 size_t row =
speciesIndex(curr_kin->kineticsSpeciesName(i));
161 netProductionRates[row] += prod_rates[i];
169 for (
int k = 0; k < dnk_dnj.outerSize(); k++) {
170 for (Eigen::SparseMatrix<double>::InnerIterator it(dnk_dnj, k); it; ++it) {
172 if (
static_cast<size_t>(it.row()) <
m_nsp) {
173 it.valueRef() = it.value() + netProductionRates[it.row()] * molarVol;
175 m_jac_trips.emplace_back(
static_cast<int>(it.row() + m_sidx),
176 static_cast<int>(it.col() + m_sidx), it.value());
183 * std::sqrt(std::numeric_limits<double>::epsilon());
185 vector<double> yCurrent(m_nv);
188 vector<double> lhsPerturbed(m_nv, 1.0), lhsCurrent(m_nv, 1.0);
189 vector<double> rhsPerturbed(m_nv, 0.0), rhsCurrent(m_nv, 0.0);
190 vector<double> yPerturbed = yCurrent;
192 yPerturbed[0] += deltaTemp;
196 eval(time, lhsPerturbed.data(), rhsPerturbed.data());
199 eval(time, lhsCurrent.data(), rhsCurrent.data());
201 for (
size_t j = 0; j < m_nv; j++) {
202 double ydotPerturbed = rhsPerturbed[j] / lhsPerturbed[j];
203 double ydotCurrent = rhsCurrent[j] / lhsCurrent[j];
205 (ydotPerturbed - ydotCurrent) / deltaTemp);
209 Eigen::VectorXd enthalpy = Eigen::VectorXd::Zero(ssize);
210 Eigen::VectorXd specificHeat = Eigen::VectorXd::Zero(ssize);
215 for (
size_t i = 0; i <
m_nsp; i++) {
216 netProductionRates[i] *=
m_vol;
218 double qdot = enthalpy.dot(netProductionRates);
221 double* moles = yCurrent.data() + m_sidx;
222 for (
size_t i = 0; i < ssize; i++) {
223 NCp += moles[i] * specificHeat[i];
225 double denom = 1 / (NCp * NCp);
226 Eigen::VectorXd hk_dnkdnj_sums = dnk_dnj.transpose() * enthalpy;
228 for (
size_t j = 0; j < ssize; j++) {
229 m_jac_trips.emplace_back(0,
static_cast<int>(j + m_sidx),
230 (specificHeat[j] * qdot - NCp * hk_dnkdnj_sums[j]) * denom);
234 Eigen::SparseMatrix<double> jac(m_nv, m_nv);
244 }
else if (nm ==
"temperature") {
253 return "temperature";
254 }
else if (k >= m_sidx && k <
neq()) {
256 if (k < m_thermo->nSpecies()) {
261 for (
auto& S : m_surfaces) {
263 if (k < th->nSpecies()) {
270 throw CanteraError(
"IdealGasConstPressureMoleReactor::componentName",
271 "Index is out of bounds.");
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 class ThermoPhase, the base class for phases with thermodynamic properties,...
Base class for exceptions thrown by Cantera classes.
void initialize(double t0=0.0) override
Initialize the reactor.
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).
void setThermoMgr(ThermoPhase &thermo) override
Specify the mixture contained in the reactor.
void eval(double t, double *LHS, double *RHS) override
Evaluate the reactor governing equations.
size_t componentIndex(const string &nm) const override
Return the index in the solution vector for this reactor of the component named nm.
Eigen::SparseMatrix< double > jacobian() override
Calculate an approximate Jacobian to accelerate preconditioned solvers.
void getState(double *y) override
Get the the current state of the reactor.
string componentName(size_t k) override
Return the name of the solution component with index i.
void updateState(double *y) override
Set the state of the reactor to correspond to the state vector y.
void initialize(double t0=0.0) override
Initialize the reactor.
vector< double > m_hk
Species molar enthalpies.
virtual void getNetProductionRates(double *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
void evalSurfaces(double *LHS, double *RHS, double *sdot) override
Evaluate terms related to surface reactions.
void getSurfaceInitialConditions(double *y) override
Get initial conditions for SurfPhase objects attached to this reactor.
void getMoles(double *y)
Get moles of the system from mass fractions stored by thermo object.
void setMassFromMoles(double *y)
Set internal mass variable based on moles given.
virtual void addSurfaceJacobian(vector< Eigen::Triplet< double > > &triplets)
For each surface in the reactor, update vector of triplets with all relevant surface jacobian derivat...
void updateSurfaceState(double *y) override
Update the state of SurfPhase objects attached to this reactor.
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.
double temperature() const
Temperature (K).
string speciesName(size_t k) const
Name of the species with index k.
const vector< double > & inverseMolecularWeights() const
Return a const reference to the internal vector of molecular weights.
virtual double density() const
Density (kg/m^3).
virtual void setMolesNoTruncate(const double *const N)
Set the state of the object with moles in [kmol].
virtual double molarVolume() const
Molar volume (m^3/kmol).
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.
double m_vol
Current volume of the reactor [m^3].
size_t m_nsp
Number of homogeneous species in the mixture.
virtual void setThermoMgr(ThermoPhase &thermo)
Specify the mixture contained in the reactor.
double time()
Current value of the simulation time [s], for reactor networks that are solved in the time domain.
Kinetics * m_kin
Pointer to the homogeneous Kinetics object that handles the reactions.
vector< double > m_wdot
Species net molar production rates.
virtual void evalWalls(double t)
Evaluate terms related to Walls.
size_t neq()
Number of equations (state variables) for this reactor.
double m_Qdot
net heat transfer into the reactor, through walls [W]
vector< Eigen::Triplet< double > > m_jac_trips
Vector of triplets representing the jacobian.
vector< double > m_sdot
Production rates of gas phase species on surfaces [kmol/s].
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 void updateConnected(bool updatePressure)
Update the state information needed by connected reactors, flow devices, and reactor walls.
Base class for a phase with thermodynamic properties.
virtual void getPartialMolarEnthalpies(double *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
virtual void setState_TP(double t, double p)
Set the temperature (K) and pressure (Pa)
virtual void getPartialMolarCp(double *cpbar) const
Return an array of partial molar heat capacities for the species in the mixture.
string type() const override
String indicating the thermodynamic model implemented.
double cp_mass() const
Specific heat at constant pressure. Units: J/kg/K.
Eigen::SparseMatrix< double > netProductionRates_ddCi()
Calculate derivatives for species net production rates with respect to species concentration at const...
Namespace for the Cantera kernel.
const size_t npos
index returned by functions to indicate "no position"
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...