Issues with padding cells and high frequencies in the analytic MT layered earth.
In [1]:
import SimPEG as simpeg
elev = 300
# 3D mesh and model
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'])
conds = [1,1e-2]
sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,[-10000,-10000,-200],[10000,10000,0],conds)
sig[M.gridCC[:,2]>elev] = 1e-8
sig[M.gridCC[:,2]<-600] = 1e-1
# Make the 1D mesh and model
mesh1d = simpeg.Mesh.TensorMesh([M.hz],np.array([M.x0[2]]))
sig1D = M.r(sig,'CC','CC','M')[0,0,:]
In [3]:
# Run for high frequency
freq = 1e4
In [4]:
# Run the analytic problem
import simpegMT as simpegmt
anaEd, anaEu, anaHd, anaHu = simpegmt.Utils.MT1Danalytic.getEHfields(mesh1d,sig1D,freq,np.array([300]))
anaE = anaEd+anaEu
anaH = anaHd+anaHu
anaZ = anaE/anaH
In [5]:
anaZ
Out[5]:
Returns nan because in the analytic solution the propagation of the fields in the layer "blows" up.
In [11]:
sig = 10
sig0 = 20
mu = 4*np.pi*1e-7
eps = 8.85*1e-12
for h in [10000,5000,1000,500,100,50,10]:
w = 2*np.pi*freq
k0 = np.sqrt(eps*mu*w**2-1j*mu*sig0*w)
k = np.sqrt(eps*mu*w**2-1j*mu*sig*w)
zp = (w*mu)/k
yp1 = k0/(w*mu)
# Convert fields to down/up going components in layer below current layer
Pj1 = np.array([[1,1],[yp1,-yp1]])
# Convert fields to down/up going components in current layer
Pjinv = 1./2*np.array([[1,zp],[1,-zp]])
# Propagate down and up components through the current layer
elamh = np.array([[np.exp(-1j*k*h),0],[0,np.exp(1j*k*h)]])
UD = elamh.dot(Pjinv.dot(Pj1)).dot([1,0])
print h, w, k
print elamh
#print Pj1, Pjinv, elamh
print UD
In [ ]: