Visualize reactor networks using `graphviz`

[Naikless]

Abstract

Include the possibility to easily visualize reactor networks by making use of the graphviz dot language.

Motivation

Having recently started to work with Cantera's zero dimensional reactor networks, I feel like there is much potential to make it easier for semi-advanced (not even to mention new) Cantera users to get into this topic.
While the creation of arbitrarily complex networks with many reactors that are connected by flow controllers and walls is possible, it is often very hard to verify that what was created is actually what was intended without a visual representation.

Such a visualization should primarily indicate connections between individual reactors. This would help to keep track of mass and energy balances. It could also indicate the state of each reactor, helping in assessing whether the obtained overall state for the network is plausible.

A visual representation of the reactor network is also one of the key elements of CHEMKIN's implementation of reactor networks. Thus including this in Cantera would strengthen the latter's abilities as an open source alternative. Here is an example I found googling.

Possible Solutions

The dot language of the open source graphviz graph visualization library allows to easily create schemes of reactor networks automatically. Each reactor is a node, each connection is an edge. The created dot file can be rendered as a vector or raster graphic or kept as a text for further manual editing.

A possible (but not flawless) implementation using the python interface of both graphviz and cantera could look like this:

from collections import defaultdict
import cantera as ct
import graphviz # obtainable by conda install -c conda-forge python-graphviz

def plot_CRN(net, **kwargs):
    
    dot = graphviz.Digraph(**kwargs)
    
    # create nodes that also display T and P for all reactors
    # find all flow controllers and walls referenced by any reactor in the net
    flow_controllers = []
    walls = []
    for reactor in net:
        dot.node(reactor.name, shape='box',
                 label=f'{reactor.name}\n\nT: {reactor.T:.2f} K\nP: {reactor.thermo.P*1e-5:.2f} bar')     
        flow_controllers.extend(reactor.inlets)
        flow_controllers.extend(reactor.outlets)
        walls.extend(reactor.walls)
    
    # filter out all flow elements that connect two reactors in the net
    connections = {}
    for fc in set(flow_controllers):
        inlet_to = [reactor for reactor in net if fc in reactor.inlets]
        outlet_to = [reactor for reactor in net if fc in reactor.outlets]
        if inlet_to and outlet_to:
            connections[fc] = (outlet_to[0], inlet_to[0])
    
    # first, group all connections that connect the same two reactors
    edges_mf = defaultdict(list)
    for key,value in connections.items():
        edges_mf[value].append(key)
    
    # now sum up the mass flows between each connected reactor pair
    for edge, flow_controllers in edges_mf.items():
        mass_flow = sum(fc.mass_flow_rate for fc in flow_controllers)
        dot.edge(edge[0].name, edge[1].name, label=f'm = {mass_flow*3600:.2f} kg/h')
    
    # same for heat fluxes through walls. However, direction can not be determined!
    heat_flows = {}
    for wall in set(walls):
        reactors = [reactor for reactor in net if wall in reactor.walls]
        if len(reactors) == 2:
            heat_flows[wall] = (reactors[0], reactors[1])
    
    edges_hf = defaultdict(list)
    for key, value in heat_flows.items():
        edges_hf[value].append(key)
    
    for edge, heat_flows  in edges_hf.items():
        heat_flow = sum([hf.qdot(0) for hf in heat_flows])
        dot.edge(edge[0].name, edge[1].name, color='red', style='dashed', label=f'q = {heat_flow*1e-3:.2f} kW')
        
    return dot

    
gas = ct.Solution('gri30.yaml')

inlet = ct.Reservoir(gas,'inlet')
outlet = ct.Reservoir(gas, 'outlet')
r1 = ct.Reactor(gas, name='R1')
r2 = ct.Reactor(gas, name='R2')

mc = ct.MassFlowController(inlet, r1)
pc = ct.PressureController(r1, r2, master=mc)
pc_out = ct.PressureController(r2, outlet, master=mc)

w = ct.Wall(r1,r2,Q=1e3)

net = [r1,r2]

sim = ct.ReactorNet(net)

sim.advance_to_steady_state()

dot = plot_CRN([*net,inlet,outlet], graph_attr={'rankdir':'LR', 'nodesep':'1'})
dot

One caveat is that the current python objects are very sparse regarding the saving of network relations:

• ReactorNet objects do not provide a list of all reactors (although their C++ counterpart certainly has it and there is even a self._reactors attribute created at initialization. Can someone enlighten me why it disappears in Python later?). For this reason, the above implementation needs one to supply a list of all involved reactors instead.
• Similarly, the only way to find the connecting elements to a flow element is by looping through the reactors and reservoirs. If the flow elements would provide these information themselves (i.e. allow access to self._upstream and self._downstream), things would become a lot easier.
• For the same reasons, wall connections also can only be detected by finding reactors connected to the same wall. Unfortunately, there is no way to determine the direction of the heat flow this way.

It would also be useful to allow for dedicated dot attributes of reactors, reservoirs, flow controllers and walls to allow to specify things like color, shape and arrow style of the visualization.

References

graphviz reference
graphviz Python interface

[ischoegl]

@Naikless ... thanks for the interesting enhancement request! I'm linking (somewhat) related issues: there was some discussion about NetworkX etc. rather than graphviz in #173, and there is a stub about a yaml summary of reactor networks in Cantera/cantera#694 that could potentially be leveraged.

[Naikless]

Thanks for the feedback!

I'll check out NetworkX to see whether it might be more suitable.

I would also be happy to extend my initial approach to an actual PR, but I believe the following questions would have to be answered before:

1. Does it even make sense to implement this at the Python interface level, or should this be done on the C++ level (which I am not very familiar with)?

2. If using the Python level, what would need to be changed to include additional information in the Python objects, such as the list of reactors in a network object or the attached reactors to a flow controller or wall. As mentioned above, from looking at the source, it seems like the respective attributes (self._reactors, self._upstream...) should already be there and I don't understand why they are not.

Edit: looking more into the cython source, it seems the answer to 2. is that the respective attributes that are currently only declared in the .pxd file would have to be included as properties in the pyx file as well. Is this correct? Any reasons not to implement this?

[bryanwweber]

Does it even make sense to implement this at the Python interface level, or should this be done on the C++ level (which I am not very familiar with)?

I think C++ would be the preferred long-term home for a feature like this so it could be used from any interface. However, something is usually better than nothing, so if you're able to implement this in Python as a first pass, I think that'd still be super helpful!

[Naikless]

Alright. I'll try to come up with something that is mergeable then, so that it can be further reviewed.

I also checked out NetworkX. From my perspective, it primarily does the job of defining the graph in terms of nodes and edges. However, the way I see it, this information is already present in the defined reactor net. Thus, implementing a NetworkX representation parallel to that feels redundant. For visualization, NetworkX also invokes graphviz, so I don't think much is gained by this additional abstraction level.

I would then go ahead and try to make the above mentioned currently hidden network attributes accessible from the Python site, if nothing speaks against this.