From 1b5cb699c60dd8718f03935632f6e0bd30c325da Mon Sep 17 00:00:00 2001 From: Gandhali Kogekar Date: Thu, 23 May 2024 12:01:21 -0400 Subject: [PATCH] Implementation of 'ReactorEdge' class --- include/cantera/zeroD/ReactorBase.h | 12 ++ include/cantera/zeroD/ReactorDelegator.h | 1 + include/cantera/zeroD/ReactorEdge.h | 81 +++++++++++++ include/cantera/zerodim.h | 3 + interfaces/cython/cantera/drawnetwork.py | 46 +++++++ interfaces/cython/cantera/reactor.pxd | 25 ++++ interfaces/cython/cantera/reactor.pyx | 112 ++++++++++++++++++ .../sourcegen/sourcegen/csharp/config.yaml | 1 + src/clib/ctreactor.cpp | 63 ++++++++++ src/zeroD/Reactor.cpp | 2 + src/zeroD/ReactorBase.cpp | 14 +++ src/zeroD/ReactorEdge.cpp | 76 ++++++++++++ 12 files changed, 436 insertions(+) create mode 100644 include/cantera/zeroD/ReactorEdge.h create mode 100644 src/zeroD/ReactorEdge.cpp diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h index 08c1e978e4..2e9b9921ef 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/include/cantera/zeroD/ReactorBase.h @@ -24,6 +24,7 @@ class FlowDevice; class WallBase; class ReactorNet; class ReactorSurface; +class ReactorEdge; class Kinetics; class ThermoPhase; class Solution; @@ -163,6 +164,16 @@ class ReactorBase return m_surfaces.size(); } + virtual void addEdge(ReactorEdge* edge); + + //! Return a reference to the *n*-th ReactorEdge connected to this + //! reactor + ReactorEdge* edge(size_t n); + //! Return the number of edges in a reactor + virtual size_t nEdges() { + return m_edges.size(); + } + /** * Initialize the reactor. Called automatically by ReactorNet::initialize. */ @@ -303,6 +314,7 @@ class ReactorBase vector m_wall; vector m_surfaces; + vector m_edges; //! Vector of length nWalls(), indicating whether this reactor is on the left (0) //! or right (1) of each wall. diff --git a/include/cantera/zeroD/ReactorDelegator.h b/include/cantera/zeroD/ReactorDelegator.h index 40c5f40870..ca4e97b28e 100644 --- a/include/cantera/zeroD/ReactorDelegator.h +++ b/include/cantera/zeroD/ReactorDelegator.h @@ -9,6 +9,7 @@ #include "Reactor.h" #include "cantera/base/Delegator.h" #include "cantera/zeroD/ReactorSurface.h" +#include "cantera/zeroD/ReactorEdge.h" #include "cantera/thermo/SurfPhase.h" namespace Cantera diff --git a/include/cantera/zeroD/ReactorEdge.h b/include/cantera/zeroD/ReactorEdge.h new file mode 100644 index 0000000000..fd7b3af8c1 --- /dev/null +++ b/include/cantera/zeroD/ReactorEdge.h @@ -0,0 +1,81 @@ +//! @file ReactorEdge.h Header file for class ReactorEdge + +// 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_REACTOR_EDGE_H +#define CT_REACTOR_EDGE_H + +#include "cantera/zeroD/ReactorBase.h" + +namespace Cantera +{ + +class Kinetics; +//class SurfPhase; +class EdgePhase; + +//! An edge where reactions can occur that is in contact with the two surfaces of a +//! Reactor. +//! @ingroup wallGroup +class ReactorEdge +{ +public: + ReactorEdge() = default; + virtual ~ReactorEdge() = default; + ReactorEdge(const ReactorEdge&) = delete; + ReactorEdge& operator=(const ReactorEdge&) = delete; + + //! Returns the edge length [m] + double length() const; + + //! Set the edge length [m] + void setLength(double a); + + //! Accessor for the EdgePhase object + EdgePhase* thermo() { + return m_thermo; + } + + //! Accessor for the InterfaceKinetics object + Kinetics* kinetics() { + return m_kinetics; + } + + //! Set the InterfaceKinetics object for the edge + void setKinetics(Kinetics* kin); + + //! Set the reactor that this edge interacts with + void setReactor(ReactorBase* reactor); + + //! 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 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 and temperature in the surface phase object to the + //! values for this surface. The temperature is set to match the bulk phase + //! of the attached Reactor. + void syncState(); + +protected: + double m_length = 1.0; + + EdgePhase* m_thermo = nullptr; + Kinetics* m_kinetics = nullptr; + ReactorBase* m_reactor = nullptr; + vector m_cov; +}; + +} + +#endif diff --git a/include/cantera/zerodim.h b/include/cantera/zerodim.h index eb80929fe6..52291e7e18 100644 --- a/include/cantera/zerodim.h +++ b/include/cantera/zerodim.h @@ -31,6 +31,9 @@ // surface #include "cantera/zeroD/ReactorSurface.h" +// edge +#include "cantera/zeroD/ReactorEdge.h" + // factories #include "cantera/zeroD/ReactorFactory.h" #include "cantera/zeroD/FlowDeviceFactory.h" diff --git a/interfaces/cython/cantera/drawnetwork.py b/interfaces/cython/cantera/drawnetwork.py index cafb980c57..adb67a409e 100644 --- a/interfaces/cython/cantera/drawnetwork.py +++ b/interfaces/cython/cantera/drawnetwork.py @@ -282,6 +282,52 @@ def draw_surface(surface, graph=None, graph_attr=None, node_attr=None, return graph +""" +def draw_edge(edge, graph=None, graph_attr=None, node_attr=None, + surface_edge_attr=None, print_state=False, **kwargs): + + Draw `~cantera.ReactorEdge` object with its connected reactor. + + :param surface: + `~cantera.ReactorEdge` object. + :param graph: + ``graphviz.graphs.BaseGraph`` object to which the connection is added. + If not provided, a new ``DiGraph`` is created. Defaults to ``None``. + :param graph_attr: + Attributes to be passed to the ``graphviz.Digraph`` function that control the + general appearance of the drawn network. + Has no effect if existing ``graph`` is provided. + See https://graphviz.org/docs/graph/ for a list of all usable attributes. + :param node_attr: + Attributes to be passed to the ``node`` method invoked to draw the reactor. + See https://graphviz.org/docs/nodes/ for a list of all usable attributes. + :param surface_edge_attr: + Attributes to be passed to the ``edge`` method invoked to draw the connection + between the surface and its reactor. + See https://graphviz.org/docs/edges/ for a list of all usable attributes. + Default is ``{"style": "dotted", "arrowhead": "none"}``. + :param print_state: + Whether state information of the reactor is printed into the node. + See ``draw_reactor`` for additional keywords to control this output. + Defaults to ``False``. + :param kwargs: + Additional keywords are passed on to `draw_reactor`. + :return: + A ``graphviz.graphs.BaseGraph`` object depicting the surface and its reactor. + + .. versionadded:: 3.1???? + + r = edge.reactor + graph = draw_reactor(r, graph, graph_attr, node_attr, print_state, **kwargs) + name = f"{r.name} edge" + edge_attr = {"style": "dotted", "arrowhead": "none", + **(surface_edge_attr or {})} + + graph.node(name, **dict(node_attr or {}, **edge.node_attr)) + graph.edge(r.name, name, **edge_attr) + + return graph +""" @_needs_graphviz def draw_flow_controllers(flow_controllers, graph=None, graph_attr=None, node_attr=None, diff --git a/interfaces/cython/cantera/reactor.pxd b/interfaces/cython/cantera/reactor.pxd index 8bb589e53d..6352fd39ad 100644 --- a/interfaces/cython/cantera/reactor.pxd +++ b/interfaces/cython/cantera/reactor.pxd @@ -27,6 +27,7 @@ cdef extern from "cantera/numerics/Integrator.h" namespace "Cantera": cdef extern from "cantera/zerodim.h" namespace "Cantera": cdef cppclass CxxWall "Cantera::Wall" cdef cppclass CxxReactorSurface "Cantera::ReactorSurface" + cdef cppclass CxxReactorEdge "Cantera::ReactorEdge" cdef cppclass CxxFlowDevice "Cantera::FlowDevice" # factories @@ -60,6 +61,7 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": CxxSparseMatrix jacobian() except +translate_exception CxxSparseMatrix finiteDifferenceJacobian() except +translate_exception void addSurface(CxxReactorSurface*) + void addEdge(CxxReactorEdge*) void setAdvanceLimit(string&, double) except +translate_exception void addSensitivityReaction(size_t) except +translate_exception void addSensitivitySpeciesEnthalpy(size_t) except +translate_exception @@ -127,6 +129,15 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": void addSensitivityReaction(size_t) except +translate_exception size_t nSensParams() + # reactor edge + + cdef cppclass CxxReactorEdge "Cantera::ReactorEdge": + CxxReactorEdge() + double length() + void setLength(double) + void setKinetics(CxxKinetics*) + void syncState() + # flow devices cdef cppclass CxxFlowDevice "Cantera::FlowDevice": @@ -223,6 +234,7 @@ cdef class ReactorBase: cdef list _outlets cdef list _walls cdef list _surfaces + cdef list _edges cdef public dict node_attr """ A dictionary containing draw attributes for the representation of the reactor as a @@ -282,6 +294,19 @@ cdef class ReactorSurface: .. versionadded:: 3.1 """ +cdef class ReactorEdge: + cdef CxxReactorEdge* edge + cdef Kinetics _kinetics + cdef ReactorBase _reactor + cdef public dict node_attr + """ + A dictionary containing draw attributes for the representation of the reactor + edge as a graphviz node. See https://graphviz.org/docs/nodes/ for a list of all + usable attributes. + + .. versionadded:: 3.1??? + """ + cdef class WallBase: cdef shared_ptr[CxxWallBase] _wall cdef CxxWallBase* wall diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index d2f383a590..059a60ad0a 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -41,6 +41,7 @@ cdef class ReactorBase: self._outlets = [] self._walls = [] self._surfaces = [] + self._edges = [] if isinstance(contents, _SolutionBase): self.insert(contents) # leave insert for the time being @@ -133,6 +134,11 @@ cdef class ReactorBase: def __get__(self): return self._surfaces + property edges: + """List of reacting edges installed on this reactor""" + def __get__(self): + return self._edges + def _add_inlet(self, inlet): """ Store a reference to ``inlet`` to prevent it from being prematurely @@ -928,6 +934,112 @@ cdef class ReactorSurface: """ self.surface.addSensitivityReaction(m) +cdef class ReactorEdge: + """ + Represents an edge in contact with the contents of a reactor. + + :param kin: + The `Kinetics` or `Interface` object representing reactions on this + surface. + :param r: + The `Reactor` into which this surface should be installed. + :param L: + The length of the reacting edge [m] + :param node_attr: + Attributes to be passed to the ``node`` method invoked to draw this surface. + See https://graphviz.org/docs/nodes/ for a list of all usable attributes. + + .. versionadded:: 3.1???? + Added the ``node_attr`` parameter. + """ + + def __cinit__(self): + self.edge = new CxxReactorEdge() + + def __dealloc__(self): + del self.edge + + def __init__(self, kin=None, Reactor r=None, *, L=None): + if kin is not None: + self.kinetics = kin + if r is not None: + self.install(r) + if L is not None: + self.length = L + #self.node_attr = node_attr or {'shape': 'underline'} + + def install(self, Reactor r): + """ + Add this `ReactorEdge` to the specified `Reactor` + """ + r._edges.append(self) + r.reactor.addEdge(self.edge) + self._reactor = r + + property length: + """ Length on which reactions can occur [m] """ + def __get__(self): + return self.edge.length() + def __set__(self, L): + self.edge.setLength(L) + + property kinetics: + """ + The `InterfaceKinetics` object used for calculating reaction rates on + this edge. + """ + def __get__(self): + self.edge.syncState() + return self._kinetics + def __set__(self, Kinetics k): + self._kinetics = k + self.edge.setKinetics(self._kinetics.kinetics) + + @property + def reactor(self): + """ + Return the `Reactor` object the edge is connected to. + + .. versionadded:: 3.1??? + """ + return self._reactor + +""" + def draw(self, graph=None, *, graph_attr=None, node_attr=None, surface_edge_attr=None, + print_state=False, **kwargs): + + Draw the edge as a ``graphviz`` ``dot`` node connected to its reactor. + The node is added to an existing ``graph`` if provided. + Optionally include current reactor state in the reactor node. + + :param graph: + ``graphviz.graphs.BaseGraph`` object to which the reactor is added. + If not provided, a new ``DiGraph`` is created. Defaults to ``None``. + :param graph_attr: + Attributes to be passed to the ``graphviz.Digraph`` function that control + the general appearance of the drawn network. + See https://graphviz.org/docs/graph/ for a list of all usable attributes. + :param node_attr: + Attributes to be passed to the ``node`` method invoked to draw the reactor. + See https://graphviz.org/docs/nodes/ for a list of all usable attributes. + :param surface_edge_attr: + Attributes to be passed to the ``edge`` method invoked to draw the + connection between the surface and its reactor. + See https://graphviz.org/docs/edges/ for a list of all usable attributes. + :param print_state: + Whether state information of the reactor is printed into its node. + Defaults to ``False`` + :param kwargs: + Additional keywords are passed on to ``~.drawnetwork.draw_reactor``. + :return: + ``graphviz.graphs.BaseGraph`` object with edge and connected + reactor. + + .. versionadded:: 3.1 + + return draw_edge(self, graph, graph_attr, node_attr, surface_edge_attr, + print_state, **kwargs) +""" cdef class WallBase: """ diff --git a/interfaces/sourcegen/sourcegen/csharp/config.yaml b/interfaces/sourcegen/sourcegen/csharp/config.yaml index be84c9c661..a2d7315b12 100644 --- a/interfaces/sourcegen/sourcegen/csharp/config.yaml +++ b/interfaces/sourcegen/sourcegen/csharp/config.yaml @@ -22,6 +22,7 @@ class_crosswalk: reactor: Reactor reactornet: ReactorNet reactorsurface: ReactorSurface + reactoredge: ReactorEdge soln: Solution surf: Surface thermo: ThermoPhase diff --git a/src/clib/ctreactor.cpp b/src/clib/ctreactor.cpp index 334c43f18e..330df05a05 100644 --- a/src/clib/ctreactor.cpp +++ b/src/clib/ctreactor.cpp @@ -25,12 +25,14 @@ typedef SharedCabinet ThermoCabinet; typedef SharedCabinet KineticsCabinet; typedef SharedCabinet SolutionCabinet; typedef SharedCabinet ReactorSurfaceCabinet; +typedef SharedCabinet ReactorEdgeCabinet; template<> ReactorCabinet* ReactorCabinet::s_storage = 0; template<> NetworkCabinet* NetworkCabinet::s_storage = 0; template<> FlowDeviceCabinet* FlowDeviceCabinet::s_storage = 0; template<> WallCabinet* WallCabinet::s_storage = 0; template<> ReactorSurfaceCabinet* ReactorSurfaceCabinet::s_storage = 0; +template<> ReactorEdgeCabinet* ReactorEdgeCabinet::s_storage = 0; template<> FuncCabinet* FuncCabinet::s_storage; // defined in ctfunc.cpp template<> ThermoCabinet* ThermoCabinet::s_storage; // defined in ct.cpp template<> KineticsCabinet* KineticsCabinet::s_storage; // defined in ct.cpp @@ -670,6 +672,66 @@ extern "C" { } } + // ReactorEdge + + int reactoredge_new(int type) + { + try { + return ReactorEdgeCabinet::add(make_shared()); + } catch (...) { + return handleAllExceptions(-1, ERR); + } + } + + int reactoredge_del(int i) + { + try { + ReactorEdgeCabinet::del(i); + return 0; + } catch (...) { + return handleAllExceptions(-1, ERR); + } + } + + int reactoredge_install(int i, int n) + { + try { + ReactorCabinet::item(n).addEdge(&ReactorEdgeCabinet::item(i)); + return 0; + } catch (...) { + return handleAllExceptions(-1, ERR); + } + } + + int reactoredge_setkinetics(int i, int n) + { + try { + ReactorEdgeCabinet::item(i).setKinetics(&KineticsCabinet::item(n)); + return 0; + } catch (...) { + return handleAllExceptions(-1, ERR); + } + } + + double reactoredge_length(int i) + { + try { + return ReactorEdgeCabinet::item(i).length(); + } catch (...) { + return handleAllExceptions(DERR, DERR); + } + } + + int reactoredge_setLength(int i, double v) + { + try { + ReactorEdgeCabinet::item(i).setLength(v); + return 0; + } catch (...) { + return handleAllExceptions(-1, ERR); + } + } + int ct_clearReactors() { try { @@ -678,6 +740,7 @@ extern "C" { FlowDeviceCabinet::clear(); WallCabinet::clear(); ReactorSurfaceCabinet::clear(); + ReactorEdgeCabinet::clear(); return 0; } catch (...) { return handleAllExceptions(-1, ERR); diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 0201222d63..7a83dde946 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -7,8 +7,10 @@ #include "cantera/zeroD/FlowDevice.h" #include "cantera/zeroD/Wall.h" #include "cantera/thermo/SurfPhase.h" +#include "cantera/thermo/EdgePhase.h" #include "cantera/zeroD/ReactorNet.h" #include "cantera/zeroD/ReactorSurface.h" +#include "cantera/zeroD/ReactorEdge.h" #include "cantera/kinetics/Kinetics.h" #include "cantera/kinetics/Reaction.h" #include "cantera/base/Solution.h" diff --git a/src/zeroD/ReactorBase.cpp b/src/zeroD/ReactorBase.cpp index 99b33688d6..f3d4a4f247 100644 --- a/src/zeroD/ReactorBase.cpp +++ b/src/zeroD/ReactorBase.cpp @@ -7,6 +7,7 @@ #include "cantera/zeroD/FlowDevice.h" #include "cantera/zeroD/ReactorNet.h" #include "cantera/zeroD/ReactorSurface.h" +#include "cantera/zeroD/ReactorEdge.h" #include "cantera/base/Solution.h" #include "cantera/thermo/ThermoPhase.h" @@ -124,11 +125,24 @@ void ReactorBase::addSurface(ReactorSurface* surf) } } +void ReactorBase::addEdge(ReactorEdge* edge) +{ + if (find(m_edges.begin(), m_edges.end(), edge) == m_edges.end()) { + m_edges.push_back(edge); + edge->setReactor(this); + } +} + ReactorSurface* ReactorBase::surface(size_t n) { return m_surfaces[n]; } +ReactorEdge* ReactorBase::edge(size_t n) +{ + return m_edges[n]; +} + void ReactorBase::restoreState() { if (!m_thermo) { throw CanteraError("ReactorBase::restoreState", "No phase defined."); diff --git a/src/zeroD/ReactorEdge.cpp b/src/zeroD/ReactorEdge.cpp new file mode 100644 index 0000000000..9e067ff601 --- /dev/null +++ b/src/zeroD/ReactorEdge.cpp @@ -0,0 +1,76 @@ +//! @file ReactorEdge.cpp + +// 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/ReactorEdge.h" +#include "cantera/zeroD/ReactorNet.h" +//#include "cantera/thermo/SurfPhase.h" +#include "cantera/thermo/EdgePhase.h" +#include "cantera/kinetics/Kinetics.h" +#include "cantera/kinetics/Reaction.h" + +namespace Cantera +{ + +double ReactorEdge::length() const +{ + return m_length; +} + +void ReactorEdge::setLength(double a) +{ + m_length = a; +} + +void ReactorEdge::setKinetics(Kinetics* kin) { + m_kinetics = kin; + if (kin == nullptr) { + m_thermo = nullptr; + return; + } + + m_thermo = dynamic_cast(&kin->thermo(0)); + if (m_thermo == nullptr) { + throw CanteraError("ReactorEdge::setKinetics", + "Specified kinetics manager does not represent an edge " + "kinetics mechanism."); + } + m_cov.resize(m_thermo->nSpecies()); + m_thermo->getCoverages(m_cov.data()); +} + +void ReactorEdge::setReactor(ReactorBase* reactor) +{ + m_reactor = reactor; +} + +void ReactorEdge::setCoverages(const double* cov) +{ + copy(cov, cov + m_cov.size(), m_cov.begin()); +} + +void ReactorEdge::setCoverages(const Composition& cov) +{ + m_thermo->setCoveragesByName(cov); + m_thermo->getCoverages(m_cov.data()); +} + +void ReactorEdge::setCoverages(const string& cov) +{ + m_thermo->setCoveragesByName(cov); + m_thermo->getCoverages(m_cov.data()); +} + +void ReactorEdge::getCoverages(double* cov) const +{ + copy(m_cov.begin(), m_cov.end(), cov); +} + +void ReactorEdge::syncState() +{ + m_thermo->setTemperature(m_reactor->temperature()); + m_thermo->setCoveragesNoNorm(m_cov.data()); +} + +}