fix(gridgen): trap duplicate vertices #1492

autotest/test_gridgen.py

@@ -4,6 +4,7 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
+import pandas as pd
 import pytest
 from matplotlib.collections import LineCollection, PathCollection, QuadMesh
 from modflow_devtools.markers import requires_exe, requires_pkg
@@ -851,3 +852,103 @@ def test_gridgen(function_tmpdir):
         "should not be (without vertical pass through activated)."
     )
     assert max(ja0) <= disu_vp.nodelay[0], msg
+
+
+@requires_exe("mf6", "gridgen")
+def test_flopy_issue_1492(function_tmpdir):
+    """
+    Submitted by David Brakenhoff in
+    https://github.com/modflowpy/flopy/issues/1492
+    """
+
+    name = "issue1492"
+    nlay = 3
+    nrow = 10
+    ncol = 10
+    delr = delc = 1.1  # <-- 1.0 converges
+    top = 1
+    bot = 0
+    dz = (top - bot) / nlay
+    botm = [top - k * dz for k in range(1, nlay + 1)]
+
+    # Create a dummy model and regular grid to use as the base grid for gridgen
+    sim = flopy.mf6.MFSimulation(
+        sim_name=name, sim_ws=function_tmpdir, exe_name="mf6"
+    )
+    gwf = flopy.mf6.ModflowGwf(sim, modelname=name)
+    dis = flopy.mf6.ModflowGwfdis(
+        gwf,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=delr,
+        delc=delc,
+        top=top,
+        botm=botm,
+    )
+    og_grid = gwf.modelgrid
+
+    # Create and build the gridgen model
+    g = Gridgen(dis, model_ws=function_tmpdir)
+    g.build()
+
+    # retrieve a dictionary of arguments to be passed
+    # directly into the flopy disv constructor
+    disv_gridprops = g.get_gridprops_disv()
+
+    # find the cell numbers for constant heads
+    chdspd = []
+    ilay = 0
+    for x, y, head in [(0, 10, 1.0), (10, 0, 0.0)]:
+        ra = g.intersect([(x, y)], "point", ilay)
+        ic = ra["nodenumber"][0]
+        chdspd.append([(ilay, ic), head])
+
+    # build run and post-process the MODFLOW 6 model
+    sim = flopy.mf6.MFSimulation(
+        sim_name=name,
+        sim_ws=function_tmpdir,
+        exe_name="mf6",
+        verbosity_level=0,
+    )
+    tdis = flopy.mf6.ModflowTdis(sim)
+    ims = flopy.mf6.ModflowIms(sim, linear_acceleration="bicgstab")
+    gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True)
+    disv = flopy.mf6.ModflowGwfdisv(gwf, **disv_gridprops)
+    ic = flopy.mf6.ModflowGwfic(gwf)
+    npf = flopy.mf6.ModflowGwfnpf(
+        gwf, xt3doptions=True, save_specific_discharge=True
+    )
+    chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chdspd)
+    budget_file = name + ".bud"
+    head_file = name + ".hds"
+    oc = flopy.mf6.ModflowGwfoc(
+        gwf,
+        budget_filerecord=budget_file,
+        head_filerecord=head_file,
+        saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
+    )
+
+    sim.write_simulation()
+    success, _ = sim.run_simulation(silent=False)
+    assert success
+
+    # debugging duplicate vertices
+    grid = gwf.modelgrid
+    og_verts = pd.DataFrame(og_grid.verts, columns=["x", "y"])
+    mg_verts = pd.DataFrame(grid.verts, columns=["x", "y"])
+
+    plot_debug = False
+    if plot_debug:
+        head = gwf.output.head().get_data()
+        bud = gwf.output.budget()
+        spdis = bud.get_data(text="DATA-SPDIS")[0]
+        pmv = flopy.plot.PlotMapView(gwf)
+        pmv.plot_array(head)
+        pmv.plot_grid(colors="white")
+        ax = plt.gca()
+        verts = grid.verts
+        ax.plot(verts[:, 0], verts[:, 1], "bo", alpha=0.25, ms=5)
+        pmv.contour_array(head, levels=[0.2, 0.4, 0.6, 0.8], linewidths=3.0)
+        pmv.plot_vector(spdis["qx"], spdis["qy"], color="white")
+        plt.show()

flopy/utils/cvfdutil.py

@@ -1,4 +1,5 @@
 import numpy as np
+import pandas as pd
 
 from .utl_import import import_optional_dependency
 
@@ -115,6 +116,7 @@ def to_cvfd(
     nodestart=None,
     nodestop=None,
     skip_hanging_node_check=False,
+    duplicate_decimals=9,
     verbose=False,
 ):
     """
@@ -135,6 +137,11 @@ def to_cvfd(
         skip the hanging node check.  this may only be necessary for quad-based
         grid refinement. (default is False)
 
+    duplicate_decimals : int
+        decimals to round duplicate vertex checks.  GRIDGEN can occasionally
+        produce very-nearly overlapping vertices, this can be used to change
+        the sensitivity for filtering out duplicates. (default is 9)
+
     verbose : bool
         print messages to the screen. (default is False)
 
@@ -176,7 +183,10 @@ def to_cvfd(
         xcyc[icell, 1] = yc
         ivertlist = []
         for p in points:
-            pt = tuple(p)
+            pt = (
+                round(p[0], duplicate_decimals),
+                round(p[1], duplicate_decimals),
+            )
             if pt in vertexdict:
                 ivert = vertexdict[pt]
             else:

flopy/utils/gridgen.py

@@ -1892,7 +1892,6 @@ def _mkvertdict(self):
         for i in range(len(shapes)):
             nodenumber = int(records[i][idx]) - 1
             self._vertdict[nodenumber] = shapes[i].points
-        return
 
     @staticmethod
     def read_qtg_nod(