From 183c0b3409ddde93ce4f5b5fa1ccf84412fcc1f1 Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Wed, 3 Mar 2004 06:31:14 +0000 Subject: [PATCH] *** empty log message *** --- Cantera/clib/src/ctreactor.cpp | 123 ++++++++++++++++------- Cantera/clib/src/ctreactor.h | 9 ++ Cantera/matlab/cantera/polynom.m | 11 +- Cantera/python/Cantera/Reactor.py | 64 +++--------- Cantera/python/Cantera/ReactorNet.py | 68 +++++++++++++ Cantera/python/src/ctreactor_methods.cpp | 62 ++++++++++++ Cantera/python/src/methods.h | 6 ++ Cantera/src/zeroD/Makefile.in | 2 +- Cantera/src/zeroD/Reactor.cpp | 20 ++-- Cantera/src/zeroD/Reactor.h | 94 ++++++++--------- Cantera/src/zeroD/ReactorNet.cpp | 90 +++++++++++++++++ Cantera/src/zeroD/ReactorNet.h | 111 ++++++++++++++++++++ 12 files changed, 513 insertions(+), 147 deletions(-) create mode 100644 Cantera/python/Cantera/ReactorNet.py create mode 100644 Cantera/src/zeroD/ReactorNet.cpp create mode 100644 Cantera/src/zeroD/ReactorNet.h diff --git a/Cantera/clib/src/ctreactor.cpp b/Cantera/clib/src/ctreactor.cpp index 1f379f9c6..c4f76a076 100755 --- a/Cantera/clib/src/ctreactor.cpp +++ b/Cantera/clib/src/ctreactor.cpp @@ -1,6 +1,7 @@ // Cantera includes #include "zeroD/Reactor.h" +#include "zeroD/ReactorNet.h" #include "zeroD/Reservoir.h" #include "zeroD/Wall.h" #include "zeroD/flowControllers.h" @@ -22,10 +23,12 @@ #define DERR -999.999 typedef ReactorBase reactor_t; -typedef FlowDevice flowdev_t; -typedef Wall wall_t; +typedef ReactorNet reactornet_t; +typedef FlowDevice flowdev_t; +typedef Wall wall_t; Cabinet* Cabinet::__storage = 0; +Cabinet* Cabinet::__storage = 0; Cabinet* Cabinet::__storage = 0; Cabinet* Cabinet::__storage = 0; @@ -33,6 +36,10 @@ inline reactor_t* _reactor(int i) { return Cabinet::cabinet()->item(i); } +inline reactornet_t* _reactornet(int i) { + return Cabinet::cabinet()->item(i); +} + inline flowdev_t* _flowdev(int i) { return Cabinet::cabinet()->item(i); } @@ -158,46 +165,94 @@ 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); +// 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 + + int DLL_EXPORT reactornet_new() { + ReactorNet* r = new ReactorNet(); + return Cabinet::cabinet()->add(r); + } + + int DLL_EXPORT reactornet_del(int i) { + try { + Cabinet::cabinet()->del(i); + return 0; + } + catch (...) { + return -1; + } + } + + int DLL_EXPORT reactornet_copy(int i) { + return Cabinet::cabinet()->newCopy(i); + } + + int DLL_EXPORT reactornet_assign(int i, int j) { + return Cabinet::cabinet()->assign(i,j); + } + + int DLL_EXPORT reactornet_setInitialTime(int i, double t) { + _reactornet(i)->setInitialTime(t); return 0; } - int DLL_EXPORT reactor_setExtTemp(int i, double t) { - reactor_t* r = _reactor(i); - if (r->type() == ReactorType) ((Reactor*)r)->setExtTemp(t); + int DLL_EXPORT reactornet_addreactor(int i, int n) { + _reactornet(i)->addReactor(_reactor(n)); 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 reactornet_advance(int i, double t) { + try { + _reactornet(i)->advance(t); + return 0; + } + catch (CanteraError) {return -1;} } - 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; + double DLL_EXPORT reactornet_step(int i, double t) { + return _reactornet(i)->step(t); } diff --git a/Cantera/clib/src/ctreactor.h b/Cantera/clib/src/ctreactor.h index 2a1124e86..92eb5aa1f 100755 --- a/Cantera/clib/src/ctreactor.h +++ b/Cantera/clib/src/ctreactor.h @@ -35,6 +35,15 @@ extern "C" { //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); + int DLL_IMPORT reactornet_assign(int i, int j); + int DLL_IMPORT reactornet_setInitialTime(int i, double t); + int DLL_IMPORT reactornet_addreactor(int i, int n); + int DLL_IMPORT reactornet_advance(int i, double t); + double DLL_IMPORT reactornet_step(int i, double t); + int DLL_IMPORT flowdev_new(int type); int DLL_IMPORT flowdev_del(int i); //int DLL_IMPORT flowdev_copy(int i); diff --git a/Cantera/matlab/cantera/polynom.m b/Cantera/matlab/cantera/polynom.m index 52a0b64f1..700748e41 100644 --- a/Cantera/matlab/cantera/polynom.m +++ b/Cantera/matlab/cantera/polynom.m @@ -1,12 +1,15 @@ function poly = polynom(coeffs) % POLY - create a polynomial Func instance -% +% +% Create a polynomial: +% polynom([-2 6 3]) 3x^2 + 6x - 2 +% [n m] = size(coeffs); if n == 1 - poly = Func('polynomial',m - 1,coeffs) + poly = Func('polynomial',m - 1,coeffs); elseif m == 1 - poly = Func('polynomial',n - 1,coeffs) + poly = Func('polynomial',n - 1,coeffs); else - error('wrong shape for coefficient array') + error('wrong shape for coefficient array'); end diff --git a/Cantera/python/Cantera/Reactor.py b/Cantera/python/Cantera/Reactor.py index 233226797..1b5257e2c 100644 --- a/Cantera/python/Cantera/Reactor.py +++ b/Cantera/python/Cantera/Reactor.py @@ -10,7 +10,6 @@ import types class ReactorBase: """Base class for reactors.""" - def __init__(self, contents = None, type = -1): """ Create a new ReactorBase instance. If 'contents' is specified, @@ -22,18 +21,14 @@ class ReactorBase: if contents: self.insert(contents) - def __del__(self): """Delete the reactor instance.""" _cantera.reactor_del(self.__reactor_id) - def reactor_id(self): - """ - The integer index used to access the kernel reactor object. - """ + """The integer index used to access the kernel reactor + object. For internal use. """ return self.__reactor_id - def insert(self, contents): """ @@ -41,12 +36,14 @@ class ReactorBase: thermodynamic properties and kinetic rates. """ self.contents = contents - self.setThermoMgr(contents) - self.setKineticsMgr(contents) + _cantera.reactor_setThermoMgr(self.__reactor_id, contents._phase_id) + _cantera.reactor_setKineticsMgr(self.__reactor_id, contents.ckin) + #self.setThermoMgr(contents) + #self.setKineticsMgr(contents) def setInitialTime(self, t0): """Set the initial time. Restarts integration from this time - using the current state as the initial condition. Default: 0.0 s""" + using the current state as the initial condition. weDefault: 0.0 s""" _cantera.reactor_setInitialTime(self.__reactor_id, t0) def setInitialVolume(self, t0): @@ -54,6 +51,8 @@ class ReactorBase: _cantera.reactor_setInitialVolume(self.__reactor_id, t0) def setEnergy(self, e): + """Turn the energy equation on or off. If off, the reactor + temperature is held constant.""" ie = 1 if e == 'off': ie = 0 @@ -68,7 +67,7 @@ class ReactorBase: return _cantera.reactor_density(self.__reactor_id) def volume(self): - """Reactor volume [m^3].""" + """Volume [m^3].""" return _cantera.reactor_volume(self.__reactor_id) def time(self): @@ -93,21 +92,14 @@ class ReactorBase: def advance(self, time): """Advance the state of the reactor in time from the current - time to time 'time'.""" + time to time 'time'. Note: this method is deprecated. See class ReactorNet.""" return _cantera.reactor_advance(self.__reactor_id, time) def step(self, time): """Advance the state of the reactor in time from the current - time to time 'time'.""" + time to time 'time'. Note: this method is deprecated. See class ReactorNet.""" return _cantera.reactor_step(self.__reactor_id, time) - def setThermoMgr(self, th): - _cantera.reactor_setThermoMgr(self.__reactor_id, th._phase_id) - - def setKineticsMgr(self, kin): - _cantera.reactor_setKineticsMgr(self.__reactor_id, kin.ckin) - #self.setThermoMgr(kin.thrm) - def massFraction(self, k): """Mass fraction of species k.""" if type(k) == types.StringType: @@ -199,12 +191,6 @@ class FlowDevice: """ return _cantera.flowdev_setpoint(self.__fdev_id) -## def reset(self): -## """ -## Reset the flow controller. Only necessary for pressure regulators. -## """ -## _cantera.flowdev_reset(self.__fdev_id) - def install(self, upstream, downstream): """ Install the device between the upstream and downstream @@ -212,27 +198,11 @@ class FlowDevice: """ _cantera.flowdev_install(self.__fdev_id, upstream.reactor_id(), downstream.reactor_id()) -## self.reset() - -## def setGains(self, gains): -## g = array(gains,'d') -## n = len(g) -## _cantera.flowdev_setGains(self.__fdev_id, n, g) - -## def getGains(self): -## n = 4 -## return _cantera.flowdev_getGains(self.__fdev_id, n) - def setParameters(self, c): params = array(c,'d') n = len(params) return _cantera.flowdev_setParameters(self.__fdev_id, n, params) -## def maxError(self): -## return _cantera.flowdev_maxError(self.__fdev_id) - -## def update(self): -## _cantera.flowdev_update(self.__fdev_id) class MassFlowController(FlowDevice): @@ -245,15 +215,6 @@ class MassFlowController(FlowDevice): self.setSetpoint(mdot) -## class PressureRegulator(FlowDevice): -## def __init__(self, upstream=None, downstream=None): -## FlowDevice.__init__(self,2) -## if upstream and downstream: -## self.install(upstream, downstream) - -## def setPressure(self, p): -## self.setSetpoint(p) - class Valve(FlowDevice): def __init__(self, upstream=None, downstream=None): @@ -345,6 +306,7 @@ class Wall: right.reactor_id()) def setKinetics(self, left, right): + """Specify surface reaction mechanisms for the left and right sides of the wall.""" ileft = 0 iright = 0 if left: diff --git a/Cantera/python/Cantera/ReactorNet.py b/Cantera/python/Cantera/ReactorNet.py new file mode 100644 index 000000000..5be02c805 --- /dev/null +++ b/Cantera/python/Cantera/ReactorNet.py @@ -0,0 +1,68 @@ +""" +Reactor networks. + +""" + +import _cantera + +class ReactorNet: + + """Networks of reactors. ReactorNet objects are used to + simultaneously advance the state of a set of coupled reactors. + + Example: + + r1 = Reactor(gas1) + r2 = Reactor(gas2) + <... install walls, inlets, outlets, etc...> + + reactor_network = ReactorNet([r1, r2]) + reactor_network.advance(time) + + """ + + + def __init__(self, reactorlist = None): + """ + Create a new ReactorNet instance. If a list of reactors is supplied, + these will be added to the network. + """ + self.__reactornet_id = _cantera.reactornet_new() + if reactorlist: + for r in reactorlist: + self.add(r) + + + def __del__(self): + """Delete the reactor network instance.""" + _cantera.reactornet_del(self.__reactornet_id) + + + def reactornet_id(self): + """ + The integer index used to access the kernel reactornet object. For internal use. + """ + return self.__reactornet_id + + + def add(self, reactor): + """ + Add a reactor to the network. + """ + _cantera.reactornet_addreactor(self.__reactornet_id, reactor.reactor_id()) + + + def setInitialTime(self, t0): + """Set the initial time. Restarts integration from this time + using the current state as the initial condition. Default: 0.0 s""" + _cantera.reactornet_setInitialTime(self.__reactornet_id, t0) + + def advance(self, time): + """Advance the state of the reactor network in time from the current + time to time 'time'.""" + return _cantera.reactornet_advance(self.__reactornet_id, time) + + def step(self, time): + """Take a single internal time step toward time 'time'. + The time after taking the step is returned.""" + return _cantera.reactornet_step(self.__reactornet_id, time) diff --git a/Cantera/python/src/ctreactor_methods.cpp b/Cantera/python/src/ctreactor_methods.cpp index 34c37a06d..d0b728746 100644 --- a/Cantera/python/src/ctreactor_methods.cpp +++ b/Cantera/python/src/ctreactor_methods.cpp @@ -450,3 +450,65 @@ py_wall_ready(PyObject *self, PyObject *args) } +static PyObject* +py_reactornet_new(PyObject *self, PyObject *args) +{ + int n = reactornet_new(); + return Py_BuildValue("i",n); +} + +static PyObject* +py_reactornet_del(PyObject *self, PyObject *args) +{ + int n; + if (!PyArg_ParseTuple(args, "i:reactornet_del", &n)) + return NULL; + int iok = reactornet_del(n); + if (iok < 0) return reportError(iok); + return Py_BuildValue("i",0); +} + +static PyObject* +py_reactornet_setInitialTime(PyObject *self, PyObject *args) +{ + int n; + double t; + if (!PyArg_ParseTuple(args, "id:reactornet_setInitialTime", &n, &t)) + return NULL; + int iok = reactornet_setInitialTime(n, t); + if (iok < 0) return reportError(iok); + return Py_BuildValue("i",0); +} + +static PyObject* +py_reactornet_addreactor(PyObject *self, PyObject *args) +{ + int n, m; + if (!PyArg_ParseTuple(args, "ii:reactornet_addreactor", &n, &m)) + return NULL; + int iok = reactornet_addreactor(n, m); + if (iok < 0) return reportError(iok); + return Py_BuildValue("i",0); +} + +static PyObject* +py_reactornet_advance(PyObject *self, PyObject *args) +{ + int n; + double t; + if (!PyArg_ParseTuple(args, "id:reactornet_advance", &n, &t)) + return NULL; + int iok = reactornet_advance(n, t); + if (iok < 0) return reportError(iok); + return Py_BuildValue("i",0); +} + +static PyObject* +py_reactornet_step(PyObject *self, PyObject *args) +{ + int n; + double t; + if (!PyArg_ParseTuple(args, "id:reactornet_step", &n, &t)) + return NULL; + return Py_BuildValue("d",reactornet_step(n, t)); +} diff --git a/Cantera/python/src/methods.h b/Cantera/python/src/methods.h index 8843915fa..80df15188 100644 --- a/Cantera/python/src/methods.h +++ b/Cantera/python/src/methods.h @@ -193,6 +193,7 @@ static PyMethodDef ct_methods[] = { {"flowdev_ready", py_flowdev_ready, 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}, @@ -206,16 +207,21 @@ static PyMethodDef ct_methods[] = { {"reactor_time", py_reactor_time, METH_VARARGS}, {"reactor_advance", py_reactor_advance, 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}, {"flowdev_setParameters", py_flowdev_setParameters, METH_VARARGS}, {"flowdev_setFunction", py_flowdev_setFunction, METH_VARARGS}, {"reactor_mass", py_reactor_mass, METH_VARARGS}, {"reactor_new", py_reactor_new, METH_VARARGS}, + {"reactornet_new", py_reactornet_new, METH_VARARGS}, {"reactor_enthalpy_mass", py_reactor_enthalpy_mass, METH_VARARGS}, {"reactor_pressure", py_reactor_pressure, METH_VARARGS}, {"reactor_setInitialVolume", py_reactor_setInitialVolume, METH_VARARGS}, {"reactor_density", py_reactor_density, METH_VARARGS}, {"reactor_setKineticsMgr", py_reactor_setKineticsMgr, METH_VARARGS}, {"reactor_del", py_reactor_del, METH_VARARGS}, + {"reactornet_del", py_reactornet_del, METH_VARARGS}, {"reactor_intEnergy_mass", py_reactor_intEnergy_mass, METH_VARARGS}, {"reactor_massFraction", py_reactor_massFraction, METH_VARARGS}, {"wall_install", py_wall_install, METH_VARARGS}, diff --git a/Cantera/src/zeroD/Makefile.in b/Cantera/src/zeroD/Makefile.in index e18f2ba8d..4c26e01de 100644 --- a/Cantera/src/zeroD/Makefile.in +++ b/Cantera/src/zeroD/Makefile.in @@ -16,7 +16,7 @@ OBJDIR = . CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT) # stirred reactors -OBJS = Reactor.o ReactorBase.o FlowDevice.o Wall.o +OBJS = Reactor.o ReactorBase.o FlowDevice.o Wall.o ReactorNet.o CXX_INCLUDES = -I.. ZEROD_LIB = @buildlib@/libzeroD.a diff --git a/Cantera/src/zeroD/Reactor.cpp b/Cantera/src/zeroD/Reactor.cpp index ae5e09324..7c17aaa70 100644 --- a/Cantera/src/zeroD/Reactor.cpp +++ b/Cantera/src/zeroD/Reactor.cpp @@ -31,19 +31,12 @@ namespace Cantera { m_maxstep(0.0), m_vdot(0.0), m_Q(0.0), - m_emis(0.0), - m_h(0.0), - m_area(1.0), - m_ext_temp(0.0), - m_ext_temp4(0.0), - m_kv(0.0), - m_p0(OneAtm), m_rtol(1.e-9), - m_trad_set(false), m_chem(true), m_energy(true) { m_integ = new CVodeInt; + // use backward differencing, with a full Jacobian computed // numerically, and use a Newton linear iterator m_integ->setMethod(BDF_Method); @@ -138,6 +131,7 @@ namespace Cantera { // The components of y are the total internal energy, // the total volume, and the mass of each species. + // Set the mass fractions and density of the mixture. doublereal u = y[0]; @@ -180,18 +174,24 @@ namespace Cantera { m_enthalpy = m_thermo->enthalpy_mass(); m_pressure = m_thermo->pressure(); m_intEnergy = m_thermo->intEnergy_mass(); + } + void Reactor::eval(doublereal time, doublereal* y, doublereal* ydot) + { + updateState(y); // synchronize the reactor state with y + evalEqs(time, y, ydot); + } /** * Called by the integrator to evaluate ydot given y at time 'time'. */ - void Reactor::eval(doublereal time, doublereal* y, doublereal* ydot) + void Reactor::evalEqs(doublereal time, doublereal* y, doublereal* ydot) { int i, k, nk; m_time = time; - updateState(y); // synchronize the reactor state with y + // updateState(y); // synchronize the reactor state with y m_vdot = 0.0; m_Q = 0.0; diff --git a/Cantera/src/zeroD/Reactor.h b/Cantera/src/zeroD/Reactor.h index d117805e5..fd46387ff 100644 --- a/Cantera/src/zeroD/Reactor.h +++ b/Cantera/src/zeroD/Reactor.h @@ -127,52 +127,52 @@ 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 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 \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; - } +// /** +// * 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 setHeatTransferCoeff(doublereal h) { +// m_h = h; +// } - void setVDotCoeff(doublereal k) { - m_kv = k; - } +// void setVDotCoeff(doublereal k) { +// m_kv = k; +// } - void setEmissivity(doublereal emis) { - m_emis = emis; - } +// void setEmissivity(doublereal emis) { +// m_emis = emis; +// } - void setExtPressure(doublereal p0) { - m_p0 = p0; - } +// void setExtPressure(doublereal p0) { +// m_p0 = p0; +// } void disableChemistry() { m_chem = false; } void enableChemistry() { m_chem = true; } @@ -205,7 +205,7 @@ namespace Cantera { //----------------------------------------------------- virtual void initialize(doublereal t0 = 0.0); - + void evalEqs(doublereal t, doublereal* y, doublereal* ydot); /** * @name Methods to specify simulation options. @@ -267,9 +267,9 @@ namespace Cantera { * vector y. */ - protected: - virtual void updateState(doublereal* y); + + protected: Kinetics* m_kin; // ReactorBase* m_env; @@ -279,14 +279,14 @@ namespace Cantera { 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; + // 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_trad_set; bool m_chem; bool m_energy; int m_nv; diff --git a/Cantera/src/zeroD/ReactorNet.cpp b/Cantera/src/zeroD/ReactorNet.cpp new file mode 100644 index 000000000..e77a8a58d --- /dev/null +++ b/Cantera/src/zeroD/ReactorNet.cpp @@ -0,0 +1,90 @@ +#include "ReactorNet.h" +#include "../CVode.h" + +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 = 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); + } + + void ReactorNet::initialize(doublereal t0) { + int n, nv; + m_nv = 0; + m_reactors.clear(); + for (n = 0; n < m_nr; n++) { + m_r[n]->initialize(t0); + if (m_r[n]->type() == ReactorType) { + Reactor* r = (Reactor*)m_r[n]; + m_reactors.push_back(r); + nv = m_reactors[n]->neq(); + m_size.push_back(nv); + m_nv += nv; + m_nreactors++; + } + } + 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); + m_init = true; + } + + void ReactorNet::advance(doublereal time) { + if (!m_init) { + m_maxstep = time; + initialize(); + } + m_integ->integrate(time); + m_time = time; + updateState(m_integ->solution()); + } + + double ReactorNet::step(doublereal time) { + if (!m_init) { + m_maxstep = time; + initialize(); + } + m_time = m_integ->step(time); + updateState(m_integ->solution()); + return m_time; + } + + void ReactorNet::eval(doublereal t, doublereal* y, doublereal* ydot) { + int n; + int start = 0; + updateState(y); + for (n = 0; n < m_nreactors; n++) { + m_reactors[n]->evalEqs(t, y + start, ydot + start); + start += m_size[n]; + } + } + + void ReactorNet::updateState(doublereal* y) { + int n; + int start = 0; + for (n = 0; n < m_nreactors; n++) { + m_reactors[n]->updateState(y + start); + start += m_size[n]; + } + } + + void ReactorNet::getInitialConditions(doublereal t0, + size_t leny, doublereal* y) { + int n; + int start = 0; + for (n = 0; n < m_nreactors; n++) { + m_reactors[n]->getInitialConditions(t0, m_size[n], y + start); + start += m_size[n]; + } + } +} + diff --git a/Cantera/src/zeroD/ReactorNet.h b/Cantera/src/zeroD/ReactorNet.h new file mode 100644 index 000000000..746ac935f --- /dev/null +++ b/Cantera/src/zeroD/ReactorNet.h @@ -0,0 +1,111 @@ +/** + * @file ReactorNet.h + */ + +/* + * $Author$ + * $Revision$ + * $Date$ + */ + +// Copyright 2004 California Institute of Technology + +#ifndef CT_REACTORNET_H +#define CT_REACTORNET_H + +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "Reactor.h" +#include "../FuncEval.h" +#include "../CVode.h" + +namespace Cantera { + + + class ReactorNet : public FuncEval { + + public: + + ReactorNet(); + virtual ~ReactorNet(){} + + //----------------------------------------------------- + + /** @name Methods to set up a simulation. */ + //@{ + + + /** + * Set initial time. Default = 0.0 s. Restarts integration + * from this time using the current mixture state as the + * initial condition. + */ + void setInitialTime(doublereal time) { + m_time = time; + m_init = false; + } + + /** + * Initialize the reactor network. + */ + void initialize(doublereal t0 = 0.0); + + /** + * Advance the state of all reactors in time. + * @param time Time to advance to (s). + */ + void advance(doublereal time); + + double step(doublereal time); + + //@} + + void addReactor(ReactorBase* r) { + m_r.push_back(r); + m_nr++; + } + + ReactorBase& reactor(int n) { + return *m_r[n]; + } + + + /// Return a reference to the integrator. + Integrator& integrator() { return *m_integ; } + + void updateState(doublereal* y); + + //----------------------------------------------------- + + // overloaded methods of class FuncEval + virtual int neq() { return m_nv; } + virtual void eval(doublereal t, doublereal* y, doublereal* ydot); + virtual void getInitialConditions(doublereal t0, size_t leny, + doublereal* y); + + + protected: + + vector m_r; + vector m_reactors; + int m_nr; + int m_nreactors; + Integrator* m_integ; + doublereal m_time; + bool m_init; + int m_nv; + vector_int m_size; + vector_fp m_atol; + doublereal m_rtol; + doublereal m_maxstep; + + private: + + }; +} + +#endif +