From 56416cba4422c034c72e40f035a25f3c7976f742 Mon Sep 17 00:00:00 2001 From: Joseph Hughes Date: Tue, 4 Aug 2026 13:33:36 -0500 Subject: [PATCH 1/2] feat(gnc): add ghost node correction data generation Ghost node correction data could not be built with FloPy, so the GNC Package was difficult to use on a quadtree grid even though the correction is what makes the control volume formulation accurate where a coarse cell meets a finer one. The new flopy.utils.gnc module computes the data from a model grid, a grid conforming array of refinement levels, and the grid connectivity, and converts it to MODFLOW 6 GNC Package input. Connectivity is taken from the grid when it is not supplied, either from the iac and ja an unstructured grid carries or from the cells that share an edge. Records with fewer contributing cells than numalphaj repeat a cell and split its contributing factor rather than pad with a cellid of zero, so that the same records are also valid MODFLOW-USG input. The computed data reproduces the ghost nodes gridgen writes to qtg.gnc.dat exactly, which is checked over five grid configurations. Also removes a duplicate MfUsgGnc entry from the flopy.mfusg __all__ list, which ruff check reports as RUF068 and which blocks the lint job. --- autotest/test_gnc.py | 329 ++++++++++++++++++++++++++++++++ flopy/mfusg/__init__.py | 1 - flopy/utils/__init__.py | 1 + flopy/utils/gnc.py | 406 ++++++++++++++++++++++++++++++++++++++++ 4 files changed, 736 insertions(+), 1 deletion(-) create mode 100644 autotest/test_gnc.py create mode 100644 flopy/utils/gnc.py diff --git a/autotest/test_gnc.py b/autotest/test_gnc.py new file mode 100644 index 0000000000..dcff6a7585 --- /dev/null +++ b/autotest/test_gnc.py @@ -0,0 +1,329 @@ +""" +Tests for ghost node correction (GNC) data computed from a grid. + +Cases: + - synthetic : a hand built grid with known contributing cells and factors. + - gridgen : the computed data must reproduce gridgen's qtg.gnc.dat. +""" + +import numpy as np +import pytest +from modflow_devtools.markers import requires_exe, requires_pkg + +import flopy +from flopy.discretization import UnstructuredGrid, VertexGrid +from flopy.utils.gnc import ( + _check_gnc, + get_gnc, + get_gnc_dtype, + get_gridprops_gnc6, + get_numalphaj, +) +from flopy.utils.gridgen import Gridgen + + +def synthetic_grid(connectivity=True): + """A coarse cell with four fine cells north and two fine cells east + + Cell 0 is a size 4 cell centered on the origin. Cells 1 and 2 are size 2 + cells on its east face and cells 3 to 6 are size 1 cells on its north + face, so the north face has a two level refinement jump. + """ + centers = [ + (0.0, 0.0, 4.0), + (3.0, 1.0, 2.0), + (3.0, -1.0, 2.0), + (-1.5, 3.0, 1.0), + (-0.5, 3.0, 1.0), + (0.5, 3.0, 1.0), + (1.5, 3.0, 1.0), + ] + vertices, iverts = [], [] + for x, y, size in centers: + half = size / 2.0 + iv = [] + for vx, vy in [ + (x - half, y - half), + (x + half, y - half), + (x + half, y + half), + (x - half, y + half), + ]: + iv.append(len(vertices)) + vertices.append([len(vertices), vx, vy]) + iverts.append(iv) + + # every cell is connected to cell 0 only + conn = {} + if connectivity: + conn["iac"] = np.array([7, 2, 2, 2, 2, 2, 2]) + conn["ja"] = np.array([0, 1, 2, 3, 4, 5, 6, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0]) + grid = UnstructuredGrid( + vertices=vertices, + iverts=iverts, + xcenters=np.array([c[0] for c in centers]), + ycenters=np.array([c[1] for c in centers]), + ncpl=np.array([len(centers)]), + **conn, + ) + level = np.array([0, 1, 1, 2, 2, 2, 2]) + return grid, level + + +def gnc_key(gnc): + """Sorted view of ghost node records for order independent comparison""" + numalphaj = get_numalphaj(gnc) + nodes = np.sort(np.column_stack([gnc[f"j{i}"] for i in range(numalphaj)]), axis=1) + alpha = np.sort( + np.column_stack([gnc[f"alpha{i}"] for i in range(numalphaj)]), axis=1 + ) + key = np.column_stack([gnc["n"], gnc["m"], nodes, alpha]) + return key[np.lexsort(key.T[::-1])] + + +def build_gridgen(ws, nlay=3, level=3, layers=None, **kwargs): + """Build a gridgen grid with a refined block in the middle""" + from shapely.geometry import Polygon + + botm = [1.0 - k * (1.0 / nlay) for k in range(1, nlay + 1)] + sim = flopy.mf6.MFSimulation(sim_name="base", sim_ws=ws) + gwf = flopy.mf6.ModflowGwf(sim, modelname="base") + flopy.mf6.ModflowGwfdis( + gwf, nlay=nlay, nrow=10, ncol=10, delr=1.0, delc=1.0, top=1.0, botm=botm + ) + g = Gridgen(gwf.modelgrid, model_ws=ws, **kwargs) + g.add_refinement_features( + [Polygon([(4, 4), (6, 4), (6, 6), (4, 6)])], + "polygon", + level, + range(nlay) if layers is None else layers, + ) + g.build() + return g + + +def test_get_gnc_dtype(): + dtype = get_gnc_dtype(3) + assert dtype.names == ( + "n", + "m", + "j0", + "j1", + "j2", + "alpha0", + "alpha1", + "alpha2", + ) + assert get_numalphaj(np.recarray((0,), dtype=dtype)) == 3 + + +@pytest.mark.parametrize("use_level", [True, False]) +def test_get_gnc_synthetic(use_level): + grid, level = synthetic_grid() + gnc = get_gnc(grid, level=level if use_level else None) + + # cell 0 has four contributing cells on its north face + assert get_numalphaj(gnc) == 4 + assert len(gnc) == 3 + + # connections without a contributing cell on the offset side are dropped + assert sorted(zip(gnc["n"], gnc["m"])) == [(0, 1), (0, 5), (0, 6)] + + records = {(rec["n"], rec["m"]): rec for rec in gnc} + + # cell 1 is offset 1.0 from cell 0 and the north cells are 3.0 away, so + # the factors total 1/3 and are shared by four contributing cells + rec = records[(0, 1)] + assert [rec[f"j{i}"] for i in range(4)] == [3, 4, 5, 6] + assert np.allclose([rec[f"alpha{i}"] for i in range(4)], 1.0 / 12.0) + + # cell 6 is offset 1.5 and the east cells are 3.0 away, so the factors + # total 0.5 shared by two contributing cells, the first of which is + # repeated three times to fill numalphaj + rec = records[(0, 6)] + assert [rec[f"j{i}"] for i in range(4)] == [1, 1, 1, 2] + assert np.allclose( + [rec[f"alpha{i}"] for i in range(4)], [1 / 12, 1 / 12, 1 / 12, 0.25] + ) + + +def factor_totals(gnc): + """Total contributing factor of each cell, keyed by ghost node""" + numalphaj = get_numalphaj(gnc) + totals = {} + for rec in gnc: + entry = totals.setdefault((rec["n"], rec["m"]), {}) + for i in range(numalphaj): + node = rec[f"j{i}"] + entry[node] = entry.get(node, 0.0) + rec[f"alpha{i}"] + return totals + + +def test_get_gnc_padding_preserves_factors(): + grid, level = synthetic_grid() + unpadded = get_gnc(grid, level=level) + padded = get_gnc(grid, level=level, numalphaj=6) + assert get_numalphaj(unpadded) == 4 + assert get_numalphaj(padded) == 6 + + # padding repeats a contributing cell and splits its factor, so the total + # factor of every cell must be unchanged + expected, actual = factor_totals(unpadded), factor_totals(padded) + assert expected.keys() == actual.keys() + for ghost, cells in expected.items(): + assert cells.keys() == actual[ghost].keys() + for node, alpha in cells.items(): + assert np.allclose(alpha, actual[ghost][node]) + assert sum(cells.values()) < 1.0 + + +def test_get_gnc_numalphaj_too_small(): + grid, level = synthetic_grid() + with pytest.raises(ValueError, match="more than numalphaj"): + get_gnc(grid, level=level, numalphaj=2) + + +def test_get_gnc_supplied_connectivity(): + grid, level = synthetic_grid() + bare, _ = synthetic_grid(connectivity=False) + assert bare.iac is None + + # the synthetic cells do not share vertices, so connectivity cannot be + # built from shared edges and has to be supplied + assert len(get_gnc(bare, level=level)) == 0 + + gnc = get_gnc(bare, level=level, iac=grid.iac, ja=grid.ja) + assert np.array_equal(gnc_key(gnc), gnc_key(get_gnc(grid, level=level))) + + +def test_check_gnc(): + dtype = get_gnc_dtype(2) + gnc = np.recarray((1,), dtype=dtype) + gnc[0] = (0, 1, 2, 2, 0.6, 0.6) + with pytest.raises(ValueError, match="must be less than one"): + _check_gnc(gnc) + + gnc[0] = (0, 1, 2, 2, 0.1, 0.1) + _check_gnc(gnc) + + # cell 0 is connected to cell 2 but not to cell 1 + ia = np.array([0, 2, 3, 5]) + ja = np.array([0, 2, 1, 2, 0]) + with pytest.raises(ValueError, match="is not connected to cell"): + _check_gnc(gnc, ia=ia, ja=ja) + + +def test_get_gridprops_gnc6_requires_ncpl(): + grid, level = synthetic_grid() + gnc = get_gnc(grid, level=level) + with pytest.raises(ValueError, match="ncpl is required"): + get_gridprops_gnc6(gnc, dis_type="disv") + with pytest.raises(ValueError, match="Unknown dis_type"): + get_gridprops_gnc6(gnc, dis_type="dis") + + gridprops = get_gridprops_gnc6(gnc, dis_type="disu") + assert gridprops["numalphaj"] == 4 + assert gridprops["numgnc"] == len(gnc) + assert gridprops["gncdata"][0][0] == (gnc["n"][0],) + + +@requires_exe("gridgen") +@requires_pkg("shapely", "geopandas") +@pytest.mark.parametrize( + "nlay,level,layers,smoothing", + [ + (3, 3, None, 1), + (1, 2, None, 1), + (1, 4, None, 1), + (3, 3, [0], 1), + (1, 3, None, 2), + ], +) +def test_get_gnc_matches_gridgen(function_tmpdir, nlay, level, layers, smoothing): + g = build_gridgen( + function_tmpdir, + nlay=nlay, + level=level, + layers=layers, + smoothing_level_horizontal=smoothing, + smoothing_level_vertical=smoothing, + ) + # gridgen writes the ghost node data whenever it exports a grid + expected = np.atleast_1d( + np.genfromtxt(function_tmpdir / "qtg.gnc.dat", dtype=get_gnc_dtype(2)) + ) + for name in ("n", "m", "j0", "j1"): + expected[name] -= 1 + assert len(expected) > 0 + + grid = UnstructuredGrid(**g.get_gridprops_unstructuredgrid()) + gnc = get_gnc(grid, numalphaj=2) + + assert len(gnc) == len(expected) + # gridgen writes the factors with six significant digits + assert np.allclose(gnc_key(gnc), gnc_key(expected), atol=2.0e-6) + + +def square_cells(centers): + """Build vertices and iverts for a list of (x, y, size) squares""" + vertices, iverts = [], [] + for x, y, size in centers: + half = size / 2.0 + iv = [] + for vx, vy in [ + (x - half, y - half), + (x + half, y - half), + (x + half, y + half), + (x - half, y + half), + ]: + iv.append(len(vertices)) + vertices.append([len(vertices), vx, vy]) + iverts.append(iv) + return vertices, iverts + + +def test_get_gnc_ihc(): + """A connection marked vertical is not corrected and does not contribute""" + grid, level = synthetic_grid() + + # cells 0 and 1 are connected at positions 1 and 8 of ja + ihc = np.ones(len(grid.ja), dtype=int) + ihc[[1, 8]] = 0 + gnc = get_gnc(grid, level=level, ihc=ihc) + + # the ghost node on the 0 to 1 connection is gone, and cell 1 is no longer + # available as a contributing cell, which leaves cell 2 on its own + assert sorted((rec["n"], rec["m"]) for rec in gnc) == [(0, 5), (0, 6)] + assert get_numalphaj(gnc) == 1 + assert np.allclose(sorted(gnc["alpha0"]), [1.0 / 6.0, 0.5]) + + +def test_get_gnc_aligned_finer_cell(): + """A finer cell centered on the face of a coarse cell needs no ghost node""" + centers = [(0.0, 0.0, 4.0), (3.0, 0.0, 2.0)] + vertices, iverts = square_cells(centers) + grid = UnstructuredGrid( + vertices=vertices, + iverts=iverts, + xcenters=np.array([c[0] for c in centers]), + ycenters=np.array([c[1] for c in centers]), + ncpl=np.array([len(centers)]), + iac=np.array([2, 2]), + ja=np.array([0, 1, 1, 0]), + ) + assert len(get_gnc(grid, level=np.array([0, 1]))) == 0 + + +def test_get_gnc_connectivity_needs_constant_ncpl(): + """Shared edges cannot give the layer layout when ncpl varies by layer""" + centers = [(0.5, 0.5, 1.0), (1.5, 0.5, 1.0), (2.5, 0.5, 1.0), (0.5, 1.5, 1.0)] + vertices, iverts = square_cells(centers) + grid = UnstructuredGrid( + vertices=vertices, + iverts=iverts, + xcenters=np.array([c[0] for c in centers]), + ycenters=np.array([c[1] for c in centers]), + ncpl=np.array([3, 1]), + ) + assert grid.iac is None + with pytest.raises(ValueError, match="different"): + get_gnc(grid) diff --git a/flopy/mfusg/__init__.py b/flopy/mfusg/__init__.py index 579db55f97..b5d4b621af 100644 --- a/flopy/mfusg/__init__.py +++ b/flopy/mfusg/__init__.py @@ -36,7 +36,6 @@ "MfUsgDpt", "MfUsgEvt", "MfUsgGnc", - "MfUsgGnc", "MfUsgLak", "MfUsgLpf", "MfUsgMdt", diff --git a/flopy/utils/__init__.py b/flopy/utils/__init__.py index 0da0e3d107..b4129b4a78 100644 --- a/flopy/utils/__init__.py +++ b/flopy/utils/__init__.py @@ -29,6 +29,7 @@ from .formattedfile import FormattedHeadFile get_modflow = get_modflow_module.run_main +from .gnc import get_gnc, get_gridprops_gnc6 from .gridintersect import GridIntersect from .hfb_util import make_hfb_array from .mflistfile import ( diff --git a/flopy/utils/gnc.py b/flopy/utils/gnc.py new file mode 100644 index 0000000000..4ea3cf2d31 --- /dev/null +++ b/flopy/utils/gnc.py @@ -0,0 +1,406 @@ +""" +Ghost node correction (GNC) data for quadtree-like grids. + +Ghost node data is computed from a grid, a grid conforming array of refinement +levels, and the grid connectivity by :func:`get_gnc`, and is converted to +MODFLOW 6 package input by :func:`get_gridprops_gnc6`. +""" + +import numpy as np + + +def get_gnc_dtype(numalphaj): + """ + Get the record dtype for ghost node data with numalphaj contributing + cells + + Parameters + ---------- + numalphaj : int + Number of contributing cells per ghost node + + Returns + ------- + dtype : np.dtype + + """ + dtype = [("n", int), ("m", int)] + dtype += [(f"j{i}", int) for i in range(numalphaj)] + dtype += [(f"alpha{i}", float) for i in range(numalphaj)] + return np.dtype(dtype) + + +def get_numalphaj(gnc): + """ + Get the number of contributing cells in a ghost node record array + + Parameters + ---------- + gnc : np.recarray + Ghost node data + + Returns + ------- + numalphaj : int + + """ + return len([name for name in gnc.dtype.names if name.startswith("j")]) + + +def _gnc_nodes(gnc): + """Return the cellid column names of a ghost node record array""" + numalphaj = get_numalphaj(gnc) + return ["n", "m"] + [f"j{i}" for i in range(numalphaj)] + + +def _get_ia(ia=None, iac=None): + """Build the zero-based ia array from ia or iac""" + if ia is not None: + return np.asarray(ia, dtype=int) + if iac is None: + return None + return np.concatenate(([0], np.cumsum(np.asarray(iac, dtype=int)))) + + +def _check_gnc(gnc, ia=None, ja=None, iac=None): + """ + Raise if ghost node records are not valid MODFLOW input + + Parameters + ---------- + gnc : np.recarray + Ghost node data + ia : array_like + Zero-based CRS row pointer. Connectivity is not checked if ja is + None or if both ia and iac are None. + ja : array_like + Zero-based CRS column indices + iac : array_like + Number of connections per cell, used if ia is None + + """ + ia = _get_ia(ia, iac) + numalphaj = get_numalphaj(gnc) + alpha = np.zeros(len(gnc)) + for i in range(numalphaj): + alpha += gnc[f"alpha{i}"] + for irec, total in enumerate(alpha): + if total >= 1.0: + raise ValueError( + f"gnc record {irec}: contributing factors sum to {total}, " + "which must be less than one" + ) + + if ia is None or ja is None: + return + ia = np.asarray(ia, dtype=int) + ja = np.asarray(ja, dtype=int) + for irec, rec in enumerate(gnc): + n = rec["n"] + neighbors = ja[ia[n] : ia[n + 1]] + # MODFLOW 6 rejects a ghost node whose n-m connection is absent + if rec["m"] not in neighbors: + raise ValueError( + f"gnc record {irec}: cell {n} is not connected to cell {rec['m']}" + ) + + +def _node_centers(modelgrid): + """Return cell center arrays with one value per node""" + xc = np.asarray(modelgrid.xcellcenters).ravel() + yc = np.asarray(modelgrid.ycellcenters).ravel() + nnodes = modelgrid.nnodes + if xc.shape[0] == nnodes: + return xc, yc + # a vertex grid stores one value per cell2d, repeated for every layer + nlay = nnodes // xc.shape[0] + return np.tile(xc, nlay), np.tile(yc, nlay) + + +def _node_layers(modelgrid): + """Return the layer index of every node""" + ncpl = modelgrid.ncpl + if np.isscalar(ncpl): + ncpl = np.full(modelgrid.nlay, ncpl, dtype=int) + return np.repeat(np.arange(len(ncpl)), ncpl) + + +def _shared_edge_connectivity(modelgrid): + """Build ia and ja from the cells that share an edge in every layer""" + neighbors = modelgrid.neighbors(method="rook") + ncpl = modelgrid.ncpl + if not np.isscalar(ncpl): + if ncpl.min() != ncpl.max(): + raise ValueError( + "Connectivity cannot be built for a grid with a different " + "number of cells in each layer, supply ia or iac and ja" + ) + ncpl = int(ncpl.min()) + nlay = modelgrid.nnodes // ncpl + + # ghost nodes are horizontal, so only the connections within a layer + # are needed and the same layout is repeated for every layer + iac, ja = [], [] + for k in range(nlay): + for icpl in range(ncpl): + conn = sorted(neighbors.get(icpl, [])) + iac.append(len(conn) + 1) + ja.append(k * ncpl + icpl) + ja.extend(k * ncpl + j for j in conn) + return _get_ia(iac=iac), np.array(ja, dtype=int) + + +def _contributing_cells(dn, naxis, d_nm, rtol): + """ + Select the cells that contribute to the ghost node on one connection + + Parameters + ---------- + dn : np.ndarray + Offset from cell n to each of its neighbors, one row per neighbor + naxis : np.ndarray + Connection axis of each neighbor of cell n + d_nm : np.ndarray + Offset from cell n to the connected cell m + rtol : float + Relative tolerance on the transverse offset + + Returns + ------- + sel : np.ndarray or None + Mask of the contributing neighbors, None when the connection needs no + ghost node + alpha : float + Total contributing factor + + """ + axis = int(np.argmax(np.abs(d_nm))) + trans = 1 - axis + offset = d_nm[trans] + # the cell centers line up across the face, so there is nothing to correct + if abs(offset) <= rtol * abs(d_nm[axis]): + return None, 0.0 + + # the contributing cells are the neighbors of n on the side of the offset + sel = (naxis == trans) & (np.sign(dn[:, trans]) == np.sign(offset)) + if not sel.any(): + return None, 0.0 + return sel, abs(offset) / np.abs(dn[sel, trans]).mean() + + +def _cell_areas(modelgrid): + """Return the area of every node from the cell vertices""" + verts = np.asarray(modelgrid.verts) + iverts = modelgrid.iverts + areas = np.empty(len(iverts)) + for i, iv in enumerate(iverts): + iv = [j for j in iv if j is not None] + x, y = verts[iv, 0], verts[iv, 1] + areas[i] = 0.5 * abs(np.dot(x, np.roll(y, -1)) - np.dot(y, np.roll(x, -1))) + if areas.shape[0] != modelgrid.nnodes: + nlay = modelgrid.nnodes // areas.shape[0] + areas = np.tile(areas, nlay) + return areas + + +def get_gnc( + modelgrid, + level=None, + ia=None, + ja=None, + iac=None, + ihc=None, + numalphaj=None, + rtol=1.0e-6, +): + """ + Compute ghost node correction data for a quadtree-like grid + + A ghost node is added in cell n for every horizontal connection to a + finer cell m whose center is offset from the center of n transverse to + the connection. The head at the ghost node is interpolated between cell + n and the neighbors of n on the side of the offset. + + Parameters + ---------- + modelgrid : flopy.discretization.UnstructuredGrid or VertexGrid + Grid the ghost nodes are computed for. Connectivity is taken from + the grid when ia, ja, and iac are None, either from the iac and ja + the grid carries or from the cells that share an edge. + level : array_like + Grid conforming array of refinement levels, where a larger value is a + finer cell. Cell areas are used if None. + ia : array_like + Zero-based CRS row pointer + ja : array_like + Zero-based CRS column indices, with the diagonal first in each row + iac : array_like + Number of connections per cell, used if ia is None + ihc : array_like + Connection type for each entry in ja, where 0 is a vertical + connection. Connections between cells in different layers are + treated as vertical if None. + numalphaj : int + Number of contributing cells written per ghost node. The largest + number found is used if None. Records with fewer contributing cells + repeat a cell rather than pad with zeros, which MODFLOW-USG requires. + rtol : float + Relative tolerance used to decide whether the center of cell m is + offset from the center of cell n (default is 1.0e-6). + + Returns + ------- + gnc : np.recarray + Record array with fields n, m, j0 to j[numalphaj-1], and alpha0 to + alpha[numalphaj-1]. Node numbers are zero-based. + + Notes + ----- + Only horizontal corrections are computed, matching gridgen. + + """ + ia = _get_ia(ia, iac) + if ia is None: + ia = _get_ia(iac=getattr(modelgrid, "iac", None)) + if ja is None: + ja = getattr(modelgrid, "ja", None) + if ia is None or ja is None: + # a vertex grid does not carry connectivity, so build it from the + # cells that share an edge + ia, ja = _shared_edge_connectivity(modelgrid) + ja = np.asarray(ja, dtype=int) + + nnodes = modelgrid.nnodes + xc, yc = _node_centers(modelgrid) + + # a larger size is a coarser cell + if level is None: + size = _cell_areas(modelgrid) + else: + size = -np.asarray(level, dtype=float).ravel() + if size.shape[0] != nnodes: + size = np.tile(size, nnodes // size.shape[0]) + + if ihc is None: + layer = _node_layers(modelgrid) + horizontal = layer[ja] == layer[np.repeat(np.arange(nnodes), np.diff(ia))] + else: + horizontal = np.asarray(ihc) != 0 + + records = [] + for n in range(nnodes): + ipos = np.arange(ia[n] + 1, ia[n + 1]) + ipos = ipos[horizontal[ipos]] + conn = ja[ipos] + if conn.size == 0: + continue + + # column 0 is the x offset and column 1 the y offset to each neighbor + d = np.column_stack((xc[conn] - xc[n], yc[conn] - yc[n])) + axis = np.argmax(np.abs(d), axis=1) + + for k, m in enumerate(conn): + if size[m] >= size[n]: + continue + sel, alpha = _contributing_cells(d, axis, d[k], rtol) + if sel is None: + continue + js = conn[sel] + records.append((n, m, js, alpha / js.size)) + + if numalphaj is None: + numalphaj = max((len(rec[2]) for rec in records), default=1) + dtype = get_gnc_dtype(numalphaj) + + gnc = np.recarray((len(records),), dtype=dtype) + for irec, (n, m, js, alpha) in enumerate(records): + if js.size > numalphaj: + raise ValueError( + f"gnc record {irec}: cell {n} has {js.size} contributing cells, " + f"which is more than numalphaj of {numalphaj}" + ) + # pad by repeating the first contributing cell and splitting its + # factor, which both MODFLOW 6 and MODFLOW-USG accumulate + nrepeat = numalphaj - js.size + 1 + nodes = np.concatenate((np.repeat(js[:1], nrepeat), js[1:])) + alphas = np.concatenate( + (np.full(nrepeat, alpha / nrepeat), np.full(js.size - 1, alpha)) + ) + gnc["n"][irec] = n + gnc["m"][irec] = m + for i in range(numalphaj): + gnc[f"j{i}"][irec] = nodes[i] + gnc[f"alpha{i}"][irec] = alphas[i] + + return gnc + + +def get_gridprops_gnc6( + gnc, dis_type="disv", ncpl=None, ia=None, ja=None, iac=None, check=True +): + """ + Get a dictionary of information needed to create a MODFLOW 6 GNC + Package. The returned dictionary can be unpacked directly into the + ModflowGwfgnc constructor. + + Parameters + ---------- + gnc : np.recarray + Ghost node data with zero-based node numbers + dis_type : str + Discretization the cellids are built for. Valid options are 'disv' + (default) and 'disu'. + ncpl : int + Number of cells per layer, required for 'disv' + ia : array_like + Zero-based CRS row pointer, used to check connectivity + ja : array_like + Zero-based CRS column indices, used to check connectivity + iac : array_like + Number of connections per cell, used if ia is None + check : bool + Verify that each n-m pair is connected and that the contributing + factors sum to less than one (default is True). + + Returns + ------- + gridprops : dict + + Notes + ----- + The correction is applied implicitly unless the explicit option is set, + so the BICGSTAB linear acceleration option should be specified in the IMS + Package. numgnc is zero for a grid without ghost nodes, in which case + the package should not be created. + + """ + if check: + _check_gnc(gnc, ia=ia, ja=ja, iac=iac) + + dis_type = dis_type.lower() + if dis_type == "disv": + if ncpl is None: + raise ValueError("ncpl is required to build disv cellids") + + def cellid(node): + return (node // ncpl, node % ncpl) + elif dis_type == "disu": + + def cellid(node): + return (node,) + else: + raise ValueError(f"Unknown dis_type {dis_type}, expected 'disv' or 'disu'") + + numalphaj = get_numalphaj(gnc) + names = _gnc_nodes(gnc) + gncdata = [ + tuple(cellid(rec[name]) for name in names) + + tuple(rec[f"alpha{i}"] for i in range(numalphaj)) + for rec in gnc + ] + + return { + "numgnc": len(gncdata), + "numalphaj": numalphaj, + "gncdata": gncdata, + } From cef0bdc3db2c2f1ce721fb3e7c2b969c0c8c7c5c Mon Sep 17 00:00:00 2001 From: Joseph Hughes Date: Tue, 4 Aug 2026 16:00:22 -0500 Subject: [PATCH 2/2] fix(gnc): report a node number that is outside the grid _check_gnc dereferenced the node numbers of a record without checking them, so a record built by hand rather than by get_gnc could fail in two ways. A node past the end of ia raised an IndexError from the slice instead of reporting the bad record, and a negative node wrapped silently and checked the connections of a different cell, which let an invalid record pass the check. Both node numbers are now range checked before they are used. --- autotest/test_gnc.py | 10 ++++++++++ flopy/utils/gnc.py | 16 ++++++++++++---- 2 files changed, 22 insertions(+), 4 deletions(-) diff --git a/autotest/test_gnc.py b/autotest/test_gnc.py index dcff6a7585..9b3b051577 100644 --- a/autotest/test_gnc.py +++ b/autotest/test_gnc.py @@ -212,6 +212,16 @@ def test_check_gnc(): _check_gnc(gnc, ia=ia, ja=ja) +@pytest.mark.parametrize("n,m", [(5, 1), (1, 5), (-3, 1), (1, -3)]) +def test_check_gnc_node_out_of_range(n, m): + """A node outside the grid is reported rather than indexed""" + gnc = np.recarray((1,), dtype=get_gnc_dtype(1)) + gnc[0] = (n, m, 0, 0.1) + # a negative node would otherwise wrap and check the wrong cell + with pytest.raises(ValueError, match="which is not a cell of a grid"): + _check_gnc(gnc, ia=np.array([0, 2, 3]), ja=np.array([0, 1, 1, 0])) + + def test_get_gridprops_gnc6_requires_ncpl(): grid, level = synthetic_grid() gnc = get_gnc(grid, level=level) diff --git a/flopy/utils/gnc.py b/flopy/utils/gnc.py index 4ea3cf2d31..9d86b174d0 100644 --- a/flopy/utils/gnc.py +++ b/flopy/utils/gnc.py @@ -95,13 +95,21 @@ def _check_gnc(gnc, ia=None, ja=None, iac=None): return ia = np.asarray(ia, dtype=int) ja = np.asarray(ja, dtype=int) + nodes = ia.shape[0] - 1 for irec, rec in enumerate(gnc): - n = rec["n"] - neighbors = ja[ia[n] : ia[n + 1]] + n, m = rec["n"], rec["m"] + # a node number outside the grid would either index past the end of ia + # or, when it is negative, silently wrap and check the wrong cell + for name, node in (("n", n), ("m", m)): + if not 0 <= node < nodes: + raise ValueError( + f"gnc record {irec}: cell {name} is {node}, which is not a " + f"cell of a grid with {nodes} cells" + ) # MODFLOW 6 rejects a ghost node whose n-m connection is absent - if rec["m"] not in neighbors: + if m not in ja[ia[n] : ia[n + 1]]: raise ValueError( - f"gnc record {irec}: cell {n} is not connected to cell {rec['m']}" + f"gnc record {irec}: cell {n} is not connected to cell {m}" )