In [1]:
import SimPEG as simpeg
from scipy.constants import mu_0
def omega(freq):
"""Change frequency to angular frequency, omega"""
return 2.*np.pi*freq
In [2]:
%pylab inline
In [3]:
# M = Mesh.TensorMesh([[(100.,32)],[(100.,34)],[(100.,18)]], x0='CCC')
M = simpeg.Mesh.TensorMesh([[(100,5,-1.5),(100.,5),(100,5,1.5)],[(100,5,-1.5),(100.,5),(100,5,1.5)],[(100,10,-1.5),(100.,10),(100,10,1.5)]], x0=['C','C','C'])
In [4]:
print M.vectorNz
In [5]:
# Set parameters
freq = 10
conds = [0.01]
elev = 300
# Setup the models
sig = np.zeros(M.nC) + 1e-8
sig[M.gridCC[:,2]<=elev] = conds[0]
# sig[M.gridCC[:,2]<-600] = 1e-1
sigBG = np.zeros(M.nC) + 1e-8
sigBG[M.gridCC[:,2]<=elev] = conds[0]
In [6]:
# Plot the models
colorbar(M.plotImage(log10(sig)))
colorbar(M.plotImage(log10(sigBG)))
Out[6]:
In [7]:
# Get the mass matrix
# The model
Msig = M.getEdgeInnerProduct(sig)
MsigBG = M.getEdgeInnerProduct(sigBG)
Mmu = M.getFaceInnerProduct(mu_0, invProp=True)
In [8]:
# Form the A matrices
C = M.edgeCurl
A = C.T*Mmu*C - 1j*omega(freq)*Msig
ABG = C.T*Mmu*C - 1j*omega(freq)*MsigBG
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 [9]:
%%time
# Solve the systems for each polarization
Ainv = simpeg.SolverLU(A)
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 = simpeg.Mesh.TensorMesh([M.hz],np.array([M.x0[2]]))
e0_1d = get1DEfields(mesh1d,M.r(sigBG,'CC','CC','M')[0,0,:],freq,sourceAmp=1).conj() # conjugate to comply with phase behavior
# 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[1]):
ex_x[i,j,:] = -e0_1d #Negative to comply with phase orientation.
eBG_x = np.vstack((simpeg.Utils.mkvc(ex_x,2),ey_x,ez_x))
# Note 100% sure why this has to be negative.
rhs_x = -ABG*eBG_x
In [11]:
# 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[1]):
ey_y[i,j,:] = e0_1d
In [12]:
eBG_y = np.r_[ex_y,simpeg.Utils.mkvc(ey_y),ez_y]
# Note 100% sure why this has to be negative.
rhs_y = -ABG*eBG_y
In [14]:
%%time
# Solve each polarization
e_x = Ainv*rhs_x
e_y = Ainv*rhs_y
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 [15]:
Meinv = M.getEdgeInnerProduct(np.ones_like(sig), invMat=True)
In [16]:
j_x = Meinv*Msig*e_x
j_y = Meinv*Msig*e_y
In [17]:
e_x_CC = M.aveE2CCV*e_x
e_y_CC = M.aveE2CCV*e_y
j_x_CC = M.aveE2CCV*j_x
j_y_CC = M.aveE2CCV*j_y
# j_x_CC = Utils.sdiag(np.r_[sig, sig, sig])*e_x_CC
In [31]:
e_x
Out[31]:
Then use "plotSlice" function, to visualize 2D sections
In [18]:
fig, ax = plt.subplots(1,2, figsize = (12, 5))
dat0 = M.plotSlice(abs(e_y_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_y_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 [19]:
# Calculate the data
rx_x, rx_y = np.meshgrid(np.arange(-250,251,50),np.arange(-250,251,50))
rx_loc = np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),elev+np.zeros((np.prod(rx_x.shape),1))))
# Get the projection matrices
Qex = M.getInterpolationMat(rx_loc,'Ex')
Qey = M.getInterpolationMat(rx_loc,'Ey')
Qez = M.getInterpolationMat(rx_loc,'Ez')
Qfx = M.getInterpolationMat(rx_loc,'Fx')
Qfy = M.getInterpolationMat(rx_loc,'Fy')
Qfz = M.getInterpolationMat(rx_loc,'Fz')
In [20]:
e_x_loc = np.hstack([simpeg.Utils.mkvc(Qex*e_x,2),simpeg.Utils.mkvc(Qey*e_x,2),simpeg.Utils.mkvc(Qez*e_x,2)])
e_y_loc = np.hstack([simpeg.Utils.mkvc(Qex*e_y,2),simpeg.Utils.mkvc(Qey*e_y,2),simpeg.Utils.mkvc(Qez*e_y,2)])
Ciw = -C/(1j*omega(freq))
b_x_loc = np.hstack([simpeg.Utils.mkvc(Qfx*Ciw*e_x,2),simpeg.Utils.mkvc(Qfy*Ciw*e_x,2),simpeg.Utils.mkvc(Qfz*Ciw*e_x,2)])
b_y_loc = np.hstack([simpeg.Utils.mkvc(Qfx*Ciw*e_y,2),simpeg.Utils.mkvc(Qfy*Ciw*e_y,2),simpeg.Utils.mkvc(Qfz*Ciw*e_y,2)])
In [21]:
from scipy.constants import mu_0
# Make a combined matrix
dt = np.dtype([('ex1',complex),('ey1',complex),('ez1',complex),('hx1',complex),('hy1',complex),('hz1',complex),('ex2',complex),('ey2',complex),('ez2',complex),('hx2',complex),('hy2',complex),('hz2',complex)])
combMat = np.empty((len(e_x_loc)),dtype=dt)
combMat['ex1'] = e_x_loc[:,0]
combMat['ey1'] = e_x_loc[:,1]
combMat['ez1'] = e_x_loc[:,2]
combMat['ex2'] = e_y_loc[:,0]
combMat['ey2'] = e_y_loc[:,1]
combMat['ez2'] = e_y_loc[:,2]
combMat['hx1'] = b_x_loc[:,0]/mu_0
combMat['hy1'] = b_x_loc[:,1]/mu_0
combMat['hz1'] = b_x_loc[:,2]/mu_0
combMat['hx2'] = b_y_loc[:,0]/mu_0
combMat['hy2'] = b_y_loc[:,1]/mu_0
combMat['hz2'] = b_y_loc[:,2]/mu_0
In [22]:
def calculateImpedance(fieldsData):
'''
Function that calculates MT impedance data from a rec array with E and H field data from both polarizations
'''
zxx = (fieldsData['ex1']*fieldsData['hy2'] - fieldsData['ex2']*fieldsData['hy1'])/(fieldsData['hx1']*fieldsData['hy2'] - fieldsData['hx2']*fieldsData['hy1'])
zxy = (-fieldsData['ex1']*fieldsData['hx2'] + fieldsData['ex2']*fieldsData['hx1'])/(fieldsData['hx1']*fieldsData['hy2'] - fieldsData['hx2']*fieldsData['hy1'])
zyx = (fieldsData['ey1']*fieldsData['hy2'] - fieldsData['ey2']*fieldsData['hy1'])/(fieldsData['hx1']*fieldsData['hy2'] - fieldsData['hx2']*fieldsData['hy1'])
zyy = (-fieldsData['ey1']*fieldsData['hx2'] + fieldsData['ey2']*fieldsData['hx1'])/(fieldsData['hx1']*fieldsData['hy2'] - fieldsData['hx2']*fieldsData['hy1'])
return zxx, zxy, zyx, zyy
zxx, zxy, zyx, zyy = calculateImpedance(combMat)
In [29]:
zxy
Out[29]:
In [28]:
ind = np.where(np.sum(np.power(rx_loc - np.array([0,0,elev]),2),axis=1)< 5)
def appResPhs(freq,z):
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
return app_res, app_phs
print appResPhs(freq,zxy[ind])
print appResPhs(freq,zyx[ind])
In [25]:
e0_1d
h0_1dC = -(mesh1d.nodalGrad*e0_1d)/(1j*omega(freq)*mu_0)
h0_1d = mesh1d.getInterpolationMat(mesh1d.vectorNx,'CC')*h0_1dC
print e0_1d, h0_1d, appResPhs(freq,e0_1d/h0_1d)
In [26]:
# Get the analytic solution for the halfspace.
import simpegMT as simpegmt
sig1DBG = M.r(sigBG,'CC','CC','M')[0,0,:]
anaEd, anaEu, anaHd, anaHu = simpegmt.Utils.MT1Danalytic.getEHfields(mesh1d,sig1DBG,freq,mesh1d.vectorNx)
anaEtemp = anaEd+anaEu
anaHtemp = anaHd+anaHu
# Scale the solution
anaE = (anaEtemp/anaEtemp[-1])#.conj()
anaH = (anaHtemp/anaEtemp[-1])#.conj()
In [27]:
anaZ = anaE/anaH
print anaZ, appResPhs(freq,anaZ)
In [27]:
In [27]: