diff --git a/Cantera/src/numerics/Integrator.h b/Cantera/src/numerics/Integrator.h index c2c0f1800..7446714b3 100755 --- a/Cantera/src/numerics/Integrator.h +++ b/Cantera/src/numerics/Integrator.h @@ -183,6 +183,9 @@ namespace Cantera { // defined in ODE_integrators.cpp Integrator* newIntegrator(std::string itype); + // defined in ODE_integrators.cpp + void deleteIntegrator(Integrator *cv); + } // namespace #endif diff --git a/Cantera/src/numerics/ODE_integrators.cpp b/Cantera/src/numerics/ODE_integrators.cpp index 53ba1980f..2846b5ace 100644 --- a/Cantera/src/numerics/ODE_integrators.cpp +++ b/Cantera/src/numerics/ODE_integrators.cpp @@ -23,4 +23,8 @@ namespace Cantera { "unknown ODE integrator: "+itype); } } + + void deleteIntegrator(Integrator *cv) { + delete cv; + } } diff --git a/Cantera/src/zeroD/ReactorNet.cpp b/Cantera/src/zeroD/ReactorNet.cpp index 1ec5f4298..62684e294 100644 --- a/Cantera/src/zeroD/ReactorNet.cpp +++ b/Cantera/src/zeroD/ReactorNet.cpp @@ -7,231 +7,261 @@ using namespace std; namespace CanteraZeroD { - ReactorNet::ReactorNet() : FuncEval(), m_nr(0), m_nreactors(0), - m_integ(0), m_time(0.0), m_init(false), - m_nv(0), m_rtol(1.0e-9), m_rtolsens(1.0e-4), - m_atols(1.0e-15), m_atolsens(1.0e-4), - m_maxstep(-1.0), - m_verbose(false), m_ntotpar(0) - { + ReactorNet::ReactorNet() : FuncEval(), m_nr(0), m_nreactors(0), + m_integ(0), m_time(0.0), m_init(false), + m_nv(0), m_rtol(1.0e-9), m_rtolsens(1.0e-4), + m_atols(1.0e-15), m_atolsens(1.0e-4), + m_maxstep(-1.0), + m_verbose(false), m_ntotpar(0) + { #ifdef DEBUG_MODE - m_verbose = true; + m_verbose = true; #endif - m_integ = newIntegrator("CVODE");// CVodeInt; + m_integ = newIntegrator("CVODE");// CVodeInt; - // use backward differencing, with a full Jacobian computed - // numerically, and use a Newton linear iterator + // 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->setMethod(BDF_Method); + m_integ->setProblemType(DENSE + NOJAC); + m_integ->setIterator(Newton_Iter); + } + + ReactorNet::~ReactorNet() { + for (int n = 0; n < m_nr; n++) { + if (m_iown[n]) { + delete m_r[n]; + } + m_r[n] = 0; + } + m_r.clear(); + m_reactors.clear(); + deleteIntegrator(m_integ); + } + + void ReactorNet::initialize(doublereal t0) { + int n, nv; + char buf[100]; + m_nv = 0; + m_reactors.clear(); + m_nreactors = 0; + if (m_verbose) { + writelog("Initializing reactor network.\n"); + } + if (m_nr == 0) + throw CanteraError("ReactorNet::initialize", + "no reactors in network!"); + for (n = 0; n < m_nr; n++) { + if (m_r[n]->type() >= ReactorType) { + m_r[n]->initialize(t0); + Reactor* r = (Reactor*)m_r[n]; + m_reactors.push_back(r); + nv = r->neq(); + m_size.push_back(nv); + m_nparams.push_back(r->nSensParams()); + m_ntotpar += r->nSensParams(); + m_nv += nv; + m_nreactors++; + + if (m_verbose) { + sprintf(buf,"Reactor %d: %d variables.\n",n,nv); + writelog(buf); + sprintf(buf," %d sensitivity params.\n", + r->nSensParams()); + writelog(buf); + } + if (m_r[n]->type() == FlowReactorType && m_nr > 1) { + throw CanteraError("ReactorNet::initialize", + "FlowReactors must be used alone."); + } + } } - void ReactorNet::initialize(doublereal t0) { - int n, nv; - char buf[100]; - m_nv = 0; - m_reactors.clear(); - m_nreactors = 0; - if (m_verbose) { - writelog("Initializing reactor network.\n"); - } - if (m_nr == 0) - throw CanteraError("ReactorNet::initialize", - "no reactors in network!"); - for (n = 0; n < m_nr; n++) { - if (m_r[n]->type() >= ReactorType) { - m_r[n]->initialize(t0); - Reactor* r = (Reactor*)m_r[n]; - m_reactors.push_back(r); - nv = r->neq(); - m_size.push_back(nv); - m_nparams.push_back(r->nSensParams()); - m_ntotpar += r->nSensParams(); - m_nv += nv; - m_nreactors++; - - if (m_verbose) { - sprintf(buf,"Reactor %d: %d variables.\n",n,nv); - writelog(buf); - sprintf(buf," %d sensitivity params.\n", - r->nSensParams()); - writelog(buf); - } - if (m_r[n]->type() == FlowReactorType && m_nr > 1) { - throw CanteraError("ReactorNet::initialize", - "FlowReactors must be used alone."); - } - } - } - - m_connect.resize(m_nr*m_nr,0); - m_ydot.resize(m_nv,0.0); - int i, j, nin, nout, nw; - ReactorBase *r, *rj; - for (i = 0; i < m_nr; i++) { - r = m_reactors[i]; - for (j = 0; j < m_nr; j++) { - if (i == j) connect(i,j); - else { - rj = m_reactors[j]; - nin = rj->nInlets(); - for (n = 0; n < nin; n++) { - if (&rj->inlet(n).out() == r) connect(i,j); - } - nout = rj->nOutlets(); - for (n = 0; n < nout; n++) { - if (&rj->outlet(n).in() == r) connect(i,j); - } - nw = rj->nWalls(); - for (n = 0; n < nw; n++) { - if (&rj->wall(n).left() == rj - && &rj->wall(n).right() == r) connect(i,j); - else if (&rj->wall(n).left() == r - && &rj->wall(n).right() == rj) connect(i,j); - } - } - } - } - - m_atol.resize(neq()); - fill(m_atol.begin(), m_atol.end(), m_atols); - m_integ->setTolerances(m_rtol, neq(), DATA_PTR(m_atol)); - m_integ->setSensitivityTolerances(m_rtolsens, m_atolsens); - m_integ->setMaxStepSize(m_maxstep); - if (m_verbose) { - sprintf(buf, "Number of equations: %d\n", neq()); - writelog(buf); - sprintf(buf, "Maximum time step: %14.6g\n", m_maxstep); - writelog(buf); - } - m_integ->initialize(t0, *this); - m_init = true; + m_connect.resize(m_nr*m_nr,0); + m_ydot.resize(m_nv,0.0); + int i, j, nin, nout, nw; + ReactorBase *r, *rj; + for (i = 0; i < m_nr; i++) { + r = m_reactors[i]; + for (j = 0; j < m_nr; j++) { + if (i == j) connect(i,j); + else { + rj = m_reactors[j]; + nin = rj->nInlets(); + for (n = 0; n < nin; n++) { + if (&rj->inlet(n).out() == r) connect(i,j); + } + nout = rj->nOutlets(); + for (n = 0; n < nout; n++) { + if (&rj->outlet(n).in() == r) connect(i,j); + } + nw = rj->nWalls(); + for (n = 0; n < nw; n++) { + if (&rj->wall(n).left() == rj + && &rj->wall(n).right() == r) connect(i,j); + else if (&rj->wall(n).left() == r + && &rj->wall(n).right() == rj) connect(i,j); + } + } + } } - void ReactorNet::advance(doublereal time) { - if (!m_init) { - if (m_maxstep < 0.0) - m_maxstep = time - m_time; - initialize(); - } - m_integ->integrate(time); - m_time = time; - updateState(m_integ->solution()); + m_atol.resize(neq()); + fill(m_atol.begin(), m_atol.end(), m_atols); + m_integ->setTolerances(m_rtol, neq(), DATA_PTR(m_atol)); + m_integ->setSensitivityTolerances(m_rtolsens, m_atolsens); + m_integ->setMaxStepSize(m_maxstep); + if (m_verbose) { + sprintf(buf, "Number of equations: %d\n", neq()); + writelog(buf); + sprintf(buf, "Maximum time step: %14.6g\n", m_maxstep); + writelog(buf); } + m_integ->initialize(t0, *this); + m_init = true; + } - double ReactorNet::step(doublereal time) { - if (!m_init) { - if (m_maxstep < 0.0) - m_maxstep = time - m_time; - initialize(); - } - m_time = m_integ->step(time); - updateState(m_integ->solution()); - return m_time; + void ReactorNet::advance(doublereal time) { + if (!m_init) { + if (m_maxstep < 0.0) + m_maxstep = time - m_time; + initialize(); } + m_integ->integrate(time); + m_time = time; + updateState(m_integ->solution()); + } -// void ReactorNet::addSensitivityParam(int n, int stype, int i) { -// m_reactors[n]->addSensitivityParam(int stype, int i); -// m_sensreactor.push_back(n); -// m_nSenseParams++; -// } + double ReactorNet::step(doublereal time) { + if (!m_init) { + if (m_maxstep < 0.0) + m_maxstep = time - m_time; + initialize(); + } + m_time = m_integ->step(time); + updateState(m_integ->solution()); + return m_time; + } -// void ReactorNet::setParameters(int np, double* p) { -// int n, nr; -// for (n = 0; n < np; n++) { -// if (n < m_nSenseParams) { -// nr = m_sensreactor[n]; -// m_reactors[nr]->setParameter(n, p[n]); -// } -// } -// } + + void ReactorNet::addReactor(ReactorBase* r, bool iown) { + if (r->type() >= ReactorType) { + m_r.push_back(r); + m_iown.push_back(iown); + m_nr++; + if (m_verbose) { + writelog("Adding reactor "+r->name()+"\n"); + } + } + else { + if (m_verbose) { + writelog("Not adding reactor "+r->name()+ + ", since type = "+int2str(r->type())+"\n"); + } + } + } + + // void ReactorNet::addSensitivityParam(int n, int stype, int i) { + // m_reactors[n]->addSensitivityParam(int stype, int i); + // m_sensreactor.push_back(n); + // m_nSenseParams++; + // } + + // void ReactorNet::setParameters(int np, double* p) { + // int n, nr; + // for (n = 0; n < np; n++) { + // if (n < m_nSenseParams) { + // nr = m_sensreactor[n]; + // m_reactors[nr]->setParameter(n, p[n]); + // } + // } + // } - void ReactorNet::eval(doublereal t, doublereal* y, - doublereal* ydot, doublereal* p) { - int n; - int start = 0; - int pstart = 0; - // use a try... catch block, since exceptions are not passed - // through CVODE, since it is C code - try { - updateState(y); - for (n = 0; n < m_nreactors; n++) { - m_reactors[n]->evalEqs(t, y + start, - ydot + start, p + pstart); - start += m_size[n]; - pstart += m_nparams[n]; - } - } - catch (...) { - showErrors(); - error("Terminating execution."); - } + void ReactorNet::eval(doublereal t, doublereal* y, + doublereal* ydot, doublereal* p) { + int n; + int start = 0; + int pstart = 0; + // use a try... catch block, since exceptions are not passed + // through CVODE, since it is C code + try { + updateState(y); + for (n = 0; n < m_nreactors; n++) { + m_reactors[n]->evalEqs(t, y + start, + ydot + start, p + pstart); + start += m_size[n]; + pstart += m_nparams[n]; + } } + catch (...) { + showErrors(); + error("Terminating execution."); + } + } - void ReactorNet::evalJacobian(doublereal t, doublereal* y, - doublereal* ydot, doublereal* p, Array2D* j) { - int n, m; - doublereal ysave, dy; - Array2D& jac = *j; + void ReactorNet::evalJacobian(doublereal t, doublereal* y, + doublereal* ydot, doublereal* p, Array2D* j) { + int n, m; + doublereal ysave, dy; + Array2D& jac = *j; - // use a try... catch block, since exceptions are not passed - // through CVODE, since it is C code - try { - //evaluate the unperturbed ydot - eval(t, y, ydot, p); - for (n = 0; n < m_nv; n++) { + // use a try... catch block, since exceptions are not passed + // through CVODE, since it is C code + try { + //evaluate the unperturbed ydot + eval(t, y, ydot, p); + for (n = 0; n < m_nv; n++) { - // perturb x(n) - ysave = y[n]; - dy = m_atol[n] + fabs(ysave)*m_rtol; - y[n] = ysave + dy; - dy = y[n] - ysave; + // perturb x(n) + ysave = y[n]; + dy = m_atol[n] + fabs(ysave)*m_rtol; + y[n] = ysave + dy; + dy = y[n] - ysave; - // calculate perturbed residual - eval(t, y, DATA_PTR(m_ydot), p); + // calculate perturbed residual + eval(t, y, DATA_PTR(m_ydot), p); - // compute nth column of Jacobian - for (m = 0; m < m_nv; m++) { - jac(m,n) = (m_ydot[m] - ydot[m])/dy; - } - y[n] = ysave; - } - } - catch (...) { - showErrors(); - error("Terminating execution."); - } + // compute nth column of Jacobian + for (m = 0; m < m_nv; m++) { + jac(m,n) = (m_ydot[m] - ydot[m])/dy; + } + y[n] = ysave; + } } - - 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]; - } + catch (...) { + showErrors(); + error("Terminating execution."); } + } - 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]; - } + 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]; } + } - int ReactorNet::globalComponentIndex(string species, int reactor) { - int start = 0; - int n; - for (n = 0; n < reactor; n++) start += m_size[n]; - return start + m_reactors[n]->componentIndex(species); + 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]; } + } + + int ReactorNet::globalComponentIndex(string species, int reactor) { + int start = 0; + int n; + for (n = 0; n < reactor; n++) start += m_size[n]; + return start + m_reactors[n]->componentIndex(species); + } } diff --git a/Cantera/src/zeroD/ReactorNet.h b/Cantera/src/zeroD/ReactorNet.h index 3ed3f6c9f..1cfdf9ae3 100644 --- a/Cantera/src/zeroD/ReactorNet.h +++ b/Cantera/src/zeroD/ReactorNet.h @@ -26,153 +26,144 @@ namespace CanteraZeroD { - class ReactorNet : public FuncEval { + class ReactorNet : public FuncEval { - public: + public: - ReactorNet(); - virtual ~ReactorNet(){} + //! Constructor + ReactorNet(); - //----------------------------------------------------- + //! Destructor + 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; - } + /** @name Methods to set up a simulation. */ + //@{ - /// Set the maximum time step. - void setMaxTimeStep(double maxstep) { - m_maxstep = maxstep; - m_init = false; - } + /** + * 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; + } + + /// Set the maximum time step. + void setMaxTimeStep(double maxstep) { + m_maxstep = maxstep; + m_init = false; + } - void setTolerances(doublereal rtol, doublereal atol) { - if (rtol >= 0.0) m_rtol = rtol; - if (atol >= 0.0) m_atols = atol; - m_init = false; - } + void setTolerances(doublereal rtol, doublereal atol) { + if (rtol >= 0.0) m_rtol = rtol; + if (atol >= 0.0) m_atols = atol; + m_init = false; + } - void setSensitivityTolerances(doublereal rtol, doublereal atol) { - if (rtol >= 0.0) m_rtolsens = rtol; - if (atol >= 0.0) m_atolsens = atol; - m_init = false; - } + void setSensitivityTolerances(doublereal rtol, doublereal atol) { + if (rtol >= 0.0) m_rtolsens = rtol; + if (atol >= 0.0) m_atolsens = atol; + m_init = false; + } - /// Current value of the simulation time. - doublereal time() { return m_time; } + /// Current value of the simulation time. + doublereal time() { return m_time; } - /// Relative tolerance. - doublereal rtol() { return m_rtol; } - doublereal atol() { return m_atols; } + /// Relative tolerance. + doublereal rtol() { return m_rtol; } + doublereal atol() { return m_atols; } - /** - * Initialize the reactor network. - */ - void initialize(doublereal t0 = 0.0); + /** + * 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); + /** + * Advance the state of all reactors in time. + * @param time Time to advance to (s). + */ + void advance(doublereal time); - double step(doublereal time); + double step(doublereal time); - //@} + //@} - void addReactor(ReactorBase* r) { - if (r->type() >= ReactorType) { - m_r.push_back(r); - m_nr++; - if (m_verbose) { - writelog("Adding reactor "+r->name()+"\n"); - } - } - else { - if (m_verbose) { - writelog("Not adding reactor "+r->name()+ - ", since type = "+int2str(r->type())+"\n"); - } - } - } + void addReactor(ReactorBase* r, bool iown = false); - ReactorBase& reactor(int n) { - return *m_r[n]; - } + ReactorBase& reactor(int n) { + return *m_r[n]; + } - bool verbose() const { return m_verbose; } - void setVerbose(bool v = true) { m_verbose = v; } + bool verbose() const { return m_verbose; } + void setVerbose(bool v = true) { m_verbose = v; } - /// Return a reference to the integrator. - Integrator& integrator() { return *m_integ; } + /// Return a reference to the integrator. + Integrator& integrator() { return *m_integ; } - void updateState(doublereal* y); + void updateState(doublereal* y); - double sensitivity(int k, int p) { - return m_integ->sensitivity(k, p)/m_integ->solution(k); - } + double sensitivity(int k, int p) { + return m_integ->sensitivity(k, p)/m_integ->solution(k); + } - double sensitivity(std::string species, int p, int reactor=0) { - int k = globalComponentIndex(species, reactor); - return sensitivity(k, p); - } + double sensitivity(std::string species, int p, int reactor=0) { + int k = globalComponentIndex(species, reactor); + return sensitivity(k, p); + } - void evalJacobian(doublereal t, doublereal* y, - doublereal* ydot, doublereal* p, Array2D* j); + void evalJacobian(doublereal t, doublereal* y, + doublereal* ydot, doublereal* p, Array2D* j); - //----------------------------------------------------- + //----------------------------------------------------- - // overloaded methods of class FuncEval - virtual int neq() { return m_nv; } - virtual void eval(doublereal t, doublereal* y, - doublereal* ydot, doublereal* p); - virtual void getInitialConditions(doublereal t0, size_t leny, - doublereal* y); - virtual int nparams() { return m_ntotpar; } + // overloaded methods of class FuncEval + virtual int neq() { return m_nv; } + virtual void eval(doublereal t, doublereal* y, + doublereal* ydot, doublereal* p); + virtual void getInitialConditions(doublereal t0, size_t leny, + doublereal* y); + virtual int nparams() { return m_ntotpar; } - int globalComponentIndex(std::string species, int reactor=0); + int globalComponentIndex(std::string species, int reactor=0); - void connect(int i, int j) { - m_connect[j*m_nr + i] = 1; - m_connect[i*m_nr + j] = 1; - } + void connect(int i, int j) { + m_connect[j*m_nr + i] = 1; + m_connect[i*m_nr + j] = 1; + } - bool connected(int i, int j) { - return (m_connect[m_nr*i + j] == 1); - } + bool connected(int i, int j) { + return (m_connect[m_nr*i + j] == 1); + } - protected: + protected: - std::vector m_r; - std::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, m_rtolsens; - doublereal m_atols, m_atolsens; - doublereal m_maxstep; - bool m_verbose; - int m_ntotpar; - vector_int m_nparams; - vector_int m_connect; - vector_fp m_ydot; + std::vector m_r; + std::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, m_rtolsens; + doublereal m_atols, m_atolsens; + doublereal m_maxstep; + bool m_verbose; + int m_ntotpar; + vector_int m_nparams; + vector_int m_connect; + vector_fp m_ydot; - private: + std::vector m_iown; - }; + private: + + }; } #endif