From e3a2ec6c8dfa5e2a77c205a014ae9a53fea48eaf Mon Sep 17 00:00:00 2001 From: GudniRos Date: Wed, 24 Jun 2015 12:33:08 -0700 Subject: [PATCH] JTvec is working but not converging. --- notebooks/Derivative test MT1D.ipynb | 92 ++++++++-------------------- simpegMT/BaseMT.py | 2 +- simpegMT/SurveyMT.py | 14 +++-- 3 files changed, 35 insertions(+), 73 deletions(-) diff --git a/notebooks/Derivative test MT1D.ipynb b/notebooks/Derivative test MT1D.ipynb index a0ffef2c..9e8e36c3 100644 --- a/notebooks/Derivative test MT1D.ipynb +++ b/notebooks/Derivative test MT1D.ipynb @@ -237,11 +237,11 @@ "==================== checkDerivative ====================\n", "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", "---------------------------------------------------------\n", - " 0 1.00e-01 3.893e-05 5.627e-08 nan\n", - " 1 1.00e-02 3.898e-06 5.632e-10 2.000\n", - " 2 1.00e-03 3.898e-07 5.633e-12 2.000\n", - " 3 1.00e-04 3.898e-08 5.633e-14 2.000\n", - " 4 1.00e-05 3.898e-09 5.637e-16 2.000\n", + " 0 1.00e-01 9.779e-05 3.219e-06 nan\n", + " 1 1.00e-02 1.007e-05 3.145e-08 2.010\n", + " 2 1.00e-03 1.010e-06 3.137e-10 2.001\n", + " 3 1.00e-04 1.010e-07 3.136e-12 2.000\n", + " 4 1.00e-05 1.010e-08 3.136e-14 2.000\n", "========================= PASS! =========================\n", "Yay passed!\n", "\n" @@ -277,7 +277,7 @@ { "data": { "text/plain": [ - "array(-0.005683351098616839)" + "array([ 0.00017028])" ] }, "execution_count": 8, @@ -299,7 +299,7 @@ { "data": { "text/plain": [ - "array([[ 0.00374478]])" + "array([[ 0.0014539]])" ] }, "execution_count": 9, @@ -338,13 +338,13 @@ "---------------------------------------------------------\n", "Project at freq: 1.000e+02\n", "Project at freq: 1.000e+02\n", - " 0 1.00e-01 7.283e-06 9.326e-09 nan\n", + " 0 1.00e-01 8.466e-06 5.557e-08 nan\n", "Project at freq: 1.000e+02\n", - " 1 1.00e-02 7.290e-07 9.338e-11 1.999\n", + " 1 1.00e-02 8.416e-07 5.519e-10 2.003\n", "Project at freq: 1.000e+02\n", - " 2 1.00e-03 7.290e-08 9.339e-13 2.000\n", + " 2 1.00e-03 8.411e-08 5.515e-12 2.000\n", "========================= PASS! =========================\n", - "Testing is important.\n", + "Not just a pretty face Gudni\n", "\n" ] }, @@ -439,22 +439,19 @@ "name": "stdout", "output_type": "stream", "text": [ - "Adjoint E formulation - e\n" + "Adjoint E formulation - e\n", + "-4.65700174631e-05 3.27055690471e-05 -7.92755865102e-05 1e-08 False\n" ] }, { - "ename": "IndexError", - "evalue": "tuple index out of range", - "output_type": "error", - "traceback": [ - "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[1;31mIndexError\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[0madjointTest\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'E'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'e'\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;36madjointTest\u001b[1;34m(fdemType, comp)\u001b[0m\n\u001b[0;32m 16\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 17\u001b[0m \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[0mJvec\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mw\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mu\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\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[0mJtvec\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mv\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mu\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 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[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[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/BaseMT.pyc\u001b[0m in \u001b[0;36mJtvec\u001b[1;34m(self, m, v, f)\u001b[0m\n\u001b[0;32m 102\u001b[0m \u001b[0mPTv\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mrx\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mprojectFieldsDeriv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrc\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[0mf\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mv\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mrx\u001b[0m\u001b[1;33m]\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;31m# wrt u, need possibility wrt m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 103\u001b[0m \u001b[1;31m# Get the\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 104\u001b[1;33m \u001b[0mdA_duIT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mATinv\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mPTv\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 105\u001b[0m \u001b[0mdA_dmT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\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[0mu_src\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdA_duIT\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[0m\n\u001b[0;32m 106\u001b[0m \u001b[0mdRHS_dmT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetRHSDeriv_m\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdA_duIT\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[0m\n", - "\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Utils/SolverUtils.pyc\u001b[0m in \u001b[0;36m__mul__\u001b[1;34m(self, b)\u001b[0m\n\u001b[0;32m 35\u001b[0m \u001b[1;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Can only multiply by a numpy array.'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 36\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 37\u001b[1;33m \u001b[1;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m1\u001b[0m \u001b[1;32mor\u001b[0m \u001b[0mb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 38\u001b[0m \u001b[0mb\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mflatten\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 39\u001b[0m \u001b[1;31m# Just one RHS\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;31mIndexError\u001b[0m: tuple index out of range" - ] + "data": { + "text/plain": [ + "False" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ @@ -463,34 +460,14 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 17, "metadata": { "collapsed": false, "scrolled": true }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "> \u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Utils/SolverUtils.py\u001b[0m(37)\u001b[0;36m__mul__\u001b[1;34m()\u001b[0m\n", - "\u001b[1;32m 36 \u001b[1;33m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[0m\u001b[1;32m---> 37 \u001b[1;33m \u001b[1;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m1\u001b[0m \u001b[1;32mor\u001b[0m \u001b[0mb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[0m\u001b[1;32m 38 \u001b[1;33m \u001b[0mb\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mflatten\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[0m\n", - "ipdb> u\n", - "> \u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/BaseMT.py\u001b[0m(104)\u001b[0;36mJtvec\u001b[1;34m()\u001b[0m\n", - "\u001b[1;32m 103 \u001b[1;33m \u001b[1;31m# Get the\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[0m\u001b[1;32m--> 104 \u001b[1;33m \u001b[0mdA_duIT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mATinv\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mPTv\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[0m\u001b[1;32m 105 \u001b[1;33m \u001b[0mdA_dmT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\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[0mu_src\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdA_duIT\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[0m\n", - "\u001b[0m\n", - "ipdb> PTv\n", - "array(-6.9926641874700175e-06)\n" - ] - } - ], + "outputs": [], "source": [ - "%debug" + "# %debug" ] }, { @@ -500,10 +477,7 @@ "collapsed": false }, "outputs": [], - "source": [ - "v = np.random.randn(1,1) + 1j* np.random.randn(1,1)\n", - "rx.projectFieldsDeriv(src, problem.mesh, fields, v, adjoint=True)" - ] + "source": [] }, { "cell_type": "code", @@ -512,21 +486,7 @@ "collapsed": false }, "outputs": [], - "source": [ - "dMf_dsig = problem.mesh.getFaceInnerProductDeriv(problem.curModel.sigma)(u0) * problem.curModel.sigmaDeriv\n", - "dsig_dm = self.curModel.sigmaDeriv\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "problem.mesh.getFaceInnerProductDeriv(problem.curModel.sigma)(u0)" - ] + "source": [] }, { "cell_type": "code", diff --git a/simpegMT/BaseMT.py b/simpegMT/BaseMT.py index 370f6476..bc630a5d 100644 --- a/simpegMT/BaseMT.py +++ b/simpegMT/BaseMT.py @@ -103,7 +103,7 @@ class BaseMTProblem(BaseFDEMProblem): # Get the dA_duIT = ATinv * PTv dA_dmT = self.getADeriv_m(freq, u_src, dA_duIT, adjoint=True) - dRHS_dmT = self.getRHSDeriv_m(src, dA_duIT, adjoint=True) + dRHS_dmT = self.getRHSDeriv_m(freq, dA_duIT, adjoint=True) # Make du_dmT if dRHS_dmT is None: du_dmT = - dA_dmT diff --git a/simpegMT/SurveyMT.py b/simpegMT/SurveyMT.py index 96f6a902..3256015b 100644 --- a/simpegMT/SurveyMT.py +++ b/simpegMT/SurveyMT.py @@ -146,9 +146,10 @@ class RxMT(Survey.BaseRx): Pbx = mesh.getInterpolationMat(self.locs,'Ex') # ex = Pex*mkvc(f[src,'e_1d'],2) # bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0 - dP_de = Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pex*v) - dP_db = - 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) - PDeriv_complex = np.sum((dP_de,dP_db)) + 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 'Z3D': raise NotImplementedError('Has not be implement for full impedance tensor') @@ -160,9 +161,10 @@ class RxMT(Survey.BaseRx): # ex = Pex*mkvc(f[src,'e_1d'],2) # bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0 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 - dP_dbTv = -mkvc(f._bDeriv_u(src,db_duv,adjoint=True)*v,2) - PDeriv_complex = np.sum((dP_deTv,dP_dbTv)) + 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') elif self.projType is 'Z3D': raise NotImplementedError('must be real or imag') Pv = np.array(getattr(PDeriv_complex, real_or_imag))