Skip to content

Commit

Permalink
Implements MoleReactor and ConstPressureMoleReactor, tests, and docs.
Browse files Browse the repository at this point in the history
This commit implements MoleReactor and ConstPressureMoleReactor classes,
it adds the appropriate python interfaces for them, it also adds a test
comparing surface chemistry of mole reactors to traditional reactors.

test_component_names had a bug because the network was not initialized
so self.net.n_vars was 0 and the loop was never entered. Adding a line
to initialize the network reveal that the test failed for
IdealGasMoleReactor so the appropriate lines were added to fix it.

Add wrappers for Extensible...MoleReactor, update docs, comment fixes.
  • Loading branch information
anthony-walker authored and speth committed Sep 9, 2022
1 parent e3080ed commit a999ead
Show file tree
Hide file tree
Showing 12 changed files with 562 additions and 70 deletions.
24 changes: 24 additions & 0 deletions doc/sphinx/cython/zerodim.rst
Expand Up @@ -42,6 +42,10 @@ Reactor
^^^^^^^
.. autoclass:: Reactor

MoleReactor
^^^^^^^^^^^
.. autoclass:: MoleReactor(contents=None, *, name=None, energy='on')

IdealGasReactor
^^^^^^^^^^^^^^^
.. autoclass:: IdealGasReactor(contents=None, *, name=None, energy='on')
Expand All @@ -54,6 +58,10 @@ ConstPressureReactor
^^^^^^^^^^^^^^^^^^^^
.. autoclass:: ConstPressureReactor(contents=None, *, name=None, energy='on')

ConstPressureMoleReactor
^^^^^^^^^^^^^^^^^^^^^^^^
.. autoclass:: ConstPressureMoleReactor(contents=None, *, name=None, energy='on')

IdealGasConstPressureReactor
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. autoclass:: IdealGasConstPressureReactor(contents=None, *, name=None, energy='on')
Expand Down Expand Up @@ -82,6 +90,22 @@ ExtensibleIdealGasConstPressureReactor
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. autoclass:: ExtensibleIdealGasConstPressureReactor(contents=None, *, name=None, energy='on')

ExtensibleMoleReactor
^^^^^^^^^^^^^^^^^^^^^
.. autoclass:: ExtensibleMoleReactor(contents=None, *, name=None, energy='on')

ExtensibleIdealGasMoleReactor
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. autoclass:: ExtensibleIdealGasMoleReactor(contents=None, *, name=None, energy='on')

ExtensibleConstPressureMoleReactor
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. autoclass:: ExtensibleConstPressureMoleReactor(contents=None, *, name=None, energy='on')

ExtensibleIdealGasConstPressureMoleReactor
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. autoclass:: ExtensibleIdealGasConstPressureMoleReactor(contents=None, *, name=None, energy='on')

Walls
-----

Expand Down
45 changes: 45 additions & 0 deletions include/cantera/zeroD/ConstPressureMoleReactor.h
@@ -0,0 +1,45 @@
//! @file ConstPressureMoleReactor.h

// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.

#ifndef CT_CONSTPRESSMOLE_REACTOR_H
#define CT_CONSTPRESSMOLE_REACTOR_H

#include "cantera/zeroD/MoleReactor.h"

namespace Cantera
{

/*!
* ConstPressureMoleReactor is a class for constant-pressure reactors
* which use a state of moles.
*/
class ConstPressureMoleReactor : public MoleReactor
{
public:
ConstPressureMoleReactor() {}

virtual std::string type() const {
return "ConstPressureMoleReactor";
};

virtual size_t componentIndex(const std::string& nm) const;

virtual std::string componentName(size_t k);

virtual void getState(double* y);

virtual void initialize(double t0 = 0.0);

virtual void eval(double t, double* LHS, double* RHS);

virtual void updateState(double* y);

protected:
const size_t m_sidx = 1;
};

}

#endif
6 changes: 2 additions & 4 deletions include/cantera/zeroD/IdealGasConstPressureMoleReactor.h
Expand Up @@ -6,7 +6,7 @@
#ifndef CT_IDEALGASCONSTPRESSMOLE_REACTOR_H
#define CT_IDEALGASCONSTPRESSMOLE_REACTOR_H

#include "cantera/zeroD/MoleReactor.h"
#include "cantera/zeroD/ConstPressureMoleReactor.h"

namespace Cantera
{
Expand All @@ -15,7 +15,7 @@ namespace Cantera
* IdealGasConstPressureMoleReactor is a class for ideal gas constant-pressure reactors
* which use a state of moles.
*/
class IdealGasConstPressureMoleReactor : public MoleReactor
class IdealGasConstPressureMoleReactor : public ConstPressureMoleReactor
{
public:
IdealGasConstPressureMoleReactor() {}
Expand Down Expand Up @@ -47,8 +47,6 @@ class IdealGasConstPressureMoleReactor : public MoleReactor

protected:
vector_fp m_hk; //!< Species molar enthalpies

const size_t m_sidx = 1;
};

}
Expand Down
27 changes: 16 additions & 11 deletions include/cantera/zeroD/MoleReactor.h
Expand Up @@ -13,8 +13,7 @@ namespace Cantera

/*!
* MoleReactor is meant to serve the same purpose as the reactor class but with a state
* vector composed of moles. It is currently not functional and only serves as a base
* class for other reactors.
* vector composed of moles. It also serves as the base class for other mole reactors.
*/
class MoleReactor : public Reactor
{
Expand All @@ -27,19 +26,25 @@ class MoleReactor : public Reactor

virtual void initialize(double t0 = 0.0);

virtual void eval(double t, double* LHS, double* RHS) {
throw NotImplementedError("MoleReactor::eval()");
}
virtual void getState(double* y);

size_t componentIndex(const std::string& nm) const {
throw NotImplementedError("MoleReactor::componentIndex");
}
virtual void updateState(double* y);

std::string componentName(size_t k) {
throw NotImplementedError("MoleReactor::componentName");
}
virtual void eval(double t, double* LHS, double* RHS);

size_t componentIndex(const std::string& nm) const;

std::string componentName(size_t k);

protected:
//! Get moles of the system from mass fractions stored by thermo object
//! @param y vector for moles to be put into
virtual void getMoles(double* y);

//! Set internal mass variable based on moles given
//! @param y vector of moles of the system
virtual void setMassFromMoles(double* y);

virtual void evalSurfaces(double* LHS, double* RHS, double* sdot);

virtual void updateSurfaceState(double* y);
Expand Down
6 changes: 6 additions & 0 deletions interfaces/cython/cantera/reactor.pxd
Expand Up @@ -193,12 +193,18 @@ cdef class Reactor(ReactorBase):
cdef CxxReactor* reactor
cdef object _kinetics

cdef class MoleReactor(Reactor):
pass

cdef class Reservoir(ReactorBase):
pass

cdef class ConstPressureReactor(Reactor):
pass

cdef class ConstPressureMoleReactor(Reactor):
pass

cdef class IdealGasReactor(Reactor):
pass

Expand Down
47 changes: 47 additions & 0 deletions interfaces/cython/cantera/reactor.pyx
Expand Up @@ -369,6 +369,14 @@ cdef class Reactor(ReactorBase):
limit = -1.
self.reactor.setAdvanceLimit(stringify(name), limit)

cdef class MoleReactor(Reactor):
"""
A homogeneous zero-dimensional reactor with a mole based state vector. By default,
they are closed (no inlets or outlets), have fixed volume, and have adiabatic,
chemically-inert walls. These properties may all be changed by adding
appropriate components such as `Wall`, `MassFlowController` and `Valve`.
"""
reactor_type = "MoleReactor"

cdef class Reservoir(ReactorBase):
"""
Expand All @@ -386,6 +394,13 @@ cdef class ConstPressureReactor(Reactor):
"""
reactor_type = "ConstPressureReactor"

cdef class ConstPressureMoleReactor(Reactor):
"""A homogeneous, constant pressure, zero-dimensional reactor with a mole based
state vector. The volume of the reactor changes as a function of time in order to
keep the pressure constant.
"""
reactor_type = "ConstPressureMoleReactor"


cdef class IdealGasReactor(Reactor):
""" A constant volume, zero-dimensional reactor for ideal gas mixtures. """
Expand Down Expand Up @@ -613,6 +628,38 @@ cdef class ExtensibleIdealGasConstPressureReactor(ExtensibleReactor):
reactor_type = "ExtensibleIdealGasConstPressureReactor"


cdef class ExtensibleMoleReactor(ExtensibleReactor):
"""
A variant of `ExtensibleReactor` where the base behavior corresponds to the
`MoleReactor` class.
"""
reactor_type = "ExtensibleMoleReactor"


cdef class ExtensibleIdealGasMoleReactor(ExtensibleReactor):
"""
A variant of `ExtensibleReactor` where the base behavior corresponds to the
`IdealGasMoleReactor` class.
"""
reactor_type = "ExtensibleIdealGasMoleReactor"


cdef class ExtensibleConstPressureMoleReactor(ExtensibleReactor):
"""
A variant of `ExtensibleReactor` where the base behavior corresponds to the
`ConstPressureMoleReactor` class.
"""
reactor_type = "ExtensibleConstPressureMoleReactor"


cdef class ExtensibleIdealGasConstPressureMoleReactor(ExtensibleReactor):
"""
A variant of `ExtensibleReactor` where the base behavior corresponds to the
`IdealGasConstPressureMoleReactor` class.
"""
reactor_type = "ExtensibleIdealGasConstPressureMoleReactor"


cdef class ReactorSurface:
"""
Represents a surface in contact with the contents of a reactor.
Expand Down
146 changes: 146 additions & 0 deletions src/zeroD/ConstPressureMoleReactor.cpp
@@ -0,0 +1,146 @@
//! @file ConstPressureMoleReactor.cpp A constant pressure
//! zero-dimensional reactor with moles as the state

// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.

#include "cantera/zeroD/Wall.h"
#include "cantera/zeroD/FlowDevice.h"
#include "cantera/zeroD/ConstPressureMoleReactor.h"
#include "cantera/base/utilities.h"
#include "cantera/thermo/SurfPhase.h"
#include "cantera/kinetics/Kinetics.h"

using namespace std;

namespace Cantera
{

void ConstPressureMoleReactor::getState(double* y)
{
if (m_thermo == 0) {
throw CanteraError("ConstPressureMoleReactor::getState",
"Error: reactor is empty.");
}
m_thermo->restoreState(m_state);
// set mass to be used in getMoles function
m_mass = m_thermo->density() * m_vol;
// set the first array element to enthalpy
y[0] = m_thermo->enthalpy_mass() * m_thermo->density() * m_vol;
// get moles of species in remaining state
getMoles(y + m_sidx);
// set the remaining components to the surface species moles on the walls
getSurfaceInitialConditions(y+m_nsp+m_sidx);
}

void ConstPressureMoleReactor::initialize(double t0)
{
MoleReactor::initialize(t0);
m_nv -= 1; // const pressure system loses 1 more variable from MoleReactor
}

void ConstPressureMoleReactor::updateState(double* y)
{
// the components of y are: [0] the enthalpy, [1...K+1) are the
// moles of each species, and [K+1...] are the moles of surface
// species on each wall.
setMassFromMoles(y + m_sidx);
m_thermo->setMolesNoTruncate(y + m_sidx);
if (m_energy) {
m_thermo->setState_HP(y[0] / m_mass, m_pressure);
} else {
m_thermo->setPressure(m_pressure);
}
m_vol = m_mass / m_thermo->density();
updateConnected(false);
updateSurfaceState(y + m_nsp + m_sidx);
}

void ConstPressureMoleReactor::eval(double time, double* LHS, double* RHS)
{
double* dndt = RHS + m_sidx; // kmol per s

evalWalls(time);

m_thermo->restoreState(m_state);

const vector_fp& imw = m_thermo->inverseMolecularWeights();

if (m_chem) {
m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
}

// evaluate reactor surfaces
evalSurfaces(LHS + m_nsp + m_sidx, RHS + m_nsp + m_sidx, m_sdot.data());

// external heat transfer
double dHdt = m_Qdot;

for (size_t n = 0; n < m_nsp; n++) {
// production in gas phase and from surfaces
dndt[n] = m_wdot[n] * m_vol + m_sdot[n];
}

// add terms for outlets
for (auto outlet : m_outlet) {
// determine enthalpy contribution
dHdt -= outlet->massFlowRate() * m_enthalpy;
// flow of species into system and dilution by other species
for (size_t n = 0; n < m_nsp; n++) {
dndt[n] -= outlet->outletSpeciesMassFlowRate(n) * imw[n];
}
}

// add terms for inlets
for (auto inlet : m_inlet) {
// enthalpy contribution from inlets
dHdt += inlet->enthalpy_mass() * inlet->massFlowRate();
// flow of species into system and dilution by other species
for (size_t n = 0; n < m_nsp; n++) {
dndt[n] += inlet->outletSpeciesMassFlowRate(n) * imw[n];
}
}

if (m_energy) {
RHS[0] = dHdt;
} else {
RHS[0] = 0.0;
}
}

size_t ConstPressureMoleReactor::componentIndex(const string& nm) const
{
size_t k = speciesIndex(nm);
if (k != npos) {
return k + m_sidx;
} else if (nm == "enthalpy") {
return 0;
} else {
return npos;
}
}

std::string ConstPressureMoleReactor::componentName(size_t k) {
if (k == 0) {
return "enthalpy";
} else if (k >= m_sidx && k < neq()) {
k -= m_sidx;
if (k < m_thermo->nSpecies()) {
return m_thermo->speciesName(k);
} else {
k -= m_thermo->nSpecies();
}
for (auto& S : m_surfaces) {
ThermoPhase* th = S->thermo();
if (k < th->nSpecies()) {
return th->speciesName(k);
} else {
k -= th->nSpecies();
}
}
}
throw CanteraError("ConstPressureMoleReactor::componentName",
"Index is out of bounds.");
}

}

0 comments on commit a999ead

Please sign in to comment.