From 960cb0a3e57f616db141b9bc04d2bd5d534569dc Mon Sep 17 00:00:00 2001 From: GudniRos Date: Mon, 2 Nov 2015 11:23:06 -0800 Subject: [PATCH] Added a Derivative test notebook --- notebooks/MT3D derivatives test.ipynb | 436 ++++++++++++++++++++++++++ simpegMT/SurveyMT.py | 34 +- 2 files changed, 453 insertions(+), 17 deletions(-) create mode 100644 notebooks/MT3D derivatives test.ipynb diff --git a/notebooks/MT3D derivatives test.ipynb b/notebooks/MT3D derivatives test.ipynb new file mode 100644 index 00000000..15335015 --- /dev/null +++ b/notebooks/MT3D derivatives test.ipynb @@ -0,0 +1,436 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import SimPEG as simpeg\n", + "import simpegEM as simpegem, simpegMT as simpegmt\n", + "from SimPEG.Utils import meshTensor\n", + "import numpy as np\n", + "import simpegMT.Tests.test_Problem3D_againstAnalytic as t3Dmt" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "sigmaHalf = 0.01\n", + "survey, problem = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.halfSpace(sigmaHalf),comp='zxyi',singleFreq=1)\n", + "if False:\n", + " fields = problem.fields(problem.curModel.sigma)\n", + " data = survey.dpred(problem.curModel.sigma)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "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.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", + "========================= PASS! =========================\n", + "Happy little convergence test!\n", + "\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "def testDeriv_dA_dm(prb,cond):\n", + " TOL = 1e-4\n", + " FLR = 1e-20\n", + " x0 = np.log(np.ones(prb.mesh.nC)*cond)\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", + " freq = src.freq\n", + " v1 = np.random.randn(prb.mesh.nE,1)\n", + " v2 = np.random.randn(prb.mesh.nE,1)\n", + " v = np.hstack(( simpeg.mkvc(v1,2), simpeg.mkvc(v2,2)))\n", + " u_0 = prb.fieldsPair(prb.mesh,prb.survey)\n", + " u_0[src,'e_pxSolution'] = v1\n", + " u_0[src,'e_pySolution'] = v2\n", + " def fun(x):\n", + " prb.curModel = x\n", + " A = prb.getA(freq) #\n", + "# return simpeg.mkvc(A*v1)+simpeg.mkvc(A*v2), lambda t: simpeg.mkvc(prb.getADeriv_m(freq, u_0[src], t))\n", + " return A*v, lambda t: (prb.getADeriv_m(freq, u_0[src], t))\n", + " return simpeg.Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)\n", + "testDeriv_dA_dm(problem,0.1)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "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.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", + "========================= PASS! =========================\n", + "Yay passed!\n", + "\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "def testDeriv_dRHS_dm(prb,cond):\n", + " TOL = 1e-4\n", + " FLR = 1e-20\n", + " x0 = np.log(np.ones(prb.mesh.nC)*cond)\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", + " u0 = prb.fields(x0)\n", + " freq = src.freq\n", + " A = prb.getA(freq) #\n", + " A_I = prb.Solver(A, **prb.solverOpts)\n", + "\n", + " ftype = prb._fieldType + 'Solution'\n", + " u0_src = u0[src, ftype]\n", + " v = np.random.randn(prb.mesh.nE,1)\n", + " \n", + " def fun(x):\n", + " prb.curModel = x\n", + " return simpeg.mkvc(np.sum(prb.getRHS(freq))), lambda x: simpeg.mkvc(prb.getRHSDeriv_m(freq, x))\n", + "# return simpeg.mkvc(prb.fields(x)[src,prb._fieldType + 'Solution']), lambda x: simpeg.mkvc(prb.getADeriv_m(freq, u0_src, x))\n", + " return simpeg.Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)\n", + "testDeriv_dA_dm(problem,0.1)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "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 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", + "========================= PASS! =========================\n", + "That was easy!\n", + "\n" + ] + } + ], + "source": [ + "def testDeriv_S_e(prb,cond):\n", + " # Initate things for the derivs Test\n", + " TOL = 1e-4\n", + " FLR = 1e-20\n", + " \n", + " src = prb.survey.srcList[0]\n", + " rx = src.rxList[0]\n", + "\n", + " x0 = np.log(np.ones(prb.mesh.nC)*cond)\n", + "# prb.mapping = simpeg.Maps.ExpMap(prb.mesh)\n", + " if True:\n", + " x0 = x0 + np.random.randn(prb.mesh.nC)*cond*1e-1 \n", + " def fun(x):\n", + " prb.curModel = x\n", + " return src.S_e(prb), lambda t: src.S_eDeriv_m(prb,t)\n", + " simpeg.Tests.checkDerivative(fun,x0,num=3,plotIt=False)\n", + "survey, problem = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.random(1),comp='All',singleFreq=1)\n", + "testDeriv_S_e(problem,0.1)" + ] + }, + { + "cell_type": "code", + "execution_count": 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", + "output_type": "stream", + "text": [ + "==================== 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", + "\n" + ] + } + ], + "source": [ + "def testDeriv_ProjFields(prb):\n", + " # Initate things for the derivs Test\n", + " src = 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 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.random(.1),comp='zyxr',singleFreq=.01)\n", + "testDeriv_ProjFields(problem)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "src = survey.srcList[0]\n", + "rx = src.rxList[0]" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "rx.locs = np.array([[0, 0, 0],[1,1,1]])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "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", + "src = survey.srcList[0]\n", + "rx = src.rxList[0]\n", + "f0 = problem.fieldsPair(survey.mesh,survey)\n", + "u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))\n", + "f0[src,'e_pxSolution'] = u0x\n", + "f0[src,'e_pySolution'] = u0y" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(2, 1)" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "rx.projectFieldsDeriv(src,survey.mesh,f0,u0).shape" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "simpeg.Utils.spzeros?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "rx.projectFields(src,survey.mesh,f0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "%debug" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/simpegMT/SurveyMT.py b/simpegMT/SurveyMT.py index f76cac0c..3a13b172 100644 --- a/simpegMT/SurveyMT.py +++ b/simpegMT/SurveyMT.py @@ -188,31 +188,31 @@ class RxMT(Survey.BaseRx): Pby = mesh.getInterpolationMat(bFLocs,'Fy') # Get the fields at location # px: x-polaration and py: y-polaration. - ex_px = Utils.sdiag(mkvc(Pex*f[src,'e_px'],2)) - ey_px = Utils.sdiag(mkvc(Pey*f[src,'e_px'],2)) - ex_py = Utils.sdiag(mkvc(Pex*f[src,'e_py'],2)) - ey_py = Utils.sdiag(mkvc(Pey*f[src,'e_py'],2)) - hx_px = Utils.sdiag(mkvc(Pbx*f[src,'b_px']/mu_0,2)) - hy_px = Utils.sdiag(mkvc(Pby*f[src,'b_px']/mu_0,2)) - hx_py = Utils.sdiag(mkvc(Pbx*f[src,'b_py']/mu_0,2)) - hy_py = Utils.sdiag(mkvc(Pby*f[src,'b_py']/mu_0,2)) + ex_px = Pex*f[src,'e_px'] + ey_px = Pey*f[src,'e_px'] + ex_py = Pex*f[src,'e_py'] + ey_py = Pey*f[src,'e_py'] + hx_px = Pbx*f[src,'b_px']/mu_0 + hy_px = Pby*f[src,'b_px']/mu_0 + hx_py = Pbx*f[src,'b_py']/mu_0 + hy_py = Pby*f[src,'b_py']/mu_0 # Derivatives as lambda functions spPe = Utils.spzeros(self.nD,mesh.nE) spPb = Utils.spzeros(self.nD,mesh.nF) # The size of the diratives should be nD,nU - ex_px_u = lambda vec: Utils.sdiag(mkvc(sp.hstack((Pex,spPe))*f._e_pxDeriv_u(src,vec),2)) - ey_px_u = lambda vec: Utils.sdiag(mkvc(sp.hstack((Pey,spPe))*f._e_pxDeriv_u(src,vec),2)) - ex_py_u = lambda vec: Utils.sdiag(mkvc(sp.hstack((spPe,Pex))*f._e_pyDeriv_u(src,vec),2)) - ey_py_u = lambda vec: Utils.sdiag(mkvc(sp.hstack((spPe,Pey))*f._e_pyDeriv_u(src,vec),2)) + ex_px_u = lambda vec: sp.hstack((Pex,spPe))*f._e_pxDeriv_u(src,vec) + ey_px_u = lambda vec: sp.hstack((Pey,spPe))*f._e_pxDeriv_u(src,vec) + ex_py_u = lambda vec: sp.hstack((spPe,Pex))*f._e_pyDeriv_u(src,vec) + ey_py_u = lambda vec: sp.hstack((spPe,Pey))*f._e_pyDeriv_u(src,vec) # NOTE: Think b_p?Deriv_u should return a 2*nF size matrix - hx_px_u = lambda vec: Utils.sdiag(mkvc(sp.hstack((Pbx,spPb))*f._b_pxDeriv_u(src,vec)/mu_0,2)) - hy_px_u = lambda vec: Utils.sdiag(mkvc(sp.hstack((Pby,spPb))*f._b_pxDeriv_u(src,vec)/mu_0,2)) - hx_py_u = lambda vec: Utils.sdiag(mkvc(sp.hstack((spPb,Pbx))*f._b_pyDeriv_u(src,vec)/mu_0,2)) - hy_py_u = lambda vec: Utils.sdiag(mkvc(sp.hstack((spPb,Pby))*f._b_pyDeriv_u(src,vec)/mu_0,2)) + hx_px_u = lambda vec: sp.hstack((Pbx,spPb))*f._b_pxDeriv_u(src,vec)/mu_0 + hy_px_u = lambda vec: sp.hstack((Pby,spPb))*f._b_pxDeriv_u(src,vec)/mu_0 + hx_py_u = lambda vec: sp.hstack((spPb,Pbx))*f._b_pyDeriv_u(src,vec)/mu_0 + hy_py_u = lambda vec: sp.hstack((spPb,Pby))*f._b_pyDeriv_u(src,vec)/mu_0 # Update the input vector # v = mkvc(v,2) # Make v into a column vector # Define the components of the derivative - Hd = Utils.sdiag(mkvc(1./(hx_px*hy_py - hx_py*hy_px).data,2)) + Hd = 1./(hx_px*hy_py - hx_py*hy_px) Hd_uV = hx_px_u(v)*hy_py + hx_px*hy_py_u(v) - hx_py*hy_px_u(v) - hx_py_u(v)*hy_px # Calculate components if 'zxx' in self.rxType: