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,:]


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

            python setup.py build_ext --inplace
    

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]:
array([ nan+nanj])

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


10000 62831.8530718 (0.628318548187-0.628318513249j)
[[  0. -0.j   0. +0.j]
 [  0. +0.j  inf+infj]]
[  0. +0.j  nan+nanj]
5000 62831.8530718 (0.628318548187-0.628318513249j)
[[  0. -0.j   0. +0.j]
 [  0. +0.j  inf+infj]]
[  0. +0.j  nan+nanj]
1000 62831.8530718 (0.628318548187-0.628318513249j)
[[  1.33271357e-273 -2.32814399e-278j   0.00000000e+000 +0.00000000e+000j]
 [  0.00000000e+000 +0.00000000e+000j   7.50348781e+272 +1.31079930e+268j]]
[  1.60872758e-273 -2.81162844e-278j  -1.55402321e+272 -2.70737840e+267j]
500 62831.8530718 (0.628318548187-0.628318513249j)
[[  3.65063497e-137 -3.18868363e-142j   0.00000000e+000 +0.00000000e+000j]
 [  0.00000000e+000 +0.00000000e+000j   2.73924950e+136 +2.39262488e+131j]]
[  4.40670622e-137 -3.85267017e-142j  -5.67317146e+135 -4.92836188e+130j]
100 62831.8530718 (0.628318548187-0.628318513249j)
[[  5.15790907e-28 -9.01045454e-34j   0.00000000e+00 +0.00000000e+00j]
 [  0.00000000e+00 +0.00000000e+00j   1.93877012e+27 +3.38687631e+21j]]
[  6.22614702e-28 -1.09272824e-33j  -4.01532439e+26 -6.82387175e+20j]
50 62831.8530718 (0.628318548187-0.628318513249j)
[[  2.27110305e-14 -1.98371768e-20j   0.00000000e+00 +0.00000000e+00j]
 [  0.00000000e+00 +0.00000000e+00j   4.40314674e+13 +3.84597256e+07j]]
[  2.74146389e-14 -2.41688373e-20j  -9.11921548e+12 -7.53244599e+06j]
10 62831.8530718 (0.628318548187-0.628318513249j)
[[  1.86744306e-03 -3.26227365e-10j   0.00000000e+00 +0.00000000e+00j]
 [  0.00000000e+00 +0.00000000e+00j   5.35491562e+02 +9.35460925e-05j]]
[  2.25420318e-03 -4.12148003e-10j  -1.10903934e+02 -1.41102131e-05j]

Is there a smart way to "fix" this so that the 1D layering can be used for the analytic solution?


In [ ]: