research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

*** empty log message ***

Browse source
This commit is contained in:
Dave Goodwin 2004-03-03 06:31:14 +00:00
parent f957c14af4
commit 183c0b3409
12 changed files with 513 additions and 147 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

123
Cantera/clib/src/ctreactor.cpp
View file

@ -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<reactor_t>* Cabinet<reactor_t>::__storage = 0;
Cabinet<reactornet_t>* Cabinet<reactornet_t>::__storage = 0;
Cabinet<flowdev_t>* Cabinet<flowdev_t>::__storage = 0;
Cabinet<wall_t>* Cabinet<wall_t>::__storage = 0;
@ -33,6 +36,10 @@ inline reactor_t* _reactor(int i) {
return Cabinet<reactor_t>::cabinet()->item(i);
}
inline reactornet_t* _reactornet(int i) {
return Cabinet<reactornet_t>::cabinet()->item(i);
}
inline flowdev_t* _flowdev(int i) {
return Cabinet<flowdev_t>::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<reactornet_t>::cabinet()->add(r);
}
int DLL_EXPORT reactornet_del(int i) {
try {
Cabinet<reactornet_t>::cabinet()->del(i);
return 0;
}
catch (...) {
return -1;
}
}
int DLL_EXPORT reactornet_copy(int i) {
return Cabinet<reactornet_t>::cabinet()->newCopy(i);
}
int DLL_EXPORT reactornet_assign(int i, int j) {
return Cabinet<reactornet_t>::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);
}

9
Cantera/clib/src/ctreactor.h
View file

@ -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);

11
Cantera/matlab/cantera/polynom.m
View file

@ -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

64
Cantera/python/Cantera/Reactor.py
View file

@ -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:

68
Cantera/python/Cantera/ReactorNet.py Normal file
View file

@ -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)

62
Cantera/python/src/ctreactor_methods.cpp
View file

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

6
Cantera/python/src/methods.h
View file

@ -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},

2
Cantera/src/zeroD/Makefile.in
View file

@ -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

20
Cantera/src/zeroD/Reactor.cpp
View file

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

94
Cantera/src/zeroD/Reactor.h
View file

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

90
Cantera/src/zeroD/ReactorNet.cpp Normal file
View file

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

111
Cantera/src/zeroD/ReactorNet.h Normal file
View file

@ -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<ReactorBase*> m_r;
vector<Reactor*> 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