23 throw CanteraError(
"IdealGasConstPressureMoleReactor::getState",
24 "Error: reactor is empty.");
39 if (m_thermo->
type() !=
"ideal-gas") {
40 throw CanteraError(
"IdealGasConstPressureMoleReactor::initialize",
41 "Incompatible phase type '{}' provided", m_thermo->
type());
62 double& mcpdTdt = RHS[0];
63 double* dndt = RHS + m_sidx;
82 for (
size_t n = 0; n <
m_nsp; n++) {
91 for (
auto outlet : m_outlet) {
92 for (
size_t n = 0; n <
m_nsp; n++) {
99 for (
auto inlet : m_inlet) {
102 for (
size_t n = 0; n <
m_nsp; n++) {
106 mcpdTdt -=
m_hk[n] * imw[n] * mdot_spec;
120 throw CanteraError(
"IdealGasConstPressureMoleReactor::jacobian",
121 "Reactor must be initialized first.");
129 size_t ssize = m_nv - m_sidx;
132 if (!m_surfaces.empty()) {
133 vector<Eigen::Triplet<double>> species_trips(dnk_dnj.nonZeros());
134 for (
int k = 0; k < dnk_dnj.outerSize(); k++) {
135 for (Eigen::SparseMatrix<double>::InnerIterator it(dnk_dnj, k); it; ++it) {
136 species_trips.emplace_back(
static_cast<int>(it.row()),
137 static_cast<int>(it.col()), it.value());
141 dnk_dnj.resize(ssize, ssize);
142 dnk_dnj.setFromTriplets(species_trips.begin(), species_trips.end());
145 Eigen::VectorXd netProductionRates = Eigen::VectorXd::Zero(ssize);
149 for (
auto &S: m_surfaces) {
150 auto curr_kin = S->kinetics();
151 vector<double> prod_rates(curr_kin->nTotalSpecies());
152 curr_kin->getNetProductionRates(prod_rates.data());
153 for (
size_t i = 0; i < curr_kin->nTotalSpecies(); i++) {
154 size_t row =
speciesIndex(curr_kin->kineticsSpeciesName(i));
156 netProductionRates[row] += prod_rates[i];
164 for (
int k = 0; k < dnk_dnj.outerSize(); k++) {
165 for (Eigen::SparseMatrix<double>::InnerIterator it(dnk_dnj, k); it; ++it) {
167 if (
static_cast<size_t>(it.row()) <
m_nsp) {
168 it.valueRef() = it.value() + netProductionRates[it.row()] * molarVol;
170 m_jac_trips.emplace_back(
static_cast<int>(it.row() + m_sidx),
171 static_cast<int>(it.col() + m_sidx), it.value());
178 * std::sqrt(std::numeric_limits<double>::epsilon());
180 vector<double> yCurrent(m_nv);
183 vector<double> lhsPerturbed(m_nv, 1.0), lhsCurrent(m_nv, 1.0);
184 vector<double> rhsPerturbed(m_nv, 0.0), rhsCurrent(m_nv, 0.0);
185 vector<double> yPerturbed = yCurrent;
187 yPerturbed[0] += deltaTemp;
191 eval(time, lhsPerturbed.data(), rhsPerturbed.data());
194 eval(time, lhsCurrent.data(), rhsCurrent.data());
196 for (
size_t j = 0; j < m_nv; j++) {
197 double ydotPerturbed = rhsPerturbed[j] / lhsPerturbed[j];
198 double ydotCurrent = rhsCurrent[j] / lhsCurrent[j];
200 (ydotPerturbed - ydotCurrent) / deltaTemp);
204 Eigen::VectorXd enthalpy = Eigen::VectorXd::Zero(ssize);
205 Eigen::VectorXd specificHeat = Eigen::VectorXd::Zero(ssize);
210 for (
size_t i = 0; i <
m_nsp; i++) {
211 netProductionRates[i] *=
m_vol;
213 double qdot = enthalpy.dot(netProductionRates);
216 double* moles = yCurrent.data() + m_sidx;
217 for (
size_t i = 0; i < ssize; i++) {
218 NCp += moles[i] * specificHeat[i];
220 double denom = 1 / (NCp * NCp);
221 Eigen::VectorXd hk_dnkdnj_sums = dnk_dnj.transpose() * enthalpy;
223 for (
size_t j = 0; j < ssize; j++) {
224 m_jac_trips.emplace_back(0,
static_cast<int>(j + m_sidx),
225 (specificHeat[j] * qdot - NCp * hk_dnkdnj_sums[j]) * denom);
229 Eigen::SparseMatrix<double> jac(m_nv, m_nv);
239 }
else if (nm ==
"temperature") {
248 return "temperature";
249 }
else if (k >= m_sidx && k <
neq()) {
251 if (k < m_thermo->nSpecies()) {
256 for (
auto& S : m_surfaces) {
258 if (k < th->nSpecies()) {
265 throw CanteraError(
"IdealGasConstPressureMoleReactor::componentName",
266 "Index is out of bounds.");
272 return 1.5 * m_thermo->
maxTemp();
281 return 0.5 * m_thermo->
minTemp();
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.
double lowerBound(size_t k) const override
Get the lower bound on the k-th component of the local state vector.
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).
double upperBound(size_t k) const override
Get the upper bound on the k-th component of the local state vector.
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.
double lowerBound(size_t k) const override
Get the lower bound on the k-th component of the local state vector.
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].
double m_mass
Current mass of the reactor [kg].
size_t m_nsp
Number of homogeneous species in the mixture.
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 double minTemp(size_t k=npos) const
Minimum temperature for which the thermodynamic data for the species or phase are valid.
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.
virtual double maxTemp(size_t k=npos) const
Maximum temperature for which the thermodynamic data for the species are valid.
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"
const double BigNumber
largest number to compare to inf.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...