Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
# control file
# Add OC package to the MODFLOW model
options = ['PRINT HEAD', 'PRINT DRAWDOWN', 'PRINT BUDGET',
'SAVE HEAD', 'SAVE DRAWDOWN', 'SAVE BUDGET',
'SAVE IBOUND', 'DDREFERENCE']
idx = 0
spd = dict()
for sp in mf.dis.nstp.array:
stress_period = idx
step = sp - 1
ke = (stress_period, step)
idx = idx + 1
spd[ke] = [options[3], options[2], options[5]]
mf.remove_package('OC')
oc = flopy.modflow.ModflowOc(mf, stress_period_data=spd, cboufm='(20i5)')
oc.write_file()
# write hob files
hob_df = pd.read_csv(r"misc_hob.csv", header=None)
start_time = hob_df.values[0][0]
end_time = hob_df.values[1][0]
#hobs
mf_name2 = r".\mf_EnS\pp_model.nam"
mf2 = flopy.modflow.Modflow.load(mf_name2, load_only=['DIS', 'BAS6','HOB'])
obs_data = mf2.hob.obs_data
if len(steady) > 1:
time_index = np.arange(start_time+1,end_time+1)
else:
nam_file = "freyberg.nam"
mo = flopy.modflow.Modflow.load(nam_file, model_ws=org_model_ws, check=False,forgive=False)
perlen = np.ones((365))
m = flopy.modflow.Modflow("freyberg_transient",model_ws=os.path.join("da","freyberg","truth"),version="mfnwt",
external_path=".")
flopy.modflow.ModflowDis(m,nrow=mo.nrow,ncol=mo.ncol,nlay=1,delr=mo.dis.delr,
delc=mo.dis.delc,top=mo.dis.top,botm=mo.dis.botm[-1],nper=len(perlen),perlen=perlen)
flopy.modflow.ModflowBas(m,ibound=mo.bas6.ibound[0],strt=mo.bas6.strt[0])
flopy.modflow.ModflowUpw(m,laytyp=mo.upw.laytyp,hk=mo.upw.hk[0],vka=mo.upw.vka[0],ss=0.00001,sy=0.01)
flopy.modflow.ModflowNwt(m)
oc_data = {}
for iper in range(m.nper):
oc_data[iper,0] = ["save head","save budget"]
flopy.modflow.ModflowOc(m,stress_period_data=oc_data)
flopy.modflow.ModflowRch(m,rech=mo.rch.rech.array[0])
wel_data = mo.wel.stress_period_data[0]
wel_data["k"][:] = 0
flopy.modflow.ModflowWel(m,stress_period_data={0:wel_data})
flopy.modflow.ModflowSfr2(m,nstrm=mo.sfr.nstrm,nss=mo.sfr.nss,istcb2=90,segment_data=mo.sfr.segment_data,reach_data=mo.sfr.reach_data)
m.write_input()
pyemu.os_utils.run("mfnwt {0}.nam".format(m.name),cwd=m.model_ws)
return m
file = '{}_ib{:02d}.ref'.format(name, kdx + 1)
files.append(file)
hfiles.append(hfile)
np.savetxt(file, ibound[kdx, :, :], fmt='%5d')
bas = flopy.modflow.ModflowBas(ml, ibound=files, strt=hfiles)
# The aquifer properties (really only the hydraulic conductivity) are defined with the
# LPF package.
lpf = flopy.modflow.ModflowLpf(ml, hk=Kh)
# Finally, we need to specify the solver we want to use (PCG with default values), and the
# output control (using the default values). Then we are ready to write all MODFLOW input
# files and run MODFLOW.
pcg = flopy.modflow.ModflowPcg(ml)
oc = flopy.modflow.ModflowOc(ml)
ml.write_input()
ml.run_model()
# change back to the starting directory
os.chdir(cwdpth)
# Once the model has terminated normally, we can read the heads file. First, a link to the heads
# file is created with `HeadFile`. The link can then be accessed with the `get_data` function, by
# specifying, in this case, the step number and period number for which we want to retrieve data.
# A three-dimensional array is returned of size `nlay, nrow, ncol`. Matplotlib contouring functions
# are used to make contours of the layers or a cross-section.
hds = flopy.utils.HeadFile(os.path.join(workspace, name + '.hds'))
h = hds.get_data(kstpkper=(0, 0))
x = y = np.linspace(0, L, N)
c = plt.contour(x, y, h[0], np.arange(90, 100.1, 0.2))
plt.clabel(c, fmt='%2.1f')
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 = {}
for kper in range(nper):
for kstp in range(nstp[kper]):
stress_period_data[(kper, kstp)] = ['save head',
'save drawdown',
'save budget',
'print head',
'print budget']
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
print('creating...', modelname)
ml = flopy.modflow.Modflow(modelname, version='mf2005', exe_name=mf_name,
model_ws=dirs[0])
discret = flopy.modflow.ModflowDis(ml, nlay=1, ncol=ncol, nrow=nrow,
delr=delr,
delc=1,
top=0, botm=[-40.0],
nper=nper, perlen=perlen, nstp=nstp)
bas = flopy.modflow.ModflowBas(ml, ibound=ibound, strt=0.05)
bcf = flopy.modflow.ModflowBcf(ml, laycon=0, tran=2 * 40)
swi = flopy.modflow.ModflowSwi2(ml, iswizt=55, nsrf=nsurf, istrat=1,
toeslope=0.2, tipslope=0.2,
nu=[0, 0.0125, 0.025],
zeta=z, ssz=ssz, isource=isource,
nsolver=1)
oc = flopy.modflow.ModflowOc(ml,
stress_period_data={(0, 999): ['save head']})
pcg = flopy.modflow.ModflowPcg(ml)
ml.write_input()
# run stratified model
if not skipRuns:
m = ml.run_model(silent=False)
# read stratified results
zetafile = os.path.join(dirs[0], '{}.zta'.format(modelname))
zobj = flopy.utils.CellBudgetFile(zetafile)
zkstpkper = zobj.get_kstpkper()
zeta = zobj.get_data(kstpkper=zkstpkper[-1], text='ZETASRF 1')[0]
zeta2 = zobj.get_data(kstpkper=zkstpkper[-1], text='ZETASRF 2')[0]
#
# vd model
modelname = 'swiex2_vd'
print('creating...', modelname)
bot = np.linspace(-H/Nlay,-H,Nlay)
delrow = delcol = L/(N-1)
dis = fmf.ModflowDis(ml,nlay=Nlay,nrow=N,ncol=N,delr=delrow,delc=delcol,top=0.0,botm=bot,laycbd=0)
Nhalf = (N-1)/2
ibound = np.ones((Nlay,N,N))
ibound[:,0,:] = -1; ibound[:,-1,:] = -1; ibound[:,:,0] = -1; ibound[:,:,-1] = -1
ibound[0,Nhalf,Nhalf] = -1
start = h1 * np.ones((N,N))
start[Nhalf,Nhalf] = h2
bas = fmf.ModflowBas(ml,ibound=ibound,strt=start)
print "o"
lpf = fmf.ModflowLpf(ml, hk=k)
pcg = fmf.ModflowPcg(ml)
oc = fmf.ModflowOc(ml)
ml.write_input()
ml.run_model()
head_file = '/home/kiruba/PycharmProjects/area_of_curve/hydrology/hydrology/mf_files/lake_example.hds'
hds = fut.HeadFile(head_file)
h = hds.get_data(kstpkper=(1,1))
print len(h)
x = y = np.linspace(0, L, N)
c = plt.contour(x, y, h[0], np.arange(90,100.1,0.2))
plt.clabel(c, fmt='%2.1f')
plt.axis('scaled')
plt.show()
c = plt.contour(x,y,h[-1],np.arange(90,100.1,0.2))
plt.clabel(c,fmt='%1.1f')
plt.axis('scaled')
plt.show()
(1,30): [],
(1,39): ['save head', 'save budget', 'print budget'],
(1,40): [],
(1,49): ['save head', 'save budget', 'print budget'],
(1,50): [],
(1,59): ['save head', 'save budget', 'print budget'],
(1,60): [],
(1,69): ['save head', 'save budget', 'print budget'],
(1,70): [],
(1,79): ['save head', 'save budget', 'print budget'],
(1,80): [],
(1,89): ['save head', 'save budget', 'print budget'],
(1,90): [],
(1,99): ['save head', 'save budget', 'print budget'],
(1,100): []}
oc = flopy.modflow.ModflowOc(ml, stress_period_data=stress_period_data)
# write the input files
ml.write_input()
else:
if not os.path.exists(cf_pths[idx]):
os.makedirs(cf_pths[idx])
filelist = [f for f in os.listdir(cf_pths[0])]
sys.stdout.write('copying files from {} to {}\n'.format(cf_pths[0], cf_pths[idx]))
for f in filelist:
if os.path.splitext(f)[1].lower() in exclude:
continue
src = os.path.join(cf_pths[0], f)
dst = os.path.join(cf_pths[idx], f)
shutil.copyfile(src, dst)
return ml, cf_pths