82 lines
2.3 KiB
Python
82 lines
2.3 KiB
Python
import numpy as np
|
|
import flopy
|
|
|
|
# Assign name and create modflow model object
|
|
modelname = "tutorial1"
|
|
mf = flopy.modflow.Modflow(modelname, exe_name="mf2005")
|
|
|
|
# Model domain and grid definition
|
|
Lx = 1000.0
|
|
Ly = 1000.0
|
|
ztop = 0.0
|
|
zbot = -50.0
|
|
nlay = 1
|
|
nrow = 10
|
|
ncol = 10
|
|
delr = Lx / ncol
|
|
delc = Ly / nrow
|
|
delv = (ztop - zbot) / nlay
|
|
botm = np.linspace(ztop, zbot, nlay + 1)
|
|
|
|
# Create the discretization object
|
|
dis = flopy.modflow.ModflowDis(
|
|
mf, nlay, nrow, ncol, delr=delr, delc=delc, top=ztop, botm=botm[1:]
|
|
)
|
|
|
|
# Variables for the BAS package
|
|
ibound = np.ones((nlay, nrow, ncol), dtype=np.int32)
|
|
ibound[:, :, 0] = -1
|
|
ibound[:, :, -1] = -1
|
|
strt = np.ones((nlay, nrow, ncol), dtype=np.float32)
|
|
strt[:, :, 0] = 10.0
|
|
strt[:, :, -1] = 0.0
|
|
bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt)
|
|
|
|
# Add LPF package to the MODFLOW model
|
|
lpf = flopy.modflow.ModflowLpf(mf, hk=10.0, vka=10.0, ipakcb=53)
|
|
|
|
# Add OC package to the MODFLOW model
|
|
spd = {(0, 0): ["print head", "print budget", "save head", "save budget"]}
|
|
oc = flopy.modflow.ModflowOc(mf, stress_period_data=spd, compact=True)
|
|
|
|
# Add PCG package to the MODFLOW model
|
|
pcg = flopy.modflow.ModflowPcg(mf)
|
|
|
|
# Write the MODFLOW model input files
|
|
mf.write_input()
|
|
|
|
# Run the MODFLOW model
|
|
success, buff = mf.run_model()
|
|
|
|
# Post process the results
|
|
import matplotlib.pyplot as plt
|
|
import flopy.utils.binaryfile as bf
|
|
|
|
plt.subplot(1, 1, 1, aspect="equal")
|
|
hds = bf.HeadFile(modelname + ".hds")
|
|
head = hds.get_data(totim=1.0)
|
|
levels = np.arange(1, 10, 1)
|
|
extent = (delr / 2.0, Lx - delr / 2.0, Ly - delc / 2.0, delc / 2.0)
|
|
plt.contour(head[0, :, :], levels=levels, extent=extent)
|
|
plt.savefig("tutorial1a.png")
|
|
|
|
fig = plt.figure(figsize=(10, 10))
|
|
ax = fig.add_subplot(1, 1, 1, aspect="equal")
|
|
|
|
hds = bf.HeadFile(modelname + ".hds")
|
|
times = hds.get_times()
|
|
head = hds.get_data(totim=times[-1])
|
|
levels = np.linspace(0, 10, 11)
|
|
|
|
cbb = bf.CellBudgetFile(modelname + ".cbc")
|
|
kstpkper_list = cbb.get_kstpkper()
|
|
frf = cbb.get_data(text="FLOW RIGHT FACE", totim=times[-1])[0]
|
|
fff = cbb.get_data(text="FLOW FRONT FACE", totim=times[-1])[0]
|
|
|
|
modelmap = flopy.plot.ModelMap(model=mf, layer=0)
|
|
qm = modelmap.plot_ibound()
|
|
lc = modelmap.plot_grid()
|
|
cs = modelmap.contour_array(head, levels=levels)
|
|
quiver = modelmap.plot_discharge(frf, fff, head=head)
|
|
plt.savefig("tutorial1b.png")
|