From 708bebadf471da1ac22dfd0e4eeb9113a300182c Mon Sep 17 00:00:00 2001 From: Joseph Hughes Date: Wed, 5 Aug 2026 12:40:21 -0500 Subject: [PATCH 1/3] docs(lakpak_utils): add embedded lake connection examples get_lak_connections builds the Lake Package connection data for a lake embedded in the model domain, and support for a vertex grid was added in #2786, but no notebook showed either one. Two notebooks are added. The first embeds a lake in a structured grid, and the second embeds the same lake in the equivalent vertex grid and checks that the two give the same connections. Both run the model and put the lake stage into the head array so that a single map shows the water surface across the lake and the aquifer rather than leaving a hole where the lake is. --- .../Notebooks/dis_lake_connections_example.py | 190 +++++++++++++++ .../disv_lake_connections_example.py | 224 ++++++++++++++++++ 2 files changed, 414 insertions(+) create mode 100644 .docs/Notebooks/dis_lake_connections_example.py create mode 100644 .docs/Notebooks/disv_lake_connections_example.py diff --git a/.docs/Notebooks/dis_lake_connections_example.py b/.docs/Notebooks/dis_lake_connections_example.py new file mode 100644 index 000000000..51949613e --- /dev/null +++ b/.docs/Notebooks/dis_lake_connections_example.py @@ -0,0 +1,190 @@ +# --- +# jupyter: +# jupytext: +# notebook_metadata_filter: all +# text_representation: +# extension: .py +# format_name: light +# format_version: '1.5' +# jupytext_version: 1.14.5 +# kernelspec: +# display_name: Python 3 (ipykernel) +# language: python +# name: python3 +# metadata: +# section: mf6 +# authors: +# - name: Joseph Hughes +# --- + +# # Embedded Lake Connections for a Structured Grid +# +# A lake embedded in the model domain replaces the cells it occupies, so the Lake (LAK) Package needs a connection to every active cell that touches it. `get_lak_connections()` builds those connections from an array of lake numbers, and returns the idomain with the lake cells deactivated, the number of connections in each lake, and the connectiondata block for the package. +# +# This notebook builds a lake embedded in one layer of a structured grid. The companion notebook does the same thing on a vertex grid. + +# + +import sys +from pathlib import Path +from tempfile import TemporaryDirectory + +import matplotlib as mpl +import matplotlib.pyplot as plt +import numpy as np + +import flopy +from flopy.discretization import StructuredGrid +from flopy.mf6.utils import get_lak_connections + +print(sys.version) +print(f"numpy version: {np.__version__}") +print(f"matplotlib version: {mpl.__version__}") +print(f"flopy version: {flopy.__version__}") +# - + +# ## Build the grid and the lake +# +# The grid is two layers of 7 rows and 7 columns. The lake occupies a three by three block in the upper layer, and is given lake number 0. Cells that are not part of a lake are masked. + +# + +nlay, nrow, ncol = 2, 7, 7 +delr = delc = np.full(7, 100.0) +top = np.full((nrow, ncol), 10.0) +botm = np.array([np.full((nrow, ncol), 0.0), np.full((nrow, ncol), -10.0)]) + +modelgrid = StructuredGrid(delr=delr, delc=delc, top=top, botm=botm, nlay=nlay) + +lake_map = np.full((nlay, nrow, ncol), -1, dtype=np.int32) +lake_map[0, 2:5, 2:5] = 0 +lake_map = np.ma.masked_where(lake_map < 0, lake_map) + +print(f"Lake cells in layer 1: {(~lake_map.mask[0]).sum()}") +# - + +# ## Build the lake connections +# +# The lake map is the same shape as the model grid, so the lake is embedded rather than sitting on top of the model. A bed leakance is given for every connection. + +# + +idomain, connection_dict, connectiondata = get_lak_connections( + modelgrid, + lake_map, + idomain=np.ones((nlay, nrow, ncol), dtype=int), + bedleak=0.1, +) + +print(f"Connections in each lake: {connection_dict}") +print(f"connectiondata rows: {len(connectiondata)}") +# - + +# A connection is horizontal where the lake meets a cell in the same layer, and vertical where it sits on the cell below. The three by three lake has twelve cells around its perimeter and nine beneath it. + +# + +claktype = [row[3] for row in connectiondata] +for kind in ("horizontal", "vertical"): + print(f"{kind:12s}{claktype.count(kind)}") + +print() +print("lakeno iconn cellid claktype bedleak connlen connwidth") +for row in connectiondata[:4]: + lakeno, iconn, cellid, kind, leak, _, _, connlen, connwidth = row + print( + f"{lakeno:6d}{iconn:6d} {cellid!s:15s} {kind:12s}" + f"{leak:8.2f}{connlen:9.1f}{connwidth:10.1f}" + ) +# - + +# The cells the lake occupies are deactivated in the returned idomain, so the lake replaces them rather than sharing the domain with them. + +print(f"idomain in the lake cells: {np.unique(idomain[0, 2:5, 2:5])}") +print(f"active cells remaining : {(idomain > 0).sum()} of {idomain.size}") + +# ## Run the model +# +# The connectiondata and the idomain go straight into the model. Constant heads +# on the left and right edges drive flow across the lake. + +# + +temp_dir = TemporaryDirectory() +name = "dis_lake" + +sim = flopy.mf6.MFSimulation(sim_name=name, sim_ws=temp_dir.name, exe_name="mf6") +flopy.mf6.ModflowTdis(sim) +flopy.mf6.ModflowIms( + sim, + linear_acceleration="bicgstab", + outer_dvclose=1e-9, + inner_dvclose=1e-10, + outer_maximum=200, +) +gwf = flopy.mf6.ModflowGwf(sim, modelname=name, newtonoptions="newton under_relaxation") +flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + idomain=idomain, +) +flopy.mf6.ModflowGwfic(gwf, strt=8.0) +flopy.mf6.ModflowGwfnpf(gwf, icelltype=1, k=10.0) +flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=[[(0, i, 0), 9.0] for i in range(nrow)] + + [[(0, i, ncol - 1), 6.0] for i in range(nrow)], +) +lak = flopy.mf6.ModflowGwflak( + gwf, + stage_filerecord=f"{name}.lak.stage.bin", + nlakes=1, + packagedata=[[0, 7.5, connection_dict[0]]], + connectiondata=connectiondata, + perioddata={0: [[0, "rainfall", 0.001]]}, +) +flopy.mf6.ModflowGwfoc(gwf, head_filerecord=f"{name}.hds", saverecord=[("HEAD", "ALL")]) + +sim.write_simulation(silent=True) +success, buff = sim.run_simulation(silent=True) +assert success, "\n".join(buff[-15:]) +print("model converged") +# - + +# ## Put the lake stage into the head array +# +# The lake cells are inactive, so the head array has no value there. Filling +# those cells with the stage of the lake that occupies them gives a single +# array of the water surface that can be plotted in one pass. + +# + +head = gwf.output.head().get_data() +stage = lak.output.stage().get_data().flatten() +for lake_number in np.unique(lake_map.compressed()): + head[lake_map == lake_number] = stage[lake_number] + +print(f"lake stage: {stage[0]:.3f}") +print(f"cells with no head value: {(np.abs(head) > 1e29).sum()}") +print(f"head range: {head.min():.3f} to {head.max():.3f}") +# - + +# ## Head map +# +# The lake reads as part of the water surface rather than as a hole in it. + +# + +fig, ax = plt.subplots(figsize=(6, 6)) +ax.set_aspect("equal") +pmv = flopy.plot.PlotMapView(modelgrid=modelgrid, ax=ax, layer=0) +cb = pmv.plot_array(head[0]) +pmv.plot_grid(colors="0.5", lw=0.5) +pmv.contour_array(head[0], colors="white", linewidths=1.0) +ax.set_title("Head with the lake stage embedded") +fig.colorbar(cb, ax=ax, shrink=0.7, label="head") +# - + +try: + temp_dir.cleanup() +except (PermissionError, NotADirectoryError): + pass diff --git a/.docs/Notebooks/disv_lake_connections_example.py b/.docs/Notebooks/disv_lake_connections_example.py new file mode 100644 index 000000000..0ef631cd1 --- /dev/null +++ b/.docs/Notebooks/disv_lake_connections_example.py @@ -0,0 +1,224 @@ +# --- +# jupyter: +# jupytext: +# notebook_metadata_filter: all +# text_representation: +# extension: .py +# format_name: light +# format_version: '1.5' +# jupytext_version: 1.14.5 +# kernelspec: +# display_name: Python 3 (ipykernel) +# language: python +# name: python3 +# metadata: +# section: mf6 +# authors: +# - name: Joseph Hughes +# --- + +# # Embedded Lake Connections for a Vertex Grid +# +# `get_lak_connections()` builds Lake (LAK) Package connection data for a lake embedded in a vertex grid as well as a structured one. The cells a lake touches are found from the cells that share an edge with it, and the width of a connection is the length of the shared edge. +# +# This notebook puts the same lake used in the [structured grid example](https://flopy.readthedocs.io/en/latest/Notebooks/dis_lake_connections_example.html) on a vertex grid, and checks that the two agree. + +# + +import sys +from tempfile import TemporaryDirectory + +import matplotlib as mpl +import matplotlib.pyplot as plt +import numpy as np + +import flopy +from flopy.discretization import StructuredGrid, VertexGrid +from flopy.mf6.utils import get_lak_connections + +print(sys.version) +print(f"numpy version: {np.__version__}") +print(f"matplotlib version: {mpl.__version__}") +print(f"flopy version: {flopy.__version__}") +# - + +# ## Build a vertex grid and the equivalent structured grid +# +# The vertex grid discretizes the same domain as the structured grid, one cell2d for every row and column, so the two must give the same lake connections. + +# + +nlay, nrow, ncol = 2, 7, 7 +delr = delc = np.full(7, 100.0) +top = np.full((nrow, ncol), 10.0) +botm = np.array([np.full((nrow, ncol), 0.0), np.full((nrow, ncol), -10.0)]) +ncpl = nrow * ncol + +structured = StructuredGrid(delr=delr, delc=delc, top=top, botm=botm, nlay=nlay) + +xv = np.concatenate(([0.0], np.cumsum(delr))) +yv = delc.sum() - np.concatenate(([0.0], np.cumsum(delc))) +vertices, ivert = [], {} +for i in range(nrow + 1): + for j in range(ncol + 1): + ivert[(i, j)] = len(vertices) + vertices.append((len(vertices), float(xv[j]), float(yv[i]))) + +cell2d = [ + ( + i * ncol + j, + 0.5 * (xv[j] + xv[j + 1]), + 0.5 * (yv[i] + yv[i + 1]), + 4, + ivert[(i, j)], + ivert[(i, j + 1)], + ivert[(i + 1, j + 1)], + ivert[(i + 1, j)], + ) + for i in range(nrow) + for j in range(ncol) +] + +vertex = VertexGrid( + vertices=vertices, + cell2d=cell2d, + top=top.flatten(), + botm=botm.reshape(nlay, ncpl), + nlay=nlay, +) +print(f"Vertex grid: {vertex.nlay} layers of {vertex.ncpl} cells") +# - + +# ## Build the lake connections +# +# The lake occupies the same three by three block in the upper layer. On the vertex grid the lake map has one value for every cell2d rather than a row and a column. + +# + +lake_map_dis = np.full((nlay, nrow, ncol), -1, dtype=np.int32) +lake_map_dis[0, 2:5, 2:5] = 0 +lake_map_disv = lake_map_dis.reshape(nlay, ncpl) + +dis_idomain, dis_conn, dis_data = get_lak_connections( + structured, + np.ma.masked_where(lake_map_dis < 0, lake_map_dis), + idomain=np.ones((nlay, nrow, ncol), dtype=int), + bedleak=0.1, +) +disv_idomain, disv_conn, disv_data = get_lak_connections( + vertex, + np.ma.masked_where(lake_map_disv < 0, lake_map_disv), + idomain=np.ones((nlay, ncpl), dtype=int), + bedleak=0.1, +) + +print(f"structured connections: {dis_conn}") +print(f"vertex connections : {disv_conn}") +print(f"first vertex record : {disv_data[0]}") +# - + +# The two grids give the same connections. A structured cellid is a layer, row, and column while a vertex cellid is a layer and a cell2d number, so the structured cellids are flattened before the two are compared. + + +# + +def normalize(row): + lakeno, _, cellid, claktype, _, _, _, connlen, connwidth = row + if len(cellid) == 3: + k, i, j = cellid + cellid = (k, i * ncol + j) + return (lakeno, *cellid, claktype, connlen, connwidth) + + +assert sorted(map(normalize, dis_data)) == sorted(map(normalize, disv_data)) +assert np.array_equal(dis_idomain.reshape(nlay, ncpl), disv_idomain) +print("The vertex grid reproduces the structured grid connections.") +# - + +# ## Run the model +# +# The connectiondata and the idomain from the vertex grid go straight into a +# DISV model. Constant heads on the left and right edges drive flow across the +# lake. + +# + +temp_dir = TemporaryDirectory() +name = "disv_lake" + +sim = flopy.mf6.MFSimulation(sim_name=name, sim_ws=temp_dir.name, exe_name="mf6") +flopy.mf6.ModflowTdis(sim) +flopy.mf6.ModflowIms( + sim, + linear_acceleration="bicgstab", + outer_dvclose=1e-9, + inner_dvclose=1e-10, + outer_maximum=200, +) +gwf = flopy.mf6.ModflowGwf(sim, modelname=name, newtonoptions="newton under_relaxation") +flopy.mf6.ModflowGwfdisv( + gwf, + nlay=nlay, + ncpl=ncpl, + top=top.flatten(), + botm=botm.reshape(nlay, ncpl), + vertices=vertices, + cell2d=cell2d, + idomain=disv_idomain, +) +flopy.mf6.ModflowGwfic(gwf, strt=8.0) +flopy.mf6.ModflowGwfnpf(gwf, icelltype=1, k=10.0) +flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=[[(0, i * ncol), 9.0] for i in range(nrow)] + + [[(0, i * ncol + ncol - 1), 6.0] for i in range(nrow)], +) +lak = flopy.mf6.ModflowGwflak( + gwf, + stage_filerecord=f"{name}.lak.stage.bin", + nlakes=1, + packagedata=[[0, 7.5, disv_conn[0]]], + connectiondata=disv_data, + perioddata={0: [[0, "rainfall", 0.001]]}, +) +flopy.mf6.ModflowGwfoc(gwf, head_filerecord=f"{name}.hds", saverecord=[("HEAD", "ALL")]) + +sim.write_simulation(silent=True) +success, buff = sim.run_simulation(silent=True) +assert success, "\n".join(buff[-15:]) +print("model converged") +# - + +# ## Put the lake stage into the head array +# +# The lake cells are inactive, so the head array has no value there. Filling +# those cells with the stage of the lake that occupies them gives a single +# array of the water surface that can be plotted in one pass. + +# + +lake_map_masked = np.ma.masked_where(lake_map_disv < 0, lake_map_disv) + +head = gwf.output.head().get_data().reshape(nlay, ncpl) +stage = lak.output.stage().get_data().flatten() +for lake_number in np.unique(lake_map_masked.compressed()): + head[lake_map_masked == lake_number] = stage[lake_number] + +print(f"lake stage: {stage[0]:.3f}") +print(f"cells with no head value: {(np.abs(head) > 1e29).sum()}") +print(f"head range: {head.min():.3f} to {head.max():.3f}") +# - + +# ## Head map +# +# The lake reads as part of the water surface rather than as a hole in it. + +# + +fig, ax = plt.subplots(figsize=(6, 6)) +ax.set_aspect("equal") +pmv = flopy.plot.PlotMapView(modelgrid=vertex, ax=ax, layer=0) +cb = pmv.plot_array(head[0]) +pmv.plot_grid(colors="0.5", lw=0.5) +pmv.contour_array(head[0], colors="white", linewidths=1.0) +ax.set_title("Head with the lake stage embedded") +fig.colorbar(cb, ax=ax, shrink=0.7, label="head") +# - + +try: + temp_dir.cleanup() +except (PermissionError, NotADirectoryError): + pass From eb8d20c3f495e4924ddac71fc06b814211ae30e3 Mon Sep 17 00:00:00 2001 From: Joseph Hughes Date: Wed, 5 Aug 2026 12:46:34 -0500 Subject: [PATCH 2/3] docs(lakpak_utils): match the voice of the other example notebooks The two notebooks were written in an impersonal voice and never named MODFLOW 6, unlike the notebooks added alongside them. They now use the first person for what the notebook does and name the program whose package input is being built. --- .../Notebooks/dis_lake_connections_example.py | 18 +++++++++--------- .../Notebooks/disv_lake_connections_example.py | 18 +++++++++--------- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/.docs/Notebooks/dis_lake_connections_example.py b/.docs/Notebooks/dis_lake_connections_example.py index 51949613e..6b6653ca3 100644 --- a/.docs/Notebooks/dis_lake_connections_example.py +++ b/.docs/Notebooks/dis_lake_connections_example.py @@ -19,9 +19,9 @@ # # Embedded Lake Connections for a Structured Grid # -# A lake embedded in the model domain replaces the cells it occupies, so the Lake (LAK) Package needs a connection to every active cell that touches it. `get_lak_connections()` builds those connections from an array of lake numbers, and returns the idomain with the lake cells deactivated, the number of connections in each lake, and the connectiondata block for the package. +# A lake embedded in the model domain replaces the cells it occupies, so the MODFLOW 6 Lake (LAK) Package needs a connection to every active cell that touches it. `get_lak_connections()` builds those connections from an array of lake numbers, and returns the idomain with the lake cells deactivated, the number of connections in each lake, and the connectiondata block for the package. # -# This notebook builds a lake embedded in one layer of a structured grid. The companion notebook does the same thing on a vertex grid. +# We build a lake embedded in one layer of a structured grid, and put the simulated lake stage into the head array so that a single map shows the water surface. The companion notebook does the same thing on a vertex grid. # + import sys @@ -44,7 +44,7 @@ # ## Build the grid and the lake # -# The grid is two layers of 7 rows and 7 columns. The lake occupies a three by three block in the upper layer, and is given lake number 0. Cells that are not part of a lake are masked. +# We use a grid of two layers with 7 rows and 7 columns. The lake occupies a three by three block in the upper layer, and is given lake number 0. Cells that are not part of a lake are masked. # + nlay, nrow, ncol = 2, 7, 7 @@ -63,7 +63,7 @@ # ## Build the lake connections # -# The lake map is the same shape as the model grid, so the lake is embedded rather than sitting on top of the model. A bed leakance is given for every connection. +# The lake map is the same shape as the model grid, so the lake is embedded rather than sitting on top of the model. We give the same bed leakance to every connection. # + idomain, connection_dict, connectiondata = get_lak_connections( @@ -101,8 +101,8 @@ # ## Run the model # -# The connectiondata and the idomain go straight into the model. Constant heads -# on the left and right edges drive flow across the lake. +# The connectiondata and the idomain go straight into the MODFLOW 6 model. We +# set constant heads on the left and right edges to drive flow across the lake. # + temp_dir = TemporaryDirectory() @@ -154,9 +154,9 @@ # ## Put the lake stage into the head array # -# The lake cells are inactive, so the head array has no value there. Filling -# those cells with the stage of the lake that occupies them gives a single -# array of the water surface that can be plotted in one pass. +# The lake cells are inactive, so the head array has no value there. We fill +# those cells with the stage of the lake that occupies them, which gives a +# single array of the water surface that can be plotted in one pass. # + head = gwf.output.head().get_data() diff --git a/.docs/Notebooks/disv_lake_connections_example.py b/.docs/Notebooks/disv_lake_connections_example.py index 0ef631cd1..7d6267061 100644 --- a/.docs/Notebooks/disv_lake_connections_example.py +++ b/.docs/Notebooks/disv_lake_connections_example.py @@ -19,9 +19,9 @@ # # Embedded Lake Connections for a Vertex Grid # -# `get_lak_connections()` builds Lake (LAK) Package connection data for a lake embedded in a vertex grid as well as a structured one. The cells a lake touches are found from the cells that share an edge with it, and the width of a connection is the length of the shared edge. +# `get_lak_connections()` builds MODFLOW 6 Lake (LAK) Package connection data for a lake embedded in a vertex grid as well as a structured one. The cells a lake touches are found from the cells that share an edge with it, and the width of a connection is the length of the shared edge. # -# This notebook puts the same lake used in the [structured grid example](https://flopy.readthedocs.io/en/latest/Notebooks/dis_lake_connections_example.html) on a vertex grid, and checks that the two agree. +# We put the same lake used in the [structured grid example](https://flopy.readthedocs.io/en/latest/Notebooks/dis_lake_connections_example.html) on a vertex grid, check that the two agree, and put the simulated lake stage into the head array so that a single map shows the water surface. # + import sys @@ -43,7 +43,7 @@ # ## Build a vertex grid and the equivalent structured grid # -# The vertex grid discretizes the same domain as the structured grid, one cell2d for every row and column, so the two must give the same lake connections. +# We build a vertex grid that discretizes the same domain as the structured grid, one cell2d for every row and column, so the two must give the same lake connections. # + nlay, nrow, ncol = 2, 7, 7 @@ -114,7 +114,7 @@ print(f"first vertex record : {disv_data[0]}") # - -# The two grids give the same connections. A structured cellid is a layer, row, and column while a vertex cellid is a layer and a cell2d number, so the structured cellids are flattened before the two are compared. +# The two grids give the same connections. A structured cellid is a layer, row, and column while a vertex cellid is a layer and a cell2d number, so we flatten the structured cellids before comparing the two. # + @@ -134,8 +134,8 @@ def normalize(row): # ## Run the model # # The connectiondata and the idomain from the vertex grid go straight into a -# DISV model. Constant heads on the left and right edges drive flow across the -# lake. +# MODFLOW 6 DISV model. We set constant heads on the left and right edges to +# drive flow across the lake. # + temp_dir = TemporaryDirectory() @@ -186,9 +186,9 @@ def normalize(row): # ## Put the lake stage into the head array # -# The lake cells are inactive, so the head array has no value there. Filling -# those cells with the stage of the lake that occupies them gives a single -# array of the water surface that can be plotted in one pass. +# The lake cells are inactive, so the head array has no value there. We fill +# those cells with the stage of the lake that occupies them, which gives a +# single array of the water surface that can be plotted in one pass. # + lake_map_masked = np.ma.masked_where(lake_map_disv < 0, lake_map_disv) From df277c14d2659cbd074e76d529534b93ff9506fe Mon Sep 17 00:00:00 2001 From: Joseph Hughes Date: Wed, 5 Aug 2026 13:36:55 -0500 Subject: [PATCH 3/3] docs(lakpak_utils): drop an unused import and compare sizes with a tolerance The structured notebook imported Path without using it. The vertex notebook compared the connection lengths and widths of the two grids for exact equality, which is brittle because the two grids compute them differently. The cells and the connection types are still compared exactly, and the lengths and widths to within rounding. --- .docs/Notebooks/dis_lake_connections_example.py | 1 - .docs/Notebooks/disv_lake_connections_example.py | 10 ++++++++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/.docs/Notebooks/dis_lake_connections_example.py b/.docs/Notebooks/dis_lake_connections_example.py index 6b6653ca3..b27f17967 100644 --- a/.docs/Notebooks/dis_lake_connections_example.py +++ b/.docs/Notebooks/dis_lake_connections_example.py @@ -25,7 +25,6 @@ # + import sys -from pathlib import Path from tempfile import TemporaryDirectory import matplotlib as mpl diff --git a/.docs/Notebooks/disv_lake_connections_example.py b/.docs/Notebooks/disv_lake_connections_example.py index 7d6267061..237fbbdd6 100644 --- a/.docs/Notebooks/disv_lake_connections_example.py +++ b/.docs/Notebooks/disv_lake_connections_example.py @@ -123,10 +123,16 @@ def normalize(row): if len(cellid) == 3: k, i, j = cellid cellid = (k, i * ncol + j) - return (lakeno, *cellid, claktype, connlen, connwidth) + return (lakeno, *cellid, claktype), (connlen, connwidth) -assert sorted(map(normalize, dis_data)) == sorted(map(normalize, disv_data)) +dis_sorted = sorted(map(normalize, dis_data)) +disv_sorted = sorted(map(normalize, disv_data)) + +# the cells and the connection types must match exactly, and the lengths and +# widths to within rounding, because the two grids compute them differently +assert [cells for cells, _ in dis_sorted] == [cells for cells, _ in disv_sorted] +assert np.allclose([size for _, size in dis_sorted], [size for _, size in disv_sorted]) assert np.array_equal(dis_idomain.reshape(nlay, ncpl), disv_idomain) print("The vertex grid reproduces the structured grid connections.") # -