feat(lakpak_utils): add support for VertexGrid embedded lakes by afc98 · Pull Request #2786 · modflowpy/flopy

afc98 · 2026-07-28T21:10:19Z

Summary

add embedded lake connection generation for VertexGrids
calculate horizontal connection widths and lengths from shared edges
support horizontal and vertical lake connections across multiple layers
respect inactive neighbouring cells and update idomain
add unit and MODFLOW 6 integration tests
cover asymmetric cells and lakes below the top layer
reformat lakpak_utils.py to pass ruff format --check

Testing

pytest -q autotest/test_lake_connections.py
cd autotest; pytest -v -n auto
ruff check autotest/test_lake_connections.py flopy/mf6/utils/lakpak_utils.py
ruff format --check autotest/test_lake_connections.py

Added a __vertex_lake_connections function that provides the same outputs as __structured_lake_connections but for vertex grids.

Calculate elevation, vertex, and cell-center arrays once when generating embedded DISV lake connections. Reuse cached vertex indices to find shared edges instead of repeatedly rebuilding VertexGrid geometry for every neighboring cell. Add regression coverage to ensure expensive grid properties are accessed only once.

codecov · 2026-07-31T22:20:15Z

Codecov Report

❌ Patch coverage is 92.50000% with 6 lines in your changes missing coverage. Please review.
✅ Project coverage is 72.6%. Comparing base (556c088) to head (c162cf6).
⚠️ Report is 182 commits behind head on develop.

Files with missing lines	Patch %	Lines
flopy/mf6/utils/lakpak_utils.py	92.5%	6 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           develop    #2786      +/-   ##
===========================================
+ Coverage     55.5%    72.6%   +17.1%     
===========================================
  Files          644      660      +16     
  Lines       124135   132441    +8306     
===========================================
+ Hits         68947    96275   +27328     
+ Misses       55188    36166   -19022

Files with missing lines	Coverage Δ
flopy/mf6/utils/lakpak_utils.py	`91.6% <92.5%> (+4.2%)`	⬆️

... and 580 files with indirect coverage changes

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

jdhughes-dev

Verified the implementation against __structured_lake_connections: for a vertex grid that discretizes the same domain as a structured grid (irregular delr/delc, interior/corner/edge lake cells, lakes in layer 1 and spanning layers, inactive neighbors, multiple lakes) it reproduces the structured connectiondata, connection_dict, and idomain exactly. Triangles, closed cell2d rings, gridgen-style hanging nodes, and grid offset/rotation are also handled correctly.

Two things need fixing before this goes in: connections are silently dropped when two cells share more than two vertices, and a lake cell that produces no connections passes without warning. Both are noted inline.

Test coverage needs to grow with the feature. Every DISV test grid is an axis-aligned unit square apart from one 3x1 cell, so none of the unstructured geometry this exists for is exercised. Please add a DIS/DISV equivalence test - build a structured grid and the vertex grid that discretizes the same domain, then assert identical connectiondata, connection_dict, and idomain - plus non-rectangular cells, two lakes including adjacent ones, and the degenerate cases below.

jdhughes-dev/flopy@fix/disv-lake-shared-edges has one way to do both, for reference: jdhughes-dev/flopy@13dd04a...jdhughes-dev:flopy:fix/disv-lake-shared-edges

The reformatting is correct and worth a line in the description - lakpak_utils.py is the only file in the repo that currently fails ruff format --check, so this brings it into line rather than adding churn.

jdhughes-dev · 2026-07-31T23:20:05Z

+
+        shared = tuple(cell_iverts & set(iverts[nicpl]))
+
+        if shared is None or len(shared) != 2:


Connections are dropped when two cells share more than two vertices. Grid._set_neighbors(method="rook") matches edges, so a cell boundary split by a vertex that both cells carry makes the two cells rook neighbors with three shared vertices, and this test discards the connection. The lake then reports too few connections, and when it is the only neighbor iconn stays 0, so the lake cell is left active in idomain.

vertices = [(0,0.,0.),(1,1.,0.),(2,1.,1.),(3,1.,2.),(4,0.,2.),(5,2.,0.),(6,2.,2.)] cell2d = [(0,0.5,1.0,5,0,1,2,3,4), (1,1.5,1.0,5,1,5,6,3,2)] get_lak_connections(grid, np.array([[0, -1]]), bedleak=1.0) # connectiondata [] nconn {0: 0} idomain [[1 1]] # expected: one horizontal connection, connwidth 2.0, connlen 0.5

Intersecting the cells' edges rather than their vertices fixes this and cannot disagree with the connectivity neighbors() reports: sum connwidth over the shared edges and take connlen as the distance to the nearest.

tuple(...) is never None, so shared is None is also dead.

I replaced the shared-vertex count with shared-edge intersection matching the edge representation used by rook connectivity. Boundaries split by hanging vertices now contribute all shared segments: their lengths are summed for connwidth, and the nearest segment determines connlen. I also removed the unreachable shared is None check and added the provided regression case, including the expected connection geometry and idomain update.

jdhughes-dev · 2026-07-31T23:20:05Z

+
+    k, icpl = cell_index
+    node = k * ncpl + icpl
+    neighbors = modelgrid.neighbors(node=node, method="rook")


A DISV grid whose cells do not carry hanging nodes at a refinement boundary has no rook neighbors, so an embedded lake cell produces no connections, connection_dict reports 0, and idomain is not deactivated. Nothing warns, and an nlakeconn of 0 reaches the LAK package. Gridgen and Voronoi output are unaffected. Please warn when an embedded lake cell yields no connections.

Added a UserWarning when an embedded lake has zero connections after all its cells have been processed. I placed the warning at lake level so legitimate interior cells in a multi-cell lake do not generate false warnings. The warning regression test confirms zero connections are reported and idomain remains unchanged for incompatible refinement topology.

jdhughes-dev · 2026-07-31T23:20:05Z

    return cellids, claktypes, belevs, televs, connlens, connwidths
+
+
+def __vertex_lake_connections(


__structured_lake_connections wraps its body in if idomain[cell_index] > 0; this function has no equivalent. Unreachable today because lake_map[idomain < 1] = -1 runs first, but please keep the two branches symmetric.

__vertex_lake_connections now has the same explicit idomain[cell_index] > 0 guard as the structured implementation.

jdhughes-dev · 2026-07-31T23:20:05Z

+
+    for neighbor in neighbors:
+        # Convert global node number to (layer, cell-per-layer)
+        nk = neighbor // ncpl


neighbors() builds 2D adjacency and re-offsets by layer, so nk is always k. The decomposition implies cross-layer handling that does not exist - drop it or note the invariant.

The comment above is also the only capitalized one in the file.

The code now uses k explicitly for horizontal neighbors and only derives the cell-per-layer index from the returned node. The comment documents that neighbors() returns 2D adjacency offset to the requested layer, and its capitalization has been normalised.

jdhughes-dev · 2026-07-31T23:20:05Z

+        return np.linalg.norm(P - A)
+
+    t = np.dot(AP, AB) / denom
+    t = np.clip(t, 0.0, 1.0)


Clamping t returns the distance to an edge endpoint when the neighbor centroid projects off the shared edge, not the face-normal length LAK uses for conductance. Fine to keep, but say so in the comment.

The helper now describes this as the shortest distance to a line segment and explicitly documents that projections outside the segment are clamped to the nearest endpoint.

jdhughes-dev · 2026-07-31T23:20:05Z

+
+    idomain = np.ones((1, modelgrid.ncpl), dtype=int)
+
+    idomain_out, pakdata, connectiondata = get_lak_connections(


idomain_out is unused.

The unused return value is now assigned to _ .

jdhughes-dev · 2026-07-31T23:20:05Z

+    for name in properties:
+        monkeypatch.setattr(VertexGrid, name, counted_property(name))
+
+    def unexpected_get_shared_edge(*args, **kwargs):


get_shared_edge() is never called by __vertex_lake_connections, so this guard cannot fire.

jdhughes-dev · 2026-07-31T23:20:06Z

+        return property(getter)
+
+    for name in properties:
+        monkeypatch.setattr(VertexGrid, name, counted_property(name))


xcellcenters and ycellcenters are defined on Grid, not VertexGrid, so monkeypatch's undo sets the inherited property objects directly on VertexGrid. The exact == 1 counts also break on any refactor that reads a property twice and depend on neighbors() internals, hence the pre-warm above. <= would be enough.

Inherited cell-center properties are now patched on Grid, while VertexGrid-specific properties are patched on VertexGrid. The assertion now checks that each property is accessed at most once instead of requiring an exact internal access pattern.

jdhughes-dev · 2026-07-31T23:20:06Z

+
+
+@requires_exe("mf6")
+def test_disv_lake_run(function_tmpdir):


Cells 1-4 touch only cell 0, which is deactivated, so this is four isolated cells at steady state with rainfall 0.0 and the test asserts only that mf6 ran. A DISV twin of test_embedded_lak_ex01 comparing heads and lake stage against the DIS run would exercise the connection geometry.

The test now builds equivalent connected DIS and DISV models with boundary-driven groundwater flow and lake rainfall, runs both, and compares the resulting aquifer heads and lake stage.

Build lake connections from shared edges so boundaries split by hanging vertices retain complete widths and distances. Warn when an embedded lake has no active-cell connections instead of silently emitting zero connections. Expand DIS/DISV equivalence and MODFLOW 6 result tests to cover unstructured geometry, inactive cells, adjacent lakes, and multilayer cases.

afc98 added 5 commits July 25, 2026 13:17

feat(get_lak_connections): add support for embedded disv modelgrids

6a4cd32

Added a __vertex_lake_connections function that provides the same outputs as __structured_lake_connections but for vertex grids.

test: tests for get_lak_connections with vertex grids

9cab1ed

fix(lakpak_utils): correct embedded DISV lake connection geometry

423f3cb

test(lakpak_utils): cover embedded DISV connection geometry

95c5cb8

Merge branch 'modflowpy:develop' into fix/disv-lake-connection-geometry

db285a0

afc98 marked this pull request as draft July 31, 2026 13:12

afc98 marked this pull request as ready for review July 31, 2026 15:28

jdhughes-dev requested changes Jul 31, 2026

View reviewed changes

jdhughes-dev mentioned this pull request Jul 31, 2026

fix(lakpak_utils): use shared edges for embedded DISV lake connections afc98/flopy#1

Open

afc98 added 2 commits August 1, 2026 21:43

style(lakpak_utils): lowercase warning message

c162cf6


		shared = tuple(cell_iverts & set(iverts[nicpl]))

		if shared is None or len(shared) != 2:

		return cellids, claktypes, belevs, televs, connlens, connwidths


		def __vertex_lake_connections(


		idomain = np.ones((1, modelgrid.ncpl), dtype=int)

		idomain_out, pakdata, connectiondata = get_lak_connections(



		@requires_exe("mf6")
		def test_disv_lake_run(function_tmpdir):

Uh oh!

Conversation

afc98 commented Jul 28, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Summary

Testing

Uh oh!

codecov Bot commented Jul 31, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

jdhughes-dev left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

afc98 commented Jul 28, 2026 •

edited

Loading

codecov Bot commented Jul 31, 2026 •

edited

Loading