34 void Reactor::getInitialConditions(
double t0,
size_t leny,
double* y)
38 cout <<
"Error: reactor is empty." << endl;
45 doublereal mass = m_thermo->
density() * m_vol;
61 size_t loc =
m_nsp + 2;
63 for (
size_t m = 0; m < m_nwalls; m++) {
64 surf = m_wall[m]->surface(m_lr[m]);
66 m_wall[m]->getCoverages(m_lr[m], y + loc);
68 loc += surf->nSpecies();
79 m_sdot.resize(
m_nsp, 0.0);
81 for (
size_t w = 0; w < m_nwalls; w++)
82 if (m_wall[w]->surface(m_lr[w])) {
83 m_nv += m_wall[w]->surface(m_lr[w])->nSpecies();
90 size_t nt = 0, maxnt = 0;
91 for (
size_t m = 0; m < m_nwalls; m++) {
92 if (m_wall[m]->kinetics(m_lr[m])) {
93 nt = m_wall[m]->kinetics(m_lr[m])->nTotalSpecies();
97 if (m_wall[m]->kinetics(m_lr[m])) {
99 &m_wall[m]->kinetics(m_lr[m])->thermo(0)) {
101 "First phase of all kinetics managers must be"
107 m_work.resize(maxnt);
111 size_t Reactor::nSensParams()
113 if (m_nsens ==
npos) {
116 m_nsens = m_pnum.size();
117 for (m = 0; m < m_nwalls; m++) {
118 ns = m_wall[m]->nSensParams(m_lr[m]);
119 m_nsens_wall.push_back(ns);
138 doublereal* mss = y + 2;
139 doublereal mass = accumulate(y+2, y+2+
m_nsp, 0.0);
144 doublereal temp = temperature();
157 size_t loc =
m_nsp + 2;
159 for (
size_t m = 0; m < m_nwalls; m++) {
160 surf = m_wall[m]->surface(m_lr[m]);
164 m_wall[m]->setCoverages(m_lr[m], y+loc);
180 void Reactor::evalEqs(doublereal time, doublereal* y,
181 doublereal* ydot, doublereal* params)
188 size_t npar = m_pnum.size();
189 for (
size_t n = 0; n < npar; n++) {
194 for (
size_t m = 0; m < m_nwalls; m++) {
195 if (m_nsens_wall[m] > 0) {
196 m_wall[m]->setSensitivityParameters(m_lr[m], params + ploc);
197 ploc += m_nsens_wall[m];
206 size_t loc =
m_nsp+2;
207 fill(m_sdot.begin(), m_sdot.end(), 0.0);
208 for (
size_t i = 0; i < m_nwalls; i++) {
209 int lr = 1 - 2*m_lr[i];
210 double vdot = lr*m_wall[i]->vdot(time);
212 m_Q += lr*m_wall[i]->Q(time);
213 Kinetics* kin = m_wall[i]->kinetics(m_lr[i]);
214 SurfPhase* surf = m_wall[i]->surface(m_lr[i]);
216 double rs0 = 1.0/surf->siteDensity();
217 size_t nk = surf->nSpecies();
219 surf->setTemperature(m_state[0]);
220 m_wall[i]->syncCoverages(m_lr[i]);
221 kin->getNetProductionRates(
DATA_PTR(m_work));
222 size_t ns = kin->surfacePhaseIndex();
223 size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
224 for (
size_t k = 1; k < nk; k++) {
225 ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
226 sum -= ydot[loc + k];
231 double wallarea = m_wall[i]->area();
232 for (
size_t k = 0; k <
m_nsp; k++) {
233 m_sdot[k] += m_work[k]*wallarea;
250 fill(ydot + 2, ydot + 2 +
m_nsp, 0.0);
252 for (
size_t n = 0; n <
m_nsp; n++) {
254 ydot[n+2] += m_sdot[n];
266 ydot[0] = - m_thermo->
pressure() * m_vdot - m_Q;
277 for (
size_t i = 0; i < m_nOutlets; i++) {
278 double mdot_out = m_outlet[i]->massFlowRate(time);
279 for (
size_t n = 0; n <
m_nsp; n++) {
280 ydot[2+n] -= mdot_out * mf[n];
288 for (
size_t i = 0; i < m_nInlets; i++) {
289 double mdot_in = m_inlet[i]->massFlowRate(time);
290 for (
size_t n = 0; n <
m_nsp; n++) {
291 ydot[2+n] += m_inlet[i]->outletSpeciesMassFlowRate(n);
294 ydot[0] += mdot_in * m_inlet[i]->enthalpy_mass();
301 size_t npar = m_pnum.size();
302 for (
size_t n = 0; n < npar; n++) {
307 for (
size_t m = 0; m < m_nwalls; m++) {
308 if (m_nsens_wall[m] > 0) {
309 m_wall[m]->resetSensitivityParameters(m_lr[m]);
310 ploc += m_nsens_wall[m];
316 void Reactor::addSensitivityReaction(
size_t rxn)
319 throw CanteraError(
"Reactor::addSensitivityReaction",
320 "Reaction number out of range ("+
int2str(rxn)+
")");
321 m_pnum.push_back(rxn);
323 m_mult_save.push_back(1.0);
327 size_t Reactor::componentIndex(
string nm)
const
342 size_t walloffset = 0, kp = 0;
344 for (
size_t m = 0; m < m_nwalls; m++) {
345 if (m_wall[m]->kinetics(m_lr[m])) {
346 kp = m_wall[m]->kinetics(m_lr[m])->reactionPhaseIndex();
347 th = &m_wall[m]->kinetics(m_lr[m])->thermo(kp);
348 k = th->speciesIndex(nm);
350 return k + 2 + m_nsp + walloffset;
352 walloffset += th->nSpecies();