feat(lakpak_utils): add support for VertexGrid embedded lakes - #2786
feat(lakpak_utils): add support for VertexGrid embedded lakes#2786afc98 wants to merge 6 commits into
Conversation
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 Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## develop #2786 +/- ##
===========================================
+ Coverage 55.5% 72.6% +17.1%
===========================================
Files 644 660 +16
Lines 124135 132425 +8290
===========================================
+ Hits 68947 96257 +27310
+ Misses 55188 36168 -19020
🚀 New features to boost your workflow:
|
jdhughes-dev
left a comment
jdhughes-dev
left a comment
There was a problem hiding this comment.
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.
|
|
||
| shared = tuple(cell_iverts & set(iverts[nicpl])) | ||
|
|
||
| if shared is None or len(shared) != 2: |
There was a problem hiding this comment.
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.5Intersecting 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.
|
|
||
| k, icpl = cell_index | ||
| node = k * ncpl + icpl | ||
| neighbors = modelgrid.neighbors(node=node, method="rook") |
There was a problem hiding this comment.
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.
| return cellids, claktypes, belevs, televs, connlens, connwidths | ||
|
|
||
|
|
||
| def __vertex_lake_connections( |
There was a problem hiding this comment.
__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.
|
|
||
| for neighbor in neighbors: | ||
| # Convert global node number to (layer, cell-per-layer) | ||
| nk = neighbor // ncpl |
There was a problem hiding this comment.
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.
| return np.linalg.norm(P - A) | ||
|
|
||
| t = np.dot(AP, AB) / denom | ||
| t = np.clip(t, 0.0, 1.0) |
There was a problem hiding this comment.
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.
|
|
||
| idomain = np.ones((1, modelgrid.ncpl), dtype=int) | ||
|
|
||
| idomain_out, pakdata, connectiondata = get_lak_connections( |
There was a problem hiding this comment.
idomain_out is unused.
| for name in properties: | ||
| monkeypatch.setattr(VertexGrid, name, counted_property(name)) | ||
|
|
||
| def unexpected_get_shared_edge(*args, **kwargs): |
There was a problem hiding this comment.
get_shared_edge() is never called by __vertex_lake_connections, so this guard cannot fire.
| return property(getter) | ||
|
|
||
| for name in properties: | ||
| monkeypatch.setattr(VertexGrid, name, counted_property(name)) |
There was a problem hiding this comment.
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.
|
|
||
|
|
||
| @requires_exe("mf6") | ||
| def test_disv_lake_run(function_tmpdir): |
There was a problem hiding this comment.
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.
2 participants
Summary
idomainTesting
pytest -q autotest/test_lake_connections.pycd autotest; pytest -v -n autoruff check autotest/test_lake_connections.py flopy/mf6/utils/lakpak_utils.pyruff format --check autotest/test_lake_connections.py