In [143]:
from SimPEG import *
import simpegEM as EM
from simpegem1d import Utils1D
%pylab inline
In [144]:
import matplotlib
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('png')
matplotlib.rcParams['savefig.dpi'] = 100
In [145]:
matplotlib.rcParams.update({'font.size': 14})
In [146]:
meshType = 'CYL'
cs, ncx, ncz, npad = 20., 25, 30, 20
hx = [(cs,ncx), (cs,npad,1.3)]
hz = [(cs,npad,-1.4), (cs,ncz), (cs,npad,1.4)]
mesh = Mesh.CylMesh([hx,1,hz], '00C')
In [148]:
print mesh2D
In [112]:
mesh2D = Mesh.TensorMesh([hx,1,hz], '00C')
In [113]:
active = mesh.vectorCCz<0.
layer1 = (mesh.vectorCCz<0.) & (mesh.vectorCCz>=-60.)
layer2 = (mesh.vectorCCz<-60) & (mesh.vectorCCz>=-160.)
actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap
In [114]:
sig_half = 1e-3
sig_air = 1e-8
sig_layer1 = 1./300
sig_layer2 = 1./5.
sigma = np.ones(mesh.nCz)*sig_air
sigma_half = sigma.copy()
sigma[active] = sig_half
sigma[layer1] = sig_layer1
sigma[layer2] = sig_layer2
sigma_half[active] = sig_half
sigma_half[layer1] = sig_layer1
sigma_half[layer2] = 1./500.
mtrue = np.log(sigma[active])
In [149]:
sig_sea = sigma[active]
sig_half = sigma_half[active]
In [150]:
sig_sea.shape
Out[150]:
In [116]:
# t = -(np.arange(8)+1)
# dummy = []
# for i in range(8):
# dummy.append()
In [117]:
print vmin, vmax
In [118]:
print sigact.shape, zact.shape
print sig.shape, z.shape
In [119]:
mesh.vectorCCz[active].shape
sigma[active].shape
Out[119]:
In [120]:
print test
In [140]:
test = mapping*np.log(sig_sea)
dat = mesh2D.plotSlice(np.log10(test), normal='Y', grid=True)
ylim(-400, 200)
xlim(0, 200)
xlabel("Easting (m)", fontsize = 14)
ylabel("Depth (m)", fontsize = 14)
title("$\mathcal{M}_{exp}[\mathcal{M}_{1D}[\mathcal{M}_{surj}[m]]]$", fontsize = 18)
vmin, vmax = np.log10(test).min(), np.log10(test).max()
cb = plt.colorbar(dat[0], ticks = np.linspace(vmin, vmax, 5), format="10$^{%3.1f}$")
cb.set_label("Conductivity (S/m)")
figsize(5, 5)
In [75]:
sigact = np.log(np.repeat(sigma[active], 2, axis=0))
zact = np.repeat(mesh.vectorCCz[active][1:], 2, axis=0)
zact = np.r_[mesh.vectorCCz[active][0], zact, mesh.vectorCCz[active][-1]]
In [141]:
sig = np.log(np.repeat(sigma, 2, axis=0))
z = np.repeat(mesh.vectorCCz[1:], 2, axis=0)
z = np.r_[mesh.vectorCCz[0], z, mesh.vectorCCz[-1]]
vmin, vmax = sig.min(), sig.max()
fig = plt.figure(figsize = (9, 4))
ax1 = plt.subplot(121)
ax2 = plt.subplot(122)
for locz in mesh.vectorCCz[active]:
ax1.plot(np.linspace(vmin, vmax, 100), np.ones(100)*locz, 'b--', lw = 0.5)
for locz in mesh.vectorCCz:
ax2.plot(np.linspace(vmin, vmax, 100), np.ones(100)*locz, 'b--', lw = 0.5)
ax1.plot(sigact, zact, 'k')
ax1.set_xlabel("$log(\sigma)$", fontsize = 18)
ax1.set_ylabel("Depth (m)")
ax1.set_xlim(-19, -1)
ax1.set_ylim(-400, 200)
ax1.set_title("$m$", fontsize = 18)
ax2.plot(sig, z, 'k')
ax2.set_xlim(-19, -1)
ax2.set_ylim(-400, 200)
ax2.set_xlabel("$log(\sigma)$", fontsize = 18)
ax2.set_title("$\mathcal{M}_{surj}[m]$", fontsize = 18)
Out[141]:
In [24]:
fig, ax = plt.subplots(1,1, figsize = (4, 5))
Utils1D.plotLayer(sigma, mesh.vectorCCz, showlayers=True, ax = ax)
ax.set_xlabel("Log conductivity")
ax.set_ylim(-400., 200.)
ax.text(2e-6, -110, 'Brine saturated layer',fontsize = 12)
ax.text(2e-6, -200, 'Basement',fontsize = 12)
Out[24]:
In [8]:
fig, ax = plt.subplots(1,1, figsize = (3, 5))
ax.plot(1000, 1000, 'k')
ax.plot(1000, 1000, 'b')
ax.legend(('With intrusion', 'Without intrusion'), loc=4,fontsize = 10)
Utils1D.plotLayer(sig_half, mesh.vectorCCz[active], ax = ax, **{'color':'b'})
Utils1D.plotLayer(sig_sea, mesh.vectorCCz[active], showlayers=True, ax = ax)
ax.set_ylim(-400., 0.)
ax.text(2e-3, -110, 'Brine saturated layer',fontsize = 12)
ax.text(2e-3, -200, 'Basement',fontsize = 12)
Out[8]:
In [9]:
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping, verbose=False)
In [10]:
prb.Solver = SolverLU
prb.timeSteps = [(1e-4/10, 15), (1e-4/10, 15), (1e-2/10, 15), (1e1/10, 15)]
In [11]:
time = np.logspace(-4, -1, 51)
In [12]:
x = np.linspace(0, 250, 26)
xyz_rx = Utils.ndgrid(x, np.r_[0.], np.r_[0.])
In [15]:
rx = EM.TDEM.RxTDEM(xyz_rx, time, 'bz')
tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'CircularLoop_MVP', [rx])
tx.radius = 250.
survey = EM.TDEM.SurveyTDEM([tx])
if prb.ispaired:
prb.unpair()
if survey.ispaired:
survey.unpair()
prb.pair(survey)
In [16]:
bz_sea = survey.dpred(np.log(sig_sea))
In [17]:
bz_half = survey.dpred(np.log(sig_half))
In [18]:
Bz_sea = bz_sea.reshape((x.size, time.size), order='F')
Bz_half = bz_half.reshape((x.size, time.size), order='F')
In [19]:
TimeX = Utils.ndgrid(x, time)
X = TimeX[:,0].reshape((x.size, time.size), order='F')
Time = TimeX[:,1].reshape((x.size, time.size), order='F')
In [20]:
print ('$10^{%3.2f}$') % (3.0)
In [26]:
test = (np.arange(7)*0.5-4).tolist()
dummy = []
for itest in test:
dummy.append('10$^{'+str(itest+3)+'}$')
In [26]:
In [27]:
ms = 1e3
fig, ax = plt.subplots(1,1, figsize = (8,6))
dat = ax.contourf(X, np.log10(Time), (Bz_sea/Bz_half), 100, cmap = cm.RdPu) #, vmin=1, vmax=1.14, clim=(1, 1.14))
ax.set_xlabel('Distance (m)', fontsize = 16)
ax.set_ylabel('Time (ms)', fontsize = 16)
ax.set_yticklabels(dummy)
cb = colorbar(dat, ax = ax)
levels = np.r_[100., 200., 350.]
CS = ax.contour(X, np.log10(Time), (Bz_sea/Bz_half), levels, colors='k')
ax.clabel(CS, inline=1, fmt = '%3.1f')
# cb = colorbar(dat, ax = ax, format='%5.2f', ticks=np.linspace(1, 1.14, 5))
cb.set_label('Amplitude ratio', fontsize = 14)
fig.savefig('figures/ampratio_time.png', dpi=200)
In [24]:
fig, ax = plt.subplots(1,1, figsize = (8,6))
dat = ax.contourf(X, np.log10(Time), (Bz_sea), 100) #, vmin=1, vmax=1.14, clim=(1, 1.14))
ax.set_xlabel('Distance (m)', fontsize = 16)
ax.set_ylabel('Time (s)', fontsize = 16)
ax.set_yticklabels(dummy)
cb = colorbar(dat, ax = ax)
# cb = colorbar(dat, ax = ax, format='%5.2f', ticks=np.linspace(1, 1.14, 5))
cb.set_label('bz (T)', fontsize = 14)
In [25]:
fig, ax = plt.subplots(1,1, figsize = (8,6))
dat = ax.contourf(X, np.log10(Time), (Bz_half), 100) #, vmin=1, vmax=1.14, clim=(1, 1.14))
ax.set_xlabel('Distance (m)', fontsize = 16)
ax.set_ylabel('Time (s)', fontsize = 16)
ax.set_yticklabels(dummy)
cb = colorbar(dat, ax = ax)
# cb = colorbar(dat, ax = ax, format='%5.2f', ticks=np.linspace(1, 1.14, 5))
cb.set_label('Amplitude ratio', fontsize = 14)
Discussions: