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

[Reactor] Replace ReactorNet::m_size with m_start

Browse source 
A small simplification since the offset into the global solution vector is
usually what we're interested in.
This commit is contained in:
Ray Speth 2014-06-03 16:51:47 +00:00
parent fb588d873b
commit 2d1003271d
2 changed files with 12 additions and 18 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

5
include/cantera/zeroD/ReactorNet.h
View file

@ -222,7 +222,10 @@ protected:
doublereal m_time; doublereal m_time;
bool m_init; bool m_init;
size_t m_nv; size_t m_nv;
std::vector<size_t> m_size;
//! m_start[n] is the starting point in the state vector for reactor n
std::vector<size_t> m_start;
vector_fp m_atol; vector_fp m_atol;
doublereal m_rtol, m_rtolsens; doublereal m_rtol, m_rtolsens;
doublereal m_atols, m_atolsens; doublereal m_atols, m_atolsens;

25
src/zeroD/ReactorNet.cpp
View file

@ -56,13 +56,13 @@ void ReactorNet::initialize()
throw CanteraError("ReactorNet::initialize", throw CanteraError("ReactorNet::initialize",
"no reactors in network!"); "no reactors in network!");
size_t sensParamNumber = 0; size_t sensParamNumber = 0;
m_start.assign(1, 0);
for (n = 0; n < m_nr; n++) { for (n = 0; n < m_nr; n++) {
if (m_r[n]->type() >= ReactorType) { if (m_r[n]->type() >= ReactorType) {
m_r[n]->initialize(m_time); m_r[n]->initialize(m_time);
Reactor* r = (Reactor*)m_r[n]; Reactor* r = (Reactor*)m_r[n];
m_reactors.push_back(r); m_reactors.push_back(r);
nv = r->neq(); nv = r->neq();
m_size.push_back(nv);
m_nparams.push_back(r->nSensParams()); m_nparams.push_back(r->nSensParams());
std::vector<std::pair<void*, int> > sens_objs = r->getSensitivityOrder(); std::vector<std::pair<void*, int> > sens_objs = r->getSensitivityOrder();
for (size_t i = 0; i < sens_objs.size(); i++) { for (size_t i = 0; i < sens_objs.size(); i++) {
@ -75,6 +75,7 @@ void ReactorNet::initialize()
} }
} }
m_nv += nv; m_nv += nv;
m_start.push_back(m_nv);
m_nreactors++; m_nreactors++;
if (m_verbose) { if (m_verbose) {
@ -189,14 +190,12 @@ void ReactorNet::eval(doublereal t, doublereal* y,
doublereal* ydot, doublereal* p) doublereal* ydot, doublereal* p)
{ {
size_t n; size_t n;
size_t start = 0;
size_t pstart = 0; size_t pstart = 0;
updateState(y); updateState(y);
for (n = 0; n < m_nreactors; n++) { for (n = 0; n < m_nreactors; n++) {
m_reactors[n]->evalEqs(t, y + start, m_reactors[n]->evalEqs(t, y + m_start[n],
ydot + start, p + pstart); ydot + m_start[n], p + pstart);
start += m_size[n];
pstart += m_nparams[n]; pstart += m_nparams[n];
} }
} }
@ -230,20 +229,17 @@ void ReactorNet::evalJacobian(doublereal t, doublereal* y,
void ReactorNet::updateState(doublereal* y) void ReactorNet::updateState(doublereal* y)
{ {
size_t start = 0;
for (size_t n = 0; n < m_nreactors; n++) { for (size_t n = 0; n < m_nreactors; n++) {
m_reactors[n]->updateState(y + start); m_reactors[n]->updateState(y + m_start[n]);
start += m_size[n];
} }
} }
void ReactorNet::getInitialConditions(doublereal t0, void ReactorNet::getInitialConditions(doublereal t0,
size_t leny, doublereal* y) size_t leny, doublereal* y)
{ {
size_t start = 0;
for (size_t n = 0; n < m_nreactors; n++) { for (size_t n = 0; n < m_nreactors; n++) {
m_reactors[n]->getInitialConditions(t0, m_size[n], y + start); m_reactors[n]->getInitialConditions(t0, m_start[n+1]-m_start[n],
start += m_size[n]; y + m_start[n]);
} }
} }
@ -252,12 +248,7 @@ size_t ReactorNet::globalComponentIndex(const string& component, size_t reactor)
if (!m_init) { if (!m_init) {
initialize(); initialize();
} }
size_t start = 0; return m_start[reactor] + m_reactors[reactor]->componentIndex(component);
size_t n;
for (n = 0; n < reactor; n++) {
start += m_size[n];
}
return start + m_reactors[n]->componentIndex(component);
} }
void ReactorNet::registerSensitivityReaction(void* reactor, void ReactorNet::registerSensitivityReaction(void* reactor,