Added a note book with MT1D derivative testing

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, "cLine truncated