In Problem P5.1 from pages 244-245 in Anderson, Woessner and Hunt (2015), we are asked to construct an areal 2D model to assess impacts of grid sensitivity on pumping. Develop a 2D areal model using Flopy to solve for heads in the upper right-hand quadrant of the island shown in Fig. P5.1. The aquifer is confined, homogeneous, and isotropic with transmissivity, T, equal to 10,000 ft2/day. Recharge, R, occurs uniformly through a leaky confining bed at a rate of 0.00305 ft/day. The half-width of the island, l, is 12,000 ft. The head at the perimeter of the island is at sea level (use h = 0 ft). The heads are symmetric across the groundwater divides that separate the island into four quadrants (Fig. P5.1). Use a point-centered FD grid so that the node at the observation well in the center of the island (Fig. P5.1) is located directly on the groundwater divides that form the lefthand side and lower boundary of the quadrant model.
Include a water budget in your model. The inflow to the model is the volume of water entering from recharge. In this notebook, we will work through the problem again using MODFLOW and the Python tool set Flopy. Notice how much code is reused from previous examples. Note also, P5.1 gives directions for calculating the waterbudget; because in this exercise we are using MODFLOW it handles all water budget calculations for us.
Write the mathematical model for this problem including the governing equation and boundary conditions of the quadrant model.
Solve the model using an error tolerance equal to 1E-4 ft and test two designs for the nodal network: (1) a 4 x 7 array of nodes (delta x = delta y = 4000 ft); (2) a 13 x 25 array of nodes (delta x = delta y= 1000 ft). For each nodal network, print out the head solution to the fourth decimal place.
Below is an iPython Notebook that builds a Python MODFLOW model for this problem and plots results. See the Github wiki associated with this Chapter for information on one suggested installation and setup configuration for Python and iPython Notebook.
[Acknowledgements: This tutorial was created by Randy Hunt and all failings are mine. The exercise here has benefited greatly from the online Flopy tutorial and example notebooks developed by Chris Langevin and Joe Hughes for the USGS Spring 2015 Python Training course GW1774]
In this example, we will create a simple groundwater flow modelusing the Flopy website approach. Visit the tutorial website here.
In [1]:
%matplotlib inline
import sys
import os
import shutil
import numpy as np
from subprocess import check_output
# Import flopy
import flopy
For this tutorial, we will work in a new subdirectory underneath the directory where the notebook is located. We can use some fancy Python tools to help us manage the directory creation. Note that if you encounter path problems with this workbook, you can stop and then restart the kernel and the paths will be reset.
In [2]:
# Set the name of the path to the model working directory
dirname = "P5-1_Island_recharge"
datapath = os.getcwd()
modelpath = os.path.join(datapath, dirname)
print 'Name of model path: ', modelpath
# Now let's check if this directory exists. If not, then we will create it.
if os.path.exists(modelpath):
print 'Model working directory already exists.'
else:
print 'Creating model working directory.'
os.mkdir(modelpath)
In [3]:
# model domain and grid definition
# for clarity, user entered variables are all caps; python syntax are lower case or mixed case
# This is an areal 2D model that uses island symmetry to reduce the grid size.
LX = 16000. # half width of island + one node for constant head boundary condition
LY = 28000. # half height of island + one node for constant head boundary condition
ZTOP = 0. # the system is confined
ZBOT = -50.
NLAY = 1
NROW = 7
NCOL = 4
DELR = LX / NCOL # recall that MODFLOW convention is DELR is along a row, thus has items = NCOL; see page XXX in AW&H (2015)
DELC = LY / NROW # recall that MODFLOW convention is DELC is along a column, thus has items = NROW; see page XXX in AW&H (2015)
DELV = (ZTOP - ZBOT) / NLAY
BOTM = np.linspace(ZTOP, ZBOT, NLAY + 1)
HK = 200.
VKA = 1.
RCH = 0.00305
# WELLQ = 0. #not needed for Problem 5.1
print "DELR =", DELR, " DELC =", DELC, ' DELV =', DELV
print "BOTM =", BOTM
print "Recharge =", RCH
#print "Pumping well rate =", WELLQ
Create a flopy MODFLOW object: flopy.modflow.Modflow.
In [4]:
# Assign name and create modflow model object
modelname = 'P5-1'
#exe_name = os.path.join(datapath, 'mf2005.exe') # for Windows OS
exe_name = os.path.join(datapath, 'mf2005') # for Mac OS
print 'Model executable: ', exe_name
MF = flopy.modflow.Modflow(modelname, exe_name=exe_name, model_ws=modelpath)
Create a flopy discretization package object: flopy.modflow.ModflowDis.
In [5]:
# Create the discretization object
TOP = ZTOP * np.ones((NROW, NCOL),dtype=np.float)
In [6]:
DIS_PACKAGE = flopy.modflow.ModflowDis(MF, NLAY, NROW, NCOL, delr=DELR, delc=DELC,
top=TOP, botm=BOTM[1:], laycbd=0)
# print DIS_PACKAGE #uncomment this on far left to see information about the flopy object
Create a flopy basic package object: flopy.modflow.ModflowBas.
In [7]:
# Variables for the BAS package
IBOUND = np.ones((NLAY, NROW, NCOL), dtype=np.int32) # all nodes are active (IBOUND = 1)
# make the top of the profile specified head by setting the IBOUND = -1
IBOUND[:, 0, :] = -1 #don't forget arrays are zero-based! Sets first row
IBOUND[:, :, 0] = -1 # Sets first column
print IBOUND
In [8]:
STRT = 1 * np.ones((NLAY, NROW, NCOL), dtype=np.float32) # set starting head to 1` through out model domain
STRT[:, 0, :] = 0. # top row ocean elevation for setting constant head
STRT[:, :, 0] = 0. # first column ocean elevation for setting constant head
print STRT
In [9]:
BAS_PACKAGE = flopy.modflow.ModflowBas(MF, ibound=IBOUND, strt=STRT)
# print BAS_PACKAGE # uncomment this at far left to see the information about the flopy BAS object
Create a flopy layer property flow package object: flopy.modflow.ModflowLpf.
In [10]:
LPF_PACKAGE = flopy.modflow.ModflowLpf(MF, laytyp=1, hk=HK, vka=VKA) # we defined the K and anisotropy at top of file
# print LPF_PACKAGE # uncomment this at far left to see the information about the flopy LPF object
In [11]:
#WEL_PACKAGE = flopy.modflow.ModflowWel(MF, stress_period_data=[0,0,0,WELLQ]) # remember python 0 index, layer 0 = layer 1 in MF
#print WEL_PACKAGE # uncomment this at far left to see the information about the flopy WEL object
Create a flopy output control object: flopy.modflow.ModflowOc.
In [12]:
OC_PACKAGE = flopy.modflow.ModflowOc(MF) # we'll use the defaults for the model output
# print OC_PACKAGE # uncomment this at far left to see the information about the flopy OC object
Create a flopy pcg package object: flopy.modflow.ModflowPcg.
In [13]:
PCG_PACKAGE = flopy.modflow.ModflowPcg(MF, mxiter=500, iter1=100, hclose=1e-04, rclose=1e-03, relax=0.98, damp=0.5)
# print PCG_PACKAGE # uncomment this at far left to see the information about the flopy PCG object
Create a flopy pcg package object: flopy.modflow.ModflowRch.
In [14]:
RCH_PACKAGE = flopy.modflow.ModflowRch(MF, rech=RCH)
# print RCH_PACKAGE # uncomment this at far left to see the information about the flopy RCH object
In [15]:
#Before writing input, destroy all files in folder to prevent reusing old files
#Here's the working directory
print modelpath
#Here's what's currently in the working directory
modelfiles = os.listdir(modelpath)
print modelfiles
In [16]:
#delete these files to prevent us from reading old results
modelfiles = os.listdir(modelpath)
for filename in modelfiles:
f = os.path.join(modelpath, filename)
if modelname in f:
try:
os.remove(f)
print 'Deleted: ', filename
except:
print 'Unable to delete: ', filename
In [17]:
#Now write the model input files
MF.write_input()
The model datasets are written using a single command (mf.write_input).
Check in the model working directory and verify that the input files have been created. Or if you might just add another cell, right after this one, that prints a list of all the files in our model directory. The path we are working in is returned from this next block.
In [18]:
# return current working directory
print "You can check the newly created files in", modelpath
In [19]:
silent = False #Print model output to screen?
pause = False #Require user to hit enter? Doesn't mean much in Ipython notebook
report = True #Store the output from the model in buff
success, buff = MF.run_model(silent=silent, pause=pause, report=report)
To read heads from the MODFLOW binary output file, we can use the flopy.utils.binaryfile module. Specifically, we can use the HeadFile object from that module to extract head data arrays.
In [20]:
#imports for plotting and reading the MODFLOW binary output file
import matplotlib.pyplot as plt
import flopy.utils.binaryfile as bf
#Create the headfile object and grab the results for last time.
headfile = os.path.join(modelpath, modelname + '.hds')
headfileobj = bf.HeadFile(headfile)
#Get a list of times that are contained in the model
times = headfileobj.get_times()
print 'Headfile (' + modelname + '.hds' + ') contains the following list of times: ', times
In [21]:
#Get a numpy array of heads for totim = 1.0
#The get_data method will extract head data from the binary file.
HEAD = headfileobj.get_data(totim=1.0)
#Print statistics on the head
print 'Head statistics'
print ' min: ', HEAD.min()
print ' max: ', HEAD.max()
print ' std: ', HEAD.std()
In [22]:
#Create a contour plot of heads
FIG = plt.figure(figsize=(12,10))
#setup contour levels and plot extent
LEVELS = np.arange(0., 26., 5.)
EXTENT = (DELR/2., LX - DELR/2., DELC/2., LY - DELC/2.)
print 'Contour Levels: ', LEVELS
print 'Extent of domain: ', EXTENT
#Make a contour plot on the first axis
AX1 = FIG.add_subplot(1, 2, 1, aspect='equal')
AX1.set_xlabel("x")
AX1.set_ylabel("y")
YTICKS = np.arange(0, 28000, 4000)
AX1.set_yticks(YTICKS)
AX1.set_title("P5.1 Island Recharge Problem")
AX1.text(2250, 10000, r"side ocean boundary condition", fontsize=10, color="blue", rotation='vertical')
AX1.text(5500, 25000, r"top ocean boundary condition", fontsize=10, color="blue")
AX1.contour(np.flipud(HEAD[0, :, :]), levels=LEVELS, extent=EXTENT)
#Make a color flood on the second axis
AX2 = FIG.add_subplot(1, 2, 2, aspect='equal')
AX2.set_xlabel("x")
AX2.set_ylabel("y")
AX2.set_yticks(YTICKS)
AX2.set_title("P5.1 color flood")
AX2.text(5500, 25000, r"top ocean boundary condition", fontsize=10, color="white")
AX2.text(2500, 10500, r"side ocean boundary condition", fontsize=10, color="white", rotation='vertical')
cax = AX2.imshow(HEAD[0, :, :], extent=EXTENT, interpolation='nearest')
cbar = FIG.colorbar(cax, orientation='vertical', shrink=0.45)
In [23]:
print HEAD
Look at the bottom of the MODFLOW output file (ending with a *.list) and write down the water balance reported to compare to the less coarse grid size done in the next part of the problem.
In [24]:
LX = 13000. # same as before - half width of island + one node for constant head boundary condition
LY = 25000. # same as before - half height of island + one node for constant head boundary condition
NLAY = 1
NROW = 25
NCOL = 13
DELR = LX / NCOL # recall that MODFLOW convention is DELR is along a row, thus has items = NCOL; see page XXX in AW&H (2015)
DELC = LY / NROW # recall that MODFLOW convention is DELC is along a column, thus has items = NROW; see page XXX in AW&H (2015)
DELV = (ZTOP - ZBOT) / NLAY
BOTM = np.linspace(ZTOP, ZBOT, NLAY + 1)
HK = 200.
VKA = 1.
RCH = 0.00305
# WELLQ = 0. #not needed for Problem 5.1
print "DELR =", DELR, " DELC =", DELC, ' DELV =', DELV
print "BOTM =", BOTM
print "Recharge =", RCH
In [25]:
# Create the discretization object
TOP = ZTOP * np.ones((NROW, NCOL),dtype=np.float)
DIS_PACKAGE = flopy.modflow.ModflowDis(MF, NLAY, NROW, NCOL, delr=DELR, delc=DELC,
top=TOP, botm=BOTM[1:], laycbd=0)
In [26]:
# Variables for the BAS package
IBOUND = np.ones((NLAY, NROW, NCOL), dtype=np.int32) # all nodes are active (IBOUND = 1)
# make the top of the profile specified head by setting the IBOUND = -1
IBOUND[:, 0, :] = -1 #don't forget arrays are zero-based! Sets first row
IBOUND[:, :, 0] = -1 # Sets first column
print IBOUND
In [27]:
STRT = 1 * np.ones((NLAY, NROW, NCOL), dtype=np.float32) # set starting head to 1` through out model domain
STRT[:, 0, :] = 0. # top row ocean elevation for setting constant head
STRT[:, :, 0] = 0. # first column ocean elevation for setting constant head
print STRT
In [28]:
BAS_PACKAGE = flopy.modflow.ModflowBas(MF, ibound=IBOUND, strt=STRT)
# print BAS_PACKAGE # uncomment this at far left to see the information about the flopy BAS object
In [29]:
#delete earlier files to prevent us from reading old results
modelfiles = os.listdir(modelpath)
for filename in modelfiles:
f = os.path.join(modelpath, filename)
if modelname in f:
try:
os.remove(f)
print 'Deleted: ', filename
except:
print 'Unable to delete: ', filename
In [30]:
#Now write the model input files
MF.write_input()
# return current working directory
print "You can check the newly created files in", modelpath
In [31]:
#Run MODFLOW
silent = False #Print model output to screen?
pause = False #Require user to hit enter? Doesn't mean much in Ipython notebook
report = True #Store the output from the model in buff
success, buff = MF.run_model(silent=silent, pause=pause, report=report)
In [32]:
#imports for plotting and reading the MODFLOW binary output file
import matplotlib.pyplot as plt
import flopy.utils.binaryfile as bf
#Create the headfile object and grab the results for last time.
headfile = os.path.join(modelpath, modelname + '.hds')
headfileobj = bf.HeadFile(headfile)
#Get a list of times that are contained in the model
times = headfileobj.get_times()
print 'Headfile (' + modelname + '.hds' + ') contains the following list of times: ', times
#Get a numpy array of heads for totim = 1.0
#The get_data method will extract head data from the binary file.
HEAD = headfileobj.get_data(totim=1.0)
#Print statistics on the head
print 'Head statistics'
print ' min: ', HEAD.min()
print ' max: ', HEAD.max()
print ' std: ', HEAD.std()
In [33]:
#Create a contour plot of heads
FIG = plt.figure(figsize=(12,10))
#setup contour levels and plot extent
LEVELS = np.arange(0., 26., 5.)
EXTENT = (DELR/2., LX - DELR/2., DELC/2., LY - DELC/2.)
print 'Contour Levels: ', LEVELS
print 'Extent of domain: ', EXTENT
#Make a contour plot on the first axis
AX1 = FIG.add_subplot(1, 2, 1, aspect='equal')
AX1.set_xlabel("x")
AX1.set_ylabel("y")
YTICKS = np.arange(0, 28000, 4000)
AX1.set_yticks(YTICKS)
AX1.set_title("P5.1 Island Recharge Problem")
AX1.text(900, 10000, r"side ocean boundary condition", fontsize=10, color="blue", rotation='vertical')
AX1.text(4000, 23900, r"top ocean boundary condition", fontsize=10, color="blue")
AX1.contour(np.flipud(HEAD[0, :, :]), levels=LEVELS, extent=EXTENT)
#Make a color flood on the second axis
AX2 = FIG.add_subplot(1, 2, 2, aspect='equal')
AX2.set_xlabel("x")
AX2.set_ylabel("y")
AX2.set_yticks(YTICKS)
AX2.set_title("P5.1 color flood")
AX2.text(4000, 23900, r"top ocean boundary condition", fontsize=10, color="white")
AX2.text(600, 10500, r"side ocean boundary condition", fontsize=10, color="white", rotation='vertical')
cax = AX2.imshow(HEAD[0, :, :], extent=EXTENT, interpolation='nearest')
cbar = FIG.colorbar(cax, orientation='vertical', shrink=0.45)
In [34]:
print HEAD
Write down the new mass balance information for this 1000 ft grid size run (or rename the *.list file). Compare your previous result obtained with the 4000 foot grid size.
In [ ]:
In [35]:
# Same as above but now 500 foot nodal spacing
LX = 12500. # same as before - half width of island + one node for constant head boundary condition
LY = 24500. # same as before - half height of island + one node for constant head boundary condition
NLAY = 1
NROW = 49
NCOL =25
DELR = LX / NCOL # recall that MODFLOW convention is DELR is along a row, thus has items = NCOL; see page XXX in AW&H (2015)
DELC = LY / NROW # recall that MODFLOW convention is DELC is along a column, thus has items = NROW; see page XXX in AW&H (2015)
DELV = (ZTOP - ZBOT) / NLAY
BOTM = np.linspace(ZTOP, ZBOT, NLAY + 1)
HK = 200.
VKA = 1.
RCH = 0.00305
# WELLQ = 0. #not needed for Problem 5.1
print "DELR =", DELR, " DELC =", DELC, ' DELV =', DELV
print "BOTM =", BOTM
print "Recharge =", RCH
In [36]:
# Create the discretization object
TOP = ZTOP * np.ones((NROW, NCOL),dtype=np.float)
DIS_PACKAGE = flopy.modflow.ModflowDis(MF, NLAY, NROW, NCOL, delr=DELR, delc=DELC,
top=TOP, botm=BOTM[1:], laycbd=0)
# Variables for the BAS package
IBOUND = np.ones((NLAY, NROW, NCOL), dtype=np.int32) # all nodes are active (IBOUND = 1)
# make the top of the profile specified head by setting the IBOUND = -1
IBOUND[:, 0, :] = -1 #don't forget arrays are zero-based! Sets first row
IBOUND[:, :, 0] = -1 # Sets first column
print IBOUND
In [37]:
STRT = 1 * np.ones((NLAY, NROW, NCOL), dtype=np.float32) # set starting head to 1` through out model domain
STRT[:, 0, :] = 0. # top row ocean elevation for setting constant head
STRT[:, :, 0] = 0. # first column ocean elevation for setting constant head
print STRT
In [38]:
BAS_PACKAGE = flopy.modflow.ModflowBas(MF, ibound=IBOUND, strt=STRT)
# print BAS_PACKAGE # uncomment this at far left to see the information about the flopy BAS object
#delete earlier files to prevent us from reading old results
modelfiles = os.listdir(modelpath)
for filename in modelfiles:
f = os.path.join(modelpath, filename)
if modelname in f:
try:
os.remove(f)
print 'Deleted: ', filename
except:
print 'Unable to delete: ', filename
#Now write the model input files
MF.write_input()
# return current working directory
print "New files written. You can check them in", modelpath
#Run MODFLOW
silent = False #Print model output to screen?
pause = False #Require user to hit enter? Doesn't mean much in Ipython notebook
report = True #Store the output from the model in buff
success, buff = MF.run_model(silent=silent, pause=pause, report=report)
#imports for plotting and reading the MODFLOW binary output file
import matplotlib.pyplot as plt
import flopy.utils.binaryfile as bf
#Create the headfile object and grab the results for last time.
headfile = os.path.join(modelpath, modelname + '.hds')
headfileobj = bf.HeadFile(headfile)
#Get a list of times that are contained in the model
times = headfileobj.get_times()
print 'Headfile (' + modelname + '.hds' + ') contains the following list of times: ', times
#Get a numpy array of heads for totim = 1.0
#The get_data method will extract head data from the binary file.
HEAD = headfileobj.get_data(totim=1.0)
#Print statistics on the head
print 'Head statistics'
print ' min: ', HEAD.min()
print ' max: ', HEAD.max()
print ' std: ', HEAD.std()
Note that the max head is getting closer to the analytical solution of 20 ft.
In [39]:
#Create a contour plot of heads
FIG = plt.figure(figsize=(12,10))
#setup contour levels and plot extent
LEVELS = np.arange(0., 26., 5.)
EXTENT = (DELR/2., LX - DELR/2., DELC/2., LY - DELC/2.)
print 'Contour Levels: ', LEVELS
print 'Extent of domain: ', EXTENT
#Make a contour plot on the first axis
AX1 = FIG.add_subplot(1, 2, 1, aspect='equal')
AX1.set_xlabel("x")
AX1.set_ylabel("y")
YTICKS = np.arange(0, 28000, 4000)
AX1.set_yticks(YTICKS)
AX1.set_title("P5.1 Island Recharge Problem")
AX1.text(500, 10000, r"side ocean boundary condition", fontsize=10, color="blue", rotation='vertical')
AX1.text(4000, 23900, r"top ocean boundary condition", fontsize=10, color="blue")
AX1.contour(np.flipud(HEAD[0, :, :]), levels=LEVELS, extent=EXTENT)
#Make a color flood on the second axis
AX2 = FIG.add_subplot(1, 2, 2, aspect='equal')
AX2.set_xlabel("x")
AX2.set_ylabel("y")
AX2.set_yticks(YTICKS)
AX2.set_title("P5.1 color flood")
AX2.text(4000, 23900, r"top ocean boundary condition", fontsize=10, color="white")
AX2.text(400, 10500, r"side ocean boundary condition", fontsize=10, color="white", rotation='vertical')
cax = AX2.imshow(HEAD[0, :, :], extent=EXTENT, interpolation='nearest')
cbar = FIG.colorbar(cax, orientation='vertical', shrink=0.45)
Again write down the mass balance information from the *.list file for comparison to other grid sizes.
In [40]:
# Same as above but now 250 foot nodal spacing
LX = 12250. # same as before - half width of island + one node for constant head boundary condition
LY = 24250. # same as before - half height of island + one node for constant head boundary condition
NLAY = 1
NROW = 97
NCOL =49
DELR = LX / NCOL # recall that MODFLOW convention is DELR is along a row, thus has items = NCOL; see page XXX in AW&H (2015)
DELC = LY / NROW # recall that MODFLOW convention is DELC is along a column, thus has items = NROW; see page XXX in AW&H (2015)
DELV = (ZTOP - ZBOT) / NLAY
BOTM = np.linspace(ZTOP, ZBOT, NLAY + 1)
HK = 200.
VKA = 1.
RCH = 0.00305
# WELLQ = 0. #not needed for Problem 5.1
print "DELR =", DELR, " DELC =", DELC, ' DELV =', DELV
print "BOTM =", BOTM
print "Recharge =", RCH
In [41]:
# Create the discretization object
TOP = ZTOP * np.ones((NROW, NCOL),dtype=np.float)
DIS_PACKAGE = flopy.modflow.ModflowDis(MF, NLAY, NROW, NCOL, delr=DELR, delc=DELC,
top=TOP, botm=BOTM[1:], laycbd=0)
# Variables for the BAS package
IBOUND = np.ones((NLAY, NROW, NCOL), dtype=np.int32) # all nodes are active (IBOUND = 1)
# make the top of the profile specified head by setting the IBOUND = -1
IBOUND[:, 0, :] = -1 #don't forget arrays are zero-based! Sets first row
IBOUND[:, :, 0] = -1 # Sets first column
STRT = 1 * np.ones((NLAY, NROW, NCOL), dtype=np.float32) # set starting head to 1` through out model domain
STRT[:, 0, :] = 0. # top row ocean elevation for setting constant head
STRT[:, :, 0] = 0. # first column ocean elevation for setting constant head
BAS_PACKAGE = flopy.modflow.ModflowBas(MF, ibound=IBOUND, strt=STRT)
# print BAS_PACKAGE # uncomment this at far left to see the information about the flopy BAS object
#delete earlier files to prevent us from reading old results
modelfiles = os.listdir(modelpath)
for filename in modelfiles:
f = os.path.join(modelpath, filename)
if modelname in f:
try:
os.remove(f)
print 'Deleted: ', filename
except:
print 'Unable to delete: ', filename
#Now write the model input files
MF.write_input()
# return current working directory
print "New files written. You can check them in", modelpath
#Run MODFLOW
silent = False #Print model output to screen?
pause = False #Require user to hit enter? Doesn't mean much in Ipython notebook
report = True #Store the output from the model in buff
success, buff = MF.run_model(silent=silent, pause=pause, report=report)
#imports for plotting and reading the MODFLOW binary output file
import matplotlib.pyplot as plt
import flopy.utils.binaryfile as bf
#Create the headfile object and grab the results for last time.
headfile = os.path.join(modelpath, modelname + '.hds')
headfileobj = bf.HeadFile(headfile)
#Get a list of times that are contained in the model
times = headfileobj.get_times()
print 'Headfile (' + modelname + '.hds' + ') contains the following list of times: ', times
#Get a numpy array of heads for totim = 1.0
#The get_data method will extract head data from the binary file.
HEAD = headfileobj.get_data(totim=1.0)
#Print statistics on the head
print 'Head statistics'
print ' min: ', HEAD.min()
print ' max: ', HEAD.max()
print ' std: ', HEAD.std()
Note that the max head is again getting closer to the analytical solution of 20 ft.
In [42]:
#Create a contour plot of heads
FIG = plt.figure(figsize=(12,10))
#setup contour levels and plot extent
LEVELS = np.arange(0., 26., 5.)
EXTENT = (DELR/2., LX - DELR/2., DELC/2., LY - DELC/2.)
print 'Contour Levels: ', LEVELS
print 'Extent of domain: ', EXTENT
#Make a contour plot on the first axis
AX1 = FIG.add_subplot(1, 2, 1, aspect='equal')
AX1.set_xlabel("x")
AX1.set_ylabel("y")
YTICKS = np.arange(0, 28000, 4000)
AX1.set_yticks(YTICKS)
AX1.set_title("P5.1 Island Recharge Problem")
AX1.text(220, 10000, r"side ocean boundary condition", fontsize=10, color="blue", rotation='vertical')
AX1.text(4000, 23800, r"top ocean boundary condition", fontsize=10, color="blue")
AX1.contour(np.flipud(HEAD[0, :, :]), levels=LEVELS, extent=EXTENT)
#Make a color flood on the second axis
AX2 = FIG.add_subplot(1, 2, 2, aspect='equal')
AX2.set_xlabel("x")
AX2.set_ylabel("y")
AX2.set_yticks(YTICKS)
AX2.set_title("P5.1 color flood")
AX2.text(4000, 23800, r"top ocean boundary condition", fontsize=10, color="white")
AX2.text(200, 10000, r"side ocean boundary condition", fontsize=10, color="white",rotation='vertical')
cax = AX2.imshow(HEAD[0, :, :], extent=EXTENT, interpolation='nearest')
cbar = FIG.colorbar(cax, orientation='vertical', shrink=0.45)
In [ ]: