added constant pressure reactor

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;
+    }
+
+}

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
+

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;
     }

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

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(){ 
+        virtual ~Reactor(){}
-#ifdef INCL_REACTOR_INTEG
-            delete m_integ; 
-#endif
-}
         
         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;

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_nsp = m_thermo->nSpecies();
-        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();

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;