diff --git a/notebooks/Derivative test MT1D.ipynb b/notebooks/Derivative test MT1D.ipynb index 41b0c6d2..751e9879 100644 --- a/notebooks/Derivative test MT1D.ipynb +++ b/notebooks/Derivative test MT1D.ipynb @@ -84,12 +84,13 @@ " srcList.append(simpegmt.SurveyMT.srcMT_polxy_1DhomotD(rxList,freq))\n", "else:\n", " for freq in freqs:\n", - " srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq,sigma))\n", + " srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))\n", "# Make the survey\n", "survey = simpegmt.SurveyMT.SurveyMT(srcList)\n", "\n", "# Set the problem\n", "problem = simpegmt.ProblemMT1D.eForm_psField(m1d)\n", + "problem.sigmaPrimary = sigma\n", "problem.pair(survey)\n", "\n", "# Get the fields\n", @@ -101,12 +102,22 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": { "collapsed": false }, - "outputs": [], - "source": [] + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "ERROR: No traceback has been produced, nothing to debug.\n" + ] + } + ], + "source": [ + "%debug\n" + ] }, { "cell_type": "markdown", @@ -136,7 +147,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 5, "metadata": { "collapsed": true }, @@ -168,7 +179,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 6, "metadata": { "collapsed": true }, @@ -192,7 +203,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 7, "metadata": { "collapsed": false }, @@ -216,7 +227,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 8, "metadata": { "collapsed": false, "scrolled": true @@ -229,13 +240,13 @@ "==================== checkDerivative ====================\n", "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", "---------------------------------------------------------\n", - " 0 1.00e-01 1.884e-05 1.227e-07 nan\n", - " 1 1.00e-02 1.873e-06 1.265e-09 1.987\n", - " 2 1.00e-03 1.872e-07 1.269e-11 1.999\n", - " 3 1.00e-04 1.872e-08 1.269e-13 2.000\n", - " 4 1.00e-05 1.872e-09 1.269e-15 2.000\n", + " 0 1.00e-01 3.089e-04 2.730e-05 nan\n", + " 1 1.00e-02 2.841e-05 2.582e-07 2.024\n", + " 2 1.00e-03 2.818e-06 2.568e-09 2.002\n", + " 3 1.00e-04 2.816e-07 2.567e-11 2.000\n", + " 4 1.00e-05 2.816e-08 2.567e-13 2.000\n", "========================= PASS! =========================\n", - "That was easy!\n", + "The test be workin!\n", "\n" ] }, @@ -245,7 +256,7 @@ "True" ] }, - "execution_count": 7, + "execution_count": 8, "metadata": {}, "output_type": "execute_result" } @@ -261,7 +272,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 9, "metadata": { "collapsed": false }, @@ -269,10 +280,10 @@ { "data": { "text/plain": [ - "array([ 0.00052762])" + "array([-0.00285083])" ] }, - "execution_count": 8, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" } @@ -283,7 +294,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 10, "metadata": { "collapsed": false }, @@ -291,10 +302,10 @@ { "data": { "text/plain": [ - "array([[ 0.00124017]])" + "array([[-0.00130618]])" ] }, - "execution_count": 9, + "execution_count": 10, "metadata": {}, "output_type": "execute_result" } @@ -305,7 +316,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 11, "metadata": { "collapsed": false }, @@ -316,7 +327,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 12, "metadata": { "collapsed": false }, @@ -328,12 +339,12 @@ "==================== checkDerivative ====================\n", "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", "---------------------------------------------------------\n", - " 0 1.00e-01 4.417e-08 4.873e-09 nan\n", - " 1 1.00e-02 4.132e-09 4.832e-11 2.004\n", - " 2 1.00e-03 4.105e-10 4.828e-13 2.000\n", - " 3 1.00e-04 4.103e-11 4.827e-15 2.000\n", + " 0 1.00e-01 1.537e-05 1.005e-07 nan\n", + " 1 1.00e-02 1.541e-06 1.008e-09 1.999\n", + " 2 1.00e-03 1.541e-07 1.008e-11 2.000\n", + " 3 1.00e-04 1.541e-08 1.008e-13 2.000\n", "========================= PASS! =========================\n", - "You get a gold star!\n", + "The test be workin!\n", "\n" ] }, @@ -343,7 +354,7 @@ "True" ] }, - "execution_count": 11, + "execution_count": 12, "metadata": {}, "output_type": "execute_result" } @@ -368,7 +379,26 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 13, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "ERROR: No traceback has been produced, nothing to debug.\n" + ] + } + ], + "source": [ + "%debug" + ] + }, + { + "cell_type": "code", + "execution_count": 14, "metadata": { "collapsed": true }, @@ -388,7 +418,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 15, "metadata": { "collapsed": false }, @@ -398,7 +428,7 @@ "output_type": "stream", "text": [ "Adjoint e formulation - projectFieldsDeriv\n", - "-2.26989762698e-05 -2.26989762698e-05 0.0 1e-08 True\n" + "0.000177264079929 0.000177264079929 3.52365706058e-19 1e-07 True\n" ] }, { @@ -407,7 +437,7 @@ "True" ] }, - "execution_count": 13, + "execution_count": 15, "metadata": {}, "output_type": "execute_result" } @@ -439,7 +469,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 21, "metadata": { "collapsed": false }, @@ -448,19 +478,25 @@ "name": "stdout", "output_type": "stream", "text": [ - "Adjoint test e formulation - getADeriv_m\n", - "(-1977540.36505+2093781.70221j) (-1977540.36505+2093781.70221j) (-1.86264514923e-09+2.79396772385e-09j) 100.0 True\n" + "Adjoint test e formulation - getADeriv_m\n" ] }, { - "data": { - "text/plain": [ - "True" - ] - }, - "execution_count": 14, - "metadata": {}, - "output_type": "execute_result" + "ename": "KeyError", + "evalue": "'Source does not have a uid: e_1dSolution'", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mKeyError\u001b[0m Traceback (most recent call last)", + "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[0;32m 20\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[0mvJw\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mwJtv\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvJw\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mwJtv\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtol\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvJw\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mwJtv\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m<\u001b[0m \u001b[0mtol\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 21\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvJw\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mwJtv\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m<\u001b[0m \u001b[0mtol\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 22\u001b[1;33m \u001b[0mgetADeriv_mAdjointTest\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[1;32m\u001b[0m in \u001b[0;36mgetADeriv_mAdjointTest\u001b[1;34m()\u001b[0m\n\u001b[0;32m 15\u001b[0m \u001b[0mw\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrandn\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm1d\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnC\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;31m#+np.random.randn(m1d.nN)*1j\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 16\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 17\u001b[1;33m \u001b[0mvJw\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mv\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mproblem\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetADeriv_m\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mf0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mw\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 18\u001b[0m \u001b[0mwJtv\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mw\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mproblem\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetADeriv_m\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mf0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0madjoint\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 19\u001b[0m \u001b[0mtol\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmax\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mTOL\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m10\u001b[0m\u001b[1;33m**\u001b[0m\u001b[0mint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlog10\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvJw\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mFLR\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/ProblemMT1D/Problems.py\u001b[0m in \u001b[0;36mgetADeriv_m\u001b[1;34m(self, freq, u, v, adjoint)\u001b[0m\n\u001b[0;32m 67\u001b[0m \u001b[0mMeMui\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmesh\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetEdgeInnerProduct\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1.0\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mmu_0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 68\u001b[0m \u001b[1;31m#\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 69\u001b[1;33m \u001b[0mu_src\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mu\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'e_1dSolution'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 70\u001b[0m \u001b[0mdMf_dsig\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmesh\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetFaceInnerProductDeriv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcurModel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msigma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mu_src\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcurModel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msigmaDeriv\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 71\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0madjoint\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Fields.pyc\u001b[0m in \u001b[0;36m__getitem__\u001b[1;34m(self, key)\u001b[0m\n\u001b[0;32m 150\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 151\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m__getitem__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkey\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 152\u001b[1;33m \u001b[0mind\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mname\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_indexAndNameFromKey\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'get'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 153\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mname\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 154\u001b[0m \u001b[0mout\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m{\u001b[0m\u001b[1;33m}\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Fields.pyc\u001b[0m in \u001b[0;36m_indexAndNameFromKey\u001b[1;34m(self, key, accessType)\u001b[0m\n\u001b[0;32m 131\u001b[0m \u001b[0msrcTestList\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mname\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mkey\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 132\u001b[0m \u001b[0mname\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_nameIndex\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maccessType\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 133\u001b[1;33m \u001b[0mind\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_srcIndex\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrcTestList\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 134\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mind\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mname\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 135\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Fields.pyc\u001b[0m in \u001b[0;36m_srcIndex\u001b[1;34m(self, srcTestList)\u001b[0m\n\u001b[0;32m 100\u001b[0m \u001b[0mind\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msrcTestList\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 101\u001b[0m \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 102\u001b[1;33m \u001b[0mind\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msurvey\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetSourceIndex\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrcTestList\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 103\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mind\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 104\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Survey.pyc\u001b[0m in \u001b[0;36mgetSourceIndex\u001b[1;34m(self, sources)\u001b[0m\n\u001b[0;32m 339\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0msrc\u001b[0m \u001b[1;32min\u001b[0m \u001b[0msources\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 340\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'uid'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 341\u001b[1;33m \u001b[1;32mraise\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Source does not have a uid: %s'\u001b[0m\u001b[1;33m%\u001b[0m\u001b[0mstr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 342\u001b[0m \u001b[0minds\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmap\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;32mlambda\u001b[0m \u001b[0msrc\u001b[0m\u001b[1;33m:\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_sourceOrder\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0muid\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msources\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 343\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mNone\u001b[0m \u001b[1;32min\u001b[0m \u001b[0minds\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[1;31mKeyError\u001b[0m: 'Source does not have a uid: e_1dSolution'" + ] } ], "source": [ @@ -480,8 +516,8 @@ " # print prb.PropMap.PropModel.nP\n", " w = np.random.randn(m1d.nC)#+np.random.randn(m1d.nN)*1j\n", "\n", - " vJw = v.dot(problem.getADeriv_m(freq,u0,w))\n", - " wJtv = w.dot(problem.getADeriv_m(freq,u0,v,adjoint=True))\n", + " vJw = v.dot(problem.getADeriv_m(freq,f0,w))\n", + " wJtv = w.dot(problem.getADeriv_m(freq,f0,v,adjoint=True))\n", " tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) \n", " print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol\n", " return np.abs(vJw - wJtv) < tol\n", @@ -490,16 +526,226 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 17, "metadata": { - "collapsed": true + "collapsed": false }, - "outputs": [], - "source": [] + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "> \u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/ProblemMT1D/Problems.py\u001b[0m(69)\u001b[0;36mgetADeriv_m\u001b[1;34m()\u001b[0m\n", + "\u001b[1;32m 68 \u001b[1;33m \u001b[1;31m#\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[0m\u001b[1;32m---> 69 \u001b[1;33m \u001b[0mu_src\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mu\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'e_1dSolution'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[0m\u001b[1;32m 70 \u001b[1;33m \u001b[0mdMf_dsig\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmesh\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetFaceInnerProductDeriv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcurModel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msigma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mu_src\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcurModel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msigmaDeriv\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[0m\n", + "ipdb> u\n", + "> \u001b[1;32m\u001b[0m(17)\u001b[0;36mgetADeriv_mAdjointTest\u001b[1;34m()\u001b[0m\n", + "\u001b[1;32m 16 \u001b[1;33m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[0m\u001b[1;32m---> 17 \u001b[1;33m \u001b[0mvJw\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mv\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mproblem\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetADeriv_m\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mu0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mw\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[0m\u001b[1;32m 18 \u001b[1;33m \u001b[0mwJtv\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mw\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mproblem\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetADeriv_m\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mu0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0madjoint\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[0m\n", + "ipdb> d\n", + "> \u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/ProblemMT1D/Problems.py\u001b[0m(69)\u001b[0;36mgetADeriv_m\u001b[1;34m()\u001b[0m\n", + "\u001b[1;32m 68 \u001b[1;33m \u001b[1;31m#\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[0m\u001b[1;32m---> 69 \u001b[1;33m \u001b[0mu_src\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mu\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'e_1dSolution'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[0m\u001b[1;32m 70 \u001b[1;33m \u001b[0mdMf_dsig\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmesh\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetFaceInnerProductDeriv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcurModel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msigma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mu_src\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcurModel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msigmaDeriv\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[0m\n", + "ipdb> p u\n", + "array([ 7.09826705e-01 -3.42993400e-01j,\n", + " 2.06495724e-02 -1.13236591e+00j,\n", + " 7.19940065e-01 +6.18299116e-01j,\n", + " -3.76573540e-01 +9.09687373e-01j,\n", + " -5.96166731e-01 +2.23309643e-01j,\n", + " 3.30533731e-01 +7.98068064e-01j,\n", + " -1.41356795e+00 -1.74137927e-01j,\n", + " 1.51554175e+00 -2.16363375e-01j,\n", + " 1.59469627e+00 -1.87971566e+00j,\n", + " -5.54637698e-01 -1.01936458e+00j,\n", + " 1.01869865e+00 +8.33443055e-01j,\n", + " -4.92616671e-01 -1.97036224e+00j,\n", + " 1.00389102e+00 -1.52567261e-01j,\n", + " 7.90254779e-01 -5.68065046e-01j,\n", + " 9.17329964e-01 -1.18879097e+00j,\n", + " -1.38317362e+00 -2.17789273e+00j,\n", + " 2.58324368e-01 -1.51003598e+00j,\n", + " -2.04803634e+00 -1.44299943e-01j,\n", + " 3.33011285e+00 +1.52905384e+00j,\n", + " 2.91429972e-01 -1.91849967e+00j,\n", + " -1.10273784e+00 -8.65794290e-01j,\n", + " -2.80103156e-01 -1.42653192e+00j,\n", + " -9.81404885e-01 +9.06845144e-02j,\n", + " 6.37579045e-01 +3.84148006e-01j,\n", + " -6.83597988e-01 +6.91187965e-01j,\n", + " -8.98575981e-01 +3.80347476e-02j,\n", + " 1.02647319e+00 -1.13876879e-01j,\n", + " 8.28006255e-01 -9.95903406e-01j,\n", + " 3.53931083e-01 -2.80188697e-01j,\n", + " -1.09339690e+00 +1.43646692e+00j,\n", + " 2.59276626e+00 +4.09662561e-01j,\n", + " -5.97873868e-01 -5.23852065e-01j,\n", + " 8.87054287e-01 -1.69872744e+00j,\n", + " 2.29767481e-01 -5.88143928e-01j,\n", + " -6.97610315e-01 +5.68810360e-01j,\n", + " -4.47049016e-01 +3.29354463e-01j,\n", + " -3.31498311e-01 -8.27120276e-01j,\n", + " 1.57689998e-01 +6.68378872e-01j,\n", + " 9.47018631e-02 +2.60180446e-01j,\n", + " -4.59177257e-01 -7.66896049e-01j,\n", + " -1.16791204e+00 +4.20585804e-01j,\n", + " 2.72405147e-01 +3.48079787e-01j,\n", + " -2.45880514e-02 +1.90226778e-01j,\n", + " -7.03428363e-01 +1.20531568e+00j,\n", + " 1.60362217e+00 +6.88472259e-01j,\n", + " -5.49144335e-01 +6.96208292e-01j,\n", + " -9.83200346e-01 -1.15726575e+00j,\n", + " 1.22826524e+00 -8.58070855e-01j,\n", + " -5.99020439e-01 -2.19182124e+00j,\n", + " -5.69992871e-01 +1.68586900e+00j,\n", + " 9.26752644e-01 +1.05767829e+00j,\n", + " 8.64247444e-01 +1.51282356e-01j,\n", + " -6.24547862e-01 -1.54517937e+00j,\n", + " -1.72043526e+00 +1.91659834e+00j,\n", + " -1.84501951e+00 -2.40165623e-01j,\n", + " 1.15210121e+00 -1.19627222e+00j,\n", + " -3.93036395e-01 -1.94918396e+00j,\n", + " -1.18266273e+00 -5.24652722e-01j,\n", + " 2.67461628e-01 -8.07279138e-01j,\n", + " -1.56363641e+00 +1.59257615e-03j,\n", + " -3.17211883e-01 -3.23379690e+00j,\n", + " -1.08658939e-01 -8.66079039e-01j,\n", + " -9.30350013e-01 -1.22668219e+00j,\n", + " -1.07263137e+00 -5.56847000e-02j,\n", + " 1.27878013e+00 -1.34173537e+00j,\n", + " 7.91441492e-01 +6.42111414e-01j,\n", + " -2.60121575e+00 -1.08140886e+00j,\n", + " -1.15027808e+00 -4.47110142e-01j,\n", + " -5.69747800e-01 +2.34065916e-01j,\n", + " 1.52155089e+00 -5.80340775e-01j,\n", + " 2.50052750e-01 -6.58793179e-01j,\n", + " 2.96944135e-01 +1.27455206e+00j,\n", + " 9.18939766e-01 -4.58122761e-01j,\n", + " 4.90015531e-01 -8.86288179e-01j,\n", + " -1.65155083e+00 -6.57733967e-01j,\n", + " 9.72148169e-01 -1.19446820e+00j,\n", + " 2.08365350e+00 -1.53585595e-01j,\n", + " -4.11101072e-01 +1.01418536e+00j,\n", + " 9.72580562e-01 -2.20470840e-01j,\n", + " 1.72922147e-01 -6.86360934e-01j,\n", + " 3.80255034e-01 +1.96862619e-01j,\n", + " -1.53628572e+00 +5.31946011e-01j,\n", + " -4.65939167e-01 +9.11725611e-01j,\n", + " -1.06390365e+00 -3.07658816e-02j,\n", + " -2.40734296e+00 -7.60530882e-01j,\n", + " 1.28242272e+00 +3.32085973e+00j,\n", + " 3.17004751e+00 +1.62654811e+00j,\n", + " -8.49758476e-01 +9.86658359e-01j,\n", + " 4.99578557e-01 +4.01392488e-01j,\n", + " -9.48695998e-01 +4.04322381e-01j,\n", + " 1.75371116e-01 +7.85458813e-01j,\n", + " 4.90232886e-01 -1.39280042e-01j,\n", + " -2.31430252e-01 -8.44065839e-01j,\n", + " -1.02319133e+00 +1.90560784e+00j,\n", + " 1.08940667e+00 +1.63046186e+00j,\n", + " -4.34233324e-01 +2.90439312e-01j,\n", + " 4.47622371e-01 +1.26507976e+00j,\n", + " -2.07317470e+00 -2.21558532e+00j,\n", + " -3.19553487e-01 -9.19210427e-01j,\n", + " 6.09870390e-01 +3.05692898e-01j,\n", + " 4.92081860e-01 -2.18326059e-01j,\n", + " 1.48544627e+00 -2.02054490e-01j,\n", + " -1.59612367e-01 +7.90853777e-02j,\n", + " 2.85854411e-01 +2.16567814e-01j,\n", + " 5.28063071e-01 +2.08838949e+00j,\n", + " 6.17580558e-01 +9.90748539e-01j,\n", + " -2.64539211e-01 +1.16598665e+00j,\n", + " -1.07357542e+00 -7.58002064e-01j,\n", + " 6.13399713e-02 +2.52338250e+00j,\n", + " -4.84180135e-01 -1.45046020e+00j,\n", + " -6.57519374e-01 -1.52719059e+00j,\n", + " -1.54949084e-01 -7.64796563e-01j,\n", + " -2.01172915e-02 +1.07377565e+00j,\n", + " 2.19840636e-01 +6.89622678e-01j,\n", + " -2.64578500e+00 +1.46933966e+00j,\n", + " 1.22168696e+00 +2.33458454e-01j,\n", + " 2.71303177e+00 +2.77858685e-01j,\n", + " 4.39426665e-01 +1.79924695e-01j,\n", + " -1.73858667e+00 -3.03473636e-01j,\n", + " 2.68736300e-01 +8.32359634e-01j,\n", + " 2.57517784e+00 +1.88097576e-01j,\n", + " 2.76775469e-01 -6.49891901e-01j,\n", + " -9.87790301e-01 +4.50328269e-01j,\n", + " 9.19033669e-01 +1.01453748e+00j,\n", + " -1.35706747e+00 +5.90240590e-01j,\n", + " 6.60118738e-01 +5.10361208e-01j,\n", + " -2.17652642e-01 -4.26632365e-01j,\n", + " -1.39163540e+00 +9.12739752e-01j,\n", + " -1.76241829e+00 -1.04835908e+00j,\n", + " -7.28279529e-01 +2.75961750e-01j,\n", + " -1.32043513e+00 -3.28790183e-01j,\n", + " 5.88227093e-01 +9.41865970e-01j,\n", + " 1.97695607e+00 +1.04535450e+00j,\n", + " -3.02599128e+00 +1.02031744e+00j,\n", + " -1.34308763e+00 -4.48222977e-01j,\n", + " -1.86687861e+00 +2.26767568e+00j,\n", + " 5.61472023e-01 -5.28652462e-01j,\n", + " -5.26747946e-01 +2.75545054e-01j,\n", + " 6.16116514e-01 +3.41765793e-01j,\n", + " 1.38709750e+00 +5.87537003e-01j,\n", + " -2.27923203e-01 +7.19317677e-01j,\n", + " -8.66404713e-01 -1.80932081e-01j,\n", + " -4.58146668e-01 +1.75955171e+00j,\n", + " -5.12526647e-01 -1.91939928e-01j,\n", + " -6.26532699e-01 -5.71301355e-01j,\n", + " -3.35566040e-01 -5.06173544e-01j,\n", + " -3.86124497e-01 -8.00560567e-01j,\n", + " 2.47252265e+00 +3.83928843e-01j,\n", + " 3.35561169e-01 +1.19309445e+00j,\n", + " -6.19075444e-01 -1.19021628e+00j,\n", + " -1.17335297e+00 +5.80500783e-01j,\n", + " -2.35111248e-01 -1.59785495e+00j,\n", + " -1.24819583e-01 -8.05281256e-01j,\n", + " -9.18606131e-01 -1.25312332e-01j,\n", + " 3.14400972e-01 +1.18650079e+00j,\n", + " 6.50630049e-01 -6.99853497e-01j,\n", + " 1.28269722e+00 -1.29832850e+00j,\n", + " 2.07880202e-01 +2.56385692e-02j,\n", + " -1.45089461e+00 -9.35093103e-01j,\n", + " 6.84146951e-01 +7.74422758e-01j,\n", + " 9.45439196e-01 -3.50695577e-01j,\n", + " 1.00057580e+00 -1.27412936e+00j,\n", + " -2.01665427e-01 +1.73591867e+00j,\n", + " 1.76623768e+00 -1.42964638e+00j,\n", + " 1.12282181e+00 +9.02767375e-01j,\n", + " -2.22663798e+00 +7.34672692e-02j,\n", + " 8.26230180e-02 -7.46601441e-01j,\n", + " -4.92810641e-01 -5.49544995e-01j,\n", + " 1.06703148e-01 +1.70461156e+00j,\n", + " -6.84949746e-01 -1.46950185e-01j,\n", + " 1.47509281e+00 +8.77421590e-01j,\n", + " -9.83722735e-03 -6.30597166e-01j,\n", + " 2.36006569e-03 +2.12883413e-01j,\n", + " -3.67899785e-01 -9.95508241e-02j,\n", + " 1.26752410e+00 +2.49370098e-01j,\n", + " -1.19439605e+00 +2.94935402e-01j,\n", + " -2.28617888e+00 +5.70564449e-02j,\n", + " 3.95617530e-02 -2.80995512e-01j,\n", + " -7.51027805e-01 -8.14567233e-01j,\n", + " -1.24160227e+00 -4.13553574e-01j, -4.76626111e-01 +4.16480574e-01j])\n", + "ipdb> u_src\n", + "*** NameError: name 'u_src' is not defined\n", + "ipdb> c\n" + ] + } + ], + "source": [ + "%debug" + ] }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 18, "metadata": { "collapsed": false }, @@ -509,7 +755,7 @@ "output_type": "stream", "text": [ "Adjoint test e formulation - getRHSDeriv_m\n", - "(13201.2196403+13827.5790776j) (13201.2196403+13827.5790776j) (-5.45696821064e-12+3.63797880709e-12j) 1.0 True\n" + "(-10077.7224119-28916.4723594j) (-10077.7224119-28916.4723594j) (1.81898940355e-12-3.63797880709e-12j) 1.0 True\n" ] }, { @@ -518,7 +764,7 @@ "True" ] }, - "execution_count": 15, + "execution_count": 18, "metadata": {}, "output_type": "execute_result" } @@ -550,19 +796,11 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": null, "metadata": { "collapsed": false }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "2\n" - ] - } - ], + "outputs": [], "source": [ "simpeg.mkvc(np.random.randn(survey.nD)+np.random.randn(survey.nD)*1j,2)\n", "\n", @@ -571,7 +809,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 19, "metadata": { "collapsed": false }, @@ -602,7 +840,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 20, "metadata": { "collapsed": false }, @@ -612,7 +850,7 @@ "output_type": "stream", "text": [ "Adjoint e formulation - Jvec\n", - "1.96695386678e-05 1.96695386678e-05 3.38813178902e-21 1e-08 True\n" + "1.09508355274e-05 1.09508355274e-05 -1.01643953671e-20 1e-08 True\n" ] }, { @@ -621,7 +859,7 @@ "True" ] }, - "execution_count": 18, + "execution_count": 20, "metadata": {}, "output_type": "execute_result" } @@ -629,6 +867,15 @@ "source": [ "JvecAdjointTest()" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] } ], "metadata": { @@ -636,18 +883,6 @@ "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.9" } }, "nbformat": 4, diff --git a/simpegMT/BaseMT.py b/simpegMT/BaseMT.py index 12261a01..8b31fbe5 100644 --- a/simpegMT/BaseMT.py +++ b/simpegMT/BaseMT.py @@ -92,8 +92,8 @@ class BaseMTProblem(BaseFDEMProblem): # Calculate the projection derivatives for rx in src.rxList: # Get the projection derivative - # v should be of size nE,2 (each column for 2 polarizations) - PDeriv_u = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v) # wrt u, we don't have have PDeriv wrt m + # v should be of size 2*nE (for 2 polarizations) + PDeriv_u = lambda t: rx.projectFieldsDeriv(src, self.mesh, u, t) # wrt u, we don't have have PDeriv wrt m Jv[src, rx] = PDeriv_u(mkvc(du_dm)) # Return the vectorized sensitivities return mkvc(Jv) diff --git a/simpegMT/FieldsMT.py b/simpegMT/FieldsMT.py index 3d2b23c5..2d3a9e1b 100644 --- a/simpegMT/FieldsMT.py +++ b/simpegMT/FieldsMT.py @@ -233,13 +233,13 @@ class FieldsMT_3D(FieldsMT): # NOTE: v needs to be length 2*nE to account for both polarizations def _b_pxSecondaryDeriv_u(self, src, v, adjoint = False): - C = sp.kron(self.mesh.edgeCurl,[[1,0],[0,1]]) + C = sp.kron(self.mesh.edgeCurl,[[1,0],[0,0]]) if adjoint: return - 1./(1j*omega(src.freq)) * (C.T * v) return - 1./(1j*omega(src.freq)) * (C * v) def _b_pySecondaryDeriv_u(self, src, v, adjoint = False): - C = sp.kron(self.mesh.edgeCurl,[[1,0],[0,1]]) + C = sp.kron(self.mesh.edgeCurl,[[0,0],[0,1]]) if adjoint: return - 1./(1j*omega(src.freq)) * (C.T * v) return - 1./(1j*omega(src.freq)) * (C * v) diff --git a/simpegMT/ProblemMT1D/Problems.py b/simpegMT/ProblemMT1D/Problems.py index d7c2cdb0..f237a086 100644 --- a/simpegMT/ProblemMT1D/Problems.py +++ b/simpegMT/ProblemMT1D/Problems.py @@ -65,8 +65,9 @@ class eForm_psField(BaseMTProblem): dsig_dm = self.curModel.sigmaDeriv MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0) - # Need to make the dMf_dsig symmetirc (nN,nN), don't know how to do this - dMf_dsig = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv + # + u_src = u['e_1dSolution'] + dMf_dsig = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv if adjoint: return 1j * omega(freq) * ( dMf_dsig.T * v ) # Note: output has to be nN/nF, not nC/nE. diff --git a/simpegMT/SurveyMT.py b/simpegMT/SurveyMT.py index a99e9212..f76cac0c 100644 --- a/simpegMT/SurveyMT.py +++ b/simpegMT/SurveyMT.py @@ -188,46 +188,48 @@ class RxMT(Survey.BaseRx): Pby = mesh.getInterpolationMat(bFLocs,'Fy') # Get the fields at location # px: x-polaration and py: y-polaration. - ex_px = Pex*mkvc(f[src,'e_px'],2) - ey_px = Pey*mkvc(f[src,'e_px'],2) - ex_py = Pex*mkvc(f[src,'e_py'],2) - ey_py = Pey*mkvc(f[src,'e_py'],2) - hx_px = Pbx*mkvc(f[src,'b_px']/mu_0,2) - hy_px = Pby*mkvc(f[src,'b_px']/mu_0,2) - hx_py = Pbx*mkvc(f[src,'b_py']/mu_0,2) - hy_py = Pby*mkvc(f[src,'b_py']/mu_0,2) + 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)) # Derivatives as lambda functions - ex_px_u = lambda vec: sp.hstack((Pex,Pex))*f._e_pxDeriv_u(src,vec) - ey_px_u = lambda vec: sp.hstack((Pey,Pey))*f._e_pxDeriv_u(src,vec) - ex_py_u = lambda vec: sp.hstack((Pex,Pex))*f._e_pyDeriv_u(src,vec) - ey_py_u = lambda vec: sp.hstack((Pey,Pey))*f._e_pyDeriv_u(src,vec) + 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)) # NOTE: Think b_p?Deriv_u should return a 2*nF size matrix - hx_px_u = lambda vec: sp.hstack((Pbx,Pbx))*f._b_pxDeriv_u(src,vec)/mu_0 - hy_px_u = lambda vec: sp.hstack((Pby,Pby))*f._b_pxDeriv_u(src,vec)/mu_0 - hx_py_u = lambda vec: sp.hstack((Pbx,Pbx))*f._b_pyDeriv_u(src,vec)/mu_0 - hy_py_u = lambda vec: sp.hstack((Pby,Pby))*f._b_pyDeriv_u(src,vec)/mu_0 - + 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)) # Update the input vector - v = mkvc(v,2) # Make v into a column vector + # v = mkvc(v,2) # Make v into a column vector # Define the components of the derivative - Hd = Utils.sdiag(1/(hx_px*hy_py - hx_py*hy_px)) - Hd_uV = hx_px_u(hy_py*v) - hx_py*hy_px_u(v) + hx_px*hy_py_u(v) - hx_py_u(hy_px*v) + Hd = Utils.sdiag(mkvc(1./(hx_px*hy_py - hx_py*hy_px).data,2)) + 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: Zij = ( ex_px*hy_py - ex_py*hy_px)*Hd - ZijN_uV = ex_px_u(hy_py*v) - ex_py*hy_px_u(v) + ex_px*hy_py_u(v) - ex_py_u(hy_px*v) + ZijN_uV = ex_px_u(v)*hy_py + ex_px*hy_py_u(v) - ex_py*hy_px_u(v) - ex_py_u(v)*hy_px elif 'zxy' in self.rxType: Zij = (-ex_px*hx_py + ex_py*hx_px)*Hd - ZijN_uV = -ex_px_u(hx_py*v) + ex_py*hx_px_u(v) - ex_px*hx_py_u(v) + ex_py_u(hx_px*v) + ZijN_uV = -ex_px_u(v)*hx_py - ex_px*hx_py_u(v) + ex_py*hx_px_u(v) + ex_py_u(v)*hx_px elif 'zyx' in self.rxType: Zij = ( ey_px*hy_py - ey_py*hy_px)*Hd - ZijN_uV = ey_px_u(hy_py*v) - ey_py*hy_px_u(v) + ey_px*hy_py_u(v) - ey_py_u(hy_px*v) + ZijN_uV = ey_px_u(v)*hy_py + ey_px*hy_py_u(v) - ey_py*hy_px_u(v) - ey_py_u(v)*hy_px elif 'zyy' in self.rxType: Zij = (-ey_px*hx_py + ey_py*hx_px)*Hd - ZijN_uV = -ey_px_u(hx_py*v) + ey_py*hx_px_u(v) - ey_px*hx_py_u(v) +ey_py_u(hx_px*v) + ZijN_uV = -ey_px_u(v)*hx_py - ey_px*hx_py_u(v) + ey_py*hx_px_u(v) + ey_py_u(v)*hx_px # Calculate the complex derivative - PDeriv_complex = ZijN_uV*Hd.toarray() - Zij * (Hd_uV*Hd.toarray()) + PDeriv_complex = Hd * (ZijN_uV - Zij * Hd_uV ) # Extract the real number for the real/imag components. Pv = np.array(getattr(PDeriv_complex, real_or_imag)) diff --git a/simpegMT/Tests/test_Problem3D_againstAnalytic.py b/simpegMT/Tests/test_Problem3D_againstAnalytic.py index cca91ef8..cbf3ed11 100644 --- a/simpegMT/Tests/test_Problem3D_againstAnalytic.py +++ b/simpegMT/Tests/test_Problem3D_againstAnalytic.py @@ -173,18 +173,22 @@ def test_DerivProjfields(inputSetup,comp='All',freq=False): survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp,freq) print 'Derivative test of data projection for eFormulation primary/secondary\n\n' # problem.mapping = simpeg.Maps.ExpMap(problem.mesh) + # Initate things for the derivs Test + src = survey.srcList[0] + rx = src.rxList[0] - # Define a src and rx - src = survey.srcList[-1] - rx = src.rxList[1] - u0 = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j + u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j + u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j + u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2))) f0 = problem.fieldsPair(survey.mesh,survey) - f0[src,'e_pxSolution'] = u0 - f0[src,'e_pySolution'] = u0 + # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2))) + f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x + f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y + def fun(u): f = problem.fieldsPair(survey.mesh,survey) - f[src,'e_pxSolution'] = u.ravel() - f[src,'e_pySolution'] = u.ravel() + f[src,'e_pxSolution'] = u[:len(u)/2] + f[src,'e_pySolution'] = u[len(u)/2::] return rx.projectFields(src,survey.mesh,f), lambda t: rx.projectFieldsDeriv(src,survey.mesh,f0,simpeg.mkvc(t,2)) return simpeg.Tests.checkDerivative(fun, u0, num=3, plotIt=False, eps=FLR)