diff --git a/notebooks/Derivative test MT1D.ipynb b/notebooks/Derivative test MT1D.ipynb index cde75b75..f75403d9 100644 --- a/notebooks/Derivative test MT1D.ipynb +++ b/notebooks/Derivative test MT1D.ipynb @@ -237,13 +237,13 @@ "==================== checkDerivative ====================\n", "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", "---------------------------------------------------------\n", - " 0 1.00e-01 1.076e-05 3.966e-07 nan\n", - " 1 1.00e-02 1.040e-06 3.875e-09 2.010\n", - " 2 1.00e-03 1.037e-07 3.866e-11 2.001\n", - " 3 1.00e-04 1.037e-08 3.865e-13 2.000\n", - " 4 1.00e-05 1.037e-09 3.865e-15 2.000\n", + " 0 1.00e-01 1.094e-05 3.087e-08 nan\n", + " 1 1.00e-02 1.097e-06 3.093e-10 1.999\n", + " 2 1.00e-03 1.097e-07 3.093e-12 2.000\n", + " 3 1.00e-04 1.097e-08 3.093e-14 2.000\n", + " 4 1.00e-05 1.097e-09 3.095e-16 2.000\n", "========================= PASS! =========================\n", - "The test be workin!\n", + "You are awesome.\n", "\n" ] }, @@ -277,7 +277,7 @@ { "data": { "text/plain": [ - "array([ -2.07544500e-05])" + "array([ 0.00127341])" ] }, "execution_count": 8, @@ -299,7 +299,7 @@ { "data": { "text/plain": [ - "array([[-0.00014517]])" + "array([[ 0.00173178]])" ] }, "execution_count": 9, @@ -338,15 +338,15 @@ "---------------------------------------------------------\n", "Project at freq: 1.000e+02\n", "Project at freq: 1.000e+02\n", - " 0 1.00e-01 1.142e-06 2.072e-08 nan\n", + " 0 1.00e-01 2.646e-06 2.277e-08 nan\n", "Project at freq: 1.000e+02\n", - " 1 1.00e-02 1.161e-07 2.039e-10 2.007\n", + " 1 1.00e-02 2.629e-07 2.246e-10 2.006\n", "Project at freq: 1.000e+02\n", - " 2 1.00e-03 1.162e-08 2.036e-12 2.001\n", + " 2 1.00e-03 2.627e-08 2.243e-12 2.001\n", "Project at freq: 1.000e+02\n", - " 3 1.00e-04 1.163e-09 2.036e-14 2.000\n", + " 3 1.00e-04 2.627e-09 2.242e-14 2.000\n", "========================= PASS! =========================\n", - "Testing is important.\n", + "Yay passed!\n", "\n" ] }, @@ -411,7 +411,7 @@ "output_type": "stream", "text": [ "Adjoint e formulation - projectFieldsDeriv\n", - "-9.0191451349e-06 -9.0191451349e-06 3.38813178902e-21 1e-09 True\n" + "-0.000991536643367 -0.000991536643367 2.16840434497e-19 1e-07 True\n" ] }, { @@ -440,10 +440,10 @@ " u = problem.fields(m)\n", " v = np.random.randn(1)#+np.random.randn(1)*1j\n", " # print prb.PropMap.PropModel.nP\n", - " w = np.random.randn(m1d.nN)#+np.random.randn(m1d.nN)*1j\n", + " w = np.random.randn(m1d.nN)+np.random.randn(m1d.nN)*1j\n", "\n", " vJw = v.dot(rx.projectFieldsDeriv(src,m1d,f0,w))\n", - " wJtv = w.dot(rx.projectFieldsDeriv(src,m1d,f0,v,adjoint=True))\n", + " wJtv = w.dot(rx.projectFieldsDeriv(src,m1d,f0,v,adjoint=True)).real\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", @@ -456,13 +456,32 @@ "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": 15, + "metadata": { + "collapsed": false + }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Adjoint test e formulation - getADeriv_m\n", - "(508851.474801-63591.6147417j) (508851.474801-63591.6147417j) (-1.7462298274e-10-5.82076609135e-11j) 10.0 True\n" + "(340193.379835-398622.996348j) (340193.379835-398622.996348j) (-2.32830643654e-10-3.49245965481e-10j) 10.0 True\n" ] }, { @@ -471,7 +490,7 @@ "True" ] }, - "execution_count": 14, + "execution_count": 15, "metadata": {}, "output_type": "execute_result" } @@ -512,7 +531,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 16, "metadata": { "collapsed": false }, @@ -522,7 +541,7 @@ "output_type": "stream", "text": [ "Adjoint test e formulation - getRHSDeriv_m\n", - "(-29427.4754714-17701.4445852j) (-29427.4754714-17701.4445852j) (2.91038304567e-11+0j) 1.0 True\n" + "(-12351.1349263+433.4655562j) (-12351.1349263+433.4655562j) (-9.09494701773e-12-1.22781784739e-11j) 1.0 True\n" ] }, { @@ -531,7 +550,7 @@ "True" ] }, - "execution_count": 15, + "execution_count": 16, "metadata": {}, "output_type": "execute_result" } @@ -563,7 +582,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 17, "metadata": { "collapsed": false }, @@ -584,7 +603,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 18, "metadata": { "collapsed": false }, @@ -593,8 +612,8 @@ "TOL = 1e-4\n", "FLR = 1e-20\n", "\n", - "def JvecAdjointTest(fdemType, comp):\n", - " print 'Adjoint %s formulation - %s' % (fdemType, comp)\n", + "def JvecAdjointTest():\n", + " print 'Adjoint e formulation - Jvec' \n", "\n", " m = np.log(np.ones(problem.mesh.nC)*0.01)\n", " if True:\n", @@ -615,7 +634,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 19, "metadata": { "collapsed": false }, @@ -624,28 +643,28 @@ "name": "stdout", "output_type": "stream", "text": [ - "Adjoint E formulation - e\n", - "-1.7867035849e-05 1.42750293876e-05 -3.21420652366e-05 1e-08 False\n" + "Adjoint e formulation - Jvec\n", + "-3.61480355369e-05 -3.61480355369e-05 -2.71050543121e-20 1e-08 True\n" ] }, { "data": { "text/plain": [ - "False" + "True" ] }, - "execution_count": 18, + "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "JvecAdjointTest('E','e')" + "JvecAdjointTest()" ] }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 20, "metadata": { "collapsed": false, "scrolled": true @@ -665,7 +684,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 21, "metadata": { "collapsed": false, "scrolled": true @@ -684,7 +703,7 @@ "array([ 9.80523303e-06, -1.98372645e-03])" ] }, - "execution_count": 20, + "execution_count": 21, "metadata": {}, "output_type": "execute_result" } @@ -695,32 +714,34 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": 21, + "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { - "data": { - "text/plain": [ - "<180x180 sparse matrix of type ''\n", - "\twith 180 stored elements in Compressed Sparse Row format>" - ] - }, - "execution_count": 21, - "metadata": {}, - "output_type": "execute_result" + "ename": "NameError", + "evalue": "name 'r' is not defined", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mNameError\u001b[0m Traceback (most recent call last)", + "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mr\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[1;31mNameError\u001b[0m: name 'r' is not defined" + ] } ], + "source": [ + "r" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], "source": [ "problem.mesh.getEdgeInnerProductDeriv(problem.curModel.sigma)(u0[1::])" ] @@ -740,18 +761,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 f1092f03..5ec79aab 100644 --- a/simpegMT/BaseMT.py +++ b/simpegMT/BaseMT.py @@ -61,16 +61,6 @@ class BaseMTProblem(BaseFDEMProblem): du_dm = dA_duI * ( - dA_dm + dRHS_dm ) # Calculate the projection derivatives for rx in src.rxList: - # Get the stacked derivative - # df_duFun = getattr(f, '_fDeriv_u', None) - # df_dmFun = getattr(f, '_fDeriv_m', None) - # df_dm = df_dmFun(src,v,adjoint=False) - # if df_dm is None: - # fDeriv_m = df_duFun(src, du_dm, adjoint=False) - # else: - # fDeriv_m = df_duFun(src, du_dm, adjoint=False) + df_dm - # Not needed for now. Since PDeriv does this currently. - # Get the projection derivative PDeriv = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m Jv[src, rx] = PDeriv(du_dm) diff --git a/simpegMT/SurveyMT.py b/simpegMT/SurveyMT.py index c440a812..4175efcd 100644 --- a/simpegMT/SurveyMT.py +++ b/simpegMT/SurveyMT.py @@ -149,12 +149,14 @@ class RxMT(Survey.BaseRx): dP_de = mkvc(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pex*v),2) dP_db = mkvc(- Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Pbx*f._bDeriv_u(src,v)/mu_0),2) PDeriv_complex = np.sum(np.hstack((dP_de,dP_db)),1) - # raise Exception('Debug error') - # elif self.projType is 'Z2D + elif self.projType is 'Z2D': + raise NotImplementedError('Has not be implement for 2D impedance tensor') elif self.projType is 'Z3D': - raise NotImplementedError('Has not be implement for full impedance tensor') + raise NotImplementedError('Has not be implement for full 3D impedance tensor') + # Extract the real number for the real/imag components. Pv = np.array(getattr(PDeriv_complex, real_or_imag)) elif adjoint: + # Note: The v vector is real and the return should be complex if self.projType is 'Z1D': Pex = mesh.getInterpolationMat(self.locs,'Fx') Pbx = mesh.getInterpolationMat(self.locs,'Ex') @@ -163,11 +165,17 @@ class RxMT(Survey.BaseRx): dP_deTv = mkvc(Pex.T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*v,2) db_duv = -Pbx.T/mu_0*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2)).T*v dP_dbTv = mkvc(f._bDeriv_u(src,db_duv,adjoint=True),2) - PDeriv_complex = np.sum(np.hstack((dP_deTv,dP_dbTv)),1) - # raise Exception('Debug error') + PDeriv_real = np.sum(np.hstack((dP_deTv,dP_dbTv)),1) + elif self.projType is 'Z2D': + raise NotImplementedError('Has not be implement for 2D impedance tensor') elif self.projType is 'Z3D': - raise NotImplementedError('must be real or imag') - Pv = np.array(getattr(PDeriv_complex, real_or_imag)) + raise NotImplementedError('Has not be implement for full 3D impedance tensor') + # Extract the data + if real_or_imag == 'imag': + Pv = 1j*PDeriv_real + elif real_or_imag == 'real': + Pv = PDeriv_real.astype(complex) + return Pv