[Cython] Full implementation of zero-dimensional reactor networks

include/cantera/zeroD/Reactor.h

@@ -97,6 +97,11 @@ public:
         }
     }
 
+    /// Returns 'true' if solution of the energy equation is enabled.
+    bool energyEnabled() const {
+        return m_energy;
+    }
+
     // overloaded methods of class FuncEval
     virtual size_t neq() {
         return m_nv;

include/cantera/zeroD/ReactorBase.h

@@ -119,7 +119,14 @@ public:
 
     //@}
 
-    void resetState();
+    //! Set the state of the Phase object associated with this reactor to the
+    //! reactor's current state.
+    void restoreState() {
+        if (!m_thermo) {
+            throw CanteraError("ReactorBase::restoreState", "No phase defined.");
+        }
+        m_thermo->restoreState(m_state);
+    }
 
     /// return a reference to the contents.
     thermo_t& contents() {

include/cantera/zeroD/Wall.h

@@ -51,6 +51,11 @@ public:
         m_area = a;
     }
 
+    //! Get the area [m^2]
+    double getArea() const {
+        return m_area;
+    }
+
     void setThermalResistance(doublereal Rth) {
         m_rrth = 1.0/Rth;
     }
@@ -60,6 +65,11 @@ public:
         m_rrth = U;
     }
 
+    //! Get the overall heat transfer coefficient [W/m^2/K].
+    double getHeatTransferCoeff() const {
+        return m_rrth;
+    }
+
     /// Set the emissivity.
     void setEmissivity(doublereal epsilon) {
         if (epsilon > 1.0 || epsilon < 0.0)
@@ -68,6 +78,10 @@ public:
         m_emiss = epsilon;
     }
 
+    double getEmissivity() const {
+        return m_emiss;
+    }
+
     /** Set the piston velocity to a specified function. */
     void setVelocity(Cantera::Func1* f=0) {
         if (f) {
@@ -82,6 +96,10 @@ public:
         m_k = k;
     }
 
+    //! Get the expansion rate coefficient
+    double getExpansionRateCoeff() const {
+        return m_k;
+    }
 
     /// Specify the heat flux function \f$ q_0(t) \f$.
     void setHeatFlux(Cantera::Func1* q) {

interfaces/cython/cantera/_cantera.pxd

@@ -206,32 +206,90 @@ cdef extern from "cantera/equil/vcs_MultiPhaseEquil.h" namespace "Cantera":
 
 cdef extern from "cantera/zeroD/ReactorBase.h" namespace "Cantera":
     cdef cppclass CxxWall "Cantera::Wall"
+    cdef cppclass CxxFlowDevice "Cantera::FlowDevice"
 
     cdef cppclass CxxReactorBase "Cantera::ReactorBase":
         CxxReactorBase()
-        void setThermoMgr(CxxThermoPhase&)
+        void setThermoMgr(CxxThermoPhase&) except +
+        void restoreState() except +
         double volume()
+        string name()
+        void setName(string)
+        void setInitialVolume(double)
+
 
 cdef extern from "cantera/zeroD/Reactor.h":
-    cdef cppclass CxxReactor "Cantera::Reactor":
+    cdef cppclass CxxReactor "Cantera::Reactor" (CxxReactorBase):
         CxxReactor()
         void setKineticsMgr(CxxKinetics&)
+        void setEnergy(int)
+        cbool energyEnabled()
+
+
+cdef extern from "cantera/zeroD/FlowReactor.h":
+    cdef cppclass CxxFlowReactor "Cantera::FlowReactor" (CxxReactor):
+        CxxFlowReactor()
+        void setMassFlowRate(double) except +
+        double speed()
+        double distance()
+
 
 cdef extern from "cantera/zeroD/Wall.h":
     cdef cppclass CxxWall "Cantera::Wall":
         CxxWall()
         cbool install(CxxReactorBase&, CxxReactorBase&)
         void setExpansionRateCoeff(double)
+        double getExpansionRateCoeff()
         double area()
         void setArea(double)
+        double getArea()
+        void setHeatTransferCoeff(double)
+        double getHeatTransferCoeff()
+        void setEmissivity(double) except +
+        double getEmissivity()
+        void setVelocity(CxxFunc1*)
+        void setHeatFlux(CxxFunc1*)
+        void setKinetics(CxxKinetics*, CxxKinetics*)
+        void setCoverages(int, double*)
+        void syncCoverages(int)
+        double vdot(double)
+        double Q(double)
+
+
+cdef extern from "cantera/zeroD/flowControllers.h":
+    cdef cppclass CxxFlowDevice "Cantera::FlowDevice":
+        CxxFlowDevice()
+        double massFlowRate(double)
+        cbool install(CxxReactorBase&, CxxReactorBase&)
+        void setFunction(CxxFunc1*)
+        void setParameters(int, double*)
+
+    cdef cppclass CxxMassFlowController "Cantera::MassFlowController" (CxxFlowDevice):
+        CxxMassFlowController()
+
+    cdef cppclass CxxValve "Cantera::Valve" (CxxFlowDevice):
+        CxxValve()
+
+    cdef cppclass CxxPressureController "Cantera::PressureController" (CxxFlowDevice):
+        CxxPressureController()
+        void setMaster(CxxFlowDevice*)
+
 
 cdef extern from "cantera/zeroD/ReactorNet.h":
     cdef cppclass CxxReactorNet "Cantera::ReactorNet":
         CxxReactorNet()
         void addReactor(CxxReactorBase*)
-        void initialize(double)
-        void advance(double)
-        double step(double)
+        void advance(double) except +
+        double step(double) except +
+        double time()
+        void setInitialTime(double)
+        void setTolerances(double, double)
+        double rtol()
+        double atol()
+        void setMaxTimeStep(double)
+        cbool verbose()
+        void setVerbose(cbool)
+
 
 cdef extern from "cantera/thermo/ThermoFactory.h" namespace "Cantera":
     cdef CxxThermoPhase* newPhase(string, string) except +
@@ -356,6 +414,9 @@ cdef class Mixture:
     cdef CxxMultiPhase* mix
     cdef list _phases
 
+cdef class Kinetics(_SolutionBase):
+    pass
+
 cdef class Func1:
     cdef CxxFunc1* func
     cdef object callable
@@ -363,7 +424,24 @@ cdef class Func1:
 
 cdef class ReactorBase:
     cdef CxxReactorBase* rbase
+    cdef object _thermo
+    cdef list _inlets
+    cdef list _outlets
+    cdef list _walls
+
+cdef class Reactor(ReactorBase):
+    cdef CxxReactor* reactor
+    cdef object _kinetics
 
 cdef class Wall:
     cdef CxxWall* wall
-    cdef double _expansionRateCoeff
+    cdef object _velocityFunc
+    cdef object _heatFluxFunc
+    cdef str name
+    cdef Kinetics _leftKinetics
+    cdef Kinetics _rightKinetics
+
+cdef class FlowDevice:
+    cdef CxxFlowDevice* dev
+    cdef Func1 _rateFunc
+    cdef str name

interfaces/cython/cantera/reactor.pyx

@@ -1,49 +1,202 @@
+from collections import defaultdict
+import numbers
+
+reactor_counts = defaultdict(int)
+
 cdef class ReactorBase:
     reactorType = "None"
     def __cinit__(self, *args, **kwargs):
         self.rbase = newReactor(stringify(self.reactorType))
 
-    def __init__(self, *args, **kwargs):
-        if args and isinstance(args[0], _SolutionBase):
-            self.insert(args[0])
+    def __init__(self, ThermoPhase contents=None, name=None, **kwargs):
+        self._inlets = []
+        self._outlets = []
+        self._walls = []
+        if isinstance(contents, ThermoPhase):
+            self.insert(contents)
+
+        if name is not None:
+            self.name = name
+        else:
+            reactor_counts[self.reactorType] += 1
+            n = reactor_counts[self.reactorType]
+            self.name = '{0}_{1}'.format(self.reactorType, n)
 
     def __dealloc__(self):
         del self.rbase
 
     def insert(self, _SolutionBase solution):
+        self._thermo = solution
         self.rbase.setThermoMgr(deref(solution.thermo))
 
-    def volume(self):
-        return self.rbase.volume()
+    property name:
+        def __get__(self):
+            return pystr(self.rbase.name())
+
+        def __set__(self, name):
+            self.rbase.setName(stringify(name))
+
+    property thermo:
+        def __get__(self):
+            self.rbase.restoreState()
+            return self._thermo
+
+    property volume:
+        def __get__(self):
+            return self.rbase.volume()
+
+        def __set__(self, double value):
+            self.rbase.setInitialVolume(value)
+
+    property T:
+        def __get__(self):
+            return self.thermo.T
+
+    property density:
+        def __get__(self):
+            return self.thermo.density
+
+    property Y:
+        def __get__(self):
+            return self.thermo.Y
+
+    # Flow devices & walls
+    property inlets:
+        """List of flow devices installed as inlets to this reactor"""
+        def __get__(self):
+            return self._inlets
+
+    property outlets:
+        """List of flow devices installed as outlets to this reactor"""
+        def __get__(self):
+            return self._outlets
+
+    property walls:
+        """List of walls installed on this reactor"""
+        def __get__(self):
+            return self._walls
+
+    def _addInlet(self, inlet):
+        """
+        Store a reference to *inlet* to prevent it from being prematurely
+        garbage collected.
+        """
+        self._inlets.append(inlet)
+
+    def _addOutlet(self, outlet):
+        """
+        Store a reference to *outlet* to prevent it from being prematurely
+        garbage collected.
+        """
+        self._outlets.append(outlet)
+
+    def _addWall(self, wall):
+        """
+        Store a reference to *wall* to prevent it from being prematurely
+        garbage collected.
+        """
+        self._walls.append(wall)
 
 
 cdef class Reactor(ReactorBase):
     reactorType = "Reactor"
-    cdef CxxReactor* reactor
 
     def __cinit__(self, *args, **kwargs):
         self.reactor = <CxxReactor*>(self.rbase)
 
+    def __init__(self, contents=None, *, name=None, energy='on', **kwargs):
+        super().__init__(contents, **kwargs)
+
+        if energy == 'off':
+            self.energyEnabled = False
+        elif energy != 'on':
+            raise ValueError("'energy' must be either 'on' or 'off'")
+
     def insert(self, _SolutionBase solution):
         ReactorBase.insert(self, solution)
+        self._kinetics = solution
         self.reactor.setKineticsMgr(deref(solution.kinetics))
 
+    property kinetics:
+        def __get__(self):
+            self.rbase.restoreState()
+            return self._kinetics
+
+    property energyEnabled:
+        def __get__(self):
+            return self.reactor.energyEnabled()
+
+        def __set__(self, pybool value):
+            self.reactor.setEnergy(int(value))
+
+
+cdef class Reservoir(ReactorBase):
+    reactorType = "Reservoir"
+
+
+cdef class ConstPressureReactor(Reactor):
+    reactorType = "ConstPressureReactor"
+
+
+cdef class FlowReactor(Reactor):
+    reactorType = "FlowReactor"
+
+    property massFlowRate:
+        def __set__(self, double value):
+            (<CxxFlowReactor*>self.reactor).setMassFlowRate(value)
+
+    property speed:
+        def __get__(self):
+            return (<CxxFlowReactor*>self.reactor).speed()
+
+    property distance:
+        def __get__(self):
+            return (<CxxFlowReactor*>self.reactor).distance()
+
 
 cdef class Wall:
     def __cinit__(self, *args, **kwargs):
         self.wall = new CxxWall()
 
-    def __init__(self, *args, **kwargs):
-        self._expansionRateCoeff = 0.0
+    def __init__(self, left, right, *, name=None, A=None, K=None, U=None,
+                 Q=None, velocity=None, kinetics=(None,None)):
+        self._velocityFunc = None
+        self._heatFluxFunc = None
+        self._leftKinetics = None
+        self._rightKinetics = None
 
-    def install(self, ReactorBase left, ReactorBase right):
+        self._install(left, right)
+        if name is not None:
+            self.name = name
+        else:
+            reactor_counts['Wall'] += 1
+            n = reactor_counts['Wall']
+            self.name = 'Wall_{0}'.format(n)
+
+        if A is not None:
+            self.area = A
+        if K is not None:
+            self.expansionRateCoeff = K
+        if U is not None:
+            self.heatTransferCoeff = U
+        if Q is not None:
+            self.setHeatFlux(Q)
+        if velocity is not None:
+            self.setVelocity(velocity)
+        if kinetics[0] is not None:
+            self.leftKinetics = kinetics[0]
+        if kinetics[1] is not None:
+            self.leftKinetics = kinetics[1]
+
+    def _install(self, ReactorBase left, ReactorBase right):
+        left._addWall(self)
+        right._addWall(self)
         self.wall.install(deref(left.rbase), deref(right.rbase))
 
     property expansionRateCoeff:
         def __get__(self):
-            return self._expansionRateCoeff
+            return self.wall.getExpansionRateCoeff()
         def __set__(self, double val):
-            self._expansionRateCoeff = val
             self.wall.setExpansionRateCoeff(val)
 
     property area:
@@ -52,19 +205,235 @@ cdef class Wall:
         def __set__(self, double value):
             self.wall.setArea(value)
 
+    property heatTransferCoeff:
+        def __get__(self):
+            return self.wall.getHeatTransferCoeff()
+        def __set__(self, double value):
+            self.wall.setHeatTransferCoeff(value)
+
+    property emissivity:
+        def __get__(self):
+            return self.wall.getEmissivity()
+        def __set__(self, double value):
+            self.wall.setEmissivity(value)
+
+    def setVelocity(self, v):
+        cdef Func1 f
+        if isinstance(v, Func1):
+            f = v
+        else:
+            f = Func1(v)
+
+        self._velocityFunc = f
+        self.wall.setVelocity(f.func)
+
+    def setHeatFlux(self, q):
+        cdef Func1 f
+        if isinstance(q, Func1):
+            f = q
+        else:
+            f = Func1(q)
+
+        self._heatFluxFunc = f
+        self.wall.setHeatFlux(f.func)
+
+    def vdot(self, double t):
+        return self.wall.vdot(t)
+
+    def qdot(self, double t):
+        return self.wall.Q(t)
+
+    property leftKinetics:
+        def __get__(self):
+            return self._leftKinetics
+        def __set__(self, Kinetics k):
+            self._leftKinetics = k
+            self._setKinetics()
+
+    property rightKinetics:
+        def __get__(self):
+            return self._rightKinetics
+        def __set__(self, Kinetics k):
+            self._rightKinetics = k
+            self._setKinetics()
+
+    def _setKinetics(self):
+        cdef CxxKinetics* L = (self._leftKinetics.kinetics
+                               if self._leftKinetics else NULL)
+        cdef CxxKinetics* R = (self._rightKinetics.kinetics
+                               if self._rightKinetics else NULL)
+        self.wall.setKinetics(L, R)
+
+    property leftCoverages:
+        def __get__(self):
+            if self._leftKinetics is None:
+                raise Exception('No kinetics manager present')
+            self.wall.syncCoverages(0)
+            return self.leftKinetics.coverages
+        def __set__(self, coverages):
+            if self._leftKinetics is None:
+                raise Exception("Can't set coverages before assigning kinetics manager.")
+            if len(coverages) != self._leftKinetics.nSpecies:
+                raise ValueError('Incorrect number of site coverages specified')
+            cdef np.ndarray[np.double_t, ndim=1] data = \
+                np.ascontiguousarray(coverages, dtype=np.double)
+            self.wall.setCoverages(0, &data[0])
+
+    property rightCoverages:
+        def __get__(self):
+            if self._rightKinetics is None:
+                raise Exception('No kinetics manager present')
+            self.wall.syncCoverages(1)
+            return self._rightKinetics.coverages
+        def __set__(self, coverages):
+            if self._rightKinetics is None:
+                raise Exception("Can't set coverages before assigning kinetics manager.")
+            if len(coverages) != self._rightKinetics.nSpecies:
+                raise ValueError('Incorrect number of site coverages specified')
+            cdef np.ndarray[np.double_t, ndim=1] data = \
+                np.ascontiguousarray(coverages, dtype=np.double)
+            self.wall.setCoverages(1, &data[0])
+
+
+cdef class FlowDevice:
+    def __cinit__(self, *args, **kwargs):
+        # Children of this abstract class are responsible for allocating dev
+        self.dev = NULL
+
+    def __init__(self, upstream, downstream, *, name=None):
+        assert self.dev != NULL
+        self._rateFunc = None
+
+        if name is not None:
+            self.name = name
+        else:
+            reactor_counts[self.__class__.__name__] += 1
+            n = reactor_counts[self.__class__.__name__]
+            self.name = '{0}_{1}'.format(self.__class__.__name__, n)
+
+        self._install(upstream, downstream)
+
+    def __dealloc__(self):
+        del self.dev
+
+    def _install(self, ReactorBase upstream, ReactorBase downstream):
+        upstream._addOutlet(self)
+        downstream._addInlet(self)
+        self.dev.install(deref(upstream.rbase), deref(downstream.rbase))
+
+    def mdot(self, double t):
+        return self.dev.massFlowRate(t)
+
+
+cdef class MassFlowController(FlowDevice):
+    def __cinit__(self, *args, **kwargs):
+        self.dev = new CxxMassFlowController()
+
+    def __init__(self, upstream, downstream, *, name=None, mdot=None):
+        super().__init__(upstream, downstream, name=name)
+        if mdot is not None:
+            self.setMassFlowRate(mdot)
+
+    def setMassFlowRate(self, m):
+        cdef Func1 f
+        if isinstance(m, Func1):
+            f = m
+        else:
+            f = Func1(m)
+
+        self._rateFunc = f
+        self.dev.setFunction(f.func)
+
+
+cdef class Valve(FlowDevice):
+    def __cinit__(self, *args, **kwargs):
+        self.dev = new CxxValve()
+
+    def __init__(self, upstream, downstream, *, name=None, K=None):
+        super().__init__(upstream, downstream, name=name)
+        if K is not None:
+            self.setValveCoeff(K)
+
+    def setValveCoeff(self, k):
+        cdef double kv
+        cdef Func1 f
+        if isinstance(k, numbers.Real):
+            kv = k
+            self.dev.setParameters(1, &kv)
+            return
+
+        if isinstance(k, Func1):
+            f = k
+        else:
+            f = Func1(k)
+        self._rateFunc = f
+        self.dev.setFunction(f.func)
+
+
+cdef class PressureController(FlowDevice):
+    def __cinit__(self, *args, **kwargs):
+        self.dev = new CxxPressureController()
+
+    def __init__(self, upstream, downstream, *, name=None, master=None, K=None):
+        super().__init__(upstream, downstream, name=name)
+        if master is not None:
+            self.setMaster(master)
+        if K is not None:
+            self.setPressureCoeff(K)
+
+    def setPressureCoeff(self, double k):
+        self.dev.setParameters(1, &k)
+
+    def setMaster(self, FlowDevice d):
+        (<CxxPressureController*>self.dev).setMaster(d.dev)
+
+
 cdef class ReactorNet:
     cdef CxxReactorNet* net
+    cdef list _reactors
+
     def __cinit__(self, *args, **kwargs):
         self.net = new CxxReactorNet()
 
-    def addReactor(self, ReactorBase r):
-        self.net.addReactor(r.rbase)
+    def __init__(self, reactors=()):
+        self._reactors = []  # prevents premature garbage collection
+        for R in reactors:
+            self.addReactor(R)
 
-    def initialize(self, double t=0.0):
-        self.net.initialize(t)
+    def addReactor(self, ReactorBase r):
+        self._reactors.append(r)
+        self.net.addReactor(r.rbase)
 
     def advance(self, double t):
         self.net.advance(t)
 
     def step(self, double t):
         return self.net.step(t)
+
+    property time:
+        def __get__(self):
+            return self.net.time()
+
+    def setInitialTime(self, double t):
+        self.net.setInitialTime(t)
+
+    def setMaxTimeStep(self, double t):
+        self.net.setMaxTimeStep(t)
+
+    property rtol:
+        def __get__(self):
+            return self.net.rtol()
+        def __set__(self, tol):
+            self.net.setTolerances(tol, -1)
+
+    property atol:
+        def __get__(self):
+            return self.net.atol()
+        def __set__(self, tol):
+            self.net.setTolerances(-1, tol)
+
+    property verbose:
+        def __get__(self):
+            return pybool(self.verbose())
+        def __set__(self, pybool v):
+            self.net.setVerbose(v)

src/zeroD/ReactorBase.cpp

@@ -29,14 +29,6 @@ ReactorBase::ReactorBase(string name) : m_nsp(0),
     m_name = name;
 }
 
-//     void ReactorBase::resetState() {
-//         m_thermo->saveState(m_state);
-//         m_enthalpy = m_thermo->enthalpy_mass();
-//         m_intEnergy = m_thermo->intEnergy_mass();
-//         m_pressure = m_thermo->pressure();
-//         m_init = false;
-//     }
-
 void ReactorBase::setThermoMgr(thermo_t& thermo)
 {
     m_thermo = &thermo;