diff --git a/notebooks/MT3D derivatives test.ipynb b/notebooks/MT3D derivatives test.ipynb index 15335015..91f11d1b 100644 --- a/notebooks/MT3D derivatives test.ipynb +++ b/notebooks/MT3D derivatives test.ipynb @@ -44,11 +44,11 @@ "==================== checkDerivative ====================\n", "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", "---------------------------------------------------------\n", - " 0 1.00e-01 3.058e+11 3.238e+10 nan\n", - " 1 1.00e-02 3.096e+10 3.246e+08 1.999\n", - " 2 1.00e-03 3.102e+09 3.251e+06 1.999\n", + " 0 1.00e-01 3.111e+11 3.481e+10 nan\n", + " 1 1.00e-02 3.112e+10 3.421e+08 2.008\n", + " 2 1.00e-03 3.116e+09 3.420e+06 2.000\n", "========================= PASS! =========================\n", - "Happy little convergence test!\n", + "Once upon a time, a happy little test passed.\n", "\n" ] }, @@ -103,11 +103,11 @@ "==================== checkDerivative ====================\n", "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", "---------------------------------------------------------\n", - " 0 1.00e-01 3.151e+11 3.512e+10 nan\n", - " 1 1.00e-02 3.103e+10 3.483e+08 2.004\n", - " 2 1.00e-03 3.101e+09 3.486e+06 2.000\n", + " 0 1.00e-01 3.187e+11 3.398e+10 nan\n", + " 1 1.00e-02 3.176e+10 3.369e+08 2.004\n", + " 2 1.00e-03 3.178e+09 3.372e+06 2.000\n", "========================= PASS! =========================\n", - "Yay passed!\n", + "That was easy!\n", "\n" ] }, @@ -164,11 +164,11 @@ "==================== checkDerivative ====================\n", "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", "---------------------------------------------------------\n", - " 0 1.00e-01 2.848e+10 2.846e+09 nan\n", - " 1 1.00e-02 2.837e+09 2.801e+07 2.007\n", - " 2 1.00e-03 2.838e+08 2.800e+05 2.000\n", + " 0 1.00e-01 2.873e+10 3.211e+09 nan\n", + " 1 1.00e-02 2.879e+09 3.204e+07 2.001\n", + " 2 1.00e-03 2.882e+08 3.207e+05 2.000\n", "========================= PASS! =========================\n", - "That was easy!\n", + "You are awesome.\n", "\n" ] } @@ -200,64 +200,6 @@ "metadata": { "collapsed": false }, - "outputs": [], - "source": [ - "src = survey.srcList[0]\n", - "u = problem.fields(problem.curModel)\n", - "u_src = u[src,:]\n", - "w = np.random.rand(problem.mesh.nC)\n" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "dA_dm = problem.getADeriv_m(src.freq,u_src,w)" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "(20536, 2)" - ] - }, - "execution_count": 8, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "dA_dm.shape" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "a, b, sig, d, e = t3Dmt.halfSpace(.01)" - ] - }, - { - "cell_type": "code", - "execution_count": 21, - "metadata": { - "collapsed": false - }, "outputs": [ { "name": "stdout", @@ -266,26 +208,28 @@ "==================== checkDerivative ====================\n", "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", "---------------------------------------------------------\n", - " 0 1.00e-01 3.699e-07 6.864e-07 nan\n", - " 1 1.00e-02 3.693e-08 6.856e-08 1.001\n", - " 2 1.00e-03 3.693e-09 6.855e-09 1.000\n", - " 3 1.00e-04 3.693e-10 6.855e-10 1.000\n", - " 4 1.00e-05 3.693e-11 6.855e-11 1.000\n", - " 5 1.00e-06 3.693e-12 6.855e-12 1.000\n", - " 6 1.00e-07 3.693e-13 6.855e-13 1.000\n", - " 7 1.00e-08 3.693e-14 6.855e-14 1.000\n", - "*********************************************************\n", - "<<<<<<<<<<<<<<<<<<<<<<<<< FAIL! >>>>>>>>>>>>>>>>>>>>>>>>>\n", - "*********************************************************\n", - "You had so much promise Gudni, oh well...\n", + " 0 1.00e-01 1.436e+01 1.792e-11 nan\n", + " 1 1.00e-02 1.436e+00 1.792e-11 0.000\n", + " 2 1.00e-03 1.436e-01 1.730e-11 0.015\n", + "========================= PASS! =========================\n", + "That was easy!\n", + "\n", + "==================== checkDerivative ====================\n", + "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", + "---------------------------------------------------------\n", + " 0 1.00e-01 1.433e+01 1.792e-11 nan\n", + " 1 1.00e-02 1.433e+00 1.776e-11 0.004\n", + " 2 1.00e-03 1.433e-01 1.763e-11 0.003\n", + "========================= PASS! =========================\n", + "The test be workin!\n", "\n" ] } ], "source": [ - "def testDeriv_ProjFields(prb):\n", + "def testDeriv_epx(prb):\n", " # Initate things for the derivs Test\n", - " src = survey.srcList[0]\n", + " src = prb.survey.srcList[0]\n", " rx = src.rxList[0]\n", " rx.locs = np.array([[0.,0.,0,],[1.,1.,1.]])\n", " u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n", @@ -303,10 +247,379 @@ " f[src,'e_pxSolution'] = u[:len(u)/2]\n", " f[src,'e_pySolution'] = u[len(u)/2::]\n", " # f[src,'b_1d'] = -(m1d.nodalGrad*u)/(1j*simpegem.Utils.EMUtils.omega(src.freq))\n", - " return rx.projectFields(src,survey.mesh,f), lambda t: rx.projectFieldsDeriv(src,survey.mesh,f0,t)\n", - " simpeg.Tests.checkDerivative(fun,u0,num=8,plotIt=False)\n", + " return f._e_px(f[src,'e_pxSolution'],[src]).ravel(), lambda t: f._e_pxDeriv_u([src],t)[:len(u0)/2].ravel()\n", + " simpeg.Tests.checkDerivative(fun,u0,num=3,plotIt=False)\n", + " \n", + "def testDeriv_epy(prb):\n", + " # Initate things for the derivs Test\n", + " src = prb.survey.srcList[0]\n", + " rx = src.rxList[0]\n", + " rx.locs = np.array([[0.,0.,0,],[1.,1.,1.]])\n", + " u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n", + " u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n", + "# u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))\n", + " u0 = np.r_[u0x, u0y]\n", + " f0 = prb.fieldsPair(survey.mesh,survey)\n", + " # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))\n", + " f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x\n", + " f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y\n", + " # f0[src,'b_1d'] = -1/(1j*simpegem.Utils.EMUtils.omega(src.freq))*m1d.nodalGrad*u0\n", + " # Run a test\n", + " def fun(u):\n", + " f = prb.fieldsPair(survey.mesh,survey)\n", + " f[src,'e_pxSolution'] = u[:len(u)/2]\n", + " f[src,'e_pySolution'] = u[len(u)/2::]\n", + " return f._e_py(f[src,'e_pySolution'],[src]).ravel(), lambda t: f._e_pyDeriv_u([src],t).ravel()\n", + " simpeg.Tests.checkDerivative(fun,u0,num=3,plotIt=False)\n", + " \n", + "survey, problem = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.random(.1),comp='zxyr',singleFreq=.01)\n", + "testDeriv_epx(problem)\n", + "testDeriv_epy(problem)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "==================== checkDerivative ====================\n", + "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", + "---------------------------------------------------------\n", + " 0 1.00e-01 3.995e-01 1.107e-13 nan\n", + " 1 1.00e-02 3.995e-02 1.105e-13 0.001\n", + " 2 1.00e-03 3.995e-03 1.110e-13 -0.002\n", + "========================= PASS! =========================\n", + "The test be workin!\n", + "\n", + "==================== checkDerivative ====================\n", + "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", + "---------------------------------------------------------\n", + " 0 1.00e-01 4.075e-01 1.086e-13 nan\n", + " 1 1.00e-02 4.075e-02 1.049e-13 0.015\n", + " 2 1.00e-03 4.075e-03 1.025e-13 0.010\n", + "========================= PASS! =========================\n", + "Happy little convergence test!\n", + "\n" + ] + } + ], + "source": [ + "def testDeriv_bpx(prb):\n", + " # Initate things for the derivs Test\n", + " src = prb.survey.srcList[0]\n", + " rx = src.rxList[0]\n", + " rx.locs = np.array([[0.,0.,0,],[1.,1.,1.]])\n", + " u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n", + " u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n", + "# u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))\n", + " u0 = np.r_[u0x, u0y]\n", + " f0 = prb.fieldsPair(survey.mesh,survey)\n", + " # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))\n", + " f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x\n", + " f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y\n", + " # f0[src,'b_1d'] = -1/(1j*simpegem.Utils.EMUtils.omega(src.freq))*m1d.nodalGrad*u0\n", + " # Run a test\n", + " def fun(u):\n", + " f = prb.fieldsPair(survey.mesh,survey)\n", + " f[src,'e_pxSolution'] = u[:len(u)/2]\n", + " f[src,'e_pySolution'] = u[len(u)/2::]\n", + " # f[src,'b_1d'] = -(m1d.nodalGrad*u)/(1j*simpegem.Utils.EMUtils.omega(src.freq))\n", + " return f._b_px(f[src,'e_pxSolution'],[src]).ravel(), lambda t: f._b_pxDeriv_u(src,t).ravel()\n", + " simpeg.Tests.checkDerivative(fun,u0,num=3,plotIt=False)\n", + " \n", + "def testDeriv_bpy(prb):\n", + " # Initate things for the derivs Test\n", + " src = prb.survey.srcList[0]\n", + " rx = src.rxList[0]\n", + " rx.locs = np.array([[0.,0.,0,],[1.,1.,1.]])\n", + " u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n", + " u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n", + "# u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))\n", + " u0 = np.r_[u0x, u0y]\n", + " f0 = prb.fieldsPair(survey.mesh,survey)\n", + " # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))\n", + " f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x\n", + " f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y\n", + " # f0[src,'b_1d'] = -1/(1j*simpegem.Utils.EMUtils.omega(src.freq))*m1d.nodalGrad*u0\n", + " # Run a test\n", + " def fun(u):\n", + " f = prb.fieldsPair(survey.mesh,survey)\n", + " f[src,'e_pxSolution'] = u[:len(u)/2]\n", + " f[src,'e_pySolution'] = u[len(u)/2::]\n", + " return f._b_py(f[src,'e_pySolution'],[src]).ravel(), lambda t: f._b_pyDeriv_u(src,t).ravel()\n", + " simpeg.Tests.checkDerivative(fun,u0,num=3,plotIt=False)\n", + " \n", "survey, problem = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.random(.1),comp='zyxr',singleFreq=.01)\n", - "testDeriv_ProjFields(problem)" + "testDeriv_bpx(problem)\n", + "testDeriv_bpy(problem)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# %debug" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "==================== checkDerivative ====================\n", + "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", + "---------------------------------------------------------\n", + " 0 1.00e-01 5.819e+11 7.675e+10 nan\n", + " 1 1.00e-02 5.440e+10 8.852e+08 1.938\n", + " 2 1.00e-03 5.403e+09 4.494e+07 1.294\n", + " 3 1.00e-04 5.399e+08 4.428e+06 1.006\n", + " 4 1.00e-05 5.399e+07 4.428e+05 1.000\n", + " 5 1.00e-06 5.399e+06 4.428e+04 1.000\n", + "*********************************************************\n", + "<<<<<<<<<<<<<<<<<<<<<<<<< FAIL! >>>>>>>>>>>>>>>>>>>>>>>>>\n", + "*********************************************************\n", + "You break it, you fix it.\n", + "\n" + ] + } + ], + "source": [ + "from scipy.constants import mu_0\n", + "def testDeriv_Hd(prb):\n", + " src = prb.survey.srcList[0]\n", + " rx = src.rxList[0]\n", + " rx.locs = np.array([[0.,0.,0,],[1.,1.,1.]])\n", + " u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n", + " u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n", + "# u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))\n", + " u0 = np.r_[u0x, u0y]\n", + " f0 = prb.fieldsPair(survey.mesh,survey)\n", + " # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))\n", + " f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x\n", + " f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y\n", + " # Run a testdef testDeriv_ZijN(rx):\n", + " def getHdcomp(rx,f,v=None):\n", + " if rx.locs.ndim == 3:\n", + " eFLocs = rx.locs[:,:,0]\n", + " bFLocs = rx.locs[:,:,1]\n", + " else:\n", + " eFLocs = rx.locs\n", + " bFLocs = rx.locs\n", + " # Get the projection\n", + " Pbx = prb.mesh.getInterpolationMat(bFLocs,'Fx')\n", + " Pby = prb.mesh.getInterpolationMat(bFLocs,'Fy')\n", + " # Get the fields at location\n", + " # px: x-polaration and py: y-polaration.\n", + " hx_px = Pbx*f[src,'b_px']/mu_0\n", + " hy_px = Pby*f[src,'b_px']/mu_0\n", + " hx_py = Pbx*f[src,'b_py']/mu_0\n", + " hy_py = Pby*f[src,'b_py']/mu_0\n", + " # Derivatives as lambda functions\n", + " spPe = simpeg.Utils.spzeros(rx.nD,prb.mesh.nE)\n", + " spPb = simpeg.Utils.spzeros(rx.nD,prb.mesh.nF)\n", + "\n", + " # NOTE: Think b_p?Deriv_u should return a 2*nF size matrix\n", + " hx_px_u = lambda vec: Pbx*f._b_pxDeriv_u(src,vec)/mu_0\n", + " hy_px_u = lambda vec: Pby*f._b_pxDeriv_u(src,vec)/mu_0\n", + " hx_py_u = lambda vec: Pbx*f._b_pyDeriv_u(src,vec)/mu_0\n", + " hy_py_u = lambda vec: Pby*f._b_pyDeriv_u(src,vec)/mu_0\n", + " # Update the input vector\n", + " sdiag = lambda t: simpeg.Utils.sdiag(simpeg.mkvc(t,2))\n", + " # Define the components of the derivative\n", + " if v is not None:\n", + " return (sdiag(hy_py)*hx_px_u(v)) + (sdiag(hx_px)*hy_py_u(v)) - (sdiag(hx_py)*hy_px_u(v)) - (sdiag(hy_px)*hx_py_u(v))\n", + " else:\n", + " return (sdiag(hx_px)*hy_py) - (sdiag(hx_py)*hy_px)\n", + " def fun(u):\n", + " f = prb.fieldsPair(survey.mesh,survey)\n", + " f[src,'e_pxSolution'] = u[:len(u)/2]\n", + " f[src,'e_pySolution'] = u[len(u)/2::]\n", + " return getHdcomp(rx,f), lambda t: getHdcomp(rx,f0,t)\n", + " simpeg.Tests.checkDerivative(fun,u0,num=6,plotIt=False)\n", + "survey, problem = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.random(.01),comp='zxyr',singleFreq=.0001)\n", + "testDeriv_Hd(problem)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false, + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "==================== checkDerivative ====================\n", + "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", + "---------------------------------------------------------\n", + " 0 1.00e-01 1.733e+05 7.391e+02 nan\n", + " 1 1.00e-02 1.734e+04 6.892e+01 1.030\n", + " 2 1.00e-03 1.734e+03 6.896e+00 1.000\n", + " 3 1.00e-04 1.734e+02 6.896e-01 1.000\n", + " 4 1.00e-05 1.734e+01 6.897e-02 1.000\n", + " 5 1.00e-06 1.734e+00 6.897e-03 1.000\n", + "*********************************************************\n", + "<<<<<<<<<<<<<<<<<<<<<<<<< FAIL! >>>>>>>>>>>>>>>>>>>>>>>>>\n", + "*********************************************************\n", + "You break it, you fix it.\n", + "\n", + "==================== checkDerivative ====================\n", + "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", + "---------------------------------------------------------\n", + " 0 1.00e-01 3.075e+05 3.072e+03 nan\n", + " 1 1.00e-02 3.057e+04 1.202e+02 1.407\n", + " 2 1.00e-03 3.055e+03 1.169e+01 1.012\n", + " 3 1.00e-04 3.055e+02 1.168e+00 1.000\n", + " 4 1.00e-05 3.055e+01 1.168e-01 1.000\n", + " 5 1.00e-06 3.055e+00 1.168e-02 1.000\n", + "*********************************************************\n", + "<<<<<<<<<<<<<<<<<<<<<<<<< FAIL! >>>>>>>>>>>>>>>>>>>>>>>>>\n", + "*********************************************************\n", + "Did you put your clever trousers on today?\n", + "\n", + "==================== checkDerivative ====================\n", + "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", + "---------------------------------------------------------\n", + " 0 1.00e-01 4.726e+05 1.575e+04 nan\n", + " 1 1.00e-02 4.791e+04 1.217e+03 1.112\n", + " 2 1.00e-03 4.797e+03 1.216e+02 1.001\n", + " 3 1.00e-04 4.798e+02 1.216e+01 1.000\n", + " 4 1.00e-05 4.798e+01 1.216e+00 1.000\n", + " 5 1.00e-06 4.798e+00 1.216e-01 1.000\n", + "*********************************************************\n", + "<<<<<<<<<<<<<<<<<<<<<<<<< FAIL! >>>>>>>>>>>>>>>>>>>>>>>>>\n", + "*********************************************************\n", + "Coffee break?\n", + "\n", + "==================== checkDerivative ====================\n", + "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", + "---------------------------------------------------------\n", + " 0 1.00e-01 4.323e+05 2.689e+04 nan\n", + " 1 1.00e-02 4.326e+04 4.536e+02 1.773\n", + " 2 1.00e-03 4.326e+03 3.671e+01 1.092\n", + " 3 1.00e-04 4.327e+02 3.660e+00 1.001\n", + " 4 1.00e-05 4.327e+01 3.660e-01 1.000\n", + " 5 1.00e-06 4.327e+00 3.660e-02 1.000\n", + "*********************************************************\n", + "<<<<<<<<<<<<<<<<<<<<<<<<< FAIL! >>>>>>>>>>>>>>>>>>>>>>>>>\n", + "*********************************************************\n", + "No gold star for you.\n", + "\n" + ] + } + ], + "source": [ + "def testDeriv_ZijN(prb):\n", + " src = prb.survey.srcList[0]\n", + " rx = src.rxList[0]\n", + " rx.locs = np.array([[0.,0.,0,],[1.,1.,1.]])\n", + " u0x = np.random.randn(survey.mesh.nE)+np.random.randn(prb.mesh.nE)*1j\n", + " u0y = np.random.randn(survey.mesh.nE)+np.random.randn(prb.mesh.nE)*1j\n", + "# u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))\n", + " u0 = np.r_[u0x, u0y]\n", + " f0 = prb.fieldsPair(prb.mesh,prb.survey)\n", + " # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))\n", + " f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x\n", + " f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y\n", + " def getZijNcomp(rx,f,v=None):\n", + " if rx.locs.ndim == 3:\n", + " eFLocs = rx.locs[:,:,0]\n", + " bFLocs = rx.locs[:,:,1]\n", + " else:\n", + " eFLocs = rx.locs\n", + " bFLocs = rx.locs\n", + " # Get the projection\n", + " Pex = prb.mesh.getInterpolationMat(eFLocs,'Ex')\n", + " Pey = prb.mesh.getInterpolationMat(eFLocs,'Ey')\n", + " Pbx = prb.mesh.getInterpolationMat(bFLocs,'Fx')\n", + " Pby = prb.mesh.getInterpolationMat(bFLocs,'Fy')\n", + " # Get the fields at location\n", + " # px: x-polaration and py: y-polaration.\n", + " ex_px = Pex*f[src,'e_px']\n", + " ey_px = Pey*f[src,'e_px']\n", + " ex_py = Pex*f[src,'e_py']\n", + " ey_py = Pey*f[src,'e_py']\n", + " hx_px = Pbx*f[src,'b_px']/mu_0\n", + " hy_px = Pby*f[src,'b_px']/mu_0\n", + " hx_py = Pbx*f[src,'b_py']/mu_0\n", + " hy_py = Pby*f[src,'b_py']/mu_0\n", + " # Derivatives as lambda functions\n", + " # The size of the deratives should be nD,nU\n", + " ex_px_u = lambda vec: Pex*f._e_pxDeriv_u(src,vec)\n", + " ey_px_u = lambda vec: Pey*f._e_pxDeriv_u(src,vec)\n", + " ex_py_u = lambda vec: Pex*f._e_pyDeriv_u(src,vec)\n", + " ey_py_u = lambda vec: Pey*f._e_pyDeriv_u(src,vec)\n", + " # NOTE: Think b_p?Deriv_u should return a 2*nF size matrix\n", + " hx_px_u = lambda vec: Pbx*f._b_pxDeriv_u(src,vec)/mu_0\n", + " hy_px_u = lambda vec: Pby*f._b_pxDeriv_u(src,vec)/mu_0\n", + " hx_py_u = lambda vec: Pbx*f._b_pyDeriv_u(src,vec)/mu_0\n", + " hy_py_u = lambda vec: Pby*f._b_pyDeriv_u(src,vec)/mu_0\n", + "\n", + " # Update the input vector\n", + " # Define the components of the derivative# Calculate components\n", + " # Update the input vector\n", + " sdiag = lambda t: simpeg.Utils.sdiag(simpeg.mkvc(t,2))\n", + " # Define the components of the derivative\n", + " if 'zxx' in rx.rxType:\n", + " if v is not None:\n", + " return sdiag(hy_py)*ex_px_u(v) + sdiag(ex_px)*hy_py_u(v) - sdiag(ex_py)*hy_px_u(v) - sdiag(hy_px)*ex_py_u(v)\n", + " return ( sdiag(ex_px)*hy_py - sdiag(ex_py)*hy_px)\n", + " elif 'zxy' in rx.rxType:\n", + " if v is not None:\n", + " return -sdiag(hx_py)*ex_px_u(v) - sdiag(ex_px)*hx_py_u(v) + sdiag(ex_py)*hx_px_u(v) + sdiag(hx_px)*ex_py_u(v)\n", + " return (-sdiag(ex_px)*hx_py + sdiag(ex_py)*hx_px)\n", + " elif 'zyx' in rx.rxType:\n", + " if v is not None:\n", + " return ey_px_u(v)*hy_py + ey_px*hy_py_u(v) - ey_py*hy_px_u(v) - ey_py_u(v)*hy_px\n", + " return ( ey_px*hy_py - ey_py*hy_px)\n", + " elif 'zyy' in rx.rxType:\n", + " if v is not None:\n", + " return -ey_px_u(v)*hx_py - ey_px*hx_py_u(v) + ey_py*hx_px_u(v) + ey_py_u(v)*hx_px\n", + " return (-ey_px*hx_py + ey_py*hx_px)\n", + " \n", + " def fun(u):\n", + " f = prb.fieldsPair(survey.mesh,prb.survey)\n", + " f[src,'e_pxSolution'] = u[:len(u)/2]\n", + " f[src,'e_pySolution'] = u[len(u)/2::]\n", + " return getZijNcomp(rx,f), lambda t: getZijNcomp(rx,f0,t)\n", + " simpeg.Tests.checkDerivative(fun,u0,num=6,plotIt=False)\n", + " \n", + "surveyxx, problemxx = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.random(.01),comp='zxxr',singleFreq=.001)\n", + "testDeriv_ZijN(problemxx)\n", + "surveyxy, problemxy = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.random(.01),comp='zxyr',singleFreq=.0001)\n", + "testDeriv_ZijN(problemxy)\n", + "surveyyx, problemyx = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.random(.01),comp='zyxr',singleFreq=.0001)\n", + "testDeriv_ZijN(problemyx)\n", + "surveyyy, problemyy = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.random(.01),comp='zyyr',singleFreq=.0001)\n", + "testDeriv_ZijN(problemyy)" ] }, { @@ -316,6 +629,127 @@ "collapsed": false }, "outputs": [], + "source": [ + "# %debug" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# dA_dm = problem.getADeriv_m(src.freq,u_src,w)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# dA_dm.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# a, b, sig, d, e = t3Dmt.halfSpace(.01)" + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "==================== checkDerivative ====================\n", + "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", + "---------------------------------------------------------\n", + " 0 1.00e-01 5.050e-06 1.906e-09 nan\n", + " 1 1.00e-02 5.049e-07 1.949e-11 1.990\n", + " 2 1.00e-03 5.048e-08 1.969e-13 1.996\n", + " 3 1.00e-04 5.048e-09 1.377e-15 2.155\n", + " 4 1.00e-05 5.048e-10 1.073e-15 0.108\n", + " 5 1.00e-06 5.049e-11 4.179e-15 -0.590\n", + " 6 1.00e-07 5.059e-12 1.079e-14 -0.412\n", + " 7 1.00e-08 5.092e-13 4.411e-15 0.388\n", + "========================= PASS! =========================\n", + "Yay passed!\n", + "\n" + ] + } + ], + "source": [ + "def testDeriv_ProjFields(prb):\n", + " # Initate things for the derivs Test\n", + " src = prb.survey.srcList[0]\n", + " rx = src.rxList[0]\n", + "# rx.locs = np.array([[0.,0.,0,],[1.,1.,1.]])\n", + "# u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n", + "# u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n", + "# u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))\n", + "# u0 = np.r_[u0x, u0y]\n", + "# f0 = prb.fieldsPair(survey.mesh,survey)\n", + " # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))\n", + "# f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x\n", + "# f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y\n", + "\n", + "# prb.mapping = simpeg.Maps.ExpMap(prb.mesh)\n", + "# if True:\n", + "# x0 = x0 + np.random.randn(prb.mesh.nC)*cond*1e-1 \n", + " survey = prb.survey\n", + " src = survey.srcList[0]\n", + " \n", + " f0 = prb.fields(prb.curModel)\n", + " u0 = np.r_[f0[src,'e_pxSolution'],f0[src,'e_pySolution']]\n", + "# def fun(x):\n", + "# prb.curModel = x\n", + " # f0[src,'b_1d'] = -1/(1j*simpegem.Utils.EMUtils.omega(src.freq))*m1d.nodalGrad*u0\n", + " # Run a test\n", + " def fun(u):\n", + " f = prb.fieldsPair(survey.mesh,survey)\n", + " f[src,'e_pxSolution'] = u[:len(u)/2]\n", + " f[src,'e_pySolution'] = u[len(u)/2::]\n", + " # f[src,'b_1d'] = -(m1d.nodalGrad*u)/(1j*simpegem.Utils.EMUtils.omega(src.freq))\n", + " return rx.projectFields(src,survey.mesh,f), lambda t: rx.projectFieldsDeriv(src,survey.mesh,f0,t)\n", + " simpeg.Tests.checkDerivative(fun,u0,num=8,plotIt=False)\n", + "survey, problem = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.halfSpace(1),comp='zxyr',singleFreq=.1)\n", + "testDeriv_ProjFields(problem)" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# %debug" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": false + }, + "outputs": [], "source": [ "src = survey.srcList[0]\n", "rx = src.rxList[0]" @@ -323,7 +757,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 18, "metadata": { "collapsed": false }, @@ -334,16 +768,30 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 30, "metadata": { "collapsed": false }, - "outputs": [], - "source": [] + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 19.30432593, 8.59301637, 7.75725432, ..., 2.18228759,\n", + " 16.72231401, 7.67029076])" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "problem.curModel.sigma" + ] }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 19, "metadata": { "collapsed": false }, @@ -361,7 +809,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 20, "metadata": { "collapsed": false }, @@ -369,47 +817,89 @@ { "data": { "text/plain": [ - "(2, 1)" + "(20536, 1)" ] }, - "execution_count": 17, + "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "rx.projectFieldsDeriv(src,survey.mesh,f0,u0).shape" + "f0._e_px(f0[src,'e_pxSolution'],[src]).shape" ] }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 21, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(0, 1)" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "f0._e_pxDeriv_u([src],u0)[len(u0)/2::].shape" + ] + }, + { + "cell_type": "code", + "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [], "source": [ - "simpeg.Utils.spzeros?" + "simpeg.Utils.spzeros" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 23, "metadata": { "collapsed": false }, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 0.00031425],\n", + " [ 0.00030875]])" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "rx.projectFields(src,survey.mesh,f0)" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 24, "metadata": { "collapsed": false }, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "ERROR: No traceback has been produced, nothing to debug.\n" + ] + } + ], "source": [ "%debug" ] @@ -429,6 +919,18 @@ "display_name": "Python 2", "language": "python", "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.10" } }, "nbformat": 4,