diff --git a/autotest/test_gridgen.py b/autotest/test_gridgen.py index eb1c53358..521835c4e 100644 --- a/autotest/test_gridgen.py +++ b/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() diff --git a/flopy/utils/cvfdutil.py b/flopy/utils/cvfdutil.py index 4491f71f7..3eb33d6e8 100644 --- a/flopy/utils/cvfdutil.py +++ b/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: diff --git a/flopy/utils/gridgen.py b/flopy/utils/gridgen.py index 2aa1eea18..9fb46f4fc 100644 --- a/flopy/utils/gridgen.py +++ b/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(