From 6314e7a563c4bafd8a74fcfdfe8e2c5d25e89290 Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Sat, 6 May 2006 14:56:57 +0000 Subject: [PATCH] added constant pressure reactor --- Cantera/src/zeroD/ConstPressureReactor.cpp | 312 +++++++++++++++++++++ Cantera/src/zeroD/ConstPressureReactor.h | 73 +++++ Cantera/src/zeroD/FlowReactor.cpp | 25 +- Cantera/src/zeroD/Reactor.cpp | 56 +--- Cantera/src/zeroD/Reactor.h | 99 +------ Cantera/src/zeroD/ReactorBase.cpp | 9 +- Cantera/src/zeroD/ReactorBase.h | 8 +- 7 files changed, 426 insertions(+), 156 deletions(-) create mode 100644 Cantera/src/zeroD/ConstPressureReactor.cpp create mode 100644 Cantera/src/zeroD/ConstPressureReactor.h diff --git a/Cantera/src/zeroD/ConstPressureReactor.cpp b/Cantera/src/zeroD/ConstPressureReactor.cpp new file mode 100644 index 000000000..077bed44e --- /dev/null +++ b/Cantera/src/zeroD/ConstPressureReactor.cpp @@ -0,0 +1,312 @@ +/** + * @file Reactor.cpp + * + * A zero-dimensional reactor + */ + +// Copyright 2001 California Institute of Technology + + +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "ConstPressureReactor.h" +#include "FlowDevice.h" +#include "Wall.h" +#include "../InterfaceKinetics.h" +#include "../SurfPhase.h" + +using namespace Cantera; + +namespace CanteraZeroD { + + ConstPressureReactor::ConstPressureReactor() : Reactor() {} + + void ConstPressureReactor:: + getInitialConditions(double t0, size_t leny, double* y) + { + m_init = true; + if (m_thermo == 0) { + throw CanteraError("getInitialConditions", + "Error: reactor is empty."); + } + m_time = t0; + m_thermo->restoreState(m_state); + + // total mass + doublereal mass = m_thermo->density() * m_vol; + + // set components y + 2 ... y + K + 1 to the + // mass M_k of each species + m_thermo->getMassFractions(y+2); + scale(y + 2, y + m_nsp + 2, y + 2, mass); + + // set the first component to the total enthalpy + y[0] = m_thermo->enthalpy_mass() * mass; + + // set the second component to the total volume + y[1] = m_vol; + + // set the remaining components to the surface species + // coverages on the walls + int loc = m_nsp + 2; + SurfPhase* surf; + for (int m = 0; m < m_nwalls; m++) { + surf = m_wall[m]->surface(m_lr[m]); + if (surf) { + m_wall[m]->getCoverages(m_lr[m], y + loc); + loc += surf->nSpecies(); + } + } + } + + void ConstPressureReactor::initialize(doublereal t0) { + m_thermo->restoreState(m_state); + m_sdot.resize(m_nsp, 0.0); + m_nv = m_nsp + 2; + for (int w = 0; w < m_nwalls; w++) + if (m_wall[w]->surface(m_lr[w])) + m_nv += m_wall[w]->surface(m_lr[w])->nSpecies(); + m_enthalpy = m_thermo->enthalpy_mass(); + m_pressure = m_thermo->pressure(); + m_intEnergy = m_thermo->intEnergy_mass(); + + int m, nt = 0, maxnt = 0; + for (m = 0; m < m_nwalls; m++) { + if (m_wall[m]->kinetics(m_lr[m])) { + nt = m_wall[m]->kinetics(m_lr[m])->nTotalSpecies(); + if (nt > maxnt) maxnt = nt; + if (m_wall[m]->kinetics(m_lr[m])) { + if (&m_kin->thermo(0) != + &m_wall[m]->kinetics(m_lr[m])->thermo(0)) { + throw CanteraError("ConstPressureReactor::initialize", + "First phase of all kinetics managers must be" + " the gas."); + } + } + } + } + m_work.resize(maxnt); + m_init = true; + } + + void ConstPressureReactor::updateState(doublereal* y) { + + // The components of y are the total enthalpy, + // the total volume, and the mass of each species. + + doublereal h = y[0]; + doublereal* mss = y + 2; + doublereal mass = accumulate(y+2, y+2+m_nsp, 0.0); + m_thermo->setMassFractions(mss); + + if (m_energy) { + m_thermo->setState_HP(h/mass, m_pressure, 1.0e-4); + } + else { + m_thermo->setPressure(m_pressure); + } + m_vol = mass / m_thermo->density(); + + int loc = m_nsp + 2; + SurfPhase* surf; + for (int m = 0; m < m_nwalls; m++) { + surf = m_wall[m]->surface(m_lr[m]); + if (surf) { + m_wall[m]->setCoverages(m_lr[m], y+loc); + loc += surf->nSpecies(); + } + } + + // save parameters needed by other connected reactors + m_enthalpy = m_thermo->enthalpy_mass(); + m_intEnergy = m_thermo->intEnergy_mass(); + m_thermo->saveState(m_state); + } + + + /* + * Called by the integrator to evaluate ydot given y at time 'time'. + */ + void ConstPressureReactor::evalEqs(doublereal time, doublereal* y, + doublereal* ydot, doublereal* params) + { + int i, k, nk; + m_time = time; + m_thermo->restoreState(m_state); + + Kinetics* kin; + int m, n, npar, ploc; + double mult; + + // process sensitivity parameters + if (params) { + + npar = m_pnum.size(); + for (n = 0; n < npar; n++) { + mult = m_kin->multiplier(m_pnum[n]); + m_kin->setMultiplier(m_pnum[n], mult*params[n]); + } + ploc = npar; + for (m = 0; m < m_nwalls; m++) { + if (m_nsens_wall[m] > 0) { + m_wall[m]->setSensitivityParameters(m_lr[m], params + ploc); + ploc += m_nsens_wall[m]; + } + } + } + + m_vdot = 0.0; + m_Q = 0.0; + + // compute wall terms + doublereal rs0, sum, wallarea; + + SurfPhase* surf; + int lr, ns, loc = m_nsp+2, surfloc; + fill(m_sdot.begin(), m_sdot.end(), 0.0); + for (i = 0; i < m_nwalls; i++) { + lr = 1 - 2*m_lr[i]; + m_Q += lr*m_wall[i]->Q(time); + kin = m_wall[i]->kinetics(m_lr[i]); + surf = m_wall[i]->surface(m_lr[i]); + if (surf && kin) { + rs0 = 1.0/surf->siteDensity(); + nk = surf->nSpecies(); + sum = 0.0; + surf->setTemperature(m_state[0]); + m_wall[i]->syncCoverages(m_lr[i]); + kin->getNetProductionRates(DATA_PTR(m_work)); + ns = kin->surfacePhaseIndex(); + surfloc = kin->kineticsSpeciesIndex(0,ns); + for (k = 1; k < nk; k++) { + ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k); + sum -= ydot[loc + k]; + } + ydot[loc] = sum; + loc += nk; + + wallarea = m_wall[i]->area(); + for (k = 0; k < m_nsp; k++) { + m_sdot[k] += m_work[k]*wallarea; + } + } + } + + // dummy equation + ydot[1] = 0.0; + + /* species equations + * Equation is: + * \dot M_k = \hat W_k \dot\omega_k + \dot m_{in} Y_{k,in} + * - \dot m_{out} Y_{k} + A \dot s_k. + */ + const doublereal* mw = DATA_PTR(m_thermo->molecularWeights()); + if (m_chem) { + m_kin->getNetProductionRates(ydot+2); // "omega dot" + } + else { + fill(ydot + 2, ydot + 2 + m_nsp, 0.0); + } + for (n = 0; n < m_nsp; n++) { + ydot[n+2] *= m_vol; // moles/s/m^3 -> moles/s + ydot[n+2] += m_sdot[n]; + ydot[n+2] *= mw[n]; + } + + + /* + * Energy equation. + * \f[ + * \dot U = -P\dot V + A \dot q + \dot m_{in} h_{in} + * - \dot m_{out} h. + * \f] + */ + if (m_energy) { + ydot[0] = - m_Q; + } + else { + ydot[0] = 0.0; + } + + // add terms for open system + if (m_open) { + + const doublereal* mf = m_thermo->massFractions(); + doublereal enthalpy = m_thermo->enthalpy_mass(); + + // outlets + + int n; + doublereal mdot_out; + for (i = 0; i < m_nOutlets; i++) { + mdot_out = m_outlet[i]->massFlowRate(time); + for (n = 0; n < m_nsp; n++) { + ydot[2+n] -= mdot_out * mf[n]; + } + if (m_energy) { + ydot[0] -= mdot_out * enthalpy; + } + } + + + // inlets + + doublereal mdot_in; + for (i = 0; i < m_nInlets; i++) { + mdot_in = m_inlet[i]->massFlowRate(time); + for (n = 0; n < m_nsp; n++) { + ydot[2+n] += m_inlet[i]->outletSpeciesMassFlowRate(n); + } + if (m_energy) { + ydot[0] += mdot_in * m_inlet[i]->enthalpy_mass(); + } + } + } + + // reset sensitivity parameters + if (params) { + npar = m_pnum.size(); + for (n = 0; n < npar; n++) { + mult = m_kin->multiplier(m_pnum[n]); + m_kin->setMultiplier(m_pnum[n], mult/params[n]); + } + ploc = npar; + for (m = 0; m < m_nwalls; m++) { + if (m_nsens_wall[m] > 0) { + m_wall[m]->resetSensitivityParameters(m_lr[m]); + ploc += m_nsens_wall[m]; + } + } + } + } + + int ConstPressureReactor::componentIndex(string nm) const { + if (nm == "H") return 0; + if (nm == "V") return 1; + // check for a gas species name + int k = m_thermo->speciesIndex(nm); + if (k >= 0) return k + 2; + + // check for a wall species + int walloffset = 0, kp = 0; + thermo_t* th; + for (int m = 0; m < m_nwalls; m++) { + if (m_wall[m]->kinetics(m_lr[m])) { + kp = m_wall[m]->kinetics(m_lr[m])->reactionPhaseIndex(); + th = &m_wall[m]->kinetics(m_lr[m])->thermo(kp); + k = th->speciesIndex(nm); + if (k >= 0) { + return k + 2 + m_nsp + walloffset; + } + else { + walloffset += th->nSpecies(); + } + } + } + return -1; + } + +} diff --git a/Cantera/src/zeroD/ConstPressureReactor.h b/Cantera/src/zeroD/ConstPressureReactor.h new file mode 100644 index 000000000..74624bf47 --- /dev/null +++ b/Cantera/src/zeroD/ConstPressureReactor.h @@ -0,0 +1,73 @@ +/** + * @file Reactor.h + * + * $Author$ + * $Revision$ + * $Date$ + */ + +// Copyright 2001 California Institute of Technology + +#ifndef CT_CONSTP_REACTOR_H +#define CT_CONSTP_REACTOR_H + +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "Reactor.h" + +namespace CanteraZeroD { + + /** + * Class ConstPressureReactor is a class for constant-pressure + * reactors. The reactor may have an arbitrary number of inlets + * and outlets, each of which may be connected to a "flow device" + * such as a mass flow controller, a pressure regulator, + * etc. Additional reactors may be connected to the other end of + * the flow device, allowing construction of arbitrary reactor + * networks. + * + */ + class ConstPressureReactor : public Reactor { + + public: + + /** + * Default constructor. + */ + ConstPressureReactor(); + + /** + * Destructor. Deletes the integrator. + */ + virtual ~ConstPressureReactor(){} + + + virtual int type() const { return ConstPressureReactorType; } + + //----------------------------------------------------- + + virtual int neq() { return m_nv; } + + virtual void getInitialConditions(doublereal t0, size_t leny, + doublereal* y); + + virtual void initialize(doublereal t0 = 0.0); + virtual void evalEqs(doublereal t, doublereal* y, + doublereal* ydot, doublereal* params); + + virtual void updateState(doublereal* y); + + virtual int componentIndex(string nm) const; + + protected: + + private: + + }; +} + +#endif + diff --git a/Cantera/src/zeroD/FlowReactor.cpp b/Cantera/src/zeroD/FlowReactor.cpp index 4b1b1a507..26c3ea778 100644 --- a/Cantera/src/zeroD/FlowReactor.cpp +++ b/Cantera/src/zeroD/FlowReactor.cpp @@ -1,5 +1,5 @@ /** - * @file ReactorZND.cpp + * @file FlowReactor.cpp * * A zero-dimensional reactor */ @@ -27,14 +27,14 @@ namespace CanteraZeroD { void FlowReactor::getInitialConditions(double t0, size_t leny, double* y) { m_init = true; - if (m_mix == 0) { + if (m_thermo == 0) { writelog("Error: reactor is empty.\n"); return; } m_time = t0; - m_mix->restoreState(m_state); + m_thermo->restoreState(m_state); - m_mix->getMassFractions(y+2); + m_thermo->getMassFractions(y+2); y[0] = 0.0; // distance @@ -46,7 +46,7 @@ namespace CanteraZeroD { * Must be called before calling method 'advance' */ void FlowReactor::initialize(doublereal t0) { - m_mix->restoreState(m_state); + m_thermo->restoreState(m_state); m_nv = m_nsp + 2; m_init = true; } @@ -58,7 +58,7 @@ namespace CanteraZeroD { m_speed = y[1]; doublereal* mss = y + 2; // doublereal mass = accumulate(y+2, y+2+m_nsp, 0.0); - m_mix->setMassFractions(mss); + m_thermo->setMassFractions(mss); doublereal rho = m_rho0 * m_speed0/m_speed; @@ -74,7 +74,7 @@ namespace CanteraZeroD { else { m_thermo->setState_TP(m_T, pmom); } - m_mix->saveState(m_state); + m_thermo->saveState(m_state); } @@ -85,10 +85,9 @@ namespace CanteraZeroD { doublereal* ydot, doublereal* params) { m_time = time; - m_mix->restoreState(m_state); + m_thermo->restoreState(m_state); double mult; - Kinetics* kin; int n, npar; // process sensitivity parameters @@ -105,10 +104,10 @@ namespace CanteraZeroD { // speed equation. Set m_fctr to a large value, so that rho*u is // held fixed - ydot[1] = m_fctr*(m_speed0 - m_mix->density()*m_speed/m_rho0); + ydot[1] = m_fctr*(m_speed0 - m_thermo->density()*m_speed/m_rho0); /* species equations */ - const doublereal* mw = DATA_PTR(m_mix->molecularWeights()); + const doublereal* mw = DATA_PTR(m_thermo->molecularWeights()); if (m_chem) { m_kin->getNetProductionRates(ydot+2); // "omega dot" @@ -116,7 +115,7 @@ namespace CanteraZeroD { else { fill(ydot + 2, ydot + 2 + m_nsp, 0.0); } - doublereal rrho = 1.0/m_mix->density(); + doublereal rrho = 1.0/m_thermo->density(); for (n = 0; n < m_nsp; n++) { ydot[n+2] *= mw[n]*rrho; } @@ -138,7 +137,7 @@ namespace CanteraZeroD { if (nm == "X") return 0; if (nm == "U") return 1; // check for a gas species name - int k = m_mix->speciesIndex(nm); + int k = m_thermo->speciesIndex(nm); if (k >= 0) return k + 2; else return -1; } diff --git a/Cantera/src/zeroD/Reactor.cpp b/Cantera/src/zeroD/Reactor.cpp index a6cb7cb68..26da3f3e2 100644 --- a/Cantera/src/zeroD/Reactor.cpp +++ b/Cantera/src/zeroD/Reactor.cpp @@ -26,9 +26,6 @@ namespace CanteraZeroD { doublereal quadInterp(doublereal x0, doublereal* x, doublereal* y); Reactor::Reactor() : ReactorBase(), -#ifdef INCL_REACTOR_INTEG - FuncEval(), -#endif m_kin(0), m_temp_atol(1.e-11), m_maxstep(0.0), @@ -37,36 +34,26 @@ namespace CanteraZeroD { m_rtol(1.e-9), m_chem(true), m_energy(true), m_nsens(-1) - { -#ifdef INCL_REACTOR_INTEG - m_integ = new CVodeInt; - - // use backward differencing, with a full Jacobian computed - // numerically, and use a Newton linear iterator - m_integ->setMethod(BDF_Method); - m_integ->setProblemType(DENSE + NOJAC); - m_integ->setIterator(Newton_Iter); -#endif - } + {} // overloaded method of FuncEval. Called by the integrator to // get the initial conditions. void Reactor::getInitialConditions(double t0, size_t leny, double* y) { m_init = true; - if (m_mix == 0) { + if (m_thermo == 0) { cout << "Error: reactor is empty." << endl; return; } m_time = t0; - m_mix->restoreState(m_state); + m_thermo->restoreState(m_state); // total mass - doublereal mass = m_mix->density() * m_vol; + doublereal mass = m_thermo->density() * m_vol; // set components y + 2 ... y + K + 1 to the // mass M_k of each species - m_mix->getMassFractions(y+2); + m_thermo->getMassFractions(y+2); scale(y + 2, y + m_nsp + 2, y + 2, mass); // set the first component to the total internal @@ -94,19 +81,13 @@ namespace CanteraZeroD { * Must be called before calling method 'advance' */ void Reactor::initialize(doublereal t0) { - m_mix->restoreState(m_state); + m_thermo->restoreState(m_state); m_sdot.resize(m_nsp, 0.0); m_nv = m_nsp + 2; for (int w = 0; w < m_nwalls; w++) if (m_wall[w]->surface(m_lr[w])) m_nv += m_wall[w]->surface(m_lr[w])->nSpecies(); -#ifdef INCL_REACTOR_INTEG - m_atol.resize(neq()); - fill(m_atol.begin(), m_atol.end(), 1.e-15); - m_integ->setTolerances(m_rtol, neq(), DATA_PTR(m_atol)); - m_integ->setMaxStepSize(m_maxstep); - m_integ->initialize(t0, *this); -#endif + m_enthalpy = m_thermo->enthalpy_mass(); m_pressure = m_thermo->pressure(); m_intEnergy = m_thermo->intEnergy_mass(); @@ -146,7 +127,7 @@ namespace CanteraZeroD { void Reactor::updateState(doublereal* y) { - phase_t& mix = *m_mix; // define for readability + phase_t& mix = *m_thermo; // define for readability // The components of y are the total internal energy, // the total volume, and the mass of each species. @@ -158,9 +139,9 @@ namespace CanteraZeroD { m_vol = y[1]; doublereal* mss = y + 2; doublereal mass = accumulate(y+2, y+2+m_nsp, 0.0); - m_mix->setMassFractions(mss); + m_thermo->setMassFractions(mss); - m_mix->setDensity(mass/m_vol); + m_thermo->setDensity(mass/m_vol); doublereal temp = temperature(); mix.setTemperature(temp); @@ -187,16 +168,9 @@ namespace CanteraZeroD { m_enthalpy = m_thermo->enthalpy_mass(); m_pressure = m_thermo->pressure(); m_intEnergy = m_thermo->intEnergy_mass(); - m_mix->saveState(m_state); + m_thermo->saveState(m_state); } -#ifdef INCL_REACTOR_INTEG - void Reactor::eval(doublereal time, doublereal* y, doublereal* ydot) - { - updateState(y); // synchronize the reactor state with y - evalEqs(time, y, ydot); - } -#endif /* * Called by the integrator to evaluate ydot given y at time 'time'. @@ -206,7 +180,7 @@ namespace CanteraZeroD { { int i, k, nk; m_time = time; - m_mix->restoreState(m_state); + m_thermo->restoreState(m_state); Kinetics* kin; int m, n, npar, ploc; @@ -279,7 +253,7 @@ namespace CanteraZeroD { * \dot M_k = \hat W_k \dot\omega_k + \dot m_{in} Y_{k,in} * - \dot m_{out} Y_{k} + A \dot s_k. */ - const doublereal* mw = DATA_PTR(m_mix->molecularWeights()); + const doublereal* mw = DATA_PTR(m_thermo->molecularWeights()); if (m_chem) { m_kin->getNetProductionRates(ydot+2); // "omega dot" } @@ -310,7 +284,7 @@ namespace CanteraZeroD { // add terms for open system if (m_open) { - const doublereal* mf = m_mix->massFractions(); + const doublereal* mf = m_thermo->massFractions(); doublereal enthalpy = m_thermo->enthalpy_mass(); // outlets @@ -374,7 +348,7 @@ namespace CanteraZeroD { if (nm == "U") return 0; if (nm == "V") return 1; // check for a gas species name - int k = m_mix->speciesIndex(nm); + int k = m_thermo->speciesIndex(nm); if (k >= 0) return k + 2; // check for a wall species diff --git a/Cantera/src/zeroD/Reactor.h b/Cantera/src/zeroD/Reactor.h index 542a0ceb7..368681fc6 100644 --- a/Cantera/src/zeroD/Reactor.h +++ b/Cantera/src/zeroD/Reactor.h @@ -17,11 +17,8 @@ #endif #include "ReactorBase.h" -#include "../FuncEval.h" -#include "../Integrator.h" #include "../Kinetics.h" -#undef INCL_REACTOR_INTEG namespace CanteraZeroD { @@ -57,11 +54,9 @@ namespace CanteraZeroD { * flow rate. Class FuncEval is the class used to define a system * of ODE's to be integrated. */ -#ifdef INCL_REACTOR_INTEG - class Reactor : public ReactorBase, public FuncEval { -#else + class Reactor : public ReactorBase { -#endif + public: /** @@ -72,56 +67,10 @@ namespace CanteraZeroD { /** * Destructor. Deletes the integrator. */ - virtual ~Reactor(){ -#ifdef INCL_REACTOR_INTEG - delete m_integ; -#endif -} + virtual ~Reactor(){} virtual int type() const { return ReactorType; } - /** - * Advance the state of the reactor in time. On the first - * call, internal method 'initialize' is called, and the maximum - * integrator step size is set. By default, this is set to - * 'time'. To specify a different maximum step size, precede the - * call to advance with a call to setMaxStep. Note that this - * cannot be reset after advance has been called. - * - * @param time Final time (s). - */ - virtual void advance(doublereal time) { -#ifdef INCL_REACTOR_INTEG - if (!m_init) { - setMaxStep(time); - initialize(); - } - m_integ->integrate(time); - m_time = time; - updateState(m_integ->solution()); - m_mix->saveState(m_state); -#else - throw CanteraError("Reactor::advance", - "Reactor::advance is deprecated. Use ReactorNet::advance"); -#endif - } - - virtual double step(doublereal time) { -#ifdef INCL_REACTOR_INTEG - if (!m_init) { - setMaxStep(time); - initialize(); - } - m_time = m_integ->step(time); - updateState(m_integ->solution()); - m_mix->saveState(m_state); - return m_time; -#else - throw CanteraError("Reactor::step", - "Reactor::step is deprecated. Use ReactorNet::step"); -#endif - } - /** * Insert something into the reactor. The 'something' must * belong to a class that is a subclass of both ThermoPhase @@ -138,15 +87,6 @@ namespace CanteraZeroD { if (m_kin->nReactions() == 0) disableChemistry(); } -#ifdef INCL_REACTOR_INTEG - /** - * Set the maximum step size for integration. - */ - void setMaxStep(doublereal maxstep) { - m_maxstep = maxstep; - } -#endif - void disableChemistry() { m_chem = false; } void enableChemistry() { m_chem = true; } @@ -156,45 +96,22 @@ namespace CanteraZeroD { else m_energy = false; } - - //----------------------------------------------------- - - /** @name References to internal objects */ - //@{ -#ifdef INCL_REACTOR_INTEG - /// Return a reference to the integrator. - Integrator& integrator() { return *m_integ; } - - //@} -#endif - - //----------------------------------------------------- - // overloaded methods of class FuncEval virtual int neq() { return m_nv; } -#ifdef INCL_REACTOR_INTEG - virtual void eval(doublereal t, doublereal* y, doublereal* ydot); -#endif + virtual void getInitialConditions(doublereal t0, size_t leny, doublereal* y); - - //----------------------------------------------------- - virtual void initialize(doublereal t0 = 0.0); - virtual void evalEqs(doublereal t, doublereal* y, doublereal* ydot, doublereal* params); + virtual void evalEqs(doublereal t, doublereal* y, + doublereal* ydot, doublereal* params); /** * Set the mixture to a state consistent with solution * vector y. */ - virtual void updateState(doublereal* y); - // virtual void addSensitivityParam(int stype, int i); - - - virtual int nSensParams(); virtual void addSensitivityReaction(int rxn); @@ -206,9 +123,7 @@ namespace CanteraZeroD { protected: Kinetics* m_kin; -#ifdef INCL_REACTOR_INTEG - Integrator* m_integ; // pointer to integrator -#endif + doublereal m_temp_atol; // tolerance on T doublereal m_maxstep; // max step size doublereal m_vdot, m_Q; diff --git a/Cantera/src/zeroD/ReactorBase.cpp b/Cantera/src/zeroD/ReactorBase.cpp index 41621e394..880755176 100644 --- a/Cantera/src/zeroD/ReactorBase.cpp +++ b/Cantera/src/zeroD/ReactorBase.cpp @@ -21,7 +21,7 @@ namespace CanteraZeroD { ReactorBase::ReactorBase(string name) : m_nsp(0), - m_mix(0), + m_thermo(0), m_time(0.0), m_vol(1.0), m_vol0(1.0), @@ -38,7 +38,7 @@ namespace CanteraZeroD { } // void ReactorBase::resetState() { -// m_mix->saveState(m_state); +// m_thermo->saveState(m_state); // m_enthalpy = m_thermo->enthalpy_mass(); // m_intEnergy = m_thermo->intEnergy_mass(); // m_pressure = m_thermo->pressure(); @@ -46,10 +46,9 @@ namespace CanteraZeroD { // } void ReactorBase::setThermoMgr(thermo_t& thermo){ - m_mix = &thermo; m_thermo = &thermo; - m_nsp = m_mix->nSpecies(); - m_mix->saveState(m_state); + m_nsp = m_thermo->nSpecies(); + m_thermo->saveState(m_state); m_enthalpy = m_thermo->enthalpy_mass(); m_intEnergy = m_thermo->intEnergy_mass(); m_pressure = m_thermo->pressure(); diff --git a/Cantera/src/zeroD/ReactorBase.h b/Cantera/src/zeroD/ReactorBase.h index f6073a3f7..f1c27f42c 100644 --- a/Cantera/src/zeroD/ReactorBase.h +++ b/Cantera/src/zeroD/ReactorBase.h @@ -33,6 +33,7 @@ namespace CanteraZeroD { const int ReservoirType = 1; const int ReactorType = 2; const int FlowReactorType = 3; + const int ConstPressureReactorType = 4; /** * Base class for stirred reactors. @@ -112,9 +113,9 @@ namespace CanteraZeroD { void resetState(); /// return a reference to the contents. - thermo_t& contents() { return *m_mix; } + thermo_t& contents() { return *m_thermo; } - const thermo_t& contents() const { return *m_mix; } + const thermo_t& contents() const { return *m_thermo; } doublereal residenceTime(); @@ -146,12 +147,9 @@ namespace CanteraZeroD { return 1; } - // virtual void addSensitivityParam(int stype, int i) {} - protected: int m_nsp; - thermo_t* m_mix; thermo_t* m_thermo; doublereal m_time; doublereal m_vol, m_vol0;