diff --git a/docs/_quarto.yml b/docs/_quarto.yml index 0a69bdbe5..ec52e0725 100644 --- a/docs/_quarto.yml +++ b/docs/_quarto.yml @@ -43,6 +43,9 @@ website: - guide/index.qmd - guide/examples.ipynb - guide/coupling.qmd + - section: "Coupling guides" + contents: + - guide/delwaq.ipynb - title: "Concepts" contents: diff --git a/docs/guide/coupling.qmd b/docs/guide/coupling.qmd index 95a5019b9..c1ac43ecb 100644 --- a/docs/guide/coupling.qmd +++ b/docs/guide/coupling.qmd @@ -10,7 +10,7 @@ Ribasim can also be (online) coupled to other kernels with the help of iMOD Coup ## Water quality -Ribasim can be offline coupled to Delwaq, the Deltares Water Quality model. Note that this functionality is still in active development. +Ribasim can be offline coupled to Delwaq, the Deltares Water Quality model. ```{mermaid} flowchart LR @@ -21,5 +21,4 @@ flowchart LR Delwaq can calculate the concentration of substances in [Basin](/reference/node/basin.qmd) nodes over time, based on initial concentrations, and of [FlowBoundary](/reference/node/flow-boundary.qmd) nodes. Ribasim exposes the `Basin / concentration`, `Basin / concentration_state`, `FlowBoundary / concentration`, and `LevelBoundary / concentration` tables to setup these substances and concentrations. -When a Ribasim model ran with the above tables, one can use the utilities in the `coupling/delwaq` folder to generate the input required for Delwaq to run, as well as to parse the output from Delwaq into a Ribasim compatible format. -For more information see the `README.md` in the same folder. +When a Ribasim model ran with the above tables, one can use the utilities in the `delwaq` namespace of the Ribasim Python API to generate the input required for Delwaq to run, as well as to parse the output from Delwaq into a Ribasim compatible format. For more information see the [guide](/guide/delwaq.ipynb). diff --git a/docs/guide/delwaq.ipynb b/docs/guide/delwaq.ipynb new file mode 100644 index 000000000..e81472512 --- /dev/null +++ b/docs/guide/delwaq.ipynb @@ -0,0 +1,208 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "title: \"Ribasim Delwaq coupling\"\n", + "---" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In order to generate the Delwaq input files, we need a completed Ribasim simulation (typically one with a results folder) that ideally also includes some substances and initial concentrations. Let's take the basic test model for example, which already has set some initial concentrations.\n", + "\n", + "All testmodels can be [downloaded from here](/install.qmd#sec-download)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "\n", + "toml_path = Path(\"../../generated_testmodels/basic/ribasim.toml\")\n", + "assert toml_path.is_file()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This Ribasim model already has substance concentrations for `Cl` and `Tracer` in the input tables, and we will use these to generate the Delwaq input files." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from ribasim import Model\n", + "\n", + "model = Model.read(toml_path)\n", + "\n", + "display(model.basin.concentration_state) # basin initial state\n", + "display(model.basin.concentration) # basin boundaries\n", + "display(model.flow_boundary.concentration) # flow boundaries\n", + "display(model.level_boundary.concentration) # level boundaries\n", + "model.plot(); # for later comparison" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# | include: false\n", + "from subprocess import run\n", + "\n", + "run(\n", + " [\n", + " \"julia\",\n", + " \"--project=../../core\",\n", + " \"--eval\",\n", + " f'using Ribasim; Ribasim.main(\"{toml_path.as_posix()}\")',\n", + " ],\n", + " check=True,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Given the path to a completed Ribasim simulation, we can call `ribasim.delwaq.generate` for generating the required input files for Delwaq from scratch." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from ribasim.delwaq import generate\n", + "\n", + "output_path = Path(\n", + " \"../../generated_testmodels/basic/delwaq\"\n", + ") # set a path where we store the Delwaq input files\n", + "graph, substances = generate(toml_path, output_path)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This call produces a handful of files in the user defined folder. Let's take a look at them:\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "list(output_path.iterdir())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "These files form a complete Delwaq simulation, and can be run by either pointing DIMR to the `dimr_config.xml` file or pointing Delwaq to the `delwaq.inp` file." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that the call to `generate` produces two output variables; `graph` and `substances` that are required for parsing the results of the Delwaq model later on. Nonetheless, we can also inspect them here, and inspect the created Delwaq network." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "substances # list of substances, as will be present in the Delwaq netcdf output" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As you can see, the complete substances list is a combination of user input (`Cl` and `Tracer` in the input tables), a `Continuity` tracer, and tracers for all nodetypes in the Ribasim model. The latter tracers allow for deeper inspection of the Ribasim model, such as debugging the mass balance by plotting fraction graphs. Let's inspect the `graph` next, which is the Delwaq network that was created from the Ribasim model:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import networkx as nx\n", + "\n", + "# Let's draw the graph\n", + "fig, ax = plt.subplots(1, 2, figsize=(10, 5))\n", + "nx.draw(\n", + " graph,\n", + " pos={k: v[\"pos\"] for k, v in graph.nodes(data=True)},\n", + " with_labels=True,\n", + " labels={k: k for k, v in graph.nodes(data=True)},\n", + " ax=ax[0],\n", + ")\n", + "ax[0].set_title(\"Delwaq node IDs\")\n", + "nx.draw(\n", + " graph,\n", + " pos={k: v[\"pos\"] for k, v in graph.nodes(data=True)},\n", + " with_labels=True,\n", + " labels={k: v[\"id\"] for k, v in graph.nodes(data=True)},\n", + " ax=ax[1],\n", + ")\n", + "ax[1].set_title(\"Ribasim node IDs\")\n", + "fig.suptitle(\"Delwaq network\");" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we plotted the Delwaq network twice, with the node IDs as used by Delwaq on the left hand side, and the corresponding Ribasim node IDs on the right hand side.\n", + "As you can see, the Delwaq network is very similar to the Ribasim network, with some notable changes:\n", + "\n", + "- All non-Basin or non-boundary types are removed (e.g. no more Pumps or TabulatedRatingCurves)\n", + "- Basin boundaries are split into separate nodes and links (drainage, precipitation, and evaporation, as indicated by the duplicated Basin IDs on the right hand side)\n", + "- All node IDs have been renumbered, with boundaries being negative, and Basins being positive." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "default", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.4" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/python/ribasim/ribasim/delwaq/generate.py b/python/ribasim/ribasim/delwaq/generate.py index fdfbd0c6e..e3ae67f92 100644 --- a/python/ribasim/ribasim/delwaq/generate.py +++ b/python/ribasim/ribasim/delwaq/generate.py @@ -31,6 +31,7 @@ ) delwaq_dir = Path(__file__).parent +output_folder = delwaq_dir / "model" env = jinja2.Environment( autoescape=True, loader=jinja2.FileSystemLoader(delwaq_dir / "template") @@ -41,20 +42,53 @@ USE_EVAP = True -def generate(toml_path: Path) -> tuple[nx.DiGraph, set[str]]: - """Generate a Delwaq model from a Ribasim model and results.""" - - # Read in model and results - model = ribasim.Model.read(toml_path) - basins = pd.read_feather(toml_path.parent / "results" / "basin.arrow") - flows = pd.read_feather(toml_path.parent / "results" / "flow.arrow") +def _boundary_name(id, type): + # Delwaq has a limit of 12 characters for the boundary name + return type[:9] + "_" + str(id) + + +def _quote(value): + return f"'{value}'" + + +def _make_boundary(data, boundary_type): + """ + Create a Delwaq boundary definition with the given data and boundary type. + Pivot our data from long to wide format, and convert the time to a string. + + Specifically, we go from a table: + `node_id, substance, time, concentration` + to + ``` + ITEM 'Drainage_6' + CONCENTRATIONS 'Cl' 'Tracer' + ABSOLUTE TIME + LINEAR DATA 'Cl' 'Tracer' + '2020/01/01-00:00:00' 0.0 1.0 + '2020/01/02-00:00:00' 1.0 -999 + ``` + """ + bid = _boundary_name(data.node_id.iloc[0], boundary_type) + piv = ( + data.pivot_table(index="time", columns="substance", values="concentration") + .reset_index() + .reset_index(drop=True) + ) + piv.time = piv.time.dt.strftime("%Y/%m/%d-%H:%M:%S") + boundary = { + "name": bid, + "substances": list(map(_quote, piv.columns[1:])), + "df": piv.to_string( + formatters={"time": _quote}, header=False, index=False, na_rep=-999 + ), + } + substances = data.substance.unique() + return boundary, substances - output_folder = delwaq_dir / "model" - output_folder.mkdir(exist_ok=True) - # Setup flow network +def _setup_graph(nodes, edge, use_evaporation=True): G = nx.DiGraph() - nodes = model.node_table() + assert nodes.df is not None for row in nodes.df.itertuples(): if row.node_type not in ribasim.geometry.edge.SPATIALCONTROLNODETYPES: @@ -66,8 +100,8 @@ def generate(toml_path: Path) -> tuple[nx.DiGraph, set[str]]: y=row.geometry.y, pos=(row.geometry.x, row.geometry.y), ) - assert model.edge.df is not None - for row in model.edge.df.itertuples(): + assert edge.df is not None + for row in edge.df.itertuples(): if row.edge_type == "flow": G.add_edge( f"{row.from_node_type} #{row.from_node_id}", @@ -184,7 +218,7 @@ def generate(toml_path: Path) -> tuple[nx.DiGraph, set[str]]: boundary=(node["id"], "precipitation"), ) - if USE_EVAP: + if use_evaporation: boundary_id -= 1 G.add_node( boundary_id, @@ -206,6 +240,57 @@ def generate(toml_path: Path) -> tuple[nx.DiGraph, set[str]]: for edge_id in d["id"]: edge_mapping[edge_id] = i + assert len(basin_mapping) == basin_id + + return G, merge_edges, node_mapping, edge_mapping, basin_mapping + + +def _setup_boundaries(model): + boundaries = [] + substances = set() + + if model.level_boundary.concentration.df is not None: + for _, rows in model.level_boundary.concentration.df.groupby(["node_id"]): + boundary, substance = _make_boundary(rows, "LevelBoundary") + boundaries.append(boundary) + substances.update(substance) + + if model.flow_boundary.concentration.df is not None: + for _, rows in model.flow_boundary.concentration.df.groupby("node_id"): + boundary, substance = _make_boundary(rows, "FlowBoundary") + boundaries.append(boundary) + substances.update(substance) + + if model.basin.concentration.df is not None: + for _, rows in model.basin.concentration.df.groupby(["node_id"]): + for boundary_type in ("Drainage", "Precipitation"): + nrows = rows.rename(columns={boundary_type.lower(): "concentration"}) + boundary, substance = _make_boundary(nrows, boundary_type) + boundaries.append(boundary) + substances.update(substance) + + return boundaries, substances + + +def generate( + toml_path: Path, + output_folder=output_folder, + use_evaporation=USE_EVAP, +) -> tuple[nx.DiGraph, set[str]]: + """Generate a Delwaq model from a Ribasim model and results.""" + + # Read in model and results + model = ribasim.Model.read(toml_path) + basins = pd.read_feather(toml_path.parent / "results" / "basin.arrow") + flows = pd.read_feather(toml_path.parent / "results" / "flow.arrow") + + output_folder.mkdir(exist_ok=True) + + # Setup flow network + G, merge_edges, node_mapping, edge_mapping, basin_mapping = _setup_graph( + model.node_table(), model.edge, use_evaporation=use_evaporation + ) + # Plot # plt.figure(figsize=(18, 18)) # nx.draw( @@ -227,7 +312,7 @@ def generate(toml_path: Path) -> tuple[nx.DiGraph, set[str]]: pointer.to_csv(output_folder / "network.csv", index=False) # not needed write_pointer(output_folder / "ribasim.poi", pointer) - total_segments = basin_id + total_segments = len(basin_mapping) total_exchanges = len(pointer) # Write attributes template @@ -309,69 +394,7 @@ def generate(toml_path: Path) -> tuple[nx.DiGraph, set[str]]: write_flows(output_folder / "ribasim.len", lengths, timestep) # Find all boundary substances and concentrations - boundaries = [] - substances = set() - - def boundary_name(id, type): - # Delwaq has a limit of 12 characters for the boundary name - return type[:9] + "_" + str(id) - - def quote(value): - return f"'{value}'" - - def make_boundary(data, boundary_type): - """ - Create a Delwaq boundary definition with the given data and boundary type. - Pivot our data from long to wide format, and convert the time to a string. - - Specifically, we go from a table: - `node_id, substance, time, concentration` - to - ``` - ITEM 'Drainage_6' - CONCENTRATIONS 'Cl' 'Tracer' - ABSOLUTE TIME - LINEAR DATA 'Cl' 'Tracer' - '2020/01/01-00:00:00' 0.0 1.0 - '2020/01/02-00:00:00' 1.0 -999 - ``` - """ - bid = boundary_name(data.node_id.iloc[0], boundary_type) - piv = ( - data.pivot_table(index="time", columns="substance", values="concentration") - .reset_index() - .reset_index(drop=True) - ) - piv.time = piv.time.dt.strftime("%Y/%m/%d-%H:%M:%S") - boundary = { - "name": bid, - "substances": list(map(quote, piv.columns[1:])), - "df": piv.to_string( - formatters={"time": quote}, header=False, index=False, na_rep=-999 - ), - } - substances = data.substance.unique() - return boundary, substances - - if model.level_boundary.concentration.df is not None: - for _, rows in model.level_boundary.concentration.df.groupby(["node_id"]): - boundary, substance = make_boundary(rows, "LevelBoundary") - boundaries.append(boundary) - substances.update(substance) - - if model.flow_boundary.concentration.df is not None: - for _, rows in model.flow_boundary.concentration.df.groupby("node_id"): - boundary, substance = make_boundary(rows, "FlowBoundary") - boundaries.append(boundary) - substances.update(substance) - - if model.basin.concentration.df is not None: - for _, rows in model.basin.concentration.df.groupby(["node_id"]): - for boundary_type in ("Drainage", "Precipitation"): - nrows = rows.rename(columns={boundary_type.lower(): "concentration"}) - boundary, substance = make_boundary(nrows, boundary_type) - boundaries.append(boundary) - substances.update(substance) + boundaries, substances = _setup_boundaries(model) # Write boundary data with substances and concentrations template = env.get_template("B5_bounddata.inc.j2") @@ -427,7 +450,7 @@ def make_boundary(data, boundary_type): bnd.sort_values(by="bid", ascending=False, inplace=True) bnd["node_type"] = [G.nodes(data="type")[bid] for bid in bnd["bid"]] bnd["node_id"] = [G.nodes(data="id")[bid] for bid in bnd["bid"]] - bnd["fid"] = list(map(boundary_name, bnd["node_id"], bnd["node_type"])) + bnd["fid"] = list(map(_boundary_name, bnd["node_id"], bnd["node_type"])) bnd["comment"] = "" bnd = bnd[["fid", "comment", "node_type"]] bnd.to_csv( diff --git a/python/ribasim/ribasim/delwaq/parse.py b/python/ribasim/ribasim/delwaq/parse.py index 7c0001e8c..3c527d02a 100644 --- a/python/ribasim/ribasim/delwaq/parse.py +++ b/python/ribasim/ribasim/delwaq/parse.py @@ -17,7 +17,9 @@ output_folder = delwaq_dir / "model" -def parse(toml_path: Path, graph, substances) -> ribasim.Model: +def parse( + toml_path: Path, graph, substances, output_folder=output_folder +) -> ribasim.Model: model = ribasim.Model.read(toml_path) # Output of Delwaq