DISV model produced by gridgen not converging when original delr/delc contain decimals

[dbrakenhoff]

When the delr and delc of the original structured grid are not whole numbers, the resulting DISV grid produced by gridgen can cause convergence issues when trying to run the MF6 model. I couldn't find out if there were requirements for the original grid for gridgen, but since it a produces a seemingly valid DISV grid, this seems like a bug?

This example showcases the issue, adapted from the example notebook flopy3_gridgen.ipynb. When delr/delc are set to [1.1, 1.2, 1.3, 1.4, 1.6, 1.7, 1.8, 1.9] the model fails to converge, but [1.25, 1.5, 1.75] work fine.

# %%
import sys

import flopy
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
from flopy.utils.gridgen import Gridgen
from shapely.geometry import Polygon

print(sys.version)
print("numpy version: {}".format(np.__version__))
print("matplotlib version: {}".format(mpl.__version__))
print("flopy version: {}".format(flopy.__version__))

# %%
gridgen_ws = "./model/gridgen"
model_ws = "./model"

name = "mymodel"
nlay = 3
nrow = 10
ncol = 10
delr = delc = 1.1  #                                 <-- change delr/dec here, e.g. from 1.0 to 1.1
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=gridgen_ws, 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,
)

# Create and build the gridgen model with a refined area in the middle
g = Gridgen(dis, model_ws=gridgen_ws)
# polys = [Polygon([(4, 4), (6, 4), (6, 6), (4, 6)])]
# g.add_refinement_features(polys, "polygon", 3, range(nlay))
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=model_ws, 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()
sim.run_simulation(silent=True)

# %%
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")
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()

[christianlangevin]

@dbrakenhoff, thanks for this. I'll take a look. I'm guessing that due to rounding errors the code is creating two vertices very close to one another rather than generating a single vertex. Just a hunch. We also noticed that the MODFLOW 6 DISV package will accept a line as a valid cell, and so we are planning to add a trap for that.

[wpbonelli]

This test case sees occasional failures which may be related:

https://github.com/modflowpy/flopy/blob/9c42c373d73e886c9fba6b9a915f51508c002977/autotest/test_grid.py#L881

The number of cells per layer in the generated grid is sometimes off by 3.

[wpbonelli]

With #2076 flopy should no longer produce invalid voronoi cell geometries, though mf6 does not yet have a corresponding trap to my knowledge.

[wpbonelli]

Ignore my prior comments, they refer to unrelated issues. Reproduced in #2119

[wpbonelli]

I'm guessing that due to rounding errors the code is creating two vertices very close to one another rather than generating a single vertex.

Confirmed @langevin-usgs' theory — flopy's gridgen driver should catch these duplicates and remove them