Offline coupling with Delwaq

[evetion]

Fixes #649

This PR adds a folder coupling/delwaq for offline coupling with Delwaq, with a more extensive description provided in README.md in the same folder.

A test is provided in test.py, that is also run on Teamcity on Windows, using the weekly DIMR releases. The basic testmodel (the only with built-in substance concentration values) is used for this. In order to speed-up the CI one can now do pixi run ribasim-core-testmodels basic to only run the basic model (but any or multiple model names can be used).

[evetion]

So this is up to date with main again (I think), and not working, given the non consecutive boundaries ids. Next up is to remove non-basin/boundary nodes, it was nice to copy the Ribasim network 1-1 to Delwaq, but it seems to introduce some error in the results. The continituity tracer gives the worst mass balance at non-basin nodes (expected, as we set some non-zero storage), so let's remove them.

However, pruning the network is not that straightforward on a pure topology basis. While it's easy to remove a single node between two basins, there are some nodes that require recalculating flows based on flow directions/rules. See our basic model for example, where I crossed the easy ones, but removing the circled tabulated rating curve #4 requires more work.

I would propose to introduce a dependency on networkx to handle the graph related logic to prune the network. This is also a good intermediate issue that can be tackled by others.

[SouthEndMusic]

Ah, you're starting to see why FractionalFlow is my nemesis.

If I understand correctly, you want a routine that takes in the Ribasim model and spits out a new graph with edges only between basins and boundary nodes. You also want to know the flows over these edges. This might work (looking only at flow edges of course):

1. Start with the original graph and make each edge its own representative edge;
2. Loop over all nodes in the network. If it is a connector node, remove it and its connected edges, and create new edges from the upstream node to all downstream nodes. The representative edge of the new edges is the edge between the connector node and the downstream node;
3. Now you should be left with the graph you want. It can happen that there are multiple edges between the same pair of nodes (so the procedure above should allow this), so in that case the flow is given by the sum of the flows over their representative edges.

Using this representative edge system, I don't think any additional calculations or checks are needed.

[evetion]

Ha, yes I understand your struggles.

What you propose is indeed what I want, with the addition that the nodes need to be renumbered so they're consecutive. Same for the boundary nodes, but then negatively numbered.

Is representative edge a name you came up with? You mean as an object to reason about (like done in the networkx package)?

This approach requires directed edges, and I wonder whether it always holds. I can imagine that if the fractionalflows here would be nodes that allow bidirectional (negative) flow, a solution would be morw complex. Not sure if such configurations are possible with our current nodetypes.

[SouthEndMusic]

Yes, representative edge is a term I came up with, I haven't looked into the networkx package yet. What I mean is that the edges you have in your delwaq graph are generally not in the original graph, so the flow over the edges in the delwaq graph is given by the sum of the flows over the representative edges in the original graph.

I'm not sure why the solution would be more complex with bidirectional edges. Do you want to convert to edges with only positive flow, so that you need edges in both directions?

I wrote this function which I think does what you want:

function create_graph_delwaq(graph::MetaGraph)::MetaGraph
    graph_delwaq = MetaGraph(
        DiGraph();
        label_type = NodeID,
        vertex_data_type = Int32,
        # Representative edges
        edge_data_type = Vector{Tuple{NodeID, NodeID}},
        graph_data = nothing,
    )

    # Add all flow edges to the delwaq graph
    for edge_metadata in values(graph.edge_data)
        (; type, from_id, to_id) = edge_metadata
        if type == EdgeType.flow

            # These zeros are the delwaq IDs to be set later
            graph_delwaq[from_id] = 0
            graph_delwaq[to_id] = 0

            edge = (from_id, to_id)
            graph_delwaq[edge...] = [edge]
        end
    end

    # Remove connector nodes from the delwaq graph
    for node_id in collect(values(graph_delwaq.vertex_labels))
        if node_id.type ∉ [
            NodeType.Basin,
            NodeType.Terminal,
            NodeType.LevelBoundary,
            NodeType.FlowBoundary,
            NodeType.UserDemand,
        ]
            inneighbor_id = only(inneighbor_labels(graph_delwaq, node_id))
            outneighbor_ids = collect(outneighbor_labels(graph_delwaq, node_id))

            # From MetagraphsNext
            rem_vertex!(graph_delwaq, code_for(graph_delwaq, node_id))

            for outneighbor_id in outneighbor_ids
                edge = (inneighbor_id, outneighbor_id)
                representative_edge = (node_id, outneighbor_id)
                # If the new edge already exists, add a representative edge,
                # otherwise create the new edge
                if haskey(graph_delwaq, edge...)
                    push!(graph_delwaq[edge...], representative_edge)
                else
                    graph_delwaq[edge...] = [representative_edge]
                end
            end
        end
    end

    # Remap IDs
    basin_id = 0
    boundary_id = 0

    for node_id in values(graph_delwaq.vertex_labels)
        if node_id.type == NodeType.Basin
            basin_id += 1
            graph_delwaq[node_id] = basin_id
        elseif node_id.type in [
            NodeType.Terminal,
            NodeType.UserDemand,
            NodeType.LevelBoundary,
            NodeType.FlowBoundary,
        ]
            boundary_id -= 1
            graph_delwaq[node_id] = boundary_id
        else
            error("Found unexpected node $node_id in delwaq graph.")
        end
    end

    return graph_delwaq
end

For the basic model this gives the following representative edges:

julia> graph_delwaq.edge_data
Dict{Tuple{Ribasim.NodeID, Ribasim.NodeID}, Vector{Tuple{Ribasim.NodeID, Ribasim.NodeID}}} with 9 entries:
  (Basin #3, Basin #9)          => [(FractionalFlow #8, Basin #9)]
  (Basin #3, Basin #6)          => [(FractionalFlow #5, Basin #6)]
  (Basin #9, LevelBoundary #17) => [(LinearResistance #10, LevelBoundary #17)]
  (Basin #6, Basin #9)          => [(Pump #7, Basin #9)]
  (Basin #1, Basin #3)          => [(ManningResistance #2, Basin #3)]
  (LevelBoundary #11, Basin #3) => [(LinearResistance #12, Basin #3)]
  (FlowBoundary #16, Basin #1)  => [(FlowBoundary #16, Basin #1)]
  (Basin #3, Terminal #14)      => [(FractionalFlow #13, Terminal #14)]
  (FlowBoundary #15, Basin #6)  => [(FlowBoundary #15, Basin #6)]

and for a simple example where there are multiple paths between 2 basins, testmodel tabulated_rating_curve:

julia> graph_delwaq.edge_data
Dict{Tuple{Ribasim.NodeID, Ribasim.NodeID}, Vector{Tuple{Ribasim.NodeID, Ribasim.NodeID}}} with 1 entry:
  (Basin #1, Basin #4) => [(TabulatedRatingCurve #2, Basin #4), (TabulatedRatingCurve #3, Basin #4)]

Edit:

Now you can also get the delwaq_id with graph_delwaq[node_id], or view all via graph_delwaq.vertex_properties. For example the basic model:

julia> graph_delwaq.vertex_properties
Dict{Ribasim.NodeID, Tuple{Int64, Int32}} with 9 entries:
  Terminal #14      => (8, -2)
  LevelBoundary #17 => (6, -3)
  FlowBoundary #16  => (3, -4)
  Basin #9          => (7, 3)
  Basin #6          => (9, 2)
  Basin #1          => (1, 1)
  Basin #3          => (4, 4)
  FlowBoundary #15  => (5, -1)
  LevelBoundary #11 => (2, -5)

Don't mind the first integer in each tuple, that is the code created internally in MetaGraphsNext for each node.

[visr]

Left a bunch of comments, mostly small. If they are addressed to the degree that they are applicable, it is good to merge from my side.

[visr]

Quality work!

[evetion]

Thanks, took me way longer that I would've liked, but very happy now that I've got it to work, including the CI.