From 8ed4d41b2d4622f73caf567b6a122b46a0a97a13 Mon Sep 17 00:00:00 2001 From: GudniRos Date: Fri, 19 Jun 2015 09:58:03 -0700 Subject: [PATCH] Added a note book with MT1D derivative testing --- notebooks/Derivative test MT1D.ipynb | 1 + 1 file changed, 1 insertion(+) create mode 100644 notebooks/Derivative test MT1D.ipynb diff --git a/notebooks/Derivative test MT1D.ipynb b/notebooks/Derivative test MT1D.ipynb new file mode 100644 index 00000000..33a308bc --- /dev/null +++ b/notebooks/Derivative test MT1D.ipynb @@ -0,0 +1 @@ +{"nbformat_minor": 0, "cells": [{"source": "Testing derivaties for 1D MT problem.\n\nEspecially the rx.projectFieldsDeriv", "cell_type": "markdown", "metadata": {}}, {"execution_count": 1, "cell_type": "code", "source": "import SimPEG as simpeg\nimport simpegEM as simpegem, simpegMT as simpegmt\nfrom SimPEG.Utils import meshTensor\nimport numpy as np", "outputs": [], "metadata": {"collapsed": true, "trusted": true}}, {"execution_count": 3, "cell_type": "code", "source": "# Setup the problem\nsigmaHalf = 1e-2\n# Frequency\nnFreq = 33\nfreqs = np.logspace(3,-3,nFreq)\n# Make the mesh\nct = 5\nair = meshTensor([(ct,25,1.3)])\n# coreT0 = meshTensor([(ct,15,1.2)])\n# coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))\ncore = np.concatenate( ( np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) )\nbot = meshTensor([(core[0],10,-1.3)])\nx0 = -np.array([np.sum(np.concatenate((core,bot)))])\nm1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)\n# Make the model\nsigma = np.zeros(m1d.nC) + sigmaHalf\nsigma[ m1d.gridCC > 0 ] = 1e-8\n\nrxList = []\nfor rxType in ['z1dr','z1di']:\n rxList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(np.array([0.0]),2).T,rxType))\n# Source list\nsrcList =[]\ntD = False\nif tD:\n for freq in freqs:\n srcList.append(simpegmt.SurveyMT.srcMT_polxy_1DhomotD(rxList,freq))\nelse:\n for freq in freqs:\n srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq,sigma))\n# Make the survey\nsurvey = simpegmt.SurveyMT.SurveyMT(srcList)\n\n# Set the problem\nproblem = simpegmt.ProblemMT1D.eForm_psField(m1d)\nproblem.pair(survey)\n\n# Get the fields\nfields = problem.fields(sigma)\n\n# Project the data\ndata = survey.projectFields(fields)\n", "outputs": [{"output_type": "stream", "name": "stdout", "text": "Project at freq: 1.000e+03\nProject at freq: 6.494e+02\nProject at freq: 4.217e+02\nProject at freq: 2.738e+02\nProject at freq: 1.778e+02\nProject at freq: 1.155e+02\nProject at freq: 7.499e+01\nProject at freq: 4.870e+01\nProject at freq: 3.162e+01\nProject at freq: 2.054e+01\nProject at freq: 1.334e+01\nProject at freq: 8.660e+00\nProject at freq: 5.623e+00\nProject at freq: 3.652e+00\nProject at freq: 2.371e+00\nProject at freq: 1.540e+00\nProject at freq: 1.000e+00\nProject at freq: 6.494e-01\nProject at freq: 4.217e-01\nProject at freq: 2.738e-01\nProject at freq: 1.778e-01\nProject at freq: 1.155e-01\nProject at freq: 7.499e-02\nProject at freq: 4.870e-02\nProject at freq: 3.162e-02\nProject at freq: 2.054e-02\nProject at freq: 1.334e-02\nProject at freq: 8.660e-03\nProject at freq: 5.623e-03\nProject at freq: 3.652e-03\nProject at freq: 2.371e-03\nProject at freq: 1.540e-03\nProject at freq: 1.000e-03\n"}], "metadata": {"collapsed": false, "trusted": true}}, {"source": "We need calculate this derivative. \n\\begin{align}\n\\underbrace{\\frac{\\partial P(f(u(m)),m^{fix})}{\\partial f}}_{Rx}\n\\end{align}\n\nUse the rule\n\\begin{align}\n\\frac{d}{dx}\\left( \\frac{a(x)}{b(x)} \\right) = \\frac{\\frac{d }{dx} a(x) b(x) - a(x)\\frac{d }{dx} b(x) }{ b(x)^2 }\n\\end{align}\n\nIn the case of the 1D MT problem the data is calculated as \n\\begin{align}\nMT1Ddata = P(f(m)) &= \\frac{P_{ex} f_e(src,m)}{P_{bx} f_b(src,m) \\frac{1}{\\mu_0}} = \\frac{P_e u}{P_b f_b(u)} \\\\\n\\frac{\\partial P(f(m))}{\\partial u} v &= \\frac{P_e}{P_b \\frac{1}{mu_0} f_b(u)}v - \\frac{P_e u}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)^2} P_b \\frac{1}{mu_0} \\frac{d f_b}{du} v\n\\end{align}\nwhere u is the fields that we solve for. \n\\begin{align}\n\\frac{d f_b}{du} = - \\frac{1}{i \\omega} \\nabla \n\\end{align}\n\n", "cell_type": "markdown", "metadata": {}}, {"execution_count": 4, "cell_type": "code", "source": "# Unused code &= \\frac{ P_{ex} P_{bx} \\frac{1}{\\mu_0} \\left( f_b(src,m) - f_e(src,m) \\right) } { \\left(P_{bx}f_b(src,m) \\frac{1}{\\mu_0} \\right)^2 }", "outputs": [], "metadata": {"collapsed": true, "trusted": true}}, {"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\\end{align}\n\n", "cell_type": "markdown", "metadata": {}}, {"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\\end{align}\n", "cell_type": "markdown", "metadata": {}}, {"execution_count": 6, "cell_type": "code", "source": "# def projectFields(self, src, mesh, u):\n# '''\n# Project the fields and return the\n# '''\n\n# if self.projType is 'Z1D':\n# Pex = mesh.getInterpolationMat(self.locs,'Fx')\n# Pbx = mesh.getInterpolationMat(self.locs,'Ex')\n# ex = Pex*mkvc(u[src,'e_1d'],2)\n# bx = Pbx*mkvc(u[src,'b_1d'],2)/mu_0\n# f_part_complex = ex/bx\n# real_or_imag = self.projComp\n# f_part = getattr(f_part_complex, real_or_imag)\n# return f_part", "outputs": [], "metadata": {"collapsed": true, "trusted": true}}, {"execution_count": 7, "cell_type": "code", "source": "# Initate things for the derivs Test\nsrc = survey.srcList[0]\nrx = src.rxList[0]\nu0 = np.random.randn(m1d.nN)+np.random.randn(m1d.nN)*1j\nf0 = problem.fieldsPair(m1d,survey)\nf0[src,'e_1d'] = u0\nf0[src,'b_1d'] = -1/(1j*simpegem.Utils.EMUtils.omega(src.freq))*m1d.nodalGrad*u0", "outputs": [], "metadata": {"collapsed": false, "trusted": true}}, {"execution_count": 8, "cell_type": "code", "source": "u0", "outputs": [{"execution_count": 8, "output_type": "execute_result", "data": {"text/plain": "array([ 1.31527642e-01+1.45769269j, 5.95564562e-01+0.42028906j,\n -9.59349955e-01-1.11757166j, -1.40416822e+00-0.41446777j,\n 3.13493260e+00+0.03428491j, -1.51344609e+00-1.00877232j,\n 3.24356486e-01+0.65648507j, 2.67095094e-01+1.13805012j,\n 1.73664285e+00-0.65728832j, -8.97426739e-01+0.2288435j ,\n 4.80610771e-01+0.49953842j, 5.13958825e-01+1.0780965j ,\n -4.99094522e-01-1.4264564j , 1.70175950e-01-1.02378638j,\n -7.24199015e-01-0.8156263j , -1.80125549e+00+0.24188448j,\n -1.34681899e+00+0.30820186j, 5.48408346e-01-0.9591145j ,\n 3.94350988e-01+0.82075162j, -4.82329564e-01-0.0788764j ,\n 1.86240416e+00+0.52052826j, -1.29718885e+00+1.18797306j,\n -2.97912489e-01+1.46517189j, 3.70500990e-01-1.38665192j,\n 6.63769486e-01-0.37204627j, -6.73249451e-01+0.20464947j,\n 8.01428477e-01-0.04983355j, -3.24457843e-01-1.19101563j,\n 7.19176651e-01-0.65664744j, 1.15011866e+00-0.32210478j,\n 4.88403249e-01+1.13171804j, -2.21457288e+00+1.18869299j,\n 1.88662905e+00+0.15674997j, 8.58485055e-01-1.99845108j,\n -1.05007147e+00-0.85845451j, 8.84874438e-01+1.30883555j,\n -1.04745228e+00+0.88262128j, 3.11484698e-02-0.58782814j,\n 1.05803715e+00-1.16282716j, 6.13611734e-01-0.47413795j,\n 1.10311284e+00+0.72486708j, 6.35852621e-02+0.58318784j,\n 9.12335549e-01-0.93814765j, -1.20727292e-01-0.61301401j,\n -2.21251950e-01-0.26118816j, -7.34294352e-02+1.20919403j,\n 1.34397029e-01-0.73505721j, 5.47378206e-01+1.08679626j,\n -7.48448764e-01+1.68429954j, -2.87751855e-01+0.78277938j,\n 1.76210804e+00+0.39999875j, 6.35728507e-03+0.16879203j,\n 1.27904302e-01-0.97850792j, -1.53435472e+00-0.19451526j,\n 5.85361570e-01-1.95740681j, -3.67260588e-01-0.84619661j,\n 2.01779544e-01+0.13567872j, -8.22524831e-01+0.01858287j,\n 1.19370940e-02-0.23221477j, 2.07497023e-01-1.25269657j,\n 5.93470657e-01-0.0706178j , -1.51051227e+00-3.47733507j,\n -6.65438952e-02+0.89626112j, -4.93838758e-01-0.91112609j,\n -8.95717984e-01+0.27883278j, -8.80481345e-01-1.22713707j,\n -7.39150759e-01-1.54478248j, -1.93704026e-01-0.4538569j ,\n -7.72571123e-01-2.00964865j, 1.80349016e-01-0.89268155j,\n 1.71229837e-01-1.63128839j, -3.37939238e-01+1.44228476j,\n 9.35523779e-01+0.2496153j , -1.05481115e-02+0.08276521j,\n -9.39380532e-01+1.72983636j, 1.59082067e-01-0.42370157j,\n -4.55304817e-01+0.5073936j , 6.71939434e-01-0.22122061j,\n 1.07996679e-01+0.87859609j, 4.35551829e-01+0.70888035j,\n -7.19405065e-01-0.83330376j, 2.47144817e-01+0.90729062j,\n -7.34323684e-01+0.29189076j, 2.37552747e-01+0.04423161j,\n -1.00953226e+00-0.05547298j, 3.92231406e-01+1.38975434j,\n 1.66953476e+00-0.50548695j, -1.53410698e+00+2.2006748j ,\n 1.45602553e+00+0.67836125j, -7.79227347e-01-1.17293762j,\n -1.82992402e+00+0.94358606j, -3.57000836e-01-0.80496331j,\n -2.05189452e-02-0.35254849j, 1.58440553e+00+0.115809j ,\n 7.87537349e-03+0.57223829j, -1.46079714e-01-0.29568146j,\n -7.20989615e-01-0.44016013j, -8.08962609e-04-1.02039032j,\n 6.38218081e-02+1.07763744j, -1.57387983e+00-0.05508794j,\n 1.87524429e-01-0.54253528j, 3.26409352e-01+0.97574936j,\n 1.04886615e-02-0.32455216j, -6.68064344e-01+0.78639857j,\n 1.68785052e+00+0.19551552j, 7.86879974e-03+0.59927633j,\n 1.45602001e-01+0.9511198j , -1.08954477e+00+0.90133763j,\n 2.60243913e-01-1.32960299j, 3.79728167e-01+0.58389849j,\n 1.07168795e+00-0.69766928j, 9.98178067e-02-0.43534991j,\n -5.48376859e-02+0.89162356j, 2.43827158e+00-1.63902447j,\n -2.99686569e-01+0.70279198j, -1.46196259e+00+0.43685617j,\n 1.11822416e+00-0.00921865j, -3.70268989e-01+0.24181441j,\n -4.07579074e-01-0.94996173j, 8.33593718e-01+1.41780647j,\n 1.93613225e+00-1.3899132j , 1.24514067e+00-0.47260864j,\n 1.04710865e+00-0.91725898j, -6.72785198e-01+1.49358496j,\n 5.95522184e-01-0.16215232j, -3.62338353e-01+0.71558586j,\n -2.82429978e-01-0.98185035j, -4.45829454e-01-1.13774383j,\n 8.64789272e-01+2.18406621j, -2.12000531e-01-0.22423491j,\n 1.12026961e+00-0.99741678j, 2.46342968e-01+1.64555958j,\n -1.76229423e+00-2.57884939j, -9.07286746e-01+0.22959098j,\n 1.77930515e+00-0.37427123j, 2.31337484e-01-0.30702036j,\n 3.89419058e-01+1.56140629j, 1.37904088e+00-0.50908746j,\n 4.14744086e-01-1.02556973j, -7.65016201e-01-1.26797832j,\n 1.57798539e-01+2.76709369j, -1.46174568e+00+0.58091253j,\n 9.55595531e-01-0.46923076j, 1.42420672e+00-0.41777273j,\n 3.21406204e-01+0.19545789j, 3.74682200e-02+0.50264904j,\n -1.00595741e+00+0.22848755j, -1.49770231e+00-0.96416412j,\n 1.08378696e-01-1.24860747j, -9.89889158e-01+1.03405002j,\n -1.24741563e-01-0.33418193j, -3.50749355e-03-1.46353074j,\n -1.53007259e+00+0.35321057j, 3.21799404e-01-0.21669782j,\n 6.23328316e-01+0.75442761j, -4.20368367e-02+2.09823917j,\n -2.25095360e-01-1.80373873j, 3.98125417e-01-0.52445374j,\n 3.01321183e-01+0.13731416j, 1.16309011e-01+0.79141987j,\n -1.81622577e+00-0.89752232j, 5.18055667e-01-0.37565116j,\n -1.16274875e-01+0.66532497j, -6.34496238e-02+0.69392704j,\n -3.62778442e-01-0.22030796j, -8.58544163e-02-0.54046723j,\n 1.18322965e-01+0.31753671j, -6.41195322e-01+1.06919709j,\n -3.15898746e-01-0.72329593j, 3.27550922e-01-0.92963881j,\n -1.02195394e+00+0.03164951j, -4.80697404e-01+0.61650807j,\n -1.12332321e+00+0.4259311j , -5.54074221e-01-0.10316683j,\n 3.38637263e-01+0.48615485j, 9.99811343e-01-0.02015324j,\n 4.97948459e-02-1.70052266j, 5.46396769e-01+1.14898578j,\n 5.95207612e-01+0.04977481j, -1.00540168e+00-0.25387621j,\n -1.69969413e+00+1.43918166j, 9.53952015e-01+0.89194531j,\n 8.89441769e-01-0.30392719j, -8.99616738e-01+0.32580947j,\n -2.95405092e-01+1.03106057j, -4.66840877e-01+0.44651271j,\n -1.12363993e+00-0.75633937j, -9.43353751e-02-0.14289243j,\n -1.43414286e+00-0.42348461j, 1.53244295e+00+0.2365174j ,\n 9.41763652e-01-0.70837202j, -1.88556648e-01+0.04127635j,\n 8.38348832e-01-2.01673488j, 1.51513292e-01+1.4245365j ,\n 1.76371601e+00-0.24490174j, -1.04350373e+00-0.50657075j,\n 6.41590202e-01-0.49815238j, -2.31204323e-01-0.21431594j,\n -8.63829119e-02+0.02111506j, -3.16823128e-01-0.92694377j,\n -2.81923631e-01+1.14107373j, 4.65817962e-01-1.7604962j ,\n 1.30020748e-02+0.49822947j, -4.04600847e-01+0.18272931j,\n 1.06840959e-02+0.32165821j, 4.31796239e-01-1.15288956j])"}, "metadata": {}}], "metadata": {"scrolled": true, "collapsed": false, "trusted": true}}, {"source": "#rx.projectFieldsDeriv(src,m1d,fields,v)\n", "cell_type": "raw", "metadata": {}}, {"execution_count": 10, "cell_type": "code", "source": "# Run a test\ndef 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)\nsimpeg.Tests.checkDerivative(fun,u0,num=3,plotIt=False)", "outputs": [{"output_type": "stream", "name": "stdout", "text": "==================== checkDerivative ====================\niter h |ft-f0| |ft-f0-h*J0*dx| Order\n---------------------------------------------------------\n 0 1.00e-01 6.154e-03 8.309e-03 nan\n 1 1.00e-02 5.944e-04 8.099e-04 1.011\n 2 1.00e-03 5.915e-05 8.070e-05 1.002\n*********************************************************\n<<<<<<<<<<<<<<<<<<<<<<<<< FAIL! >>>>>>>>>>>>>>>>>>>>>>>>>\n*********************************************************\nDid you put your clever trousers on today?\n\n"}, {"execution_count": 10, "output_type": "execute_result", "data": {"text/plain": "False"}, "metadata": {}}], "metadata": {"collapsed": false, "trusted": true}}, {"execution_count": null, "cell_type": "code", "source": "print rx.projectFieldsDeriv(src,m1d,f0,u0)\nprint m1d.nF\nprint m1d.nN", "outputs": [], "metadata": {"collapsed": false, "trusted": true}}, {"execution_count": null, "cell_type": "code", "source": "fields._b_1dDeriv_u(src,u0)", "outputs": [], "metadata": {"collapsed": false, "trusted": true}}, {"execution_count": null, "cell_type": "code", "source": "%debug", "outputs": [], "metadata": {"collapsed": false, "trusted": true}}, {"execution_count": null, "cell_type": "code", "source": "", "outputs": [], "metadata": {"collapsed": true, "trusted": true}}], "nbformat": 4, "metadata": {"kernelspec": {"display_name": "Python 2", "name": "python2", "language": "python"}}} \ No newline at end of file