DISV model produced by gridgen not converging when original delr/delc contain decimals #1492
Description
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()