In [3]:
from SimPEG import *
from scipy.constants import mu_0
def omega(freq):
"""Change frequency to angular frequency, omega"""
return 2.*np.pi*freq
In [4]:
%pylab inline
In [5]:
# M = Mesh.TensorMesh([[(100.,32)],[(100.,34)],[(100.,18)]], x0='CCC')
M = Mesh.TensorMesh([[(100.,20)],[(100.,20)],[(100.,18)]], x0='CCC')
In [6]:
print M
In [7]:
# Setup the model
conds = [1e-2,1]
sig = Utils.ModelBuilder.defineBlock(M.gridCC,[0,0,-300],[500,500,-100],conds)
sig[M.gridCC[:,2]>0] = 1e-8
sigBG = np.zeros(M.nC) + conds[0]
sigBG[M.gridCC[:,2]>0] = 1e-8
colorbar(M.plotImage(log10(sig)))
Out[7]:
In [8]:
# Get the mass matrix
# The model
Msig = M.getEdgeInnerProduct(sig)
MsigBG = M.getEdgeInnerProduct(sigBG)
Mmu = M.getFaceInnerProduct(mu_0, invProp=True)
In [9]:
freq = 1e1
C = M.edgeCurl
A = C.T*Mmu*C - 1j*omega(freq)*Msig
ABG = C.T*Mmu*C - 1j*omega(freq)*MsigBG
In [10]:
# Need to solve x and y polarizations of the source.
from simpegMT.Utils import get1DEfields
# Get a 1d solution for a halfspace background
mesh1d = Mesh.TensorMesh([M.hz],np.array([M.x0[2]]))
e0_1d = get1DEfields(mesh1d,M.r(sigBG,'CC','CC','M')[0,0,:],freq)
# Setup x (east) polarization (_x)
ex_x = np.zeros(M.vnEx,dtype=complex)
ey_x = np.zeros((M.nEy,1),dtype=complex)
ez_x = np.zeros((M.nEz,1),dtype=complex)
# Assign the source to ex_x
for i in arange(M.vnEx[0]):
for j in arange(M.vnEx[2]):
ex_x[i,j,:] = e0_1d
eBG_x = np.vstack((Utils.mkvc(M.r(ex_x,'Ex','Ex','V'),2),ey_x,ez_x))
rhs_x = ABG.dot(eBG_x)
In [11]:
M.vnEy
Out[11]:
In [12]:
# Setup y (north) polarization (_y)
ex_y = np.zeros(M.nEx, dtype='complex128')
ey_y = np.zeros((M.vnEy), dtype='complex128')
ez_y = np.zeros(M.nEz, dtype='complex128')
# Assign the source to ex_x
for i in arange(M.vnEy[0]):
for j in arange(M.vnEy[2]):
ey_y[i,j,:] = e0_1d
# eBG_y = np.vstack((ex_y,Utils.mkvc(M.r(ey_y,'Ey','Ey','V'),2),ez_y))
# rhs_y = ABG.dot(eBG_y)
In [13]:
eBG_y = np.r_[ex_y,Utils.mkvc(ey_y),ez_y]
rhs_y = ABG.dot(eBG_y)
We are using splu in scipy package. This is bit slow, but on the cluster you can use mumps, which might a lot faster. We can think about having better iterative solver.
In [14]:
%%time
# Solve the systems for each polarization
Ainv = SolverLU(A)
In [23]:
%%time
e_x = Ainv*rhs_x
I want to visualize electrical field, which is a vector, so I average them on to cell center. Also I want to see current density ($\vec{j} = \sigma \vec{e}$).
In [34]:
Meinv = M.getEdgeInnerProduct(np.ones_like(sig), invMat=True)
In [35]:
j_x = Meinv*Msig*e_x
In [36]:
e_x_CC = M.aveE2CCV*e_x
j_x_CC = M.aveE2CCV*j_x
# j_x_CC = Utils.sdiag(np.r_[sig, sig, sig])*e_x_CC
Then use "plotSlice" function, to visualize 2D sections
In [37]:
fig, ax = plt.subplots(1,2, figsize = (12, 5))
dat0 = M.plotSlice(abs(e_x_CC), vType='CCv', view='vec', streamOpts={'color': 'k'}, normal='X', ax = ax[0])
cb0 = plt.colorbar(dat0[0], ax = ax[0])
dat1 = M.plotSlice(abs(j_x_CC), vType='CCv', view='vec', streamOpts={'color': 'k'}, normal='X', ax = ax[1])
cb1 = plt.colorbar(dat1[0], ax = ax[1])
Is it reasonable?: Based on that you put resistive target that makes sense to me; current does not want to flow on resistive target so they just do roundabout:). And see air interface. It is continuous on current but not on electric field, which looks reasonable.
In [33]:
# # Get a 1d solution for a halfspace background
# mesh1d = simpeg.Mesh.TensorMesh(M.hz,M.x0[2])
# e0_1d = get1DEfields(mesh1d,M.r(sigBG,'CC','CC','M')[0,0,:],freq)
# # Setup x (east) polarization (_x)
# ex_x = np.zeros(M.vnEx)
# ey_x = np.zeros(M.vnEy)
# ez_x = np.zeros(M.vnEz)
# # Assign the source to ex_x
# for i in arange(M.vnEx[0]):
# for j in arange(M.vnEx[2]):
# ex_x[i,j,:] = e0_1d
# eBG_x = np.vstack(M.r(ex_x,'Ex','Ex','V'),M.r(ex_y)