This page was generated from
gridgen_example.py.
It's also available as a notebook.
Creating Layered Quadtree Grids with GRIDGEN
FloPy has a module that can be used to drive the GRIDGEN program. This notebook shows how it works.
Import Modules and Locate Gridgen Executable
[1]:
import os
import sys
from pprint import pformat
from tempfile import TemporaryDirectory
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
# run installed version of flopy or add local path
import flopy
from flopy.utils import flopy_io
from flopy.utils.gridgen import Gridgen
print(sys.version)
print(f"numpy version: {np.__version__}")
print(f"matplotlib version: {mpl.__version__}")
print(f"flopy version: {flopy.__version__}")
3.12.2 | packaged by conda-forge | (main, Feb 16 2024, 20:50:58) [GCC 12.3.0]
numpy version: 1.26.4
matplotlib version: 3.8.4
flopy version: 3.7.0.dev0
The Flopy GRIDGEN module requires that the gridgen executable can be called using subprocess (i.e., gridgen is in your path).
[2]:
gridgen_exe = flopy.which("gridgen")
if gridgen_exe is None:
msg = (
"Warning, gridgen is not in your path. "
"When you create the griden object you will need to "
"provide a full path to the gridgen binary executable."
)
print(msg)
else:
print(
f"gridgen executable was found at: {flopy_io.relpath_safe(gridgen_exe)}"
)
gridgen executable was found at: ../../../../../.local/bin/modflow/gridgen
[3]:
# temporary directory
temp_dir = TemporaryDirectory()
model_ws = temp_dir.name
gridgen_ws = os.path.join(model_ws, "gridgen")
if not os.path.exists(gridgen_ws):
os.makedirs(gridgen_ws, exist_ok=True)
print(f"Model workspace is : {flopy_io.scrub_login(model_ws)}")
print(f"Gridgen workspace is : {flopy_io.scrub_login(gridgen_ws)}")
Model workspace is : /tmp/tmplpq06fke
Gridgen workspace is : /tmp/tmplpq06fke/gridgen
Basic Gridgen Operations
Setup Base MODFLOW Grid
GRIDGEN works off of a base MODFLOW grid. The following information defines the basegrid.
[4]:
Lx = 100.0
Ly = 100.0
nlay = 2
nrow = 51
ncol = 51
delr = Lx / ncol
delc = Ly / nrow
h0 = 10
h1 = 5
top = h0
botm = np.zeros((nlay, nrow, ncol), dtype=np.float32)
botm[1, :, :] = -10.0
[5]:
ms = flopy.modflow.Modflow(rotation=-20.0)
dis = flopy.modflow.ModflowDis(
ms,
nlay=nlay,
nrow=nrow,
ncol=ncol,
delr=delr,
delc=delc,
top=top,
botm=botm,
)
Create the Gridgen Object
[6]:
g = Gridgen(ms.modelgrid, model_ws=gridgen_ws)
Add an Optional Active Domain
Cells outside of the active domain will be clipped and not numbered as part of the final grid. If this step is not performed, then all cells will be included in the final grid.
[7]:
# setup the active domain
adshp = os.path.join(gridgen_ws, "ad0")
adpoly = [[[(0, 0), (0, 60), (40, 80), (60, 0), (0, 0)]]]
# g.add_active_domain(adpoly, range(nlay))
Refine the Grid
[8]:
x = Lx * np.random.random(10)
y = Ly * np.random.random(10)
wells = list(zip(x, y))
g.add_refinement_features(wells, "point", 3, range(nlay))
rf0shp = os.path.join(gridgen_ws, "rf0")
[9]:
river = [[(-20, 10), (60, 60)]]
g.add_refinement_features(river, "line", 3, range(nlay))
rf1shp = os.path.join(gridgen_ws, "rf1")
[10]:
g.add_refinement_features(adpoly, "polygon", 1, range(nlay))
rf2shp = os.path.join(gridgen_ws, "rf2")
Plot the Gridgen Input
[11]:
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(1, 1, 1, aspect="equal")
mm = flopy.plot.PlotMapView(model=ms)
mm.plot_grid()
flopy.plot.plot_shapefile(rf2shp, ax=ax, facecolor="yellow", edgecolor="none")
flopy.plot.plot_shapefile(rf1shp, ax=ax, linewidth=10)
flopy.plot.plot_shapefile(rf0shp, ax=ax, facecolor="red", radius=1)
[11]:
<matplotlib.collections.PatchCollection at 0x7fbb20a99ee0>
Build the Grid
[12]:
g.build(verbose=False)
Plot the Grid
[13]:
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(1, 1, 1, aspect="equal")
g.plot(ax, linewidth=0.5)
flopy.plot.plot_shapefile(
rf2shp, ax=ax, facecolor="yellow", edgecolor="none", alpha=0.2
)
flopy.plot.plot_shapefile(rf1shp, ax=ax, linewidth=10, alpha=0.2)
flopy.plot.plot_shapefile(rf0shp, ax=ax, facecolor="red", radius=1, alpha=0.2)
[13]:
<matplotlib.collections.PatchCollection at 0x7fbb20a9a2a0>
Create a Flopy ModflowDisu Object
[14]:
mu = flopy.mfusg.MfUsg(model_ws=gridgen_ws, modelname="mfusg")
disu = g.get_disu(mu)
disu.write_file()
# print(disu)
Intersect Features with the Grid
[15]:
adpoly_intersect = g.intersect(adpoly, "polygon", 0)
print(adpoly_intersect.dtype.names)
print(adpoly_intersect)
print(adpoly_intersect.nodenumber)
('nodenumber', 'polyid', 'totalarea', 'SHAPEID')
[( 322, 0, 0.961169 , 0) ( 382, 0, 0.961169 , 0) ( 325, 0, 0.961169 , 0)
... (5991, 0, 0.0977309, 0) (5993, 0, 0.348072 , 0)
(5994, 0, 0.0428516, 0)]
[ 322 382 325 ... 5991 5993 5994]
[16]:
well_intersect = g.intersect(wells, "point", 0)
print(well_intersect.dtype.names)
print(well_intersect)
print(well_intersect.nodenumber)
('nodenumber', 'pointid', 'SHAPEID')
[(2533, 0, 0) (4759, 1, 1) (1206, 2, 2) (2703, 3, 3) (3987, 4, 4)
(2794, 5, 5) (2759, 6, 6) (3321, 7, 7) (1855, 8, 8)]
[2533 4759 1206 2703 3987 2794 2759 3321 1855]
[17]:
river_intersect = g.intersect(river, "line", 0)
print(river_intersect.dtype.names)
# print(river_intersect)
# print(river_intersect.nodenumber)
('nodenumber', 'arcid', 'length', 'starting_distance', 'ending_distance', 'SHAPEID')
Plot Intersected Features
[18]:
a = np.zeros((g.nodes), dtype=int)
a[adpoly_intersect.nodenumber] = 1
a[well_intersect.nodenumber] = 2
a[river_intersect.nodenumber] = 3
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(1, 1, 1, aspect="equal")
g.plot(ax, a=a, masked_values=[0], edgecolor="none", cmap="jet")
flopy.plot.plot_shapefile(rf2shp, ax=ax, facecolor="yellow", alpha=0.25)
[18]:
<matplotlib.collections.PatchCollection at 0x7fbb0ee38c80>
Use Gridgen to Build MODFLOW 6 DISV Model
In this section, we will reproduce the MODFLOW 6 Quick Start example that is shown on the main page of the flopy repository (https://github.com/modflowpy/flopy).
A main idea for DISV in MODFLOW 6 is that each layer much have the same spatial grid. Gridgen allows the creation of a grid that has a different number of cells within each layer. This type of grid cannot be used with the DISV Package. To make sure that the resulting grid is the same for each layer, refinement should be added to all layers when using the flopy Gridgen wrapper.
[19]:
from shapely.geometry import Polygon
name = "dummy"
nlay = 3
nrow = 10
ncol = 10
delr = delc = 1.0
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(gwf.modelgrid, 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()
[20]:
# Create and plot a flopy VertexGrid object
# Note that this is not necessary, because it can
# be created automatically after the disv package
# is added to the gwf model (gwf.modelgrid should exist)
gridprops_vg = g.get_gridprops_vertexgrid()
vgrid = flopy.discretization.VertexGrid(**gridprops_vg)
vgrid.plot()
[20]:
<matplotlib.collections.LineCollection at 0x7fbb0eefb5c0>
[21]:
# retrieve a dictionary of arguments to be passed
# directly into the flopy disv constructor
disv_gridprops = g.get_gridprops_disv()
disv_gridprops.keys()
[21]:
dict_keys(['nlay', 'ncpl', 'top', 'botm', 'nvert', 'vertices', 'cell2d'])
[22]:
# 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])
chdspd
[22]:
[[(0, 0), 1.0], [(0, 435), 0.0]]
[23]:
# build uun and post-process the MODFLOW 6 model
ws = os.path.join(model_ws, "gridgen_disv")
name = "mymodel"
sim = flopy.mf6.MFSimulation(
sim_name=name, sim_ws=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 = f"{name}.bud"
head_file = f"{name}.hds"
oc = flopy.mf6.ModflowGwfoc(
gwf,
budget_filerecord=budget_file,
head_filerecord=head_file,
saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
)
sim.write_simulation()
success, buff = sim.run_simulation(silent=True)
assert success, pformat(buff)
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")
[23]:
<matplotlib.quiver.Quiver at 0x7fbb0e068140>
Use Gridgen to Build MODFLOW 6 DISU Model
In this section, we will reproduce the MODFLOW 6 Quick Start example that is shown on the main page of the flopy repository (https://github.com/modflowpy/flopy).
DISU is the most general grid form that can be used with MODFLOW 6. It does not have the requirement that each layer must use the same grid
[24]:
from shapely.geometry import Polygon
name = "dummy"
nlay = 3
nrow = 10
ncol = 10
delr = delc = 1.0
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(gwf.modelgrid, model_ws=gridgen_ws)
polys = [Polygon([(4, 4), (6, 4), (6, 6), (4, 6)])]
g.add_refinement_features(polys, "polygon", 3, layers=[0])
g.build()
[25]:
# retrieve a dictionary of arguments to be passed
# directly into the flopy disu constructor
disu_gridprops = g.get_gridprops_disu6()
disu_gridprops.keys()
[25]:
dict_keys(['nodes', 'top', 'bot', 'area', 'iac', 'nja', 'ja', 'cl12', 'ihc', 'hwva', 'angldegx', 'nvert', 'vertices', 'cell2d'])
[26]:
# Create and plot a flopy UnstructuredGrid object
gridprops_ug = g.get_gridprops_unstructuredgrid()
ugrid = flopy.discretization.UnstructuredGrid(**gridprops_ug)
f = plt.figure(figsize=(10, 10))
for ilay in range(g.nlay):
ax = plt.subplot(1, g.nlay, ilay + 1)
ugrid.plot(layer=ilay, ax=ax)
[27]:
# find the cell numbers for constant heads
chdspd = []
for x, y, head in [(0, 10, 1.0), (10, 0, 0.0)]:
ra = g.intersect([(x, y)], "point", 0)
ic = ra["nodenumber"][0]
chdspd.append([(ic,), head])
chdspd
[27]:
[[(0,), 1.0], [(435,), 0.0]]
[28]:
# build run and post-process the MODFLOW 6 model
ws = os.path.join(model_ws, "gridgen_disu")
name = "mymodel"
sim = flopy.mf6.MFSimulation(
sim_name=name, sim_ws=ws, exe_name="mf6", verbosity_level=1
)
tdis = flopy.mf6.ModflowTdis(sim)
ims = flopy.mf6.ModflowIms(sim, linear_acceleration="bicgstab")
gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True)
disu = flopy.mf6.ModflowGwfdisu(gwf, **disu_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 = f"{name}.bud"
head_file = f"{name}.hds"
oc = flopy.mf6.ModflowGwfoc(
gwf,
budget_filerecord=budget_file,
head_filerecord=head_file,
saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
)
sim.write_simulation()
success, buff = sim.run_simulation(silent=True)
assert success, pformat(buff)
head = gwf.output.head().get_data()
bud = gwf.output.budget()
spdis = bud.get_data(text="DATA-SPDIS")[0]
gwf.modelgrid.set_coord_info(angrot=15)
f = plt.figure(figsize=(10, 10))
vmin = head.min()
vmax = head.max()
for ilay in range(gwf.modelgrid.nlay):
ax = plt.subplot(1, g.nlay, ilay + 1)
pmv = flopy.plot.PlotMapView(gwf, layer=ilay, ax=ax)
ax.set_aspect("equal")
pmv.plot_array(head.flatten(), cmap="jet", vmin=vmin, vmax=vmax)
pmv.plot_grid(colors="k", alpha=0.1)
pmv.contour_array(
head, levels=[0.2, 0.4, 0.6, 0.8], linewidths=3.0, vmin=vmin, vmax=vmax
)
ax.set_title(f"Layer {ilay + 1}")
pmv.plot_vector(spdis["qx"], spdis["qy"], color="white")
writing simulation...
writing simulation name file...
writing simulation tdis package...
writing solution package ims_-1...
writing model mymodel...
writing model name file...
writing package disu...
writing package ic...
writing package npf...
writing package chd_0...
INFORMATION: maxbound in ('gwf6', 'chd', 'dimensions') changed to 2 based on size of stress_period_data
writing package oc...
Use Gridgen to Build MODFLOW-USG DISU Model
In this section, we will reproduce the MODFLOW 6 Quick Start example that is shown on the main page of the flopy repository (https://github.com/modflowpy/flopy).
In this last example, MODFLOW-USG will be used to simulate the problem.
[29]:
# build run and post-process the MODFLOW 6 model
ws = os.path.join(model_ws, "gridgen_mfusg")
name = "mymodel"
chdspd = []
for x, y, head in [(0, 10, 1.0), (10, 0, 0.0)]:
ra = g.intersect([(x, y)], "point", 0)
ic = ra["nodenumber"][0]
chdspd.append([ic, head, head])
gridprops = g.get_gridprops_disu5()
# create the mfusg modoel
m = flopy.mfusg.MfUsg(
modelname=name,
model_ws=ws,
version="mfusg",
exe_name="mfusg",
structured=False,
)
disu = flopy.mfusg.MfUsgDisU(m, **gridprops)
bas = flopy.modflow.ModflowBas(m)
lpf = flopy.mfusg.MfUsgLpf(m)
chd = flopy.modflow.ModflowChd(m, stress_period_data=chdspd)
sms = flopy.mfusg.MfUsgSms(m)
oc = flopy.modflow.ModflowOc(m)
m.write_input()
success, buff = m.run_model(silent=True, report=True)
assert success, pformat(buff)
# head is returned as a list of head arrays for each layer
head_file = os.path.join(ws, f"{name}.hds")
head = flopy.utils.HeadUFile(head_file).get_data()
# MODFLOW-USG does not have vertices, so we need to create
# and unstructured grid and then assign it to the model. This
# will allow plotting and other features to work properly.
gridprops_ug = g.get_gridprops_unstructuredgrid()
ugrid = flopy.discretization.UnstructuredGrid(**gridprops_ug)
m.modelgrid = ugrid
f = plt.figure(figsize=(10, 10))
vmin = 0.0
vmax = 1.0
for ilay in range(disu.nlay):
ax = plt.subplot(1, g.nlay, ilay + 1)
pmv = flopy.plot.PlotMapView(m, layer=ilay, ax=ax)
ax.set_aspect("equal")
pmv.plot_array(head[ilay], cmap="jet", vmin=vmin, vmax=vmax)
pmv.plot_grid(colors="k", alpha=0.1)
pmv.contour_array(head[ilay], levels=[0.2, 0.4, 0.6, 0.8], linewidths=3.0)
ax.set_title(f"Layer {ilay + 1}")
[30]:
try:
# ignore PermissionError on Windows
temp_dir.cleanup()
except:
pass