In [1]:
%pylab inline --no-import-all
In [2]:
import os
import shutil
from subprocess import Popen
import math
import numpy as np
#import pylab as pl
import scipy.signal as signal
import scipy.fftpack as fftpack
import matplotlib.pyplot as mplt
import flopy.modflow as mf
import flopy.utils as fu
In [3]:
#--default matplotlib
mplt.rcParams['font.sans-serif'] = 'Arial'
mplt.rcParams['font.serif'] = 'Times'
mplt.rcParams['font.cursive'] = 'Zapf Chancery'
mplt.rcParams['font.fantasy'] = 'Comic Sans MS'
mplt.rcParams['font.monospace'] = 'Courier New'
mplt.rcParams['pdf.compression'] = 0
mplt.rcParams['pdf.fonttype'] = 42
mplt.rcParams['mathtext.default'] = 'regular'
ticksize = 8
mplt.rcParams['legend.fontsize'] = 8
mplt.rcParams['axes.labelsize'] = 8
mplt.rcParams['xtick.labelsize'] = ticksize
mplt.rcParams['ytick.labelsize'] = ticksize
In [4]:
home = os.path.expanduser("~")
working_dir = os.path.join(home, 'TidalModflow')
#working_dir = os.path.join('D:/','Users', 'jdhughes', '_Data', 'Projects', 'Tidal', 'TidalModflow')
if not os.path.exists(working_dir):
print 'creating...{0}'.format(working_dir)
os.makedirs(working_dir)
print 'Module working directory: {0}'.format(working_dir)
os.chdir(working_dir)
In [5]:
a = 1.0
lbda = 0.125
days = 7.
ctime = '15 min'
#ctime = '1 min'
#ctime = '1 hr'
if ctime == '1 hr':
tmult = 24
elif ctime == '30 min':
tmult = 24 * 2
elif ctime == '15 min':
tmult = 24 * 4
elif ctime == '10 min':
tmult = 24 * 6
elif ctime == '1 min':
tmult = 24 * 60
delti = float(tmult)
delt = 1. / delti
ntimes = tmult * int(days)
times = np.linspace(0, days, ntimes)
rad = np.copy(times) / lbda
tide = a * np.sin(rad)
In [6]:
fig = mplt.figure()
mplt.plot(times, tide, color='blue', linewidth=0.5)
mplt.xlim((0.,days))
mplt.ylim((-a * 1.25,a * 1.25))
mplt.show()
In [7]:
nlay, nrow, ncol = 1, 1, 250
dy = 1.
dx = np.empty(ncol, np.float)
xcell = np.empty(ncol, np.float)
xmax = 100000.
#dx.fill(xmax / float(ncol))
dxmult = 1.05
dx[0] = xmax * ((dxmult-1) / (math.pow(dxmult, float(ncol)) - 1))
xcell[0] = dx[0] / 2.
for idx in xrange(1, ncol):
dx[idx] = dx[idx-1] * dxmult
xcell[idx] = xcell[idx-1] + 0.5 * (dx[idx-1] + dx[idx])
xanal = np.copy(xcell) - dx[0] / 2.
In [8]:
exe_name='mf2005.exe'
model_name = 'TidalModflow'
nper, itmuni, lenuni = ntimes-1, 4, 0
perlen = np.empty(nper, np.float)
perlen.fill(delt)
nstp = np.ones(nper, np.int)
tsmult = np.ones(nper, np.float)
steady = False
hnoflo = 1.0e+30
hdry = 1.0e+31
top = 10.
botm = -30.
ibound = np.ones((nlay,nrow,ncol), np.int)
ibound[:,:,0] = -1
ibound[:,:,-1] = -1
h0, hl = 0., 2.5
strt = (np.sqrt((h0 * h0) - (xcell / xmax) * (h0 * h0 - hl * hl))).reshape((nlay, nrow, ncol))
laytype = np.ones(nlay, dtype=np.int)
layavg = np.zeros(nlay, dtype=np.int)
#--solver iteration data
nouter = 100
ninner = 50
#--solver convergence data
hclose = 1.0e-6
rclose = 0.1
relax = 0.97
Build CHD from tides
In [9]:
chd_data = []
tidal_data = []
for idx in xrange(len(times)-1):
t = [1,1,1,tide[idx],tide[idx+1]]
chd_data.append([t])
tidal_data.append(tide[idx+1])
tidal_data = np.array(tidal_data)
In [10]:
#--build list of observation locations
ts_list = []
for jcol in xrange(ncol):
ts_list.append((1, 1, jcol+1))
In [11]:
#--storage parameters
ss = 1.0e-5
sy = 0.15
#--calculate average storativity
xt = np.arange(0,xmax+1.,1.)
b = (np.sqrt((h0 * h0) - (xt / xmax) * (h0 * h0 - hl * hl))) - botm
S = np.average(b) * sy
In [12]:
all_hd = np.array([0.01, 0.1, 1., 10., 100., 1000., 10000.])
#all_hd = np.array([1., 10.])
all_kh = np.copy(all_hd) * S
In [13]:
a_sim = np.zeros(ncol, np.float)
off_sim = np.zeros(ncol, np.float)
In [14]:
fig = mplt.figure()
ax1 = mplt.subplot(2,1,1)
ax2 = mplt.subplot(2,1,2)
ax1min = None
ax1max = None
cl = []
m = None
m = mf.Modflow(modelname=model_name, exe_name=exe_name)
for idx, [hk,hd] in enumerate(zip(all_kh,all_hd)):
if idx == 0:
dis = mf.ModflowDis(m, nlay=nlay, nrow=nrow, ncol=ncol, nper=nper,
delr=dx, delc=dy, laycbd=0, top=top, botm=botm,
perlen=perlen, nstp=nstp, tsmult=tsmult, itmuni=itmuni, lenuni=lenuni, steady=steady)
bas = mf.ModflowBas(m, hnoflo=hnoflo, ibound=ibound, strt=strt, stoper=0.01)
lpf = mf.ModflowLpf(m, laytyp=laytype, layavg=layavg, layvka=0, laywet=0,
iwdflg=0, wetfct=0, ihdwet=0, hdry=hdry,
hk=hk, hani=1.0, vka=hk, ss=ss, sy=sy)
chd = mf.ModflowChd(m, layer_row_column_shead_ehead=chd_data)
pcg = mf.ModflowPcg(m, hclose=hclose, rclose=rclose, mxiter=nouter, iter1=ninner, relax=relax, iprpcg=1, mutpcg=1)
oc = mf.ModflowOc(m, ihedfm=-4, iddnfm=-4, save_head_every=1)
m.write_input()
else:
packlist = 'LPF'
m.remove_package(packlist)
lpf = mf.ModflowLpf(m, laytyp=laytype, layavg=layavg, layvka=0, laywet=0,
iwdflg=0, wetfct=0, ihdwet=0, hdry=hdry,
hk=hk, hani=1.0, vka=hk, ss=ss, sy=sy)
m.write_input(SelPackList=[packlist])
#--run the model
Popen([exe_name, model_name + '.nam']).communicate()[0]
#m.run_model2(stdout=False)
#--read the data for each column
hds = fu.HeadFile('{0}.hds'.format(model_name))
v = hds.get_ts(ts_list)
nsamples = v.shape[0]
dt = np.arange(1-nsamples, nsamples)
hds.file.close()
#--calculate offset
tt = v[:, 1] - np.average(v[:, 1])
A = fftpack.fft(tt)
jdx = np.abs(A).argmax()
AO = np.angle(A) / (360. * delt)
a0 = AO[jdx]
for jcol in xrange(ncol):
tp = v[:, jcol+1] - np.average(v[:, jcol+1])
#--fft
A = fftpack.fft(tp)
#--amplitude
ta = np.abs(A).max()
jdx = np.abs(A).argmax()
#--offset
AO = np.angle(A) / (360. * delt)
to = AO[jdx]
#xcorr = signal.correlate(tp, tt)
#to = float(dt[xcorr.argmax()] * delhrs)
#--fill arrays
off_sim[jcol] = (a0 - to) * 24.
a_sim[jcol] = ta
#--
jdx = 0
d0 = 0.0
for jcol in xrange(0, ncol-1):
d1 = off_sim[jcol] - off_sim[jcol+1]
if jcol > 0:
sd0 = np.sign(d0)
sd1 = np.sign(d1)
if abs(d0) > 0.1:
if sd1 != sd0:
jdx = jcol + 1
if ax1min is None:
ax1min = off_sim[:jdx].min()
else:
ax1min = max(ax1min, off_sim[:jdx].min())
if ax1max is None:
ax1max = off_sim[:jdx].max()
else:
ax1max = min(ax1max, off_sim[:jdx].max())
#print jcol, sd0, sd1, d0, d1
break
d0 = d1
#if a_sim[jcol] < 0.1:
# jdx = jcol
# off_base = off_sim[jdx-1]
# break
#off_sim[jdx:] = off_base
off_sim[jdx:] = np.nan
#--plot values
ax1.semilogx(xcell, off_sim, linewidth=0.5, label='{0:5.2g}'.format(hd))
p2, = ax2.semilogx(xcell, a*a_sim/np.abs(a_sim).max(), linewidth=0.5, label='{0:5.2g}'.format(hd))
cl.append(p2.get_color())
for idx, hd in enumerate(all_hd):
tm = np.sqrt(lbda/((delti/2.) * float(hd)))
xx = np.multiply(xanal, -1.0) * tm
aa = np.exp(xx)
tm = np.sqrt(lbda/(float(hd)))
xx = np.copy(xanal) * tm
#--plot values
ax1.semilogx(xanal, xx, linewidth=0.5, linestyle=':', color=cl[idx], label='_None')
ax2.semilogx(xanal, aa, linewidth=0.5, linestyle=':', color=cl[idx], label='_None')
l = ax1.legend(loc='best',ncol=2, title='Hydraulic\ndiffusivity')
l.draw_frame(False)
l.get_title().set_ha("center")
ax1.set_xlim((1.e-2,xmax))
ax1.set_ylim((ax1min,ax1max))
ax1.set_ylabel('Time offset, minutes')
ax2.set_xlim((1.e-2,xmax))
ax2.set_ylabel(r'$\Delta$ amplitude, meters')
ax2.set_xlabel('Distance from coast, meters')
mplt.show()
Get out of the working directory and delete it
In [15]:
os.chdir(home)
shutil.rmtree(working_dir, ignore_errors=True, onerror=None)
In [15]: