MTanalytic1D_layerIssue.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:4f51688cd2ee8a11dad3df1928925d3c9cad0da43a3f6a3c3c840024caae5fe1"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Issues with padding cells and high frequencies in the analytic MT layered earth."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import SimPEG as simpeg\n",
      "elev = 300\n",
      "# 3D mesh and model\n",
      "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'])\n",
      "conds = [1,1e-2]\n",
      "sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,[-10000,-10000,-200],[10000,10000,0],conds)\n",
      "sig[M.gridCC[:,2]>elev] = 1e-8\n",
      "sig[M.gridCC[:,2]<-600] = 1e-1\n",
      "# Make the 1D mesh and model\n",
      "mesh1d = simpeg.Mesh.TensorMesh([M.hz],np.array([M.x0[2]]))\n",
      "sig1D = M.r(sig,'CC','CC','M')[0,0,:]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Efficiency Warning: Interpolation will be slow, use setup.py!\n",
        "\n",
        "            python setup.py build_ext --inplace\n",
        "    \n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Run for high frequency\n",
      "freq = 1e4"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Run the analytic problem\n",
      "import simpegMT as simpegmt\n",
      "anaEd, anaEu, anaHd, anaHu = simpegmt.Utils.MT1Danalytic.getEHfields(mesh1d,sig1D,freq,np.array([300]))\n",
      "anaE = anaEd+anaEu\n",
      "anaH = anaHd+anaHu\n",
      "anaZ = anaE/anaH\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "anaZ"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "array([ nan+nanj])"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Returns nan because in the analytic solution the propagation of the fields in the layer \"blows\" up."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sig = 10\n",
      "sig0 = 20\n",
      "mu = 4*np.pi*1e-7\n",
      "eps = 8.85*1e-12\n",
      "for h in [10000,5000,1000,500,100,50,10]:\n",
      "    w = 2*np.pi*freq\n",
      "    k0 = np.sqrt(eps*mu*w**2-1j*mu*sig0*w)\n",
      "    k = np.sqrt(eps*mu*w**2-1j*mu*sig*w)\n",
      "    zp = (w*mu)/k\n",
      "    yp1 = k0/(w*mu)\n",
      "    # Convert fields to down/up going components in layer below current layer\n",
      "    Pj1 = np.array([[1,1],[yp1,-yp1]])\n",
      "    # Convert fields to down/up going components in current layer\n",
      "    Pjinv = 1./2*np.array([[1,zp],[1,-zp]])\n",
      "    # Propagate down and up components through the current layer\n",
      "    elamh = np.array([[np.exp(-1j*k*h),0],[0,np.exp(1j*k*h)]])\n",
      "    UD = elamh.dot(Pjinv.dot(Pj1)).dot([1,0])\n",
      "    print h, w, k \n",
      "    print elamh\n",
      "    #print Pj1, Pjinv, elamh\n",
      "    print UD"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 62831.8530718 (0.628318548187-0.628318513249j)\n",
        "[[  0. -0.j   0. +0.j]\n",
        " [  0. +0.j  inf+infj]]\n",
        "[  0. +0.j  nan+nanj]\n",
        "5000 62831.8530718 (0.628318548187-0.628318513249j)\n",
        "[[  0. -0.j   0. +0.j]\n",
        " [  0. +0.j  inf+infj]]\n",
        "[  0. +0.j  nan+nanj]\n",
        "1000 62831.8530718 (0.628318548187-0.628318513249j)\n",
        "[[  1.33271357e-273 -2.32814399e-278j   0.00000000e+000 +0.00000000e+000j]\n",
        " [  0.00000000e+000 +0.00000000e+000j   7.50348781e+272 +1.31079930e+268j]]\n",
        "[  1.60872758e-273 -2.81162844e-278j  -1.55402321e+272 -2.70737840e+267j]\n",
        "500 62831.8530718 (0.628318548187-0.628318513249j)\n",
        "[[  3.65063497e-137 -3.18868363e-142j   0.00000000e+000 +0.00000000e+000j]\n",
        " [  0.00000000e+000 +0.00000000e+000j   2.73924950e+136 +2.39262488e+131j]]\n",
        "[  4.40670622e-137 -3.85267017e-142j  -5.67317146e+135 -4.92836188e+130j]\n",
        "100 62831.8530718 (0.628318548187-0.628318513249j)\n",
        "[[  5.15790907e-28 -9.01045454e-34j   0.00000000e+00 +0.00000000e+00j]\n",
        " [  0.00000000e+00 +0.00000000e+00j   1.93877012e+27 +3.38687631e+21j]]\n",
        "[  6.22614702e-28 -1.09272824e-33j  -4.01532439e+26 -6.82387175e+20j]\n",
        "50 62831.8530718 (0.628318548187-0.628318513249j)\n",
        "[[  2.27110305e-14 -1.98371768e-20j   0.00000000e+00 +0.00000000e+00j]\n",
        " [  0.00000000e+00 +0.00000000e+00j   4.40314674e+13 +3.84597256e+07j]]\n",
        "[  2.74146389e-14 -2.41688373e-20j  -9.11921548e+12 -7.53244599e+06j]\n",
        "10 62831.8530718 (0.628318548187-0.628318513249j)\n",
        "[[  1.86744306e-03 -3.26227365e-10j   0.00000000e+00 +0.00000000e+00j]\n",
        " [  0.00000000e+00 +0.00000000e+00j   5.35491562e+02 +9.35460925e-05j]]\n",
        "[  2.25420318e-03 -4.12148003e-10j  -1.10903934e+02 -1.41102131e-05j]\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Is there a smart way to \"fix\" this so that the 1D layering can be used for the analytic solution?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}