*** empty log message ***

Cantera/clib/src/ctfunc.cpp

@@ -59,6 +59,9 @@ extern "C" {
                         "not enough Arrhenius coefficients");
                 r = new Arrhenius1(n, params);
             }
+            else if (type == PeriodicFuncType) {
+                r = new PeriodicFunc(*_func(n), params[0]);
+            }
             else if (type == SumFuncType) {
                 r = new Func1Sum(*_func(n), *_func(m));
             }

Cantera/clib/src/ctreactor.cpp

@@ -94,10 +94,10 @@ extern "C" {
         return 0;
     }
 
-    int DLL_EXPORT reactor_setInitialTime(int i, double t) {
+    //int DLL_EXPORT reactor_setInitialTime(int i, double t) {
-        _reactor(i)->setInitialTime(t);
+    //    _reactor(i)->setInitialTime(t);
-        return 0;
+    //    return 0;
-    }
+    //}
 
     int DLL_EXPORT reactor_setThermoMgr(int i, int n) {
         _reactor(i)->setThermoMgr(*_th(n));
@@ -111,17 +111,17 @@ extern "C" {
         return 0;
     }
 
-    int DLL_EXPORT reactor_advance(int i, double t) {
+    //int DLL_EXPORT reactor_advance(int i, double t) {
-        try {
+    //    try {
-            _reactor(i)->advance(t);
+    //        _reactor(i)->advance(t);
-            return 0;
+    //        return 0;
-        }
+    //    }
-        catch (CanteraError) {return -1;}
+    //    catch (CanteraError) {return -1;}
-    }
+    //}
 
-    double DLL_EXPORT reactor_step(int i, double t) {
+    //    double DLL_EXPORT reactor_step(int i, double t) {
-        return _reactor(i)->step(t);
+    //    return _reactor(i)->step(t);
-    }
+    //}
 
     double DLL_EXPORT reactor_time(int i) {
         return _reactor(i)->time();
@@ -165,48 +165,6 @@ extern "C" {
         return 0;
     }
 
-//     int DLL_EXPORT reactor_setArea(int i, double a) {
-//         reactor_t* r = _reactor(i);
-//         if (r->type() == ReactorType) ((Reactor*)r)->setArea(a);
-//         return 0;
-//     }
-
-//     int DLL_EXPORT reactor_setExtTemp(int i, double t) {
-//         reactor_t* r = _reactor(i);
-//         if (r->type() == ReactorType) ((Reactor*)r)->setExtTemp(t);
-//         return 0;
-//     }
-
-//     int DLL_EXPORT reactor_setExtRadTemp(int i, double t) {
-//         reactor_t* r = _reactor(i);
-//         if (r->type() == ReactorType) ((Reactor*)r)->setExtRadTemp(t);
-//         return 0;
-//     }
-
-//     int DLL_EXPORT reactor_setVDotCoeff(int i, double v) {
-//         reactor_t* r = _reactor(i);
-//         if (r->type() == ReactorType) ((Reactor*)r)->setVDotCoeff(v);
-//         return 0;
-//     }
-
-//     int DLL_EXPORT reactor_setHeatTransferCoeff(int i, double h) {
-//         reactor_t* r = _reactor(i);
-//         if (r->type() == ReactorType) ((Reactor*)r)->setHeatTransferCoeff(h);
-//         return 0;
-//     }
-
-//     int DLL_EXPORT reactor_setEmissivity(int i, double eps) {
-//         reactor_t* r = _reactor(i);
-//         if (r->type() == ReactorType) ((Reactor*)r)->setEmissivity(eps);
-//         return 0;
-//     }
-
-//     int DLL_EXPORT reactor_setExtPressure(int i, double p) {
-//         reactor_t* r = _reactor(i);
-//         if (r->type() == ReactorType) ((Reactor*)r)->setExtPressure(p);
-//         return 0;
-//     }
-
 
     // reactor networks
 
@@ -248,11 +206,16 @@ extern "C" {
             _reactornet(i)->advance(t);
             return 0;
         }
-        catch (CanteraError) {return -1;}
+        catch (...) {return -1;}
     }
 
     double DLL_EXPORT reactornet_step(int i, double t) {
-        return _reactornet(i)->step(t);
+        try {
+            return _reactornet(i)->step(t);
+        }
+        catch (...) {
+            return DERR;
+        }
     }
 
 
@@ -263,8 +226,8 @@ extern "C" {
         switch (type) {
         case MFC_Type:
             r = new MassFlowController(); break;
-        case PressureReg_Type:
+        case PressureController_Type:
-            r = new PressureRegulator(); break;
+            r = new PressureController(); break;
         case Valve_Type:
             r = new Valve(); break;
         default:
@@ -278,14 +241,6 @@ extern "C" {
         return 0;
     }
 
-    int DLL_EXPORT flowdev_copy(int i) {
-        return Cabinet<flowdev_t>::cabinet()->newCopy(i);
-    }
-
-    int DLL_EXPORT flowdev_assign(int i, int j) {
-        return Cabinet<flowdev_t>::cabinet()->assign(i,j);
-    }
-
     int DLL_EXPORT flowdev_install(int i, int n, int m) {
         try {
             bool ok = _flowdev(i)->install(*_reactor(n), *_reactor(m) );
@@ -297,26 +252,19 @@ extern "C" {
         }
     }
 
-    double DLL_EXPORT flowdev_massFlowRate(int i) {
+    int DLL_EXPORT flowdev_setMaster(int i, int n) {
-        return _flowdev(i)->massFlowRate();
+        if (_flowdev(i)->type() == PressureController_Type) {
-    }
+            ((PressureController*)_flowdev(i))->setMaster(_flowdev(n));
-
+        }
-    double DLL_EXPORT flowdev_setpoint(int i) {
-        return _flowdev(i)->setpoint();
-    }
-
-    int DLL_EXPORT flowdev_setSetpoint(int i, double v) {
-        _flowdev(i)->setSetpoint(v);
         return 0;
     }
 
-    int DLL_EXPORT flowdev_setGains(int i, int n, double* gains) {
+    double DLL_EXPORT flowdev_massFlowRate(int i, double time) {
-        _flowdev(i)->setGains(n, gains);
+        return _flowdev(i)->massFlowRate(time);
-        return 0;
     }
 
-    int DLL_EXPORT flowdev_getGains(int i, int n, double* gains) {
+    int DLL_EXPORT flowdev_setMassFlowRate(int i, double mdot) {
-        _flowdev(i)->getGains(n, gains);
+        _flowdev(i)->setMassFlowRate(mdot);
         return 0;
     }
 
@@ -330,20 +278,6 @@ extern "C" {
         return 0;
     }
 
-    int DLL_EXPORT flowdev_reset(int i) {
-        _flowdev(i)->reset();
-        return 0;
-    }
-
-    int DLL_EXPORT flowdev_update(int i) {
-        _flowdev(i)->update();
-        return 0;
-    }
-
-    double DLL_EXPORT flowdev_maxError(int i) {
-        return _flowdev(i)->maxError();
-    }    
-
     int DLL_EXPORT flowdev_ready(int i) {
         bool ok = _flowdev(i)->ready();
         if (ok) return 1;
@@ -427,8 +361,13 @@ extern "C" {
         return 0;
     }
 
-    int DLL_EXPORT wall_setExpansionRate(int i, int n) {
+    int DLL_EXPORT wall_setVelocity(int i, int n) {
-        _wall(i)->setExpansionRate(_func(n));
+        _wall(i)->setVelocity(_func(n));
+        return 0;
+    }
+
+    int DLL_EXPORT wall_setEmissivity(int i, double epsilon) {
+        _wall(i)->setEmissivity(epsilon);
         return 0;
     }
 

Cantera/clib/src/ctreactor.h

@@ -10,12 +10,12 @@ extern "C" {
     int DLL_IMPORT reactor_copy(int i);
     int DLL_IMPORT reactor_assign(int i, int j);
     int DLL_IMPORT reactor_setInitialVolume(int i, double v);
-    int DLL_IMPORT reactor_setInitialTime(int i, double t);
+    //    int DLL_IMPORT reactor_setInitialTime(int i, double t);
     int DLL_IMPORT reactor_setEnergy(int i, int eflag);
     int DLL_IMPORT reactor_setThermoMgr(int i, int n);
     int DLL_IMPORT reactor_setKineticsMgr(int i, int n);
-    int DLL_IMPORT reactor_advance(int i, double t);
+    //int DLL_IMPORT reactor_advance(int i, double t);
-    double DLL_IMPORT reactor_step(int i, double t);
+    //double DLL_IMPORT reactor_step(int i, double t);
     double DLL_IMPORT reactor_time(int i);
     double DLL_IMPORT reactor_mass(int i);
     double DLL_IMPORT reactor_volume(int i);
@@ -26,15 +26,6 @@ extern "C" {
     double DLL_IMPORT reactor_pressure(int i);
     double DLL_IMPORT reactor_massFraction(int i, int k);
 
-    //int DLL_IMPORT reactor_setArea(int i, double a);
-    //int DLL_IMPORT reactor_setExtTemp(int i, double t);
-    //int DLL_IMPORT reactor_setExtRadTemp(int i, double t);
-    //int DLL_IMPORT reactor_setVDotCoeff(int i, double v);
-    //int DLL_IMPORT reactor_setHeatTransferCoeff(int i, double h);
-    //int DLL_IMPORT reactor_setEmissivity(int i, double eps);
-    //int DLL_IMPORT reactor_setExtPressure(int i, double p);
-    //int DLL_IMPORT reactor_setEnergy(int i, int eflag);
-
     int DLL_IMPORT reactornet_new();
     int DLL_IMPORT reactornet_del(int i);
     int DLL_IMPORT reactornet_copy(int i);
@@ -46,19 +37,12 @@ extern "C" {
 
     int DLL_IMPORT flowdev_new(int type);
     int DLL_IMPORT flowdev_del(int i);
-    //int DLL_IMPORT flowdev_copy(int i);
-    //int DLL_IMPORT flowdev_assign(int i, int j);
     int DLL_IMPORT flowdev_install(int i, int n, int m);
-    double DLL_IMPORT flowdev_massFlowRate(int i);
+    int DLL_IMPORT flowdev_setMaster(int i, int n);
-    double DLL_IMPORT flowdev_setpoint(int i);
+    double DLL_IMPORT flowdev_massFlowRate(int i, double time);
-    int DLL_IMPORT flowdev_setSetpoint(int i, double v);
+    int DLL_IMPORT flowdev_setMassFlowRate(int i, double mdot);
-    //int DLL_IMPORT flowdev_setGains(int i, int n, double* gains);
-    //int DLL_IMPORT flowdev_getGains(int i, int n, double* gains);
     int DLL_IMPORT flowdev_setParameters(int i, int n, double* v);
     int DLL_IMPORT flowdev_setFunction(int i, int n);
-    //int DLL_IMPORT flowdev_reset(int i);
-    //int DLL_IMPORT flowdev_update(int i);
-    //double DLL_IMPORT flowdev_maxError(int i);
     int DLL_IMPORT flowdev_ready(int i);
 
      int DLL_IMPORT wall_new(int type);
@@ -75,7 +59,8 @@ extern "C" {
      int DLL_IMPORT wall_setHeatTransferCoeff(int i, double u);
      int DLL_IMPORT wall_setHeatFlux(int i, int n);
      int DLL_IMPORT wall_setExpansionRateCoeff(int i, double k);
-     int DLL_IMPORT wall_setExpansionRate(int i, int n);
+     int DLL_IMPORT wall_setVelocity(int i, int n);
+    int DLL_IMPORT wall_setEmissivity(int i, double epsilon);
      int DLL_IMPORT wall_ready(int i);
 
 }

Cantera/matlab/Makefile.in

@@ -18,6 +18,7 @@ SRCS = cantera/private/ctmethods.cpp \
        cantera/private/kineticsmethods.cpp \
        cantera/private/transportmethods.cpp \
        cantera/private/reactormethods.cpp \
+       cantera/private/reactornetmethods.cpp \
        cantera/private/wallmethods.cpp \
        cantera/private/flowdevicemethods.cpp \
        cantera/private/onedimmethods.cpp \
@@ -60,6 +61,7 @@ install:
 	@INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/@Solution
 	@INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/@XML_Node/private
 	@INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/@Reactor/private
+	@INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/@ReactorNet/private
 	@INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/@Wall/private
 	@INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/@FlowDevice/private
 	@INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/@Func/private
@@ -94,6 +96,10 @@ install:
             @prefix@/matlab/toolbox/cantera/cantera/@Reactor
 	cd cantera/@Reactor/private; @INSTALL@ *.m \
             @prefix@/matlab/toolbox/cantera/cantera/@Reactor/private
+	cd cantera/@ReactorNet; @INSTALL@ *.m \
+            @prefix@/matlab/toolbox/cantera/cantera/@ReactorNet
+	cd cantera/@ReactorNet/private; @INSTALL@ *.m \
+            @prefix@/matlab/toolbox/cantera/cantera/@ReactorNet/private
 	cd cantera/@Wall; @INSTALL@ *.m \
             @prefix@/matlab/toolbox/cantera/cantera/@Wall
 	cd cantera/@Wall/private; @INSTALL@ *.m \

Cantera/matlab/cantera/private/ctmethods.cpp

@@ -24,6 +24,7 @@ const int PHASE_CLASS = 30;
 const int KINETICS_CLASS = 40;
 const int TRANSPORT_CLASS = 50;
 const int REACTOR_CLASS = 60;
+const int REACTORNET_CLASS = 65;
 const int WALL_CLASS = 70;
 const int FLOWDEVICE_CLASS = 80;
 const int ONEDIM_CLASS = 90;
@@ -54,6 +55,9 @@ void transportmethods( int nlhs, mxArray *plhs[], int nrhs,
 void reactormethods( int nlhs, mxArray *plhs[], int nrhs, 
     const mxArray *prhs[] );
 
+void reactornetmethods( int nlhs, mxArray *plhs[], int nrhs, 
+    const mxArray *prhs[] );
+
 void wallmethods( int nlhs, mxArray *plhs[], int nrhs, 
     const mxArray *prhs[] );
 
@@ -93,6 +97,8 @@ extern "C" {
             transportmethods(nlhs, plhs, nrhs, prhs); break;
         case REACTOR_CLASS:
             reactormethods(nlhs, plhs, nrhs, prhs); break;
+        case REACTORNET_CLASS:
+            reactornetmethods(nlhs, plhs, nrhs, prhs); break;
         case WALL_CLASS:
             wallmethods(nlhs, plhs, nrhs, prhs); break;
         case FLOWDEVICE_CLASS:

Cantera/matlab/cantera/private/reactormethods.cpp

@@ -46,18 +46,18 @@
             case 4:
                 iok = reactor_setInitialVolume(i, v);
                 break;
-            case 5:
+                //            case 5:
-                iok = reactor_setInitialTime(i, v);
+                //iok = reactor_setInitialTime(i, v);
-                break;
+                //break;
             case 6:
                 iok = reactor_setThermoMgr(i, int(v));
                 break;
             case 7:
                 iok = reactor_setKineticsMgr(i, int(v));
                 break;
-            case 8:
+                //case 8:
-                iok = reactor_advance(i, v);
+                //iok = reactor_advance(i, v);
-                break;
+                //break;
             case 9:
                 iok = reactor_setEnergy(i, int(v));
                 break;
@@ -76,9 +76,9 @@
 
         else if (job < 40) {
             switch (job) {
-            case 21:
+                //case 21:
-                r = reactor_step(i, v);
+                //r = reactor_step(i, v);
-                break;
+                //break;
             case 22:
                 r = reactor_time(i);
                 break;

Cantera/python/Cantera/Func.py

@@ -13,7 +13,11 @@ import types
 
 
 class Func1:
-    """Base class for functions of one variable."""
+    """A class for functors of one variable.
+
+    A Functor is an object that behaves like a function. Class 'Func1'
+    is the base class from which several functor classes derive. These
+    classes are designed to be used with the Cantera kernel.  """
     
     def __init__(self, typ, n, coeffs=[]):
         self.n = n
@@ -69,44 +73,76 @@ class Func1:
         return RatioFunction(other, self)            
 
     def func_id(self):
+        """Return the integer index used internally to access the kernel-level object."""
         return self._func_id
 
 
 class Polynomial(Func1):
-    """A polynomial. The degree is determined by the number of coefficients
+    """A polynomial.
-    supplied. Examples:
+    Instances of class 'Polynomial' evaluate
-    
+    \f[
-    p1 = Polynomial([1.0, -2.0, 3.0])   #    3t^2 - 2t + 1
+    f(t) = \sum_{n = 0}^N a_n t^n.
-    p2 = Polynomial([6.0, 8.0])         #    8t + 6    
+    \f]
+    The coefficients are supplied as a list, beginning with \f$a_N\f$ and ending with \f$a_0\f$.
+    >>> p1 = Polynomial([1.0, -2.0, 3.0])   #    3t^2 - 2t + 1
+    >>> p2 = Polynomial([6.0, 8.0])         #    8t + 6
     """
     def __init__(self, coeffs=[]):
         Func1.__init__(self, 2, len(coeffs)-1, coeffs)
 
 
 class Gaussian(Func1):
-    """A Gaussian pulse.
+    """A Gaussian pulse. Instances of class 'Gaussian' evaluate 
+    \f[
+    f(t) = A \exp[-(t - t_0) / \tau]
+    \f]
+    where
+    \f[
+    \tau = \frac{\mbox{FWHM}}{2.0\sqrt{\log(2.0)}}
+    \f]
+    Here FWHM denotes the full width at half maximum. 
     """
-    def __init__(self, A = 0.0, t0 = 0.0, FWHM = 0.0):
+    def __init__(self, A = 0.0, t0 = 0.0, FWHM = 1.0):
-        coeffs = array([A, t0, 0.5*FWHM], 'd')
+        """
+
+        A - Peak value.
+
+        t0 - time at which pulse is centered.
+
+        FWHM - full width at half-maximum.
+
+        """
+        coeffs = array([A, t0, FWHM], 'd')
         Func1.__init__(self, 4, 0, coeffs)
+
         
 class Fourier(Func1):
+    """Fourier series. Instances of class 'Fourier' evaluate the Fourier series
+    \f[
+    f(t) = \frac{a_0}{2} + \sum_{n=1}^N [a_n \cos(n\omega t) + b_n \sin(n \omega t)]
+    \f]
+    where
+    \f[
+    a_n = \int_{-\pi/\omega}^{\pi/\omega} f(t) \cos(n \omega t) dt
+    \f]
+    and
+    \f[
+    b_n = \int_{-\pi/\omega}^{\pi/\omega} f(t) \sin(n \omega t) dt.
+    \f]
+    The function \f$ f(t) \f$ must be periodic, with period \f$ T = 2\pi/\omega \f$.
+    >>> coeffs = [(a0, b0), (a1, b1), (a2, b2)]
+    >>> f = Fourier(omega, coeffs)
+    Note that b0 must be specified, but is not
+    used. The value of b0 is arbitrary.
     """
-    Fourier series.
+    def __init__(self, omega, coefficients):
-    
+        """
-    f(t) = a[0]/2 + sum_{i=1}^n [a[i]*cos(n*omega*t) + b[i]*sin(n*omega*t)]
+        omega - fundamental frequency [radians/sec].
-
+        
-    Note that b[0] must be specified for symmetry with 'a', but is not
+        coefficients - List of (a,b) pairs, beginning with \f$n = 0\f$.
-    used.
+        
-
+        """
-    Example:
+        cc = asarray(coefficients,'d')       
-
-    coeffs = [(a0, b0), (a1, b1), (a2, b2)]
-    f = Fourier(omega, coeffs)
-    """
-    
-    def __init__(self, omega, c):
-        cc = asarray(c,'d')       
         n, m = cc.shape
         if m <> 2:
             raise CanteraError('provide (a, b) for each term')
@@ -115,17 +151,21 @@ class Fourier(Func1):
 
 
 class Arrhenius(Func1):
-    """Sum of modified Arrhenius terms.
+    """Sum of modified Arrhenius terms. Instances of class 'Arrhenius' evaluate
-
+    \f[
     f(T) = \sum_{i=1}^n A_n T^{b_n}\exp(-E_n/T)
-
+    \f]
     Example:
 
-    f = Arrhenius([(a0, b0, e0), (a1, b1, e1)])
+    >>> f = Arrhenius([(a0, b0, e0), (a1, b1, e1)])
     
     """
-    def __init__(self, c):
+    def __init__(self, coefficients):
-        cc = asarray(c,'d')
+        """
+        coefficients - sequence of \f$(A, b, E)\f$ triplets.
+        
+        """
+        cc = asarray(coefficients,'d')
         n, m = cc.shape
         if m <> 3:
             raise CanteraError('Three Arrhenius parameters (A, b, E) required.')
@@ -134,35 +174,87 @@ class Arrhenius(Func1):
 
 
 def Const(value):
-    """Constant function."""
+    """Constant function.
+    >>> f = Const(4.0)  # evaluates f(t) = 4.0.
+    """
     return Polynomial([value])
 
 
+class PeriodicFunction(Func1):
+    def __init__(self, func, T):
+        Func1.__init__(self, 50, func._func_id(), array([T],'d'))
+        
 
 # functions that combine two functions
 
 class SumFunction(Func1):
-    """f = f1 + f2"""
+    """Sum of two functions.
+    Instances of class SumFunction evaluate the sum of two supplied functors.
+    It is not necessary to explicitly create an instance of SumFunction, since
+    the addition operator of the base class is overloaded to return a SumFunction
+    instance.
+    >>> f1 = Polynomial([2.0, 1.0])
+    >>> f2 = Polynomial([3.0, -5.0])
+    >>> f3 = f1 + f2     # functor to evaluate (2t + 1) + (3t - 5)
+    In this example, object 'f3' is a functor of class'SumFunction' that calls f1 and f2
+    and returns their sum.
+    """
+    
     def __init__(self, f1, f2):
+        """
+        f1 - first functor.
+        
+        f2 - second functor.
+        """
         self.f1 = f1
         self.f2 = f2
         self.n = -1
-        self._func_id = _cantera.func_newcombo(20, f1.func_id(), f2.func_id())
+        self._func_id = _cantera.func_newcombo(20, f1._func_id(), f2._func_id())
 
 class ProdFunction(Func1):
-    """f = f1 * f2"""  
+    """Product of two functions.
+    Instances of class ProdFunction evaluate the product of two supplied functors.
+    It is not necessary to explicitly create an instance of 'ProdFunction', since 
+    the multiplication operator of the base class is overloaded to return a 'ProdFunction'
+    instance.
+    >>> f1 = Polynomial([2.0, 1.0])
+    >>> f2 = Polynomial([3.0, -5.0])
+    >>> f3 = f1 * f2     # functor to evaluate (2t + 1)*(3t - 5)
+    In this example, object 'f3' is a functor of class'ProdFunction' that calls f1 and f2
+    and returns their product.
+    """    
     def __init__(self, f1, f2):
+        """
+        f1 - first functor.
+        
+        f2 - second functor.
+        """        
         self.f1 = f1
         self.f2 = f2        
         self.n = -1
-        self._func_id = _cantera.func_newcombo(30, f1.func_id(), f2.func_id())
+        self._func_id = _cantera.func_newcombo(30, f1._func_id(), f2._func_id())
 
 class RatioFunction(Func1):
-    """f = f1 / f2"""    
+    """Ratio of two functions.
+    Instances of class RatioFunction evaluate the ratio of two supplied functors.
+    It  is not necessary to explicitly create an instance of 'RatioFunction', since
+    the division operator of the base class is overloaded to return a RatioFunction
+    instance.
+    >>> f1 = Polynomial([2.0, 1.0])
+    >>> f2 = Polynomial([3.0, -5.0])
+    >>> f3 = f1 / f2     # functor to evaluate (2t + 1)/(3t - 5)
+    In this example, object 'f3' is a functor of class'RatioFunction' that calls f1 and f2
+    and returns their ratio.
+    """        
     def __init__(self, f1, f2):
+        """
+        f1 - first functor.
+        
+        f2 - second functor.
+        """                
         self.f1 = f1
         self.f2 = f2        
         self.n = -1
-        self._func_id = _cantera.func_newcombo(40, f1.func_id(), f2.func_id())
+        self._func_id = _cantera.func_newcombo(40, f1._func_id(), f2._func_id())
     
         

Cantera/python/Cantera/Phase.py

@@ -36,32 +36,29 @@ class Phase:
     def phase_id(self):
         return self._phase_id
     
-##     def addElement(self, symbol, atomicWeight):
-##         """Add an element."""
-##         _cantera.phase_addelement(self._phase_id,symbol,atomicWeight)
-
-##     def addSpecies(self, name, atoms = [],
-##                    thermoType = 'poly',
-##                    thermoCoeffs = []):
-##         _cantera.phase_addSpecies(self._phase_id, name, array(atoms,'d'),
-##                            thermoType, array(thermoCoeffs,'d'),
-##                            minTemp, maxTemp, refPressure)
-    
     def nElements(self):
         """Number of elements."""
         return _cantera.phase_nelements(self._phase_id)
 
-    def atomicWeights(self):
+    def atomicWeights(self, elements = []):
         """Array of element molar masses [kg/kmol]."""
-        return _cantera.phase_getarray(self._phase_id,1)
+        atw = _cantera.phase_getarray(self._phase_id,1)
-
+        if elements:
+            ae = []
+            m = 0
+            for e in elements:
+                m = self.elementIndex(e)
+                ae.append(atw[m])
+            return Numeric.asarray(ae)
+        else:
+            return atw
+        
     def nSpecies(self):
         """Number of species."""
         return _cantera.phase_nspecies(self._phase_id)
 
     def nAtoms(self, species = -1, element = -1):
         """Number of atoms of element 'element' in species 'species'.
-
         The  element and species may be specified by name or by number."""
         try:
             m = self.elementIndex(element)
@@ -97,18 +94,20 @@ class Phase:
         """Mean molar mass [kg/kmol]."""
         return _cantera.phase_meanmolwt(self._phase_id)
 
-    def molecularWeights(self):
+    def molarMasses(self, species = []):
+        """Array of species molar masses [kg/kmol]."""
+        mm  = _cantera.phase_getarray(self._phase_id,22)
+        return self.selectSpecies(mm, species)
+
+    def molecularWeights(self, species = []):
         """Array of species molar masses [kg/kmol].
         DEPRECATED: use molarMasses"""
-        return _cantera.phase_getarray(self._phase_id,22)
+        return self.molarMasses(species)
 
-    def molarMasses(self):
+    def moleFractions(self, species = []):
-        """Array of species molar masses [kg/kmol]."""
-        return _cantera.phase_getarray(self._phase_id,22)    
-
-    def moleFractions(self):
         """Species mole fraction array."""
-        return _cantera.phase_getarray(self._phase_id,20)
+        x = _cantera.phase_getarray(self._phase_id,20)
+        return self.selectSpecies(x, species)
 
     def moleFraction(self, species=-1):
         """Mole fraction of a species, referenced by name or
@@ -119,12 +118,15 @@ class Phase:
         k = self.speciesIndex(species)
         return _cantera.phase_molefraction(self._phase_id,k)    
 
-    def massFractions(self):
+
+    def massFractions(self, species = []):
         """Species mass fraction array."""
-        return _cantera.phase_getarray(self._phase_id,21)
+        y = _cantera.phase_getarray(self._phase_id,21)
+        return self.selectSpecies(y, species)
+
 
     def massFraction(self, species=-1):
-        """Mass fraction of one or more species, referenced by name or
+        """Mass fraction of one species, referenced by name or
         index number.
         >>> ph.massFraction(4)
         >>> ph.massFraction('CH4')
@@ -219,7 +221,7 @@ class Phase:
         if type(x) == types.StringType:
             _cantera.phase_setstring(self._phase_id,2,x)
         else:        
-            _cantera.phase_setarray(self._phase_id,2,norm,x)        
+            _cantera.phase_setarray(self._phase_id,2,norm,Numeric.asarray(x))        
         
     def setState_TRX(self, t, rho, x):
         """Set the temperature, density, and mole fractions."""
@@ -230,7 +232,7 @@ class Phase:
     def setState_TRY(self, t, rho, y):
         """Set the temperature, density, and mass fractions."""        
         self.setTemperature(t)
-        self.setMassFractions(x)
+        self.setMassFractions(y)
         self.setDensity(rho)        
 
     def setState_TR(self, t, rho):
@@ -238,3 +240,14 @@ class Phase:
         self.setTemperature(t)
         self.setDensity(rho)
     
+    def selectSpecies(self, f, sp):
+        if sp:
+            fs = []
+            k = 0
+            for s in sp:
+                k = self.speciesIndex(s)
+                fs.append(f[k])
+            return Numeric.asarray(fs)
+        else:
+            return f
+        

Cantera/python/Cantera/Reactor.py

@@ -20,14 +20,23 @@ class ReactorBase:
         2 = Reservoir).
         """
         self.__reactor_id = _cantera.reactor_new(type)
+        self._type = type
         self._inlets = []
         self._outlets = []
         self._walls = []
+        self._reservoirs = []
         self._name = name
         self._verbose = verbose
         self.insert(contents)
         self._setInitialVolume(volume)
         self._setEnergy(energy)
+        if self._verbose:
+            print 'Created '+self._name
+            print '    Volume = ',volume,' m^3'
+            if energy <> 'on':
+                print '    Temperature will be held constant'
+            print '    Initial State:'
+            print contents
 
 
     def __del__(self):
@@ -38,12 +47,14 @@ class ReactorBase:
 
     def __str__(self):
         s = self._name
+        s += ':\n     Volume = '+`self.volume()`
         if self._contents:
-            s += ": \n"+`self._contents`
+            s += "\n"+`self._contents`
         return s
     
     def __repr__(self):
         s = self._name
+        s += ':\n     Volume = '+`self.volume()`        
         if self._contents:
             s += ": \n"+`self._contents`
         return s
@@ -206,27 +217,33 @@ class ReactorBase:
         """
         return self._walls
     
-    def _addInlet(self, inlet):
+    def _addInlet(self, inlet, other):
         """For internal use. Store a reference to 'inlet'
         so that it will not be deleted before this object."""        
         self._inlets.append(inlet)
+        if self._type == 1 and other._type == 2:
+            self._reservoirs.append(other)        
 
-    def _addOutlet(self, outlet):
+    def _addOutlet(self, outlet, other):
         """For internal use. Store a reference to 'outlet'
         so that it will not be deleted before this object."""        
         self._outlets.append(outlet)
+        if self._type == 1 and other._type == 2:
+            self._reservoirs.append(other)        
 
-    def _addWall(self, wall):
+    def _addWall(self, wall, other):
         """For internal use. Store a reference to 'wall'
         so that it will not be deleted before this object."""
         self._walls.append(wall)
+        if self._type == 1 and other._type == 2:
+            self._reservoirs.append(other)
 
     def syncContents(self):
         """Set the state of the object representing the reactor contents
         to the current reactor state.
         >>> r = Reactor(gas)
         >>> (statements that change the state of object 'gas')
-        >>> r.syncContents(self)
+        >>> r.syncContents()
         After this statement, the state of object 'gas' is synchronized
         with the reactor state.
         See 'contents'.
@@ -356,6 +373,7 @@ class Reservoir(ReactorBase):
 
 
 
+
 #------------------ FlowDevice ---------------------------------
 
 class FlowDevice:
@@ -388,21 +406,9 @@ class FlowDevice:
         """
         return _cantera.flowdev_ready(self.__fdev_id)
 
-    def massFlowRate(self):
+    def massFlowRate(self, time = -999.0):
         """Mass flow rate (kg/s). """
-        return _cantera.flowdev_massFlowRate(self.__fdev_id)
+        return _cantera.flowdev_massFlowRate(self.__fdev_id, time)
-
-    def setSetpoint(self, v):
-        """
-        Deprecated. Set the set point.
-        """
-        _cantera.flowdev_setSetpoint(self.__fdev_id, v)
-
-    def setpoint(self):
-        """
-        Deprecated. The setpoint value.
-        """
-        return _cantera.flowdev_setpoint(self.__fdev_id)
 
     def install(self, upstream, downstream):
         """
@@ -413,8 +419,8 @@ class FlowDevice:
         if self._verbose:
             print
             print self._name+': installing between '+upstream.name()+' and '+downstream.name()
-        upstream._addOutlet(self)
+        upstream._addOutlet(self, downstream)
-        downstream._addInlet(self)
+        downstream._addInlet(self, upstream)
         _cantera.flowdev_install(self.__fdev_id, upstream.reactor_id(),
                                   downstream.reactor_id())
     def _setParameters(self, c):
@@ -422,7 +428,12 @@ class FlowDevice:
         n = len(params)
         return _cantera.flowdev_setParameters(self.__fdev_id, n, params)    
 
+    def setFunction(self, f):
+        _cantera.flowdev_setFunction(self.__fdev_id, f.func_id())
 
+    def flowdev_id(self):
+        return self.__fdev_id
+        
 _mfccount = 0
 
 class MassFlowController(FlowDevice):
@@ -487,14 +498,18 @@ class MassFlowController(FlowDevice):
         """
         if self._verbose:
             print self._name+': setting mdot to '+`mdot`+' kg/s'
-        self.setSetpoint(mdot)
+        if type(mdot) == types.InstanceType:
+            self.setFunction(mdot)
+        else:
+            _cantera.flowdev_setMassFlowRate(self.flowdev_id(), mdot)
+
 
     def set(self, mdot = 0.0):
         """Set the mass flow rate [kg/s].
 
         >>> mfc.set(mdot = 0.2)
         """
-        self.setSetpoint(mdot)
+        self._setMassFlowRate(mdot)
 
 
 _valvecount = 0
@@ -535,7 +550,7 @@ class Valve(FlowDevice):
     
     """
     def __init__(self, upstream=None, downstream=None,
-                 name='', Kv = 0.0, verbose=0):
+                 name='', Kv = 0.0, mdot0 = 0.0, verbose=0):
         """
         upstream - upstream reactor or reservoir.
         
@@ -559,19 +574,89 @@ class Valve(FlowDevice):
         FlowDevice.__init__(self,3,name,verbose)
         if upstream and downstream:
             self.install(upstream, downstream)
-        self.setValveCoeff(Kv)
+        self.setValveCoeff(Kv, mdot0)
 
 
-    def setValveCoeff(self, v):
+    def setValveCoeff(self, Kv, mdot0 = 0.0):
         """Set or reset the valve coefficient \f$ K_v \f$."""
-        vv = zeros(1,'d')
+        vv = zeros(2,'d')
-        vv[0] = v
+        vv[0] = Kv
+        vv[1] = mdot0
         if self._verbose:
             print
-            print self._name+': setting valve coefficient to '+`v`+' kg/Pa-s'
+            print self._name+': setting valve coefficient to '+`Kv`+' kg/Pa-s'
         self._setParameters(vv)
 
+    def _setValveCharacteristic(self, f):
+        """Set or reset the valve characteristics.
+        """
+        if type(f) == types.InstanceType:
+            self.setFunction(f)
+        else:
+            raise CanteraError("Wrong type for valve characteristic function.")
 
+    def set(self, Kv = -1.0, mdot = 0.0, F = None):
+        if F:
+            self.setFunction(F)
+        if Kv > 0.0:
+            self.setValveCoeff(Kv, mdot0 = mdot)
+            
+
+
+_pccount = 0
+
+class PressureController(FlowDevice):
+
+    def __init__(self, upstream=None, downstream=None,
+                 name='', master = None, Kv = 0.0, verbose=0):
+        """
+        upstream - upstream reactor or reservoir.
+        
+        downstream - downstream reactor or reservoir.
+        
+        name - name used to identify the valve in output.
+        If no name is specified, it defaults to 'Valve_n', where n is an
+        integer assigned in the order the Valve object
+        was created.
+
+        Kv - the constant in the mass flow rate equation.
+
+        verbose - if set to a positive integer, additional diagnostic
+        information will be printed.
+        
+        """        
+        global _pccount
+        if name == '':
+            name = 'PressureController_'+`_pccount`
+        _pccount += 1
+        FlowDevice.__init__(self,2,name,verbose)
+        if upstream and downstream:
+            self.install(upstream, downstream)
+        self.setPressureCoeff(Kv)
+        self.setMaster(master)
+
+
+    def setPressureCoeff(self, Kv):
+        """Set or reset the pressure coefficient \f$ K_v \f$."""
+        vv = zeros(1,'d')
+        vv[0] = Kv
+        if self._verbose:
+            print
+            print self._name+': setting pressure coefficient to '+`Kv`+' kg/Pa-s'
+        self._setParameters(vv)
+
+    def setMaster(self, master):
+        _cantera.flowdev_setMaster(self.flowdev_id(),
+                                   master.flowdev_id())
+            
+    def set(self, Kv = -1.0, master = None):
+        if master:
+            self.setMaster(master)
+        if Kv > 0.0:
+            self.setPressureCoeff(Kv)
+            
+
+            
 #------------- Wall ---------------------------
 
 _wallcount = 0
@@ -583,7 +668,7 @@ class Wall:
     """
     def __init__(self, left=None, right=None, name = '',
                  A = 1.0, K = 0.0, U = 0.0,
-                 Q = None, Vdot = None,
+                 Q = None, velocity = None,
                  kinetics = [None, None]):
         typ = 0
         self.__wall_id = _cantera.wall_new(typ)
@@ -601,7 +686,7 @@ class Wall:
             raise CanteraError('both left and right reactors must be specified.')
         self.setArea(A)
         self.setExpansionRateCoeff(K)
-        self.setExpansionRate(Vdot)        
+        self.setVelocity(velocity)        
         self.setHeatTransferCoeff(U)
         self.setHeatFlux(Q)
 
@@ -641,6 +726,13 @@ class Wall:
         """
         return _cantera.wall_setHeatTransferCoeff(self.__wall_id, u)
 
+    def setEmissivity(self, epsilon):
+        """
+        Set the emissivity.
+        The radiative heat flux through the wall is computed from
+        \f[ q_r = \epsion \sigma (T_\ell^4 - T_r^4) \f]
+        """
+        
     def setHeatFlux(self, qfunc=None):
         """
         Specify the time-dependent heat flux function [W/m2].
@@ -655,13 +747,13 @@ class Wall:
         resulting from a unit pressure drop."""
         _cantera.wall_setExpansionRateCoeff(self.__wall_id, k)        
         
-    def setExpansionRate(self, vfunc=None):
+    def setVelocity(self, vfunc=None):
         """
-        Specify the volumetric expansion rate function [m^3/s].
+        Specify the velocity function [m/s].
         """
         n = 0
         if vfunc: n = vfunc.func_id()
-        _cantera.wall_setExpansionRate(self.__wall_id, n)
+        _cantera.wall_setVelocity(self.__wall_id, n)
 
     def vdot(self):
         """Rate of volume change [m^3]. A positive value corresponds
@@ -676,8 +768,8 @@ class Wall:
         return _cantera.wall_Q(self.__wall_id)
     
     def install(self, left, right):
-        left._addWall(self)
+        left._addWall(self, right)
-        right._addWall(self)
+        right._addWall(self, left)
         _cantera.wall_install(self.__wall_id, left.reactor_id(),
                                right.reactor_id())
 

Cantera/python/Cantera/ThermoPhase.py

@@ -106,48 +106,53 @@ class ThermoPhase(Phase):
         """ The pressure [Pa]."""
         return _cantera.thermo_getfp(self._phase_id,7)        
 
-    def chemPotentials(self):
+    def chemPotentials(self, species = []):
         """Species chemical potentials.
         
         This method returns an array containing the species
         chemical potentials [J/kmol]. The expressions used to
         compute these depend on the model implemented by the
         underlying kernel thermo manager."""
-        return _cantera.thermo_getarray(self._phase_id,20)
+        mu = _cantera.thermo_getarray(self._phase_id,20)
+        return self.selectSpecies(mu, species)
     
-    def enthalpies_RT(self):
+    def enthalpies_RT(self, species = []):
         """Pure species non-dimensional enthalpies.
         
         This method returns an array containing the pure-species
         standard-state enthalpies divided by RT. For gaseous species,
         these values are ideal gas enthalpies."""
-        return _cantera.thermo_getarray(self._phase_id,23)
+        hrt = _cantera.thermo_getarray(self._phase_id,23)
-
+        return self.selectSpecies(hrt, species)
-    def entropies_R(self):
+    
+    def entropies_R(self, species = []):
         """Pure species non-dimensional entropies.
         
         This method returns an array containing the pure-species
         standard-state entropies divided by R. For gaseous species,
         these values are ideal gas entropies."""
-        return _cantera.thermo_getarray(self._phase_id,24)    
+        sr = _cantera.thermo_getarray(self._phase_id,24)
+        return self.selectSpecies(sr, species)
 
-    def gibbs_RT(self):
+    def gibbs_RT(self, species = []):
         """Pure species non-dimensional Gibbs free energies.
         
         This method returns an array containing the pure-species
         standard-state Gibbs free energies divided by R.
         For gaseous species, these are ideal gas values."""        
-        return (_cantera.thermo_getarray(self._phase_id,23)
+        grt = (_cantera.thermo_getarray(self._phase_id,23)
                 - _cantera.thermo_getarray(self._phase_id,24))
+        return self.selectSpecies(grt, species)
     
-    def cp_R(self):
+    def cp_R(self, species = []):
         """Pure species non-dimensional heat capacities
         at constant pressure.
         
         This method returns an array containing the pure-species
         standard-state heat capacities divided by R. For gaseous
         species, these values are ideal gas heat capacities."""
-        return _cantera.thermo_getarray(self._phase_id,25)
+        cpr = _cantera.thermo_getarray(self._phase_id,25)
+        return self.selectSpecies(cpr, species)
     
 
     def setPressure(self, p):

Cantera/python/src/ctreactor_methods.cpp

@@ -32,17 +32,17 @@ py_reactor_setInitialVolume(PyObject *self, PyObject *args)
     return Py_BuildValue("i",0);
 }
 
-static PyObject*
+// static PyObject*
-py_reactor_setInitialTime(PyObject *self, PyObject *args)
+// py_reactor_setInitialTime(PyObject *self, PyObject *args)
-{
+// {
-    int n;
+//     int n;
-    double t;
+//     double t;
-    if (!PyArg_ParseTuple(args, "id:reactor_setInitialTime", &n, &t))
+//     if (!PyArg_ParseTuple(args, "id:reactor_setInitialTime", &n, &t))
-        return NULL;
+//         return NULL;
-    int iok = reactor_setInitialTime(n, t);
+//     int iok = reactor_setInitialTime(n, t);
-    if (iok < 0) return reportError(iok);
+//     if (iok < 0) return reportError(iok);
-    return Py_BuildValue("i",0);
+//     return Py_BuildValue("i",0);
-}
+// }
 
 static PyObject*
 py_reactor_setEnergy(PyObject *self, PyObject *args)
@@ -79,27 +79,27 @@ py_reactor_setKineticsMgr(PyObject *self, PyObject *args)
     return Py_BuildValue("i",0);
 }
 
-static PyObject*
+// static PyObject*
-py_reactor_advance(PyObject *self, PyObject *args)
+// py_reactor_advance(PyObject *self, PyObject *args)
-{
+// {
-    int n;
+//     int n;
-    double t;
+//     double t;
-    if (!PyArg_ParseTuple(args, "id:reactor_advance", &n, &t))
+//     if (!PyArg_ParseTuple(args, "id:reactor_advance", &n, &t))
-        return NULL;
+//         return NULL;
-    int iok = reactor_advance(n, t);
+//     int iok = reactor_advance(n, t);
-    if (iok < 0) return reportError(iok);
+//     if (iok < 0) return reportError(iok);
-    return Py_BuildValue("i",0);
+//     return Py_BuildValue("i",0);
-}
+// }
 
-static PyObject*
+// static PyObject*
-py_reactor_step(PyObject *self, PyObject *args)
+// py_reactor_step(PyObject *self, PyObject *args)
-{
+// {
-    int n;
+//     int n;
-    double t;
+//     double t;
-    if (!PyArg_ParseTuple(args, "id:reactor_step", &n, &t))
+//     if (!PyArg_ParseTuple(args, "id:reactor_step", &n, &t))
-        return NULL;
+//         return NULL;
-    return Py_BuildValue("d",reactor_step(n, t));
+//     return Py_BuildValue("d",reactor_step(n, t));
-}
+// }
 
 static PyObject*
 py_reactor_time(PyObject *self, PyObject *args)
@@ -224,34 +224,46 @@ py_flowdev_install(PyObject *self, PyObject *args)
     return Py_BuildValue("i",0);
 }
 
+static PyObject*
+py_flowdev_setMaster(PyObject *self, PyObject *args)
+{
+    int n, m;
+    if (!PyArg_ParseTuple(args, "ii:flowdev_setMaster", &n, &m))
+        return NULL;
+    int iok = flowdev_setMaster(n, m);
+    if (iok < 0) return reportError(iok);
+    return Py_BuildValue("i",0);
+}
+
 static PyObject*
 py_flowdev_massFlowRate(PyObject *self, PyObject *args)
 {
     int n;
-    if (!PyArg_ParseTuple(args, "i:flowdev_massFlowRate", &n))
+    double t;
+    if (!PyArg_ParseTuple(args, "id:flowdev_massFlowRate", &n, &t))
         return NULL;
-    double mdot = flowdev_massFlowRate(n);
+    double mdot = flowdev_massFlowRate(n, t);
     return Py_BuildValue("d",mdot);
 }
 
-static PyObject*
+// static PyObject*
-py_flowdev_setpoint(PyObject *self, PyObject *args)
+// py_flowdev_setpoint(PyObject *self, PyObject *args)
-{
+// {
-    int n;
+//     int n;
-    if (!PyArg_ParseTuple(args, "i:flowdev_setpoint", &n))
+//     if (!PyArg_ParseTuple(args, "i:flowdev_setpoint", &n))
-        return NULL;
+//         return NULL;
-    double v = flowdev_setpoint(n);
+//     double v = flowdev_setpoint(n);
-    return Py_BuildValue("d",v);
+//     return Py_BuildValue("d",v);
-}
+// }
 
 static PyObject*
-py_flowdev_setSetpoint(PyObject *self, PyObject *args)
+py_flowdev_setMassFlowRate(PyObject *self, PyObject *args)
 {
     int n;
-    double v;
+    double mdot;
-    if (!PyArg_ParseTuple(args, "id:flowdev_setSetpoint", &n, &v))
+    if (!PyArg_ParseTuple(args, "id:flowdev_setMassFlowRate", &n, &mdot))
         return NULL;
-    int iok = flowdev_setSetpoint(n, v);
+    int iok = flowdev_setMassFlowRate(n, mdot);
     if (iok < 0) return reportError(iok);
     return Py_BuildValue("i",0);
 }
@@ -404,6 +416,18 @@ py_wall_setHeatTransferCoeff(PyObject *self, PyObject *args)
     return Py_BuildValue("i",0);
 }
 
+static PyObject*
+py_wall_setEmissivity(PyObject *self, PyObject *args)
+{
+    int n;
+    double epsilon;
+    if (!PyArg_ParseTuple(args, "id:wall_setEmissivity", &n, &epsilon))
+        return NULL;
+    int iok = wall_setEmissivity(n,epsilon);
+    if (iok < 0) return reportError(iok);
+    return Py_BuildValue("i",0);
+}
+
 static PyObject*
 py_wall_setExpansionRateCoeff(PyObject *self, PyObject *args)
 {
@@ -417,12 +441,12 @@ py_wall_setExpansionRateCoeff(PyObject *self, PyObject *args)
 }
 
 static PyObject*
-py_wall_setExpansionRate(PyObject *self, PyObject *args)
+py_wall_setVelocity(PyObject *self, PyObject *args)
 {
     int n, m;
-    if (!PyArg_ParseTuple(args, "ii:wall_setExpansionRate", &n, &m))
+    if (!PyArg_ParseTuple(args, "ii:wall_setVelocity", &n, &m))
         return NULL;
-    int iok = wall_setExpansionRate(n,m);
+    int iok = wall_setVelocity(n,m);
     if (iok < 0) return reportError(iok);
     return Py_BuildValue("i",0);
 }

Cantera/python/src/methods.h

@@ -192,21 +192,22 @@ static PyMethodDef ct_methods[] = {
     {"bndry_setmdot", py_bndry_setmdot, METH_VARARGS},
 
     {"flowdev_ready", py_flowdev_ready, METH_VARARGS},
-    {"reactor_setInitialTime", py_reactor_setInitialTime, METH_VARARGS},
+    //{"reactor_setInitialTime", py_reactor_setInitialTime, METH_VARARGS},
     {"reactornet_setInitialTime", py_reactornet_setInitialTime, METH_VARARGS},
     {"flowdev_new", py_flowdev_new, METH_VARARGS},
     {"flowdev_massFlowRate", py_flowdev_massFlowRate, METH_VARARGS},
     {"flowdev_del", py_flowdev_del, METH_VARARGS},
-    {"flowdev_setpoint", py_flowdev_setpoint, METH_VARARGS},
+    //    {"flowdev_setpoint", py_flowdev_setpoint, METH_VARARGS},
     {"reactor_temperature", py_reactor_temperature, METH_VARARGS},
-    {"flowdev_setSetpoint", py_flowdev_setSetpoint, METH_VARARGS},
+    {"flowdev_setMassFlowRate", py_flowdev_setMassFlowRate, METH_VARARGS},
     {"flowdev_install", py_flowdev_install, METH_VARARGS},
+    {"flowdev_setMaster", py_flowdev_setMaster, METH_VARARGS},
     {"reactor_setThermoMgr", py_reactor_setThermoMgr, METH_VARARGS},
     {"reactor_setEnergy", py_reactor_setEnergy, METH_VARARGS},
     {"reactor_volume", py_reactor_volume, METH_VARARGS},
     {"reactor_time", py_reactor_time, METH_VARARGS},
-    {"reactor_advance", py_reactor_advance, METH_VARARGS},
+    //    {"reactor_advance", py_reactor_advance, METH_VARARGS},
-    {"reactor_step", py_reactor_step, METH_VARARGS},
+    //{"reactor_step", py_reactor_step, METH_VARARGS},
     {"reactornet_addreactor", py_reactornet_addreactor, METH_VARARGS},
     {"reactornet_advance", py_reactornet_advance, METH_VARARGS},
     {"reactornet_step", py_reactornet_step, METH_VARARGS},
@@ -235,7 +236,8 @@ static PyMethodDef ct_methods[] = {
     {"wall_new", py_wall_new, METH_VARARGS},
     {"wall_vdot", py_wall_vdot, METH_VARARGS},
     {"wall_del", py_wall_del, METH_VARARGS},
-    {"wall_setExpansionRate", py_wall_setExpansionRate, METH_VARARGS},
+    {"wall_setVelocity", py_wall_setVelocity, METH_VARARGS},
+    {"wall_setEmissivity", py_wall_setEmissivity, METH_VARARGS},
     {"wall_setExpansionRateCoeff", py_wall_setExpansionRateCoeff, METH_VARARGS},
     {"wall_ready", py_wall_ready, METH_VARARGS},
 

Cantera/src/Func1.h

@@ -28,7 +28,8 @@ namespace Cantera {
     const int DiffFuncType = 25;
     const int ProdFuncType = 30;
     const int RatioFuncType = 40;
- 
+    const int PeriodicFuncType = 50;
+
     /**
      * Base class for 'functor' classes that evaluate a function of
      * one variable.
@@ -47,10 +48,10 @@ namespace Cantera {
 
     class Gaussian : public Func1 {
     public:
-        Gaussian(double A, double t0, double tau) {
+        Gaussian(double A, double t0, double fwhm) {
             m_A = A;
             m_t0 = t0;
-            m_tau = tau;
+            m_tau = fwhm/(2.0*sqrt(log(2.0)));
         }
         virtual ~Gaussian() {}
         virtual doublereal eval(doublereal t) {
@@ -177,6 +178,29 @@ namespace Cantera {
         vector_fp m_A, m_b, m_E;
     };
 
+    /**
+     *  Periodic function. Takes any function and makes it
+     *  periodic with period T.
+     */
+    class PeriodicFunc : public Func1 {
+    public:
+        PeriodicFunc(Func1& f, doublereal T) {
+            m_func = &f;
+            m_period = T;
+        }
+        virtual ~PeriodicFunc() {}
+        virtual doublereal eval(doublereal t) {
+            int np = int(t/m_period);
+            doublereal time = t - np*m_period;
+            return m_func->eval(time);
+        }
+    protected:
+        Func1* m_func;
+        doublereal m_period;
+
+    private:
+    };
+            
 
     /**
      * Sum of two functions.

Cantera/src/IdealGasPhase.h

@@ -228,6 +228,14 @@ namespace Cantera {
         }
 
 
+        virtual doublereal isothermalCompressibility() {
+            return -1.0/pressure();
+        }
+
+        virtual doublereal thermalExpansionCoeff() {
+            return 1.0/temperature();
+        }
+
         // new methods defined here
 
         const array_fp& enthalpy_RT() const {

Cantera/src/NasaThermo.h

@@ -86,11 +86,7 @@ namespace Cantera {
             vector_fp chigh(7);
             copy(c + 8, c + 15, chigh.begin());
 
-            // Make cp exactly continuous at tmid by offsetting the
+            checkContinuity(tmid, clow, chigh.begin());
-            // high temp fit.
-            doublereal cplow = poly4(tmid, clow+2);
-            doublereal cphigh = poly4(tmid, chigh.begin()+2);
-            chigh[0] += (cplow - cphigh);
 
             m_high[igrp-1].push_back(NasaPoly1(index, tmid, thigh, 
                                          pref, chigh.begin()));
@@ -105,7 +101,6 @@ namespace Cantera {
             m_low_map[index] = &m_low[igrp-1].back();
         }
 
-
         /** 
          * update the properties for only one species.
          */
@@ -193,6 +188,52 @@ namespace Cantera {
         doublereal                         m_p0;
         int                                m_ngroups;
         mutable vector_fp                  m_t;
+
+    private:
+
+
+        /// Check the continuity of properties at the midpoint
+        /// temperature, and adjust the high-T coefficients to
+        /// make the properties exactly continuous at Tmid.
+        void checkContinuity(double tmid, const doublereal* clow,
+            doublereal* chigh) {
+
+            // heat capacity
+            doublereal cplow = poly4(tmid, clow+2);
+            doublereal cphigh = poly4(tmid, chigh+2);
+            doublereal delta = cplow - cphigh;
+            if (fabs(delta/cplow) > 0.001) {
+                writelog("WARNING: discontinuity in cp/R detected at Tmid = "
+                    +fp2str(tmid)+"\n");
+                writelog("  Adjusting high-temperature coefficient to fix.\n");
+            }
+            chigh[2] += cplow - cphigh;
+
+            // enthalpy
+            doublereal hrtlow = enthalpy_RT(tmid, clow);
+            doublereal hrthigh = enthalpy_RT(tmid, chigh);
+            chigh[0] += tmid*(hrtlow - hrthigh);
+
+            // entropy
+            doublereal srlow = entropy_R(tmid, clow);
+            doublereal srhigh = entropy_R(tmid, chigh);
+            chigh[1] += srlow - srhigh;
+        }
+
+        /// for internal use by checkContinuity
+        doublereal enthalpy_RT(double t, const doublereal* c) {
+            return c[2] + 0.5*c[3]*t + OneThird*c[4]*t*t 
+                + 0.25*c[5]*t*t*t + 0.2*c[6]*t*t*t*t
+                + c[0]/t;
+        }
+
+        /// for internal use by checkContinuity
+        doublereal entropy_R(double t, const doublereal* c) {
+            return c[2]*log(t) + c[3]*t + 0.5*c[4]*t*t 
+                + OneThird*c[5]*t*t*t + 0.25*c[6]*t*t*t*t
+                + c[1];
+        }
+
     };
 
 }
@@ -200,7 +241,10 @@ namespace Cantera {
 #endif
 
 // $Log$
-// Revision 1.2  2003-11-01 04:50:35  dggoodwin
+// Revision 1.3  2004-04-22 21:44:36  dggoodwin
+// *** empty log message ***
+//
+// Revision 1.2  2003/11/01 04:50:35  dggoodwin
 // *** empty log message ***
 //
 // Revision 1.1.1.1  2003/04/14 17:57:51  dggoodwin

Cantera/src/PureFluidPhase.h

@@ -149,6 +149,13 @@ namespace Cantera {
             mu[0] = gibbs_mole();
         }
 
+        virtual doublereal isothermalCompressibility() {
+            return m_sub->isothermalCompressibility();
+        }
+
+        virtual doublereal thermalExpansionCoeff() {
+            return m_sub->thermalExpansionCoeff();
+        }
 
         tpx::Substance& TPX_Substance() { return *m_sub; }
 

Cantera/src/ThermoPhase.cpp

@@ -115,13 +115,15 @@ namespace Cantera {
         for (int n = 0; n < 50; n++) {
             dt = (u - intEnergy_mass())/cv_mass();
             if (dt > 100.0) dt = 100.0;
-            else if (dt < -100.0) dt = -100.0; 
+            else if (dt < -100.0) dt = -100.0;
             setTemperature(temperature() + dt);
             if (fabs(dt) < tol) {
                 return;
             }
         }
-        throw CanteraError("setState_UV","no convergence. dt = " + fp2str(dt));
+        throw CanteraError("setState_UV",
+            "no convergence. dt = " + fp2str(dt)+"\n"
+            +"u = "+fp2str(u)+" v = "+fp2str(v)+"\n");
     }
 
     void ThermoPhase::setState_SP(doublereal s, doublereal p, 

Cantera/src/ThermoPhase.h

@@ -569,7 +569,14 @@ namespace Cantera {
          */
         virtual void setParameters(int n, doublereal* c) {}
         
-        
+        virtual doublereal isothermalCompressibility() {
+            err("isothermalCompressibility"); return -1.0;
+        }
+
+        virtual doublereal thermalExpansionCoeff() {
+            err("thermalExpansionCoeff()"); return -1.0;
+        }
+
         //---------------------------------------------------------
         /// @name Critical state properties.
         

Cantera/src/ctexceptions.h

@@ -24,6 +24,7 @@ namespace Cantera {
         CanteraError() {}
         CanteraError(string proc, string msg) {
             setError(proc, msg);
+            writelog("Throwing CanteraError. "+proc+" "+msg+"\n");
             //m_msg = msg;
         }
         virtual ~CanteraError(){}

Cantera/src/zeroD/FlowDevice.cpp

@@ -34,14 +34,16 @@ namespace Cantera {
         return true; 
     }
 
-    void FlowDevice::setFunction(Func1* f) {}
+    void FlowDevice::setFunction(Func1* f) {
+        m_func = f;
+    }
 
 
     /**
      * Mass flow rate of outlet species k.  Returns zero if this
      * species is not present in the upstream mixture.
      */
-    doublereal FlowDevice::massFlowRate(int k) {
+    doublereal FlowDevice::outletSpeciesMassFlowRate(int k) {
         if (k < 0 || k >= m_nspout) return 0.0;
         int ki = m_out2in[k];
         if (ki < 0) return 0.0;

Cantera/src/zeroD/FlowDevice.h

@@ -24,7 +24,7 @@ namespace Cantera {
     class Func1;
 
     const int MFC_Type = 1;
-    const int PressureReg_Type = 2;
+    const int PressureController_Type = 2;
     const int Valve_Type = 3;
 
     /**
@@ -43,79 +43,92 @@ namespace Cantera {
     public:
 
         /// Constructor
-        FlowDevice() : m_mdot(0.0), 
+        FlowDevice() : m_mdot(0.0), m_func(0), m_type(0),
                        m_nspin(0), m_nspout(0),
                        m_in(0), m_out(0) {}
 
         /// Destructor (does nothing)
         virtual ~FlowDevice(){}
 
-        /// Copy constructor.
+ //        /// Copy constructor.
-        FlowDevice(const FlowDevice& a) : m_in(a.m_in), m_out(a.m_out) {}
+//         FlowDevice(const FlowDevice& a) : m_in(a.m_in), m_out(a.m_out) {}
 
-        /// Assignment operator
+//         /// Assignment operator
-        FlowDevice& operator=(const FlowDevice& a) {
+//         FlowDevice& operator=(const FlowDevice& a) {
-            if (this == &a) return *this;
+//             if (this == &a) return *this;
-            m_in = a.m_in;
+//             m_in = a.m_in;
-            m_out = a.m_out;
+//             m_out = a.m_out;
-            return *this;
+//             return *this;
-        }
+//         }
-            
-        /**
-         * Mass flow rate (kg/s). May be overloaded in derived
-         * classes.
-         */
-        virtual doublereal massFlowRate() {return m_mdot;}
-        
-        doublereal massFlowRate(int k);
-        virtual doublereal enthalpy_mass();
 
-        /** 
+        int type() { return m_type; }
-         * Setpoint. Default = 0.0.
-         */
-        virtual doublereal setpoint() { warn("setpoint"); return 0.0; }
-
-        /* Update the internal state, if necessary. By default this method
-         * does nothing, but may be overloaded for devices that have a
-         * state.
-         */
-        virtual void update() {warn("update");}
-
-        /* Reset the device. By default this method does nothing, but
-         * may be overloaded for devices that have a state that depends on
-         * past history.
-         */
-        virtual void reset() {warn("reset");}
 
         /**
-         * Set the setpoint. May be changed at any time. By default, 
+         * Mass flow rate (kg/s).
-         * this does nothing.
          */
-        virtual void setSetpoint(doublereal value) {warn("setSetpoint");}
+        doublereal massFlowRate(double time = -999.0) {
-
+            if (time != -999.0) updateMassFlowRate(time);
-        /**
+            return m_mdot;
-         * Set the controller gains. Returns false if the number of
-         * gains is too small, or if an illegal value is specified.
-         */
-        virtual bool setGains(int n, const doublereal* gains) {
-            warn("setGains");
-            return true;
         }
 
-        /** 
+        // Update the mass flow rate at time 'time'. This must be 
-         * Get the controller gains. Returns false if the 'gains'
+        // overloaded in subclassess to update m_mdot.
-         * array is too small.
+        virtual void updateMassFlowRate(doublereal time) {}
-         */
-        virtual bool getGains(int n, doublereal* gains) {
-            warn("getGains");
-            return true;
-        }
 
-        /** 
+        // mass flow rate of outlet species k
-         * Maximum difference between input and setpoint since 
+        doublereal outletSpeciesMassFlowRate(int k);
-         * last call to 'reset'.
+
-         */
+        // specific enthalpy
-        virtual doublereal maxError() {warn("maxError"); return 0.0;}
+        doublereal enthalpy_mass();
+
+//         /** 
+//          * Setpoint. Default = 0.0.
+//          */
+//         virtual doublereal setpoint() { warn("setpoint"); return 0.0; }
+
+
+//         /* Update the internal state, if necessary. By default this method
+//          * does nothing, but may be overloaded for devices that have a
+//          * state.
+//          */
+//         virtual void update() {warn("update");}
+
+
+//         /* Reset the device. By default this method does nothing, but
+//          * may be overloaded for devices that have a state that depends on
+//          * past history.
+//          */
+//         virtual void reset() {warn("reset");}
+
+//         /**
+//          * Set the setpoint. May be changed at any time. By default, 
+//          * this does nothing.
+//          */
+//         virtual void setSetpoint(doublereal value) {warn("setSetpoint");}
+
+//         /**
+//          * Set the controller gains. Returns false if the number of
+//          * gains is too small, or if an illegal value is specified.
+//          */
+//         virtual bool setGains(int n, const doublereal* gains) {
+//             warn("setGains");
+//             return true;
+//         }
+
+//         /** 
+//          * Get the controller gains. Returns false if the 'gains'
+//          * array is too small.
+//          */
+//         virtual bool getGains(int n, doublereal* gains) {
+//             warn("getGains");
+//             return true;
+//         }
+
+//         /** 
+//          * Maximum difference between input and setpoint since 
+//          * last call to 'reset'.
+//          */
+//         virtual doublereal maxError() {warn("maxError"); return 0.0;}
 
         /**
          * Install a flow device between two reactors.
@@ -125,7 +138,6 @@ namespace Cantera {
         bool install(ReactorBase& in, ReactorBase& out);
 
         virtual bool ready() { return (m_in != 0 && m_out != 0); }
-        doublereal m_mdot;
 
         /// Return a reference to the upstream reactor.
         ReactorBase& in() const { return *m_in; }
@@ -139,11 +151,16 @@ namespace Cantera {
             copy(coeffs, coeffs + n, m_coeffs.begin());
         }
 
-        virtual void setFunction(Func1* f);
+        void setFunction(Func1* f);
+        void setMassFlowRate(doublereal mdot) {m_mdot = mdot;}
+
 
     protected:
 
+        doublereal m_mdot;
+        Func1* m_func;
         vector_fp m_coeffs;
+        int m_type;
 
     private:
 

Cantera/src/zeroD/Reactor.cpp

@@ -26,7 +26,6 @@ namespace Cantera {
     Reactor::Reactor() : ReactorBase(), 
                          FuncEval(),
                          m_kin(0),
-                         m_integ(0),
                          m_temp_atol(1.e-11), 
                          m_maxstep(0.0),
                          m_vdot(0.0), 
@@ -35,13 +34,15 @@ namespace Cantera {
                          m_chem(true),
                          m_energy(true)
     {
+#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);        
+        m_integ->setIterator(Newton_Iter);
+#endif
     }
 
 
@@ -79,7 +80,8 @@ namespace Cantera {
         for (int m = 0; m < m_nwalls; m++) {
             surf = m_wall[m]->surface(m_lr[m]);
             if (surf) {
-                surf->getCoverages(y+loc);
+                m_wall[m]->getCoverages(m_lr[m], y + loc);
+                //surf->getCoverages(y+loc);
                 loc += surf->nSpecies();
             }
         }
@@ -96,12 +98,13 @@ namespace Cantera {
         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(), m_atol.begin());
         m_integ->setMaxStep(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();
@@ -140,33 +143,26 @@ namespace Cantera {
         doublereal* mss = y + 2;
         doublereal mass = accumulate(y+2, y+2+m_nsp, 0.0);
         m_mix->setMassFractions(mss);
+
         m_mix->setDensity(mass/m_vol);
 
         doublereal temp = temperature();
         mix.setTemperature(temp);
 
         if (m_energy) {
-            doublereal u_mass = u/mass;       // specific int. energy
+            m_thermo->setState_UV(u/mass,m_vol/mass); 
-            doublereal delta;
+            temp = mix.temperature(); //mix.setTemperature(temp);
-
-            do {
-                delta = -(m_thermo->intEnergy_mass() 
-                    - u_mass)/m_thermo->cv_mass();
-                temp += delta;
-                mix.setTemperature(temp);
-            }
-            while (fabs(delta) > m_temp_atol);
         }
-        mix.setTemperature(temp);
+        //m_state[0] = temp;
-        m_state[0] = temp;
 
         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) {
-                surf->setTemperature(temp);
+                //                surf->setTemperature(temp);
-                surf->setCoverages(y+loc);
+                //surf->setCoverages(y+loc);
+                m_wall[m]->setCoverages(m_lr[m], y+loc);
                 loc += surf->nSpecies();
             }
         }
@@ -175,7 +171,7 @@ namespace Cantera {
         m_enthalpy = m_thermo->enthalpy_mass();
         m_pressure = m_thermo->pressure();
         m_intEnergy = m_thermo->intEnergy_mass();
-
+        //m_kappa = m_thermo->isothermalCompressibility();
         m_mix->saveState(m_state);
     }
 
@@ -217,6 +213,8 @@ namespace Cantera {
                 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(m_work.begin());
                 ns = kin->surfacePhaseIndex();
                 surfloc = kin->kineticsSpeciesIndex(0,ns);
@@ -281,12 +279,16 @@ namespace Cantera {
             int n;
             doublereal mdot_out;
             for (i = 0; i < m_nOutlets; i++) {
-                mdot_out = m_outlet[i]->massFlowRate();
+                mdot_out = m_outlet[i]->massFlowRate(time);
                 for (n = 0; n < m_nsp; n++) {
                     ydot[2+n] -= mdot_out * mf[n];
                 }
-                if (m_energy)
+                if (m_energy) {
+                    //                    cout << "before = " << ydot[0] << endl;
                     ydot[0] -= mdot_out * enthalpy;
+                    //cout << mdot_out << " " << enthalpy << endl;
+                    //cout << "after = " << ydot[0] << endl;
+                }
             }
 
 
@@ -294,12 +296,13 @@ namespace Cantera {
 
             doublereal mdot_in;
             for (i = 0; i < m_nInlets; i++) {
-                mdot_in = m_inlet[i]->massFlowRate();
+                mdot_in = m_inlet[i]->massFlowRate(time);
                 for (n = 0; n < m_nsp; n++) {
-                    ydot[2+n] += m_inlet[i]->massFlowRate(n);
+                    ydot[2+n] += m_inlet[i]->outletSpeciesMassFlowRate(n);
                 }
-                if (m_energy)
+                if (m_energy) {
                     ydot[0] += mdot_in * m_inlet[i]->enthalpy_mass();
+                }
             }
         }
     }

Cantera/src/zeroD/Reactor.h

@@ -69,7 +69,11 @@ namespace Cantera {
         /**
          * Destructor. Deletes the integrator.
          */
-        virtual ~Reactor(){ delete m_integ; }
+        virtual ~Reactor(){ 
+#ifdef INCL_REACTOR_INTEG
+            delete m_integ; 
+#endif
+}
         
         virtual int type() const { return ReactorType; }
 
@@ -84,6 +88,7 @@ namespace Cantera {
          * @param time Final time (s).
          */
         virtual void advance(doublereal time) {
+#ifdef INCL_REACTOR_INTEG
             if (!m_init) {
                 setMaxStep(time);
                 initialize();
@@ -92,9 +97,14 @@ namespace Cantera {
             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();
@@ -103,6 +113,10 @@ namespace Cantera {
             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
         }
 
         /**
@@ -128,53 +142,6 @@ namespace Cantera {
             m_maxstep = maxstep;
         }
 
- //        /**
-//          * Set the reactor surface area [m$^2$]. Can be changed at any time.
-//          */
-//         void setArea(doublereal area) {
-//             m_area = area;
-//         }
-
-//         /**
-//          * Set the external temperature \f$ T_0 \f$ 
-//          * used for heat loss calculations.
-//          * The heat loss rate is calculated from
-//          * \f[
-//          * \dot Q_{out} = h A (T - T_0) + \epsilon A (T^4 - T_{0,R}^4).
-//          * \f]
-//          * @see setArea, setEmissivity, setExtRadTemp
-//          */
-//         void setExtTemp(doublereal ts) {
-//             m_ext_temp = ts;
-//             if (!m_trad_set) m_ext_temp4 = ts*ts*ts*ts;
-//         }
-
-//         /**
-//          * Set the external temperature for radiation. By default, this
-//          * is the same as the temperature set by setExtTemp. But if 
-//          * setExtRadTemp is called, then subsequent of calls to 
-//          * setExtTemp do not modify the value set here.
-//          */
-//         void setExtRadTemp(doublereal tr) {        
-//             m_ext_temp4 = tr*tr*tr*tr;
-//         }
-
-//         void setHeatTransferCoeff(doublereal h) {
-//             m_h = h;
-//         }
-
-//         void setVDotCoeff(doublereal k) {
-//             m_kv = k;
-//         }
-
-//         void setEmissivity(doublereal emis) {
-//             m_emis = emis;
-//         }
-
-//         void setExtPressure(doublereal p0) {
-//             m_p0 = p0;
-//         }
-
         void disableChemistry() { m_chem = false; }
         void enableChemistry() { m_chem = true; }
 
@@ -208,61 +175,6 @@ namespace Cantera {
         virtual void initialize(doublereal t0 = 0.0);
 	void evalEqs(doublereal t, doublereal* y, doublereal* ydot);
 
-        /**
-         * @name Methods to specify simulation options.
-         * These virtual methods may be overloaded in 
-         * derived classes to implement models for heat gain/loss,
-         * surface chemistry, and compression/expansion.
-         */
-        //@{
-
-        /**
-         * Initialize the boundary conditions, if necessary. This
-         * method does nothing, but may be overloaded in derived classes if
-         * initialization is needed.
-         */
-        //        virtual void initBC() {}
-
-        /**
-         * Evaluate the reactor boundary conditions. This procedure is
-         * called during integration to evaluate the rate of volume
-         * change \f$ dV/dt \f$ [m^3/s], the heat loss rate [W], and
-         * the species production rates due to surface chemistry.
-         *
-         * It may be overloaded in derived classes to implement other
-         * boundary conditions. If not overloaded, this routine 
-         * implements the following boundary conditions.
-         *
-         * The rate of volume change is
-         * \f[
-         * dV/dt = K ( P - P_{ext})
-         * \f]
-         * where K is set in procedure setVDotCoeff.
-         *
-         * 
-         * The heat loss rate is calculated from
-         * \f[
-         * \dot Q_{out} = h A (T - T_0) + \epsilon A (T^4 - T_{0,R}^4).
-         * \f]
-         * @see setArea, setEmissivity, setExtRadTemp
-         */
-
-//         virtual void evalBC(doublereal& vdot, 
-//             doublereal& heatLossRate, doublereal* sdot) {
-//             doublereal t = m_mix->temperature();
-
-//             //m_p0 = m_env->pressure();
-//             vdot = m_kv * (m_thermo->pressure()/m_p0 - 1.0)*m_vol0;
-//             heatLossRate = m_area * (
-//                 m_h * (t - m_ext_temp)
-//                 + m_emis * StefanBoltz * (t*t*t*t - m_ext_temp4)
-//                 );
-//         }
-
-        //@}
-
-        //-----------------------------------------------------
-
         /**
          * Set the mixture to a state consistent with solution
          * vector y.
@@ -273,21 +185,15 @@ namespace Cantera {
     protected:
         
         Kinetics*   m_kin;
-        //        ReactorBase*     m_env;
-        //        Thermo*     m_thermo;
 
         Integrator* m_integ;         // pointer to integrator
         doublereal m_temp_atol;      // tolerance on T
         doublereal m_maxstep;        // max step size
         doublereal m_vdot, m_Q;
-        //        doublereal m_emis, m_h, m_area;
-        //doublereal m_ext_temp, m_ext_temp4;
-        //doublereal m_kv, m_p0;
         vector_fp m_atol;
         doublereal m_rtol;
         vector_fp m_work;
         vector_fp m_sdot;            // surface production rates
-        //bool m_trad_set;
         bool m_chem;
         bool m_energy;
         int m_nv;

Cantera/src/zeroD/ReactorBase.cpp

@@ -32,6 +32,7 @@ namespace Cantera {
                                  m_enthalpy(0.0),
                                  m_intEnergy(0.0), 
                                  m_pressure(0.0),
+                                 m_kappa(0.0),
                                  m_nwalls(0)
     {
         m_name = name;
@@ -50,9 +51,11 @@ namespace Cantera {
         m_thermo = &thermo;
         m_nsp = m_mix->nSpecies();
         m_mix->saveState(m_state);
+        m_rho0 = m_thermo->density();
         m_enthalpy = m_thermo->enthalpy_mass();
         m_intEnergy = m_thermo->intEnergy_mass();
         m_pressure = m_thermo->pressure();
+        m_kappa = m_thermo->isothermalCompressibility();
     }
 
     void ReactorBase::addInlet(FlowDevice& inlet) {

Cantera/src/zeroD/ReactorBase.h

@@ -134,6 +134,8 @@ namespace Cantera {
         doublereal mass() const { return m_vol * density(); }
         const doublereal* massFractions() const { return m_state.begin() + 2; }
         doublereal massFraction(int k) const { return m_state[k+2]; }
+        doublereal compressibility() const { return m_kappa; }
+
         //@}
 
         int error(string msg) const {
@@ -142,7 +144,7 @@ namespace Cantera {
         }
 
     protected:
-
+ 
         int m_nsp;
         thermo_t* m_mix;
         thermo_t*  m_thermo;
@@ -160,7 +162,9 @@ namespace Cantera {
         vector_int m_lr;
         int m_nwalls;
         string m_name;
-        
+        double m_rho0;
+        double m_kappa;
+
     private:
 
         void tilt(string method="") const { 

Cantera/src/zeroD/ReactorNet.cpp

@@ -4,7 +4,8 @@
 namespace Cantera {
 
     ReactorNet::ReactorNet() : FuncEval(), m_nr(0), m_nreactors(0),
-                               m_integ(0), m_init(false), m_nv(0), m_rtol(1.0e-6),
+                               m_integ(0), m_init(false), 
+                               m_nv(0), m_rtol(1.0e-6),
                                m_verbose(false)
     {
         m_integ = new CVodeInt;
@@ -23,6 +24,9 @@ namespace Cantera {
         m_nv = 0;
         m_reactors.clear();
         m_nreactors = 0;
+        if (m_verbose) {
+            writelog("Initializing reactor network.\n");
+        }
         for (n = 0; n < m_nr; n++) {
             if (m_r[n]->type() == ReactorType) {
                 m_r[n]->initialize(t0);
@@ -54,7 +58,7 @@ namespace Cantera {
 
     void ReactorNet::advance(doublereal time) {
         if (!m_init) {
-            m_maxstep = time;
+            m_maxstep = time - m_time;
             initialize();
         }
         m_integ->integrate(time);
@@ -64,7 +68,7 @@ namespace Cantera {
 
     double ReactorNet::step(doublereal time) {
         if (!m_init) {
-            m_maxstep = time;
+            m_maxstep = time - m_time;
             initialize();
         }
         m_time = m_integ->step(time);
@@ -75,10 +79,15 @@ namespace Cantera {
     void ReactorNet::eval(doublereal t, doublereal* y, doublereal* ydot) {
         int n;
         int start = 0;
-        updateState(y);
+        try {
-        for (n = 0; n < m_nreactors; n++) {
+            updateState(y);
-            m_reactors[n]->evalEqs(t, y + start, ydot + start);
+            for (n = 0; n < m_nreactors; n++) {
-            start += m_size[n];
+                m_reactors[n]->evalEqs(t, y + start, ydot + start);
+                start += m_size[n];
+            }
+        }
+        catch (CanteraError) {
+            showErrors(cout);
         }
     }
 

Cantera/src/zeroD/Wall.cpp

@@ -9,7 +9,7 @@ namespace Cantera {
 
         
     Wall::Wall() : m_left(0), m_right(0),  
-                   m_area(0.0), m_k(0.0), m_rrth(0.0), 
+                   m_area(0.0), m_k(0.0), m_rrth(0.0), m_emiss(0.0), 
                    m_vf(0), m_qf(0) {
         for (int n = 0; n < 2; n++) {
             m_chem[n] = 0;
@@ -41,6 +41,8 @@ namespace Cantera {
             if (ileft >= 0) {
                 m_surf[0] = (SurfPhase*)&left->thermo(ileft);
                 m_nsp[0] = m_surf[0]->nSpecies();
+                m_leftcov.resize(m_nsp[0]);
+                m_surf[0]->getCoverages(m_leftcov.begin());
             }
         }
         if (right) {
@@ -48,6 +50,8 @@ namespace Cantera {
             if (iright >= 0) {
                 m_surf[1] = (SurfPhase*)&right->thermo(iright);
                 m_nsp[1] = m_surf[1]->nSpecies();
+                m_rightcov.resize(m_nsp[1]);
+                m_surf[1]->getCoverages(m_rightcov.begin());
             }
         }
         if (ileft < 0 || iright < 0) {
@@ -61,11 +65,14 @@ namespace Cantera {
      * The volume rate of change is given by 
      * \f[ \dot V = K A (P_{left} - P_{right}) + F(t) \f]
      * where \f$ F(t) \f$ is a specified function of time.
+     *
+     * This method is used by class Reactor to compute the
+     * rate of volume change of the reactor.
      */
     doublereal Wall::vdot(doublereal t) {
         double rate1 = m_k * m_area * 
                        (m_left->pressure() - m_right->pressure()); 
-        if (m_vf) rate1 += m_vf->eval(t);
+        if (m_vf) rate1 += m_area * m_vf->eval(t);
         return rate1;
     }
 
@@ -78,8 +85,33 @@ namespace Cantera {
     doublereal Wall::Q(doublereal t) {
         double q1 = (m_area * m_rrth) *
                     (m_left->temperature() - m_right->temperature());
+        if (m_emiss > 0.0) {
+            double tl = m_left->temperature();
+            double tr = m_right->temperature();
+            q1 += m_area * StefanBoltz * (tl*tl*tl*tl - tr*tr*tr*tr);
+        }
         if (m_qf) q1 += m_area * m_qf->eval(t);
         return q1;
     }
 
+    void Wall::setCoverages(int leftright, const doublereal* cov) {
+        if (leftright == 0)
+            copy(cov, cov + m_nsp[0], m_leftcov.begin());
+        else
+            copy(cov, cov + m_nsp[1], m_rightcov.begin());
+    }
+
+    void Wall::getCoverages(int leftright, doublereal* cov) {
+        if (leftright == 0)
+            copy(m_leftcov.begin(), m_leftcov.end(), cov);
+        else
+            copy(m_rightcov.begin(), m_rightcov.end(), cov);
+    }
+
+    void Wall::syncCoverages(int leftright) {
+        if (leftright == 0)
+            m_surf[0]->setCoverages(m_leftcov.begin());
+        else
+            m_surf[1]->setCoverages(m_rightcov.begin());
+    }
 }        

Cantera/src/zeroD/Wall.h

@@ -57,8 +57,13 @@ namespace Cantera {
         /// Set the overall heat transfer coefficient [W/m^2/K].
         void setHeatTransferCoeff(doublereal U) { m_rrth = U; }
 
-        /** Set the rate of volume change to a specified function.*/
+        void setEmissivity(doublereal epsilon) { m_emiss = epsilon; }
-        void setExpansionRate(Func1* f=0) {if (f) m_vf = f;}
+
+        //  /** Set the rate of volume change to a specified function.*/
+        //        void setExpansionRate(Func1* f=0) {if (f) m_vf = f;}
+
+        /** Set the piston velocity to a specified function. */
+        void setVelocity(Func1* f=0) {if (f) m_vf = f;}
 
         /** 
          * Set the expansion rate coefficient.
@@ -98,6 +103,13 @@ namespace Cantera {
             return m_chem[leftright];
         }
 
+        void setCoverages(int leftright, const doublereal* cov);
+
+        void getCoverages(int leftright, doublereal* cov);
+
+        void syncCoverages(int leftright);
+
+
     protected:
 
         vector_fp m_coeffs;
@@ -108,8 +120,10 @@ namespace Cantera {
         SurfPhase* m_surf[2];
         int m_nsp[2];
         doublereal m_area, m_k, m_rrth;
+        doublereal m_emiss;
         Func1 *m_vf;
         Func1 *m_qf;
+        vector_fp m_leftcov, m_rightcov;
 
     private:
 

Cantera/src/zeroD/flowControllers.h

@@ -21,141 +21,107 @@
 
 #include "FlowDevice.h"
 #include "ReactorBase.h"
-#include "PID_Controller.h"
+//#include "PID_Controller.h"
 #include "../Func1.h"
 
 namespace Cantera {
 
 
-    ///////////////////////////////////////////////////////////////
-
-
-    /** 
-     * A base class for devices that do not use closed-loop control.
-     * This is defined only for convenience, in order to overload
-     * virtual methods of FlowDevice that print warnings with ones
-     * that do nothing.
-     */
-    class NoController : public FlowDevice {
-    public:
-
-        NoController() {}
-        virtual ~NoController() {}
-
-        NoController(const NoController& a) 
-            : FlowDevice(a) {}
-
-        NoController& operator=(const NoController& a) {
-            return *this;
-        }
-
-        // unneeded methods 
-        virtual void update() {}
-        virtual void reset() {}
-        virtual bool setGains(int n, const doublereal* gains) {return true;}
-        virtual bool getGains(int n, doublereal* gains) {return true;}
-        virtual doublereal maxError() { return 0.0; }
-        virtual doublereal setpoint() { return 0.0; }
-        virtual void setSetpoint(doublereal mdot) { }
-        virtual bool ready() {
-            return FlowDevice::ready();
-        }
-
-    protected:
-    private:
-    };
-
-
-    //////////////////////////////////////////////////////////
-
-
     /**
      * A class for mass flow controllers. The mass flow rate is constant,
      * independent of any other parameters.
      */
-    class MassFlowController : public NoController {
+    class MassFlowController : public FlowDevice {
     public:
 
-        MassFlowController() {}
+        MassFlowController() : FlowDevice() {
-        virtual ~MassFlowController() {}
+            m_type = MFC_Type;
-
-        MassFlowController(const MassFlowController& a) 
-            : NoController(a) {}
-
-        MassFlowController& operator=(const MassFlowController& a) {
-            if (this == &a) return *this;
-            m_mdot = a.m_mdot;
-            return *this;
         }
 
-        virtual doublereal setpoint() { return m_mdot; }
+        virtual ~MassFlowController() {}
-        virtual void setSetpoint(doublereal mdot) { m_mdot = mdot; }
+
         virtual bool ready() {
             return FlowDevice::ready() && m_mdot >= 0.0;
         }
 
-    protected:
+        /// If a function of time has been specified for
-    private:
+        /// mdot, then update the stored mass flow rate.
-    };
+        /// Otherwise, mdot is a constant, and does not need
-
+        /// updating.
-
+        virtual void updateMassFlowRate(doublereal time) {
-    class UserValve : public NoController {
+            if (m_func) m_mdot = m_func->eval(time);
-    public:
+            if (m_mdot < 0.0) m_mdot = 0.0;
-
-        UserValve() : m_func(0) {}
-        virtual ~UserValve() {}
-
-        UserValve(const UserValve& a) : NoController(a) {}
-
-        UserValve& operator=(const UserValve& a) {
-            if (this == &a) return *this;
-            m_func = a.m_func;
-            return *this;
-        }
-
-        virtual bool ready() {
-            return FlowDevice::ready() && m_func != 0;
-        }
-
-        virtual void setFunction(Func1* f) { m_func = f; }
-
-        virtual doublereal massFlowRate() {
-            return m_func->eval(in().pressure() - out().pressure());
         }
 
     protected:
-        Func1* m_func;
 
     private:
     };
 
-
     /**
      * A class for mass flow controllers. The mass flow rate is constant,
      * independent of any other parameters.
      */
-    class Valve : public NoController {
+    class PressureController : public FlowDevice {
     public:
 
-        Valve() {}
+        PressureController() : FlowDevice(), m_master(0) {
-        virtual ~Valve() {}
+            m_type = PressureController_Type;
-
-        Valve(const Valve& a) : NoController(a) {}
-
-        Valve& operator=(const Valve& a) {
-            if (this == &a) return *this;
-            m_mdot = a.m_mdot;
-            return *this;
         }
 
+        virtual ~PressureController() {}
+
+        virtual bool ready() {
+            return FlowDevice::ready() && m_master != 0;
+        }
+
+        void setMaster(FlowDevice* master) {
+            m_master = master;
+        }
+
+        virtual void updateMassFlowRate(doublereal time) {
+            doublereal master_mdot = m_master->massFlowRate(time);
+            m_mdot = master_mdot + m_coeffs[0]*(in().pressure() - 
+                out().pressure());
+            if (m_mdot < 0.0) m_mdot = 0.0;
+        }
+
+    protected:
+        FlowDevice* m_master;
+
+    private:
+    };
+
+
+    /// Valve objects supply a mass flow rate that is a function of the
+    /// pressure drop across the valve. The default behavior is a linearly
+    /// proportional to the pressure difference. Note that
+    /// real valves do not have this behavior, so this class
+    /// does not model real, physical valves.
+    class Valve : public FlowDevice {
+    public:
+
+        Valve() : FlowDevice() {
+            m_type = Valve_Type;
+        }
+
+        virtual ~Valve() {}
+
         virtual bool ready() {
             return FlowDevice::ready() && m_coeffs.size() >= 1;
         }
 
-        virtual doublereal massFlowRate() {
+        /// Compute the currrent mass flow rate, based on
-            m_mdot = m_coeffs[0]* (in().pressure() - out().pressure());
+        /// the pressure difference.
+        virtual void updateMassFlowRate(doublereal time) {
+            double delta_P = in().pressure() - out().pressure();
+            if (m_func) {
+                m_mdot = m_func->eval(delta_P);
+            }
+            else {
+                m_mdot = m_coeffs[0]*delta_P;
+            }
             if (m_mdot < 0.0) m_mdot = 0.0;
-            return m_mdot;
         }
 
     protected:
@@ -163,57 +129,6 @@ namespace Cantera {
     private:
     };
 
-
-
-    /**
-     * A PressureRegulator is a device that controls the pressure 
-     * of the upstream reactor by regulating the mass flow rate.
-     */ 
-    class PressureRegulator : public FlowDevice  {
-
-    public:
-
-        PressureRegulator() {}
-        virtual ~PressureRegulator() {}
-        PressureRegulator(const PressureRegulator& p) : m_pid(p.m_pid) {}
-        PressureRegulator& operator=(const PressureRegulator& p) {
-            if (this == &p) return *this;
-            m_pid = p.m_pid;
-            return *this;
-        }
-
-        // overloaded virtual methods
-        virtual void setSetpoint(doublereal p0) { m_pid.setpoint(-p0); }
-        virtual doublereal setpoint() { return -m_pid.setpoint(); }
-        virtual bool ready() { 
-            return FlowDevice::ready() && m_pid.ready(); }
-        virtual void reset() {
-            m_pid.reset(in().time()-1.e-12, -in().pressure());
-        }
-        virtual void update() {
-            m_pid.update(in().time(), -in().pressure());
-        }
-
-        virtual bool setGains(int n, const doublereal* gains) {
-            return m_pid.setGains(n, gains);
-        }
-
-        virtual bool getGains(int n, doublereal* gains) {
-            return m_pid.getGains(n, gains);
-        }
-
-        virtual doublereal maxError() { return m_pid.maxError(); }
-
-        virtual doublereal massFlowRate() {
-            m_mdot =  m_pid.output(-in().pressure());
-            return m_mdot;
-        }
-
-    protected:
-
-    private:
-        PID_Controller m_pid;
-    };
 }
 #endif
 

configure

@@ -183,7 +183,7 @@ LAPACK_FTN_STRING_LEN_AT_END='y'
 CXX=${CXX:=g++}
 
 # C++ compiler flags
-CXXFLAGS=${CXXFLAGS:="-O0 -g -Wall"}
+CXXFLAGS=${CXXFLAGS:="-O0 -Wall"}
 
 # the C++ flags required for linking
 #LCXX_FLAGS=

data/inputs/ptcombust.cti

@@ -1,5 +1,6 @@
 #
-# see http://reaflow.iwr.uni-heidelberg.de/~Olaf.Deutschmann/ for more about this mechanism
+# see http://reaflow.iwr.uni-heidelberg.de/~Olaf.Deutschmann/ for 
+# more about this mechanism
 #
 #---------------------------------------------------------------------!
 #***********************************************************************
@@ -28,9 +29,17 @@
 
 units(length = "cm", time = "s", quantity = "mol", act_energy = "J/mol")
 
+
+#
+# Define a gas mixture with species imported from GRI-Mech.
+# Reactions will be imported from GRI-Mech 3.0, as long as they
+# don't involve species not declared here. Transport properties
+# will be computed using a mixture-averaged model.
+# 
 ideal_gas(name = "gas",
          elements = "O H C N Ar",
-         species = """gri30: H2      H       O       O2      OH      H2O     HO2     H2O2 
+         species = """gri30: H2      H       O       O2      OH      
+                             H2O     HO2     H2O2 
          C       CH      CH2     CH2(S)  CH3     CH4     CO      CO2     
          HCO     CH2O    CH2OH   CH3O    CH3OH   C2H     C2H2    C2H3    
          C2H4    C2H5    C2H6    HCCO    CH2CO   HCCOH AR N2""",
@@ -55,12 +64,12 @@ ideal_interface(name = "Pt_surf",
                                       coverages = 'O(S):0.0, PT(S):0.5, H(S):0.5')
                 )
 
-#-------------------------------------------------------------------------------
+#-----------------------------------------------------------------------------
 #  Species data
 #
 # Note that reactions 12-14 are reversible, and therefore require thermo
 # data
-#-------------------------------------------------------------------------------
+#-----------------------------------------------------------------------------
 
 species(name = "PT(S)",
     atoms = " Pt:1 ",
@@ -202,22 +211,23 @@ species(name = "O(S)",
 
 #  Reaction 1
 surface_reaction("H2 + 2 PT(S) => 2 H(S)", [4.45790E+10, 0.5, 0], 
-		  order = "PT(S):1")
+                 order = "PT(S):1")
 
 #  Reaction 2
 surface_reaction( "2 H(S) => H2 + 2 PT(S)",   
-		  Arrhenius(3.70000E+21, 0, 67400, coverage = ['H(S)', 0.0, 0.0, -6000.0]))
+		  Arrhenius(3.70000E+21, 0, 67400,
+                            coverage = ['H(S)', 0.0, 0.0, -6000.0]))
 
 #  Reaction 3
 surface_reaction( "H + PT(S) => H(S)",   stick(1.00000E+00, 0, 0))
 
 #  Reaction 4
 surface_reaction( "O2 + 2 PT(S) => 2 O(S)",   Arrhenius(1.80000E+21, -0.5, 0),
-         options = 'duplicate')
+                  options = 'duplicate')
 
 #  Reaction 5
 surface_reaction( "O2 + 2 PT(S) => 2 O(S)",   stick(2.30000E-02, 0, 0),
-         options = 'duplicate')
+                  options = 'duplicate')
 
 #  Reaction 6
 surface_reaction( "2 O(S) => O2 + 2 PT(S)",   
@@ -266,7 +276,8 @@ surface_reaction( "CH4 + 2 PT(S) => CH3(S) + H(S)",   [4.63340E+20, 0.5, 0],
 		  order = "PT(S):2.3")
 
 #  Reaction 20
-surface_reaction( "CH3(S) + PT(S) => CH2(S)s + H(S)",   [3.70000E+21, 0, 20000])
+surface_reaction( "CH3(S) + PT(S) => CH2(S)s + H(S)",
+                  [3.70000E+21, 0, 20000])
 
 #  Reaction 21
 surface_reaction( "CH2(S)s + PT(S) => CH(S) + H(S)",   [3.70000E+21, 0, 20000])

ext/tpx/Sub.h

@@ -133,7 +133,7 @@ namespace tpx {
             return T*(s2 - s1)/(2.0*dt);
         }
 
-        virtual double thermExpCoeff() {
+        virtual double thermalExpansionCoeff() {
             double Tsave = T, dt = 1.e-4*T;
             double p0 = P();
             Set(TP, Tsave - dt, p0);
@@ -144,6 +144,16 @@ namespace tpx {
             return (v2 - v1)/((v2 + v1)*dt);
         }
 
+        virtual double isothermalCompressibility() {
+            double Psave = P(), dp = 1.e-4*Psave;
+            Set(TP, T, Psave - dp);
+            double v1 = v();
+            Set(TP, T, Psave + dp);
+            double v2 = v();
+            Set(TP, T, Psave);
+            return -(v2 - v1)/((v2 + v1)*dp);
+        }
+
 
         // saturation properties