Fix examples folder

This commit is contained in:
seogi committed 2014-02-27 09:15:41 -08:00
1 parent 3cba7845ca
commit 751d7098b6
4 files changed
+486 -7

No files matched your search

+1 -1
View File
@@ -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 <http://www.ngdc.noaa.gov/geomag-web/>`_).
+47
View File
@@ -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()
+38 -6
View File
@@ -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()
+400
View File
@@ -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": [
"[<matplotlib.lines.Line2D at 0x301f910>]"
]
},
{
"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+VlZWl5uZmTZgwQaLine truncated
"text": [
"<matplotlib.figure.Figure at 0x2d7f0d0>"
]
}
],
"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": {}
}
]
}