Closed rabbl closed 6 years ago
You are missing a number of MODFLOW packages.
The following runs:
pth = os.path.join(cpth, 'simple')
modelname = 'hob_simple'
nlay, nrow, ncol = 1, 11, 11
shape3d = (nlay, nrow, ncol)
shape2d = (nrow, ncol)
ib = np.ones(shape3d, dtype=np.int)
ib[0, 0, 0] = -1
m = flopy.modflow.Modflow(modelname=modelname, model_ws=pth,
verbose=True, exe_name=exe_name,)
dis = flopy.modflow.ModflowDis(m, nlay=1, nrow=11, ncol=11, nper=2,
perlen=[1, 1])
bas = flopy.modflow.ModflowBas(m, ibound=ib, strt=10.)
lpf = flopy.modflow.ModflowLpf(m)
pcg = flopy.modflow.ModflowPcg(m)
obs = flopy.modflow.HeadObservation(m, layer=0, row=5, column=5,
time_series_data=[[1., 54.4],
[2., 55.2]])
hob = flopy.modflow.ModflowHob(m, iuhobsv=51, hobdry=-9999.,
obs_data=[obs])
# Write the model input files
m.write_input()
# run the modflow-2005 model
if run:
success, buff = m.run_model(silent=False)
assert success, 'could not run simple MODFLOW-2005 model'
FloPy is using the following executable to run the model: /Users/jdhughes/Documents/Development/bin/mac/mf2005
MODFLOW-2005
U.S. GEOLOGICAL SURVEY MODULAR FINITE-DIFFERENCE GROUND-WATER FLOW MODEL
Version 1.12.00 2/3/2017
Using NAME file: hob_simple.nam
Run start date and time (yyyy/mm/dd hh:mm:ss): 2017/08/01 9:48:47
Solving: Stress period: 1 Time step: 1 Ground-Water Flow Eqn.
Solving: Stress period: 2 Time step: 1 Ground-Water Flow Eqn.
Run end date and time (yyyy/mm/dd hh:mm:ss): 2017/08/01 9:48:47
Elapsed run time: 0.004 Seconds
Normal termination of simulation
Does this resolve your issue?
Thanks for the response resolving this issue. Your example is running well with my configuration.
I have tries hob in the tutorial-files and when I try to add observationPoints to the Tuturial_2, inserting it into line 45 it brings the following error:
At line 907 of file obs2bas7.f (unit = 51, file = 'tutorial2.hds')
Fortran runtime error: Format present for UNFORMATTED data transfer
Here is the complete script:
import numpy as np
import flopy
# Model domain and grid definition
Lx = 1000.
Ly = 1000.
ztop = 10.
zbot = -50.
nlay = 1
nrow = 10
ncol = 10
delr = Lx / ncol
delc = Ly / nrow
delv = (ztop - zbot) / nlay
botm = np.linspace(ztop, zbot, nlay + 1)
hk = 1.
vka = 1.
sy = 0.1
ss = 1.e-4
laytyp = 1
# Variables for the BAS package
# Note that changes from the previous tutorial!
ibound = np.ones((nlay, nrow, ncol), dtype=np.int32)
strt = 10. * np.ones((nlay, nrow, ncol), dtype=np.float32)
# Time step parameters
nper = 3
perlen = [1, 100, 100]
nstp = [1, 100, 100]
steady = [True, False, False]
# Flopy objects
modelname = 'tutorial2'
mf = flopy.modflow.Modflow(modelname, exe_name='mf2005')
dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, delr=delr, delc=delc,
top=ztop, botm=botm[1:],
nper=nper, perlen=perlen, nstp=nstp,
steady=steady)
bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt)
lpf = flopy.modflow.ModflowLpf(mf, hk=hk, vka=vka, sy=sy, ss=ss, laytyp=laytyp, ipakcb=53)
pcg = flopy.modflow.ModflowPcg(mf)
obs = flopy.modflow.HeadObservation(mf, layer=0, row=5, column=5, time_series_data=[[1., 54.4], [2., 55.2]])
hob = flopy.modflow.ModflowHob(mf, iuhobsv=51, hobdry=-9999., obs_data=[obs])
# Make list for stress period 1
stageleft = 10.
stageright = 10.
bound_sp1 = []
for il in range(nlay):
condleft = hk * (stageleft - zbot) * delc
condright = hk * (stageright - zbot) * delc
for ir in range(nrow):
bound_sp1.append([il, ir, 0, stageleft, condleft])
bound_sp1.append([il, ir, ncol - 1, stageright, condright])
print('Adding ', len(bound_sp1), 'GHBs for stress period 1.')
# Make list for stress period 2
stageleft = 10.
stageright = 0.
condleft = hk * (stageleft - zbot) * delc
condright = hk * (stageright - zbot) * delc
bound_sp2 = []
for il in range(nlay):
for ir in range(nrow):
bound_sp2.append([il, ir, 0, stageleft, condleft])
bound_sp2.append([il, ir, ncol - 1, stageright, condright])
print('Adding ', len(bound_sp2), 'GHBs for stress period 2.')
# We do not need to add a dictionary entry for stress period 3.
# Flopy will automatically take the list from stress period 2 and apply it
# to the end of the simulation
stress_period_data = {0: bound_sp1, 1: bound_sp2}
# Create the flopy ghb object
ghb = flopy.modflow.ModflowGhb(mf, stress_period_data=stress_period_data)
# Create the well package
# Remember to use zero-based layer, row, column indices!
pumping_rate = -500.
wel_sp1 = [[0, nrow/2 - 1, ncol/2 - 1, 0.]]
wel_sp2 = [[0, nrow/2 - 1, ncol/2 - 1, 0.]]
wel_sp3 = [[0, nrow/2 - 1, ncol/2 - 1, pumping_rate]]
stress_period_data = {0: wel_sp1, 1: wel_sp2, 2: wel_sp3}
wel = flopy.modflow.ModflowWel(mf, stress_period_data=stress_period_data)
# Output control
stress_period_data = {(0, 0): ['save head',
'save drawdown',
'save budget',
'print head',
'print budget']}
save_head_every = 1
oc = flopy.modflow.ModflowOc(mf, stress_period_data=stress_period_data,
compact=True)
# Write the model input files
mf.write_input()
# Run the model
success, mfoutput = mf.run_model(silent=False, pause=False)
if not success:
raise Exception('MODFLOW did not terminate normally.')
# Imports
import matplotlib.pyplot as plt
import flopy.utils.binaryfile as bf
# Create the headfile and budget file objects
headobj = bf.HeadFile(modelname+'.hds')
times = headobj.get_times()
cbb = bf.CellBudgetFile(modelname+'.cbc')
# Setup contour parameters
levels = np.linspace(0, 10, 11)
extent = (delr/2., Lx - delr/2., delc/2., Ly - delc/2.)
print('Levels: ', levels)
print('Extent: ', extent)
# Well point
wpt = ((float(ncol/2)-0.5)*delr, (float(nrow/2-1)+0.5)*delc)
wpt = (450., 550.)
# Make the plots
mytimes = [1.0, 101.0, 201.0]
for iplot, time in enumerate(mytimes):
print('*****Processing time: ', time)
head = headobj.get_data(totim=time)
#Print statistics
print('Head statistics')
print(' min: ', head.min())
print(' max: ', head.max())
print(' std: ', head.std())
# Extract flow right face and flow front face
frf = cbb.get_data(text='FLOW RIGHT FACE', totim=time)[0]
fff = cbb.get_data(text='FLOW FRONT FACE', totim=time)[0]
#Create the plot
#plt.subplot(1, len(mytimes), iplot + 1, aspect='equal')
plt.subplot(1, 1, 1, aspect='equal')
plt.title('stress period ' + str(iplot + 1))
modelmap = flopy.plot.ModelMap(model=mf, layer=0)
qm = modelmap.plot_ibound()
lc = modelmap.plot_grid()
qm = modelmap.plot_bc('GHB', alpha=0.5)
cs = modelmap.contour_array(head, levels=levels)
plt.clabel(cs, inline=1, fontsize=10, fmt='%1.1f', zorder=11)
quiver = modelmap.plot_discharge(frf, fff, head=head)
mfc = 'None'
if (iplot+1) == len(mytimes):
mfc='black'
plt.plot(wpt[0], wpt[1], lw=0, marker='o', markersize=8,
markeredgewidth=0.5,
markeredgecolor='black', markerfacecolor=mfc, zorder=9)
plt.text(wpt[0]+25, wpt[1]-25, 'well', size=12, zorder=12)
plt.show()
plt.show()
# Plot the head versus time
idx = (0, nrow/2 - 1, ncol/2 - 1)
ts = headobj.get_ts(idx)
plt.subplot(1, 1, 1)
ttl = 'Head at cell ({0},{1},{2})'.format(idx[0] + 1, idx[1] + 1, idx[2] + 1)
plt.title(ttl)
plt.xlabel('time')
plt.ylabel('head')
plt.plot(ts[:, 0], ts[:, 1])
plt.show()
Try changing the unit number for HOB output to another unit number:
hob = flopy.modflow.ModflowHob(mf, iuhobsv=1051, hobdry=-9999., obs_data=[obs])
The default unit for binary head output is 51. The unit numbers are still a little messy in FloPy. Sorry.
Yes, that works. Thanks!
Im not familiar with Fortran.
Why you can't use a file-unit twice? Works it like a link or internal reference to another file?
So we can close this issue I think. Thanks again.
Multiple packages can write to the same unit number provided the files are opened for writing the same type of data (formatted or unformatted). By default the head file is opened for binary output ('wb' in python) and the HOB file is opened for formatted ascii output ('w' in python).
I will close this issue.
I have tried the hob-package-implementation and got some Fortran-runtime errors.
Fortran runtime error: Missing format for FORMATTED data transfer
Here a small sample (as in the documentation):
Returns: