In [1]:

    

import simpegSP as SP
from SimPEG import Mesh, np
from SimPEG.EM import Static
%pylab inline


    

    




Efficiency Warning: Interpolation will be slow, use setup.py!

            python setup.py build_ext --inplace
    
Populating the interactive namespace from numpy and matplotlib

In [2]:

    

# workdir = "/Documents and Settings/"
workdir = "C:/00PLAN_UBC/seepageModeling/Cheong_dam/"
fname = "drawdown(0sec).csv"
fluiddata = SP.Utils.readSeepageModel(workdir+fname)


    

In [3]:

    

xyz =  fluiddata["xyz"]
h =  fluiddata["h"]
Sw =  fluiddata["Sw"]
Kx =  fluiddata["Kx"]
Ky =  fluiddata["Ky"]
P = fluiddata["P"]
Ux =  fluiddata["Ux"]
Uy =  fluiddata["Uy"]
Uy[Uy>1e10] = np.nan
mesh = fluiddata["mesh"]
xsurf, ysurf, yup = fluiddata["xsurf"], fluiddata["ysurf"], fluiddata["yup"]
hcc = fluiddata["hcc"]
actind = fluiddata["actind"]
waterind = fluiddata["waterind"]


    

In [4]:

    

mesh.plotGrid()
plt.gca().set_aspect('equal', adjustable='box')


    

In [5]:

    

fig = plt.figure(figsize = (12, 5))
ax = plt.subplot(111)
dat = Static.Utils.plot2Ddata(xyz, h, ax=ax, ncontour=30, contourOpts={"cmap":"viridis"})
ax.fill_between(xsurf, ysurf, yup, facecolor='white', edgecolor="white")

ax.set_xlim(50, 130)
ax.set_ylim(25, 44)
cb = plt.colorbar(dat[1], orientation="horizontal")
cb.set_label("Total head (m)")


    

In [6]:

    

rho = np.ones(mesh.nC)*1e8
rho[actind] = 200.
rho[waterind] = 0.33


    

In [7]:

    

tempx = np.r_[0. , 82, 125, 0., 0. ]
tempy = np.r_[40., 40, 27.5 , 27.5, 40.]
inds = SP.Utils.PolygonInd(mesh, np.c_[tempx, tempy])
satinds = (inds) & (actind)
# L = np.zeros(mesh.nC)
L = np.ones(mesh.nC)*0.
L[satinds] = 2.5*1e-4


    

In [8]:

    

# coreind = SP.Utils.PolygonInd(mesh, np.c_[tempx, tempy])
# rho[coreind] = res_core


    

In [9]:

    

plt.plot(tempx, tempy, 'k-')


    

    
Out[9]:





[<matplotlib.lines.Line2D at 0xb7297b8>]






    

In [10]:

    

fig = plt.figure(figsize = (12, 5))
ax = plt.subplot(111)
dat = Static.Utils.plot2Ddata(xyz, Kx*Ux, ax=ax, ncontour=30, contourOpts={"cmap":"viridis"})
ax.fill_between(xsurf, ysurf, yup, facecolor='white', edgecolor="white")

ax.set_xlim(50, 130)
ax.set_ylim(25, 44)
cb = plt.colorbar(dat[1], orientation="horizontal")
cb.set_label("Total head (m)")


    

In [11]:

    

fig = plt.figure(figsize = (12, 5))
ax = plt.subplot(111)
dat = Static.Utils.plot2Ddata(xyz, Ky*Uy, ax=ax, ncontour=30, contourOpts={"cmap":"viridis"})
ax.fill_between(xsurf, ysurf, yup, facecolor='white', edgecolor="white")

ax.set_xlim(50, 130)
ax.set_ylim(25, 44)
cb = plt.colorbar(dat[1], orientation="horizontal")
cb.set_label("Total head (m)")


    

In [12]:

    

fig = plt.figure(figsize = (12, 5))
ax = plt.subplot(111)
dat = Static.Utils.plot2Ddata(xyz, np.log10(Ux*h), ax=ax, ncontour=30, contourOpts={"cmap":"viridis"})
ax.fill_between(xsurf, ysurf, yup, facecolor='white', edgecolor="white")

ax.set_xlim(50, 130)
ax.set_ylim(25, 44)
cb = plt.colorbar(dat[1], orientation="horizontal")
cb.set_label("Total head (m)")


    

In [13]:

    

rho = np.ones(mesh.nC)*1e8
rho[actind] = 200.
rho[waterind] = 0.33


    

In [14]:

    

mesh.plotImage(np.log10(1./rho), pcolorOpts={"cmap":"viridis"})


    

    
Out[14]:





(<matplotlib.collections.QuadMesh at 0xc874c50>,)






    

In [15]:

    

tempx = np.r_[0. , 82, 125, 0., 0. ]
tempy = np.r_[40., 40, 27.5 , 27.5, 40.]
inds = SP.Utils.PolygonInd(mesh, np.c_[tempx, tempy])
satinds = (inds) & (actind)
# L = np.zeros(mesh.nC)
L = np.ones(mesh.nC)*0.
L[satinds] = 2.5*1e-4


    

In [16]:

    

mesh.plotImage(satinds, pcolorOpts={"cmap":"viridis"})


    

    
Out[16]:





(<matplotlib.collections.QuadMesh at 0xcea3c88>,)






    

In [17]:

    

from SimPEG import Utils


    

In [18]:

    

airind = ~np.logical_or(actind, waterind)
topo = np.c_[xsurf, ysurf]
pts = mesh.vectorCCx[np.logical_and(mesh.vectorCCx>85, mesh.vectorCCx<120)]
locs = SP.Utils.drapeTopotoLoc(mesh, topo, pts)
locM = locs.copy()
locN = locs.copy()
locN = np.c_[np.ones(locs.shape[0])*locM[0,0], np.ones(locs.shape[0])*locM[0,1]]


    

In [19]:

    

fig = plt.figure(figsize = (12*1.5, 5*1.5))
ax = plt.subplot(111)
sigma = 1./rho
sigma[airind] = np.nan
dat = mesh.plotImage(np.log10(sigma), grid=True, pcolorOpts={"cmap":"viridis"}, ax=ax, clim=(-3, -2))
ax.plot(locM[:,0], locM[:,1], 'go')
ax.plot(locN[:,0], locN[:,1], 'ro', ms=5)
ax.plot(topo[:,0], topo[:,1], 'w-')
ax.plot(tempx, tempy, 'w-', lw=2)
ax.set_xlim(50, 130)
ax.set_ylim(25, 44)
# ax.plot(128, 27, 'wo', ms=5)
plt.gca().set_aspect('equal', adjustable='box')
cb = plt.colorbar(dat[0], orientation="horizontal")


    

In [20]:

    

fig = plt.figure(figsize = (12*1.5, 5*1.5))
ax = plt.subplot(111)
sigma = 1./rho
sigma[airind] = np.nan
mesh.plotImage(L, grid=True, pcolorOpts={"cmap":"viridis"}, ax=ax)
ax.plot(locM[:,0], locM[:,1], 'go')
ax.plot(locN[:,0], locN[:,1], 'r.', ms=5)
ax.plot(topo[:,0], topo[:,1], 'w-')
ax.plot(tempx, tempy, 'w-', lw=2)
ax.set_xlim(50, 130)
ax.set_ylim(25, 44)
plt.gca().set_aspect('equal', adjustable='box')


    

In [21]:

    

prob = SP.Problem_CC(mesh, rho=rho)
rx = SP.Rx.Dipole(locM, locN)
src = SP.Src.StreamingCurrents([rx], L=L, mesh=mesh)
survey = SP.Survey([src])
survey.pair(prob)


    

In [22]:

    

f = prob.fields(hcc)
data = survey.dpred(hcc, f=f)


    

In [23]:

    

fig = plt.figure(figsize = (12, 5))
ax = plt.subplot(111)
temp = hcc.copy()
temp[airind] = np.nan
dat = mesh.plotImage(temp, grid=False, pcolorOpts={"cmap":"viridis"}, ax=ax, clim=(h.min(), h.max()) )
ax.plot(locM[:,0], locM[:,1], 'go')
ax.plot(locN[:,0], locN[:,1], 'r.', ms=5)
ax.plot(topo[:,0], topo[:,1], 'r-')
ax.set_xlim(50, 130)
ax.set_ylim(25, 44)
plt.gca().set_aspect('equal', adjustable='box')
cb = plt.colorbar(dat[0], orientation="horizontal")
cb.set_label("Total head (m)")


    

In [24]:

    

fig = plt.figure(figsize = (12, 5))
ax = plt.subplot(111)
dat = Static.Utils.plot2Ddata(mesh.gridCC, f[src,'phi'], ax=ax, ncontour=30, contourOpts={"cmap":"viridis"})
ax.fill_between(xsurf, ysurf, yup, facecolor='white', edgecolor="white")
ax.set_xlim(50, 130)
ax.set_ylim(25, 44)
cb = plt.colorbar(dat[1], orientation="horizontal")
cb.set_label("Potential (V)")


    

In [25]:

    

fig = plt.figure(figsize = (12, 5))
ax = plt.subplot(111)
sigma = 1./rho
sigma[airind] = np.nan
dat = mesh.plotImage(f[src,'e'], view="vec", vType="F",grid=False, pcolorOpts={"cmap":"viridis"}, ax=ax, streamOpts={"color":"w"})
ax.plot(locM[:,0], locM[:,1], 'go')
ax.plot(locN[:,0], locN[:,1], 'r.', ms=5)
ax.plot(topo[:,0], topo[:,1], 'r-')
ax.plot(tempx, tempy, 'w-', lw=2)
ax.set_xlim(50, 130)
ax.set_ylim(25, 44)
plt.gca().set_aspect('equal', adjustable='box')
cb = plt.colorbar(dat[0], orientation="horizontal")
cb.set_label("Potential (V)")


    

In [26]:

    

figsize(12, 3)
plt.plot(locs[:,0], data*1000, '.-')
plt.xlim(50, 130)
grid(True)


    

In [27]:

    

q = src.eval(prob)


    

In [28]:

    

fig = plt.figure(figsize = (12, 5))
ax = plt.subplot(111)

dat = mesh.plotImage(q, grid=False, pcolorOpts={"cmap":"viridis"}, ax=ax)
# ax.plot(locM[:,0], locM[:,1], 'go')
# ax.plot(locN[:,0], locN[:,1], 'r.', ms=5)
# ax.plot(topo[:,0], topo[:,1], 'r-')
# ax.plot(tempx, tempy, 'w-', lw=2)

ax.set_xlim(50, 130)
ax.set_ylim(25, 44)
plt.gca().set_aspect('equal', adjustable='box')
cb = plt.colorbar(dat[0], orientation="horizontal")
cb.set_label("Potential (V)")


    

In [29]:

    

fig = plt.figure(figsize = (12, 5))
ax = plt.subplot(111)

dat = mesh.plotImage(f[src,"charge"], grid=False, pcolorOpts={"cmap":"viridis"}, ax=ax)
# ax.plot(locM[:,0], locM[:,1], 'go')
# ax.plot(locN[:,0], locN[:,1], 'r.', ms=5)
# ax.plot(topo[:,0], topo[:,1], 'r-')
# ax.plot(tempx, tempy, 'w-', lw=2)

ax.set_xlim(50, 130)
ax.set_ylim(25, 44)
plt.gca().set_aspect('equal', adjustable='box')
cb = plt.colorbar(dat[0], orientation="horizontal")
cb.set_label("Potential (V)")


    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]: