<reponame>andrewcalderwood/flopy<filename>.docs/pysrc/tutorial1.py
import os
import sys
import numpy as np
import flopy

# Assign name and create modflow model object
modelname = "tutorial1"
exe_name = os.path.join("..", ".bin", "mf2005")
mf = flopy.modflow.Modflow(modelname, exe_name=exe_name)

# 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
mf.run_model()

import matplotlib.pyplot as plt
import flopy.utils.binaryfile as bf

plt.subplot(1, 1, 1, aspect="equal")
hds = bf.HeadFile(f"{modelname}.hds")
head = hds.get_data(totim=1.0)
levels = np.arange(1, 10, 1)
extent = (delr / 2.0, Lx - delc / 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(f"{modelname}.hds")
times = hds.get_times()
head = hds.get_data(totim=times[-1])
levels = np.linspace(0, 10, 11)

cbb = bf.CellBudgetFile(f"{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]

pmv = flopy.plot.PlotMapView(model=mf, layer=0)
qm = pmv.plot_ibound()
lc = pmv.plot_grid()
cs = pmv.contour_array(head, levels=levels)
quiver = pmv.plot_discharge(frf, fff, head=head)
plt.savefig("tutorial1b.png")
