diff --git a/notebooks/Derivative test MT1D.ipynb b/notebooks/Derivative test MT1D.ipynb index e0027572..29e5ce9c 100644 --- a/notebooks/Derivative test MT1D.ipynb +++ b/notebooks/Derivative test MT1D.ipynb @@ -13,7 +13,7 @@ "cell_type": "code", "execution_count": 1, "metadata": { - "collapsed": true + "collapsed": false }, "outputs": [], "source": [ @@ -29,6 +29,28 @@ "metadata": { "collapsed": false }, + "outputs": [ + { + "data": { + "text/plain": [ + "simpegMT.FieldsMT.FieldsMT_1D" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "simpegmt.FieldsMT.FieldsMT_1D" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, "outputs": [ { "name": "stdout", @@ -85,6 +107,15 @@ "data = survey.projectFields(fields)\n" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, { "cell_type": "markdown", "metadata": {}, @@ -113,7 +144,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 4, "metadata": { "collapsed": true }, @@ -128,7 +159,7 @@ "source": [ "As matrices the formulas above can be written as\n", "\\begin{align}\n", - "\\left[ \\frac{\\partial P(f(m))}{\\partial u} v \\right] = diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right] [P_e v] - diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right]^T diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right] \\left[ P_b \\frac{d f_b}{du}(v) \\frac{1}{mu_0} \\right]\n", + "\\left[ \\frac{\\partial P(f(m))}{\\partial u} v \\right] = diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right] [P_e v] - diag[P_e u] diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right]^T diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right] \\left[ P_b \\frac{d f_b}{du}(v) \\frac{1}{mu_0} \\right]\n", "\\end{align}\n", "\n" ] @@ -139,13 +170,13 @@ "source": [ "The adjoint problem is done simliarly\n", "\\begin{align}\n", - "\\left[ \\frac{\\partial P(f(m))}{\\partial u} v \\right]^T = [P_e v]^T diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right]^T - \\left[ P_b \\frac{d f_b}{du}(v) \\frac{1}{mu_0} \\right]^T diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right] diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right]^T \n", + "\\left[ \\frac{\\partial P(f(m))}{\\partial u} v \\right]^T = [P_e v]^T diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right]^T - \\left[ P_b \\frac{d f_b}{du}(v) \\frac{1}{mu_0} \\right]^T diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right] diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right]^T diag \\left[ P_e u \\right]^T\n", "\\end{align}\n" ] }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 5, "metadata": { "collapsed": true }, @@ -169,7 +200,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 6, "metadata": { "collapsed": false }, @@ -178,15 +209,16 @@ "# Initate things for the derivs Test\n", "src = survey.srcList[0]\n", "rx = src.rxList[0]\n", + "v = np.random.randn(m1d.nN)\n", "u0 = np.random.randn(m1d.nN)+np.random.randn(m1d.nN)*1j\n", "f0 = problem.fieldsPair(m1d,survey)\n", - "f0[src,'e_1d'] = u0\n", - "f0[src,'b_1d'] = -1/(1j*simpegem.Utils.EMUtils.omega(src.freq))*m1d.nodalGrad*u0" + "f0[src,'e_1dSolution'] = u0\n", + "# f0[src,'b_1d'] = -1/(1j*simpegem.Utils.EMUtils.omega(src.freq))*m1d.nodalGrad*u0" ] }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 7, "metadata": { "collapsed": false, "scrolled": true @@ -195,100 +227,100 @@ { "data": { "text/plain": [ - "array([ 0.36132584-1.64741907j, -1.09543652-0.5983257j ,\n", - " -0.37399249-0.17183374j, -0.29875616+0.13386035j,\n", - " -0.13497219-1.67592883j, -0.43805630+0.79351159j,\n", - " 1.10360642-1.88408115j, -1.55655836+0.09909074j,\n", - " 1.09262762-0.26644991j, 1.00456242+1.36030358j,\n", - " 0.10649574-0.2738965j , -0.93711934-0.82724463j,\n", - " 1.86172288-1.39322063j, 0.17645120-2.38628835j,\n", - " 0.68232151-0.96482639j, -1.54748216+0.02907389j,\n", - " -0.33109688+0.19681238j, -0.70235831-1.00115598j,\n", - " -1.13598355+0.21054259j, -0.09009229-0.07819883j,\n", - " -0.01097127+0.68165894j, 0.52661676+0.27378467j,\n", - " -0.70418270-1.10247881j, -0.51995781-0.20521301j,\n", - " -0.12731995-0.33125501j, -1.59215935-2.21034752j,\n", - " -0.18124185+0.32745866j, 0.67431311+1.33705062j,\n", - " 0.92167361-0.38476494j, 0.81511489-0.31198285j,\n", - " -0.34442259+0.957995j , -0.18007675-0.11694359j,\n", - " -0.72309201+0.17539006j, 0.61922884+0.50879561j,\n", - " 0.47818118+0.81750598j, -0.37937463-0.33322861j,\n", - " -0.06936763-0.64247149j, -1.19310390+0.49472846j,\n", - " 0.95134389+0.1500819j , -0.44352264+2.10893973j,\n", - " 2.09530714-0.49663377j, -1.13328736+0.46104269j,\n", - " 0.00899419+0.98561084j, 1.79288071-0.98577358j,\n", - " -1.21189456-1.69393797j, -0.79574262+0.49529776j,\n", - " -0.58983776+0.38069564j, -2.30858124+2.30767413j,\n", - " -1.58107977-0.88577507j, -0.94772249+0.1889765j ,\n", - " -0.38215882-0.31335661j, 0.38162110-1.18276337j,\n", - " -0.15421254-0.10918159j, -0.58284090-1.41758673j,\n", - " 0.85594502-0.22857923j, 1.57628220+0.9592039j ,\n", - " 0.46253320+0.18527238j, -1.80418304-1.85163236j,\n", - " -0.40546603+1.23461379j, -0.11477894+0.22189668j,\n", - " -1.16436004-0.11841764j, 0.30911061-1.32708718j,\n", - " 0.32877677+0.80977796j, 1.67348108-0.09360135j,\n", - " 0.63761689+1.22794045j, -0.34255102-0.11160981j,\n", - " 0.59182860-0.13894497j, -2.43635154-1.21416755j,\n", - " -0.17375475-0.77779331j, 1.11791248+0.43335221j,\n", - " -0.29401484-0.45669262j, 1.39676700-1.17486543j,\n", - " 1.71647097+1.6277104j , 1.95008743+1.06926207j,\n", - " -0.67552544+0.48327768j, -0.11000128-0.82649714j,\n", - " 0.96933809-0.37167971j, 1.47494336+0.24152755j,\n", - " 0.75463706+0.97596453j, 0.27418483-1.46035469j,\n", - " -2.17773096-0.16074446j, -0.27717438-1.56310869j,\n", - " -0.28038940+0.10826212j, -0.55314635+0.54203288j,\n", - " -1.09967872-0.0837615j , -0.16844443-0.85471157j,\n", - " -0.33760983-0.07555467j, -1.10668556-0.48712673j,\n", - " 1.03395454+0.28116204j, -0.48505901+0.42626281j,\n", - " 1.24278733+0.14643982j, -0.67434494-0.44708822j,\n", - " 0.88512301+0.58888926j, 0.29114641-1.21153896j,\n", - " 2.09025380+1.32932413j, 1.10848862+0.02614621j,\n", - " -1.02811260-0.8517983j , 0.26992478+0.43162846j,\n", - " 0.41252221+0.28673221j, -0.02741592-0.91702344j,\n", - " 0.17601394-1.0670725j , -1.02325192+0.01888599j,\n", - " -1.80534626-1.4886524j , 1.19712663+1.1909822j ,\n", - " 0.68362676+1.89111425j, 0.06090091+0.13260624j,\n", - " -0.75587177+1.08278499j, 0.88939421+0.86220205j,\n", - " 0.44347443-0.80363074j, 0.56780383+0.1887053j ,\n", - " -1.66790268+1.53615809j, -0.01736428-0.06506712j,\n", - " -2.25744003+1.06202019j, 0.39400372+0.3498991j ,\n", - " -0.55776393-1.08756479j, -1.11171590+0.33655697j,\n", - " -1.09279545+1.66655595j, 0.58095680+0.53565886j,\n", - " -0.93002846+0.34343169j, 0.72921615-1.8539533j ,\n", - " -0.19275461+0.41839652j, -0.13118494-0.5427038j ,\n", - " -0.42579986-0.11640398j, -2.29461643+1.02111832j,\n", - " 0.91308216-1.38100176j, -3.32212238+0.68156176j,\n", - " 0.56641453-0.95066951j, -0.16976454+0.49554783j,\n", - " 1.50014425-0.58622796j, -0.60263781-2.15510148j,\n", - " 1.21695459+1.78261241j, -0.63597912-0.8384838j ,\n", - " -3.64982757-0.55347908j, -0.97532271+0.55805176j,\n", - " -0.00551552-0.1537918j , -0.95478088+0.73917938j,\n", - " -1.48127806+0.68079515j, 1.67421015+0.25426024j,\n", - " 1.33335078-0.12101624j, 0.29645596+0.04605261j,\n", - " -0.47334868-1.05512171j, -0.26655968-1.55359388j,\n", - " 1.29559009-0.7934454j , 0.28283852-0.18507402j,\n", - " 0.23508679+2.28747714j, -0.86893735-0.00295461j,\n", - " 0.18639473+1.1307612j , -1.90052723+0.21204624j,\n", - " -0.11999826-0.27195367j, 0.29379502+0.21711147j,\n", - " 1.31071507+0.4451999j , 0.24524378+1.64849073j,\n", - " 1.68293321+0.73050439j, 0.36039213-0.90352994j,\n", - " 0.72257192+0.87566287j, 2.11195235+0.74321181j,\n", - " 0.18958287-1.06053094j, 0.46269877+1.0830655j ,\n", - " -0.32403398-1.12584818j, 0.72269683-0.36983534j,\n", - " -0.75930979+0.36833169j, -0.08659410+0.45782596j,\n", - " 0.68410073+0.64559686j, 1.08174726+0.24836564j,\n", - " 1.16991513-1.29388377j, 0.48984416+0.3837798j ,\n", - " 0.09385395-0.46595543j, 0.39693801+0.53049104j,\n", - " 0.78547311-0.61634797j, -0.30417945+1.17182948j,\n", - " -0.16075966-0.09621673j, -0.25778022-1.53597405j,\n", - " 0.56695410+2.15202438j, -0.48969409+0.11914719j,\n", - " -0.59882416-0.57579404j, -0.15237306+0.77126722j,\n", - " -1.65801751-0.60162042j, 0.30512004-2.08686648j,\n", - " 1.47493551-0.82515753j, -0.53121616-0.2578771j ,\n", - " -0.19552096-0.33023782j])" + "array([ -5.21576299e-01-0.28596337j, 1.33775588e+00-1.08159078j,\n", + " 1.02311300e+00-0.45033307j, -1.14140017e+00-0.53814712j,\n", + " -1.51095317e+00-0.0225078j , 2.29842631e-01-0.33435588j,\n", + " -6.13742564e-01+0.14777572j, 1.74709600e+00+0.07973247j,\n", + " -6.07949930e-01+0.53467919j, 7.44109732e-01-1.06383585j,\n", + " -6.22557065e-01-0.6604577j , 5.54530896e-01+0.49763768j,\n", + " -1.42655640e+00+1.03849457j, -8.79197413e-01-0.05384536j,\n", + " -5.91554616e-01-0.29187425j, -2.09918231e+00+0.65127732j,\n", + " -2.79577518e-01+0.08460458j, -8.49733033e-01-1.15162916j,\n", + " -5.61159802e-01+0.40323403j, -1.32856823e-01+0.83145628j,\n", + " 1.38730097e-01-0.64695336j, 4.74841673e-01+2.16259344j,\n", + " 5.20314820e-02+1.39337638j, 2.40669484e+00+1.44814886j,\n", + " 1.02271608e+00-0.64912151j, 6.52379785e-01+0.59394264j,\n", + " 2.50946314e-01-1.71467789j, -3.35061639e-01+1.80721483j,\n", + " 1.87006132e+00-0.18721108j, 1.85379814e+00+0.06414066j,\n", + " 7.66449796e-01+0.28757766j, 3.77939332e-01+0.5861157j ,\n", + " -2.01974930e+00+2.36530034j, -5.70969153e-01+0.16116942j,\n", + " 3.37395529e-01+0.76006325j, -1.12084482e-01-2.02591559j,\n", + " 1.36070180e-01-0.88380389j, 1.99378191e+00-0.50681261j,\n", + " 1.50662216e-01+1.76352214j, -4.33494912e-01-0.80584325j,\n", + " -5.30448843e-01+0.99832665j, -1.52568430e+00+0.46127907j,\n", + " -4.78999736e-01-0.52250698j, 5.07078441e-01+0.02687913j,\n", + " 1.52990246e+00+0.09587891j, 9.10105334e-01-0.22098742j,\n", + " 6.63562676e-01-0.45367061j, -1.08906319e+00-1.81627004j,\n", + " 3.41582990e-02+1.71027836j, -2.22532053e-01+1.06004075j,\n", + " 1.74917677e+00-1.41735414j, 1.06979706e+00-1.89094749j,\n", + " -3.35092655e-01+0.39943895j, 1.05206441e+00+0.80828213j,\n", + " -6.87168211e-01+1.72353937j, 1.40047301e+00-1.3727236j ,\n", + " -1.03495040e-01+0.26680256j, -1.12873662e+00-0.56805852j,\n", + " -8.57584940e-01+1.85418143j, 1.03116279e+00+0.06711151j,\n", + " 5.46075841e-01+1.07590359j, 2.67048944e-01+0.75324202j,\n", + " 9.60209324e-01-0.09519331j, -1.36504332e-01-1.35387957j,\n", + " -3.38824018e-01+0.70215028j, -4.82487190e-02-0.24027277j,\n", + " 1.55731408e+00-2.38970922j, -3.35634102e-03-1.11778716j,\n", + " 1.05996566e+00-1.54847977j, 1.37371396e+00+0.21763618j,\n", + " -2.34257640e+00+1.16351971j, -5.41859146e-01+1.0902067j ,\n", + " -1.79555485e+00+0.06162846j, 7.46946726e-01-0.87760098j,\n", + " -9.37453865e-02+0.26210336j, -6.00849545e-01+0.49005865j,\n", + " 2.47173181e-04-0.41345406j, 1.00840211e+00+1.23027562j,\n", + " 8.39952932e-01+0.43605496j, 8.65132399e-01+1.07984908j,\n", + " 2.04841885e+00-0.0851883j , 1.44345869e-01+0.4539564j ,\n", + " 9.19226689e-01-0.37858608j, -1.51448024e-01-1.2974253j ,\n", + " 5.32408479e-02-1.60301977j, -1.84229724e+00-0.15278222j,\n", + " -2.13350910e+00+1.87102448j, 9.03468928e-01+0.23498138j,\n", + " -9.26803345e-01+0.51213564j, 7.94552525e-01+0.54895852j,\n", + " -9.77798824e-01-0.49901395j, -6.05296444e-01-0.03209891j,\n", + " 1.24765034e+00-0.56443723j, -1.90295464e+00-0.21190697j,\n", + " 1.44761031e+00+0.14009266j, -9.23615842e-01+1.93414685j,\n", + " -2.89396133e-01+0.47920215j, 1.22494934e-01+1.33535626j,\n", + " -6.00157469e-01+0.61558122j, -1.40975425e-01+1.61693392j,\n", + " 7.08238379e-01+0.536495j , 5.90147090e-01+0.60041602j,\n", + " -7.09067864e-01+0.08489417j, -1.24030970e+00-0.73986538j,\n", + " -1.03488057e+00-0.09736817j, 2.89640823e-01-0.41632583j,\n", + " -4.94975749e-02-0.61519627j, 5.41641220e-01+0.77100395j,\n", + " 7.22415122e-01+1.59030908j, -5.13304114e-01-0.09690537j,\n", + " 5.41180474e-01+2.6467645j , -8.05664508e-01+1.44276035j,\n", + " -4.06768654e-01+1.58421981j, 1.33368455e+00+0.63576588j,\n", + " 4.29526056e-02+0.16497284j, 5.92169353e-01-0.76120082j,\n", + " 1.84734942e+00-0.29548982j, 3.25330682e-01+0.84903188j,\n", + " 1.08709120e-01-1.73236832j, -5.01403552e-01-0.083525j ,\n", + " 3.43892536e-01+0.3488206j , -2.42917431e+00-0.25699335j,\n", + " 5.53836789e-01-0.45887772j, -6.75365309e-01+0.0848306j ,\n", + " 4.78327738e-01-0.77164497j, -4.05543707e-01-0.04087569j,\n", + " 6.48445605e-01+1.00477682j, 4.50922906e-02-0.09451425j,\n", + " -2.97801954e-01-0.40438134j, -2.47619828e+00+0.25935539j,\n", + " 1.95520198e+00+0.93190891j, 7.52207479e-01+0.60075967j,\n", + " -8.95365798e-01+1.22241712j, -4.65516271e-01-0.10833261j,\n", + " -1.55272666e-01-1.20074828j, 5.78006258e-01+0.06705178j,\n", + " 9.35873450e-01-1.244475j , 3.08977268e-03-1.57055076j,\n", + " -1.69903140e+00-1.00183406j, -5.73059946e-01-1.96431137j,\n", + " 3.80255478e-01-0.58071467j, -3.27026898e-01-0.44312791j,\n", + " 2.66653209e-01+0.20442809j, -7.26239447e-02-0.73470193j,\n", + " -1.27190710e+00+0.6681173j , -1.66240293e+00-0.68067832j,\n", + " -1.30894453e-01+0.4142016j , 1.17940864e-01-0.29430427j,\n", + " 5.40779365e-01-0.1570116j , 4.59411844e-01-1.33449408j,\n", + " -2.86538969e-01-0.46958326j, -2.91593793e-02-1.25227084j,\n", + " 2.45440216e-02+0.07172228j, -8.17237335e-02-0.71911812j,\n", + " 1.74142448e+00-0.26649839j, -4.31150193e-01-0.79560216j,\n", + " -7.04116352e-01-1.53969088j, -4.88139638e-01-0.59075677j,\n", + " -1.61059835e-01-0.62180822j, -1.75918449e-01-0.31350392j,\n", + " -2.28918438e-01-1.25884104j, -3.94756005e-01-0.08464047j,\n", + " 3.23000023e-01+0.87737735j, 1.05597723e-02-0.74774027j,\n", + " -2.25509504e+00+2.81870339j, -9.95775953e-01-1.21962899j,\n", + " -8.59196834e-01+1.4407182j , -6.38861601e-01+0.67564119j,\n", + " 1.99028382e-01+2.32353483j, -1.15758435e+00-0.70250997j,\n", + " -9.18642125e-01-0.72232602j, 1.00415977e+00+0.40547797j,\n", + " -2.61807248e-01-1.27673902j, -1.39773535e+00+1.04765733j,\n", + " 8.93091651e-01-0.67650682j, -1.26058855e+00+1.24753141j,\n", + " -4.41649966e-01-0.64575466j, -1.35393159e+00+0.42326246j,\n", + " 1.35067221e-01+1.44186272j, 3.95784662e-01-0.99344302j,\n", + " 3.22900419e-01+1.13931166j])" ] }, - "execution_count": 6, + "execution_count": 7, "metadata": {}, "output_type": "execute_result" } @@ -306,7 +338,63 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 8, + "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 5.703e-05 1.929e-06 nan\n", + " 1 1.00e-02 5.877e-06 1.954e-08 1.994\n", + " 2 1.00e-03 5.894e-07 1.956e-10 1.999\n", + " 3 1.00e-04 5.896e-08 1.957e-12 2.000\n", + "========================= PASS! =========================\n", + "Not just a pretty face Gudni\n", + "\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Run a test\n", + "def fun(u):\n", + " f = problem.fieldsPair(m1d,survey)\n", + " f[src,'e_1dSolution'] = u\n", + "# f[src,'b_1d'] = -(m1d.nodalGrad*u)/(1j*simpegem.Utils.EMUtils.omega(src.freq))\n", + " return rx.projectFields(src,m1d,f), lambda t: rx.projectFieldsDeriv(src,m1d,f0,t)\n", + "simpeg.Tests.checkDerivative(fun,u0,num=4,plotIt=False)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Test the Jvec derivative." + ] + }, + { + "cell_type": "code", + "execution_count": 10, "metadata": { "collapsed": false }, @@ -318,39 +406,60 @@ "==================== checkDerivative ====================\n", "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", "---------------------------------------------------------\n", - " 0 1.00e-01 1.969e-05 8.331e-07 nan\n", - " 1 1.00e-02 2.045e-06 7.979e-09 2.019\n", - " 2 1.00e-03 2.052e-07 7.945e-11 2.002\n", - " 3 1.00e-04 2.052e-08 7.942e-13 2.000\n", - "========================= PASS! =========================\n", - "That was easy!\n", - "\n" + "Project at freq: 1.000e+02\n", + "Project at freq: 1.000e+02\n" ] }, { - "data": { - "text/plain": [ - "True" - ] - }, - "execution_count": 7, - "metadata": {}, - "output_type": "execute_result" + "ename": "ValueError", + "evalue": "dimension mismatch", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mValueError\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 13\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 14\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0msurvey\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdpred\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;32mlambda\u001b[0m \u001b[0mx\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[0mm0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 15\u001b[1;33m \u001b[0msimpeg\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mTests\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcheckDerivative\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mu0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnum\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mplotIt\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Tests/TestUtils.pyc\u001b[0m in \u001b[0;36mcheckDerivative\u001b[1;34m(fctn, x0, num, plotIt, dx, expectedOrder, tolerance, eps, ax)\u001b[0m\n\u001b[0;32m 256\u001b[0m \u001b[1;31m# 1st order Taylor\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 257\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0minspect\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0misfunction\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mJ0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 258\u001b[1;33m \u001b[0mE1\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0ml2norm\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mft\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mf0\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mh\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mJ0\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdx\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 259\u001b[0m \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 260\u001b[0m \u001b[1;31m# We assume it is a numpy.ndarray\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m(x)\u001b[0m\n\u001b[0;32m 12\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0mfun\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 13\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 14\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0msurvey\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdpred\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;32mlambda\u001b[0m \u001b[0mx\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[0mm0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 15\u001b[0m \u001b[0msimpeg\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mTests\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcheckDerivative\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mu0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnum\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mplotIt\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\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/BaseMT.pyc\u001b[0m in \u001b[0;36mJvec\u001b[1;34m(self, m, v, f)\u001b[0m\n\u001b[0;32m 64\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 65\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 66\u001b[1;33m \u001b[0mJv\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[0mP\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdu_dm\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 67\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 68\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mJv\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/BaseMT.pyc\u001b[0m in \u001b[0;36m\u001b[1;34m(v)\u001b[0m\n\u001b[0;32m 61\u001b[0m \u001b[0mdu_dm\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[0mdf_dm\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 62\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 63\u001b[1;33m \u001b[0mP\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mlambda\u001b[0m \u001b[0mv\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[1;31m# wrt u, also have wrt m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 64\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 65\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/SurveyMT.py\u001b[0m in \u001b[0;36mprojectFieldsDeriv\u001b[1;34m(self, src, mesh, f, v, adjoint)\u001b[0m\n\u001b[0;32m 141\u001b[0m \u001b[1;31m# ex = Pex*mkvc(f[src,'e_1d'],2)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 142\u001b[0m \u001b[1;31m# bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 143\u001b[1;33m \u001b[0mderiv_complex\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1.\u001b[0m\u001b[1;33m/\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPbx\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'b_1d'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mmu_0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPex\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPex\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'e_1d'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1.\u001b[0m\u001b[1;33m/\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPbx\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'b_1d'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mmu_0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1.\u001b[0m\u001b[1;33m/\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPbx\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'b_1d'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mmu_0\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[0mPbx\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_bDeriv_u\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m)\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[0m\u001b[0;32m 144\u001b[0m \u001b[1;31m# elif self.projType is 'Z2D\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 145\u001b[0m \u001b[1;32melif\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mprojType\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;34m'Z3D'\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[1;32m/home/gudni/anaconda/lib/python2.7/site-packages/scipy/sparse/base.pyc\u001b[0m in \u001b[0;36m__mul__\u001b[1;34m(self, other)\u001b[0m\n\u001b[0;32m 325\u001b[0m \u001b[1;31m# dense row or column vector\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 326\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mother\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m \u001b[1;33m!=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mN\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mand\u001b[0m \u001b[0mother\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m \u001b[1;33m!=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mN\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 327\u001b[1;33m \u001b[1;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'dimension mismatch'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 328\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 329\u001b[0m \u001b[0mresult\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_mul_vector\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mravel\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mother\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[1;31mValueError\u001b[0m: dimension mismatch" + ] } ], "source": [ - "# Run a test\n", - "def fun(u):\n", - " f = problem.fieldsPair(m1d,survey)\n", - " f[src,'e_1d'] = u\n", - " f[src,'b_1d'] = -(m1d.nodalGrad*u)/(1j*simpegem.Utils.EMUtils.omega(src.freq))\n", - " return rx.projectFields(src,m1d,f), lambda t: rx.projectFieldsDeriv(src,m1d,f0,t)\n", - "simpeg.Tests.checkDerivative(fun,u0,num=4,plotIt=False)" + "# print '%s formulation - %s' % (fdemType, comp)\n", + "CONDUCTIVITY = 0.01\n", + "m0 = np.log(np.ones(problem.mesh.nC)*CONDUCTIVITY)\n", + "# mu = np.log(np.ones(problem.mesh.nC)*MU)\n", + "\n", + "if True:\n", + " m0 = m0 + np.random.randn(problem.mesh.nC)*CONDUCTIVITY*1e-1 \n", + "# mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1\n", + "\n", + "# prb.mu = mu\n", + "# survey = prb.survey\n", + "def fun(x):\n", + " \n", + " return survey.dpred(m0), lambda x: problem.Jvec(m0, x)\n", + "simpeg.Tests.checkDerivative(fun, u0, num=3, plotIt=False)" ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "problem.getADeriv_m(freq,fields[src,'e_1dSolution'],v)" + ] + }, + { + "cell_type": "code", + "execution_count": null, "metadata": { "collapsed": false }, @@ -359,64 +468,80 @@ "name": "stdout", "output_type": "stream", "text": [ - "[ -8.67361738e-19]\n", - "181\n", - "181\n" + "> \u001b[1;32m/home/gudni/anaconda/lib/python2.7/site-packages/scipy/sparse/base.py\u001b[0m(327)\u001b[0;36m__mul__\u001b[1;34m()\u001b[0m\n", + "\u001b[1;32m 326 \u001b[1;33m \u001b[1;32mif\u001b[0m \u001b[0mother\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m \u001b[1;33m!=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mN\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mand\u001b[0m \u001b[0mother\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m \u001b[1;33m!=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mN\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\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--> 327 \u001b[1;33m \u001b[1;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'dimension mismatch'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[0m\u001b[1;32m 328 \u001b[1;33m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[0m\n", + "ipdb> u\n", + "> \u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/SurveyMT.py\u001b[0m(143)\u001b[0;36mprojectFieldsDeriv\u001b[1;34m()\u001b[0m\n", + "\u001b[1;32m 142 \u001b[1;33m \u001b[1;31m# bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[0m\u001b[1;32m--> 143 \u001b[1;33m \u001b[0mderiv_complex\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1.\u001b[0m\u001b[1;33m/\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPbx\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'b_1d'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mmu_0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPex\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPex\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'e_1d'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1.\u001b[0m\u001b[1;33m/\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPbx\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'b_1d'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mmu_0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1.\u001b[0m\u001b[1;33m/\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPbx\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'b_1d'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mmu_0\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[0mPbx\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_bDeriv_u\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m)\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[0m\u001b[1;32m 144 \u001b[1;33m \u001b[1;31m# elif self.projType is 'Z2D\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[0m\n", + "ipdb> Pex*v\n", + "*** ValueError: dimension mismatch\n", + "ipdb> v.shape\n", + "(180,)\n" ] } ], "source": [ - "print rx.projectFieldsDeriv(src,m1d,f0,u0)\n", - "print m1d.nF\n", - "print m1d.nN" + "%debug" ] }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ - "# fields._b_1dDeriv_u(src,u0)" + "problem.getA" ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "metadata": { "collapsed": false }, - "outputs": [ - { - "ename": "TypeError", - "evalue": "object of type 'FieldsMT' has no len()", - "output_type": "error", - "traceback": [ - "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[1;31mTypeError\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[0msurvey\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdpred\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", - "\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Utils/CounterUtils.pyc\u001b[0m in \u001b[0;36mwrapper\u001b[1;34m(self, *args, **kwargs)\u001b[0m\n\u001b[0;32m 81\u001b[0m \u001b[0mcounter\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'counter'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 82\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mtype\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcounter\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mCounter\u001b[0m\u001b[1;33m:\u001b[0m \u001b[0mcounter\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcount\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__class__\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__name__\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;34m'.'\u001b[0m\u001b[1;33m+\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__name__\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 83\u001b[1;33m \u001b[0mout\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 84\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 85\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mwrapper\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Utils/codeutils.pyc\u001b[0m in \u001b[0;36mrequiresVarWrapper\u001b[1;34m(self, *args, **kwargs)\u001b[0m\n\u001b[0;32m 224\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvar\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[0;32m 225\u001b[0m \u001b[1;32mraise\u001b[0m \u001b[0mException\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mextra\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 226\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 227\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 228\u001b[0m \u001b[0mdoc\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mrequiresVarWrapper\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__doc__\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;36mdpred\u001b[1;34m(self, m, u)\u001b[0m\n\u001b[0;32m 306\u001b[0m \u001b[0mWhere\u001b[0m \u001b[0mP\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0ma\u001b[0m \u001b[0mprojection\u001b[0m \u001b[0mof\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mfields\u001b[0m \u001b[0monto\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mdata\u001b[0m \u001b[0mspace\u001b[0m\u001b[1;33m.\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 307\u001b[0m \"\"\"\n\u001b[1;32m--> 308\u001b[1;33m \u001b[1;32mif\u001b[0m \u001b[0mu\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m:\u001b[0m \u001b[0mu\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mprob\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfields\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 309\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mprojectFields\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[0;32m 310\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/ProblemMT1D/Problems.pyc\u001b[0m in \u001b[0;36mfields\u001b[1;34m(self, m)\u001b[0m\n\u001b[0;32m 90\u001b[0m '''\n\u001b[0;32m 91\u001b[0m \u001b[1;31m# Set the current model\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 92\u001b[1;33m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcurModel\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mm\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 93\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 94\u001b[0m \u001b[0mF\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mFieldsMT\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[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msurvey\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/Problem.pyc\u001b[0m in \u001b[0;36mcurModel\u001b[1;34m(self, value)\u001b[0m\n\u001b[0;32m 76\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[1;31m# it is the same!\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 77\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mPropMap\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 78\u001b[1;33m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_curModel\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmapping\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvalue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 79\u001b[0m \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 80\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_curModel\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mModels\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mModel\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvalue\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmapping\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/PropMaps.pyc\u001b[0m in \u001b[0;36m__call__\u001b[1;34m(self, vec)\u001b[0m\n\u001b[0;32m 254\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 255\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m__call__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvec\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 256\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mPropModel\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvec\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 257\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 258\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m__contains__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mval\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/PropMaps.pyc\u001b[0m in \u001b[0;36m__init__\u001b[1;34m(self, propMap, vector)\u001b[0m\n\u001b[0;32m 117\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpropMap\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpropMap\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 118\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mvector\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mvector\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 119\u001b[1;33m \u001b[1;32massert\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mvector\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m==\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnP\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 120\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 121\u001b[0m \u001b[1;33m@\u001b[0m\u001b[0mproperty\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;31mTypeError\u001b[0m: object of type 'FieldsMT' has no len()" - ] - } - ], + "outputs": [], "source": [ - "survey.dpred(f0)" + "dMf_dsig = problem.mesh.getFaceInnerProductDeriv(problem.curModel.sigma)(u0) * problem.curModel.sigmaDeriv\n", + "dsig_dm = self.curModel.sigmaDeriv\n" ] }, { - "cell_type": "raw", + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "problem.mesh.getFaceInnerProductDeriv(problem.curModel.sigma)(u0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "problem.mesh.getEdgeInnerProductDeriv(problem.curModel.sigma)(u0[1::])" + ] + }, + { + "cell_type": "code", + "execution_count": null, "metadata": { "collapsed": true }, + "outputs": [], "source": [] } ], diff --git a/simpegMT/BaseMT.py b/simpegMT/BaseMT.py index fbc87187..31b28300 100644 --- a/simpegMT/BaseMT.py +++ b/simpegMT/BaseMT.py @@ -23,3 +23,100 @@ class BaseMTProblem(BaseFDEMProblem): # Use the forward and devs from BaseFDEMProblem # Might need to add more stuff here. + def Jvec(self, m, v, f=None): + if f is None: + f = self.fields(m) + + self.curModel = m + + Jv = self.dataPair(self.survey) + + for freq in self.survey.freqs: + dA_du = self.getA(freq) # + dA_duI = self.Solver(dA_du, **self.solverOpts) + + for src in self.survey.getSrcByFreq(freq): + ftype = self._fieldType + 'Solution' + u_src = f[src, ftype] + dA_dm = self.getADeriv_m(freq, u_src, v) + dRHS_dm = self.getRHSDeriv_m(freq, v) + if dRHS_dm is None: + du_dm = dA_duI * ( - dA_dm ) + else: + du_dm = dA_duI * ( - dA_dm + dRHS_dm ) + for rx in src.rxList: + # df_duFun = u.deriv_u(rx.fieldsUsed, m) + if 'e' in self._fieldType: + projField = 'b' + elif 'b' in self._fieldType: + projField = 'e' + df_duFun = getattr(f, '_%sDeriv_u'%projField, None) + df_du = df_duFun(src, du_dm, adjoint=False) + if df_du is not None: + du_dm = df_du + + df_dmFun = getattr(f, '_%sDeriv_m'%projField, None) + df_dm = df_dmFun(src, v, adjoint=False) + if df_dm is not None: + du_dm += df_dm + + P = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m + + + Jv[src, rx] = P(du_dm) + + return Utils.mkvc(Jv) + + def Jtvec(self, m, v, f=None): + if f is None: + f = self.fields(m) + + self.curModel = m + + # Ensure v is a data object. + if not isinstance(v, self.dataPair): + v = self.dataPair(self.survey, v) + + Jtv = np.zeros(m.size) + + for freq in self.survey.freqs: + AT = self.getA(freq).T + ATinv = self.Solver(AT, **self.solverOpts) + + for src in self.survey.getSrcByFreq(freq): + ftype = self._fieldType + 'Solution' + u_src = f[src, ftype] + + for rx in src.rxList: + PTv = rx.projectFieldsDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt u, need possibility wrt m + + df_duTFun = getattr(f, '_%sDeriv_u'%rx.projField, None) + df_duT = df_duTFun(src, PTv, adjoint=True) + if df_duT is not None: + dA_duIT = ATinv * df_duT + else: + 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) + + if dRHS_dmT is None: + du_dmT = - dA_dmT + else: + du_dmT = -dA_dmT + dRHS_dmT + + df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None) + dfT_dm = df_dmFun(src, PTv, adjoint=True) + if dfT_dm is not None: + du_dmT += dfT_dm + + real_or_imag = rx.projComp + if real_or_imag == 'real': + Jtv += du_dmT.real + elif real_or_imag == 'imag': + Jtv += - du_dmT.real + else: + raise Exception('Must be real or imag') + + return Jtv \ No newline at end of file diff --git a/simpegMT/FieldsMT.py b/simpegMT/FieldsMT.py index 97edd224..0fe0f3be 100644 --- a/simpegMT/FieldsMT.py +++ b/simpegMT/FieldsMT.py @@ -9,15 +9,88 @@ from simpegEM.Utils.EMUtils import omega ############## class FieldsMT(Problem.Fields): """Field Storage for a MT survey.""" - knownFields = {'b_px': 'F','b_py': 'F', 'e_px': 'E','e_py': 'E','b_1d':'E','e_1d':'F'} + knownFields = {} dtype = complex - def _b_1dDeriv_u(self,src,v,adjoint=False): - """ - The derivative of b_1d wrt u - """ - nG = self.mesh.nodalGrad - if adjoint: - return - 1./( 1j*omega(src.freq) ) * ( nG.T * v) - return - 1./( 1j*omega(src.freq) ) * ( nG * v) \ No newline at end of file +class FieldsMT_1D(FieldsMT): + """ + Fields storage for the 1D MT solution. + """ + knownFields = {'e_1dSolution':'F'} + aliasFields = { + 'e_1d' : ['e_1dSolution','F','_e'], + 'e_1dPrimary' : ['e_1dSolution','F','_ePrimary'], + 'e_1dSecondary' : ['e_1dSolution','F','_eSecondary'], + 'b_1d' : ['e_1dSolution','E','_b'], + 'b_1dPrimary' : ['e_1dSolution','E','_bPrimary'], + 'b_1dSecondary' : ['e_1dSolution','E','_bSecondary'] + } + + def __init__(self,mesh,survey,**kwargs): + FieldsMT.__init__(self,mesh,survey,**kwargs) + + def _ePrimary(self, eSolution, srcList): + ePrimary = np.zeros_like(eSolution) + for i, src in enumerate(srcList): + ep = src.ePrimary(self.survey.prob) + if ep is not None: + ePrimary[:,i] = ep[:,-1] + return ePrimary + + def _eSecondary(self, eSolution, srcList): + return eSolution + + def _e(self, eSolution, srcList): + return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList) + + def _eDeriv_u(self, src, v, adjoint = False): + return None + + def _eDeriv_m(self, src, v, adjoint = False): + # assuming primary does not depend on the model + return None + + def _bPrimary(self, eSolution, srcList): + bPrimary = np.zeros([self.survey.mesh.nE,eSolution.shape[1]], dtype = complex) + for i, src in enumerate(srcList): + bp = src.bPrimary(self.survey.prob) + if bp is not None: + bPrimary[:,i] += bp[:,-1] + return bPrimary + + def _bSecondary(self, eSolution, srcList): + C = self.mesh.nodalGrad + b = (C * eSolution) + for i, src in enumerate(srcList): + b[:,i] *= - 1./(1j*omega(src.freq)) + # There is no magnetic source in the MT problem + # S_m, _ = src.eval(self.survey.prob) + # if S_m is not None: + # b[:,i] += 1./(1j*omega(src.freq)) * S_m + return b + + def _b(self, eSolution, srcList): + return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList) + + def _bSecondaryDeriv_u(self, src, v, adjoint = False): + C = self.mesh.nodalGrad + if adjoint: + return - 1./(1j*omega(src.freq)) * (C.T * v) + return - 1./(1j*omega(src.freq)) * (C * v) + + def _bSecondaryDeriv_m(self, src, v, adjoint = False): + S_mDeriv, _ = src.evalDeriv(self.survey.prob, adjoint) + S_mDeriv = S_mDeriv(v) + if S_mDeriv is not None: + return 1./(1j * omega(src.freq)) * S_mDeriv + return None + + def _bDeriv_u(self, src, v, adjoint=False): + # Primary does not depend on u + return self._bSecondaryDeriv_u(src, v, adjoint) + + def _bDeriv_m(self, src, v, adjoint=False): + # Assuming the primary does not depend on the model + return self._bSecondaryDeriv_m(src, v, adjoint) + diff --git a/simpegMT/ProblemMT1D/Problems.py b/simpegMT/ProblemMT1D/Problems.py index 51982255..7ffe46fa 100644 --- a/simpegMT/ProblemMT1D/Problems.py +++ b/simpegMT/ProblemMT1D/Problems.py @@ -3,7 +3,7 @@ from SimPEG import mkvc from scipy.constants import mu_0 from simpegMT.BaseMT import BaseMTProblem from simpegMT.SurveyMT import SurveyMT -from simpegMT.FieldsMT import FieldsMT +from simpegMT.FieldsMT import FieldsMT_1D from simpegMT.DataMT import DataMT from simpegMT.Utils.MT1Danalytic import getEHfields import numpy as np @@ -18,19 +18,19 @@ class eForm_psField(BaseMTProblem): """ # From FDEMproblem: Used to project the fields. Currently not used for MTproblem. - _fieldType = 'e' + _fieldType = 'e_1d' _eqLocs = 'EF' def __init__(self, mesh, **kwargs): BaseMTProblem.__init__(self, mesh, **kwargs) + self.fieldsPair = FieldsMT_1D - def getA(self, freq,): + def getA(self, freq): """ Function to get the A matrix. :param float freq: Frequency - :param logic full: Return full A or the inner part :rtype: scipy.sparse.csr_matrix :return: A """ @@ -54,11 +54,13 @@ 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 if adjoint: - return 1j * omega(freq) * ( dsig_dm.T * ( dMf_dsig.T * v ) ) - - return 1j * omega(freq) * ( dMf_dsig * ( dsig_dm * v ) ) + return 1j * omega(freq) * ( dMf_dsig.T * v ) + # Note: output has to be nN/nF, not nC/nE. + return 1j * omega(freq) * ( (dMf_dsig * dMf_dsig.T)**(1/2) * v) def getRHS(self, freq): """ @@ -73,7 +75,7 @@ class eForm_psField(BaseMTProblem): S_e = Src.S_e(self) return -1j * omega(freq) * S_e - def getRHSderiv_m(self, freq, u, v, adjoint=False): + def getRHSDeriv_m(self, freq, v, adjoint=False): """ The derivative of the RHS wrt sigma """ @@ -91,7 +93,7 @@ class eForm_psField(BaseMTProblem): # Set the current model self.curModel = m - F = FieldsMT(self.mesh, self.survey) + F = FieldsMT_1D(self.mesh, self.survey) for freq in self.survey.freqs: if self.verbose: startTime = time.time() @@ -110,12 +112,12 @@ class eForm_psField(BaseMTProblem): # Store the fields # NOTE: only store - F[Src, 'e_1d'] = e[:,1] # Only storing the yx polarization as 1d + F[Src, 'e_1dSolution'] = e[:,1] # Only storing the yx polarization as 1d # F[Src, 'e_py'] = 0*e[:,0] # Note curl e = -iwb so b = -curl e /iw - b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) ) + # b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) ) # F[Src, 'b_px'] = 0*b[:,0] - F[Src, 'b_1d'] = b[:,1] + # F[Src, 'b_1d'] = b[:,1] if self.verbose: print 'Ran for {:f} seconds'.format(time.time()-startTime) sys.stdout.flush() diff --git a/simpegMT/SurveyMT.py b/simpegMT/SurveyMT.py index 77b2cfd6..491966dd 100644 --- a/simpegMT/SurveyMT.py +++ b/simpegMT/SurveyMT.py @@ -83,7 +83,7 @@ class RxMT(Survey.BaseRx): """Component projection (real/imag)""" return self.knownRxTypes[self.rxType][1] - def projectFields(self, src, mesh, u): + def projectFields(self, src, mesh, f): ''' Project the fields and return the ''' @@ -91,8 +91,8 @@ class RxMT(Survey.BaseRx): if self.projType is 'Z1D': Pex = mesh.getInterpolationMat(self.locs,'Fx') Pbx = mesh.getInterpolationMat(self.locs,'Ex') - ex = Pex*mkvc(u[src,'e_1d'],2) - bx = Pbx*mkvc(u[src,'b_1d'],2)/mu_0 + ex = Pex*mkvc(f[src,'e_1d'],2) + bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0 f_part_complex = ex/bx # elif self.projType is 'Z2D': elif self.projType is 'Z3D': @@ -103,14 +103,14 @@ class RxMT(Survey.BaseRx): Pby = mesh.getInterpolationMat(self.locs,'Fy') # Get the fields at location # px: x-polaration and py: y-polaration. - ex_px = Pex*u[src,'e_px'] - ey_px = Pey*u[src,'e_px'] - ex_py = Pex*u[src,'e_py'] - ey_py = Pey*u[src,'e_py'] - hx_px = Pbx*u[src,'b_px']/mu_0 - hy_px = Pby*u[src,'b_px']/mu_0 - hx_py = Pbx*u[src,'b_py']/mu_0 - hy_py = Pby*u[src,'b_py']/mu_0 + 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 # Make the complex data if 'zxx' in self.rxType: f_part_complex = (ex_px*hy_py - ex_py*hy_px)/(hx_px*hy_py - hx_py*hy_px) @@ -140,7 +140,7 @@ 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 - deriv_complex = Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pex*v) - 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._b_1dDeriv_u(src,v)/mu_0) + deriv_complex = Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pex*v) - 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) # elif self.projType is 'Z2D elif self.projType is 'Z3D': pass @@ -156,7 +156,8 @@ class RxMT(Survey.BaseRx): ############### ### Sources ### ############### -class srcMT(Survey.BaseSrc): +# Note: Should like inheret from FDEM +class srcMT(SrcFDEM): # Survey.BaseSrc): ''' Sources for the MT problem. Use the SimPEG BaseSrc, since the source fields share properties with the transmitters. @@ -205,7 +206,11 @@ class srcMT_polxy_1Dprimary(srcMT): def bPrimary(self,problem): # Project ePrimary to bPrimary # Satisfies the primary(background) field conditions - bBG_bp = (- problem.mesh.edgeCurl * self.ePrimary )/( 1j*omega(freq) ) + if problem.mesh.dim == 1: + C = problem.mesh.nodalGrad + elif problem.mesh.dim == 3: + C = problem.mesh.edgeCurl + bBG_bp = (- C * self.ePrimary(problem) )/( 1j*omega(self.freq) ) return bBG_bp def S_e(self,problem): @@ -231,15 +236,18 @@ class srcMT_polxy_1Dprimary(srcMT): def S_eDeriv(self, problem, v, adjoint = False): # Need to deal with if problem.mesh.dim == 1: - pass + # Need to use the faceInnerProduct + MsigmaDeriv = problem.mesh.getFaceInnerProductDeriv(problem.curModel.sigma)(self.ePrimary(problem)[:,-1]) * problem.curModel.sigmaDeriv + MsigmaDeriv = ( MsigmaDeriv * MsigmaDeriv.T)**2 if problem.mesh.dim == 2: pass if problem.mesh.dim == 3: - MesigmaDeriv = problem.MeSigmaDeriv(self.ePrimary(problem)) + MsigmaDeriv = problem.MeSigmaDeriv(self.ePrimary(problem)) if adjoint: - return MesigmaDeriv.T * v + return MsigmaDeriv.T * v else: - return MesigmaDeriv * v + # Moved the v in front to make the multi work + return MsigmaDeriv * v ##############