diff --git a/docs/api_PF.rst b/docs/api_PF.rst index 2ac74424..d045dd08 100644 --- a/docs/api_PF.rst +++ b/docs/api_PF.rst @@ -100,7 +100,7 @@ Since most materials in the earth have lower permeability than \\(\\mu_0\\), us Since we compute secondary field based on the earth field, which can be different from different locations in the world, we can expect different anomalous responses in different locations in the earth. For instance, assume we have two susceptible spheres, which are exactly same. However, anomalous responses in Seoul and Vancouver are going to be different. -.. plot :: /home/seogi/Documents/simpegpf/docs/figures/figure1.py +.. plot :: examples/Mag_back_1.py Since we can measure total fields ( \\(\\vec{B}\\)), and usually have reasonably accurate earth field (\\(\\vec{B}_0\\)), we can compute anmalous fields, \\(\\vec{B}_s\\) from our observed data. If you want to download earth magnetic fields at specific location see this website (`noaa `_). diff --git a/docs/examples/Mag_back_1.py b/docs/examples/Mag_back_1.py new file mode 100644 index 00000000..b4dd1459 --- /dev/null +++ b/docs/examples/Mag_back_1.py @@ -0,0 +1,47 @@ +from simpegPF.MagAnalytics import MagSphereAnalFunA, IDTtoxyz +from SimPEG import * +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable + +xr = np.linspace(-300, 300, 41) +yr = np.linspace(-300, 300, 41) +X, Y = np.meshgrid(xr, yr) +Z = np.ones((np.size(xr), np.size(yr)))*150 + + +# Bz component at Korea +inckr = -8. + 3/60 +deckr = 54. + 9/60 +btotkr = 50898.6 +Bokr = IDTtoxyz(inckr, deckr, btotkr) + +bx,by,bz = MagSphereAnalFunA(X, Y, Z,100.,0.,0.,0.,0.01,Bokr,'secondary') +Bzkr = np.reshape(bz, (np.size(xr), np.size(yr)), order='F') + +# Bz component at Canada +incca = 16. + 49/60 +decca = 70. + 19/60 +btotca = 54692.1 +Boca = IDTtoxyz(incca, decca, btotca) + +bx,by,bz = MagSphereAnalFunA(X, Y, Z,100.,0.,0.,0.,0.01,Boca,'secondary') +Bzca = np.reshape(bz, (np.size(xr), np.size(yr)), order='F') + +fig = plt.figure( figsize = (14,5) ) + +ax1 = plt.subplot(121) +dat1 = plt.imshow(Bzkr, extent=[min(xr), max(xr), min(yr), max(yr)]); +divider = make_axes_locatable(ax1) +cax1 = divider.append_axes("right", size="5%", pad=0.05) +ax1.set_xlabel('East-West (m)'); ax1.set_ylabel('South-North (m)') +plt.colorbar(dat1, cax=cax1) +ax1.set_title('$B_z$ field at Seoul, South Korea') + +ax2 = plt.subplot(122) +dat2 = plt.imshow(Bzca, extent=[min(xr), max(xr), min(yr), max(yr)]); +divider = make_axes_locatable(ax2) +cax2 = divider.append_axes("right", size="5%", pad=0.05) +ax2.set_xlabel('East-West (m)'); ax2.set_ylabel('South-North (m)') +plt.colorbar(dat2, cax=cax2) +ax2.set_title('$B_z$ field at Vancouver, Canada') +plt.show() \ No newline at end of file diff --git a/simpegPF/Tests/test_forward_PFproblem.py b/simpegPF/Tests/test_forward_PFproblem.py index 8451fe9c..e786a34d 100644 --- a/simpegPF/Tests/test_forward_PFproblem.py +++ b/simpegPF/Tests/test_forward_PFproblem.py @@ -1,23 +1,55 @@ import unittest from SimPEG import * +import matplotlib.pyplot as plt import simpegPF as PF class MagProblemTests(unittest.TestCase): def setUp(self): - M = Mesh.TensorMesh([10,10]) - mod = Model.LogModel(M) - prob = PF.Mag.MagProblem(M, mod, None) + hxind = ((5,25,1.3),(41, 12.5),(5,25,1.3)) + hyind = ((5,25,1.3),(41, 12.5),(5,25,1.3)) + hzind = ((5,25,1.3),(40, 12.5),(5,25,1.3)) + hx, hy, hz = Utils.meshTensors(hxind, hyind, hzind) + M = Mesh.TensorMesh([hx, hy, hz], [-hx.sum()/2,-hy.sum()/2,-hz.sum()/2]) + + chibkg = 0. + chiblk = 0.01 + chi = np.ones(M.nC)*chibkg + sph_ind = PF.MagAnalytics.spheremodel(M, 0., 0., 0., 100) + chi[sph_ind] = chiblk + model = PF.BaseMag.BaseMagModel(M) + prob = PF.Magnetics.MagneticsDiffSecondary(M, model) self.prob = prob self.M = M + self.chi = chi - def test_forward(self): - passed = True - self.assertTrue(passed) + def test_anal_forward(self): + data = PF.BaseMag.BaseMagData() + data.setBackgroundField(x=1., y=1., z=0.) + xr = np.linspace(-300, 300, 41) + yr = np.linspace(-300, 300, 41) + X, Y = np.meshgrid(xr, yr) + Z = np.ones((xr.size, yr.size))*150 + rxLoc = np.c_[Utils.mkvc(X), Utils.mkvc(Y), Utils.mkvc(Z)] + data.rxLoc = rxLoc + + self.prob.pair(data) + B = self.prob.fields(self.chi) + + bxa,bya,bza = PF.MagAnalytics.MagSphereAnalFunA(rxLoc[:,0],rxLoc[:,1],rxLoc[:,2],100.,0.,0.,0.,0.01,np.array([1.,1.,0.]),'secondary') + + dpred = data.projectFieldsAsVector(B) + err = np.linalg.norm(dpred-np.r_[bxa, bya, bza])/np.linalg.norm(np.r_[bxa, bya, bza]) + + if err > 0.05: + raise Exception('Anaytic test is failed T.T') + else: + print "Anaytic test is passed" + pass if __name__ == '__main__': unittest.main() diff --git a/simpegPF/notebooks/Jacobian.ipynb b/simpegPF/notebooks/Jacobian.ipynb new file mode 100644 index 00000000..8b8994c4 --- /dev/null +++ b/simpegPF/notebooks/Jacobian.ipynb @@ -0,0 +1,400 @@ +{ + "metadata": { + "name": "" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from SimPEG import Mesh, Problem, Utils, np, sp, Tests\n", + "from simpegPF.Magnetics import BaseMag, MagneticsDiffSecondary\n", + "from scipy.constants import mu_0\n", + "from simpegPF.MagAnalytics import spheremodel, CongruousMagBC\n", + "%pylab inline" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Populating the interactive namespace from numpy and matplotlib\n" + ] + } + ], + "prompt_number": 3 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "hxind = ((5,25,1.3),(21, 25.),(5,25,1.3))\n", + "hyind = ((5,25,1.3),(21, 25.),(5,25,1.3))\n", + "hzind = ((5,25,1.3),(20, 25.),(5,25,1.3))\n", + "hx, hy, hz = Utils.meshTensors(hxind, hyind, hzind)\n", + "M3 = Mesh.TensorMesh([hx, hy, hz], [-sum(hx)/2,-sum(hy)/2,-sum(hz)/2])" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 4 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "import matplotlib.pyplot as plt\n", + "hxind = ((5,25,1.3),(41, 12.5),(5,25,1.3))\n", + "hyind = ((5,25,1.3),(41, 12.5),(5,25,1.3))\n", + "hzind = ((5,25,1.3),(40, 12.5),(5,25,1.3))\n", + "hx, hy, hz = Utils.meshTensors(hxind, hyind, hzind)\n", + "mesh = Mesh.TensorMesh([hx, hy, hz], [-hx.sum()/2,-hy.sum()/2,-hz.sum()/2])\n", + "\n", + "chibkg = 0.001\n", + "chiblk = 0.1\n", + "chi = np.ones(mesh.nC)*chibkg\n", + "sph_ind = spheremodel(mesh, 0., 0., 0., 100)\n", + "chi[sph_ind] = chiblk\n", + "model = BaseMag.BaseMagModel(mesh)\n", + "mu = (1.+chi)*mu_0\n", + "\n", + "data = BaseMag.BaseMagData()\n", + "data.setBackgroundField(x=1., y=1., z=0.)\n", + "xr = np.linspace(-300, 300, 41)\n", + "yr = np.linspace(-300, 300, 41)\n", + "X, Y = np.meshgrid(xr, yr)\n", + "Z = np.ones((xr.size, yr.size))*150\n", + "rxLoc = np.c_[Utils.mkvc(X), Utils.mkvc(Y), Utils.mkvc(Z)]\n", + "data.rxLoc = rxLoc\n", + "\n", + "prob = MagneticsDiffSecondary(mesh, model)\n", + "\n", + "prob.pair(data)\n" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 5 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "phi, B = prob.fields(chi)\n", + "# mesh.plotSlice(B, 'F', view='vec', showIt=True)\n", + "dpred = data.dpred(chi, u=B)" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 4 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "test = {'B': B, 'u': phi}" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 5 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "print test['B']" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "[ 0.001 0.00032989 0.00039923 ..., 0. 0. 0. ]\n" + ] + } + ], + "prompt_number": 7 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## ?? they are different ?" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "MfMui = mesh.getFaceInnerProduct(mu).diagonal()\n", + "MfMui_test = mesh.getFaceMass(mu).diagonal()\n", + "print np.linalg.norm(MfMui - MfMui_test)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "6.99232804893\n" + ] + } + ], + "prompt_number": 6 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "plot(MfMui/MfMui_test)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 7, + "text": [ + "[]" + ] + }, + { + "metadata": {}, + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAD9CAYAAACm2+DgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHNxJREFUeJzt3X9QVNf9//HXIsw4hIQEq4jgBBQCrMAuCm7TYLIWDYmm\nJpq0NTGMyRCHMkkTO36chkkzwXaqsY3toOYPazU2TUttZ1L101GiTNxJNFHGFGorxjFGwvLLCBEQ\nfwTZPZ8//LrfEK6IDQtr8nzM7MzuPWf3vu9R7ou95yxrM8YYAQDwJWEjXQAAIDQREAAASwQEAMAS\nAQEAsERAAAAsERAAAEsDBsTFixflcrnkdDplt9tVWlrar8+f/vQnORwOZWVl6a677tLhw4cDbZWV\nlUpLS1NKSopWr1499NUDAILGdq3PQZw/f16RkZHq7e1VXl6eXnnlFeXl5QXa33//fdntdkVHR6uy\nslJlZWU6cOCAfD6fUlNTVVVVpfj4eOXm5qqiokLp6elBPygAwFd3zUtMkZGRkqSenh75fD7FxMT0\nab/zzjsVHR0tSXK5XGpsbJQkVVdXKzk5WYmJiYqIiNDChQu1ffv2oa4fABAk1wwIv98vp9Op2NhY\nzZw5U3a7/ap9N23apDlz5kiSmpqaNHHixEBbQkKCmpqahqBkAMBwCL9Wh7CwMNXW1qqzs1MFBQXy\neDxyu939+u3du1ebN2/W/v37JUk2m21QBQy2HwCgr2D/paRBr2KKjo7W3LlzdejQoX5thw8f1pIl\nS7Rjxw7ddtttkqT4+Hh5vd5AH6/Xq4SEBMvXNsaE/O2ll14a8Rqokxqpkzqv3IbDgAHR1tamjo4O\nSdKFCxe0Z88eZWdn9+nT0NCgBQsW6I033lBycnJge05Ojo4fP676+nr19PRo69atmjdvXhAOAQAQ\nDANeYmppadHixYvl9/vl9/tVWFio/Px8bdiwQZJUXFysn//85zpz5oxKSkokSREREaqurlZ4eLjW\nr1+vgoIC+Xw+FRUVsYIJAG4g11zmGvQCbLZhe7v0VVxt7iXUUOfQuRFqlKhzqN0odQ7HuZOAAIAb\n0HCcO/lTGwAASwQEAMASAQEAsERAAAAsERAAAEsEBADAEgEBALBEQAAALBEQAABLBAQAwBIBAQCw\nREAAACwREAAASwQEAMASAQEAsERAAAAsERAAAEsEBADA0oABcfHiRblcLjmdTtntdpWWlvbr8+GH\nH+rOO+/U6NGjtWbNmj5tiYmJysrKUnZ2tqZPnz60lQMAgip8oMbRo0dr7969ioyMVG9vr/Ly8rRv\n3z7l5eUF+owZM0br1q3Ttm3b+j3fZrPJ4/EoJiZm6CsHAATVNS8xRUZGSpJ6enrk8/n6nezHjh2r\nnJwcRUREWD4/2F+qDQAIjgHfQUiS3+/X1KlTdeLECZWUlMhutw/6xW02m2bNmqVRo0apuLhYS5Ys\nsexXVlYWuO92u+V2uwe9DwD4JvB4PPJ4PMO6T5sZ5K/4nZ2dKigo0Msvv2x5Al+xYoWioqK0bNmy\nwLaWlhbFxcXp9OnTmj17ttatW6cZM2b0LcBm410GAFyn4Th3DnoVU3R0tObOnatDhw4N+sXj4uIk\nXb4MNX/+fFVXV19/hQCAETFgQLS1tamjo0OSdOHCBe3Zs0fZ2dmWfb+cZOfPn9fZs2clSefOndPu\n3buVmZk5FDUDAIbBgHMQLS0tWrx4sfx+v/x+vwoLC5Wfn68NGzZIkoqLi9Xa2qrc3Fx1dXUpLCxM\n5eXlqqur06effqoFCxZIknp7e7Vo0SLde++9wT8iAMCQGPQcRNAKYA4CAK5bSM1BAAC+WQgIAIAl\nAgIAYImAAABYIiAAAJYICACAJQICAGCJgAAAWCIgAACWCAgAgCUCAgBgiYAAAFgiIAAAlggIAIAl\nAgIAYImAAABYIiAAAJYICACAJQICAGBpwIC4ePGiXC6XnE6n7Ha7SktL+/X58MMPdeedd2r06NFa\ns2ZNn7bKykqlpaUpJSVFq1evHtrKAQBBZTPX+Nbr8+fPKzIyUr29vcrLy9Mrr7yivLy8QPvp06f1\nySefaNu2bbrtttu0bNkySZLP51NqaqqqqqoUHx+v3NxcVVRUKD09vW8Bw/DF2wDwdTMc585rXmKK\njIyUJPX09Mjn8ykmJqZP+9ixY5WTk6OIiIg+26urq5WcnKzExERFRERo4cKF2r59+xCWDgAIpvBr\ndfD7/Zo6dapOnDihkpIS2e32Qb1wU1OTJk6cGHickJCggwcPWvYtKysL3He73XK73YPaBwB8U3g8\nHnk8nmHd5zUDIiwsTLW1ters7FRBQYE8Hs+gTuA2m23QRXwxIAAA/X35l+cVK1YEfZ+DXsUUHR2t\nuXPn6tChQ4PqHx8fL6/XG3js9XqVkJBw/RUCAEbEgAHR1tamjo4OSdKFCxe0Z88eZWdnW/b98mRJ\nTk6Ojh8/rvr6evX09Gjr1q2aN2/eEJUNAAi2AS8xtbS0aPHixfL7/fL7/SosLFR+fr42bNggSSou\nLlZra6tyc3PV1dWlsLAwlZeXq66uTlFRUVq/fr0KCgrk8/lUVFTUbwUTACB0XXOZa9ALYJkrAFy3\nkFjmCgD4ZiIgAACWCAgAgCUCAgBgiYAAAFgiIAAAlggIAIAlAgIAYImAAABYIiAAAJYICACAJQIC\nAGCJgAAAWCIgAACWCAgAgCUCAgBgiYAAAFgiIAAAlggIAIClAQPi4sWLcrlccjqdstvtKi0ttez3\n7LPPKiUlRQ6HQzU1NYHtiYmJysrKUnZ2tqZPnz60lQMAgip8oMbRo0dr7969ioyMVG9vr/Ly8rRv\n3z7l5eUF+uzcuVMfffSRjh8/roMHD6qkpEQHDhyQdPlLtT0ej2JiYoJ7FACAIXfNS0yRkZGSpJ6e\nHvl8vn4n+x07dmjx4sWSJJfLpY6ODp06dSrQbowZynoBAMNkwHcQkuT3+zV16lSdOHFCJSUlstvt\nfdqbmpo0ceLEwOOEhAQ1NTUpNjZWNptNs2bN0qhRo1RcXKwlS5ZY7qOsrCxw3+12y+12/3dHAwBf\nUx6PRx6PZ1j3ec2ACAsLU21trTo7O1VQUCCPx9PvBH61dwn79u3ThAkTdPr0ac2ePVtpaWmaMWNG\nv35fDAgAQH9f/uV5xYoVQd/noFcxRUdHa+7cuTp06FCf7fHx8fJ6vYHHjY2Nio+PlyRNmDBBkjR2\n7FjNnz9f1dXVQ1EzAGAYDBgQbW1t6ujokCRduHBBe/bsUXZ2dp8+8+bN0+uvvy5JOnDggG699VbF\nxsbq/PnzOnv2rCTp3Llz2r17tzIzM4NxDACAIBjwElNLS4sWL14sv98vv9+vwsJC5efna8OGDZKk\n4uJizZkzRzt37lRycrJuuukmvfbaa5Kk1tZWLViwQJLU29urRYsW6d577w3y4QAAhorNjPAyI5vN\nxkonALhOw3Hu5JPUAABLBAQAwBIBAQCwREAAACwREAAASwQEAMASAQEAsERAAAAsERAAAEsEBADA\nEgEBALBEQAAALBEQAABLBAQAwBIBAQCwREAAACwREAAASwQEAMASAQEAsDRgQFy8eFEul0tOp1N2\nu12lpaWW/Z599lmlpKTI4XCopqYmsL2yslJpaWlKSUnR6tWrh7ZyAEBQDRgQo0eP1t69e1VbW6vD\nhw9r79692rdvX58+O3fu1EcffaTjx4/rd7/7nUpKSiRJPp9PzzzzjCorK1VXV6eKigodPXo0eEcC\nABhS4dfqEBkZKUnq6emRz+dTTExMn/YdO3Zo8eLFkiSXy6WOjg61trbq5MmTSk5OVmJioiRp4cKF\n2r59u9LT0/vt43//96seBgB8M9x+u5SVNTz7umZA+P1+TZ06VSdOnFBJSYnsdnuf9qamJk2cODHw\nOCEhQU1NTWpubu63/eDBg5b7+J//KQvcHzPGrTFj3Nd5GADw9dbe7lF7u0fJyVJu7vDs85oBERYW\nptraWnV2dqqgoEAej0dut7tPH2PMVyri2LGyr/R8APj6c/+/22UrVqwI+h4HvYopOjpac+fO1aFD\nh/psj4+Pl9frDTxubGxUQkJCv+1er1cJCQlDUDIAYDgMGBBtbW3q6OiQJF24cEF79uxRdnZ2nz7z\n5s3T66+/Lkk6cOCAbr31VsXGxionJ0fHjx9XfX29enp6tHXrVs2bNy9IhwEAGGoDXmJqaWnR4sWL\n5ff75ff7VVhYqPz8fG3YsEGSVFxcrDlz5mjnzp1KTk7WTTfdpNdee+3yC4eHa/369SooKJDP51NR\nUZHlBDUAIDTZzFedQPiqBdhsX3kOAwC+aYbj3MknqQEAlggIAIAlAgIAYImAAABYIiAAAJYICACA\nJQICAGCJgAAAWCIgAACWCAgAgCUCAgBgiYAAAFgiIAAAlggIAIAlAgIAYImAAABYIiAAAJYICACA\nJQICAGBpwIDwer2aOXOmpkyZooyMDK1du7ZfnzNnzmj+/PlyOBxyuVw6cuRIoC0xMVFZWVnKzs7W\n9OnTh756AEDQ2MwA33rd2tqq1tZWOZ1OdXd3a9q0adq2bZvS09MDfZYvX65bbrlFL774oo4dO6an\nn35aVVVVkqSkpCR98MEHiomJuXoBw/DF2wDwdTMc584B30GMHz9eTqdTkhQVFaX09HQ1Nzf36XP0\n6FHNnDlTkpSamqr6+nqdPn060M7JHwBuTOGD7VhfX6+amhq5XK4+2x0Oh958803l5eWpurpan3zy\niRobGzV27FjZbDbNmjVLo0aNUnFxsZYsWWL52mVlZYH7brdbbrf7vzoYAPi68ng88ng8w7rPAS8x\nXdHd3S23262f/exneuihh/q0nT17Vs8995xqamqUmZmpDz/8UL///e+VlZWl5uZmTZgwQadPn9bs\n2bO1bt06zZgxo28BXGICgOs2HOfOawbEpUuX9MADD+j+++/X0qVLr/mCSUlJ+ve//62oqKg+21es\nWKGoqCgtW7asbwEEBABctxGfgzDGqKioSHa7/arh0NnZqZ6eHknSxo0bdc899ygqKkrnz5/X2bNn\nJUnnzp3T7t27lZmZOcTlAwCCZcA5iP379+uNN94ILFWVpJUrV6qhoUGSVFxcrLq6Oj3xxBOy2WzK\nyMjQpk2bJEmnTp3S/PnzJUm9vb1atGiR7r333mAeCwBgCA1qDiKoBXCJCQCu24hfYgIAfHMREAAA\nSwQEAMASAQEAsERAAAAsERAAAEsEBADAEgEBALBEQAAALBEQAABLBAQAwBIBAQCwREAAACwREAAA\nSwQEAMASAQEAsERAAAAsERAAAEsEBADA0oAB4fV6NXPmTE2ZMkUZGRlau3Ztvz5nzpzR/Pnz5XA4\n5HK5dOTIkUBbZWWl0tLSlJKSotWrVw999QCAoLGZAb71urW1Va2trXI6neru7ta0adO0bds2paen\nB/osX75ct9xyi1588UUdO3ZMTz/9tKqqquTz+ZSamqqqqirFx8crNzdXFRUVfZ4rDc8XbwPA181w\nnDsHfAcxfvx4OZ1OSVJUVJTS09PV3Nzcp8/Ro0c1c+ZMSVJqaqrq6+v16aefqrq6WsnJyUpMTFRE\nRIQWLlyo7du3B+kwAABDLXywHevr61VTUyOXy9Vnu8Ph0Jtvvqm8vDxVV1frk08+UWNjo5qamjRx\n4sRAv4SEBB08eNDytcvKygL33W633G739R0FAHzNeTweeTyeYd3noAKiu7tbjzzyiMrLyxUVFdWn\n7fnnn9dzzz2n7OxsZWZmKjs7W6NGjZLNZht0EV8MCABAf1/+5XnFihVB3+c1A+LSpUt6+OGH9fjj\nj+uhhx7q137zzTdr8+bNgcdJSUmaPHmyLly4IK/XG9ju9XqVkJAwRGUDAIJtwDkIY4yKiopkt9u1\ndOlSyz6dnZ3q6emRJG3cuFH33HOPoqKilJOTo+PHj6u+vl49PT3aunWr5s2bN/RHAAAIigHfQezf\nv19vvPGGsrKylJ2dLUlauXKlGhoaJEnFxcWqq6vTE088IZvNpoyMDG3atOnyC4eHa/369SooKJDP\n51NRUVG/FUwAgNA14DLXYSmAZa4AcN1GfJkrAOCbi4AAAFgiIAAAlggIAIAlAgIAYImAAABYIiAA\nAJYICACAJQICAGCJgAAAWCIgAACWCAgAgCUCAgBgiYAAAFgiIAAAlggIAIAlAgIAYImAAABYGjAg\nvF6vZs6cqSlTpigjI0Nr167t16etrU333XefnE6nMjIytGXLlkBbYmJi4Pusp0+fPuTFAwCCZ8Dv\npG5tbVVra6ucTqe6u7s1bdo0bdu2Tenp6YE+ZWVl+vzzz7Vq1Sq1tbUpNTVVp06dUnh4uJKSkvTB\nBx8oJibm6gXwndQAcN1G/Dupx48fL6fTKUmKiopSenq6mpub+/SJi4tTV1eXJKmrq0tjxoxReHh4\noJ2TPwDcmAY9B1FfX6+amhq5XK4+25csWaIjR45owoQJcjgcKi8vD7TZbDbNmjVLOTk52rhx49BV\nDQAIuvBrd5G6u7v1yCOPqLy8XFFRUX3aVq5cKafTKY/HoxMnTmj27Nn617/+pZtvvln79+9XXFyc\nTp8+rdmzZystLU0zZszo9/plZWWB+263W263+ysdFAB83Xg8Hnk8nmHd54BzEJJ06dIlPfDAA7r/\n/vu1dOnSfu1z5szRCy+8oLvuukuSlJ+fr9WrVysnJ6dPvxUrVigqKkrLli3rWwBzEABw3UZ8DsIY\no6KiItntdstwkKS0tDRVVVVJkk6dOqVjx45p0qRJOn/+vM6ePStJOnfunHbv3q3MzMwhLh8AECwD\nvoPYt2+f7r77bmVlZclms0m6fEmpoaFBklRcXKy2tjY9+eSTamhokN/vV2lpqR577DF9/PHHWrBg\ngSSpt7dXixYtUmlpaf8CeAcBANdtOM6d17zEFGwEBABcvxG/xAQA+OYiIAAAlggIAIAlAgIAYImA\nAABYIiAAAJYICACAJQICAGCJgAAAWCIgAACWCAgAgCUCAgBgiYAAAFgiIAAAlggIAIAlAgIAYImA\nAABYIiAAAJYICACApQEDwuv1aubMmZoyZYoyMjK0du3afn3a2tp03333yel0KiMjQ1u2bAm0VVZW\nKi0tTSkpKVq9evWQFz+cPB7PSJcwKNQ5dG6EGiXqHGo3Sp3DYcCAiIiI0G9/+1sdOXJEBw4c0Kuv\nvqqjR4/26bN+/XplZ2ertrZWHo9Hy5YtU29vr3w+n5555hlVVlaqrq5OFRUV/Z57I7lR/tNQ59C5\nEWqUqHOo3Sh1DocBA2L8+PFyOp2SpKioKKWnp6u5ublPn7i4OHV1dUmSurq6NGbMGIWHh6u6ulrJ\nyclKTExURESEFi5cqO3btwfpMAAAQy18sB3r6+tVU1Mjl8vVZ/uSJUv03e9+VxMmTNDZs2f117/+\nVZLU1NSkiRMnBvolJCTo4MGDQ1Q2ACDozCCcPXvWTJs2zfz973/v1/aLX/zCPPfcc8YYYz766COT\nlJRkurq6zN/+9jfz1FNPBfr98Y9/NM8880y/50vixo0bN27/xS3YrvkO4tKlS3r44Yf1+OOP66GH\nHurX/t577+mFF16QJE2ePFlJSUk6duyYEhIS5PV6A/28Xq8SEhL6Pf9yRgAAQs2AcxDGGBUVFclu\nt2vp0qWWfdLS0lRVVSVJOnXqlI4dO6ZJkyYpJydHx48fV319vXp6erR161bNmzdv6I8AABAUNjPA\nr/D79u3T3XffraysLNlsNknSypUr1dDQIEkqLi5WW1ubnnzySTU0NMjv96u0tFSPPfaYJGnXrl1a\nunSpfD6fioqKVFpaOgyHBAAYEkG/iDWAXbt2mdTUVJOcnGxefvnlYdnn7bffbjIzM43T6TS5ubnG\nGGPa29vNrFmzTEpKipk9e7Y5c+ZMoP/KlStNcnKySU1NNW+99VZg+6FDh0xGRoZJTk42zz77bGD7\nxYsXzQ9+8AOTnJxsXC6Xqa+vH1RdTz75pBk3bpzJyMgIbBuuurZs2WJSUlJMSkqK+cMf/nBdNb70\n0ksmPj7eOJ1O43Q6zc6dO0e0RmOMaWhoMG6329jtdjNlyhRTXl4ekuN5tTpDbUwvXLhgpk+fbhwO\nh0lPTzfPP/98SI7n1eoMtfE0xpje3l7jdDrNAw88EJJjecWIBURvb6+ZPHmyOXnypOnp6TEOh8PU\n1dUFfb+JiYmmvb29z7bly5eb1atXG2OMefnll81Pf/pTY4wxR44cMQ6Hw/T09JiTJ0+ayZMnG7/f\nb4wxJjc31xw8eNAYY8z9999vdu3aZYwx5tVXXzUlJSXGGGP+8pe/mB/+8IeDquudd94x//znP/uc\nfIejrvb2djNp0iRz5swZc+bMmcD9wdZYVlZm1qxZ06/vSNVojDEtLS2mpqbGGHN5gcUdd9xh6urq\nQm48r1ZnKI7puXPnjDHGXLp0ybhcLvPuu++G3Hherc5QHM81a9aYxx57zHzve98zxoTez/oVIxYQ\n7733nikoKAg8XrVqlVm1alXQ95uYmGja2tr6bEtNTTWtra3GmMs/tKmpqcaYy8n9xXc2BQUF5v33\n3zfNzc0mLS0tsL2iosIUFxcH+hw4cMAYc/k/6be+9a1B13by5Mk+J9/hqOvPf/6z+dGPfhR4TnFx\nsamoqBh0jWVlZeaVV17p128ka/yyBx980OzZsyckx9OqzlAe03PnzpmcnBzzn//8J6TH84t1htp4\ner1ek5+fb95+++3AO4hQHcsR+1tMVp+TaGpqCvp+bTabZs2apZycHG3cuFHS5cn12NhYSVJsbKxO\nnTolSWpubu6z8upKjV/eHh8fH6j9i8cVHh6u6OhoffbZZ/9VrcGuq729/aqvdT3WrVsnh8OhoqIi\ndXR0hFSNX/z8TiiP55U6v/3tb0sKvTH1+/1yOp2KjY0N/PmdUBxPqzpDbTx/8pOf6Ne//rXCwv7/\n6TcUx1IawT/Wd2XSe7jt379fNTU12rVrl1599VW9++67/eoaqdoGEqp1lZSU6OTJk6qtrVVcXJyW\nLVs20iUFdHd36+GHH1Z5ebluvvnmPm2hNJ7d3d165JFHVF5erqioqJAc07CwMNXW1qqxsVHvvPOO\n9u7d26c9VMbzy3V6PJ6QGs9//OMfGjdunLKzs6+6xD9UxlIawYCIj48f1OckhlpcXJwkaezYsZo/\nf76qq6sVGxur1tZWSVJLS4vGjRtnWWNjY6MSEhIUHx+vxsbGftuvPOfKKq/e3l51dnYqJibmv6o1\n2HWNGTPmK/87jBs3LvAf+qmnnlJ1dXVI1Hjl8zuFhYWBz++E4nhafc4oVMdUkqKjozV37lx98MEH\nITmeX67z0KFDITWe7733nnbs2KGkpCQ9+uijevvtt1VYWBi6YzngBaggunTpkpk0aZI5efKk+fzz\nz4dlkvrcuXOmq6vLGGNMd3e3+c53vmPeeusts3z58sB1vlWrVvWbIPr888/Nxx9/bCZNmhSYIJo+\nfbo5cOCA8fv9/SaIrlznq6ioGPQktTH9r+8PR13t7e0mKSnJnDlzxnz22WeB+4Otsbm5OXD/N7/5\njXn00UdHvEa/328KCwvN0qVL+2wPtfG8Wp2hNqanT58OtJ0/f97MmDHDVFVVhdx4Xq3OlpaWkBrP\nKzweT2AOItTG8ooRXea6c+dOc8cdd5jJkyeblStXBn1/H3/8sXE4HMbhcJgpU6YE9tne3m7y8/Mt\nl5j98pe/NJMnTzapqammsrIysP3KErPJkyebH//4x4HtFy9eNN///vcDS8xOnjw5qNoWLlxo4uLi\nTEREhElISDCbN28etro2b95skpOTTXJystmyZcuga9y0aZMpLCw0mZmZJisryzz44IOBibaRqtEY\nY959911js9mMw+EILG3ctWtXyI2nVZ07d+4MuTE9fPiwyc7ONg6Hw2RmZppf/epXxpjh+7n5qnWG\n2nhe4fF4AquYQm0srxjwg3IAgG8uvlEOAGCJgAAAWCIgAACWCAgAgCUCAgBgiYAAAFj6P0f5Arjx\nzIiGAAAAAElFTkSuQmCC\n", + "text": [ + "" + ] + } + ], + "prompt_number": 7 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "def MfmuIfun(mu):\n", + " \n", + "# MfMui = mesh.getFaceMass(1./mu)\n", + " \n", + " MfMui = mesh.getFaceInnerProduct(1./mu)*1/3\n", + " # There are scaling issue for Face inner product!! and bit slower than getFaceMass\n", + " \n", + " MfMuI = 1/MfMui.diagonal()\n", + "# dMfMuI = Utils.sdiag(MfMuI**2)*mesh.aveF2CC.T*Utils.sdiag(mesh.vol*1./mu**2)\n", + " dMfMuI = Utils.sdiag(MfMuI**2)*mesh.aveF2CC.T*Utils.sdiag(mesh.vol*1./mu**2)\n", + " \n", + " return MfMuI, dMfMuI" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 8 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "d_mu = 0.8*mu\n", + "# I am not sure why it does not give us second order accuracy\n", + "Tests.checkDerivative(MfmuIfun, mu, dx=d_mu, num=5, plotIt=False)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "==================== checkDerivative ====================\n", + "iter\th\t\t|J0-Jt|\t\t|J0+h*dJ'*dx-Jt|\tOrder\n", + "---------------------------------------------------------\n", + "0\t1.00e-01\t7.240e-08\t\t1.904e-22\t\tnan" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "1\t1.00e-02\t7.240e-09\t\t1.009e-22\t\t0.276" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "2\t1.00e-03\t7.240e-10\t\t2.166e-22\t\t-0.332" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "3\t1.00e-04\t7.240e-11\t\t9.316e-23\t\t0.366" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "4\t1.00e-05\t7.240e-12\t\t2.209e-22\t\t-0.375" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "========================= PASS! =========================\n", + "You are awesome.\n", + "\n" + ] + }, + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 9, + "text": [ + "True" + ] + } + ], + "prompt_number": 9 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "def derBbc(chi):\n", + " Bbc, const = CongruousMagBC(mesh, np.array([1., 0., 0.]), chi)\n", + " bc = Utils.sdiag(prob.mesh.vol)*prob.mesh.faceDiv*prob._Pout.T*Bbc \n", + " dBbc = lambda x: Utils.sdiag(prob.mesh.vol)*prob.mesh.faceDiv*prob._Pout.T*const*Bbc*(np.inner(prob.mesh.vol, x))\n", + " return bc, dBbc" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 8 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "d_chi = 0.8*chi\n", + "# I am not sure why it does not give us second order accuracy\n", + "Tests.checkDerivative(derBbc, chi, dx=d_chi, num=5, plotIt=False)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "==================== checkDerivative ====================\n", + "iter\th\t\t|J0-Jt|\t\t|J0+h*dJ'*dx-Jt|\tOrder\n", + "---------------------------------------------------------\n", + "0\t1.00e-01\t2.212e+00\t\t8.604e-05\t\tnan" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "1\t1.00e-02\t2.212e-01\t\t8.604e-07\t\t2.000" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "2\t1.00e-03\t2.212e-02\t\t8.574e-09\t\t2.001" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "3\t1.00e-04\t2.212e-03\t\t9.857e-11\t\t1.939" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "4\t1.00e-05\t2.212e-04\t\t1.424e-11\t\t0.840" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "========================= PASS! =========================\n", + "Awesome, Seogi, just awesome.\n", + "\n" + ] + }, + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 9, + "text": [ + "True" + ] + } + ], + "prompt_number": 9 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "##################\n", + "# Test J\n", + "##################\n", + "\n", + "d_chi = 0.8*chi #np.random.rand(mesh.nCz)\n", + "d_sph_ind = spheremodel(mesh, 0., 0., -100., 50)\n", + "d_chi[d_sph_ind] = 0.02\n", + "\n", + "from SimPEG.Tests import checkDerivative\n", + "\n", + "derChk = lambda m: [prob.data.dpred(m), lambda mx: -prob.Jvec(chi, mx)]\n", + "print '\\n'\n", + "passed = checkDerivative(derChk, chi, plotIt=False, dx=d_chi, num=2)\n", + "\n", + "\n", + "# plt.pcolor(X, Y, dpred.reshape(X.shape, order='F'))\n", + "\n", + "# plt.show()" + ], + "language": "python", + "metadata": {}, + "outputs": [] + } + ], + "metadata": {} + } + ] +} \ No newline at end of file