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

[Reactor] Add ReactorSurface to simplify use of surface reactions

Browse source 
This separates the handling of interactions between reactors (mediated by
Wall objects) and surfaces on which surface reactions occur (handled by
ReactorSurface). This simplifies the implementation within reactor, and
reduces the complexity of user code involving surface reactions by
eliminating the need to set up a Reservoir object for the opposite side
of a Wall object that is only being used for surface reactions.
This commit is contained in:
Ray Speth 2016-05-15 21:17:36 -04:00
parent a8c6fe09d9
commit 3f766370b9
8 changed files with 321 additions and 195 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

8
include/cantera/zeroD/ReactorBase.h
View file

@ -13,6 +13,7 @@ namespace Cantera
class FlowDevice;
class Wall;
class ReactorNet;
class ReactorSurface;
const int ReservoirType = 1;
const int ReactorType = 2;
@ -114,6 +115,12 @@ public:
//! Return a reference to the *n*-th Wall connected to this reactor.
Wall& wall(size_t n);
void addSurface(ReactorSurface* surf);
//! Return a reference to the *n*-th ReactorSurface connected to this
//! reactor
ReactorSurface* surface(size_t n);
/**
* Initialize the reactor. Called automatically by ReactorNet::initialize.
*/
@ -230,6 +237,7 @@ protected:
vector_fp m_state;
std::vector<FlowDevice*> m_inlet, m_outlet;
std::vector<Wall*> m_wall;
std::vector<ReactorSurface*> m_surfaces;
vector_int m_lr;
std::string m_name;

91
include/cantera/zeroD/ReactorSurface.h Normal file
View file

@ -0,0 +1,91 @@
//! @file ReactorSurface.h Header file for class ReactorSurface
#ifndef CT_REACTOR_SURFACE_H
#define CT_REACTOR_SURFACE_H
#include "cantera/zeroD/ReactorBase.h"
namespace Cantera
{
class Kinetics;
class SurfPhase;
class ReactorSurface
{
public:
ReactorSurface();
//! Returns the surface area [m^2]
double area() const;
//! Set the surface area [m^2]
void setArea(double a);
//! Accessor for the SurfPhase object
SurfPhase* thermo() {
return m_thermo;
}
//! Accessor for the InterfaceKinetics object
Kinetics* kinetics() {
return m_kinetics;
}
//! Set the InterfaceKinetics object for this surface
void setKinetics(Kinetics* kin);
//! Set the reactor that this Surface interacts with
void setReactor(ReactorBase* reactor);
//! Number of sensitivity parameters associated with reactions on this
//! surface
size_t nSensParams() const {
return m_params.size();
}
//! Set the surface coverages. Array `cov` has length equal to the number of
//! surface species.
void setCoverages(const double* cov);
//! Set the surface coverages by name
void setCoverages(const Composition& cov);
//! Set the surface coverages by name
void setCoverages(const std::string& cov);
//! Get the surface coverages. Array `cov` should have length equal to the
//! number of surface species.
void getCoverages(double* cov) const;
//! Set the coverages in the surface phase object to the values for this
//! surface.
void syncCoverages();
//! Enable calculation of sensitivities with respect to the rate constant
//! for reaction `i`.
void addSensitivityReaction(size_t i);
//! Set reaction rate multipliers. `params` is the global vector of
//! sensitivity parameters. This function is called within
//! ReactorNet::eval() before the reaction rates are evaluated.
void setSensitivityParameters(const double* params);
//! Set reaction rate multipliers back to their initial values. This
//! function is called within ReactorNet::eval() after all rates have been
//! evaluated.
void resetSensitivityParameters();
protected:
double m_area;
SurfPhase* m_thermo;
Kinetics* m_kinetics;
ReactorBase* m_reactor;
vector_fp m_cov;
std::vector<SensitivityParameter> m_params;
};
}
#endif

27
include/cantera/zeroD/Wall.h
View file

@ -8,6 +8,7 @@
#include "cantera/base/ctexceptions.h"
#include "cantera/numerics/Func1.h"
#include "cantera/zeroD/ReactorBase.h"
#include "cantera/zeroD/ReactorSurface.h"
namespace Cantera
{
@ -61,6 +62,8 @@ public:
//! Set the area [m^2].
void setArea(doublereal a) {
m_area = a;
m_surf[0].setArea(a);
m_surf[1].setArea(a);
}
//! Get the area [m^2]
@ -147,13 +150,17 @@ public:
//! Return a pointer to the surface phase object for the left
//! (`leftright=0`) or right (`leftright=1`) wall surface.
SurfPhase* surface(int leftright) {
return m_surf[leftright];
return m_surf[leftright].thermo();
}
ReactorSurface* reactorSurface(int leftright) {
return &m_surf[leftright];
}
//! Return a pointer to the surface kinetics object for the left
//! (`leftright=0`) or right (`leftright=1`) wall surface.
Kinetics* kinetics(int leftright) {
return m_chem[leftright];
return m_surf[leftright].kinetics();
}
//! Set the surface coverages on the left (`leftright = 0`) or right
@ -178,11 +185,7 @@ public:
//! Number of sensitivity parameters associated with reactions on the left
//! (`lr = 0`) or right (`lr = 1`) side of the wall.
size_t nSensParams(int lr) const {
if (lr == 0) {
return m_pleft.size();
} else {
return m_pright.size();
}
return m_surf[lr].nSensParams();
}
void addSensitivityReaction(int leftright, size_t rxn);
void setSensitivityParameters(double* params);
@ -191,17 +194,13 @@ public:
protected:
ReactorBase* m_left;
ReactorBase* m_right;
Kinetics* m_chem[2];
SurfPhase* m_surf[2];
size_t m_nsp[2];
std::vector<ReactorSurface> m_surf;
doublereal m_area, m_k, m_rrth;
doublereal m_emiss;
Func1* m_vf;
Func1* m_qf;
vector_fp m_leftcov, m_rightcov;
std::vector<SensitivityParameter> m_pleft, m_pright;
vector_fp m_leftmult_save, m_rightmult_save;
};
}

15
src/zeroD/ConstPressureReactor.cpp
View file

@ -5,6 +5,7 @@
#include "cantera/zeroD/ConstPressureReactor.h"
#include "cantera/zeroD/FlowDevice.h"
#include "cantera/zeroD/Wall.h"
#include "cantera/thermo/SurfPhase.h"
using namespace std;
@ -162,16 +163,12 @@ std::string ConstPressureReactor::componentName(size_t k) {
} else {
k -= m_thermo->nSpecies();
}
for (size_t m = 0; m < m_wall.size(); m++) {
Wall& w = *m_wall[m];
if (w.kinetics(m_lr[m])) {
size_t kp = w.kinetics(m_lr[m])->reactionPhaseIndex();
ThermoPhase& th = w.kinetics(m_lr[m])->thermo(kp);
if (k < th.nSpecies()) {
return th.speciesName(k);
for (auto& S : m_surfaces) {
ThermoPhase* th = S->thermo();
if (k < th->nSpecies()) {
return th->speciesName(k);
} else {
k -= th.nSpecies();
}
k -= th->nSpecies();
}
}
}

105
src/zeroD/Reactor.cpp
View file

@ -7,6 +7,7 @@
#include "cantera/zeroD/Wall.h"
#include "cantera/thermo/SurfPhase.h"
#include "cantera/zeroD/ReactorNet.h"
#include "cantera/zeroD/ReactorSurface.h"
#include <cfloat>
@ -70,12 +71,9 @@ void Reactor::getState(double* y)
void Reactor::getSurfaceInitialConditions(double* y)
{
size_t loc = 0;
for (size_t m = 0; m < m_wall.size(); m++) {
SurfPhase* surf = m_wall[m]->surface(m_lr[m]);
if (surf) {
m_wall[m]->getCoverages(m_lr[m], y + loc);
loc += surf->nSpecies();
}
for (auto& S : m_surfaces) {
S->getCoverages(y + loc);
loc += S->thermo()->nSpecies();
}
}
@ -88,31 +86,28 @@ void Reactor::initialize(doublereal t0)
m_thermo->restoreState(m_state);
m_sdot.resize(m_nsp, 0.0);
m_wdot.resize(m_nsp, 0.0);
m_nv = m_nsp + 3;
for (size_t w = 0; w < m_wall.size(); w++) {
if (m_wall[w]->surface(m_lr[w])) {
m_nv += m_wall[w]->surface(m_lr[w])->nSpecies();
}
}
m_enthalpy = m_thermo->enthalpy_mass();
m_pressure = m_thermo->pressure();
m_intEnergy = m_thermo->intEnergy_mass();
size_t nt = 0, maxnt = 0;
for (size_t m = 0; m < m_wall.size(); m++) {
m_wall[m]->initialize();
if (m_wall[m]->kinetics(m_lr[m])) {
nt = m_wall[m]->kinetics(m_lr[m])->nTotalSpecies();
for (size_t n = 0; n < m_wall.size(); n++) {
Wall* W = m_wall[n];
W->initialize();
if (W->kinetics(m_lr[n])) {
addSurface(W->reactorSurface(m_lr[n]));
}
}
m_nv = m_nsp + 3;
size_t maxnt = 0;
for (auto& S : m_surfaces) {
m_nv += S->thermo()->nSpecies();
size_t nt = S->kinetics()->nTotalSpecies();
maxnt = std::max(maxnt, nt);
if (m_wall[m]->kinetics(m_lr[m])) {
if (&m_kin->thermo(0) !=
&m_wall[m]->kinetics(m_lr[m])->thermo(0)) {
if (&m_kin->thermo(0) != &S->kinetics()->thermo(0)) {
throw CanteraError("Reactor::initialize",
"First phase of all kinetics managers must be"
" the gas.");
}
}
"First phase of all kinetics managers must be the gas.");
}
}
m_work.resize(maxnt);
@ -121,8 +116,8 @@ void Reactor::initialize(doublereal t0)
size_t Reactor::nSensParams()
{
size_t ns = m_sensParams.size();
for (size_t m = 0; m < m_wall.size(); m++) {
ns += m_wall[m]->nSensParams(m_lr[m]);
for (auto& S : m_surfaces) {
ns += S->nSensParams();
}
return ns;
}
@ -191,12 +186,9 @@ void Reactor::updateState(doublereal* y)
void Reactor::updateSurfaceState(double* y)
{
size_t loc = 0;
for (size_t m = 0; m < m_wall.size(); m++) {
SurfPhase* surf = m_wall[m]->surface(m_lr[m]);
if (surf) {
m_wall[m]->setCoverages(m_lr[m], y+loc);
loc += surf->nSpecies();
}
for (auto& S : m_surfaces) {
S->setCoverages(y+loc);
loc += S->thermo()->nSpecies();
}
}
@ -284,15 +276,15 @@ double Reactor::evalSurfaces(double t, double* ydot)
size_t loc = 0; // offset into ydot
double mdot_surf = 0.0; // net mass flux from surface
for (size_t i = 0; i < m_wall.size(); i++) {
Kinetics* kin = m_wall[i]->kinetics(m_lr[i]);
SurfPhase* surf = m_wall[i]->surface(m_lr[i]);
if (surf && kin) {
for (auto S : m_surfaces) {
Kinetics* kin = S->kinetics();
SurfPhase* surf = S->thermo();
double rs0 = 1.0/surf->siteDensity();
size_t nk = surf->nSpecies();
double sum = 0.0;
surf->setTemperature(m_state[0]);
m_wall[i]->syncCoverages(m_lr[i]);
S->syncCoverages();
kin->getNetProductionRates(&m_work[0]);
size_t ns = kin->surfacePhaseIndex();
size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
@ -303,13 +295,12 @@ double Reactor::evalSurfaces(double t, double* ydot)
ydot[loc] = sum;
loc += nk;
double wallarea = m_wall[i]->area();
double wallarea = S->area();
for (size_t k = 0; k < m_nsp; k++) {
m_sdot[k] += m_work[k]*wallarea;
mdot_surf += m_sdot[k] * mw[k];
}
}
}
return mdot_surf;
}
@ -350,18 +341,14 @@ size_t Reactor::speciesIndex(const string& nm) const
}
// check for a wall species
size_t walloffset = 0, kp = 0;
thermo_t* th;
for (size_t m = 0; m < m_wall.size(); 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);
size_t offset = m_nsp;
for (auto& S : m_surfaces) {
ThermoPhase* th = S->thermo();
k = th->speciesIndex(nm);
if (k != npos) {
return k + m_nsp + walloffset;
return k + offset;
} else {
walloffset += th->nSpecies();
}
offset += th->nSpecies();
}
}
return npos;
@ -412,16 +399,12 @@ std::string Reactor::componentName(size_t k) {
} else {
k -= m_thermo->nSpecies();
}
for (size_t m = 0; m < m_wall.size(); m++) {
Wall& w = *m_wall[m];
if (w.kinetics(m_lr[m])) {
size_t kp = w.kinetics(m_lr[m])->reactionPhaseIndex();
ThermoPhase& th = w.kinetics(m_lr[m])->thermo(kp);
if (k < th.nSpecies()) {
return th.speciesName(k);
for (auto& S : m_surfaces) {
ThermoPhase* th = S->thermo();
if (k < th->nSpecies()) {
return th->speciesName(k);
} else {
k -= th.nSpecies();
}
k -= th->nSpecies();
}
}
}
@ -441,8 +424,8 @@ void Reactor::applySensitivity(double* params)
m_thermo->modifyOneHf298SS(p.local, p.value + params[p.global]);
}
}
for (size_t m = 0; m < m_wall.size(); m++) {
m_wall[m]->setSensitivityParameters(params);
for (auto& S : m_surfaces) {
S->setSensitivityParameters(params);
}
m_thermo->invalidateCache();
m_kin->invalidateCache();
@ -460,8 +443,8 @@ void Reactor::resetSensitivity(double* params)
m_thermo->resetHf298(p.local);
}
}
for (size_t m = 0; m < m_wall.size(); m++) {
m_wall[m]->resetSensitivityParameters();
for (auto& S : m_surfaces) {
S->resetSensitivityParameters();
}
m_thermo->invalidateCache();
m_kin->invalidateCache();

14
src/zeroD/ReactorBase.cpp
View file

@ -5,6 +5,7 @@
#include "cantera/zeroD/ReactorBase.h"
#include "cantera/zeroD/FlowDevice.h"
#include "cantera/zeroD/ReactorNet.h"
#include "cantera/zeroD/ReactorSurface.h"
using namespace std;
namespace Cantera
@ -67,6 +68,19 @@ Wall& ReactorBase::wall(size_t n)
return *m_wall[n];
}
void ReactorBase::addSurface(ReactorSurface* surf)
{
if (find(m_surfaces.begin(), m_surfaces.end(), surf) == m_surfaces.end()) {
m_surfaces.push_back(surf);
surf->setReactor(this);
}
}
ReactorSurface* ReactorBase::surface(size_t n)
{
return m_surfaces[n];
}
ReactorNet& ReactorBase::network()
{
if (m_net) {

106
src/zeroD/ReactorSurface.cpp Normal file
View file

@ -0,0 +1,106 @@
//! @file ReactorSurface.cpp
#include "cantera/zeroD/ReactorSurface.h"
#include "cantera/zeroD/ReactorNet.h"
#include "cantera/thermo/SurfPhase.h"
#include "cantera/kinetics/Kinetics.h"
namespace Cantera
{
ReactorSurface::ReactorSurface()
: m_area(1.0)
, m_thermo(nullptr)
, m_kinetics(nullptr)
, m_reactor(nullptr)
{
}
double ReactorSurface::area() const
{
return m_area;
}
void ReactorSurface::setArea(double a)
{
m_area = a;
}
void ReactorSurface::setKinetics(Kinetics* kin) {
m_kinetics = kin;
if (kin == nullptr) {
m_thermo = nullptr;
return;
}
size_t i = kin->surfacePhaseIndex();
if (i == npos) {
throw CanteraError("ReactorSurface::setKinetics",
"Specified surface kinetics manager does not represent a surface "
"kinetics mechanism.");
}
m_thermo = dynamic_cast<SurfPhase*>(&kin->thermo(i));
m_cov.resize(m_thermo->nSpecies());
m_thermo->getCoverages(m_cov.data());
}
void ReactorSurface::setReactor(ReactorBase* reactor)
{
m_reactor = reactor;
}
void ReactorSurface::setCoverages(const double* cov)
{
copy(cov, cov + m_cov.size(), m_cov.begin());
}
void ReactorSurface::setCoverages(const Composition& cov)
{
m_thermo->setCoveragesByName(cov);
m_thermo->getCoverages(m_cov.data());
}
void ReactorSurface::setCoverages(const std::string& cov)
{
m_thermo->setCoveragesByName(cov);
m_thermo->getCoverages(m_cov.data());
}
void ReactorSurface::getCoverages(double* cov) const
{
copy(m_cov.begin(), m_cov.end(), cov);
}
void ReactorSurface::syncCoverages()
{
m_thermo->setCoverages(m_cov.data());
}
void ReactorSurface::addSensitivityReaction(size_t i)
{
if (i >= m_kinetics->nReactions()) {
throw CanteraError("ReactorSurface::addSensitivityReaction",
"Reaction number out of range ({})", i);
}
size_t p = m_reactor->network().registerSensitivityParameter(
m_kinetics->reactionString(i), 1.0, 1.0);
m_params.emplace_back(
SensitivityParameter{i, p, 1.0, SensParameterType::reaction});
}
void ReactorSurface::setSensitivityParameters(const double* params)
{
for (auto& p : m_params) {
p.value = m_kinetics->multiplier(p.local);
m_kinetics->setMultiplier(p.local, p.value*params[p.global]);
}
}
void ReactorSurface::resetSensitivityParameters()
{
for (auto& p : m_params) {
m_kinetics->setMultiplier(p.local, p.value);
}
}
}

104
src/zeroD/Wall.cpp
View file

@ -8,14 +8,10 @@
namespace Cantera
{
Wall::Wall() : m_left(0), m_right(0),
m_surf(2),
m_area(1.0), m_k(0.0), m_rrth(0.0), m_emiss(0.0),
m_vf(0), m_qf(0)
{
for (int n = 0; n < 2; n++) {
m_chem[n] = 0;
m_surf[n] = 0;
m_nsp[n] = 0;
}
}
bool Wall::install(ReactorBase& rleft, ReactorBase& rright)
@ -28,37 +24,15 @@ bool Wall::install(ReactorBase& rleft, ReactorBase& rright)
m_right = &rright;
m_left->addWall(*this, 0);
m_right->addWall(*this, 1);
m_surf[0].setReactor(&rleft);
m_surf[1].setReactor(&rright);
return true;
}
void Wall::setKinetics(Kinetics* left, Kinetics* right)
{
m_chem[0] = left;
m_chem[1] = right;
size_t ileft = 0, iright = 0;
if (left) {
ileft = left->surfacePhaseIndex();
if (ileft != npos) {
m_surf[0] = (SurfPhase*)&left->thermo(ileft);
m_nsp[0] = m_surf[0]->nSpecies();
m_leftcov.resize(m_nsp[0]);
m_surf[0]->getCoverages(m_leftcov.data());
}
}
if (right) {
iright = right->surfacePhaseIndex();
if (iright != npos) {
m_surf[1] = (SurfPhase*)&right->thermo(iright);
m_nsp[1] = m_surf[1]->nSpecies();
m_rightcov.resize(m_nsp[1]);
m_surf[1]->getCoverages(m_rightcov.data());
}
}
if (ileft == npos || iright == npos) {
throw CanteraError("Wall::setKinetics",
"specified surface kinetics manager does not "
"represent a surface reaction mechanism.");
}
m_surf[0].setKinetics(left);
m_surf[1].setKinetics(right);
}
doublereal Wall::vdot(doublereal t)
@ -88,90 +62,44 @@ doublereal Wall::Q(doublereal t)
void Wall::setCoverages(int leftright, const doublereal* cov)
{
if (leftright == 0) {
copy(cov, cov + m_nsp[0], m_leftcov.begin());
} else {
copy(cov, cov + m_nsp[1], m_rightcov.begin());
}
m_surf[leftright].setCoverages(cov);
}
void Wall::setCoverages(int leftright, const compositionMap& cov)
{
m_surf[leftright]->setCoveragesByName(cov);
if (leftright == 0) {
m_surf[0]->getCoverages(&m_leftcov[0]);
} else {
m_surf[1]->getCoverages(&m_rightcov[0]);
}
m_surf[leftright].setCoverages(cov);
}
void Wall::setCoverages(int leftright, const std::string& cov)
{
m_surf[leftright]->setCoveragesByName(cov);
if (leftright == 0) {
m_surf[0]->getCoverages(&m_leftcov[0]);
} else {
m_surf[1]->getCoverages(&m_rightcov[0]);
}
m_surf[leftright].setCoverages(cov);
}
void Wall::getCoverages(int leftright, doublereal* cov)
{
if (leftright == 0) {
copy(m_leftcov.begin(), m_leftcov.end(), cov);
} else {
copy(m_rightcov.begin(), m_rightcov.end(), cov);
}
m_surf[leftright].getCoverages(cov);
}
void Wall::syncCoverages(int leftright)
{
if (leftright == 0) {
m_surf[0]->setCoverages(m_leftcov.data());
} else {
m_surf[1]->setCoverages(m_rightcov.data());
}
m_surf[leftright].syncCoverages();
}
void Wall::addSensitivityReaction(int leftright, size_t rxn)
{
if (rxn >= m_chem[leftright]->nReactions()) {
throw CanteraError("Wall::addSensitivityReaction",
"Reaction number out of range ({})", rxn);
}
if (leftright == 0) {
size_t p = m_left->network().registerSensitivityParameter(
m_chem[0]->reactionString(rxn), 1.0, 1.0);
m_pleft.emplace_back(
SensitivityParameter{rxn, p, 1.0, SensParameterType::reaction});
} else {
size_t p = m_right->network().registerSensitivityParameter(
m_chem[1]->reactionString(rxn), 1.0, 1.0);
m_pright.emplace_back(
SensitivityParameter{rxn, p, 1.0, SensParameterType::reaction});
}
m_surf[leftright].addSensitivityReaction(rxn);
}
void Wall::setSensitivityParameters(double* params)
{
// process sensitivity parameters
for (auto& p : m_pleft) {
p.value = m_chem[0]->multiplier(p.local);
m_chem[0]->setMultiplier(p.local, p.value*params[p.global]);
}
for (auto& p : m_pright) {
p.value = m_chem[1]->multiplier(p.local);
m_chem[1]->setMultiplier(p.local, p.value*params[p.global]);
}
m_surf[0].setSensitivityParameters(params);
m_surf[1].setSensitivityParameters(params);
}
void Wall::resetSensitivityParameters()
{
for (auto& p : m_pleft) {
m_chem[0]->setMultiplier(p.local, p.value);
}
for (auto& p : m_pright) {
m_chem[1]->setMultiplier(p.local, p.value);
}
m_surf[0].resetSensitivityParameters();
m_surf[1].resetSensitivityParameters();
}
}