*** empty log message ***

This commit is contained in:
Dave Goodwin 2005-11-10 15:02:30 +00:00
parent 12712e54ae
commit ed770b98c8
8 changed files with 294 additions and 42 deletions

View file

@ -18,9 +18,9 @@ do_ranlib = @DO_RANLIB@
CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT)
# stirred reactors
OBJS = Reactor.o ReactorBase.o FlowDevice.o Wall.o ReactorNet.o
OBJS = Reactor.o ReactorBase.o FlowDevice.o Wall.o ReactorNet.o FlowReactor.o
ZEROD_H = Reactor.h ReactorBase.h FlowDevice.h Wall.h ReactorNet.h \
flowControllers.h PID_Controller.h Reservoir.h
flowControllers.h PID_Controller.h Reservoir.h FlowReactor.h
CXX_INCLUDES = -I..
ZEROD_LIB = @buildlib@/libzeroD.a

View file

@ -13,7 +13,7 @@
#endif
#include "Reactor.h"
#include "../CVode.h"
//#include "../CVode.h"
#include "FlowDevice.h"
#include "Wall.h"
#include "../InterfaceKinetics.h"
@ -26,7 +26,9 @@ namespace CanteraZeroD {
doublereal quadInterp(doublereal x0, doublereal* x, doublereal* y);
Reactor::Reactor() : ReactorBase(),
#ifdef INCL_REACTOR_INTEG
FuncEval(),
#endif
m_kin(0),
m_temp_atol(1.e-11),
m_maxstep(0.0),
@ -34,7 +36,7 @@ namespace CanteraZeroD {
m_Q(0.0),
m_rtol(1.e-9),
m_chem(true),
m_energy(true)
m_energy(true), m_nsens(-1)
{
#ifdef INCL_REACTOR_INTEG
m_integ = new CVodeInt;
@ -47,7 +49,6 @@ namespace CanteraZeroD {
#endif
}
// overloaded method of FuncEval. Called by the integrator to
// get the initial conditions.
void Reactor::getInitialConditions(double t0, size_t leny, double* y)
@ -89,7 +90,6 @@ namespace CanteraZeroD {
}
}
/*
* Must be called before calling method 'advance'
*/
@ -111,8 +111,8 @@ namespace CanteraZeroD {
m_pressure = m_thermo->pressure();
m_intEnergy = m_thermo->intEnergy_mass();
int nt = 0, maxnt = 0;
for (int m = 0; m < m_nwalls; m++) {
int m, nt = 0, maxnt = 0;
for (m = 0; m < m_nwalls; m++) {
if (m_wall[m]->kinetics(m_lr[m])) {
nt = m_wall[m]->kinetics(m_lr[m])->nTotalSpecies();
if (nt > maxnt) maxnt = nt;
@ -127,10 +127,23 @@ namespace CanteraZeroD {
}
}
m_work.resize(maxnt);
m_init = true;
}
int Reactor::nSensParams() {
if (m_nsens < 0) {
// determine the number of sensitivity parameters
int m, ns;
m_nsens = m_pnum.size();
for (m = 0; m < m_nwalls; m++) {
ns = m_wall[m]->nSensParams(m_lr[m]);
m_nsens_wall.push_back(ns);
m_nsens += ns;
}
}
return m_nsens;
}
void Reactor::updateState(doublereal* y) {
phase_t& mix = *m_mix; // define for readability
@ -140,6 +153,7 @@ namespace CanteraZeroD {
// Set the mass fractions and density of the mixture.
doublereal u = y[0];
m_vol = y[1];
doublereal* mss = y + 2;
@ -176,22 +190,46 @@ namespace CanteraZeroD {
m_mix->saveState(m_state);
}
#ifdef INCL_REACTOR_INTEG
void Reactor::eval(doublereal time, doublereal* y, doublereal* ydot)
{
updateState(y); // synchronize the reactor state with y
evalEqs(time, y, ydot);
}
#endif
/*
* Called by the integrator to evaluate ydot given y at time 'time'.
*/
void Reactor::evalEqs(doublereal time, doublereal* y, doublereal* ydot)
void Reactor::evalEqs(doublereal time, doublereal* y,
doublereal* ydot, doublereal* params)
{
int i, k, nk;
m_time = time;
m_mix->restoreState(m_state);
Kinetics* kin;
int m, n, npar, ploc;
double mult;
// process sensitivity parameters
if (params) {
npar = m_pnum.size();
for (n = 0; n < npar; n++) {
//m_mult_save[n] = m_kin->multiplier(m_pnum[n]);
mult = m_kin->multiplier(m_pnum[n]);
m_kin->setMultiplier(m_pnum[n], mult*params[n]);
// m_kin->setMultiplier(m_pnum[n], m_mult_save[n]*params[n]);
}
ploc = npar;
for (m = 0; m < m_nwalls; m++) {
if (m_nsens_wall[m] > 0) {
m_wall[m]->setSensitivityParameters(m_lr[m], params + ploc);
ploc += m_nsens_wall[m];
}
}
}
// updateState(y); // synchronize the reactor state with y
m_vdot = 0.0;
@ -199,7 +237,7 @@ namespace CanteraZeroD {
// compute wall terms
doublereal vdot, rs0, sum, wallarea;
Kinetics* kin;
// Kinetics* kin;
SurfPhase* surf;
int lr, ns, loc = m_nsp+2, surfloc;
fill(m_sdot.begin(), m_sdot.end(), 0.0);
@ -242,8 +280,6 @@ namespace CanteraZeroD {
* - \dot m_{out} Y_{k} + A \dot s_k.
*/
const doublereal* mw = m_mix->molecularWeights().begin();
int n;
if (m_chem) {
m_kin->getNetProductionRates(ydot+2); // "omega dot"
}
@ -268,7 +304,7 @@ namespace CanteraZeroD {
ydot[0] = - m_thermo->pressure() * m_vdot - m_Q;
}
else {
ydot[0] = 0.0;
ydot[0] = 0.0;
}
// add terms for open system
@ -287,10 +323,7 @@ namespace CanteraZeroD {
ydot[2+n] -= mdot_out * mf[n];
}
if (m_energy) {
// cout << "before = " << ydot[0] << endl;
ydot[0] -= mdot_out * enthalpy;
//cout << mdot_out << " " << enthalpy << endl;
//cout << "after = " << ydot[0] << endl;
}
}
@ -308,5 +341,59 @@ namespace CanteraZeroD {
}
}
}
// reset sensitivity parameters
if (params) {
npar = m_pnum.size();
for (n = 0; n < npar; n++) {
mult = m_kin->multiplier(m_pnum[n]);
m_kin->setMultiplier(m_pnum[n], mult/params[n]);
//m_kin->setMultiplier(m_pnum[n], m_mult_save[n]);
}
ploc = npar;
for (m = 0; m < m_nwalls; m++) {
if (m_nsens_wall[m] > 0) {
m_wall[m]->resetSensitivityParameters(m_lr[m]);
ploc += m_nsens_wall[m];
}
}
}
}
void Reactor::addSensitivityReaction(int rxn) {
m_pnum.push_back(rxn);
m_pname.push_back(name()+": "+m_kin->reactionString(rxn));
m_mult_save.push_back(1.0);
if (rxn < 0 || rxn >= m_kin->nReactions())
throw CanteraError("Reactor::addSensitivityReaction",
"Reaction number out of range ("+int2str(rxn)+")");
}
int Reactor::componentIndex(string nm) const {
if (nm == "U") return 0;
if (nm == "V") return 1;
// check for a gas species name
int k = m_mix->speciesIndex(nm);
if (k >= 0) return k + 2;
// check for a wall species
int walloffset = 0, kp = 0;
thermo_t* th;
for (int m = 0; m < m_nwalls; m++) {
if (m_wall[m]->kinetics(m_lr[m])) {
kp = m_wall[m]->kinetics(m_lr[m])->reactionPhaseIndex();
th = &m_wall[m]->kinetics(m_lr[m])->thermo(kp);
k = th->speciesIndex(nm);
if (k >= 0) {
return k + 2 + m_nsp + walloffset;
}
else {
walloffset += th->nSpecies();
}
}
}
return -1;
}
}

View file

@ -7,7 +7,7 @@
*/
// Copyright 2001 California Institute of Technology
#ifndef CT_REACTOR_H
#define CT_REACTOR_H
@ -18,10 +18,10 @@
#include "ReactorBase.h"
#include "../FuncEval.h"
#include "../CVode.h"
#include "../Integrator.h"
#include "../Kinetics.h"
#define INCL_REACTOR_INTEG
#undef INCL_REACTOR_INTEG
namespace CanteraZeroD {
@ -57,9 +57,11 @@ namespace CanteraZeroD {
* flow rate. Class FuncEval is the class used to define a system
* of ODE's to be integrated.
*/
#ifdef INCL_REACTOR_INTEG
class Reactor : public ReactorBase, public FuncEval {
#else
class Reactor : public ReactorBase {
#endif
public:
/**
@ -67,7 +69,6 @@ namespace CanteraZeroD {
*/
Reactor();
/**
* Destructor. Deletes the integrator.
*/
@ -137,12 +138,14 @@ namespace CanteraZeroD {
if (m_kin->nReactions() == 0) disableChemistry();
}
#ifdef INCL_REACTOR_INTEG
/**
* Set the maximum step size for integration.
*/
void setMaxStep(doublereal maxstep) {
m_maxstep = maxstep;
}
#endif
void disableChemistry() { m_chem = false; }
void enableChemistry() { m_chem = true; }
@ -153,29 +156,33 @@ namespace CanteraZeroD {
else m_energy = false;
}
//-----------------------------------------------------
/** @name References to internal objects */
//@{
#ifdef INCL_REACTOR_INTEG
/// Return a reference to the integrator.
Integrator& integrator() { return *m_integ; }
//@}
#endif
//-----------------------------------------------------
// overloaded methods of class FuncEval
virtual int neq() { return m_nv; }
#ifdef INCL_REACTOR_INTEG
virtual void eval(doublereal t, doublereal* y, doublereal* ydot);
#endif
virtual void getInitialConditions(doublereal t0, size_t leny,
doublereal* y);
//-----------------------------------------------------
virtual void initialize(doublereal t0 = 0.0);
void evalEqs(doublereal t, doublereal* y, doublereal* ydot);
virtual void evalEqs(doublereal t, doublereal* y, doublereal* ydot, doublereal* params);
/**
* Set the mixture to a state consistent with solution
@ -184,11 +191,24 @@ namespace CanteraZeroD {
virtual void updateState(doublereal* y);
// virtual void addSensitivityParam(int stype, int i);
virtual int nSensParams();
virtual void addSensitivityReaction(int rxn);
virtual string sensParamID(int p) { return m_pname[p]; }
// virtual string component(int k) const;
virtual int componentIndex(string nm) const;
protected:
Kinetics* m_kin;
#ifdef INCL_REACTOR_INTEG
Integrator* m_integ; // pointer to integrator
#endif
doublereal m_temp_atol; // tolerance on T
doublereal m_maxstep; // max step size
doublereal m_vdot, m_Q;
@ -200,6 +220,12 @@ namespace CanteraZeroD {
bool m_energy;
int m_nv;
int m_nsens;
vector_int m_pnum;
vector<string> m_pname;
vector_int m_nsens_wall;
vector_fp m_mult_save;
private:
};
}

View file

@ -32,6 +32,7 @@ namespace CanteraZeroD {
const int ReactorType = 1;
const int ReservoirType = 2;
const int FlowReactorType = 3;
/**
* Base class for stirred reactors.
@ -145,6 +146,8 @@ namespace CanteraZeroD {
return 1;
}
// virtual void addSensitivityParam(int stype, int i) {}
protected:
int m_nsp;

View file

@ -1,15 +1,16 @@
#include "ReactorNet.h"
#include "../CVode.h"
#include "Integrator.h"
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_atols(1.0e-15),
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_verbose(false), m_ntotpar(0)
{
m_integ = new CVodeInt;
m_integ = newIntegrator("CVODE");// CVodeInt;
// use backward differencing, with a full Jacobian computed
// numerically, and use a Newton linear iterator
@ -29,23 +30,34 @@ namespace CanteraZeroD {
writelog("Initializing reactor network.\n");
}
for (n = 0; n < m_nr; n++) {
if (m_r[n]->type() == ReactorType) {
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_atol.resize(neq());
fill(m_atol.begin(), m_atol.end(), m_atols);
m_integ->setTolerances(m_rtol, neq(), m_atol.begin());
m_integ->setSensitivityTolerances(m_rtolsens, m_atolsens);
m_integ->setMaxStepSize(m_maxstep);
if (m_verbose) {
sprintf(buf, "Number of equations: %d\n", neq());
@ -79,20 +91,39 @@ namespace CanteraZeroD {
return m_time;
}
void ReactorNet::eval(doublereal t, doublereal* y, doublereal* ydot) {
// 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);
m_reactors[n]->evalEqs(t, y + start,
ydot + start, p + pstart);
start += m_size[n];
pstart += m_nparams[n];
}
}
catch (CanteraError) {
catch (...) {
showErrors();
error("Terminating execution.");
}
@ -116,5 +147,13 @@ namespace CanteraZeroD {
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);
}
}

View file

@ -20,7 +20,7 @@
#include "Reactor.h"
#include "../FuncEval.h"
#include "../CVode.h"
#include "../Integrator.h"
namespace CanteraZeroD {
@ -59,6 +59,12 @@ namespace CanteraZeroD {
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; }
@ -92,20 +98,33 @@ namespace CanteraZeroD {
return *m_r[n];
}
void setVerbose(bool v = true) { m_verbose = v; }
/// Return a reference to the integrator.
Integrator& integrator() { return *m_integ; }
void updateState(doublereal* y);
double sensitivity(int k, int p) {
return m_integ->sensitivity(k, p)/m_integ->solution(k);
}
double sensitivity(string species, int p, int reactor=0) {
int k = globalComponentIndex(species, reactor);
return sensitivity(k, p);
}
//-----------------------------------------------------
// overloaded methods of class FuncEval
virtual int neq() { return m_nv; }
virtual void eval(doublereal t, doublereal* y, doublereal* ydot);
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(string species, int reactor=0);
protected:
@ -119,10 +138,12 @@ namespace CanteraZeroD {
int m_nv;
vector_int m_size;
vector_fp m_atol;
doublereal m_rtol;
doublereal m_atols;
doublereal m_rtol, m_rtolsens;
doublereal m_atols, m_atolsens;
doublereal m_maxstep;
bool m_verbose;
int m_ntotpar;
vector_int m_nparams;
private:

View file

@ -114,4 +114,58 @@ namespace CanteraZeroD {
else
m_surf[1]->setCoverages(m_rightcov.begin());
}
}
void Wall::addSensitivityReaction(int leftright, int rxn) {
if (rxn < 0 || rxn >= m_chem[leftright]->nReactions())
throw CanteraError("Wall::addSensitivityReaction",
"Reaction number out of range ("+int2str(rxn)+")");
if (leftright == 0) {
m_pleft.push_back(rxn);
m_leftmult_save.push_back(1.0);
m_pname_left.push_back(m_chem[0]->reactionString(rxn));
}
else {
m_pright.push_back(rxn);
m_rightmult_save.push_back(1.0);
m_pname_right.push_back(m_chem[1]->reactionString(rxn));
}
}
void Wall::setSensitivityParameters(int lr, double* params) {
// process sensitivity parameters
int n, npar;
if (lr == 0) {
npar = m_pleft.size();
for (n = 0; n < npar; n++) {
m_leftmult_save[n] = m_chem[0]->multiplier(m_pleft[n]);
m_chem[0]->setMultiplier(m_pleft[n],
m_leftmult_save[n]*params[n]);
}
}
else {
npar = m_pright.size();
for (n = 0; n < npar; n++) {
m_rightmult_save[n] = m_chem[1]->multiplier(m_pright[n]);
m_chem[1]->setMultiplier(m_pright[n],
m_rightmult_save[n]*params[n]);
}
}
}
void Wall::resetSensitivityParameters(int lr) {
int n, npar;
if (lr == 0) {
npar = m_pleft.size();
for (n = 0; n < npar; n++) {
m_chem[0]->setMultiplier(m_pleft[n], m_leftmult_save[n]);
}
}
else {
npar = m_pright.size();
for (n = 0; n < npar; n++) {
m_chem[1]->setMultiplier(m_pright[n],
m_rightmult_save[n]);
}
}
}
}

View file

@ -135,6 +135,24 @@ namespace CanteraZeroD {
void syncCoverages(int leftright);
int nSensParams(int lr) const {
if (lr == 0)
return m_pleft.size();
else
return m_pright.size();
}
void addSensitivityReaction(int leftright, int rxn);
string sensitivityParamID(int leftright, int p) {
if (leftright == 0)
return m_pname_left[p];
else
return m_pname_right[p];
}
void setSensitivityParameters(int lr, double* params);
void resetSensitivityParameters(int lr);
// int componentIndex(string nm) const;
protected:
//vector_fp m_coeffs;
@ -150,6 +168,10 @@ namespace CanteraZeroD {
Func1 *m_qf;
vector_fp m_leftcov, m_rightcov;
vector_int m_pleft, m_pright;
vector_fp m_leftmult_save, m_rightmult_save;
vector<string> m_pname_left, m_pname_right;
private:
};