diff --git a/simpegPF/Dev/Intgrl_MAG_Inv_Driver.py b/simpegPF/Dev/Intgrl_MAG_Inv_Driver.py index 7d186f7c..39778559 100644 --- a/simpegPF/Dev/Intgrl_MAG_Inv_Driver.py +++ b/simpegPF/Dev/Intgrl_MAG_Inv_Driver.py @@ -43,29 +43,31 @@ if topofile == 'null': topo = np.c_[mkvc(Nx),mkvc(Ny),mkvc(Nz)] else: - topofile = np.genfromtxt(topofile,delimiter=' \n',dtype=np.str,skip_header=0) + topo = np.genfromtxt(topofile,skip_header=1) # Work with flat topogrphy for now -nullcell = np.ones(mesh.nC) +actv = PF.Magnetics.getActiveTopo(mesh,topo,'N') + +nC = int(sum(actv)) # Load model file if isinstance(mstart, float): - mstart = np.ones(mesh.nC) * mstart + mstart = np.ones(nC) * mstart else: mstart = Utils.meshutils.readUBCTensorModel(mstart,mesh) - + mstart = mstart[actv==1] # Get magnetization vector for MOF if magfile=='DEFAULT': - M_xyz = PF.Magnetics.dipazm_2_xyz(np.ones(mesh.nC) * M[0], np.ones(mesh.nC) * M[1]) + M_xyz = PF.Magnetics.dipazm_2_xyz(np.ones(nC) * M[0], np.ones(nC) * M[1]) else: M_xyz = np.genfromtxt(magfile,delimiter=' \n',dtype=np.str,comments='!') # Create forward operator -F = PF.Magnetics.Intrgl_Fwr_Op(mesh,B,M_xyz,rxLoc,'tmi') +F = PF.Magnetics.Intrgl_Fwr_Op(mesh,B,M_xyz,rxLoc,actv,'tmi') # Get distance weighting function wr = PF.Magnetics.get_dist_wgt(mesh,rxLoc,3.,np.min(mesh.hx)/4) diff --git a/simpegPF/Magnetics.py b/simpegPF/Magnetics.py index b3ed6cab..0a167674 100644 --- a/simpegPF/Magnetics.py +++ b/simpegPF/Magnetics.py @@ -502,7 +502,7 @@ def Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,flag): return d -def Intrgl_Fwr_Op(mesh,B,M,rxLoc,flag): +def Intrgl_Fwr_Op(mesh,B,M,rxLoc,actv,flag): """ Magnetic forward operator in integral form @@ -534,13 +534,30 @@ def Intrgl_Fwr_Op(mesh,B,M,rxLoc,flag): @author: dominiquef """ + # Find non-zero cells + inds = np.nonzero(actv)[0] + # Create active cell projector + P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), + shape=(mesh.nC, len(inds))) + + mcell = len(inds) + + # Create vectors of nodal location (lower and upper coners for each cell) xn = mesh.vectorNx; yn = mesh.vectorNy; zn = mesh.vectorNz; + + yn2,xn2,zn2 = np.meshgrid(yn[1:], xn[1:], zn[1:]) + yn1,xn1,zn1 = np.meshgrid(yn[0:-1], xn[0:-1], zn[0:-1]) + + Yn = P.T*np.c_[mkvc(yn1), mkvc(yn2)] + Xn = P.T*np.c_[mkvc(xn1), mkvc(xn2)] + Zn = P.T*np.c_[mkvc(zn1), mkvc(zn2)] ndata = rxLoc.shape[0] + # Convert Bdecination from north to cartesian D = (450.-float(B[1]))%360. @@ -549,9 +566,9 @@ def Intrgl_Fwr_Op(mesh,B,M,rxLoc,flag): # Pre-allocate space and create magnetization matrix if required if (flag=='tmi') | (flag == 'xyz'): # If assumes uniform magnetization direction - if M.shape != (mesh.nC,3): + if M.shape != (mcell,3): - print 'Magnetization vector must be 3*Ncells' + print 'Magnetization vector must be Nc x 3' return @@ -564,7 +581,7 @@ def Intrgl_Fwr_Op(mesh,B,M,rxLoc,flag): if flag == 'tmi': - F = np.zeros((ndata, mesh.nC)) + F = np.zeros((ndata, mcell)) # Projection matrix Ptmi = mkvc(np.r_[np.cos(np.deg2rad(B[0]))*np.cos(np.deg2rad(D)), @@ -572,10 +589,10 @@ def Intrgl_Fwr_Op(mesh,B,M,rxLoc,flag): np.sin(np.deg2rad(B[0]))],2).T; elif flag == 'xyz': - F = np.zeros((int(3*ndata), mesh.nC)) + F = np.zeros((int(3*ndata), mcell)) elif flag == 'full': - F = np.zeros((int(3*ndata), int(3*mesh.nC))) + F = np.zeros((int(3*ndata), int(3*mcell))) else: print """Flag must be either 'tmi' | 'xyz' | 'full', please revised""" @@ -589,7 +606,7 @@ def Intrgl_Fwr_Op(mesh,B,M,rxLoc,flag): count = -1; for ii in range(ndata): - tx, ty, tz = get_T_mat(xn,yn,zn,rxLoc[ii,:]) + tx, ty, tz = get_T_mat(Xn,Yn,Zn,rxLoc[ii,:]) if flag=='tmi': F[ii,:] = Ptmi.dot(np.vstack((tx,ty,tz)))*Mxyz @@ -609,7 +626,7 @@ def Intrgl_Fwr_Op(mesh,B,M,rxLoc,flag): count = progress(ii,count,ndata) - print "Done 100% ...forward modeling completed!!\n" + print "Done 100% ...forward operator completed!!\n" return F @@ -637,6 +654,7 @@ def get_T_mat(Xn,Yn,Zn,rxLoc): @author: dominiquef """ + eps = 1e-10 # add a small value to the locations to avoid /0 mcell = Xn.shape[0] @@ -644,15 +662,15 @@ def get_T_mat(Xn,Yn,Zn,rxLoc): Tx = np.zeros((1,3*mcell)) Ty = np.zeros((1,3*mcell)) Tz = np.zeros((1,3*mcell)) - - dz2 = rxLoc[2] - Zn[:,0] - dz1 = rxLoc[2] - Zn[:,1] + + dz2 = rxLoc[2] - Zn[:,0] + eps + dz1 = rxLoc[2] - Zn[:,1] + eps - dy2 = Yn[:,1] - rxLoc[1] - dy1 = Yn[:,0] - rxLoc[1] + dy2 = Yn[:,1] - rxLoc[1] + eps + dy1 = Yn[:,0] - rxLoc[1] + eps - dx2 = Xn[:,1] - rxLoc[0] - dx1 = Xn[:,0] - rxLoc[0] + dx2 = Xn[:,1] - rxLoc[0] + eps + dx1 = Xn[:,0] - rxLoc[0] + eps R1 = ( dy2**2 + dx2**2 ) R2 = ( dy2**2 + dx1**2 ) diff --git a/simpegPF/notebooks/SimPEG Tutorial - MAG Linear Problem.ipynb b/simpegPF/notebooks/SimPEG Tutorial - MAG Linear Problem.ipynb new file mode 100644 index 00000000..7d6ab4e0 --- /dev/null +++ b/simpegPF/notebooks/SimPEG Tutorial - MAG Linear Problem.ipynb @@ -0,0 +1,1290 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Objective:** \n", + "\n", + "In this tutorial we will create a simple magnetic problem from scratch using the SimPEG framework.\n", + "\n", + "We are using the integral form of the magnetostatic problem. In the absence of free-currents or changing magnetic field, magnetic material can give rise to a secondary magnetic field according to:\n", + "\n", + "$$\\vec b = \\frac{\\mu_0}{4\\pi} \\int_{V} \\vec M \\cdot \\nabla \\nabla \\left(\\frac{1}{r}\\right) \\; dV $$\n", + "\n", + "Where $\\mu_0$ is the magnetic permealitity of free-space, $\\vec M$ is the magnetization per unit volume and $r$ defines the distance between the observed field $\\vec b$ and the magnetized object. Assuming a purely induced response, the strenght of magnetization can be written as:\n", + "\n", + "$$ \\vec M = \\mu_0 \\kappa \\vec H_0 $$\n", + "\n", + "where $\\vec H$ is an external inducing magnetic field, and $\\kappa$ the magnetic susceptibility of matter.\n", + "As derived by Sharma 1966, the integral can be evaluated for rectangular prisms such that:\n", + "\n", + "$$ \\vec b(P) = \\mathbf{T} \\cdot \\vec H_0 \\; \\kappa $$\n", + "\n", + "Where the tensor matrix $\\bf{T}$ relates the three components of magnetization $\\vec M$ to the components of the field $\\vec b$:\n", + "\n", + "$$\\mathbf{T} =\n", + "\t \\begin{pmatrix}\n", + " \t\tT_{xx} & T_{xy} & T_{xz} \\\\\n", + "\t\tT_{yx} & T_{yy} & T_{yz} \\\\\n", + "\t\tT_{zx} & T_{zy} & T_{zz} \n", + "\t\\end{pmatrix} $$\n", + " \n", + "In general, we discretize the earth into a collection of cells, each contributing to the magnetic data such that:\n", + "\n", + "$$\\vec b(P) = \\sum_{j=1}^{nc} \\mathbf{T}_j \\cdot \\vec H_0 \\; \\kappa_j$$\n", + "\n", + "giving rise to a linear problem.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Using matplotlib backend: nbAgg\n", + "Populating the interactive namespace from numpy and matplotlib\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "C:\\Users\\dominiquef.MIRAGEOSCIENCE\\AppData\\Local\\Continuum\\Anaconda\\lib\\site-packages\\IPython\\kernel\\__init__.py:13: ShimWarning: The `IPython.kernel` package has been deprecated. You should import from ipykernel or jupyter_client instead.\n", + " \"You should import from ipykernel or jupyter_client instead.\", ShimWarning)\n" + ] + } + ], + "source": [ + "%matplotlib notebook\n", + "%pylab" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Efficiency Warning: Interpolation will be slow, use setup.py!\n", + "\n", + " python setup.py build_ext --inplace\n", + " \n" + ] + } + ], + "source": [ + "from SimPEG import *\n", + "import simpegPF as PF" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false, + "scrolled": true + }, + "outputs": [], + "source": [ + "# First we need to define the direction of the inducing field\n", + "# As a simple case, we pick a vertical inducing field of magnitude 50,000nT. \n", + "# From old convention, field orientation is given as an azimuth from North \n", + "# (positive clockwise) and dip from the horizontal (positive downward).\n", + "H0 = np.array(([90.,0.,50000.]))\n", + "\n", + "# Assume all induced so the magnetization M is also in the same direction\n", + "M = np.array([90,0])\n", + "\n", + "# Create a mesh\n", + "hxind = [(5, 20)]\n", + "hyind = [(5, 20)]\n", + "hzind = [(5, 10)]\n", + "\n", + "mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCC')\n", + "\n", + "# Assume flat topo for now, so all cells are active\n", + "nC = mesh.nC \n", + "actv = np.ones(nC)\n", + "\n", + "\n", + "# Create and array of observation points\n", + "xr = np.linspace(-20., 20., 20)\n", + "yr = np.linspace(-20., 20., 20)\n", + "X, Y = np.meshgrid(xr, yr)\n", + "Z = np.ones(X.size)*(mesh.vectorNz[-1]+1.) # Let just put the observation flat\n", + "\n", + "rxLoc = np.c_[Utils.mkvc(X.T), Utils.mkvc(Y.T), Utils.mkvc(Z.T)]\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now that we have all our spatial components, we can create our linear system. For a single location and single component of the data, the system would looks like this:\n", + "\n", + "$$ b_x =\n", + "\t\\begin{bmatrix}\n", + "\tT_{xx}^1 &... &T_{xx}^{nc} & T_{xy}^1 & ... & T_{xy}^{nc} & T_{xz}^1 & ... & T_{xz}^{nc}\\\\\n", + "\t \\end{bmatrix}\n", + "\t \\begin{bmatrix}\n", + "\t\t\\mathbf{M}_x \\\\ \\mathbf{M}_y \\\\ \\mathbf{M}_z\n", + "\t\\end{bmatrix} \\\\ $$\n", + "\n", + "where each of $T_{xx},\\;T_{xy},\\;T_{xz}$ are [nc x 1] long. For the $y$ and $z$ component, we need the two other rows of the tensor $\\mathbf{T}$.\n", + "In our simple induced case, the magnetization direction $\\mathbf{M_x,\\;M_y\\;,Mz}$ are known and assumed to be constant everywhere, so we can reduce the size of the system such that: \n", + "\n", + "$$ \\vec{\\mathbf{d}}_{\\text{pred}} = (\\mathbf{T\\cdot M})\\; \\kappa$$\n", + "\n", + "\n", + "\n", + "In most geophysical surveys, we are not collecting all three components, but rather the magnitude of the field, or $Total\\;Magnetic\\;Intensity$ (TMI) data.\n", + "Because the inducing field is really large, we will assume that the anomalous fields are parallel to $H_0$:\n", + "\n", + "$$ d^{TMI} = \\hat H_0 \\cdot \\vec d$$\n", + "\n", + "We then end up with a much smaller system:\n", + "\n", + "$$ d^{TMI} = \\mathbf{F\\; \\kappa}$$\n", + "\n", + "where $\\mathbf{F} \\in \\mathbb{R}^{nd \\times nc}$ is our $forward$ operator." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Begin calculation of forward operator: tmi\n", + "Done 0.0 %\n", + "Done 10.0 %\n", + "Done 20.0 %\n", + "Done 30.0 %\n", + "Done 40.0 %\n", + "Done 50.0 %\n", + "Done 60.0 %\n", + "Done 70.0 %\n", + "Done 80.0 %\n", + "Done 90.0 %\n", + "Done 100% ...forward operator completed!!\n", + "\n" + ] + } + ], + "source": [ + "# First, convert the magnetization direction to Cartesian\n", + "mi = np.ones(mesh.nC) * M[0]\n", + "md = np.ones(mesh.nC) * M[1]\n", + "M_xyz = PF.Magnetics.dipazm_2_xyz( mi , md ) # Ouputs an nc x 3 array\n", + "\n", + "# Create the forward model operator\n", + "F = PF.Magnetics.Intrgl_Fwr_Op(mesh,H0,M_xyz,rxLoc,actv,'tmi')" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " fig.waiting = false;\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " this.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width);\n", + " canvas.attr('height', height);\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('